# A Nonclassical Radau Collocation Method for Nonlinear Initial-Value Problems with Applications to Lane-Emden Type Equations

Authors: Mohammad Maleki; M. Tavassoli Kajani; I. Hashim; A. Kilicman; K. A. M. Atan
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

We propose a numerical method for solving nonlinear
initial-value problems of Lane-Emden type. The method is based upon
nonclassical Gauss-Radau collocation points, and weighted interpolation.
Nonclassical orthogonal polynomials, nonclassical Radau points
and weighted interpolation are introduced on arbitrary intervals. Then
they are utilized to reduce the computation of nonlinear initial-value
problems to a system of nonlinear algebraic equations. We also present
the comparison of this work with some well-known results and show
that the present solution is very accurate.


---

## Body

## 1. Introduction

Many problems in the literature of mathematical physics can be formulated as equations of the Lane-Emden type defined in the form(1.1)y′′+2xy′+g(y(x))=0,0<x<∞,
subject to
(1.2)y(0)=A,y′(0)=0,
where prime denotes differentiation with respect to x. The solution of the Lane-Emden equation, as well as those of a variety of nonlinear problems in quantum mechanics and astrophysics such as the scattering length calculations in the variable phase approach, is numerically challenging because of the singular point at the origin. Equations (1.1) and (1.2) with specializing g(y(x)) and A occur in several models of mathematical physics and astrophysics such as the theory of stellar structure, the thermal behavior of a spherical cloud of gas, and theories of thermionic currents. It has been studied widely in the literature; see, for example, [1–12]. This equation was first studied by the astrophysicists Jonathan Homer Lane and Robert Emden, which considered the thermal behavior of a spherical cloud of gas acting under the mutual attraction of its molecules and subject to the classical laws of thermodynamics [13]. For g(y)=yα, and A=1 in (1.1) and (1.2), we obtain the standard Lane-Emden equation of index α which has been the object of much study [1–3]. It was physically shown that interesting values of α lie in the interval [0,5], and this equation has analytical solutions for α=0, 1 and 5. Various alternative techniques have been developed for the solution of the Lane-Emden equation in the literature. Among others, this equation has been solved by means of perturbation methods and a (1-1) Padé approximation (Bender et al. [2]), Adomian's decomposition method [3], the quasilinearization method of [4], the homotopy analysis method [5, 6], a variational approach which uses a semi-inverse method to obtain a variational principle [7], a power series solution [8], a linearization technique [9, 10], the variational iteration method [11, 12], hybrid functions collocation method [14], Lagrangian method [15], Hermite functions collocation method [16], Sinc-Collocation method [17], rational Legendre pseudospectral approach [18], a modified Legendre-spectral method [19], and a numerical technique based on converting the Lane-Emden equations into integral equations [20].In the present paper, we first consider the nonlinear ordinary differential equations of the form(1.3)f(x,y(x),y′(x),y′′(x))=0,0⩽x<∞
with initial conditions
(1.4)y(0)=c1,y′(0)=c2.
We assume that (1.3) and (1.4) have a unique solution y(x) to be determined. We then solve a variety of Lane-Emden equations which fall into this category. Here, we introduce a new direct computational method for solving (1.3) and (1.4). This method consists of reducing the solution of (1.3) and (1.4) to a set of algebraic equations by first interpolating y(x) using weighted Lagrange interpolation based on Freud-type weights and sets of nonclassical Gauss-Radau (NGR) nodes. These nodes, which arise from nonclassical orthogonal polynomials based on Freud-type weights over interval [a,T], are presented. Equation (1.3) is then collocated at these NGR collocation points to evaluate the unknown coefficients, which are the values of the function y(x) at these collocation points.This paper is organized as follows. In Section2, we describe the generation of NGR collocation points, function approximation, and selection of weights. In Section 3, we explain our method, and in Section 4, the present method is applied to a nonlinear Lane-Emden equation as well as the standard Lane-Emden equation of index α. The numerical solutions are compared in Section 5 with available exact or approximate solutions in order to assess the accuracy of the proposed method.


## 2. Nonclassical Radau Collocation Method

### 2.1. NGR Points

In classical pseudospectral methods [14, 22–24], the classical Gauss-Lobatto and Gauss-Radau collocation points are based on Chebyshev or Legendre polynomials and lie on the closed interval [-1,1] and half-open interval [-1,1), respectively. In the present work, we consider the generation of the NGR collocation points based on nonclassical orthogonal polynomials with respect to exponential weights in the intervals [a,T] where a and T are any real numbers.LetN+1 be the number of collocation points and let PN(t) be the Nth-degree nonclassical orthogonal polynomial with respect to weight w(t) which can be obtained from the following three-term recurrence relation [25]:
(2.1)Pk+1(x)=(x-αk)Pk(x)-βkPk-1(x),k=0,1,2,…,P-1(x)=0,P0(x)=1.
The recurrence coefficients in (2.1) are given in [26] by
(2.2)αk=∫aT‍tw(x)Pk2(x)dx∫aT‍w(x)Pk2(x)dx,k=0,1,2,…,β0=∫aT‍w(x)P02(x)dx,βk=∫aT‍w(x)Pk2(x)dx∫aT‍w(x)Pk-12(x)dx.
The NGR collocation points xj for j=0,1,…,N are obtained by the method outlined by Golub [27]. The tridiagonal Jacobi-Radau matrix of order N+1 is defined by
(2.3)JN+1R=(α0β1β1α1β2β2α2β3⋱⋱⋱βN-1αN-1βNβNαNR),
where αNR is obtained from
(2.4)αNR=a-βNPN-1(a)PN(a).Theorem 2.1 (see Golub [27]).
The Gauss-Radau nodesa=x0<x1<⋯<xN<T are the eigenvalues of JN+1R, and the Gauss-Radau weights wj are given by
(2.5)wj=β0(v1j)2,j=0,1,…,N,
where vj is the normalized eigenvector of JN+1R corresponding to the eigenvalue xj (i.e., vjTvj=1) and v1j its first component.


### 2.2. Function Approximation and Differentiation Matrices

Consider the NGR collocation points(x0,x1,…,xN) defined in Section 2.1 on the interval [a,T) and additional noncollocated point xN+1=T. The function y∈L2[a,T] is approximated by weighted Lagrange interpolation as [26, 28]
(2.6)y(x)≈yN(x)=∑j=0N+1‍W(x)W(xj)Lj(x)yj,x∈[a,T],
where W is a positive and bounded weight function, and  yj=y(xj),
(2.7)Lj(x)=∏k=0k≠jN+1‍x-xkxj-xk,j=0,1,…,N+1,
is a basis of (N+1)th-degree Lagrange polynomials. Notice that the basis includes the function LN+1 corresponding to the terminal value xN+1=T. Differentiating twice the series of (2.6) and evaluating the NGR collocation points xi, 0⩽i⩽N give
(2.8)y′(xi)≈∑j=0N+1‍yjddx[W(x)W(xj)Lj(x)]x=xi=Di(1)·Y,y′′(xi)≈∑j=0N+1‍yjd2dx2[W(x)W(xj)Lj(x)]x=xi=Di(2)·Y,
where
(2.9)Y=[y0,y1,…,yN+1]T,Di(1) and Di(2) are the ith row of D(1) and D(2),
(2.10)Dij(1)=ddx[W(x)W(xj)Lj(x)]x=xi,Dij(2)=d2dx2[W(x)W(xj)Lj(x)]x=xi.
The rectangular (N+1)×(N+2) matrices D(1) and D(2) formed by the coefficients Dij(1) and Dij(2), i=0,1,…,N; j=0,1,…,N+1 are the first and second order Gauss-Radau differentiation matrices, respectively. These matrices transform the approximation of y(x) at x0,…,xN+1 to the first and second derivatives of y at the collocation points x0,…,xN.


### 2.3. Weight Selection

When studying the uniform convergence behavior of the weighted Lagrange interpolation asN→∞, a crucial role is played by the Lebesgue function(2.11)λN,W(x)=W(x)∑j=0N+1‍|Lj(x)|W(xj),
and the  Lebesgue constant(2.12)ΛN,W=∥λN,W(x)∥∞,
(see, e.g., [29–32] and references therein).In general, the orthogonal weight functionw and weight of interpolation W are chosen independently [26, 28, 33–35]. Nevertheless, if we expect a reasonable upper estimate for the Lebesgue constant, then we have to assume some connections between these two weights. The most natural assumption, as suggested in [29], is to generate the interpolation nodes (collocation points) xj with respect to the weight w=W2. In addition, to obtain uniform convergence in weighted interpolation some conditions on W should be considered. The Bernstein's approximation problem deals with the uniform convergence behavior of weighted interpolation. The problem is as followsLetW>0 be measurable. When it is true that for every continuous y:ℝ→ℝ with
(2.13)lim|x|→∞W(x)y(x)=0,
there exist a sequence of polynomials {FN}N=1∞ with
(2.14)EN(y)W=limN→∞∥W(y-FN)∥=0?
The condition lim|x|→∞W(x)y(x)=0 is essential to counteract the growth of any polynomial at infinity.As stated above, the Lebesgue constant plays an important role in answering the Bernstein's problem. Now it is known that for any set of interpolation nodes and any weightW, ΛN,W is unbounded with respect to N. A consequence of this is that there exists y such that weighted interpolation does not converge uniformly to y. However, if y is not too badly behaved (e.g., as measured by the modulus of continuity) and the ΛN,W are not too large, then uniform convergence is achieved (as positive answers to the Bernstein’s problem).The Freud-type weights, as positive answers to the Bernstein’s problem, (see, e.g., [36]) are defined as W(t)=exp(-Q(t)) with the following conditions: Q:ℝ→ℝ is even, continuous in ℝ, Q′′ is continuous in (0,∞), Q′>0 on (0,∞) and for some 1<a≤b(2.15)a≤(xQ′(x))′Q′(x)≤b,x∈(0,∞).IfW is a Freud weight we write W∈ℱ, and if, moreover, Q is differentiable in ℝ, we write w∈ℱ1. By definition, if W∈ℱ, then W2∈ℱ, too. Canonical example is as
(2.16)Wβ(x)=exp(-|x|β),β≥1,x∈ℝ.
Clearly Wβ∈ℱ1 with a=b=β. In this work, we consider the Freud weight with β=1.Vértesi in [32] has shown that given a Freud-type weight satisfying the previous conditions and any set of interpolation nodes
(2.17)logN=O(ΛN,W),N→∞.


## 2.1. NGR Points

In classical pseudospectral methods [14, 22–24], the classical Gauss-Lobatto and Gauss-Radau collocation points are based on Chebyshev or Legendre polynomials and lie on the closed interval [-1,1] and half-open interval [-1,1), respectively. In the present work, we consider the generation of the NGR collocation points based on nonclassical orthogonal polynomials with respect to exponential weights in the intervals [a,T] where a and T are any real numbers.LetN+1 be the number of collocation points and let PN(t) be the Nth-degree nonclassical orthogonal polynomial with respect to weight w(t) which can be obtained from the following three-term recurrence relation [25]:
(2.1)Pk+1(x)=(x-αk)Pk(x)-βkPk-1(x),k=0,1,2,…,P-1(x)=0,P0(x)=1.
The recurrence coefficients in (2.1) are given in [26] by
(2.2)αk=∫aT‍tw(x)Pk2(x)dx∫aT‍w(x)Pk2(x)dx,k=0,1,2,…,β0=∫aT‍w(x)P02(x)dx,βk=∫aT‍w(x)Pk2(x)dx∫aT‍w(x)Pk-12(x)dx.
The NGR collocation points xj for j=0,1,…,N are obtained by the method outlined by Golub [27]. The tridiagonal Jacobi-Radau matrix of order N+1 is defined by
(2.3)JN+1R=(α0β1β1α1β2β2α2β3⋱⋱⋱βN-1αN-1βNβNαNR),
where αNR is obtained from
(2.4)αNR=a-βNPN-1(a)PN(a).Theorem 2.1 (see Golub [27]).
The Gauss-Radau nodesa=x0<x1<⋯<xN<T are the eigenvalues of JN+1R, and the Gauss-Radau weights wj are given by
(2.5)wj=β0(v1j)2,j=0,1,…,N,
where vj is the normalized eigenvector of JN+1R corresponding to the eigenvalue xj (i.e., vjTvj=1) and v1j its first component.


## 2.2. Function Approximation and Differentiation Matrices

Consider the NGR collocation points(x0,x1,…,xN) defined in Section 2.1 on the interval [a,T) and additional noncollocated point xN+1=T. The function y∈L2[a,T] is approximated by weighted Lagrange interpolation as [26, 28]
(2.6)y(x)≈yN(x)=∑j=0N+1‍W(x)W(xj)Lj(x)yj,x∈[a,T],
where W is a positive and bounded weight function, and  yj=y(xj),
(2.7)Lj(x)=∏k=0k≠jN+1‍x-xkxj-xk,j=0,1,…,N+1,
is a basis of (N+1)th-degree Lagrange polynomials. Notice that the basis includes the function LN+1 corresponding to the terminal value xN+1=T. Differentiating twice the series of (2.6) and evaluating the NGR collocation points xi, 0⩽i⩽N give
(2.8)y′(xi)≈∑j=0N+1‍yjddx[W(x)W(xj)Lj(x)]x=xi=Di(1)·Y,y′′(xi)≈∑j=0N+1‍yjd2dx2[W(x)W(xj)Lj(x)]x=xi=Di(2)·Y,
where
(2.9)Y=[y0,y1,…,yN+1]T,Di(1) and Di(2) are the ith row of D(1) and D(2),
(2.10)Dij(1)=ddx[W(x)W(xj)Lj(x)]x=xi,Dij(2)=d2dx2[W(x)W(xj)Lj(x)]x=xi.
The rectangular (N+1)×(N+2) matrices D(1) and D(2) formed by the coefficients Dij(1) and Dij(2), i=0,1,…,N; j=0,1,…,N+1 are the first and second order Gauss-Radau differentiation matrices, respectively. These matrices transform the approximation of y(x) at x0,…,xN+1 to the first and second derivatives of y at the collocation points x0,…,xN.


## 2.3. Weight Selection

When studying the uniform convergence behavior of the weighted Lagrange interpolation asN→∞, a crucial role is played by the Lebesgue function(2.11)λN,W(x)=W(x)∑j=0N+1‍|Lj(x)|W(xj),
and the  Lebesgue constant(2.12)ΛN,W=∥λN,W(x)∥∞,
(see, e.g., [29–32] and references therein).In general, the orthogonal weight functionw and weight of interpolation W are chosen independently [26, 28, 33–35]. Nevertheless, if we expect a reasonable upper estimate for the Lebesgue constant, then we have to assume some connections between these two weights. The most natural assumption, as suggested in [29], is to generate the interpolation nodes (collocation points) xj with respect to the weight w=W2. In addition, to obtain uniform convergence in weighted interpolation some conditions on W should be considered. The Bernstein's approximation problem deals with the uniform convergence behavior of weighted interpolation. The problem is as followsLetW>0 be measurable. When it is true that for every continuous y:ℝ→ℝ with
(2.13)lim|x|→∞W(x)y(x)=0,
there exist a sequence of polynomials {FN}N=1∞ with
(2.14)EN(y)W=limN→∞∥W(y-FN)∥=0?
The condition lim|x|→∞W(x)y(x)=0 is essential to counteract the growth of any polynomial at infinity.As stated above, the Lebesgue constant plays an important role in answering the Bernstein's problem. Now it is known that for any set of interpolation nodes and any weightW, ΛN,W is unbounded with respect to N. A consequence of this is that there exists y such that weighted interpolation does not converge uniformly to y. However, if y is not too badly behaved (e.g., as measured by the modulus of continuity) and the ΛN,W are not too large, then uniform convergence is achieved (as positive answers to the Bernstein’s problem).The Freud-type weights, as positive answers to the Bernstein’s problem, (see, e.g., [36]) are defined as W(t)=exp(-Q(t)) with the following conditions: Q:ℝ→ℝ is even, continuous in ℝ, Q′′ is continuous in (0,∞), Q′>0 on (0,∞) and for some 1<a≤b(2.15)a≤(xQ′(x))′Q′(x)≤b,x∈(0,∞).IfW is a Freud weight we write W∈ℱ, and if, moreover, Q is differentiable in ℝ, we write w∈ℱ1. By definition, if W∈ℱ, then W2∈ℱ, too. Canonical example is as
(2.16)Wβ(x)=exp(-|x|β),β≥1,x∈ℝ.
Clearly Wβ∈ℱ1 with a=b=β. In this work, we consider the Freud weight with β=1.Vértesi in [32] has shown that given a Freud-type weight satisfying the previous conditions and any set of interpolation nodes
(2.17)logN=O(ΛN,W),N→∞.


## 3. Solution of Nonlinear Initial-Value Problems

Our discrete approximation to the nonlinear initial-value problem in (1.3)-(1.4) is obtained by evaluating (2.6) at N collocation points x1,…,xN and replacing y'(xi) and y′′(xi) by their discrete approximations in (2.8), and evaluating the boundary conditions in (1.4) at collocation point x0=0. Hence, the discrete approximation to the nonlinear initial-value problem is
(3.1)f(xi,yi,Di(1)·Y,Di(2)·Y)=0,i=1,2,…,N,(3.2)y(x0)=y0=c1,D0(1)·Y=c2.Using (3.1) and (3.2), we obtain a system of N+2 nonlinear algebraic equations which can be solved using the Newton's iterative method.It is well known that the initial guess for Newton's iterative method is very important especially for complicated problems. To choose the initial guess for our problem, in the first stage we setN=5 and apply the Newton's iterative method for solving N+2 nonlinear equations by choosing c1 in (1.4) as our initial guess. We then increase N by 5 and use the approximate solution in stage one as our initial guess in this stage. We continue this approach until the results are similar up to a required number of decimal places for two consecutive stages.It is worth mentioning that, in the case that the initial-value problem has a singularity atx=0 (e.g., the Lane-Emden type equations), this method avoids the singularity, because we compute (3.1) at the collocation points that are straightly more than zero.


## 4. Illustrative Examples

We applied the method presented in this paper and solved two problems. The first example is a nonlinear Lane-Emden equation considered in [10] and the second example is the standard Lane-Emden equation of index α. As stated in Section 2.3, we consider the Freud-type weight W(x)=exp(-|x|) as interpolation weight function and generate the nonclassical orthogonal polynomials and the NGR collocation points in the intervals [0,T] with respect to the weight w=W2.Example 4.1.
This example corresponds to the following singular nonlinear Lane-Emden equation [10]:
(4.1)y′′+2xy'+4(2ey+ey/2)=0,0<x<∞,y(0)=0,y′(0)=0,
which has the following exact solution:
(4.2)ye(x)=-2ln(1+x2).
Define
(4.3)EN,T=max{|ye(x)-yN(x)|:0⩽x⩽T},E~N,T=max{|ye′(x)-yN′(x)|:0⩽x⩽T},
where yN(x) and ye(x) denote the approximate solution obtained by the present method and the exact solution, respectively. In Table 1, the maximum absolute errors between approximate and exact solutions are denoted by EN,T, and maximum absolute error between the derivative of approximate and exact solutions is denoted by E~N,T, for different values of N and T are given, which show the efficiency of the present method in large interval calculation. Further, in Table 2 a comparison is made between the values of y(x) obtained using the present method for T=10 and N=30 together with the values given in [16] using Hermite functions collocation method and the exact solutions.Table 1
The maximum absolute errorsEN,T and E~N,T for Example 4.1.
T
N
E
N
,
T
E
~
N
,
T
20
20
1.80
×
10
-
4
4.40
×
10
-
4
25
1.15
×
10
-
5
1.50
×
10
-
4
30
1.25
×
10
-
6
9.20
×
10
-
6
10
20
8.50
×
10
-
6
1.25
×
10
-
4
25
4.10
×
10
-
7
1.24
×
10
-
5
30
3.65
×
10
-
8
4.50
×
10
-
7
5
20
1.85
×
10
-
7
4.30
×
10
-
6
25
5.40
×
10
-
9
1.83
×
10
-
7
30
1.60
×
10
-
10
7.50
×
10
-
9Table 2
Comparison ofy(x) between present method, method in [16] and the exact values for Example 4.1.
x
Method in [16]
Error
Present method
Error
Exact
0.00
0.0000000000
0.00
×
1
0
0
0.0000000000
0.00
×
1
0
0
0.0000000000
0.01
−0.0001970587
2.93
×
1
0
-
6
−0.0001999904
4.08
×
1
0
-
10
−0.0001999900
0.10
−0.0198967225
3.94
×
1
0
-
6
−0.0199006640
2.27
×
1
0
-
9
−0.0199006617
0.50
−0.4462840851
3.02
×
1
0
-
6
−0.4462870994
3.28
×
1
0
-
9
−0.4462871026
1.00
−1.3862934297
9.31
×
1
0
-
7
−1.3862943490
1.21
×
1
0
-
8
−1.3862943611
2.00
−3.2188763248
5.00
×
1
0
-
7
−3.2188758411
1.63
×
1
0
-
8
−3.2188758249
3.00
−4.6051709964
8.10
×
1
0
-
7
−4.6051702055
1.95
×
1
0
-
8
−4.6051701860
4.00
−5.6664274573
7.69
×
1
0
-
7
−5.6664266969
8.80
×
1
0
-
9
−5.6664266881
5.00
−6.5161937402
6.64
×
1
0
-
7
−6.5161930821
6.09
×
1
0
-
9
−6.5161930760
6.00
−7.2218363729
5.48
×
1
0
-
7
−7.2218358310
5.68
×
1
0
-
9
−7.2218358253
7.00
−7.8240461812
1.70
×
1
0
-
7
−7.8240460154
4.57
×
1
0
-
9
−7.8240460109
8.00
−8.3487734467
1.09
×
1
0
-
7
−8.3487745418
1.97
×
1
0
-
9
−8.3487745398
9.00
−8.8134506165
1.21
×
1
0
-
5
−8.8134384925
2.05
×
1
0
-
9
−8.8134384945
10.00
−9.2302027821
3.83
×
1
0
-
5
−9.2302410337
1.18
×
1
0
-
11
−9.2302410337Example 4.2.
Consider the standard Lane-Emden equation of indexα given by
(4.4)y′′+2xy′+yα=0,0<x<∞,y(0)=1,y′(0)=0,
This equation is linear for α=0 and 1, nonlinear otherwise, and exact solutions exist only for α=0,1 and 5 and are given in Bender et al. [2], respectively, by
(4.5)y(x)=1-x26,y(x)=sin(x)x,y(x)=(1+x23)-12.
Moreover, Bender et al. [2] determined the zeros of y(x) asymptotically, here denoted by ξ, and found that
(4.6)ξ=π+0.885273956δ+0.24222δ2,
for δ=-0.5,0,0.5,1.0 and 1.5 which correspond to α=0,1,2,3, and 4, respectively.We applied the method presented in this paper and solved this example and then evaluated the zeros ofy(x), which are also evaluated in [10] using linearization technique, in [15] by using Lagrangian method, and in [16] using Hermite functions collocation method.The selection ofT is crucial for the computing of zero, ξ, of y(x). In order to obtain reasonable approximations of zeros, for each value of α, in the first stage we select a sufficiently large T and set N=25 then we solve the problem using the present method to obtain ξ1 as the first approximation of the zero. Then according to the obtained ξ1 we select a value for T. Finally, we resolve the problem for different values of N.In Table3, the resulting values of the zeros of y(x) for α=0,0.5,…,4.5 using the present method for different N, with the results obtained in [10, 15, 16] together with the exact solutions of Horedt [21], are presented. Table 3 shows that the present method provides very accurate predictions of the zeros of y(x) even in large intervals. In order to demonstrate the accuracy of the proposed method, in Tables 4 and 5 we have compared the numerical results of y(x) using the present method for α=3, N=25 and α=4, N=30, respectively, with the methods in [15–17, 21]. In addition, Table 6 shows the maximum absolute errors for α=5, for which the exact solution exists, with N=30 and choosing different values of T. As can be seen from Tables 3-6, the present method provides very accurate results even for large values of T. The resulting graphs of the standard Lane-Emden equation for α=0,0.5,1,1.5,2,2.5,3,3.5,4 and 4.5 and different values of N given in Table 3 are shown in Figure 1.Table 3
Comparison of numerical results forξ (y(ξ)=0) for Example 4.2.
α
N
Present method
Method in [10]
Method in [15]
Method in [16]
Exact [21]
0
20
2.4494897427832
2.44899
—
—
2.4494897427832
0.5
27
2.752698051
—
—
—
2.75269805
1
20
3.1415926535898
3.14048
—
—
3.1415926535898
1.5
25
3.65375374
—
—
3.65375374
3.65375374
2
15
4.35287460
4.35086
4.352875
4.35287460
4.35287460
2.5
17
5.35527546
—
—
5.35527546
5.35527546
3
25
6.89684862
6.89312
6.896846
6.89684862
6.89684862
3.5
25
9.53580535
—
—
—
9.53580534
4
30
14.9715464
14.96518
14.971546
14.9715463
14.9715463
4.5
30
31.8364633
—
—
—
31.8364632Table 4
Comparison ofy(x) for present method, and other methods in the literature for α=3.
x
Present method
Error
Error in [17]
Error in [15]
Error in [16]
Exact [21]
0.000
1.0000000000
0.00×100
0.00×100
0.00×100
0.00×100
1.0000000000
0.100
0.9983358295
2.95×10-8
1.82×10-5
1.28×10-5
1.40×10-6
0.9983358000
0.500
0.9598390699
3.00×10-8
1.12×10-4
1.81×10-5
2.99×10-6
0.9598391000
1.000
0.8550575686
3.14×10-8
1.07×10-4
6.00×10-7
1.99×10-6
0.8550576000
5.000
0.1108198351
3.51×10-8
2.02×10-4
2.00×10-7
3.90×10-7
0.1108198000
6.000
0.0437379839
3.90×10-9
5.00×10-5
2.00×10-5
1.14×10-6
0.0437379800
6.800
0.0041677894
4.00×10-10
2.93×10-4
2.79×10-6
1.05×10-5
0.0041677890
6.896
0.0000360111
1.00×10-10
2.10×10-5
1.01×10-6
8.88×10-8
0.0000360112Table 5
Comparison ofy(x) for present method, and other methods in the literature for α=4.
x
Present method
Error
Error in [16]
Exact [21]
0.0
1.0000000000
0.00
×
1
0
0
0.00
×
1
0
0
1.0000000000
0.1
0.9983366596
4.04
×
1
0
-
8
2.51
×
1
0
-
4
0.9983367000
0.2
0.9933862134
1.34
×
1
0
-
8
2.48
×
1
0
-
4
0.9933862000
0.5
0.9603109024
2.36
×
1
0
-
9
2.05
×
1
0
-
4
0.9603109000
1.0
0.8608138124
1.24
×
1
0
-
8
1.93
×
1
0
-
4
0.8608138000
5.0
0.2359227314
3.14
×
1
0
-
8
8.59
×
1
0
-
5
0.2359227000
10.0
0.0596727418
1.80
×
1
0
-
9
6.22
×
1
0
-
5
0.0596727400
14.0
0.0083305270
5.30
×
1
0
-
12
2.47
×
1
0
-
5
0.0083305270
14.9
0.0005764180
9.09
×
1
0
-
10
5.00
×
1
0
-
7
0.0005764189Table 6
The maximum absolute errorsEN,T and E~N,T for α=5 and N=30.
T
E
N
,
T
E
~
N
,
T
40
8.0
×
1
0
-
6
1.0
×
1
0
-
5
30
5.5
×
1
0
-
8
6.4
×
1
0
-
7
20
6.6
×
1
0
-
9
1.1
×
1
0
-
7
10
1.4
×
1
0
-
10
3.5
×
1
0
-
9
5
1.5
×
1
0
-
13
3.2
×
1
0
-
12Figure 1
Lane-Emden graphs obtained by the present method forα=0,0.5,1,…,4.5.


## 5. Comparison with Other Methods

As stated in the  Introduction section, the Lane-Emden equations have been widely solved using both numerical and analytical methods. In this section we aim to present the advantages of our numerical method over some other existing methods in the literature.(i)
Comparison with some analytical solutions: among others, the Lane-Emden equations have been solved with the variational iteration method [11, 12], homotopy-perturbation method [5, 6], and the Adomian decomposition method [3]. However, this type of solution methods is dependent on the initial guess so that the obtained series solution is changed with changing the initial guess, whereas the present method is not dependent on the initial guess. Furthermore, in the mentioned methods the interval of convergence of the obtained series solution is limited (usually [0,1]), whereas Section 4 shows that our method provides accurate approximate solutions in larger domains. For instance, consider Example 4.1 of Section 4. The series solution for this example in [3, 6, 12] is as follows:
(5.1)y(x)=-2(x2-12x4+13x6-14x8+⋯),which converges to the exact solution only on the interval[0,1], while we have solved this example with high accuracy in the interval [0,20].Note that in the mentioned analytical solutions, Padé approximant can be implemented for manipulating a polynomial approximation into a rational function to gain more information about the approximate solution. Nevertheless, comparison between Figure 1 of Wazwaz [3] obtained using Padé approximants [6/6] with Figure 1 and Table 3 of this paper shows that our method provides much more accurate predictions of the zeros of the standard Lane-Emden equation in Example 4.2.(ii)
Comparison with some spectral methods: several spectral methods have been established for solving the Lane-Emden equations [15–20]. Methods in  [19, 20] can only be implemented in the interval [0,1], while our method can be implemented in larger interval. Further, methods in [15–18] are based on orthogonal functions on the semi-infinite interval [0,∞], such as Laguerre functions, Hermite functions, and radial basis functions. All of these methods need certain quadratures on unbounded domains, which introduce errors and so weaken the merit of spectral approximations. Moreover, for an infinitely smooth function f∈L2(0,+∞), the spectral convergence of the truncated series in these functions occurs only if f decays exponentially fast at +∞. According to these reasons, we see from Tables 2–5 that our method provides more accurate numerical results. Also, Table 3 shows that the spectral methods in [15–18] can solve the standard Lane-Emden equation up to the interval [0,15], whereas we have solved this equation with high accuracy up to the interval [0,32].(iii)
The present method has also the following advantages: first, for another type of equations the Freud-type weights described in Section2.3 can be tuned to improve the accuracy of the discrete approximation. Second, this method provides very accurate results with moderate number of collocation points even in large intervals. Third, methods in [26, 28, 33–35], by considering some arbitrary weight functions, have utilized nonclassical basis polynomials on the interval [0,1]. But, in the present work we have developed this idea over an arbitrary interval [0,T] based on the Freud-type weights.


## 6. Conclusion

The Lane-Emden equation occurs in the theory of stellar structure and describes the temperature variation of a spherical gas cloud. The difficulty in this type of equations, due to the existence of singular point atx=0, is overcome here. In the standard Lane-Emden equation, the first zero of y is an important point of the function, so we have computed y up to this zero by utilizing the nonclassical Radau collocation method. A set of nonclassical orthogonal polynomials based on Freud-type weights is proposed to provide an effective but simple way to improve the convergence of the solution by a Radau collocation method. Numerical examples demonstrate the validity and high accuracy of the technique.



---

*Source: 103205-2012-10-23.xml*