A Multi-Domain Spectral Collocation Approach for Solving Lane-Emden Type Equations

Abstract

In this work, we explore the application of a novel multi-domain spectral collocation method for solving general non-linear singular initial value differential equations of the Lane-Emden type. The proposed solution approach is a simple iterative approach that does not employ linearisation of the differential equations. Spectral collocation is used to discretise the iterative scheme to form matrix equations that are solved over a sequence of non-overlapping sub-intervals of the domain. Continuity conditions are used to advance the solution across the non-overlapping sub-intervals. Different Lane-Emden equations that have been reported in the literature have been used for numerical experimentation. The results indicate that the method is very effective in solving Lane-Emden type equations. Computational error analysis is presented to demonstrate the fast convergence and high accuracy of the method of solution.

1. Introduction

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

$y^{″}+\frac{\gamma }{x}{y}^{\prime }+f(x,y)=g(x),x\in [0,\infty ),y(0)={\alpha }_0,y^{\prime }(0)={\beta }_0$E1

where the prime denotes differentiation with respect to x, f(x, y) is a non-linear function, g(x) is a prescribed function and $\gamma$, ${\alpha }_0$, ${\beta }_0$ are known constants. In recent years, problems described by this class of differential equations have been widely investigated by many researchers because of their applications in astronomy, mathematical biology, mathematical physics, non-Newtonian fluid mechanics, and other areas of science and engineering. From a solution method viewpoint, it has been observed that, owing to the singularity at 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 $[\mathrm{0,1}]$ on 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₀, the initial value problem (1) transforms to

$y^{\prime }(x)=h(x),y(0)={\alpha }_0,$E2

$h^{\prime }(x)+\frac{\gamma }{x}h+f(x,y)=g(x),h(0)={\beta }_0,$E3

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

${y^{\prime }}_{i+1}(x)=h_i(x),y_{i+1}(0)={\alpha }_0,$E4

${h^{\prime }}_{i+1}(x)+\frac{\gamma }{x}{h}_{i+1}(x)+f(x,y_{i+1})=g(x),h_{i+1}(0)={\beta }_0\mathrm{.}$E5

Assuming that an initial approximation $h_0(x)$ is given, Eq. (4) can be solved for $y_{i+1}(x)$ and the solution can be used immediately in Eq. (5) which is, in turn, solved for $h_{i+1}$. Thus, at each iteration level 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 $p$ non-overlapping sub-intervals. Let $x\in \text{Λ}$, where $\text{Λ}=[a,b]$ is the interval where the solution of Eq. (1) exists. The sub-intervals are defined as

${\text{Λ}}_k=(x_{k-1},x_k),k=\mathrm{1,2,}\cdots ,p,\text{with},a=x_0<x_1<x_2<\cdots <x_p=b\mathrm{.}$E6

The solution procedure assumes that the solution at each sub-interval ${\text{Λ}}_k$, denoted by $y^{(k)}(x)$, can be approximated by a Lagrange interpolation polynomial of the form

$y^{(k)}(x)\sim {\displaystyle {\displaystyle \sum _{s=0}^N}}{y}^{(k)}(x_s)L_s(x),$E7

which interpolates $y^{(k)}(x)$ at selected points, chosen to be the Chebyshev-Gauss-Lobatto, defined by

${\left\{{\tau }_s\right\}}_{s=0}^N={\left\{\cos \left(\frac{\pi s}{N}\right)\right\}}_{s=0}^N\mathrm{.}$E8

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 ${\tau }_s$. The function L_s(x) is the characteristic Lagrange cardinal polynomial defined as

$L_s(x)={\displaystyle {\displaystyle \prod _{\begin{array}{c}s=0\\ s\ne k\end{array}}^M}}\frac{x-x_k}{x_s-x_k},$E9

that obey the Kronecker delta equation:

$L_s(x_k)={\delta }_{sk}=\{\begin{array}{cc}0& \text{if}s\ne k\\ 1& \text{if}s=k\mathrm{.}\end{array}\mathrm{.}$E10

The solution method seeks to solve each of the system of first-order ODEs independently in each k_th sub-interval ${\text{Λ}}_k$ using the solutions obtained in the preceding interval ${\text{Λ}}_{k-1}$ as initial conditions. In the first interval the solution is computed on $[x_0,x_1]$ and is labelled $y^{(1)}(x).$ The value of the solution at the last node $y^{(1)}(x_1)$ is used as an initial condition when seeking a solution in the second sub-interval ${\text{Λ}}_2$. This procedure is repeated in each ${\text{Λ}}_k$ interval with the solutions marched across each interval using the following continuity condition:

$y^{(k)}(x_{k-1})=y^{(k-1)}(x_{k-1}),h^{(k)}(x_{k-1})=h^{(k-1)}(x_{k-1}).$E11

Therefore, in each sub-interval ${\text{Λ}}_k$, the following system of ODEs is solved:

$\frac{\text{d}{y}_{i+1}^{(k)}(x)}{\text{d}x}=h_i^{(k)}(x),y_{i+1}^{(k)}(x_{k-1})={\alpha }_{k-1}$E12

$\frac{\text{d}{h}_{i+1}^{(k)}(x)}{\text{d}x}+\frac{\gamma }{x}{h}_{i+1}^{(k)}(x)=g(x)-f(x,y_{i+1}^{(k)}),h_{i+1}^{(k)}(x_{k-1})={\beta }_{k-1},$E13

where

${\alpha }_{k-1}=y^{(k-1)}(x_{k-1}),{\beta }_{k-1}=h^{(k-1)}(x_{k-1}),k=\mathrm{1,2,}\dots ,p\mathrm{.}$

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 $x_j$ for $j=\mathrm{0,1,2,}\cdots ,N$ for each sub-interval ${\text{Λ}}_k$. Before the spectral method is applied, each sub-interval $[x_{k-1},x_k]$ is transformed to $[-\mathrm{1,1}]$ using the linear transformation:

$x=\frac{x_k-x_{k-1}}{2}\tau +\frac{x_k+x_{k-1}}{2},\tau \in [-\mathrm{1,1}].$E14

The derivatives at the collocation points are evaluated as

${\frac{\text{d}{y}^{(k)}}{\text{d}x}|}_{x=x_j}={\displaystyle {\displaystyle \sum _{s=0}^N}}{y}^{(k)}(x_s)\frac{\text{d}{L}_s(x_j)}{\text{d}x}=\frac{2}{\text{Δ}{x}_k}{\displaystyle {\displaystyle \sum _{s=0}^N}}{y}^{(k)}({\tau }_s)\frac{\text{d}{L}_s({\tau }_j)}{\text{d}\tau }={\displaystyle {\displaystyle \sum _{s=0}^N}}{D}_{js}{y}^{(k)}({\tau }_s),$E15

where $\text{Δ}{x}_k=x_k-x_{k-1}$, $D_{js}=\frac{2}{\text{Δ}{x}_k}{D}_{js}$ with $D_{js}=\frac{d{L}_s({\tau }_j)}{d\tau }$ being the jth and sth entry of the standard first derivative Chebyshev differentiation matrix of size $(N+1)\times (N+1)$ as defined in [27]. Evaluating Eqs. (12) and (13) at the grid points ${\tau }_j$ for $j=\mathrm{0,1,2,}\cdots ,N$ gives

${\displaystyle {\displaystyle \sum _{s=0}^N}}{D}_{js}{y}_{i+1}^{(k)}({\tau }_s)=h_i^{(k)}({\tau }_j),y_{i+1}^{(k)}({\tau }_N)={\alpha }_{k-1},$E16

${\displaystyle {\displaystyle \sum _{s=0}^N}}{D}_{js}{h}_{i+1}^{(k)}({\tau }_s)+\frac{\gamma }{x_j}{h}_{i+1}^{(k)}({\tau }_j)=g(x_j)-f(x_s,y_{i+1}^{(k)}({\tau }_j)),h_{i+1}^{(k)}({\tau }_N)={\beta }_{k-1}\mathrm{.}$E17

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

${\displaystyle {\displaystyle \sum _{s=0}^{N-1}}}{D}_{js}{y}_{i+1}^{(k)}({\tau }_s)=V_i^{(k)}({\tau }_j),$E18

${\displaystyle {\displaystyle \sum _{s=0}^{N-1}}}{D}_{js}{h}_{i+1}^{(k)}({\tau }_s)+\frac{\gamma }{x_j}{h}_{i+1}^{(k)}({\tau }_j)=W_i^{(k)}({\tau }_j).$E19

where

$V_i^{(k)}({\tau }_j)=h_i^{(k)}({\tau }_j)-D_{jN}{y}_{i+1}^{(k)}({\tau }_N),$E20

$W_i^{(k)}({\tau }_j)=g(x_j)-f(x_j,y_{i+1}^{(k)}({\tau }_j))-D_{jN}{h}_{i+1}^{(k)}({\tau }_N),$E21

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

$D{Y}_{i+1}^{(k)}=V_i^{(k)},$E22

$\left(D+\zeta \right)H_{i+1}^{(k)}=W_i^{(k)},$E23

where the vectors $Y^{(k)},H^{(k)},V^{(k)}$ and $W^{(k)}$ are respectively given by

$Y^{(k)}=[y^{(k)}({\tau }_0),y^{(k)}({\tau }_1),y^{(k)}({\tau }_2),\cdots ,y^{(k)}({\tau }_{N-1}{)]}^T,$E24

$H^{(k)}=[h^{(k)}({\tau }_0),h^{(k)}({\tau }_1),h^{(k)}({\tau }_2),\cdots ,h^{(k)}({\tau }_{N-1}{)]}^T,$E25

$V^{(k)}=[v^{(k)}({\tau }_0),v^{(k)}({\tau }_1),v^{(k)}({\tau }_2),\cdots ,v^{(k)}({\tau }_{N-1}{)]}^T,$E26

$W^{(k)}=[w^{(k)}({\tau }_0),w^{(k)}({\tau }_1),w^{(k)}({\tau }_2),\cdots ,w^{(k)}({\tau }_{N-1}{)]}^T,$E27

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

$\zeta =\left[\begin{array}{cccc}\frac{\gamma }{x_0}& & & \\ & \frac{\gamma }{x_1}& & \\ & & \ddots & \\ & & & \frac{\gamma }{x_{N-1}}\end{array}\right]\mathrm{.}$E28

The approximate values of $y^{(k)}(x)$ and $\frac{d{y}^{(k)}(x)}{dx}$ are obtained by iteratively solving Eqs. (22) and (23), respectively, for $i=\mathrm{0,1,2,}\dots$, starting from suitable initial approximation.

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

$E_j=\frac{|y_e(x_j)-y_a(x_j)|}{|y_e(x_j)|}$E29

where $E_j$ is the relative error at a grid point $x_j$, $y_e(x_j)$ and $y_a(x_j)$ are the exact and approximate solutions at a grid point $x_j$, respectively.

Example 1

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

$y^{″}+\frac{2}{x}{y}^{\prime }-2(2x^2+3)y=0,x>0,\text{subject}\text{to}y(0)=1\text{and}{y}^{\prime }(0)=0,$E30

which has the exact solution:

$y(x)=e^{x^2}\mathrm{.}$

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:

$y^{″}+\frac{2}{x}{y}^{\prime }+4(2e^y+e^{y/2})=0,x>0,\text{subject}\text{to}y(0)=0\text{and}{y}^{\prime }(0)=0$E31

which has the exact solution:

$y(x)=-2\ln (1+x^2).$

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:

$y^{″}+\frac{2}{x}{y}^{\prime }-6y=4y\ln (y),x>0,\text{subject}\text{to}y(0)=1\text{and}{y}^{\prime }(0)=0$E32

which has the exact solution:

$y(x)=e^{x^2}\mathrm{.}$

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:

$y^{″}+\frac{6}{x}{y}^{\prime }+2y(7+\ln (y^2))=0,x>0,\text{subject}\text{to}y(0)=1\text{and}{y}^{\prime }(0)=0$E33

which has the exact solution [32]:

$y(x)=e^{-x^2}\mathrm{.}$

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:

$y^{″}+\frac{2}{x}{y}^{\prime }+y^3-6-x^6=0,x>0,\text{subject}\text{to}y(0)=0\text{and}{y}^{\prime }(0)=0,$E34

which has the exact solution [33]:

$y(x)=x^2\mathrm{.}$

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:

$y^{″}+\frac{2}{x}{y}^{\prime }+y^m=0,x>0,\text{subject}\text{to}y(0)=1\text{and}{y}^{\prime }(0)=0.$E35

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 $[\mathrm{0,5}]$ are most physically interesting to study. The equation is linear when m = 0 and m = 1 and non-linear for values of $m>1$. Wazwaz [35] gave the general solution of Eq. (35) in series form as

$\begin{array}{c}y(x)=1-\frac{1}{6}{x}^2+\frac{m}{120}{x}^4-\frac{m(8m-5)}{3.7!}{x}^6+\frac{m(70-183m+122m^2)}{9.9!}{x}^8\\ +\frac{m(3150-1080m+12642m^2-5032m^3)}{45.11!}{x}^{10}+\dots \mathrm{.}\end{array}$E36

Analytical solutions for $m=\mathrm{0,1}$ and $m=5$ are given as [35]

$y(x)=1-\frac{1}{3!}{x}^2,y(x)=\frac{\sin x}{x},\text{and}y(x)={\left(1+\frac{x^2}{3}\right)}^{-\frac{1}{2}},$E37

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

$y^{″}+\frac{2}{x}{y}^{\prime }+y^5=0,x>0,\text{subject}\text{to}y(0)=1\text{and}{y}^{\prime }(0)=0$E38

with the exact solution:

$y(x)={\left(1+\frac{x^2}{3}\right)}^{-1/2}\mathrm{.}$E39

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

$R_k=\frac{\left|y_e(x_k)-y_a(x_k)\right|}{\left|y_e(x_k)\right|},0\le k\le N,$E40

where N is the number of grid points, $y_a(x_k)$ is the approximate solution and $y_e(x_k)$ is the exact solution at the grid point 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.

		N = 4		N = 6
x	Exact	Approximate	Relative error	Approximate	Relative error
0.2	1.060313	1.060313	1.884728e−14	1.060313	4.251109e−014
0.4	1.140018	1.140018	3.505912e−14	1.140018	6.505415e−014
0.6	1.371416	1.371416	7.852585e−14	1.371416	1.066980e−013
0.8	1.787189	1.787189	1.413886e−13	1.787189	1.474757e−013
1.0	2.522988	2.522988	2.355112e−13	2.522988	1.864022e−013
1.2	3.858367	3.858367	3.776361e−13	3.858367	2.234047e−013
1.4	6.391979	6.391979	5.987445e−13	6.391979	2.591455e−013
1.6	11.471251	11.471251	9.227679e−13	11.471251	2.937560e−013
1.8	22.301278	22.301278	1.396995e−12	22.301278	3.316738e−013
2.0	46.966942	46.966942	2.066546e−12	46.966942	3.633883e−013
CPU Time (sec)	0.659414	0.629082		0.659414

Table 1.

Analytical, approximate solutions and relative errors for Example 1.

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⁻¹² and for N = 6, the relative error is on average 10⁻¹³. 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.

Error graph.

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 $x\in [\mathrm{0,2}]$. Table 1 and Figures 1 and 2 validate the accuracy and computational efficiency of the multi-domain pseudospectral relaxation method for Eq. (30).

Figure 2.

Comparison of analytical and approximate solutions for Example 1.

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⁻¹¹. For N = 4 the relative error is approximately 10⁻¹³. 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 $x\in [\mathrm{0,20}]$. A maximum of N = 6 grid points ensured that the numerical method converged to an error of 10⁻¹³ within a fraction of a second.

Figure 3.

Error graph.

Figure 4.

Comparison of analytical and approximate solutions for Example 2.

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 $x\in [\mathrm{0,20}]$. The match between the exact and approximate solutions in Table 2 and Figures 3 and 4 validates the accuracy and computational efficiency of the multi-domain pseudospectral relaxation method for Eq. (31).

		N = 4		N = 6
x	Exact	Approximate	Relative error	Approximate	Relative error
2	–3.271943	–3.271943	8.810376e−011	–3.271943	2.704972e−014
4	–5.697684	–5.697682	1.000287e−011	–5.697684	6.202217e−014
6	–7.243401	–7.243399	8.316853e−012	–7.243401	2.136148e−014
8	–8.365152	–8.365151	1.308867e−011	–8.365152	3.825343e−015
10	–9.243421	–9.243420	1.383787e−011	–9.243421	4.807193e−015
12	–9.964487	–9.964487	1.318604e−011	–9.964487	5.172115e−015
14	–10.575872	–10.575872	1.205994e−011	–10.575872	6.720993e−015
16	–11.106445	–11.106445	1.080384e−011	–11.106445	3.199792e−016
18	–11.575028	–11.575028	9.561315e−012	–11.575028	6.293749e−015
20	–11.927621	–11.927622	8.504839e−012	–11.927621	1.738910e−014
CPU Time (sec)	0.649135	0.602359		0.649135

Table 2.

Analytical, approximate solutions and relative errors for Example 2.

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⁻¹², while for N = 6, the relative error is approximately 10⁻¹³. 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 $x\in [\mathrm{0,2}]$. N = 6 grid points ensured that the numerical method converged to an error of 10⁻¹³ within few seconds.

		N = 4		N = 6
x	Exact	Approximate	Relative error	Approximate	Relative error
0.2	1.029322	1.029322	2.502345e−014	1.029322	1.682611e−014
0.4	1.146713	1.146713	1.316722e−013	1.146713	4.124439e−014
0.6	1.383892	1.383892	3.732052e−013	1.383892	6.915367e−014
0.8	1.809228	1.809228	8.178657e−013	1.809228	1.078787e−013
1.0	2.562286	2.562286	1.581696e−012	2.562286	1.493997e−013
1.2	3.931024	3.931024	2.814656e−012	3.931024	1.883216e−013
1.4	6.53322	6.53322	4.761036e−012	6.53322	2.300241e−013
1.6	11.762306	11.762306	7.722919e−012	11.762306	2.750095e−013
1.8	22.940410	22.940410	1.211045e−011	22.940410	3.176323e−013
2.0	48.467816	48.467816	1.845750e−011	48.467816	3.667955e−013
CPU Time (sec)	1.069773	1.010628		1.069773

Table 3.

Analytical, approximate solutions and relative errors for Example 3.

Figure 5.

Error graph.

Figure 6.

A comparison of analytical and approximate solutions in Example 3.

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 $x\in [\mathrm{0,2}]$. The match between the exact and approximate solutions in 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⁻¹³. For N = 6, the relative error is also approximately 10⁻¹³. 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⁻¹³ in a fraction of seconds.

		N = 4		N = 6
x	Exact	Approximate	Relative error	Approximate	Relative error
0.2	0.974097	0.974097	7.978218e−015	0.974097	2.803774e−014
0.4	0.877179	0.877179	1.215047e−014	0.877179	6.999178e−014
0.6	0.729173	0.729173	1.248514e−014	0.729173	1.079508e−013
0.8	0.559538	0.559538	1.805602e−014	0.559538	1.341305e−013
1.0	0.396355	0.396355	3.221242e−014	0.396355	1.606419e−013
1.2	0.259177	0.259177	6.232707e−014	0.259177	1.756296e−013
1.4	0.156446	0.156446	1.011254e−013	0.156446	1.754615e−013
1.6	0.087174	0.087174	1.525094e−013	0.087174	1.404105e−013
1.8	0.044840	0.044840	1.897188e−013	0.044840	3.528213e−014
2.0	0.021292	0.021292	1.404623e−013	0.021292	2.123230e−013
CPU Time (sec)	1.080828	1.046274		1.080828

Table 4.

Analytical, approximate solutions and relative errors for Example 4.

Figure 7.

Error graph.

Figure 8.

A comparison of analytical and approximate solutions in Example 4.

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 $x\in [\mathrm{0,2}]$. 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⁻¹⁴. For N = 6, the relative error is approximately 10⁻¹⁴. 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 $x\in [\mathrm{0,2}]$. The numerical method converged to an error of 10⁻¹⁴ in a fraction of a second.

		N = 4		N = 6
x	Exact	Approximate	Relative error	Approximate	Relative error
0.2	0.026244	0.026244	4.759186e−015	0.026244	1.692155e−014
0.4	0.131044	0.131044	1.397903e−014	0.131044	3.473577e−014
0.6	0.315844	0.315844	2.495721e−014	0.315844	5.518706e−014
0.8	0.580644	0.580644	3.537301e−014	0.580644	6.864276e−014
1.0	0.925444	0.925444	4.138845e−014	0.925444	8.445643e−014
1.2	1.350244	1.350244	3.880967e−014	1.350244	9.981979e−014
1.4	1.855044	1.855044	4.093663e−014	1.855044	1.098825e−013
1.6	2.439844	2.439844	3.749517e−014	2.439844	8.700337e−014
1.8	3.104644	3.104644	2.445989e−014	3.104644	4.706026e−014
2.0	3.849444	3.849444	4.614580e−016	3.849444	6.575777e−015
CPU Time (sec)	0.705910	0.638383		0.705910

Table 5.

Analytical, approximate solutions and relative errors for Example 5.

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 $x\in [\mathrm{0,2}]$. Table 5 and Figures 9 and 10 validate the accuracy and computational efficiency of the multi-domain pseudospectral relaxation method for Eq. (34).

Figure 9.

Error graph.

Figure 10.

A comparison of analytical and approximate solutions in Example 5.

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 ${10}^{-11}$ and 10⁻¹³ 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 $x\in [\mathrm{0,20}]$. N = 6 grid points ensured that the numerical method converged to an error of 10⁻¹³ within a few seconds as shown in Table 6.

		N = 4		N = 6
x	Exact	Approximate	Relative error	Approximate	Relative error
2	0.650926	0.650926	1.991857e−011	0.650926	3.189482e−014
4	0.395693	0.395693	2.135316e−011	0.395693	3.745701e−014
6	0.276499	0.276499	6.204217e−012	0.276499	5.842239e−014
8	0.211100	0.211100	8.877991e−012	0.211100	1.519920e−013
10	0.170333	0.170333	2.323778e−011	0.170333	2.328536e−013
12	0.142624	0.142624	3.707619e−011	0.142624	3.096185e−013
14	0.122609	0.122609	5.055364e−011	0.122609	4.021536e−013
16	0.107492	0.107492	6.378059e−011	0.107492	4.983484e−013
18	0.095677	0.095607	7.683911e−011	0.095677	5.994828e−013
20	0.087057	0.087057	8.847457e−011	0.087057	6.926360e−013
CPU Time (sec)	1.199892	1.046417		1.199892

Table 6.

Analytical, approximate solutions and relative errors for Example 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 $x\in [\mathrm{0,20}]$. 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).

Figure 11.

Plot showing solutions to Eq. (35) for some values of m.

Figure 12.

Error graph.

Figure 13.

A comparison of analytical and approximate solutions in Example 6.

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.

Acknowledgments

This work is based on the research supported in part by the National Research Foundation of South Africa (Grant No: 85596).