**2. The MSRM for the model problem**

We seek an approximate solution to

$$y'' + \frac{\chi}{x}y' + f(\mathbf{x})\mathbf{g}(\mathbf{y}) = \mathbf{0}, \mathbf{0} < \mathbf{x} \le L,\\ y(\mathbf{0}) = a, \mathbf{y}'(\mathbf{0}) = \mathbf{0} \tag{5}$$

where *γ* >0, *L*, *α* are given constants and *f* and *g* are given functions.

We follow the idea behind the solution method by Ramos [23], where the singularity at *x* ¼ 0 is isolated in a sufficiently small subinterval *I<sup>ε</sup>* ¼ ½ � 0, *ε* of *I* ¼ ½ � 0, *L* where *ε* >0. The point *x* ¼ *ε* splits interval *I* into two subintervals: *I<sup>ε</sup>* and *I* � *I<sup>ε</sup>* ¼ ½ � *ε*, *L* . A linearisation method is used to solve Eq. (5) restricted to *Iε*, i.e., *near the singularity* at *x* ¼ 0. In order to improve the accuracy of the method on the subinterval *I* � *I<sup>ε</sup>* we form a partition

$$I - I\_{\varepsilon} = \downarrow\_{m=1}^{M} I\_{m} \text{ where} \\ I\_{m} = [\mathfrak{x}\_{m-1}, \mathfrak{x}\_{m}], \mathfrak{x}\_{0} = \varepsilon \text{ and} \\ \mathfrak{x}\_{M} = L \tag{6}$$

The method by Ramos proceeds by using the same linearisation method on the subintervals *Im*,ð Þ *m* ¼ 1, ⋯, *M* of *I* � *Iε*, i.e. *away from the singularity*. To this end, at the interface *xm*�<sup>1</sup> of *Im*�<sup>1</sup> and *Im* we make use the solution of Eq. (3) restricted to *Im*�<sup>1</sup> to generate the initial conditions for Eq. (3) restricted to *Im*. This ensures continuity of the solution. In this chapter we avoid linearisation on *I* � *I<sup>ε</sup>* by making use of the SRM on the subintervals of *I* � *Iε*. However, as was done in [23] we make use of linearisation on *Iε*. This approach results in the MSRM. A detailed description of the MSRM is given in Sections 2.1 and 2.2.

#### **2.1 Near the singularity**

Let *ε*∈ð Þ 0, *L* be a sufficiently small number. Restrict problem (5) to 0, ½ � *ε* and re-arrange to get

symmetric, polytropic fluid. The equation provides a useful approximation for

<sup>þ</sup> *qy* <sup>¼</sup> *h x*, *<sup>y</sup>*, *<sup>y</sup>*<sup>0</sup> ð Þ whereðÞ<sup>0</sup> <sup>¼</sup> *<sup>d</sup>*

an equation which was discussed in [6]. An existence result for the solution is

Exact solutions are available for particular cases of Eq. (3) [7], but not for the general case according to the best of our knowledge. This is motivation enough for seeking approximate solutions. To this end several approximate analytical methods were used by other researchers to solve Eq. (3). Van Gorder [8] made use of the Homotopy analysis method (HAM) and its variant, the Optimal homotopy analysis method, to solve a boundary value problem for the Lane-Emden equation of the second kind. The two respective analytical solutions that they obtained were in strong agreement. The Homotopy perturbation method (HPM) is another variant of the HAM that was used by Chowdhury and Hashim [9] to solve an initial value problem for Eq. (3). Their analytical solutions were the same as those that were obtained by Wazwaz [10] using the Adomian decomposition method (ADM). Chowdhury and Hashim observed that the HPM was less computationally expensive than the ADM. Wazwaz [11] made use of the variational iteration method (VIM) to solve both initial value problems and boundary value problems for Eq. (3) and for some inhomogeneous Emden-Fowler equations. The results that they

Some numerical methods have been used by other researchers to approximate solutions to Eq. (3). Many of these numerical methods fall in the class of *collocation methods*. Examples of these collocation methods include but are not limited to the Chebyshev wavelet finite difference method (CWFDM) [12], the Haar wavelet collocation method (HWCM) [13], the Taylor wavelet method (TWM) [14] and the Radial basis function - differential quadrature method (RDF-DQM) [15]. One distinct feature of these collocation methods is the choice of collocation points for discretizing the problem domain. Another distinct feature is the choice of basis functions that are used either for constructing numerical solutions or for numerical differentiation. The CWFDM makes use of Chebyshev-Gauss-Lobatto collocation points and Chebyshev wavelet finite difference basis functions. The HWCM makes

*x <sup>j</sup>* ¼ ð Þ *j* � 0*:*5 *=*2*M*, *j* ¼ 1, ⋯, 2*M*

on the problem domain 0, ½ � *L* , and the method uses integrals of Haar wavelets as basis functions. The TWM uses roots of shifted Legendre polynomials as collocation points, and as basis functions the method uses Taylor wavelets which are special functions that defined in terms of Taylor polynomials. Convergence results for the Taylor wavelet solution were presented in [14]. The RDF-DQM uses collocation points

*N* � 1

, *j* ¼ 1, ⋯, *N*

<sup>1</sup> � cos *<sup>j</sup>* � <sup>1</sup>

*dx* <sup>þ</sup> *f x*ð Þ*g y*ð Þ¼ <sup>0</sup> (3)

*dx* ðÞ, (4)

A more general form for Eq. (2) is the Emden-Fowler equation

given therein under certain conditions on *p x*ð Þ, *q x*ð Þ and *h x*, *y*, *y*<sup>0</sup> ð Þ.

obtained demonstrated the reliability and effectiveness of the VIM.

*d*2 *y dx*<sup>2</sup> <sup>þ</sup> *<sup>γ</sup> x dy*

*py*<sup>0</sup> ð Þ<sup>0</sup>

stars.

which can be written as

*Recent Advances in Numerical Simulations*

use of the collocation points

**116**

*<sup>x</sup> <sup>j</sup>* <sup>¼</sup> <sup>2</sup> *L* *Recent Advances in Numerical Simulations*

$$\mathbf{y}'' = \underbrace{-\left[\frac{\mathbf{y}}{\mathbf{x}}\mathbf{y}' + f(\mathbf{x})\mathbf{g}(\mathbf{y})\right]}\_{\mathbf{u}(\mathbf{x},\mathbf{y},\mathbf{y}')}, \mathbf{0} < \mathbf{x} \le \mathbf{e}, \mathbf{y}(\mathbf{0}) = a, \mathbf{y}'(\mathbf{0}) = \mathbf{0} \tag{7}$$

The first step of the SRM for (11) is to let *v* ¼ *y* and *w* ¼ *v*<sup>0</sup> so that we obtain the

*w* þ *f x*ð Þ*g v*ð Þ¼ 0, *w*ð Þ¼ *ε α*<sup>0</sup>

*L* � *ε* <sup>2</sup> *<sup>η</sup>*

*w* þ *f x*ð Þ ð Þ*η g v*ð Þ¼ 0, *w*ð Þ¼ �1 *α*<sup>0</sup>

for �1≤ *η*≤ 1 where *β* ¼ 2*=*ð Þ L � *ε :* As described in [19] the next step of the SRM

mimicks the Gauss–Seidel method for linear systems and it yields the iteration

*wr*þ<sup>1</sup> ¼ �*f x*ð Þ ð Þ*η g v*ð Þ *<sup>r</sup>*þ<sup>1</sup> |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} *<sup>R</sup>*ð Þ<sup>2</sup> ð Þ*<sup>η</sup>*

for �1 ¼ *η*<sup>N</sup> ≤*η*≤*η*<sup>0</sup> ¼ 1 where *r* ¼ 0, 1, ⋯ and on ½ � �1, 1 we formed a grid

If the initial approximations *v*<sup>0</sup> and *w*<sup>0</sup> to *v* and *w*, respectively, are prescribed

approximations. To this end we assume that *vr* and *wr* are known at the end of the *<sup>r</sup>*th iteration. Once Eq. (16) is solved for *vr*þ<sup>1</sup> the right hand side *<sup>R</sup>*ð Þ<sup>2</sup> of Eq. (17) becomes known and we solve this equation for *wr*þ1. Since *vr*þ<sup>1</sup> and *wr*þ<sup>1</sup> are now known, we proceed in a similar manner to compute *vr*þ<sup>2</sup> and *wr*þ2, and so on. As done in [19] we solve Eqs. (16) and (17) by making use of Chebyshev differentia-

*<sup>d</sup><sup>η</sup>* <sup>¼</sup> *wr* |{z} *<sup>R</sup>*ð Þ<sup>1</sup> ð Þ*<sup>η</sup>*

*L* þ *ε* 2 þ

*<sup>w</sup>*<sup>0</sup> <sup>þ</sup> *<sup>γ</sup> x*

*DOI: http://dx.doi.org/10.5772/intechopen.96611*

*β dw dη*

*<sup>β</sup> dwr*þ<sup>1</sup> *dη*

<sup>þ</sup> *<sup>γ</sup> x*ð Þ*η*

þ *γ x*

*<sup>η</sup><sup>i</sup>* <sup>¼</sup> cos *<sup>π</sup><sup>i</sup>*

then Eqs. (16) and (17) generate sequences f g *vr* <sup>∞</sup>

*D*^ |{z} *A*1

where *<sup>D</sup>*^ <sup>¼</sup> *<sup>β</sup>D*, *<sup>D</sup>* is the *<sup>N</sup>* � *<sup>N</sup>* Chebyshev differentiation matrix,

*<sup>D</sup>*^ <sup>þ</sup> *<sup>γ</sup>*diag 1*=<sup>x</sup> <sup>η</sup>*<sup>0</sup> ½ � ð Þ, <sup>⋯</sup>, 1*=x*ð Þ *<sup>η</sup><sup>N</sup>* � � |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} *A*2

tion [24] to obtain

**119**

*<sup>β</sup> dvr*þ<sup>1</sup>

consisting of the Chebyshev-Gauss-Lobatto collocation points

*N*

for *ε*≤*x*≤*L* which upon making use of the change of variable

*A Modified Spectral Relaxation Method for Some Emden-Fowler Equations*

*x*ð Þ¼ *η*

*<sup>β</sup> dv*

*v*<sup>0</sup> ¼ *w*, *v*ð Þ¼ *ε α*0, (12)

*<sup>d</sup><sup>η</sup>* <sup>¼</sup> *<sup>w</sup>*, *<sup>v</sup>*ð Þ¼ �<sup>1</sup> *<sup>α</sup>*<sup>0</sup> (14)

, *vr*þ<sup>1</sup>ð Þ¼ *η<sup>N</sup> α*<sup>0</sup> (16)

, *wr*þ<sup>1</sup>ð Þ¼ *η<sup>N</sup> α*<sup>0</sup>

� �, *<sup>i</sup>* <sup>¼</sup> 0, 1, <sup>⋯</sup>, *<sup>N</sup>:* (18)

*<sup>r</sup>*¼<sup>1</sup> and f g *wr* <sup>∞</sup>

**<sup>v</sup>***<sup>r</sup>*þ<sup>1</sup> <sup>¼</sup> **<sup>R</sup>**ð Þ<sup>1</sup> , *vr*þ<sup>1</sup>ð Þ¼ *<sup>η</sup><sup>N</sup> <sup>α</sup>*0, (19)

**<sup>w</sup>***<sup>r</sup>*þ<sup>1</sup> <sup>¼</sup> **<sup>R</sup>**ð Þ<sup>2</sup> , *wr*þ<sup>1</sup>ð Þ¼ *<sup>η</sup><sup>N</sup> <sup>α</sup>*<sup>0</sup>

<sup>0</sup> (13)

<sup>0</sup> (15)

<sup>0</sup> (17)

*<sup>r</sup>*¼<sup>1</sup> of consecutive

<sup>0</sup> (20)

equivalent problem

becomes

If we Taylor expand *u* about *ε*<sup>0</sup> ¼ *ε*, *y*ð Þ 0 |ffl{zffl} *y*0 , *y*<sup>0</sup> ð Þ 0 |ffl{zffl} *y*0 0 0 B@ 1 CA and neglect higher order

terms we get

$$u(\mathbf{x}, \boldsymbol{y}, \boldsymbol{y}') = \underbrace{u(\mathbf{e}\_0)}\_{\boldsymbol{u}\_0} + \underbrace{u\_{\mathcal{V}}(\mathbf{e}\_0)}\_{H\_0} (\boldsymbol{y} - \boldsymbol{y}\_0) + \underbrace{u\_{\mathcal{V}}(\mathbf{e}\_0)}\_{G\_0} (\boldsymbol{y}' - \boldsymbol{y}\_0') + \underbrace{u\_{\mathbf{x}}(\mathbf{e}\_0)}\_{L\_0} (\boldsymbol{x} - \boldsymbol{e})$$

where *us* denotes *<sup>∂</sup><sup>u</sup> ∂s* . Consequently, Eq. (7) can be replaced by

$$y'' = \mathfrak{u}\_0 + H\_0(\mathfrak{y} - \mathfrak{y}\_0) + G\_0(\mathfrak{y}' - \mathfrak{y}\_0') + L\_0(\mathfrak{x} - \mathfrak{e}),\\ \mathbf{0} \le \mathfrak{x} \le \mathfrak{e},\\ y(\mathbf{0}) = a, \mathfrak{y}'(\mathbf{0}) = \mathbf{0},\tag{8}$$

where the differential equation now holds at *x* ¼ 0 because the singularity there is no longer present. If *<sup>H</sup>*<sup>0</sup> 6¼ 0 and *<sup>R</sup>*<sup>0</sup> <sup>≔</sup> ð Þ *<sup>G</sup>*0*=*<sup>2</sup> <sup>2</sup> <sup>þ</sup> *<sup>H</sup>*<sup>0</sup> <sup>&</sup>gt; 0, then problem (8) has analytical solution [23]

$$y(\mathbf{x}) = A\_0 \exp\left(\lambda\_0^+(\mathbf{x} - \boldsymbol{\varepsilon})\right) + B\_0 \exp\left(\lambda\_0^-(\mathbf{x} - \boldsymbol{\varepsilon})\right) + C\_0(\mathbf{x} - \boldsymbol{\varepsilon}) + D\_0 \tag{9}$$

for 0≤ *x*≤ *ε*, where *λ*� <sup>0</sup> <sup>¼</sup> *<sup>G</sup>*0*=*<sup>2</sup> � ffiffiffiffiffi *R*0 <sup>p</sup> ,*C*<sup>0</sup> ¼ �*L*0*=H*0,

$$D\_0 = -(G\_0 C\_0 + P\_0) / H\_0,\\ P\_0 = u\_0 - H\_0 y\_0 - G\_0 y\_0',$$

$$A\_0 = \frac{\exp\left(\lambda\_0^+ \varepsilon\right)}{\lambda\_0^+ - \lambda\_0^-} \left(y\_0' - C\_0 - \lambda\_0^- \left(y\_0 + C\_0 \varepsilon - D\_0\right)\right),$$

$$B\_0 = \frac{\exp\left(\lambda\_0^- \varepsilon\right)}{\lambda\_0^+ - \lambda\_0^-} \left(\lambda\_0^+ \left(y\_0 + C\_0 \varepsilon - D\_0\right) - \left(y\_0' - C\_0\right)\right)$$

Thus

$$\mathcal{Y}(\varepsilon) = \underbrace{A\_0 + B\_0 + D\_0}\_{a\_0} \text{ and } \mathcal{Y}(\varepsilon) = \underbrace{\lambda\_0^+ A\_0 + \lambda\_0^- B\_0 + C\_0}\_{d\_0'} \tag{10}$$

We take *α*<sup>0</sup> and *α*<sup>0</sup> <sup>0</sup> as initial values for problem (7) restricted to ½ � *ε*, *L* in the next section.

#### **2.2 Away from the singularity**

We seek *y x*ð Þ satisfying

$$y'' + \frac{\chi}{\varkappa}y' + f(\varkappa)\mathbf{g}(\mathbf{y}) = \mathbf{0}, \boldsymbol{\varepsilon} \le \boldsymbol{\varkappa} \le L,\\ y(\boldsymbol{\varepsilon}) = a\_0,\\ y'(\boldsymbol{\varepsilon}) = a\_0' \tag{11}$$

In this section we begin by describing the SRM for problem (11). In practical applications it is usually important to obtain a solution to (11) which possesses a prescibed degree of accuracy. To this end we make use of the SRM on a partition of the problem domain ½ � *ε*, *L* . This is our last task in this section.

*A Modified Spectral Relaxation Method for Some Emden-Fowler Equations DOI: http://dx.doi.org/10.5772/intechopen.96611*

The first step of the SRM for (11) is to let *v* ¼ *y* and *w* ¼ *v*<sup>0</sup> so that we obtain the equivalent problem

$$v' = w,\\ v(\varepsilon) = a\_0,\tag{12}$$

$$
\omega w' + \frac{\chi}{\varkappa} w + f(\mathfrak{x}) \mathfrak{g}(v) = \mathbf{0}, \\
w(\varepsilon) = a'\_0 \tag{13}
$$

for *ε*≤*x*≤*L* which upon making use of the change of variable

$$\mathfrak{x}(\eta) = \frac{L+\varepsilon}{2} + \frac{L-\varepsilon}{2}\eta$$

becomes

*<sup>y</sup>*<sup>00</sup> ¼ � *<sup>γ</sup>*

*Recent Advances in Numerical Simulations*

*u x*, *y*, *y*<sup>0</sup> ð Þ¼ *u*ð Þ *ε*<sup>0</sup>

where *us* denotes *<sup>∂</sup><sup>u</sup>*

*y*<sup>00</sup> ¼ *u*<sup>0</sup> þ *H*<sup>0</sup> *y* � *y*<sup>0</sup>

analytical solution [23]

Thus

section.

**118**

We take *α*<sup>0</sup> and *α*<sup>0</sup>

**2.2 Away from the singularity**

We seek *y x*ð Þ satisfying

*<sup>y</sup>*<sup>00</sup> <sup>þ</sup> *<sup>γ</sup>*

*y x*ð Þ¼ *A*<sup>0</sup> exp *λ*<sup>þ</sup>

for 0≤ *x*≤ *ε*, where *λ*�

terms we get

*<sup>x</sup> <sup>y</sup>*<sup>0</sup> <sup>þ</sup> *f x*ð Þ*g y*ð Þ h i

, 0 <*x*≤*ε*, *y*ð Þ¼ 0 *α*, *y*<sup>0</sup>

1

*y*<sup>0</sup> � *y*<sup>0</sup> 0 � � <sup>þ</sup> *ux*ð Þ *<sup>ε</sup>*<sup>0</sup>

, *y*<sup>0</sup> ð Þ 0 |ffl{zffl} *y*0 0


� � <sup>þ</sup> *<sup>L</sup>*0ð Þ *<sup>x</sup>* � *<sup>ε</sup>* , 0 <sup>≤</sup>*x*≤*ε*, *<sup>y</sup>*ð Þ¼ <sup>0</sup> *<sup>α</sup>*, *<sup>y</sup>*<sup>0</sup>


� � <sup>þ</sup> *uy*0ð Þ *<sup>ε</sup>*<sup>0</sup>

. Consequently, Eq. (7) can be replaced by

where the differential equation now holds at *x* ¼ 0 because the singularity there is no longer present. If *<sup>H</sup>*<sup>0</sup> 6¼ 0 and *<sup>R</sup>*<sup>0</sup> <sup>≔</sup> ð Þ *<sup>G</sup>*0*=*<sup>2</sup> <sup>2</sup> <sup>þ</sup> *<sup>H</sup>*<sup>0</sup> <sup>&</sup>gt; 0, then problem (8) has

*R*0

*D*<sup>0</sup> ¼ �ð Þ *G*0*C*<sup>0</sup> þ *P*<sup>0</sup> *=H*0, *P*<sup>0</sup> ¼ *u*<sup>0</sup> � *H*0*y*<sup>0</sup> � *G*0*y*<sup>0</sup>

<sup>0</sup> � *C*<sup>0</sup> � *λ*�

and *y*<sup>0</sup>

*<sup>x</sup> <sup>y</sup>*<sup>0</sup> <sup>þ</sup> *f x*ð Þ*g y*ð Þ¼ 0, *<sup>ε</sup>*<sup>≤</sup> *<sup>x</sup>*<sup>≤</sup> *<sup>L</sup>*, *<sup>y</sup>*ð Þ¼ *<sup>ε</sup> <sup>α</sup>*0, *<sup>y</sup>*<sup>0</sup>

the problem domain ½ � *ε*, *L* . This is our last task in this section.

In this section we begin by describing the SRM for problem (11). In practical applications it is usually important to obtain a solution to (11) which possesses a prescibed degree of accuracy. To this end we make use of the SRM on a partition of

<sup>0</sup> *y*<sup>0</sup> þ *C*0*ε* � *D*<sup>0</sup> � � � *<sup>y</sup>*<sup>0</sup>

<sup>p</sup> ,*C*<sup>0</sup> ¼ �*L*0*=H*0,

� � � � ,

� � � �

ð Þ¼ *ε λ*<sup>þ</sup>

<sup>0</sup> *y*<sup>0</sup> þ *C*0*ε* � *D*<sup>0</sup>

<sup>0</sup> *A*<sup>0</sup> þ *λ*�

<sup>0</sup> as initial values for problem (7) restricted to ½ � *ε*, *L* in the next


0

B@

*y* � *y*<sup>0</sup>

0

<sup>0</sup> ð Þ *<sup>x</sup>* � *<sup>ε</sup>* � � <sup>þ</sup> *<sup>B</sup>*<sup>0</sup> exp *<sup>λ</sup>*�

<sup>0</sup> <sup>¼</sup> *<sup>G</sup>*0*=*<sup>2</sup> � ffiffiffiffiffi

*y*0

*λ*þ

<sup>0</sup> *<sup>ε</sup>* � �

<sup>0</sup> *<sup>ε</sup>* � �

ð Þ¼ 0 0 (7)

ð Þ *x* � *ε*

ð Þ¼ 0 0,

(8)

(10)

CA and neglect higher order


<sup>0</sup> ð Þ *<sup>x</sup>* � *<sup>ε</sup>* � � <sup>þ</sup> *<sup>C</sup>*0ð Þþ *<sup>x</sup>* � *<sup>ε</sup> <sup>D</sup>*<sup>0</sup> (9)

<sup>0</sup> � *C*<sup>0</sup>

<sup>0</sup> *B*<sup>0</sup> þ *C*<sup>0</sup>

ð Þ¼ *ε α*<sup>0</sup>

<sup>0</sup> (11)

0,


þ *uy*ð Þ *ε*<sup>0</sup> |fflffl{zfflffl} *H*<sup>0</sup>

If we Taylor expand *u* about *ε*<sup>0</sup> ¼ *ε*, *y*ð Þ 0


*∂s*

� � <sup>þ</sup> *<sup>G</sup>*<sup>0</sup> *<sup>y</sup>*<sup>0</sup> � *<sup>y</sup>*<sup>0</sup>

*<sup>A</sup>*<sup>0</sup> <sup>¼</sup> exp *<sup>λ</sup>*<sup>þ</sup>

*<sup>B</sup>*<sup>0</sup> <sup>¼</sup> exp *<sup>λ</sup>*�

*λ*þ <sup>0</sup> � *λ*� 0

*y*ð Þ¼ *ε A*<sup>0</sup> þ *B*<sup>0</sup> þ *D*<sup>0</sup> |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} *<sup>α</sup>*<sup>0</sup>

*λ*þ <sup>0</sup> � *λ*� 0

$$
\beta \frac{dv}{d\eta} = w,\\
v(-1) = a\_0 \tag{14}
$$

$$
\rho \frac{dw}{d\eta} + \frac{\chi}{\varkappa(\eta)} w + f(\varkappa(\eta)) \mathbf{g}(v) = \mathbf{0}, w(-\mathbf{1}) = a'\_0 \tag{15}
$$

for �1≤ *η*≤ 1 where *β* ¼ 2*=*ð Þ L � *ε :* As described in [19] the next step of the SRM mimicks the Gauss–Seidel method for linear systems and it yields the iteration

$$\beta \frac{dv\_{r+1}}{d\eta} = \underbrace{w\_r}\_{\mathbb{R}^{(1)}(\eta)}, v\_{r+1}(\eta\_N) = a\_0 \tag{16}$$

$$\beta \frac{dw\_{r+1}}{d\eta} + \frac{\gamma}{\varkappa} w\_{r+1} = \underbrace{-f(\varkappa(\eta)) \lg(v\_{r+1})}\_{R^{(2)}(\eta)}, w\_{r+1}(\eta\_N) = a'\_0 \tag{17}$$

for �1 ¼ *η*<sup>N</sup> ≤*η*≤*η*<sup>0</sup> ¼ 1 where *r* ¼ 0, 1, ⋯ and on ½ � �1, 1 we formed a grid consisting of the Chebyshev-Gauss-Lobatto collocation points

$$\eta\_i = \cos\left(\frac{\pi i}{N}\right), \qquad \qquad i = 0, 1, \cdots, N. \tag{18}$$

If the initial approximations *v*<sup>0</sup> and *w*<sup>0</sup> to *v* and *w*, respectively, are prescribed then Eqs. (16) and (17) generate sequences f g *vr* <sup>∞</sup> *<sup>r</sup>*¼<sup>1</sup> and f g *wr* <sup>∞</sup> *<sup>r</sup>*¼<sup>1</sup> of consecutive approximations. To this end we assume that *vr* and *wr* are known at the end of the *<sup>r</sup>*th iteration. Once Eq. (16) is solved for *vr*þ<sup>1</sup> the right hand side *<sup>R</sup>*ð Þ<sup>2</sup> of Eq. (17) becomes known and we solve this equation for *wr*þ1. Since *vr*þ<sup>1</sup> and *wr*þ<sup>1</sup> are now known, we proceed in a similar manner to compute *vr*þ<sup>2</sup> and *wr*þ2, and so on. As done in [19] we solve Eqs. (16) and (17) by making use of Chebyshev differentiation [24] to obtain

$$\underbrace{\hat{D}\lrcorner}\_{A\_1}\mathbf{v}\_{r+1} = \mathbf{R}^{(1)}, \mathbf{v}\_{r+1}(\eta\_N) = a\_0,\tag{19}$$

$$\underbrace{\left(\hat{\mathbf{D}} + \mathbf{y} \text{diag}[\mathbf{1}/\mathbf{x}(\eta\_0), \dots, \mathbf{1}/\mathbf{x}(\eta\_N)]\right)}\_{A\_2} \mathbf{w}\_{r+1} = \mathbf{R}^{(2)}, \mathbf{w}\_{r+1}(\eta\_N) = a'\_0 \tag{20}$$

where *<sup>D</sup>*^ <sup>¼</sup> *<sup>β</sup>D*, *<sup>D</sup>* is the *<sup>N</sup>* � *<sup>N</sup>* Chebyshev differentiation matrix,

$$\text{diag}[b\_0, \dots, b\_N] = \begin{pmatrix} b\_0 \\ & \ddots \\ & & b\_N \end{pmatrix} \tag{21}$$

4.**for** *m* ¼ 1, 2, ⋯, *M* **do**

for the estimate ~*y* to *y* on *Im*.

*DOI: http://dx.doi.org/10.5772/intechopen.96611*

6. Set *α*<sup>0</sup> ¼ ~*y x*ð Þ *<sup>m</sup>* and set *α*<sup>0</sup>

computational effeciency of the MSRM.

*y*<sup>00</sup> þ 1

constructing the numerical solution.

*<sup>y</sup>*ð Þ*<sup>ε</sup>* <sup>≔</sup> *<sup>α</sup>*<sup>0</sup> <sup>¼</sup> <sup>9</sup>*:*<sup>999999926</sup> � <sup>10</sup>�<sup>1</sup> and *<sup>y</sup>*<sup>0</sup>

*<sup>i</sup>* <sup>¼</sup> *<sup>w</sup> <sup>η</sup><sup>i</sup>* ð Þ and *<sup>R</sup>*ð Þ<sup>2</sup>

analytical solution

1. Coefficients

for problem (8).

2. Initial values

for problem (11).

*R*ð Þ<sup>1</sup>

3. Entries

**121**

7.**end**

**2.3 Examples**

5. Define *Im* ≔ ½ � *xm*�1, *xm* and make use of the SRM to solve *y*<sup>00</sup> ¼ *u x*, *y*, *y*<sup>0</sup> ð Þ, *x*∈*Im*, *y x*ð Þ¼ *<sup>m</sup>*�<sup>1</sup> *α*0, *y*<sup>0</sup>

*A Modified Spectral Relaxation Method for Some Emden-Fowler Equations*

<sup>0</sup> ¼ ~*y*<sup>0</sup>

Example 1 We look for a numerical solution to the problem

ð Þ *xm* .

In this section we make use of some examples to investigate the accuracy and

*<sup>x</sup> <sup>y</sup>*<sup>0</sup> <sup>þ</sup> <sup>3</sup>*y*<sup>5</sup> � *<sup>y</sup>*<sup>3</sup> <sup>¼</sup> 0, 0 <sup>&</sup>lt;*x*≤10*y*ð Þ¼ <sup>0</sup> 1, *<sup>y</sup>*<sup>0</sup>

which Wazwaz [11] solved using the VIM and obtained the approximate.

*y x*ð Þ¼ <sup>1</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffi

When we apply the MSRM on (25) we get the following components for

ð Þ¼ *xm*�<sup>1</sup> *α*<sup>0</sup>

0

ð Þ¼ 0 0 (25)

<sup>1</sup> <sup>þ</sup> *<sup>x</sup>*<sup>2</sup> <sup>p</sup> (26)

<sup>0</sup> ¼ �1*:*<sup>264241093</sup> � <sup>10</sup>�<sup>4</sup> (28)

*<sup>u</sup>*<sup>0</sup> ¼ �2, *<sup>L</sup>*<sup>0</sup> <sup>¼</sup> 0, *<sup>H</sup>*<sup>0</sup> ¼ �12, *<sup>G</sup>*<sup>0</sup> ¼ �10<sup>4</sup> (27)

ð Þ*ε* ≔ *α*<sup>0</sup>

*<sup>r</sup>*þ<sup>1</sup> *<sup>η</sup><sup>i</sup>* ð Þ� *<sup>v</sup>*<sup>3</sup>

*<sup>r</sup>*þ<sup>1</sup> *<sup>η</sup><sup>i</sup>* ð Þ � �, *<sup>i</sup>* <sup>¼</sup> 0, <sup>⋯</sup>, *<sup>N</sup>*

*<sup>i</sup>* ¼ � <sup>3</sup>*v*<sup>5</sup>

The MSRM generates a numerical solution to problem (25) which is plotted together with the analytical solution (26) in **Figure 1(a)**. **Figure 1(a)** shows a good agreement between the numerical and analytical solutions. For a more detailed comparion of the two solutions see **Table 1**. **Table 1** and **Figure 1(b)** show that the absolute error of the numerical solution is no more than 0*:*<sup>5</sup> � <sup>10</sup>�6. Thus the numerical solution agrees with the analytical solution in the first 6 decimal places.

in the vectors on the right hand sides of Eqs. (19) and (20).

is a diagonal matrix, and **<sup>R</sup>**ð Þ*<sup>i</sup>* <sup>¼</sup> *<sup>R</sup>*ð Þ*<sup>i</sup> <sup>η</sup>*<sup>0</sup> ð Þ, <sup>⋯</sup>*R*ð Þ*<sup>i</sup>* ð Þ *<sup>η</sup><sup>N</sup>* � �*<sup>T</sup>* with *<sup>i</sup>* <sup>¼</sup> 1, 2. Moreover,

$$\mathbf{v}\_r = \begin{bmatrix} v\_r(\eta\_0), \cdots, v\_r(\eta\_N) \end{bmatrix}^T, \mathbf{w}\_r = \begin{bmatrix} w\_r(\eta\_0), \cdots, w\_r(\eta\_N) \end{bmatrix}^T, \qquad r = \mathbf{0}, \mathbf{1}, \cdots \tag{22}$$

and ½ �*<sup>T</sup>* denotes matrix transpose.

We prescribe *v*<sup>0</sup> and *w*<sup>0</sup> by requiring that

$$\nu\_0(\varepsilon) \equiv \mathcal{y}\_0(\varepsilon) = a\_0 \text{ and } \mathfrak{w}\_0(\varepsilon) \equiv \mathcal{y}\_0'(\varepsilon) = a\_0'$$

Moreover, we assume that

$$\lim\_{r \to \infty} v\_r = \nu \text{ and } \lim\_{r \to \infty} w\_r = w\_r$$

The initial conditions

$$\upsilon\_{r+1}(\eta\_N) = a\_0 \text{ and } \upsilon\_{r+1}(\eta\_N) = a\_0'$$

are included in the iterative scheme consisting of Eqs. (19) and (20) as shown below.

$$\mathbf{P}\left(\frac{A\_1}{\mathbf{0}\cdot\mathbf{\color{red}{0}}\mathbf{1}}\right)\begin{pmatrix}\begin{matrix}\boldsymbol{v}\_{r+1}(\boldsymbol{\eta}\_0)\\\vdots\\\boldsymbol{v}\_{r+1}(\boldsymbol{\eta}\_{N-1})\\\hline\boldsymbol{v}\_{r+1}(\boldsymbol{\eta}\_N)\end{matrix}\end{pmatrix} = \begin{pmatrix}\begin{matrix}\boldsymbol{R}^{(1)}(\boldsymbol{\eta}\_0)\\\vdots\\\boldsymbol{R}^{(1)}(\boldsymbol{\eta}\_{N-1})\\\hline\boldsymbol{a}\_0\end{matrix}\end{pmatrix}\tag{23}$$

$$
\chi \left( \frac{A\_2}{\mathbf{0} \cdots \mathbf{0}} \mathbf{1} \right) \begin{pmatrix} w\_{r+1}(\eta\_0) \\ \vdots \\ w\_{r+1}(\eta\_{N-1}) \\ \frac{w\_{r+1}(\eta\_N)}{w\_{r+1}(\eta\_N)} \end{pmatrix} = \begin{pmatrix} R^{(2)}(\eta\_0) \\ \vdots \\ R^{(2)}(\eta\_{N-1}) \\ \frac{R^{(2)}(\eta\_{N-1})}{a\_0'} \end{pmatrix} \tag{24}$$

The larger *L* is the less reliable is the SRM. As a workaround to this problem is we subdivide interval ½ � *ε*, *L* into a disjoint union of non-overlapping subintervals as detailed in Eq. (6). Given the the model problem on *I*1, we use the SRM to compute estimate ~*y* to *y* on *I*1. We make use of ~*y* to compute initial values for the problem on *I*2. We repeat this procedure for the problem on *I*<sup>2</sup> and continue in a similar manner until we exhaust ½ � *ε*, *L* . Shown in Algorithm 2.2 is an outline of the MSRM for problem (5).

**Algorithm 1**: Putting the MSRM together.


*A Modified Spectral Relaxation Method for Some Emden-Fowler Equations DOI: http://dx.doi.org/10.5772/intechopen.96611*

4.**for** *m* ¼ 1, 2, ⋯, *M* **do**

5. Define *Im* ≔ ½ � *xm*�1, *xm* and make use of the SRM to solve

$$\mathbf{y}^{\prime\prime} = u(\boldsymbol{\kappa}, \boldsymbol{y}, \boldsymbol{y}^{\prime}), \boldsymbol{\kappa} \in I\_m,\\ y(\boldsymbol{\kappa}\_{m-1}) = a\_0, \boldsymbol{y}^{\prime}(\boldsymbol{\kappa}\_{m-1}) = a\_0^{\prime}$$

for the estimate ~*y* to *y* on *Im*.

6. Set *α*<sup>0</sup> ¼ ~*y x*ð Þ *<sup>m</sup>* and set *α*<sup>0</sup> <sup>0</sup> ¼ ~*y*<sup>0</sup> ð Þ *xm* .

7.**end**

diag½ �¼ *b*0, ⋯, *bN*

is a diagonal matrix, and **<sup>R</sup>**ð Þ*<sup>i</sup>* <sup>¼</sup> *<sup>R</sup>*ð Þ*<sup>i</sup> <sup>η</sup>*<sup>0</sup> ð Þ, <sup>⋯</sup>*R*ð Þ*<sup>i</sup>* ð Þ *<sup>η</sup><sup>N</sup>*

**v***<sup>r</sup>* ¼ *vr η*<sup>0</sup> ½ � ð Þ, ⋯, *vr*ð Þ *η<sup>N</sup>*

*Recent Advances in Numerical Simulations*

Moreover, we assume that

The initial conditions

below.

problem (5).

**120**

to Eq. (7) on *Iε*.

3. Set *α*<sup>0</sup> ¼ *yε*ð Þ*ε* , set *α*<sup>0</sup>

and ½ �*<sup>T</sup>* denotes matrix transpose.

We prescribe *v*<sup>0</sup> and *w*<sup>0</sup> by requiring that

*A*1 0 ⋯ 0 1 � �

*A*2 0⋯0 1 � �

**Algorithm 1**: Putting the MSRM together.

<sup>0</sup> ¼ *y*<sup>0</sup>

*b*0

0

B@

*<sup>T</sup>*, **<sup>w</sup>***<sup>r</sup>* <sup>¼</sup> *wr <sup>η</sup>*<sup>0</sup> ½ � ð Þ, <sup>⋯</sup>, *wr*ð Þ *<sup>η</sup><sup>N</sup>*

lim*<sup>r</sup>*!<sup>∞</sup>*vr* <sup>¼</sup> *<sup>v</sup>* and lim*<sup>r</sup>*!<sup>∞</sup> *wr* <sup>¼</sup> *<sup>w</sup>*

*vr*þ<sup>1</sup>ð Þ¼ *η<sup>N</sup> α*<sup>0</sup> and *wr*þ<sup>1</sup>ð Þ¼ *η<sup>N</sup> α*<sup>0</sup>

are included in the iterative scheme consisting of Eqs. (19) and (20) as shown

1

1

The larger *L* is the less reliable is the SRM. As a workaround to this problem is we

1.Let 0, ½ �¼ *L Iε*∪½ � *ε*, *L* where *I<sup>ε</sup>* ≔ ½ � 0, *ε* for some sufficiently small number *ε*>0.

2.Replace nonlinear Eq. (5) with linear Eq. (8) on *I<sup>ε</sup>* and compute solution *yε*ð Þ *x*

*<sup>ε</sup>*ð Þ*ε* and let *ε* ¼ *x*<sup>0</sup> < *x*<sup>1</sup> < ⋯ <*xM* ¼ *L*.

CCA ¼

CCA ¼

0

BBB@

0

BBB@

*vr*þ<sup>1</sup> *η*<sup>0</sup> ð Þ ⋮ *vr*þ<sup>1</sup> *η<sup>N</sup>*�<sup>1</sup> ð Þ *vr*þ<sup>1</sup>ð Þ *η<sup>N</sup>*

*wr*þ<sup>1</sup> *η*<sup>0</sup> ð Þ ⋮ *wr*þ<sup>1</sup> *η<sup>N</sup>*�<sup>1</sup> ð Þ *wr*þ<sup>1</sup>ð Þ *η<sup>N</sup>*

subdivide interval ½ � *ε*, *L* into a disjoint union of non-overlapping subintervals as detailed in Eq. (6). Given the the model problem on *I*1, we use the SRM to compute estimate ~*y* to *y* on *I*1. We make use of ~*y* to compute initial values for the problem on *I*2. We repeat this procedure for the problem on *I*<sup>2</sup> and continue in a similar manner until we exhaust ½ � *ε*, *L* . Shown in Algorithm 2.2 is an outline of the MSRM for

0

BB@

0

BB@

*v*0ð Þ� *ε y*0ð Þ¼ *ε α*<sup>0</sup> and*w*0ð Þ� *ε y*<sup>0</sup>

⋱

*bN*

1

� �*<sup>T</sup>* with *<sup>i</sup>* <sup>¼</sup> 1, 2. Moreover,

<sup>0</sup>ð Þ¼ *ε α*<sup>0</sup> 0

0

*<sup>R</sup>*ð Þ<sup>1</sup> *<sup>η</sup>*<sup>0</sup> ð Þ ⋮ *<sup>R</sup>*ð Þ<sup>1</sup> *<sup>η</sup><sup>N</sup>*�<sup>1</sup> ð Þ *α*0

*<sup>R</sup>*ð Þ<sup>2</sup> *<sup>η</sup>*<sup>0</sup> ð Þ ⋮ *<sup>R</sup>*ð Þ<sup>2</sup> *<sup>η</sup><sup>N</sup>*�<sup>1</sup> ð Þ *α*0 0

1

CCCA

1

CCCA

(23)

(24)

CA (21)

*<sup>T</sup>*, *<sup>r</sup>* <sup>¼</sup> 0, 1, <sup>⋯</sup> (22)

#### **2.3 Examples**

In this section we make use of some examples to investigate the accuracy and computational effeciency of the MSRM.

Example 1 We look for a numerical solution to the problem

$$y'' + \frac{1}{x}y' + 3y^5 - y^3 = 0,\\ 0 < x \le 10\\ y(0) = 1,\\ y'(0) = 0 \tag{25}$$

which Wazwaz [11] solved using the VIM and obtained the approximate. analytical solution

$$\mathbf{y}(\mathbf{x}) = \frac{\mathbf{1}}{\sqrt{\mathbf{1} + \mathbf{x}^2}} \tag{26}$$

When we apply the MSRM on (25) we get the following components for constructing the numerical solution.

1. Coefficients

$$u\_0 = -2, L\_0 = 0, H\_0 = -12, G\_0 = -10^4 \tag{27}$$

for problem (8).

2. Initial values

$$y(\varepsilon) \coloneqq a\_0 = 9.999999926 \times 10^{-1} \text{ and } y'(\varepsilon) \coloneqq a\_0' = -1.264241093 \times 10^{-4} \tag{28}$$

for problem (11).

3. Entries

$$R\_i^{(1)} = \boldsymbol{w}(\eta\_i) \text{ and } R\_i^{(2)} = -\left(\mathfrak{d}v\_{r+1}^5(\eta\_i) - v\_{r+1}^3(\eta\_i)\right), i = 0, \cdots, N$$

in the vectors on the right hand sides of Eqs. (19) and (20).

The MSRM generates a numerical solution to problem (25) which is plotted together with the analytical solution (26) in **Figure 1(a)**. **Figure 1(a)** shows a good agreement between the numerical and analytical solutions. For a more detailed comparion of the two solutions see **Table 1**. **Table 1** and **Figure 1(b)** show that the absolute error of the numerical solution is no more than 0*:*<sup>5</sup> � <sup>10</sup>�6. Thus the numerical solution agrees with the analytical solution in the first 6 decimal places.

Applying the MSRM on (29) produces the following building blocks for

*A Modified Spectral Relaxation Method for Some Emden-Fowler Equations*

*<sup>i</sup>* ¼ �8 *e*

decimal places. For the MSRM to achieve this degree of accuracy we chose

*<sup>N</sup>* <sup>¼</sup> 60, *<sup>M</sup>* <sup>¼</sup> 10,*rmax* <sup>¼</sup> 20 and *<sup>ε</sup>* <sup>¼</sup> <sup>10</sup>�6. We observed that the MSRM stopped at

*<sup>x</sup> <sup>y</sup>*<sup>0</sup> � <sup>6</sup>*<sup>y</sup>* � <sup>4</sup>*<sup>y</sup>* ln *<sup>y</sup>* <sup>¼</sup> , 0 <sup>&</sup>lt;*x*≤1, *<sup>y</sup>*ð Þ¼ <sup>0</sup> 1, *<sup>y</sup>*<sup>0</sup>

*y x*ð Þ¼ *e*

is the exact solution [25]. We restrict the problem to a relatively small interval

*x*2

of the vectors on the right hand sides of Eqs. (19) and (20).

*<sup>u</sup>*<sup>0</sup> ¼ �24, *<sup>L</sup>*<sup>0</sup> <sup>¼</sup> 0, *<sup>H</sup>*<sup>0</sup> ¼ �16, *<sup>G</sup>*<sup>0</sup> ¼ �<sup>5</sup> � <sup>10</sup><sup>4</sup> (31)

<sup>0</sup> ¼ �4*:*<sup>767657761</sup> � <sup>10</sup>�<sup>4</sup> (32)

ð Þ¼ 0 0, (33)

(34)

ð Þ*ε* ≔ *α*<sup>0</sup>

*vr*þ<sup>1</sup> *<sup>η</sup><sup>i</sup>* ð Þ <sup>þ</sup> <sup>2</sup>*<sup>e</sup>*

A comparison of the MSRM solution to problem (29) with the analytical solution (30) is shown in **Figure 2(a)**. The graph suggests that the two solutions are exactly the same. However, a closer look at the two solutions is provided in **Table 2** and we observe that the analytical and numerical solutions are slightly different. A plot of the absolute error of the MSRM solution over a grid on the problem domain 0, 10 ½ � is shown in **Figure 2(b)**. We observe that the absolute error does not exceed <sup>0</sup>*:*<sup>5</sup> � <sup>10</sup>�7. Hence the MSRM solution and the analytical solution agree to within 7

*vr*þ<sup>1</sup> *<sup>η</sup><sup>i</sup>* ð Þ*=*<sup>2</sup> , *<sup>i</sup>* <sup>¼</sup> 0, <sup>⋯</sup>, *<sup>N</sup>*

constructing the numerical solution.

*DOI: http://dx.doi.org/10.5772/intechopen.96611*

*<sup>y</sup>*ð Þ*<sup>ε</sup>* <sup>≔</sup> *<sup>α</sup>*<sup>0</sup> ¼ �3*:*<sup>846468396</sup> � <sup>10</sup>�<sup>8</sup> and *<sup>y</sup>*<sup>0</sup>

*<sup>i</sup>* <sup>¼</sup> *<sup>w</sup> <sup>η</sup><sup>i</sup>* ð Þ and *<sup>R</sup>*ð Þ<sup>2</sup>

1. Coefficients

for problem (8).

2. Initial values

for problem (11).

*R*ð Þ<sup>1</sup>

iteration *r* ¼ 5 for problem (29).

*y*<sup>00</sup> þ 2

where

**Figure 2.**

**123**

Example 3 We seek a numerical solution to

½ � 0, 1 because the solution (34) grows rapidly.

*(a) Solutions to problem (29). (b) Absolute error of the MSRM solution to (29).*

3. Elements

**Figure 1.** *(a) Solutions to problem (25). (b) Absolute error of the MSRM solution to (25).*


#### **Table 1.**

*Comparison of numerical values of y with solution (26).*

This degree of accuracy was achieved by the MSRM upon choosing *N* ¼ 60, *M* ¼ 10, setting the maximum number of iterations as *rmax* ¼ 20 and imposing the stopping criterion

$$\max\left\{ \|\mathbf{v}\_{r+1} - \mathbf{v}\_r\|\_{\infty}, \|\mathbf{w}\_{r+1} - \mathbf{w}\_r\|\_{\infty} \right\} \le \varepsilon, r = 0, 1, \dots$$

for the iterative method, where *<sup>ε</sup>* <sup>¼</sup> <sup>10</sup>�<sup>6</sup> is the tolerance and <sup>∥</sup>**w**∥<sup>∞</sup> <sup>¼</sup> max <sup>0</sup> <sup>≤</sup>*i*<sup>≤</sup> *<sup>N</sup>*∣*wi*∣ is the infinity norm. It took only *r* ¼ 5 iterations for the MSRM to achive the given degree of accuracy.

Example 2 We consider the initial value problem

$$y'' + \frac{5}{x}y' + 8\left(e^y + 2e^{y/2}\right) = 0,\\ 0 < x \le 10,\\ y(0) = y'(0) = 0 \tag{29}$$

that was solved by Wazwaz [10] using the ADM to obtain the approximate analytical solution

$$y(\mathbf{x}) = -2\ln\left(\mathbf{1} + \mathbf{x}^2\right) \tag{30}$$

Applying the MSRM on (29) produces the following building blocks for constructing the numerical solution.

1. Coefficients

$$
\mu\_0 = -24, L\_0 = 0, H\_0 = -16, G\_0 = -5 \times 10^4 \tag{31}
$$

for problem (8).

2. Initial values

$$y(\varepsilon) \coloneqq a\_0 = -3.846468396 \times 10^{-8} \text{ and } y'(\varepsilon) \coloneqq a\_0' = -4.767657761 \times 10^{-4} \tag{32}$$

for problem (11).

3. Elements

$$R\_i^{(1)} = w(\eta\_i) \text{ and } R\_i^{(2)} = -8 \left( e^{v\_{r+1}(\eta\_i)} + 2e^{v\_{r+1}(\eta\_i)/2} \right), i = \mathbf{0}, \dots, N$$

of the vectors on the right hand sides of Eqs. (19) and (20).

A comparison of the MSRM solution to problem (29) with the analytical solution (30) is shown in **Figure 2(a)**. The graph suggests that the two solutions are exactly the same. However, a closer look at the two solutions is provided in **Table 2** and we observe that the analytical and numerical solutions are slightly different. A plot of the absolute error of the MSRM solution over a grid on the problem domain 0, 10 ½ � is shown in **Figure 2(b)**. We observe that the absolute error does not exceed <sup>0</sup>*:*<sup>5</sup> � <sup>10</sup>�7. Hence the MSRM solution and the analytical solution agree to within 7 decimal places. For the MSRM to achieve this degree of accuracy we chose *<sup>N</sup>* <sup>¼</sup> 60, *<sup>M</sup>* <sup>¼</sup> 10,*rmax* <sup>¼</sup> 20 and *<sup>ε</sup>* <sup>¼</sup> <sup>10</sup>�6. We observed that the MSRM stopped at iteration *r* ¼ 5 for problem (29).

Example 3 We seek a numerical solution to

$$y'' + \frac{2}{x}y' - 6y - 4y \ln y = \text{, } 0 < x \le 1, y(0) = 1, y'(0) = 0,\tag{33}$$

where

This degree of accuracy was achieved by the MSRM upon choosing *N* ¼ 60, *M* ¼ 10, setting the maximum number of iterations as *rmax* ¼ 20 and imposing the stopping

*x* **MSRM Solution (26) Absolute error** 1 0*:*<sup>70710678</sup> <sup>0</sup>*:*<sup>70710678</sup> <sup>3</sup>*:*<sup>1</sup> � <sup>10</sup>�<sup>9</sup> 2 0*:*<sup>44721363</sup> <sup>0</sup>*:*<sup>44721360</sup> <sup>3</sup> � <sup>10</sup>�<sup>8</sup> 3 0*:*<sup>31622782</sup> <sup>0</sup>*:*<sup>31622777</sup> <sup>5</sup> � <sup>10</sup>�<sup>8</sup> 4 0*:*<sup>24253570</sup> <sup>0</sup>*:*<sup>24253563</sup> <sup>7</sup>*:*<sup>1</sup> � <sup>10</sup>�<sup>8</sup> 5 0*:*<sup>19611623</sup> <sup>0</sup>*:*<sup>19611614</sup> <sup>9</sup>*:*<sup>4</sup> � <sup>10</sup>�<sup>8</sup> 6 0*:*<sup>16439911</sup> <sup>0</sup>*:*<sup>16439899</sup> <sup>1</sup>*:*<sup>2</sup> � <sup>10</sup>�<sup>7</sup> 7 0*:*<sup>14142151</sup> <sup>0</sup>*:*<sup>14142136</sup> <sup>1</sup>*:*<sup>5</sup> � <sup>10</sup>�<sup>7</sup> 8 0*:*<sup>12403492</sup> <sup>0</sup>*:*<sup>12403473</sup> <sup>1</sup>*:*<sup>9</sup> � <sup>10</sup>�<sup>7</sup> 9 0*:*<sup>11043175</sup> <sup>0</sup>*:*<sup>11043153</sup> <sup>2</sup>*:*<sup>3</sup> � <sup>10</sup>�<sup>7</sup> <sup>10</sup> <sup>0</sup>*:*<sup>09950399</sup> <sup>0</sup>*:*<sup>09950372</sup> <sup>2</sup>*:*<sup>7</sup> � <sup>10</sup>�<sup>7</sup>

max f g ∥**v***<sup>r</sup>*þ<sup>1</sup> � **v***r*∥∞, ∥**w***<sup>r</sup>*þ<sup>1</sup> � **w***r*∥<sup>∞</sup> ≤*ε*,*r* ¼ 0, 1, ⋯

that was solved by Wazwaz [10] using the ADM to obtain the approximate

¼ 0, 0 <*x*≤10, *y*ð Þ¼ 0 *y*<sup>0</sup>

*y x*ð Þ¼�2 ln 1 <sup>þ</sup> *<sup>x</sup>*<sup>2</sup> (30)

ð Þ¼ 0 0 (29)

for the iterative method, where *<sup>ε</sup>* <sup>¼</sup> <sup>10</sup>�<sup>6</sup> is the tolerance and <sup>∥</sup>**w**∥<sup>∞</sup> <sup>¼</sup> max <sup>0</sup> <sup>≤</sup>*i*<sup>≤</sup> *<sup>N</sup>*∣*wi*∣ is the infinity norm. It took only *r* ¼ 5 iterations for the MSRM to

criterion

**Table 1.**

**Figure 1.**

achive the given degree of accuracy.

*y*<sup>00</sup> þ 5

analytical solution

**122**

Example 2 We consider the initial value problem

*<sup>y</sup>* <sup>þ</sup> <sup>2</sup>*<sup>e</sup> <sup>y</sup>=*<sup>2</sup>

*(a) Solutions to problem (25). (b) Absolute error of the MSRM solution to (25).*

*Recent Advances in Numerical Simulations*

*<sup>x</sup> <sup>y</sup>*<sup>0</sup> <sup>þ</sup> <sup>8</sup> *<sup>e</sup>*

*Comparison of numerical values of y with solution (26).*

$$\mathcal{Y}(\mathfrak{x}) = \mathfrak{e}^{\mathfrak{x}^{\mathfrak{x}}} \tag{34}$$

is the exact solution [25]. We restrict the problem to a relatively small interval ½ � 0, 1 because the solution (34) grows rapidly.

**Figure 2.** *(a) Solutions to problem (29). (b) Absolute error of the MSRM solution to (29).*


#### **Table 2.**

*Comparison of numerical values of y with solution (30).*

The MSRM for (34) gives the following components for constructing the numerical solution.

1. Coefficients

$$
\mu\_0 = \mathfrak{G}, L\_0 = \mathbf{0}, H\_0 = \mathbf{10} \text{ and } G\_0 = -2 \times \mathbf{10}^4 \tag{35}
$$

**Figure 3(a)** shows that the numerical solution agrees well with the analytical solution (34). For a closer look at how the numerical and analytical solutions compare see **Table 3**. We observe that the MSRM solution agrees with the analytical solution (34) on the problem domain 0, 1 ½ � to within at least 7 decimal places. A plot of the absolute error at these grid points on 0, 1 ½ � is shown in **Figure 3(b)**. We observe that the absolute error in the MSRM increases as me move away from the

*x* **MSRM Solution (34) Absolute error** *:*1 1*:*<sup>01005018</sup> <sup>1</sup>*:*<sup>01005017</sup> <sup>1</sup>*:*<sup>4</sup> � <sup>10</sup>�<sup>8</sup> *:*2 1*:*<sup>04081079</sup> <sup>1</sup>*:*<sup>04081077</sup> <sup>1</sup>*:*<sup>4</sup> � <sup>10</sup>�<sup>8</sup> *:*3 1*:*<sup>09417430</sup> <sup>1</sup>*:*<sup>09417428</sup> <sup>1</sup>*:*<sup>5</sup> � <sup>10</sup>�<sup>8</sup> *:*4 1*:*<sup>17351089</sup> <sup>1</sup>*:*<sup>17351087</sup> <sup>1</sup>*:*<sup>6</sup> � <sup>10</sup>�<sup>8</sup> *:*5 1*:*<sup>28402543</sup> <sup>1</sup>*:*<sup>28402542</sup> <sup>1</sup>*:*<sup>8</sup> � <sup>10</sup>�<sup>8</sup> *:*6 1*:*<sup>43332943</sup> <sup>1</sup>*:*<sup>43332941</sup> <sup>2</sup> � <sup>10</sup>�<sup>8</sup> *:*7 1*:*<sup>63231624</sup> <sup>1</sup>*:*<sup>63231622</sup> <sup>2</sup>*:*<sup>2</sup> � <sup>10</sup>�<sup>8</sup> *:*8 1*:*<sup>89648091</sup> <sup>1</sup>*:*<sup>89648088</sup> <sup>2</sup>*:*<sup>8</sup> � <sup>10</sup>�<sup>8</sup> *:*9 2*:*<sup>24790802</sup> <sup>2</sup>*:*<sup>24790799</sup> <sup>3</sup>*:*<sup>6</sup> � <sup>10</sup>�<sup>8</sup> 1 2*:*<sup>71828187</sup> <sup>2</sup>*:*<sup>71828183</sup> <sup>4</sup>*:*<sup>6</sup> � <sup>10</sup>�<sup>8</sup>

*A Modified Spectral Relaxation Method for Some Emden-Fowler Equations*

accuracy is achieved with *<sup>N</sup>* <sup>¼</sup> 60, *<sup>M</sup>* <sup>¼</sup> 10,*rmax* <sup>¼</sup> 20 and *<sup>ε</sup>* <sup>¼</sup> <sup>10</sup>�6. Moreover, only

In this chapter we presented a modified spectral relaxation method (denoted by MSRM) for solving singular initial value problems for some Emden-Fowler equations. We made use of some examples of the model problem to demonstrate that the MSRM is reliable, accurate and computationally efficient. The method provided a reliable treatment of the singular point. The MSRM solutions were compared with analytical solutions that were obtained using other methods, i.e., the Variational iteration method and the Adomian decomposition method. There was agreement between the solutions that were compared in the first 6 decimal places. A possible way of increasing the degree of accuracy of the MSRM would be to increase the tolerance for the method. This and other ways for optimizing the method could consistute future work. In all the examples that were considered, it took at most 6 iterations for the

. This degree of

singular point *<sup>x</sup>* <sup>¼</sup> 0, but the absolute error never exceeds 0*:*<sup>5</sup> � <sup>10</sup>�<sup>7</sup>

MSRM to converge. Hence the method exhibited rapid convergence.

6 iterations were required to achieve this degree of accuracy.

*Comparison of numerical values of y with solution (34).*

*DOI: http://dx.doi.org/10.5772/intechopen.96611*

**3. Conclusions**

**125**

**Table 3.**

for problem (8).

2. Initial values

$$a\_0 = 1.000000017 \text{ and } a\_0' = 2.593994191 \times 10^{-04} \tag{36}$$

for problem (11).

3. Right hand sides

$$R\_i^{(1)} = w\_r(\eta\_i) \text{ and} \\ R\_i^{(2)} = 6v\_{r+1}(\eta\_i) + 4v\_{r+1}(\eta\_i) \ln \left( v\_{r+1}(\eta\_i) \right), i = 0, \dots, N \tag{37}$$

for linear systems (19) and (20).

**Figure 3.** *(a) Solutions to problem (33). (b) Absolute error of MSRM solution to (33).*


*A Modified Spectral Relaxation Method for Some Emden-Fowler Equations DOI: http://dx.doi.org/10.5772/intechopen.96611*

#### **Table 3.**

The MSRM for (34) gives the following components for constructing the

*x* **MSRM Solution (34) Absolute error** �1*:*<sup>38629437</sup> �1*:*<sup>38629436</sup> <sup>9</sup>*:*<sup>9</sup> � <sup>10</sup>�<sup>9</sup> �3*:*<sup>21887583</sup> �3*:*<sup>21887582</sup> <sup>1</sup>*:*<sup>6</sup> � <sup>10</sup>�<sup>9</sup> �4*:*<sup>60517018</sup> �4*:*<sup>60517019</sup> <sup>1</sup>*:*<sup>4</sup> � <sup>10</sup>�<sup>9</sup> �5*:*<sup>66642669</sup> �5*:*<sup>66642669</sup> <sup>1</sup>*:*<sup>4</sup> � <sup>10</sup>�<sup>9</sup> �6*:*<sup>51619308</sup> �6*:*<sup>51619308</sup> <sup>1</sup>*:*<sup>4</sup> � <sup>10</sup>�<sup>10</sup> �7*:*<sup>22183583</sup> �7*:*<sup>22183583</sup> <sup>2</sup>*:*<sup>3</sup> � <sup>10</sup>�<sup>10</sup> �7*:*<sup>82404601</sup> �7*:*<sup>82404601</sup> <sup>3</sup>*:*<sup>4</sup> � <sup>10</sup>�<sup>10</sup> �8*:*<sup>34877454</sup> �8*:*<sup>34877454</sup> <sup>3</sup>*:*<sup>4</sup> � <sup>10</sup>�<sup>10</sup> �8*:*<sup>81343849</sup> �8*:*<sup>81343849</sup> <sup>1</sup>*:*<sup>2</sup> � <sup>10</sup>�<sup>9</sup> �9*:*<sup>23024103</sup> �9*:*<sup>23024103</sup> <sup>1</sup>*:*<sup>9</sup> � <sup>10</sup>�<sup>9</sup>

*α*<sup>0</sup> ¼ 1*:*000000017 and *α*<sup>0</sup>

*(a) Solutions to problem (33). (b) Absolute error of MSRM solution to (33).*

*<sup>u</sup>*<sup>0</sup> <sup>¼</sup> 6, *<sup>L</sup>*<sup>0</sup> <sup>¼</sup> 0, *<sup>H</sup>*<sup>0</sup> <sup>¼</sup> 10 and *<sup>G</sup>*<sup>0</sup> ¼ �<sup>2</sup> � <sup>10</sup><sup>4</sup> (35)

*<sup>i</sup>* ¼ 6*vr*þ<sup>1</sup> *η<sup>i</sup>* ð Þþ 4*vr*þ<sup>1</sup> *η<sup>i</sup>* ð Þln *vr*þ<sup>1</sup> *η<sup>i</sup>* ð Þ ð Þ , *i* ¼ 0, ⋯, *N* (37)

<sup>0</sup> <sup>¼</sup> <sup>2</sup>*:*<sup>593994191</sup> � <sup>10</sup>�<sup>04</sup> (36)

numerical solution.

**Table 2.**

1. Coefficients

for problem (8).

2. Initial values

for problem (11).

3. Right hand sides

*<sup>i</sup>* <sup>¼</sup> *wr <sup>η</sup><sup>i</sup>* ð Þand*R*ð Þ<sup>2</sup>

for linear systems (19) and (20).

*Comparison of numerical values of y with solution (30).*

*Recent Advances in Numerical Simulations*

*R*ð Þ<sup>1</sup>

**Figure 3.**

**124**

*Comparison of numerical values of y with solution (34).*

**Figure 3(a)** shows that the numerical solution agrees well with the analytical solution (34). For a closer look at how the numerical and analytical solutions compare see **Table 3**. We observe that the MSRM solution agrees with the analytical solution (34) on the problem domain 0, 1 ½ � to within at least 7 decimal places. A plot of the absolute error at these grid points on 0, 1 ½ � is shown in **Figure 3(b)**. We observe that the absolute error in the MSRM increases as me move away from the singular point *<sup>x</sup>* <sup>¼</sup> 0, but the absolute error never exceeds 0*:*<sup>5</sup> � <sup>10</sup>�<sup>7</sup> . This degree of accuracy is achieved with *<sup>N</sup>* <sup>¼</sup> 60, *<sup>M</sup>* <sup>¼</sup> 10,*rmax* <sup>¼</sup> 20 and *<sup>ε</sup>* <sup>¼</sup> <sup>10</sup>�6. Moreover, only 6 iterations were required to achieve this degree of accuracy.

### **3. Conclusions**

In this chapter we presented a modified spectral relaxation method (denoted by MSRM) for solving singular initial value problems for some Emden-Fowler equations. We made use of some examples of the model problem to demonstrate that the MSRM is reliable, accurate and computationally efficient. The method provided a reliable treatment of the singular point. The MSRM solutions were compared with analytical solutions that were obtained using other methods, i.e., the Variational iteration method and the Adomian decomposition method. There was agreement between the solutions that were compared in the first 6 decimal places. A possible way of increasing the degree of accuracy of the MSRM would be to increase the tolerance for the method. This and other ways for optimizing the method could consistute future work. In all the examples that were considered, it took at most 6 iterations for the MSRM to converge. Hence the method exhibited rapid convergence.

*Recent Advances in Numerical Simulations*

**References**

492-497.

[1] Van Gorder R.A., Vajravelu K., Analytic and numerical solutions to the

*DOI: http://dx.doi.org/10.5772/intechopen.96611*

*A Modified Spectral Relaxation Method for Some Emden-Fowler Equations*

[10] Wazwaz A.M., Adomian decomposition method for a reliable treatment of the Emden-Fowler equation, Applied Mathematics and Computation 161 (2005) 543-560.

217 (2011) 10387-10395.

34 (2015) 178-186.

(2018) 9, 615-629.

[11] Wazwaz A.M., A reliable treatment of singular Emden-Fowler initial value problems and boundary value problems, Applied Mathematics and Computation

[12] Kazemi Nasab A., Kılıçman A., Pashazadeh Atabakan Z., Leong W.J., A numerical approach for solving singular nonlinear Lane-Emden type equations arising in astrophysics, New Astronomy

[13] Singh R., Garg H., Guleria V., Haar wavelet collocation method for Lane-Emden equations with Dirichlet, Neumann and Neumann-Robin boundary conditions, Journal of Computational and Applied Mathematics 346 (2019) 150-161.

[14] Gümgüm S., Taylor wavelet solution of linear and nonlinear Lane-Emden equations, Applied Numerical Mathematics 158 (2020) 44-53.

[15] Parand K., Hashemi S., RBF-DQ method for solving non-linear

differential equations of Lane-Emden type, Ain Shams Engineering Journal

[16] Iserles A., A first course in the numerical analysis of differential euqations, Cambridge, 2009.

[17] Motsa S.S., A new spectral

Volume 201, 2014 - Issue 2.

[18] Motsa S.S., Dlamini P.G., Khumalo M., Spectral relaxation method and spectral quasilinearization

relaxation method for similarity variable nonlinear boundary layer flow systems, Chemical Engineering Communications

[2] Šmarda Z., Khan Y., An efficient computational approach to solving singular initial value problems for Lane-Emden type equations, Journal of Computational and Applied Mathematics 290 (2015) 65-73.

[3] Van Gorder R.A., Analytic and numerical solutions to the Lane-Emden equation, New Astronomy 16 (2011)

[5] Căruntu B., Bota C., Approximate polynomial solutions of the nonlinear Lane-Emden type equations arising in

[6] Răb M., Bounds for solutions of the

Math.2, Scripta Fac.Sci.Nat.UJEP Brunessis, X1, 1975, 79-84.

[7] Torres-Córdoba R., Martínez-García E.A., Exact analytic solution of an unsolvable class of first Lane-Emden equation for polytropic gas sphere, New

Astronomy 82 (2021) 101458.

[8] Van Gorder R.A., Relation between Lane-Emden solutions and radial solutions to the elliptic Heavenly equation on a disk, New Astronomy 37

[9] Chowdhury M.S.H., Hashim I., Solutions of Emden-Fowler equations by homotopy-perturbation method, Nonlinear analysis: real world applications 10 (2009) 104-115.

þ *qu* ¼ *h x*, *u*, *u*<sup>0</sup> ð Þ, Arch.

astrophysics using the squared remainder minimization method, Computer Physics Communications 184

[4] Momoniat E., Harley C., Approximate implicit solution of a Lane-Emden equation, New Astronomy

11 (2006) 520-526.

(2013) 1643-1648.

equation *pu*<sup>0</sup> ð Þ0

(2015) 42-47.

**127**

Lane–Emden equation, Physics Letters A, 372 (2008) 6060-6065.
