# A New Tau Method for Solving Nonlinear Lane-Emden Type Equations via Bernoulli Operational Matrix of Differentiation

Authors: E. Tohidi; Kh. Erfani; M. Gachpazan; S. Shateyi
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

A new and efficient numerical approach is developed for solving nonlinear Lane-Emden type equations via Bernoulli operational matrix of differentiation. The fundamental structure of the presented method is based on the Tau method together with the Bernoulli polynomial approximations in which a new operational matrix is introduced. After implementation of our scheme, the main problem would be transformed into a system of algebraic equations such that its solutions are the unknown Bernoulli coefficients. Also, under several mild conditions the error analysis of the proposed method is provided. Several examples are included to illustrate the efficiency and accuracy of the proposed technique and also the results are compared with the different methods. All calculations are done in Maple 13.


---

## Body

## 1. Introduction

In recent years, the researches on the singular initial value problems (SIVPs) in several special second-order ordinary differential equations (ODEs) have received considerable attention among mathematicians and physicists. One of the most well-known classes of such equations are the Lane-Emden type equations which model many phenomena in mathematical physics and astrophysics. They are nonlinear ordinary differential equations which describe the equilibrium density distribution in self-gravitating sphere of polytrophic isothermal gas and have a singularity at the origin [1]. It must be noted that these equations have fundamental importance in the field of radiative cooling and modeling of clusters of galaxies. Moreover, it has been recently observed that the density profiles of dark matter halos are often modeled by the isothermal Lane-Emden equation with suitable boundary conditions at the origin [2]. Since getting the analytic solution of these equations is a difficult task in many cases, robust numerical schemes must be constructed for obtaining the approximated solutions. In this paper, we will present an efficient method for computing the numerical solution of the Lane-Emden type equations [3, 4]
(1)y′′(x)+αxy′(x)+f(x,y(x))=g(x),x∈[0,1],α>0,
with the initial conditions
(2)y(0)=a,y′(0)=0,
where the prime denotes the differentiation with respect to x, a is a constant, f and g are nonlinear continuous functions. Selecting α=2, f(x,y(x))=(y(x))n, g(x)=0 and a=1 yields [5, 6]
(3)y′′(x)+2xy′(x)+yn(x)=0,x∈[0,1],
which has another form
(4)1x2ddx(x2dydx)+yn(x)=0,x∈[0,1],
subject to the initial conditions
(5)y(0)=1,y′(0)=0.
Similarly, isothermal gas spheres equation (in the case of f(x,y(x))=ey(x)) is modeled by [5]
(6)y′′(x)+2xy′(x)+ey(x)=0,x∈[0,1],
with the zero values of y(0) and y′(0).Recently, many analytical and numerical methods have been used to solve Lane-Emden type equations (1), (3), and (6). We note that the main difficulty arises in the singularity of the equations at x=0. In [7, 8], Wazwaz employed the Adomian decomposition method with an alternate framework designed to overcome the difficulty of the singular point. The authors of [9] also applied pseudospectral method based on rational Legendre functions. Ramos [10–13] and Yousefi [14] solved Lane-Emden equation through several powerful numerical methods. Also, Yıldırım and Öziş [15] presented approximate exact solutions of a class of Lane-Emden type singular IVPs problems, by the VIM.Since the beginning of 1994, the Bernoulli, Bernstein, Legendre, Taylor, Fourier, and Bessel matrix methods have been used in the works [2, 14, 16–26] to solve high-order linear and nonlinear differential (including hyperbolic partial differential equations) Fredholm Volterra integro-differential difference delay equations and their systems. The main characteristic of these approaches is based on the operational matrices of differentiation instead of integration. However, we can use high order Gauss quadrature rules such as [27, 28] for problems in different integration forms. To the best of our knowledge, this is the first work concerning the Bernoulli matrix method for solving nonlinear SIVPs such as Lane-Emden type equations (1), (3), and (6). This partially motivated our interest in such method. Therefore, in the light of the above-mentioned methods and by means of the matrix relations between the Bernoulli polynomials and their derivatives together with the idea of the Tau scheme, we develop a new method for solving nonlinear Lane-Emden type equations.The Bernoulli polynomials play an important role in several branches of mathematical analysis, such as the theory of distributions in p-adic analysis [29], the theory of modular forms [30], the study of polynomial expansions of analytic functions [31], and so forth. Recently, new applications of the Bernoulli polynomials have also been found in mathematical physics, in connection with the theory of the Korteweg-de Vries equation [32] and Lamé equation [33], and in the study of vertex algebras [34].In this paper, we are concerned with the direct solution technique for solving (1), (3), and (6) (with the associated initial conditions) using the Tau method based on Bernoulli operational matrix, such that it can be implemented efficiently and at the same time has a good convergence property. It should be noted that, the Tau method originate from vanishing the inner product of the residual and each of the test functions. By the aid of Tau scheme, the mentioned equations would be transformed into the systems of algebraic equations.The remainder of the paper is organized as follows. In Section2 we introduce some mathematical preliminaries of Bernoulli polynomials and also their operational matrix of differentiation. Section 3 is devoted to applying the Tau method for solving Lane-Emden type equations using the Bernoulli operational matrix. In Section 4 we provide the error analysis of the proposed method. In Section 5 the proposed method is applied to several examples together with a full comparison to other recent methods. Also conclusions are given in Section 6.


## 2. Bernoulli Polynomials and Their Operational Matrix of Differentiation

Bernoulli polynomials and their applications can be found in number theory and classical analysis initially. They also appear in the integral representation of the differentiable periodic functions, since they are employed for approximating such functions in terms of polynomials. They are also used for representing the remainder term of the composite Euler-Maclaurin quadrature rule [16].In this section, we recall some properties of the Bernoulli polynomials which will be of fundamental importance in the sequel.Property 1 (differentiation: see [16]).
B
n
′
(
x
)
=
n
B
n
-
1
(
x
), n=1,2,….Property 2 (integral means conditions: see [16]).
∫
0
1
B
n
(
x
)
d
x
=
0, n=1,2,….Property 3 (differences: see [16]).
B
n
(
x
+
1
)
- Bn(x)=nxn-1,  n=1,2,….Property 4 (monomials representation: see [16]).
B
n
(
x
)
=
∑
k
=
0
n
‍
(
n
k
)
B
k
(
0
)
x
n
-
k,  n=1,2,….If we introduce the Bernoulli vectorB(x) in the form B(x)=[B0(x),B1(x),…,BN(x)], then the derivative of the B(x), with the aid of the first property, can be expressed in the matrix form by
(7)[B0′(x)B1′(x)B2′(x)⋮BN-1′(x)BN′(x)]︸B′(x)T=[000…000100…000020…000⋮⋮⋮⋱⋮⋮⋮000…N-100000…0N0]︸M[B0(x)B1(x)B2(x)⋮BN-1(x)BN(x)]︸BT(x),
where M is the (N+1)×(N+1) operational matrix of differentiation.Accordingly, thekth derivative of B(x) can be given by
(8)B′(x)T=MB(x)T⟹B(1)(x)=B(x)MT,B(2)(x)=B(1)(x)MT=B(x)(MT)2,B(3)(x)=B(1)(x)(MT)2=B(x)(MT)3,B(3)(x)⋮B(k)(x)=B(x)(MT)k,
where M is defined in (7).


## 3. Implementation of the Tau Scheme Using Bernoulli Operational Matrix

In this section, we apply the Tau method for solving the Lane-Emden type equations numerically. Therefore, we first assume that the solution of the Lane-Emden type equations can be written in terms of linear combination of Bernoulli polynomials (with unknown coefficients). Then, inner product of the residual with each of the test functions should be assumed to zero. Thus, the considered problem would be transformed into a system of nonlinear algebraic equations in which its solutions are the associated Bernoulli coefficients. Again consider the Lane-Emden equation(9)y′′(x)+αxy′(x)+f(x,y)=g(x),α,x≥0
with the initial conditions
(10)y(0)=a,y′(0)=0.
Now we approximate y(x), f(x,y(x)) and g(x) by Bernoulli polynomials in the following forms
(11)y(x)≈yN(x)=∑i=0N‍ciBi(x)=CBT(x),f(x,y)≈f(x,yN(x))=f(x,CBT(x)),g(x)≈gN(x)=∑i=0N‍giBi(x)=GBT(x),
where the unknowns are C=[c0c1⋯cN] meanwhile G=[g0g1⋯gN] are known and B(x)=[B0(x)B1(x)⋯BN(x)] are the Bernoulli polynomials.Using operational matrix of differentiation of Bernoulli polynomials, (1) can be written approximately as follows:
(12)CM2BT(x)+αxCM1BT(x)+f(x,CBT(x))≈GBT(x).
The residual RN(x) for (12) can be written as
(13)RN(x)=CM2BT(x)+αxCM1BT(x)+f(x,CBT(x))-GBT(x).
Applying the typical Tau method, which is used in the sense of a particular form of the Petrov-Galerkin method [4], (12) can be converted in N-1 nonlinear equations by applying
(14)〈RN(x),Bi(x)〉=∫01‍(RN(x)Bi(x))dx=0,∫01‍(RN(x)Bi(x))dx=01i=0,1,…,N-2.
The initial conditions are given by
(15)y(0)=CTB(0)=a,y(1)(0)=CTMB(0)=0.
Equations (14) and (15) generate N+1 sets of nonlinear equations in terms of N+1 unknown Bernoulli coefficients ci, i=0,1,…,N. These nonlinear equations can be solved using Newton’s iterative method that was implemented in fsolve command of MAPLE software for obtaining the unknown coefficients of the row vector C and hence the solution y(x) can be approximated easily by yN(x).


## 4. Error Analysis and Accuracy of the Solution

In this section, we will illustrate convergence of the proposed method assuming the known functions and also the unknown solution are in the spaceCm[0,1] (where m is the degree of smoothness of the problem) with bounded derivatives. But some lemmas should be recalled from the literature and then the main theorem of this section would be provided.Lemma 1.
Assume thatg∈L2[0,1] is an enough smooth function and also is approximated by the Bernoulli series ∑n=0∞‍gnBn(x), then the coefficients gn for all n=0,1,…,∞ can be calculated from the following relation:
(16)gn=1n!∫01‍g(n)(x)dx.Proof.
See [20].Lemma 2.
Assume that one approximates the functiong(x) on the interval [0,1] by Bernoulli polynomials as argued in Lemma 1. Then the coefficients gn decays as follows:
(17)gn≤Gnn!,
where Gn=supx∈[0,1]|g(n)(x)|.Proof.
Since it is trivial we omit the proof.The previous lemma implies that Bernoulli coefficients are decayed rapidly as increasing ofn under the condition of boundedness of all derivatives of g(x) with respect to n!. We will illustrate this fact experimentally in our numerical examples. It should be noted that, there are some examples such as g(x)=1/(2x+5)5 which has bounded derivatives, but Gn/n! does not go to zero rapidly.Theorem 3 (see [37]).
Suppose thatg(x) is in the space Cm[0,1] and is approximated by Bernoulli polynomials as done in Lemma 1. With more details assume that PN[g](x) is the approximate polynomial of g(x) in terms of Bernoulli polynomials and RN[g](x) is the remainder term. Then, the associated formulas are stated as follows:
(18)g(x)=PN[g](x)+RN[g](x),x∈[0,1],PN[g](x)=∫01‍g(x)dx+∑j=1N‍Bj(x)j!(g(j-1)(1)-g(j-1)(0)),RN[g](x)=-1N!∫01‍BN*(x-t)g(N)(t)dt,
where BN*(x)=BN(x-[x]) and [x] denotes the largest integer not greater than x.Proof.
See [37].Lemma 4.
Supposeg(x)∈C∞[0,1] (with bounded derivatives) and gN(x) is the approximated polynomial using Bernoulli polynomials. Then the error bound would be obtained as follows:
(19)∥E(gN(x))∥∞≤CG^(2π)-N,x∈[0,1],
where G^ denotes a bound for all the derivatives of function g(x) (i.e., ∥g(i)(x)∥∞≤G^, for i=0,1,…) and C is a positive constant.Proof.
See [25].Next, we will provide the main theorem of this section. We show that, if both of they(x) and g(x) are approximated by the Bernoulli polynomials in the Lane-Emden equation (1), then the error of the approximation of y(x) depends directly on the approximation of g(x). Therefore, using enough values of N may gives us high order approximation of the desired solutions.Theorem 5.
Assume thatF(x,y(x))=x-α∫0x‍tαf(t,y(t))dt and G(x)=a+Lα(g(x)), where Lα(·)=∫0x‍x-α∫0x‍tα(·)dtdx is a linear integral operator. If we approximate y(x) and G(x) by yN(x) and GN(x), respectively, using Bernoulli polynomials, then
(20)∥y(x)-yN(x)∥∞≤11-LF∥G(x)-GN(x)∥∞,
where LF is the Lipschitz constant of the function F(x,y(x)) with respect to its second variable y(x) and also LF≪1.Proof.
Consider the integral operatorLα(·)=∫0x‍x-α∫0x‍tα(·)dtdx. Applying Lα to both sides of (1) yields
(21)y(x)=a+Lα(g(x))+Lα(f(x,y(x))).
Considering the following assumptions:
(22)G(x)=a+Lα(g(x)),F(x,y(x))=x-α∫0x‍tαf(t,y(t))dt,
transforms (21) into the following Volterra integral equations:
(23)y(x)=G(x)+∫0x‍F(t,y(t))dt,x∈[0,1].
Now approximate both functions y(x) and G(x) by yN(x) and GN(x) using Bernoulli polynomials. Therefore
(24)yN(x)=GN(x)+∫0x‍F(t,yN(t))dt,x∈[0,1].
Thus
(25)∥y(x)-yN(x)∥∞=∥G(x)-GN(x)+∫0x‍(F(t,y(t))-F(t,yN(t)))dt∥∞≤∥G(x)-GN(x)∥∞+LF∥y(x)-yN(x)∥∞.
In other words
(26)∥y(x)-yN(x)∥∞≤11-LF∥G(x)-GN(x)∥∞
and this completes the proof.


## 5. Numerical Illustrations

In this section, several numerical examples are given to illustrate the accuracy and effectiveness of the proposed method. All calculations are done on a 64 bits PC laptop. As we claimed in Section2, our method which is based on the Bernoulli polynomials has more efficiency with respect to the Legendre method [4] in isothermal gas spheres equation (i.e., Example 3) and also has the same accuracy with respect to methods which use high-order Gauss quadrature rules [36]. We must recall that, the results of our method and the method of [36] are the same in Examples 5 and 6, but our results were obtained in less CPU time; meanwhile the results of [36] were obtained in more CPU time by using the same PC and equipment. Moreover in Table 5, we show the vanishing of the Bernoulli coefficients (as shown theoretically in Lemma 2) for all the examples of this section experimentally.Example 1 (see [4]).
At first we consider the equationy′′(x)+(2/x)y′(x)+1=0 with the initial conditions y(0)=1,  y′(0)=0. Trivially the exact solution of this equation is y(x)=1-(x2/6).
Note thatn=0, f(x,y)=y0(x)=1 and g(x)=0. We apply the method that was explained in Section 3 for N=2. Thus assume that
(27)y(x)≈y2(x)=CBT(x)=[c0c1c2][B0(x)B1(x)B2(x)]=c0B0(x)+c1B1(x)+c2B2(x).
Our aim is to determine the unknown Bernoulli coefficientsc0, c1, and c2 by using Tau method. Also the operational matrix of Bernoulli polynomials and its square are as follows:
(28)M1=[000100020],M2=[000000200].
By using (14) we have
(29)〈R2(x),B0(x)〉=∫01‍(R2(x)B0(x))dx=0,
where
(30)R2(x)=xy2′′(x)+2y2′(x)+x=x[c0c1c2]M2BT(x)+2[c0c1c2]MBT(x)+x=x[c0c1c2][000000200][B0(x)B1(x)B2(x)]+2[c0c1c2][000100020][B0(x)B1(x)B2(x)]+x=2c2x+2c1+4c2(x-12)+x.
Therefore
(31)∫01‍(2c2x+2c1+4c2(x-12)+x)dx=2c1+c2+12=0.
By applying the initial conditions we have
(32)c0-12c1+16c2=1,c1-c2=0.
Solving (31)-(32) yields c0=17/18, c1=-1/6, c2=-1/6. Thus
(33)y2(x)=c0B0(x)+c1B1(x)+c2B2(x)=(1718-16-16)(1x-12x2-x+16)=1-x26,
which is the exact solution.Example 2 (see [4]).
We now consider the following Lane-Emden equation:(34)y′′(x)+2xy′(x)+y(x)=0,
with initial conditions y(0)=1, y′(0)=0 which has the exact solution y(x)=sin(x)/x.
Again, we apply our proposed method for solving the above equation for different values ofN such as N=8 and N=10. In other words we assume that y(x)≈yN(x)=CBT(x) and f(x,y)=y(x)≈yN(x)=CBT(x). The numerical results of our method are depicted in Figure 1. In the Figure 1(a), the exact solution y(x) together with the numerical solutions y8(x) and y10(x) is illustrated in a larger computational interval. Moreover the error history of these approximated solutions is shown in the Figure 1(b). From this figure one can conclude that our method obtained highly accurate solutions even in large computational intervals.(a) Comparison of the exact and approximated solutions and (b) history of the errors.
(a)
(b)Example 3 (see [1, 4, 7]).
As the third example we consider the isothermal gas spheres equation as follows:(35)y′′(x)+2xy′(x)+ey(x)=0,x≥0
with the initial conditions y(0)=0 and y′(0)=0.
In this case we havef(x,y)=ey and g(x)=0. We approximate f by using the five terms of its Maclaurin expansion (i.e., f(x,y)≈1+y+y2/2+y3/6+y4/24). Again we use yN(x) instead of y(x) in the procedure of approximation. In other words
(36)f(x,y)≈1+CBT(x)+(CBT(x))22+(CBT(x))36+(CBT(x))424.
We apply our method for solving this problem by using several values of N. We provide the numerical solutions at the points x=0.0,0.1,0.2,0.5, and 1.0 in the case of N=10 in Table 1 together with the exact solution that was reported in [7]. Moreover we make an interesting comparison with the methods that are based on Legendre operational matrix of differentiation [4] and the Hermite collocation method [1] in the associated values of the errors by the assumption of N=10. From these comparisons we see that the presented method (PM) has more efficiency with respect to the above-mentioned methods.Table 1
The exact and approximated solution of Example3 together with the error comparisons for N=10.
x
Solution of PM
Exact Solution [7]
Error of PM
Error of [4]
Error of [1]
0.0
+
3.01060
e
-
31
+
0.0000000000
3.01
e
-
31
9.24
e
-
18
0.00
e
+
00
0.1
-
0.0016658339
-
0.0016658338
6.21
e
-
11
5.28
e
-
10
5.85
e
-
07
0.2
-
0.0066533671
-
0.0066533671
9.22
e
-
13
3.37
e
-
08
6.04
e
-
07
0.5
-
0.0411539573
-
0.0411539567
5.95
e
-
10
8.12
e
-
06
5.58
e
-
07
1.0
-
0.1588276831
-
0.1588272857
3.97
e
-
07
4.93
e
-
04
8.20
e
-
07Example 4 (see [1, 4, 7]).
As the fourth example we consider the following Lane-Emden equation:(37)y′′(x)+2xy′(x)+sin(y)=0,
with the initial conditions y(0)=1 and y′(0)=0.
Similar to the third example we havef(x,y)=sin(y) and g(x)=0. Approximating f by using the three terms of its Maclaurin expansion yields f(x,y)=y-y3/6+y5/120 and also by applying the Bernoulli polynomials in the procedure of approximations we have
(38)f(x,y)≈CBT(x)-(CBT(x))36+(CBT(x))5120.
Again the method that is proposed in Section 3 will be used. The numerical solutions at the points x=0.0,0.1,0.2,0.5, and 1.0 in the case of N=10 are written in Table 2 together with the exact solution that was reported in [7]. Moreover the error of the presented method (PM) together with the errors of [1, 4] by the assumption of N=10 is provided. Evidently the results of our methods are very close to [4] and are superior with respect to the results that were obtained in [1].Table 2
The exact and approximated solution of Example4 together with the error comparisons for N=10.
x
Solution of PM
Exact Solution [7]
Error of PM
Error of [4]
Error of [1]
0.0
0.9999999999
1.0000000000
1.11
e
-
16
1.11
e
-
16
0.00
e
+
00
0.1
0.9985976023
0.9985979273
3.25
e
-
07
3.25
e
-
07
7.21
e
-
06
0.2
0.9943949769
0.9943962648
1.28
e
-
06
1.28
e
-
06
1.00
e
-
05
0.5
0.9651702487
0.9651777797
7.53
e
-
06
7.53
e
-
06
1.04
e
-
05
1.0
0.8636571676
0.8636807315
2.35
e
-
05
2.35
e
-
05
7.03
e
-
06Example 5 (see [1, 35, 36]).
We now consider the Lane-Emden equationy′′(x)+(2/x)y′(x)+y3(x)=0 with the initial conditions y(0)=1, y′(0)=0.
The exact solution of this equation was reported in [35]. The numerical results of our scheme together with two other methods [1, 36] are provided in Table 3. Not only does our method needs to lower values of Bernoulli polynomials (as the test functions) with respect to the methods that are based on Gauss quadrature rules [36] for obtaining high accurate solutions, but also the CPU time of our method for solving the associated system of nonlinear algebraic is mush less than others such as [36]. The basic reason for this claim is based upon using the operational matrices of differentiation and this fact leads to sparse equations.Table 3
The exact and approximated solution of Example5 together with the error comparisons.
x
Solution of PM
Exact Solution [35]
Error of PMN=10
Error of [1]
Error of [36] N=25
0.0
1.0000000
1.0000000
0.0
e
+
0.0
0.0
e
+
0.0
0.0
e
+
0.0
0.1
0.9983358
0.9983358
2.95
e
-
08
1.40
e
-
06
2.95
e
-
08
0.5
0.9598390
0.9598391
3.00
e
-
08
2.99
e
-
06
3.00
e
-
08
1.0
0.8550575
0.8550576
3.14
e
-
08
1.99
e
-
06
3.14
e
-
08Example 6 (see [1, 35, 36]).
As the final example we consider the Lane-Emden equationy′′(x)+(2/x)y′(x)+y4(x)=0 with initial conditions y(0)=1, y′(0)=0.
Again we recall that the exact solution of this equation was reported in [35]. The numerical results of this Example are given in Table 4. Similar to the previous example we can see that the results of our approach are the same [36] by using lower values of Bernoulli polynomials with respect to the shifted Legendre polynomials and also are superior with respect to the results of [1].Table 4
The exact and approximated solution of Example6 together with the error comparisons.
x
Solution of PM
Exact Solution [35]
Error atN=10
Error of [1]
Error [36] N=25
0.0
1.0000000000
1.0000000000
0.0
e
+
0.0
0.0
e
+
0.0
0.0
e
+
0.0
0.1
0.9983366595
0.9983367000
4.05
e
-
08
2.50
e
-
04
4.04
e
-
08
0.2
0.9933862135
0.9933862000
1.35
e
-
08
2.48
e
-
04
1.35
e
-
08
0.5
0.9603109023
0.9603109000
2.33
e
-
09
2.05
e
-
04
2.34
e
-
09
1.0
0.8608138122
0.8608138000
1.22
e
-
08
1.93
e
-
04
1.22
e
-
08Table 5
Vanishing of the Bernoulli coefficients as the number of index is increased.
N
=
8
Example2
Example3
Example4
Example5
Example6
c
0
+
0.946083
+
0.948962
+
0.951426
+
0.954020
-
0.053961
c
1
-
0.158529
-
0.144942
-
0.133975
-
0.136343
-
0.158828
c
2
-
0.150584
-
0.126065
-
0.108253
-
0.132283
-
0.151451
c
3
+
0.015700
+
0.035408
+
0.046539
+
0.008232
+
0.014332
c
4
+
0.007379
+
0.012099
+
0.011798
+
0.004384
+
0.006137
c
5
-
0.000558
-
0.005891
+
0.011145
+
0.000293
-
0.001159
c
6
-
0.000174
-
0.000918
-
0.001069
+
0.000050
-
0.000277
c
7
+
0.000011
+
0.001152
+
0.003566
+
0.000041
+
0.000106
c
8
+
0.000002
-
0.000076
-
0.001424
-
0.000010
+
0.000010HintAs we claimed in Lemma2 the Bernoulli coefficients in the procedure of approximating any arbitrary continuous function must tend to zero as the number of the test functions (i.e., the Bernoulli polynomials) tend to infinity. This fact is shown experimentally in Table 5.


## 6. Conclusions

The Lane-Emden type equations describe a variety of phenomena in theoretical physics and astrophysics, including the aspects of stellar structure, the thermal history of a spherical cloud of gas, isothermal gas spheres, and thermionic currents. Lane-Emden type equations have been considered by many mathematicians as mentioned before [4]. The fundamental goal of this paper has been to construct an approximation to the solution of nonlinear Lane-Emden type equations in the computational interval [0,1]. A set of Bernoulli polynomials is proposed to provide an effective but simple way to improve the convergence of the solution by the Tau method. The validity of the method is based on the assumption that it converges by increasing the number of Bernoulli polynomials. A comparison is made among the numerical (and exact) solutions of [35] and the series solutions of [7], Legendre operational matrix [4], Hermite collocation method [1] and high-order Gauss quadrature rule [36], and the current work. It has been shown that the present work provides acceptable approach for Lane-Emden type equations.



---

*Source: 850170-2013-05-29.xml*