`Mathematical Problems in EngineeringVolume 2012 (2012), Article ID 280702, 14 pageshttp://dx.doi.org/10.1155/2012/280702`
Research Article

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

1School of Mathematical Sciences, University of KwaZulu-Natal, Private Bag X01, Pietermaritzburg, Scottsville 3209, South Africa
2Department of Mathematics, University of Venda, Privatre Bag X5050, Thohoyandou 0950, South Africa

Received 9 March 2012; Revised 16 May 2012; Accepted 29 May 2012

Copyright © 2012 S. S. Motsa and S. Shateyi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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 with initial conditions where 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 , exact analytical solutions were obtained by Chandrasekhar [1]. For all other values of , 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 [913], 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 . For , 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 is defined as the smallest positive value for which . 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 because their region of convergence is usually less than . 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 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 and . Analytical treatment of (1.3) leading to series solutions was carried out in [2830] 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 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 [3639] 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 , where is the first zero of . 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 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., [3639]). 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 where are unknown functions that are obtained by iteratively solving the linearized version of the governing equations assuming that are known from previous iterations. The algorithm starts with an initial approximation 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 () 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 An initial approximation for that was found to give accurate results for all was This value of was obtained by substituting the initial approximation given by (2.4) in (2.2) and solving the resulting equation at , the middle of the domain for .

Substituting the expansions (2.3) and (2.5) into (2.2) and neglecting nonlinear terms in and give subject to where where the dots denote differentiation with respect to the scaled variable . To solve (2.7), we use the Chebyshev collocation method in which case the functions 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 We use the Gauss-Lobatto collocation points [40, 41] to define the Chebyshev nodes in [−1,1] as where is the number of collocation points. The derivatives of at the collocation points can be written as where , is transpose, and is the Chebyshev derivative matrix of size whose entries are defined [40, 41] as with Applying the Chebyshev spectral collocation method in (2.7)-(2.8) gives subject to The matrices in (2.15) are defined as We note that the system (2.15) consists of unknowns in and an additional unknown . To solve for all 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 st rows of the coefficient matrix and vectors and deleting the first and ()st columns of the coefficient matrix. Thus, starting from the initial approximations and , the solutions and can be determined from

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

The Frank-Kamenetskii equation admits exact analytical solutions when (see e.g., [27, 42]). When , it has two solutions. The problem has no solution for and when , it has a unique solution. The solutions are given by [42] where the constants and are given by

To solve (1.3), it is convenient to introduce the boundary condition at as where is a constant. The governing equation (1.3) is then solved subject to the Dirichlet boundary conditions and using the SLM approach as described in the previous subsection using the Neumann boundary condition as an extra condition for determining the unknown . To this end, we look for a solution of the form where and 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 and give subject to where A suitable initial guess that satisfies the boundary conditions (1.4) and (2.23) is where 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 . The solutions for 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 where The equation system can be written as the matrix equation where the boundary conditions have been imposed on the first and last rows of and and the derivative boundary condition has been added in the row. Thus, starting from and given by (2.28) and (2.29), the subsequent solutions for 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 16 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 against numerical results given by Boyd [25] and Seidov [24].
Table 2: Comparison of the first zero () of the Lane-Emden equation at different iterations for different values of against numerical results given by Boyd [25].
Table 3: Comparison of the value of , for the Lane-Emden equation at different iterations for different values of against numerical results given by Boyd [25].
Table 4: SLM results for corresponding to the two roots of (1.3) when , and .
Table 5: Maximum absolute errors for solution of (1.3) (1st root).
Table 6: Maximum absolute errors for solution of (1.3) (2nd root).

In the case of the Lane-Emden equation (1.1), exact analytical solutions are only available for the case when . For any other values of , numerical methods are used to integrate the equation. The main difficulty in analyzing the Lane-Emden equation is the singularity behaviour occurring at . 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 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 collocation points for and for . 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 , five iterations for , six iterations for , and seven iterations for . 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 are given in Table 2.

In Table 3, we give a comparison of the values of which is used in the definition of the dimensionless mass of the polytropic star given by 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]).

Figure 1: Graph of the SLM solutions of the Lane-Emden equation (1.1) for , and 4 using 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 and . The initial approximations required to start the SLM algorithm were generated using (2.28)-(2.29). For and , (2.29) generates two solutions for 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 and 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 and the number of iterations improves the accuracy of the method.

In Figure 2, we plot the solution profile 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 .

#### 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 index . 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.

#### References

1. S. Chandrasekhar, An Introduction to the Study of Stellar Structure, Dover, New York, NY, USA, 1958.
2. H. T. Davis, Introduction to Nonlinear Differential and Integral Equations, Dover, New York, NY, USA, 1962.
3. A.-M. Wazwaz, “A new method for solving singular initial value problems in the second-order ordinary differential equations,” Applied Mathematics and Computation, vol. 128, no. 1, pp. 45–57, 2002.
4. A.-M. Wazwaz, “Adomian decomposition method for a reliable treatment of the Emden-Fowler equation,” Applied Mathematics and Computation, vol. 161, no. 2, pp. 543–560, 2005.
5. V. S. Ertürk, “Differential transformation method for solving differential equations of Lane-Emden type,” Mathematical & Computational Applications, vol. 12, no. 3, pp. 135–139, 2007.
6. A. Yildirim and T. Öziş, “Solutions of singular IVPs of Lane-Emden type by homotopy perturbation method,” Physics Letters, Section A, vol. 369, no. 1-2, pp. 70–76, 2007.
7. A. S. Bataineh, M. S. M. Noorani, and I. Hashim, “Homotopy analysis method for singular IVPs of Emden-Fowler type,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 4, pp. 1121–1131, 2009.
8. S. Liao, “A new analytic algorithm of Lane-Emden type equations,” Applied Mathematics and Computation, vol. 142, no. 1, pp. 1–16, 2003.
9. A. Aslanov, “Determination of convergence intervals of the series solutions of Emden-Fowler equations using polytropes and isothermal spheres,” Physics Letters. A, vol. 372, no. 20, pp. 3555–3561, 2008.
10. C. Hunter, “Series solutions for polytropes and the isothermal sphere,” Monthly Notices of the Royal Astronomical Society, vol. 328, no. 3, pp. 839–847, 2001.
11. C. Mohan and A. R. Al-Bayaty, “Power-series solutions of the Lane-Emden equation,” Astrophysics and Space Science, vol. 73, no. 1, pp. 227–239, 1980.
12. I. W. Roxburgh and L. M. Stockman, “Power series solutions of the polytrope equations,” Monthly Notices of the Royal Astronomical Society, vol. 303, no. 3, pp. 466–470, 1999.
13. Z. F. Seidov, “The power series as solution of the Lane-Emden equation with index two,” Doklady Akademii Nauk SSSR, vol. 35, pp. 21–24, 1979.
14. J. H. He, “Variational approach to the Lane-Emden equation,” Applied Mathematics and Computation, vol. 143, no. 2-3, pp. 539–541, 2003.
15. A. Yıldırım and T. Öziş, “Solutions of singular IVPs of Lane-Emden type by the variational iteration method,” Nonlinear Analysis, vol. 70, no. 6, pp. 2480–2484, 2009.
16. S. S. Motsa and S. Shateyi, “New analytic solution to the Lane-Emden equation of index 2,” Mathematical Problems in Engineering, vol. 2012, Article ID 614796, 19 pages, 2012.
17. S. Liao, Beyond Perturbation: Introduction to Homotopy Analysis Method, vol. 2 of CRC Series: Modern Mechanics and Mathematics, Chapman & Hall/CRC, Boca Raton, Fla, USA, 2004.
18. K. Parand and M. Razzaghi, “Rational legendre approximation for solving some physical problems on semi-infinite intervals,” Physica Scripta, vol. 69, no. 5, pp. 353–357, 2004.
19. K. Parand and A. Pirkhedri, “Sinc-Collocation method for solving astrophysics equations,” New Astronomy, vol. 15, no. 6, pp. 533–537, 2010.
20. K. Parand, A. R. Rezaei, and A. Taghavi, “Lagrangian method for solving LaneEmden type equation arising in astrophysics on semi-infinite domains,” Acta Astronautica, vol. 67, no. 7-8, pp. 673–680, 2010.
21. S. S. Motsa and P. Sibanda, “A new algorithm for solving singular IVPs of Lane-Emden type,” in Proceedings of the 4th International Conference on Applied Mathematics, Simulation, Modelling (ASM'10), pp. 176–180, Corfu Island, Greece, July 2010.
22. K. Parand, S. Abbasbandy, S. Kazem, and A. R. Rezaei, “An improved numerical method for a class of astrophysics problems based on radial basis functions,” Physica Scripta, vol. 83, no. 1, Article ID 015011, 2011.
23. G. P. Horedt, “Seven-digit tables of Lane-Emden functions,” Astrophysics and Space Science, vol. 126, no. 2, pp. 357–408, 1986.
24. Z. F. Seidov, “Lane-Emden equation: perturbation method,” In press, http://arxiv.org/abs/astro-ph/0402130.
25. J. P. Boyd, “Chebyshev spectral methods and the Lane-Emden problem,” Numerical Mathematics. Theory, Methods and Applications, vol. 4, no. 2, pp. 142–157, 2011.
26. D. A. Frank-Kamenetskii, Diffusion and Heat Transfer in Chemical Kinetics, Plenum Press, New York, NY, USA, 1969.
27. M. Kubíček and V. Hlavačék, Numerical Solution of Nonlinear Boundary Value Problems with Applications, Prentice-Hall, Englewood Cliffs, NJ, USA, 1983.
28. S. A. Khuri and A. Sayfy, “A novel approach for the solution of a class of singular boundary value problems arising in physiology,” Mathematical and Computer Modelling, vol. 52, no. 3-4, pp. 626–636, 2010.
29. M. Kumar and N. Singh, “Modified Adomian Decomposition Method and computer implementation for solving singular boundary value problems arising in various physical problems,” Computers and Chemical Engineering, vol. 34, no. 11, pp. 1750–1760, 2010.
30. A.-M. Wazwaz, “A new algorithm for solving differential equations of Lane-Emden type,” Applied Mathematics and Computation, vol. 118, no. 2-3, pp. 287–310, 2001.
31. E. Momoniat and C. Harley, “Approximate implicit solution of a Lane-Emden equation,” New Astronomy, vol. 11, no. 7, pp. 520–526, 2006.
32. A. Asaithambi, “A shooting method for nonlinear heat transfer using automatic differentiation,” Applied Mathematics and Computation, vol. 180, no. 1, pp. 264–269, 2006.
33. H. Çaǧlar, N. Çaǧlar, and M. Özer, “B-spline solution of non-linear singular boundary value problems arising in physiology,” Chaos, Solitons and Fractals, vol. 39, no. 3, pp. 1232–1237, 2009.
34. S. R. K. Iyengar and P. Jain, “Spline finite difference methods for singular two point boundary value problems,” Numerische Mathematik, vol. 50, no. 3, pp. 363–376, 1987.
35. M. Kumar and Y. Gupta, “Methods for solving singular boundary value problems using splines: a review,” Journal of Applied Mathematics and Computing, vol. 32, no. 1, pp. 265–278, 2010.
36. F. G. Awad, P. Sibanda, S. S. Motsa, and O. D. Makinde, “Convection from an inverted cone in a porous medium with cross-diffusion effects,” Computers & Mathematics with Applications, vol. 61, no. 5, pp. 1431–1441, 2011.
37. Z. G. Makukula, P. Sibanda, and S. S. Motsa, “A novel numerical technique for two-dimensional laminar flow between two moving porous walls,” Mathematical Problems in Engineering, vol. 2010, Article ID 528956, 15 pages, 2010.
38. S. S. Motsa and S. Shateyi, “A new approach for the solution of three-dimensional magnetohydrodynamic rotating flow over a shrinking sheet,” Mathematical Problems in Engineering, vol. 2010, Article ID 586340, 2010.
39. S. S. Motsa, “New algorithm for solving non-linear BVPs in heat transfer,” International Journal of Modeling, Simulation & Scientific Computing, vol. 2, no. 3, pp. 355–373, 2011.
40. C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer Series in Computational Physics, Springer, New York, NY, USA, 1988.
41. L. N. Trefethen, Spectral Methods in MATLAB, vol. 10, SIAM, Philadelphia, Pa, USa, 2000.
42. C. Harley and E. Momoniat, “Efficient boundary value problem solution for a Lane-Emden equation,” Mathematical & Computational Applications, vol. 15, no. 4, pp. 613–620, 2010.
43. G. P. Horedt, Polytropes Applications in Astrophysics and Related Fields, Kluwer Academic Publishers, Dordrecht, The Netherlands, 2004.