A Successive Linearization Method Approach to Solve Lane-Emden Type of Equations

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.


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 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 of y x is defined as the smallest positive value x 0 for which y x 0 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 x 0 because their region of convergence is usually less than x 0 .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 which is solved subject to the boundary conditions 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.

Numerical Solution
In this section, we present the method used to solve the governing equations 1.1 -1.4 .

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 x αt.
2.1 Substituting 2.1 into 1.1 and simplifying gives 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 Mathematical Problems in Engineering of linear differential equations which are then solved using analytical or numerical methods wherever possible.In applying the SLM on 2.2 , we set where y i are unknown functions that are obtained by iteratively solving the linearized version of the governing equations assuming that y i 0 ≤ m ≤ i − 1 are known from previous iterations.The algorithm starts with an initial approximation y 0 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 We observe that 2.2 is a second-order differential equation but has three boundary conditions.The third boundary condition y t 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 An initial approximation for α that was found to give accurate results for all n ∈ 0, 5 was

2.6
This value of α 0 was obtained by substituting the initial approximation y 0 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 y i and α i give subject to where 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 y i 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
We use the Gauss-Lobatto collocation points 40, 41 to define the Chebyshev nodes in −1,1 as where N 1 is the number of collocation points.The derivatives of y i at the collocation points can be written as where Y i y i z 0 , y i z 1 , . . ., y i z N−1 , y i z N T , T is transpose, and D is the Chebyshev derivative matrix of size N 1 × N 1 whose entries are defined 40, 41 as

Mathematical Problems in Engineering
Applying the Chebyshev spectral collocation method in 2.7 -2.8 gives subject to The matrices in 2.15 are defined as

2.18
We note that the system 2.15 consists of N 1 unknowns in Y i 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 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 Y 0 and α 0 , the solutions Y i and α i i 1, 2, 3, . . .can be determined from Mathematical Problems in Engineering 7

Solution of the Frank-Kamenetskii Equation
The Frank-Kamenetskii equation admits exact analytical solutions when k 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 where the constants c 1 and c 2 are given by To solve 1.3 , it is convenient to introduce the boundary condition at x 0 as 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 where y i 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 y i and β i give subject to where

2.27
A suitable initial guess that satisfies the boundary conditions 1.4 and 2.23 is where β 0 is the initial approximation for β which is given by 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 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.31
The equation system 2.30 can be written as the matrix equation where the boundary conditions have been imposed on the first and last rows of C i−1 and S i−1 and the derivative boundary condition has been added in the N 2 row.Thus, starting from y 0 and β 0 given by 2.28 and 2.29 , the subsequent solutions for y i , β i i 1, 2, 3, . . .can be obtained by solving the matrix system 2.32 .

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.
In 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 Table 3, 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 −α 2 y α against the numerical results of 25 .Again, it can be seen that the present SLM results rapidly converge to the numerical results of 25 .
Figure 1 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 .
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 Table 5, 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.Table 6 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 Figure 2, 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.

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 index n ∈ 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.

Table 1 :
Comparison of the first zero α of the Lane-Emden equation at different iterations for different values of n against numerical results given byBoyd 25 and Seidov 24 .

Table 2 :
Comparison of the first zero α of the Lane-Emden equation at different iterations for different values of n against numerical results given byBoyd 25 .

Table 3 :
Comparison of the value of dy/dt 1 , for the Lane-Emden equation at different iterations for different values of n against numerical results given byBoyd 25 .

Table 4 :
SLM results for y 0 corresponding to the two roots of 1.3 when k 1, δ 1, and N 50.

Table 5 :
Maximum absolute errors for solution of 1.3 1st root .

Table 6 :
Maximum absolute errors for solution of 1.3 2nd root .