# A Successive Linearization Method Approach to Solve Lane-Emden

Authors: S. S. Motsa; S. Shateyi
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

We propose a new application of the successive linearization method for solving singular initial and boundary value problems of Lane-Emden type. To demonstrate the reliability of the proposed method, a comparison is made with results from existing methods in the literature and with exact analytical solutions. It was found that the method is easy to implement, yields accurate results, and performs better than some numerical methods.


---

## Body

## 1. Introduction

In this work, we investigate the application of a novel approach of solving a class of nonlinear singular initial value and boundary value problems in the second-order ordinary differential equations. In recent years, these problems have attracted the attention of many researchers because of their useful applications in astronomy, mathematical biology, physics, and other areas of science and engineering. One of the equations describing this type of equations is the Lane-Emden type equation formulated as(1.1)y′′(x)+2xy′(x)+yn=0,
with initial conditions
(1.2)y(0)=1,y′(0)=0,
where n∈[0,5] is a constant parameter. This equation is very useful in astrophysics in the study of polytropic models and stellar structures [1, 2]. For the special case when n=0,1,5, exact analytical solutions were obtained by Chandrasekhar [1]. For all other values of n, approximate analytical methods and numerical methods are used to approximate the solution of the Lane-Emden equation. Analytical approaches that have recently been applied in solving the Lane-Emden equations include the Adomian decomposition method [3, 4], differential transformation method [5], homotopy perturbation method [6], homotopy analysis method [7, 8], power series expansions [9–13], and variational iteration method [14, 15]. Generally, when all the above cited analytical approaches are used to solve Lane-Emden equation, a truncated power series solution of the true solution is obtained. This solution converges rapidly in a very small region (0<x<1). For x>1, convergence is very slow and the solutions are inaccurate even when using a large number of terms [16]. Convergence acceleration methods such as Padé approximations may be used to improve the convergence of the resulting series or to enlarge their domain of convergence. The homotopy analysis method [8, 17] has a unique advantage over the other analytic approximation methods because it has a convergence controlling parameter that can be adjusted to improve the region of convergence of the resulting series.An important physical parameter associated with the Lane-Emden function is the location of its first positive real zero. The first zero ofy(x) is defined as the smallest positive value x0 for which y(x0)=0. This value is important because it gives the radius of a polytropic star. The analytic approaches on their own are not very useful in solving for x0 because their region of convergence is usually less than x0. Recently, there has been a surge in the number of numerical methods that have been proposed to find solutions of the Lane-Emden equations. Recent numerical methods that have been proposed include the Legendre Tau method [18], Sinc-collocation method [19], the Lagrangian approach [20], the successive linearization method [21], and a numerical method based on radial basis functions [22]. Accurate results for the Lane-Emden equation have previously been reported in [23] where the Runge-Kutta routine with self-adapting step was used to generate seven digit tables of Lane-Emden functions. These tables are now widely used as a benchmark for testing the accuracy of new methods for solving the Lane-Emden equations. More recently, numerical values giving the first zero of the Lane-Emden equation for various values of n with at least ten decimal places of accuracy were tabulated in [24] who use a numerical perturbation approach and in [25] who use the Chebyshev spectral method.Another important class of singular ODE's of Lane-Emden type is the Frank-Kamenetskii [26] boundary value problem
(1.3)y′′(x)+kxy′(x)+δey=0,
which is solved subject to the boundary conditions
(1.4)y′(0)=0,y(1)=0.
This equation is used to model the temperature distribution in a vessel before a thermal explosion takes place. A substantial amount of research has been done on this type of problem using both analytical and numerical treatments. Closed form solutions of this problem were reported in [26, 27] for the case when k=0 and 1. Analytical treatment of (1.3) leading to series solutions was carried out in [28–30] using the Adomian decomposition method and in [31] using a power series approach based on symmetry reduction. Asaithambi [32] used a shooting method based approach that uses automatic differentiation to obtain a Taylor series expansion of the solution. Numerical approaches based on B-spline functions (see e.g, [28, 33, 34]) have been widely used to develop solutions of singular boundary value problems of the type (1.3). Kumar and Gupta [35] present a survey and review of papers of spline solution of singular boundary value problems of the type (1.3).In the present paper, we present a very simple and robust method that gives very accurate solutions to the nonlinear initial and boundary value problems governed by (1.1) and (1.3), respectively. The method extends the approach of [21] who attempted to solve the Lane-Emden equation (1.1) for the case when n=2,3,4 using the successive linearization method (SLM). The successive linearization method is a new numerical iterative approach that has been recently reported and successfully utilized in solving boundary value problems [36–39] occurring in boundary layer flows in the application of fluid mechanics.The SLM approach is based on transforming an ordinary nonlinear differential equation into an iterative scheme made up of linear equations which are then solved using numerical approaches. In [21], the SLM was used to obtain solutions of (1.1) that were accurate only up to seven decimal places. In this study, a modified version of the SLM is presented which converges much faster and gives solutions that are accurate to more than fifteen decimal places. We also extend the application of the present method to singular boundary value problems governed by the Frank-Kamenetskii equation. The accuracy of the present method it tested against results obtained using other methods was found in, literature.


## 2. Numerical Solution

In this section, we present the method used to solve the governing equations (1.1)–(1.4).


### 2.1. Solution of Lane-Emden Equation

To solve the Lane-Emden equation (1.1), it is convenient to recast the problem from being an initial value problem to a boundary value problem by considering only the domain x∈[0,α], where α is the first zero of y(x). In most practical applications of the Lane-Emden equation (1.1), the goal is to integrate the governing equation from 0 to α. Since α is an unknown parameter, we rescale the problem by setting
(2.1)x=αt.
Substituting (2.1) into (1.1) and simplifying gives
(2.2)ytt+2tyt+α2yn=0,y(0)=1,y(1)=0,yt(0)=0,
which is a nonlinear eigenvalue problem with α as the eigenvalue. To solve (2.2), we use the successive linearization method (SLM) (see e.g., [36–39]). The SLM approach is based on transforming an ordinary nonlinear differential equation into an iterative scheme made up of linear differential equations which are then solved using analytical or numerical methods wherever possible. In applying the SLM on (2.2), we set
(2.3)y(t)=yi(t)+∑m=0i-1‍ym(t),i=1,2,3,…,
where yi are unknown functions that are obtained by iteratively solving the linearized version of the governing equations assuming that yi(0≤m≤i-1) are known from previous iterations. The algorithm starts with an initial approximation y0(t) which is chosen to satisfy the boundary conditions in (2.2). A suitable initial guess that satisfies all the boundary conditions in (2.2) of this example is
(2.4)y0(t)=1-t2.
We observe that (2.2) is a second-order differential equation but has three boundary conditions. The third boundary condition (yt(0)=0) is used as an extra condition that is required to solve for the eigenvalue α. Since α is an unknown parameter, we also apply the SLM approach in solving for α by setting
(2.5)α=αi+∑m=0i-1‍αm.
An initial approximation for α that was found to give accurate results for all n∈(0,5) was
(2.6)α02=6(43)n.
This value of α0 was obtained by substituting the initial approximation y0 given by (2.4) in (2.2) and solving the resulting equation at x=1/2, the middle of the domain for t.Substituting the expansions (2.3) and (2.5) into (2.2) and neglecting nonlinear terms in yi and αi give
(2.7)y¨i+2ty˙i+a1,i-1yi+a2,i-1αi=ri-1,
subject to
(2.8)yi(0)=0,yi(1)=0,y˙i(0)=0,
where
(2.9)a1,i-1=n(∑m=0i-1‍αm)2(∑m=0n-1‍ym)n-1,a2,i-1=2∑m=0i-1‍αm(∑m=0n-1‍ym)n,ri-1=-[∑m=0i-1‍y¨m+2t∑m=0i-1‍y˙m+(∑m=0i-1‍αm)2(∑m=0i-1‍ym)n],
where the dots denote differentiation with respect to the scaled variable t. To solve (2.7), we use the Chebyshev collocation method in which case the functions yi are approximated by Chebyshev interpolating polynomials in terms of their values at the collocation points. Before applying the spectral method, it is convenient to transform the physical region [0,1] of problem 11 into the domain [−1,1]. This can be achieved by using the mapping
(2.10)t=z+12,-1≤z≤1.
We use the Gauss-Lobatto collocation points [40, 41] to define the Chebyshev nodes in [−1,1] as
(2.11)zj=cos(πjN),j=0,1,…N,
where (N+1) is the number of collocation points. The derivatives of yi at the collocation points can be written as
(2.12)y˙i(zj)=∑k=0N‍Djkyi(zk)=DYi,y¨i(zj)=∑k=0N‍Djk2yi(zk)=D2Yi,
where Yi=[yi(z0),yi(z1),…,yi(zN-1),yi(zN)]T, T is transpose, and D is the Chebyshev derivative matrix of size (N+1)×(N+1) whose entries are defined [40, 41] as
(2.13)Djk=cjck(-1)j+kzj-zkj≠k;j,k=0,1,…,N,Dkk=-zk2(1-zk2)k=1,2,…,N-1,D00=2N2+16=-DNN,
with
(2.14)ck={2,k=0,N,1,-1≤k≤N-1.
Applying the Chebyshev spectral collocation method in (2.7)-(2.8) gives
(2.15)Ai-1Yi+Bi-1αi=Ri-1,
subject to
(2.16)yi(zn)=0,yi(z0)=0,(2.17)∑k=0N‍DNkyi(zk)=0.
The matrices in (2.15) are defined as
(2.18)Ri-1=[ri-1(z0),ri-1(z1),…,ri-1(zN-1),ri-1(zN)]T,Ai-1=4D2+[4t]dD+a1,i-1,a1,i-1=[a1,i-1(t0)0⋯00a1,i-1(t1)⋯0⋮⋮⋮⋮00⋯a1,i-1(tN)],[4t]d=[4t00⋯004t1⋯0⋮⋮⋮⋮00⋯4tN],Bi-1=a2,i-1,witha2,i-1=[a2,i-1(t0)a2,i-1(t1)⋮a2,i-1(tN)].
We note that the system (2.15) consists of (N+1) unknowns in Yi and an additional unknown αi. To solve for all (N+2) unknowns, we increase the dimension of the system (2.15) by introducing an additional row on which we impose the Neumann boundary condition (2.17). The resulting system takes the form
(2.19)[Ai-1Bi-1DN,0⋯DN,N0][Yiαi]=[Ri-10].The boundary conditions (2.16) can be imposed on system (2.19) by deleting the first and (N+1)st rows of the coefficient matrix and vectors and deleting the first and (N+1)st columns of the coefficient matrix. Thus, starting from the initial approximations Y0 and α0, the solutions Yi and αi(i=1,2,3,…) can be determined from
(2.20)[Yiαi]=[Ai-1Bi-1DN,1⋯DN,N-10]-1[Ri-10].


### 2.2. Solution of the Frank-Kamenetskii Equation

The Frank-Kamenetskii equation admits exact analytical solutions whenk=1 (see e.g., [27, 42]). When δ<2, it has two solutions. The problem has no solution for δ>2 and when δ=2, it has a unique solution. The solutions are given by [42]
(2.21)y=log[16ec1(2δ+ec1x2)2],y=log[16ec2(1+2δec2x2)2],
where the constants c1 and c2 are given by
(2.22)c1=log[2(4-δ)±42(2-δ)],c2=log[4-δ±22(2-δ)2δ2].To solve (1.3), it is convenient to introduce the boundary condition at x=0 as
(2.23)y(0)=β,
where β is a constant. The governing equation (1.3) is then solved subject to the Dirichlet boundary conditions y(0)=β and y(1)=0 using the SLM approach as described in the previous subsection using the Neumann boundary condition y'(0)=0 as an extra condition for determining the unknown β. To this end, we look for a solution of the form
(2.24)y(x)=yi(x)+∑m=0i-1‍ym(x),β=βi+∑m=0i-1‍βm,
where yi and βi are obtained iteratively by solving the linearized equations that result from substituting (2.24) into the governing equation (1.3). Substituting (2.24) into (1.3)-(1.4) and neglecting nonlinear terms in yi and βi give
(2.25)yi′′+kxyi′+b1,i-1yi=si-1,
subject to
(2.26)yi(0)-βi=pi-1,yi(1)=0,yi′(0)=0,
where
(2.27)b1,i-1=δexp(∑m=0i-1‍ym),pi-1=∑m=0i-1‍βm-∑m=0i-1‍ym(0),si-1=-[∑m=0i-1‍ym′′+kx∑m=0i-1‍ym′+δexp(∑m=0i-1‍ym)].
A suitable initial guess that satisfies the boundary conditions (1.4) and (2.23) is
(2.28)y0(x)=β0(1-x2),
where β0 is the initial approximation for β which is given by
(2.29)δexp(3β04)=2β0(1+k).
Equation (2.29) is obtained by substituting (2.28) into the governing equation (1.3) and setting x=1/2. The solutions for β0 can be obtained using nonlinear equation solvers in scientific computing software such as the fsolve routine that is available in MAPLE or MATLAB.Equations (2.25)-(2.26) are solved using the Chebyshev spectral collocation method. The details of the implementation of this approach are similar to those given for the Lane-Emden equation in the previous section. The discretized equation system to be solved is
(2.30)Ci-1Yi=Si-1,∑k=0N‍DNkyi(zk)=0,yi(z0)=0,yi(zn)-βi=pi-1,
where
(2.31)Si-1=[si-1(z0),si-1(z1),…,si-1(zN)]T,Ci-1=4D2+[2kx]dD+b1,i-1,b1,i-1=[b1,i-1(x0)0⋯00b1,i-1(x1)⋯0⋮⋮⋮⋮00⋯b1,i-1(xN)],[2kx]d=2k[1x00⋯001x1⋯0⋮⋮⋮⋮00⋯1xN].
The equation system  can be written as the matrix equation
(2.32)[10⋯0000Ci-1⋮00⋯01-1DN,0DN,1⋯DN,N-1DN,N0][yi(z0)yi(z1)⋮yi(zN-1)yi(zN)βi]=[0si-1(z1)⋮si-1(zN-1)pi-10],
where the boundary conditions have been imposed on the first and last rows of Ci-1 and Si-1 and the derivative boundary condition has been added in the (N+2) row. Thus, starting from y0 and β0 given by (2.28) and (2.29), the subsequent solutions for yi,βi(i=1,2,3,…) can be obtained by solving the matrix system (2.32).


## 2.1. Solution of Lane-Emden Equation

To solve the Lane-Emden equation (1.1), it is convenient to recast the problem from being an initial value problem to a boundary value problem by considering only the domain x∈[0,α], where α is the first zero of y(x). In most practical applications of the Lane-Emden equation (1.1), the goal is to integrate the governing equation from 0 to α. Since α is an unknown parameter, we rescale the problem by setting
(2.1)x=αt.
Substituting (2.1) into (1.1) and simplifying gives
(2.2)ytt+2tyt+α2yn=0,y(0)=1,y(1)=0,yt(0)=0,
which is a nonlinear eigenvalue problem with α as the eigenvalue. To solve (2.2), we use the successive linearization method (SLM) (see e.g., [36–39]). The SLM approach is based on transforming an ordinary nonlinear differential equation into an iterative scheme made up of linear differential equations which are then solved using analytical or numerical methods wherever possible. In applying the SLM on (2.2), we set
(2.3)y(t)=yi(t)+∑m=0i-1‍ym(t),i=1,2,3,…,
where yi are unknown functions that are obtained by iteratively solving the linearized version of the governing equations assuming that yi(0≤m≤i-1) are known from previous iterations. The algorithm starts with an initial approximation y0(t) which is chosen to satisfy the boundary conditions in (2.2). A suitable initial guess that satisfies all the boundary conditions in (2.2) of this example is
(2.4)y0(t)=1-t2.
We observe that (2.2) is a second-order differential equation but has three boundary conditions. The third boundary condition (yt(0)=0) is used as an extra condition that is required to solve for the eigenvalue α. Since α is an unknown parameter, we also apply the SLM approach in solving for α by setting
(2.5)α=αi+∑m=0i-1‍αm.
An initial approximation for α that was found to give accurate results for all n∈(0,5) was
(2.6)α02=6(43)n.
This value of α0 was obtained by substituting the initial approximation y0 given by (2.4) in (2.2) and solving the resulting equation at x=1/2, the middle of the domain for t.Substituting the expansions (2.3) and (2.5) into (2.2) and neglecting nonlinear terms in yi and αi give
(2.7)y¨i+2ty˙i+a1,i-1yi+a2,i-1αi=ri-1,
subject to
(2.8)yi(0)=0,yi(1)=0,y˙i(0)=0,
where
(2.9)a1,i-1=n(∑m=0i-1‍αm)2(∑m=0n-1‍ym)n-1,a2,i-1=2∑m=0i-1‍αm(∑m=0n-1‍ym)n,ri-1=-[∑m=0i-1‍y¨m+2t∑m=0i-1‍y˙m+(∑m=0i-1‍αm)2(∑m=0i-1‍ym)n],
where the dots denote differentiation with respect to the scaled variable t. To solve (2.7), we use the Chebyshev collocation method in which case the functions yi are approximated by Chebyshev interpolating polynomials in terms of their values at the collocation points. Before applying the spectral method, it is convenient to transform the physical region [0,1] of problem 11 into the domain [−1,1]. This can be achieved by using the mapping
(2.10)t=z+12,-1≤z≤1.
We use the Gauss-Lobatto collocation points [40, 41] to define the Chebyshev nodes in [−1,1] as
(2.11)zj=cos(πjN),j=0,1,…N,
where (N+1) is the number of collocation points. The derivatives of yi at the collocation points can be written as
(2.12)y˙i(zj)=∑k=0N‍Djkyi(zk)=DYi,y¨i(zj)=∑k=0N‍Djk2yi(zk)=D2Yi,
where Yi=[yi(z0),yi(z1),…,yi(zN-1),yi(zN)]T, T is transpose, and D is the Chebyshev derivative matrix of size (N+1)×(N+1) whose entries are defined [40, 41] as
(2.13)Djk=cjck(-1)j+kzj-zkj≠k;j,k=0,1,…,N,Dkk=-zk2(1-zk2)k=1,2,…,N-1,D00=2N2+16=-DNN,
with
(2.14)ck={2,k=0,N,1,-1≤k≤N-1.
Applying the Chebyshev spectral collocation method in (2.7)-(2.8) gives
(2.15)Ai-1Yi+Bi-1αi=Ri-1,
subject to
(2.16)yi(zn)=0,yi(z0)=0,(2.17)∑k=0N‍DNkyi(zk)=0.
The matrices in (2.15) are defined as
(2.18)Ri-1=[ri-1(z0),ri-1(z1),…,ri-1(zN-1),ri-1(zN)]T,Ai-1=4D2+[4t]dD+a1,i-1,a1,i-1=[a1,i-1(t0)0⋯00a1,i-1(t1)⋯0⋮⋮⋮⋮00⋯a1,i-1(tN)],[4t]d=[4t00⋯004t1⋯0⋮⋮⋮⋮00⋯4tN],Bi-1=a2,i-1,witha2,i-1=[a2,i-1(t0)a2,i-1(t1)⋮a2,i-1(tN)].
We note that the system (2.15) consists of (N+1) unknowns in Yi and an additional unknown αi. To solve for all (N+2) unknowns, we increase the dimension of the system (2.15) by introducing an additional row on which we impose the Neumann boundary condition (2.17). The resulting system takes the form
(2.19)[Ai-1Bi-1DN,0⋯DN,N0][Yiαi]=[Ri-10].The boundary conditions (2.16) can be imposed on system (2.19) by deleting the first and (N+1)st rows of the coefficient matrix and vectors and deleting the first and (N+1)st columns of the coefficient matrix. Thus, starting from the initial approximations Y0 and α0, the solutions Yi and αi(i=1,2,3,…) can be determined from
(2.20)[Yiαi]=[Ai-1Bi-1DN,1⋯DN,N-10]-1[Ri-10].


## 2.2. Solution of the Frank-Kamenetskii Equation

The Frank-Kamenetskii equation admits exact analytical solutions whenk=1 (see e.g., [27, 42]). When δ<2, it has two solutions. The problem has no solution for δ>2 and when δ=2, it has a unique solution. The solutions are given by [42]
(2.21)y=log[16ec1(2δ+ec1x2)2],y=log[16ec2(1+2δec2x2)2],
where the constants c1 and c2 are given by
(2.22)c1=log[2(4-δ)±42(2-δ)],c2=log[4-δ±22(2-δ)2δ2].To solve (1.3), it is convenient to introduce the boundary condition at x=0 as
(2.23)y(0)=β,
where β is a constant. The governing equation (1.3) is then solved subject to the Dirichlet boundary conditions y(0)=β and y(1)=0 using the SLM approach as described in the previous subsection using the Neumann boundary condition y'(0)=0 as an extra condition for determining the unknown β. To this end, we look for a solution of the form
(2.24)y(x)=yi(x)+∑m=0i-1‍ym(x),β=βi+∑m=0i-1‍βm,
where yi and βi are obtained iteratively by solving the linearized equations that result from substituting (2.24) into the governing equation (1.3). Substituting (2.24) into (1.3)-(1.4) and neglecting nonlinear terms in yi and βi give
(2.25)yi′′+kxyi′+b1,i-1yi=si-1,
subject to
(2.26)yi(0)-βi=pi-1,yi(1)=0,yi′(0)=0,
where
(2.27)b1,i-1=δexp(∑m=0i-1‍ym),pi-1=∑m=0i-1‍βm-∑m=0i-1‍ym(0),si-1=-[∑m=0i-1‍ym′′+kx∑m=0i-1‍ym′+δexp(∑m=0i-1‍ym)].
A suitable initial guess that satisfies the boundary conditions (1.4) and (2.23) is
(2.28)y0(x)=β0(1-x2),
where β0 is the initial approximation for β which is given by
(2.29)δexp(3β04)=2β0(1+k).
Equation (2.29) is obtained by substituting (2.28) into the governing equation (1.3) and setting x=1/2. The solutions for β0 can be obtained using nonlinear equation solvers in scientific computing software such as the fsolve routine that is available in MAPLE or MATLAB.Equations (2.25)-(2.26) are solved using the Chebyshev spectral collocation method. The details of the implementation of this approach are similar to those given for the Lane-Emden equation in the previous section. The discretized equation system to be solved is
(2.30)Ci-1Yi=Si-1,∑k=0N‍DNkyi(zk)=0,yi(z0)=0,yi(zn)-βi=pi-1,
where
(2.31)Si-1=[si-1(z0),si-1(z1),…,si-1(zN)]T,Ci-1=4D2+[2kx]dD+b1,i-1,b1,i-1=[b1,i-1(x0)0⋯00b1,i-1(x1)⋯0⋮⋮⋮⋮00⋯b1,i-1(xN)],[2kx]d=2k[1x00⋯001x1⋯0⋮⋮⋮⋮00⋯1xN].
The equation system  can be written as the matrix equation
(2.32)[10⋯0000Ci-1⋮00⋯01-1DN,0DN,1⋯DN,N-1DN,N0][yi(z0)yi(z1)⋮yi(zN-1)yi(zN)βi]=[0si-1(z1)⋮si-1(zN-1)pi-10],
where the boundary conditions have been imposed on the first and last rows of Ci-1 and Si-1 and the derivative boundary condition has been added in the (N+2) row. Thus, starting from y0 and β0 given by (2.28) and (2.29), the subsequent solutions for yi,βi(i=1,2,3,…) can be obtained by solving the matrix system (2.32).


## 3. Results

In this section, we present the approximate solutions of the governing equations (1.1)–(1.4) using the successive linearization method (SLM). In order to assess the performance and reliability of the present method of solution, the results are tabulated in Tables 1–6 and compared with accurate results from the literature and exact analytical solutions for the cases where analytical solutions are available.Table 1
Comparison of the first zero (α) of the Lane-Emden equation at different iterations for different values of n against numerical results given by Boyd [25] and Seidov [24].
i
n
n
=
1
n
=
2
n
=
3
n
=
4
1
3.103790966018498
4.035737655330261
5.557521240863945
8.688676175146224
2
3.141338575299632
4.330814323152805
6.592739714660735
11.771827678456009
3
3.141592635714092
4.352749271550132
6.877461391335034
13.886237798409815
4
3.141592653589793
4.352874592230670
6.896766842576802
14.817460536386396
5
3.141592653589793
4.352874595946125
6.896848617928429
14.967975594173165
6
3.141592653589793
4.352874595946125
6.896848619376960
14.971544376671904
7
3.141592653589793
4.352874595946125
6.896848619376960
14.971546348837493
8
3.141592653589793
4.352874595946125
6.896848619376960
14.971546348838095
Reference [25]
3.141592653589793
4.352874595946125
6.896848619376960
14.971546348838095
Reference [24]
3.141592653589793
4.352874595946124
6.896848619376960
14.971546348838093Table 2
Comparison of the first zero (α) of the Lane-Emden equation at different iterations for different values of n against numerical results given by Boyd [25].
Iteration
n
n
=
0.5
n
=
1.5
n
=
2.5
n
=
3.5
1
2.748427032416
3.524247954347
4.685837835928
6.799314270510
2
2.752696116339
3.650158732919
5.264679963832
8.585590955706
3
2.752698054064
3.653750290899
5.353307428656
9.386863426519
4
2.752698054065
3.653753736216
5.355274556453
9.531659161527
5
2.752698054065
3.653753736219
5.355275459011
9.535802061278
6
2.752698054065
3.653753736219
5.355275459011
9.535805344243
7
2.752698054065
3.653753736219
5.355275459011
9.535805344245
8
2.752698054065
3.653753736219
5.355275459011
9.535805344245
Reference [25]
2.752698054065
3.653753736219
5.355275459011
9.535805344245Table 3
Comparison of the value ofdy/dt(1), for the Lane-Emden equation at different iterations for different values of n against numerical results given by Boyd [25].
i
n
n
=
1
n
=
2
n
=
3
n
=
4
1
-
1.032146158998257
-
0.614886114475695
-
0.316160520906889
-
0.009814787486055
2
-
1.000237675841930
-
0.558083545440036
-
0.306609857182704
-
0.146741321451808
3
-
1.000000012692448
-
0.553918189951109
-
0.293448460329348
-
0.128608843668195
4
-
1.000000000000000
-
0.553897421500438
-
0.292634978984266
-
0.121103290612396
5
-
1.000000000000000
-
0.553897420877304
-
0.292631615214246
-
0.120066832081784
6
-
1.000000000000000
-
0.553897420877304
-
0.292631615154770
-
0.120043051612478
7
-
1.000000000000000
-
0.553897420877304
-
0.292631615154770
-
0.120043038478707
8
-
1.000000000000000
-
0.553897420877304
-
0.292631615154770
-
0.120043038478703
Reference [25]
-
1.000000000000000
-
0.553897420877304
-
0.292631615154770
-
0.120043038478703Table 4
SLM results fory(0) corresponding to the two roots of (1.3) when k=1, δ=1, and N=50.
Iteration
1st root
2nd root
1
0.31669044212643981117
4.12191170750356255708
2
0.31669436763842469104
3.85609488634651084646
3
0.31669436764074987779
3.84275755885741291875
4
0.31669436764074987779
3.84218880566051216841
5
0.31669436764074987779
3.84218871571893306220
6
0.31669436764074987779
3.84218871571892197872
7
0.31669436764074987779
3.84218871571892197872
Exact
0.31669436764074987779
3.84218871571892197872Table 5
Maximum absolute errors for solution of (1.3) (1st root).
N
Iter.
1
2
3
Reference [33]
20
3.926
×
10
-
6
2.325
×
10
-
12
2.488
×
10
-
22
3.161
×
10
-
5
40
3.926
×
10
-
6
2.325
×
10
-
12
5.976
×
10
-
25
7.874
×
10
-
6
60
3.926
×
10
-
6
2.325
×
10
-
12
5.976
×
10
-
25
3.501
×
10
-
6
90
3.926
×
10
-
6
2.325
×
10
-
12
5.976
×
10
-
25
1.552
×
10
-
6
111
3.926
×
10
-
6
2.325
×
10
-
12
5.976
×
10
-
25
1.039
×
10
-
6Table 6
Maximum absolute errors for solution of (1.3) (2nd root).
N
Iter.
2
4
6
8
20
2.460
×
10
-
2
1.094
×
10
-
7
6.407
×
10
-
10
6.407
×
10
-
10
40
2.482
×
10
-
2
1.093
×
10
-
7
8.611
×
10
-
19
8.611
×
10
-
19
60
2.480
×
10
-
2
1.094
×
10
-
7
5.013
×
10
-
28
5.515
×
10
-
28
90
2.482
×
10
-
2
1.094
×
10
-
7
5.014
×
10
-
29
3.300
×
10
-
38
111
2.482
×
10
-
2
1.094
×
10
-
7
5.014
×
10
-
29
3.300
×
10
-
38In the case of the Lane-Emden equation (1.1), exact analytical solutions are only available for the case when n=0,1,5. For any other values of n, numerical methods are used to integrate the equation. The main difficulty in analyzing the Lane-Emden equation is the singularity behaviour occurring at x=0. Most researchers revisit the Lane-Emden equation repeatedly to test the viability of their new numerical methods against it. Accurate numerical results to seven decimal places are tabulated in Horedt [23, 43]. More accurate results including more decimal places have recently been reported in [24, 25]. In this work, we compare the results of the SLM approximate results to at least twelve decimal places for the dimensionless radius and mass of the polytropic stars with the results of [24, 25]. In Table 1, we present the SLM solution of the Lane-Emden equation (1.1) for the first zero α. The first zero is the smallest root α such that y(α)=0 and gives the dimensionless radius of the polytropic star [1]. The present results, at different iterations, are compared with the numerical results of [24, 25] to fifteen decimal places. The results presented in Table 1 we generated using N=40 collocation points for n=1,2 and N=70 for n=3,4. It can be seen from the table that the present SLM approach gives very accurate results which rapidly converge to the numerical result. Accuracy to fifteen decimal places is obtained after only four iterations for n=1, five iterations for n=2, six iterations for n=3, and seven iterations for n=4. We remark that the original SLM solution of the Lane-Emden equation reported in [21] gave results which were only accurate to seven decimal places even after increasing the number of iterations. This shows that the present modification of the original SLM approach results in a more robust method of solution. Similar results for nonintegral values of n are given in Table 2.In Table3, we give a comparison of the values of y'(1) which is used in the definition of the dimensionless mass of the polytropic star given by -α2y'(α) against the numerical results of [25]. Again, it can be seen that the present SLM results rapidly converge to the numerical results of [25].Figure1 shows the Lane-Emden equation solution curves for different indices obtained using eight iterations of the SLM solution series. These solution curves are qualitatively similar to those obtained using other numerical methods (see e.g., [22]).Figure 1
Graph of the SLM solutions of the Lane-Emden equation (1.1) for n=0.5,1,1.5,2,2.5,3,3.5, and 4 using i=8 iterations.The Frank-Kamenetskii equation (1.3) was solved using the SLM iteration scheme given by the matrix equation (2.32), and the present results were compared against the exact analytical results given by (2.21)-(2.22) for the case when k=1 and δ=1. The initial approximations required to start the SLM algorithm were generated using (2.28)-(2.29). For k=1 and δ, (2.29) generates two solutions for β0 which then result in two SLM solutions. Table 4 gives a comparison between the SLM approximate results to twenty decimal places against the exact analytical solutions. It can be seen from the table that for the first root, convergence to twenty decimal places is achieved after only three iterations. For the second root, convergence to twelve decimal places of accuracy is achieved after only five iterations and convergence to twenty decimal places is achieved after six iterations. From this observation, it can be seen that the present algorithm is very robust.In Table5, we give results for the maximum absolute errors for (1.3) when k=1 and δ=1 for different collocation points (grid points) at different iterations. The present results are compared with the recent B-Spline numerical results of [33]. It is worth noting that the B-Spline approach of [33] was meant to be an improvement on other numerical approaches that have previously been used to solve (1.3). For example, the results reported in [33] are more accurate than those reported in [34] where a three point finite difference approach was used to solve the same problem. Table 5 indicates that the present results are far more superior than the B-Spline results of [33] when the same number of grid points are used. We remark that results for only the first root were reported in [33] whereas our approach gives two solutions of (1.3). It can also be seen from Table 5 that the convergence of the method improves with an increase in the number of iterations.Table6 gives the results for the second root. Again, it can be seen from this table that the present method gives very accurate results which rapidly converge to the exact solution. Increasing N and the number of iterations improves the accuracy of the method.In Figure2, we plot the solution profile y(x) for different values of the Frank-Kamenetskii parameter δ after 6 iterations. It can be seen from the figure that there is good agreement between our approximate results and the exact analytical solutions.Figure 2
Comparison between the analytical solution (in diamonds) and the SLM approximate solutions for different values ofδ when k=1.


## 4. Conclusion

In this paper, we presented a new application of the successive linearization method (SLM) in solving Lane-Emden equations that model polytropic models of arbitrary indexn∈(0,5). The method is also tested on a singular boundary value Frank-Kamenetskii problem. The governing nonlinear equations were transformed to eigenvalue problems and subsequently solved using the SLM approach. A comparison was made between exact analytical solutions, numerical results from the literature and the present approximate solutions. The comparison indicates that the present method is robust and gives very accurate results and performs better that other numerical methods that have previously been applied to solve the same problems. The results presented in this paper can easily be extended to other initial and boundary value problems which are difficult to solve using other numerical methods.



---

*Source: 280702-2012-10-21.xml*