## 1. Introduction

In the most general form, Lane-Emden type equations are given as

where the prime denotes differentiation with respect to *x*, *f*(*x, y*) is a non-linear function, *g*(*x*) is a prescribed function and *x* = 0, Lane-Emden type equations are not trivial to solve. For this reason, the equations are normally used as benchmark equations for testing the effectiveness and robustness of new analytical and numerical methods of solution.

Analytical approaches that have recently been used in solving the Lane-Emden equations are mostly based on truncated series expansions. Examples include the Adomian decomposition method [1–3], differential transformation method [4, 5], Laplace transform [6, 7], homotopy analysis method [8–10], power series expansions [11–14] and variational iteration method [15–17]. Being power series based, the above methods have a small region of convergence and are not suitable for generating solutions in very large values of *x*. For this reason, most analytical approaches have only reported solutions of Lane-Emden type equations on small interval *x* axis. Despite this limitation, analytical approaches have been found to be desirable because they easily overcome the difficulty caused by the singularity at *x* = 0.

To overcome the limitations of analytical solution methods, several numerical approaches have been proposed for the solution of Lane-Emden type equations. Numerical methods based on spectral collocation have been found to be particularly effective. Collocation methods that have been reported recently for the solution of Lane-Emden type equations include the Bessel collocation method [18, 19], Jacobi-Gauss collocation method [20], Legendre Tau method [21], Sinc-collocation method [22], Chebyshev spectral methods [23], quasi-linearisation based Chebyshev pseudo-spectral method [24, 25] and a collocation method based on radial basis functions [26]. The discretisation scheme of collocation based methods is only implemented on interior nodes of the discretised domain. This property makes it possible for these collocation methods to overcome the difficulty of dealing with the singular point.

In this work, we present a multi-domain spectral collocation method for solving Lane-Emden equations. The method is based on the innovative idea of reducing the governing non-linear differential equations to a system of first-order equations which are solved iteratively using a Gauss-Seidel–like relaxation approach. The domain of the problem is divided into smaller non-overlapping sub-intervals on which the Chebyshev spectral collocation method is used to solve the iteration scheme. The continuity condition is used to advance the solution across neighbouring sub-intervals. The advantage of the approach is that it does not use Taylor-series based linearisation methods to simplify the non-linear differential equations. The method is free of errors associated with series truncation. The algorithm is also very easy to develop and yields very accurate results using only a few discretisation nodes. The accuracy of the method is validated against known results from the literature. The aim of the study is to explore the applicability of the multi-domain spectral collocation method to Lane-Emden type equations over semi-infinite domains. The results confirm that the method is suitable for solving all types of Lane-Emden equations.

## 2. A multi-domain pseudospectral relaxation method

In this section, we describe the development of the multi-domain collocation algorithm for the solution of Lane-Emden type equations. This algorithm addresses some limitations of standard collocation methods. Without loss of generality, we express the second-order generalised Lane-Emden Eq. (1) as a system of first-order ordinary differential equations (ODEs). If we let *h* – *y*_{0}, the initial value problem (1) transforms to

where *y′*(*x*) = *h*(*x*). The following iterative scheme for solving (2) and (3) is introduced:

Assuming that an initial approximation *i* + 1, Eqs. (4) and (5) form a pair of linear decoupled first-order differential equations. The solution procedure is discussed below.

Below, we describe the development of the multi-domain approach for solving the system of first-order Eqs. (4) and (5). The multi-domain technique approach assumes that the main interval can be decomposed into

The solution procedure assumes that the solution at each sub-interval

which interpolates

The choice of the interpolation function (7) with Lagrangre polynomial basis and grid points (8) enables the use of simple formulas for converting the continuous derivatives to discrete matrix-vector form at the collocation points *L*_{s}(x) is the characteristic Lagrange cardinal polynomial defined as

that obey the Kronecker delta equation:

The solution method seeks to solve each of the system of first-order ODEs independently in each *k*_{th} sub-interval

Therefore, in each sub-interval

where

Eqs. (12) and (13) cannot be solved exactly because of the non-linear function *f*(*x, y*). Accordingly, we employ the spectral collocation at the grid points

The derivatives at the collocation points are evaluated as

where *j*th and *s*th entry of the standard first derivative Chebyshev differentiation matrix of size

Including the relevant initial conditions in Eqs. (16) and (17) yields

where

Eqs. (18) and (19) can be expressed in matrix form as follows:

where the vectors

and *T* denotes the transpose of the vector. The *N* × *N* diagonal matrix ζ is given by

The approximate values of

## 3. Numerical experiments

In this section, we discuss the numerical experiments to be used to demonstrate the accuracy and general performance of the multi-domain pseudospectral relaxation method. In these numerical experiments, we have selected equations with known exact solutions, and to determine the accuracy of the method, we find the relative error. The relative error is defined as

where

### Example 1

We first consider the linear, homogeneous Lane-Emden equation, with variable coefficients:

which has the exact solution:

Eq. (30) has been solved by various researchers using different techniques such as the variational iteration method and the homotopy-pertubation method [16, 28, 29].

### Example 2

Secondly, the non-linear, homogeneous Lane-Emden equation, with variable coefficients is considered:

which has the exact solution:

Eq. (31) has been solved by [29, 30, 31] using the lie symmetry, the homotopy-pertubation method and the variational iteration method.

### Example 3

In this example, we consider the non-linear, variable coefficients, homogeneous Lane-Emden equation:

which has the exact solution:

Ramos [28] solved Eq. (32) using the variation iteration method, while Yildirim [16] solved the same equation using the linearisation method.

### Example 4

We consider the non-linear, homogeneous Lane-Emden equation:

which has the exact solution [32]:

This example was solved by Wazwaz [32] using the Adomian decomposition method.

### Example 5

We consider the non-linear, homogenous Lane-Emden equation, which represent an infinite circular cylinder in astrophysics:

which has the exact solution [33]:

This example was solved by Yin et al. [33] using the modified Laplace decomposition method.

### Example 6

Lastly, we consider the Lane-Emden equation with the general form:

Eq. (35) is the standard Lane-Emden equation that models the thermal behaviour of a spherical cloud of glass acting under the mutual attraction of its molecules and subject to classical laws of thermodynamics [34]. The values of *m* in the interval *m* = 0 and *m* = 1 and non-linear for values of

(36) |

Analytical solutions for

respectively. In this example, we consider the Lane-Emden Eq. (35) for *m* = 5. We therefore consider the equation:

with the exact solution:

## 4. Results and discussion

In this section, we discuss and present the results obtained using the proposed algorithm. We used the six examples in the previous section. The results were generated using MATLAB 2013. To validate the accuracy, computational time and general performance of the method, we computed relative errors and computational time of each of the numerical examples. The level of accuracy of the algorithm at a particular level is determined by the relative error *R _{k}* defined by

where *N* is the number of grid points, *x _{k}*. The graphs were all generated using

*N*= 4. For each numerical experiment, we present the relative error and the corresponding approximate solution for

*N*= 4 and

*N*= 6 in a tabular form. The central processing unit computational time is displayed. Graphs showing an excellent agreement between the analytical and approximate solutions are presented for each numerical experiment. These graphs validate the accuracy of the method. Error graphs showing the distribution of the relative errors are also presented. These error graphs are in excellent agreement with the results presented in the tables for all the numerical experiments used in this chapter.

**Table 1** shows the exact, approximate solution and the relative error for Eq. (30). For *N* = 4, the multi-domain spectral relaxation method gives a relative error of approximately 10^{−12} and for *N* = 6, the relative error is on average 10^{−13}. Increasing the number of grid points results in a more accurate solution. The results obtained from the multi-domain spectral relaxation method are remarkable since few grid points give accurate results in a large domain. Using a few grid points ensures that the numerical method converges within few seconds.

**Figure 1** shows the relative error displayed in **Table 1** for *N* = 4. The results in **Figure 1** are in excellent agreement with those in **Table 1**. **Figure 2** shows the analytical and approximate solutions. Since the approximate solution is superimposed on the exact solutions, this implies that the multi-domain pseudospectral relaxation method converged to the exact solution over the domain **Table 1** and **Figures 1** and **2** validate the accuracy and computational efficiency of the multi-domain pseudospectral relaxation method for Eq. (30).

The results obtained from approximating the solution to Eq. (31) using the multi-domain pseudospectral relaxation method are shown in **Table 1** and **Figures 3** and **4**. In **Table 2**, we display the exact solution, approximate solution and the relative error of Eq. (31) in Example 2. For *N* = 4, the multi-domain spectral relaxation method gives a relative error of approximately 10^{−11}. For *N* = 4 the relative error is approximately 10^{−13}. We observe that increasing the number of grid points decreases the relative error. The multi-domain pseudospectral relaxation method uses a few grid points to achieve accurate results in the domain *N* = 6 grid points ensured that the numerical method converged to an error of 10^{−13} within a fraction of a second.

**Figure 3** shows the relative error displayed in **Table 2** for *N* = 4. The results in **Figure 3** are in excellent agreement with those in **Table 2**. **Figure 4** shows the analytical and approximate solutions of Eq. (31). The approximate solution superimposed on the exact solutions shows that the multi-domain pseudospectral relaxation method converged to the exact solution over the domain **Table 2** and **Figures 3** and **4** validates the accuracy and computational efficiency of the multi-domain pseudospectral relaxation method for Eq. (31).

The results obtained from approximating the solution to Eq. (32) are given in **Table 3** and **Figures 5** and **6**. **Table 3** shows the exact solution, the approximate solution and the relative error of Eq. (32). For *N* = 4, the multi-domain spectral relaxation method gives a relative error of approximately 10^{−12}, while for *N* = 6, the relative error is approximately 10^{−13}. We observe that increasing the number of grid points decreases the relative error and hence increases the accuracy of the method. This pseudospectral method uses a few grid points to achieve accurate results in the domain *N* = 6 grid points ensured that the numerical method converged to an error of 10^{−13} within few seconds.

The relative error shown in **Table 3** is displayed in **Figure 5**. The results in **Figure 5** are in excellent agreement with those in **Table 3**. **Figure 6** shows the analytical and approximate solutions of Eq. (32). The approximate solution being superimposed on the exact solutions implies that the multi-domain pseudospectral relaxation method converged to the exact solution over the domain **Table 3** and **Figures 5** and **6** validates the accuracy and computational efficiency of the multi-domain pseudospectral relaxation method for Eq. (32).

The results obtained from approximating the solution to Eq. (33) are given in **Table 4** and **Figures 7** and **8**. **Table 4** shows the exact solution, the approximate solution and the relative error of Eq. (33). For *N* = 4, the multi-domain spectral relaxation method gives a relative error of approximately 10^{−13}. For *N* = 6, the relative error is also approximately 10^{−13}. Increasing the number of grid points decreases the relative error. Thus, a maximum of *N* = 6 grid points ensures convergence of the method. *N* = 6 grid points ensured that the numerical method converged to an error of 10^{−13} in a fraction of seconds.

**Figure 7** shows the relative error displayed in **Table 4** for *N* = 4. The results in **Figure 7** are in excellent agreement with those in **Table 4**. **Figure 8** shows the analytical and approximate solutions. Since the approximate solution is superimposed on the exact solutions, this implies that the multi-domain pseudospectral relaxation method converged to the exact solution over the domain **Table 4** and **Figures 7** and **8** validate the accuracy and computational efficiency of the multi-domain pseudospectral relaxation method for Eq. (33).

The results obtained from approximating the solution to Eq. (34) are given in **Table 5** and **Figures 9** and **10**. **Table 5** shows the exact solution, the approximate solution and the relative error of Eq. (34). The multi-domain pseudospectral relaxation method gives a relative error of approximately 10^{−14}. For *N* = 6, the relative error is approximately 10^{−14}. Increasing the number of grid points decreases the relative error and thus implying that the numerical method converged to the exact solution. The pseudospectral collocation method uses a few grid points to achieve accurate results in the domain ^{−14} in a fraction of a second.

**Figure 9** shows the relative error displayed in **Table 5** for *N* = 4. The results in **Figure 9** are in excellent agreement with those in **Table 5**. **Figure 10** shows the analytical and approximate solutions. Since the approximate solution is superimposed on the exact solutions, this implies that the multi-domain pseudospectral relaxation method converged to the exact solution over the domain **Table 5** and **Figures 9** and **10** validate the accuracy and computational efficiency of the multi-domain pseudospectral relaxation method for Eq. (34).

The results obtained from approximating the solution to Eq. (38) are given in **Table 6** and **Figures 11**, **12** and **13**. **Table 6** shows values obtained for the exact and approximate solution together with the respective relative error values of Eq. (38). For *N* = 4 the multi-domain pseudospectral collocation method gives a relative error of approximately ^{−13} for *N* = 6. An increase in the number of grid points results in a decrease in the relative error. This implies that the numerical method converges to the exact solution of Eq. (38). The multi-domain pseudospectral relaxation method used a few grid points to achieve accurate results in the domain *N* = 6 grid points ensured that the numerical method converged to an error of 10^{−13} within a few seconds as shown in **Table 6**.

**Figure 11** shows the plot for different values of *m*. The results are in good agreement with those obtained by [25, 26]. **Figure 12** displays the relative error graph of Eq. (38). The results in **Figure 12** are in excellent agreement with those obtained in **Table 6**. The comparison between the exact solution and approximate solution of Eq. (38) is shown in **Figure 13**. The superimposition of the approximate solution on the exact solution implies that the multi-domain pseudospectral relaxation method converged to the exact solution over the domain **Table 6** and **Figures 11**, **12** and **13** give a validation of the accuracy and computational efficiency of the multi-domain pseudospectral relaxation method for Eq. (38).

## 5. Conclusion

In this work, we presented a multi-domain spectral collocation method for solving Lane-Emden equations. This numerical method was used to solve six Lane-Emden equations. The results obtained were remarkable in the sense that using few grid points, we were able to achieve accurate results. We were able to compute the results using minimal computational time. The approximate solutions were in excellent agreement with the exact solutions in all the numerical experiments. We presented error graphs, approximate and exact solution graphs to show accuracy of the method. Tables showing relative errors were also generated to show accuracy and computational efficiency of the numerical method presented.

This approach is useful for solving other nonlinear, singular initial value problems. This approach is an alternative to the already existing list of numerical methods that can be used to solve such equations. This numerical approach can be extended to solve time-dependent Lane-Emden equations and other singular value type of equations.