**2. A multi-domain pseudospectral relaxation method**

0 0 *y y f xy gx x y y* ( , ) = ( ), [0, ), (0) = , (0) = *<sup>x</sup>*

where the prime denotes differentiation with respect to *x*, *f*(*x, y*) is a non-linear function, *g*(*x*) is a prescribed function and *γ*, *α*0, *β*<sup>0</sup> 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

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 0,1 on *x* axis. Despite this limitation, analytical approaches have been found to be desirable because

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 colloca‐

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 nonoverlapping 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 neigh‐ bouring 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

a

 b

¢ (1)

¢¢ ¢ + + Î ¥

they easily overcome the difficulty caused by the singularity at *x* = 0.

tion methods to overcome the difficulty of dealing with the singular point.

g

144 Numerical Simulation - From Brain Imaging to Turbulent Flows

numerical methods of solution.

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

$$
\chi'(\mathbf{x}) \equiv h(\mathbf{x}), \qquad \chi(0) \equiv \alpha\_0,\tag{2}
$$

$$h'(\mathbf{x}) + \frac{\mathcal{I}}{\mathcal{X}}h + f(\mathbf{x}, \mathbf{y}) = \mathbf{g}(\mathbf{x}), \quad h(\mathbf{0}) = \mathcal{J}\_0,\tag{3}$$

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

$$\mathbf{y}'\_{i+1}(\mathbf{x}) \equiv h\_i(\mathbf{x}), \qquad \mathbf{y}\_{i+1}(0) \equiv \alpha\_0,\tag{4}$$

$$h\_{l+1}'(\mathbf{x}) + \frac{\gamma}{\chi} h\_{l+1}(\mathbf{x}) + f(\mathbf{x}, \mathbf{y}\_{l+1}) = \mathbf{g}(\mathbf{x}), \qquad h\_{l+1}(0) = \beta\_0. \tag{5}$$

Assuming that an initial approximation *h*0(*x*) is given, Eq. (4) can be solved for *yi*+1(*x*) and the solution can be used immediately in Eq. (5) which is, in turn, solved for *hi*+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* ∈Λ, where Λ= *a*, *b* is the interval where the solution of Eq. (1) exists. The sub-intervals are defined as

$$\mathbf{A}\_{k} = (\mathbf{x}\_{k-1}, \mathbf{x}\_{k}), \quad k = 1, 2, \dots, p, \quad \text{with}, \quad a = \mathbf{x}\_{0} < \mathbf{x}\_{1} < \mathbf{x}\_{2} < \dots < \mathbf{x}\_{p} = b. \tag{6}$$

The solution procedure assumes that the solution at each sub-interval Λ*<sup>k</sup>* , denoted by *<sup>y</sup>* (*<sup>k</sup>* ) (*x*), can be approximated by a Lagrange interpolation polynomial of the form

$$y^{(k)}(\mathbf{x}) \sim \sum\_{s=0}^{N} y^{(k)}(\mathbf{x}\_s) L\_s(\mathbf{x}) \tag{7}$$

which interpolates *y* (*<sup>k</sup>* ) (*x*) at selected points, chosen to be the Chebyshev-Gauss-Lobatto, defined by

$$\{\{\pi\_s\}\_{s=0}^N = \left\{\cos\left(\frac{\pi s}{N}\right)\right\}\_{s=0}^N. \tag{8}$$

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 *τs*. The function *L*s(x) is the characteristic Lagrange cardinal polynomial defined as

$$L\_s(\mathbf{x}) = \prod\_{s=0 \atop s \neq k}^{M} \frac{\mathbf{x} - \mathbf{x}\_k}{\mathbf{x}\_s - \mathbf{x}\_k},\tag{9}$$

that obey the Kronecker delta equation:

$$L\_s(\mathbf{x}\_k) = \mathcal{S}\_{st} = \begin{cases} 0 & \text{if} \quad s \neq k \\ 1 & \text{if} \quad s = k \end{cases} \tag{10}$$

The solution method seeks to solve each of the system of first-order ODEs independently in each *k*th sub-interval Λ*k* using the solutions obtained in the preceding interval Λ*<sup>k</sup>* <sup>−</sup>1 as initial conditions. In the first interval the solution is computed on *x*0, *<sup>x</sup>*1 and is labelled *<sup>y</sup>* (1)(*x*). The value of the solution at the last node *<sup>y</sup>* (1)(*x*1) is used as an initial condition when seeking a solution in the second sub-interval Λ2. This procedure is repeated in each Λ*<sup>k</sup>* interval with the solutions marched across each interval using the following continuity condition:

$$h^{(k)}(\mathbf{x}\_{k-1}) = \mathbf{y}^{(k-1)}(\mathbf{x}\_{k-1}), \quad h^{(k)}(\mathbf{x}\_{k-1}) \equiv h^{(k-1)}(\mathbf{x}\_{k-1}).\tag{11}$$

Therefore, in each sub-interval Λ*<sup>k</sup>* , the following system of ODEs is solved:

$$\frac{d\mathbf{y}\_{i+1}^{(k)}(\mathbf{x})}{d\mathbf{x}} = h\_i^{(k)}(\mathbf{x}), \quad \mathbf{y}\_{i+1}^{(k)}(\mathbf{x}\_{k-1}) = \alpha\_{k-1} \tag{12}$$

A Multi-Domain Spectral Collocation Approach for Solving Lane-Emden Type Equations http://dx.doi.org/10.5772/63016 147

$$\frac{d h\_{i+1}^{(k)}(\mathbf{x})}{d \mathbf{x}} + \frac{\mathcal{Y}}{\mathbf{x}} h\_{i+1}^{(k)}(\mathbf{x}) = \mathbf{g}(\mathbf{x}) - f(\mathbf{x}, \mathbf{y}\_{i+1}^{(k)}), \quad h\_{i+1}^{(k)}(\mathbf{x}\_{k-1}) = \mathcal{B}\_{k-1}, \tag{13}$$

where

 ∑ ( ) ( ) =0 ( ) ( ) ( ),

=0

*N s s*

*s*

d

solutions marched across each interval using the following continuity condition:

( ) ( 1) ( ) ( 1) 1 11 1 ( ) = ( ), ( ) = ( ). *kk kk k kk k yx y x hx h x* - -

Therefore, in each sub-interval Λ*<sup>k</sup>* , the following system of ODEs is solved:

1 ( ) ( )

*i k k*

*dy x hx yx dx*

( ) ( ), ( )

*i ik k*

( )

*k*

*L x*

t

{ } = cos .

=0 ( )= , *M*

*x x L x x x* ¹

*<sup>s</sup> s k s k*

0 if ( )= = . 1 if = . *s k sk*

The solution method seeks to solve each of the system of first-order ODEs independently in each *k*th sub-interval Λ*k* using the solutions obtained in the preceding interval Λ*<sup>k</sup>* <sup>−</sup>1 as initial conditions. In the first interval the solution is computed on *x*0, *<sup>x</sup>*1 and is labelled *<sup>y</sup>* (1)(*x*). The value of the solution at the last node *<sup>y</sup>* (1)(*x*1) is used as an initial condition when seeking a solution in the second sub-interval Λ2. This procedure is repeated in each Λ*<sup>k</sup>* interval with the

<sup>ì</sup> <sup>¹</sup> <sup>í</sup>

*s k*

*s k*


11 1

a<sup>+</sup> = = +- - (12)

which interpolates *y* (*<sup>k</sup>* )

146 Numerical Simulation - From Brain Imaging to Turbulent Flows

cardinal polynomial defined as

that obey the Kronecker delta equation:

defined by

*s*

*s s*

=0

*s*

*k*

*N*

*s N* p

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 *τs*. The function *L*s(x) is the characteristic Lagrange

*y x y xLx* (7)

ì ü æ ö í ý ç ÷ î þ è ø (8)


<sup>î</sup> (10)

(*x*) at selected points, chosen to be the Chebyshev-Gauss-Lobatto,

*<sup>N</sup> k k*

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

Eqs. (12) and (13) cannot be solved exactly because of the non-linear function *f*(*x, y*). Accord‐ ingly, we employ the spectral collocation at the grid points *xj* for *j* =0,1,2,⋯, *N* for each subinterval Λ*<sup>k</sup>* . Before the spectral method is applied, each sub-interval *xk* <sup>−</sup>1, *xk* is transformed to −1,1 using the linear transformation:

$$\mathbf{x} = \frac{\mathbf{x}\_k - \mathbf{x}\_{k-1}}{2}\boldsymbol{\tau} + \frac{\mathbf{x}\_k + \mathbf{x}\_{k-1}}{2}, \quad \boldsymbol{\tau} \in [-1, 1]. \tag{14}$$

The derivatives at the collocation points are evaluated as

$$\left. \frac{d\mathbf{y}^{(k)}}{d\mathbf{x}} \right|\_{\mathbf{x}=\mathbf{x}\_{j}} = \sum\_{s=0}^{N} \mathbf{y}^{(k)}(\mathbf{x}\_{s}) \frac{d L\_{i}(\mathbf{x}\_{j})}{d\mathbf{x}} = \frac{2}{\Delta \mathbf{x}\_{k}} \sum\_{s=0}^{N} \mathbf{y}^{(k)}(\boldsymbol{\pi}\_{s}) \frac{d L\_{i}(\boldsymbol{\pi}\_{j})}{d\boldsymbol{\pi}} = \sum\_{s=0}^{N} D\_{js} \mathbf{y}^{(k)}(\boldsymbol{\pi}\_{s}), \tag{15}$$

where Δ*xk* <sup>=</sup> *xk* <sup>−</sup> *xk* <sup>−</sup>1, **<sup>D</sup>** *js* <sup>=</sup> <sup>2</sup> <sup>Δ</sup>*xk Djs* with *Djs* <sup>=</sup> *d L <sup>s</sup>*(*τ<sup>j</sup>* ) *<sup>d</sup>τ* being the *j*th and *s*th entry of the standard first derivative Chebyshev differentiation matrix of size (*N* + 1)×(*N* + 1) as defined in [27]. Evaluating Eqs. (12) and (13) at the grid points *τ<sup>j</sup>* for *j* =0,1,2,⋯, *N* gives

$$\sum\_{s=0}^{N} \mathbf{D}\_{js} \mathbf{y}\_{i+1}^{(k)}(\boldsymbol{\tau}\_s) = h\_i^{(k)}(\boldsymbol{\tau}\_j), \quad \mathbf{y}\_{i+1}^{(k)}(\boldsymbol{\tau}\_N) = \boldsymbol{\alpha}\_{k-1}, \tag{16}$$

$$\sum\_{s=0}^{N} \mathbf{D}\_{ji} h\_{i+1}^{(k)}(\boldsymbol{\tau}\_{\cdot}) + \frac{\mathcal{T}}{\mathbf{x}\_{j}} h\_{i+1}^{(k)}(\boldsymbol{\tau}\_{\cdot}) = \mathbf{g}(\mathbf{x}\_{j}) - f(\mathbf{x}\_{s}, \mathbf{y}\_{i+1}^{(k)}(\boldsymbol{\tau}\_{\cdot})), \quad h\_{i+1}^{(k)}(\boldsymbol{\tau}\_{\cdot}) = \mathcal{B}\_{k-1}.\tag{17}$$

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

$$\sum\_{s=0}^{N-1} \mathbf{D}\_{js} \mathbf{y}\_{i+1}^{(k)}(\boldsymbol{\tau}\_s) = \boldsymbol{V}\_i^{(k)}(\boldsymbol{\tau}\_j),\tag{18}$$

$$\sum\_{s=0}^{N-1} \mathbf{D}\_{ja} h\_{i+1}^{(k)}(\boldsymbol{\tau}\_s) + \frac{\mathcal{Y}}{\boldsymbol{\varkappa}\_j} h\_{i+1}^{(k)}(\boldsymbol{\tau}\_j) = W\_i^{(k)}(\boldsymbol{\tau}\_j). \tag{19}$$

where

$$W\_{i}^{(k)}(\\\tau\_{j}) \equiv h\_{i}^{(k)}(\\\tau\_{j}) - \mathbf{D}\_{jN} \mathbf{y}\_{i+1}^{(k)}(\\\tau\_{N}),\tag{20}$$

$$W\_i^{(k)}(\\\tau\_j) = \mathbf{g}(\mathbf{x}\_j) - f(\mathbf{x}\_j, \mathbf{y}\_{i+1}^{(k)}(\tau\_j)) - \mathbf{D}\_{jN} h\_{i+1}^{(k)}(\tau\_N),\tag{21}$$

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

$$\mathbf{D} \mathbf{Y}\_{i+1}^{(k)} = \mathbf{V}\_i^{(k)},\tag{22}$$

$$\mathbf{H}\left(\mathbf{D} + \boldsymbol{\zeta}\right)\mathbf{H}\_{i+1}^{(k)} = \mathbf{W}\_{i}^{(k)},\tag{23}$$

where the vectors **Y**(*<sup>k</sup>* ) , **H**(*<sup>k</sup>* ) , **V**(*<sup>k</sup>* ) and **W**(*<sup>k</sup>* ) are respectively given by

$$\mathbf{Y}^{(k)} = \left[ \mathbf{y}^{(k)}(\tau\_0), \mathbf{y}^{(k)}(\tau\_1), \mathbf{y}^{(k)}(\tau\_2), \dots, \mathbf{y}^{(k)}(\tau\_{N-1}) \right]^{\mathrm{r}},\tag{24}$$

$$\mathbf{H}^{(k)} = \left[ h^{(k)}(\tau\_0), h^{(k)}(\tau\_1), h^{(k)}(\tau\_2), \dots, h^{(k)}(\tau\_{N-1}) \right]^\top,\tag{25}$$

$$\mathbf{V}^{(k)} = \left[ \mathbf{v}^{(k)}(\tau\_0), \mathbf{v}^{(k)}(\tau\_1), \mathbf{v}^{(k)}(\tau\_2), \dots, \mathbf{v}^{(k)}(\tau\_{N-1}) \right]^r,\tag{26}$$

$$\mathbf{W}^{(k)} = \left[ \boldsymbol{\omega}^{(k)}(\tau\_0), \boldsymbol{\omega}^{(k)}(\tau\_1), \boldsymbol{\omega}^{(k)}(\tau\_2), \dots, \boldsymbol{\omega}^{(k)}(\tau\_{N-1}) \right]^{\mathsf{T}},\tag{27}$$

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

$$
\zeta = \begin{bmatrix}
\frac{\mathcal{Y}}{\mathcal{X}\_0} & & & \\
& \ddots & & \\
& & \mathcal{Y} & \\
& & \ddots & \\
& & & \mathcal{Y} \\
& & & & \mathcal{X}\_{N-1}
\end{bmatrix}.
\tag{28}
$$

The approximate values of *y* (*<sup>k</sup>* ) (*x*) and *d y* (*<sup>k</sup>* ) (*x*) *dx* are obtained by iteratively solving Eqs. (22) and (23), respectively, for *i* =0,1,2,…, starting from suitable initial approximation.
