Open access peer-reviewed chapter

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

Written By

Gerald Tendayi Marewo

Submitted: December 14th, 2020 Reviewed: February 13th, 2021 Published: April 27th, 2021

DOI: 10.5772/intechopen.96611

From the Edited Volume

## Recent Advances in Numerical Simulations

Edited by Francisco Bulnes and Jan Peter Hessling

Chapter metrics overview

View Full Metrics

## Abstract

In this chapter, we present a modified version of the spectral relaxation method for solving singular initial value problems for some Emden-Fowler equations. This study was motivated by the several applications that these equations have in Science. The first step of the method of solution makes use of linearisation to solve the model problem on a small subinterval of the problem domain. This subinterval contains a singularity at the initial instant. The first step is combined with using the spectral relaxation method to recursively solve the model problem on the rest of the problem domain. We make use of examples to demonstrate that the method is reliable, accurate and computationally efficient. The numerical solutions that are obtained in this chapter are in good agreement with other solutions in the literature.

### Keywords

• Emden-Fowler equations
• Lane-Emden equations
• singular initial value problem
• spectral relaxation method
• numerical method

## 1. Introduction

The singular initial value problem

d2ydx2+γxdydx+rxy=sx,0<x,γ>0y0=α,dydx0=0,αRE1

for the Lane-Emden Eq. (1) models several phenomena such as the thermal behaviour of a spherical cloud of gas acting under the mutual attraction of its molecules [1], the temperature variation of a self gravitating star, the kinetics of combustion [2], thermal explosion in a rectangular slab [3] and the density distribution in isothermal gas spheres [4]. Moreover, Eq. (1) has been used many a time as a benchmark for new methods.

A particular case of Eq. (1) is the Emden-Fowler equation of the first kind:

d2ydx2+2xdydx+ym=0,y0=1,dydx0=0,mNE2

As mentioned in [5], Eq. (2) represents the dimensionless form of the governing equation for the gravitational potential of a Newtonian self-gravitating, spherically symmetric, polytropic fluid. The equation provides a useful approximation for stars.

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

d2ydx2+γxdydx+fxgy=0E3

which can be written as

py+qy=hxyy where=ddx,E4

an equation which was discussed in [6]. An existence result for the solution is given therein under certain conditions on px,qx and hxyy.

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 obtained demonstrated the reliability and effectiveness of the VIM.

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 use of the collocation points

xj=j0.5/2M,j=1,,2M

on the problem domain 0L, 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

xj=2L1cosj1N1,j=1,,N

on the problem domain 0L and the method makes use of Radial basis functions. Unlike making use of collocation methods for solving Eq. (3) Van Gorder and Vajravelu [1] used the Runge–Kutta-Felhberg 4-5 (RKF45) method to validate the analytical solutions that they obtained from using the HAM and from using the traditional power series method. The RKF45 method is an embedded Runge–Kutta-pair which makes use of an adaptive stepsize to control the method and to ensure stability properties such as A-stability. See [16] for more details on the RKF45 method.

In this chapter we make use of a modified version of the spectral relation method (SRM) to solve an initial value problem for Eq. (3). We denote our method by MSRM. The SRM was successfully used to solve fluid flow problems by for example Motsa [17], Motsa et al. [18], Shateyi et al. [19] and Gangadhar et al. [20] to mention a few. In [17, 18, 19, 20] the SRM was shown to be accurate, computationally efficient and easy to implement. Moreover, the SRM was applied only after transforming the governing partial differential differential equations to ordinary differential equations. However, not so long ago the SRM was modified in such a way that it was directly applicable to partial differential equations. See for example [21]. The SRM was used to solve other types of problems such as hyperchaotic systems [22]. It is to the best of our knowledge that the MSRM has not been used in existing literature. We chose the MSRM because it is not computationally intensive and it is easy to implement.

In Section 2 we describe the MSRM for the model problem. In Section 3 we make use of examples to demonstrate the accuracy and computational efficiency of the MSRM. Section 4 concludes this chapter.

## 2. The MSRM for the model problem

We seek an approximate solution to

y+γxy+fxgy=0,0<xL,y0=α,y0=0E5

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ε=0ε of I=0L where ε>0. The point x=ε splits interval I into two subintervals: Iε and IIε=ε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 IIε we form a partition

IIε=m=1MImwhereIm=xm1xm,x0=εandxM=LE6

The method by Ramos proceeds by using the same linearisation method on the subintervals Im,m=1M of IIε, i.e. away from the singularity. To this end, at the interface xm1 of Im1 and Im we make use the solution of Eq. (3) restricted to Im1 to generate the initial conditions for Eq. (3) restricted to Im. This ensures continuity of the solution. In this chapter we avoid linearisation on IIε by making use of the SRM on the subintervals of IIε. 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 ε0L be a sufficiently small number. Restrict problem (5) to 0ε and re-arrange to get

y=γxy+fxgyuxyy,0<xε,y0=α,y0=0E7

If we Taylor expand u about ε0=εy0y0y0y0 and neglect higher order terms we get

uxyy=uε0u0+uyε0H0yy0+uyε0G0yy0+uxε0L0xε

where us denotes us. Consequently, Eq. (7) can be replaced by

y=u0+H0yy0+G0yy0+L0xε,0xε,y0=α,y0=0,E8

where the differential equation now holds at x=0 because the singularity there is no longer present. If H00 and R0G0/22+H0>0, then problem (8) has analytical solution [23]

yx=A0expλ0+xε+B0expλ0xε+C0xε+D0E9

for 0xε,whereλ0±=G0/2±R0,C0=L0/H0,

D0=G0C0+P0/H0,P0=u0H0y0G0y0,A0=expλ0+ελ0+λ0y0C0λ0y0+C0εD0,B0=expλ0ελ0+λ0λ0+y0+C0εD0y0C0

Thus

yε=A0+B0+D0α0andyε=λ0+A0+λ0B0+C0α0E10

We take α0 and α0 as initial values for problem (7) restricted to εL in the next section.

### 2.2 Away from the singularity

We seek yx satisfying

y+γxy+fxgy=0,εxL,yε=α0,yε=α0E11

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.

The first step of the SRM for (11) is to let v=y and w=v so that we obtain the equivalent problem

v=w,vε=α0,E12
w+γxw+fxgv=0,wε=α0E13

for εxL which upon making use of the change of variable

xη=L+ε2+Lε2η

becomes

βdv=w,v1=α0E14
βdw+γxηw+fxηgv=0,w1=α0E15

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

βdvr+1=wrR1η,vr+1ηN=α0E16
βdwr+1+γxwr+1=fxηgvr+1R2η,wr+1ηN=α0E17

for 1=ηNηη0=1 where r=0,1, and on 11 we formed a grid consisting of the Chebyshev-Gauss-Lobatto collocation points

ηi=cosπiN,i=0,1,,N.E18

If the initial approximations v0 and w0 to v and w, respectively, are prescribed then Eqs. (16) and (17) generate sequences vrr=1 and wrr=1 of consecutive approximations. To this end we assume that vr and wr are known at the end of the rth iteration. Once Eq. (16) is solved for vr+1 the right hand side R2 of Eq. (17) becomes known and we solve this equation for wr+1. Since vr+1 and wr+1 are now known, we proceed in a similar manner to compute vr+2 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

D̂A1vr+1=R1,vr+1ηN=α0,E19
D̂+γdiag1/xη01/xηNA2wr+1=R2,wr+1ηN=α0E20

where D̂=βD,D is the N×N Chebyshev differentiation matrix,

diagb0bN=b0bNE21

is a diagonal matrix, and Ri=Riη0RiηNT with i=1,2. Moreover,

vr=vrη0vrηNT,wr=wrη0wrηNT,r=0,1,E22

and T denotes matrix transpose.

We prescribe v0 and w0 by requiring that

v0εy0ε=α0andw0εy0ε=α0

Moreover, we assume that

limrvr=vandlimrwr=w

The initial conditions

vr+1ηN=α0andwr+1ηN=α0

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

A1001vr+1η0vr+1ηN1vr+1ηN=R1η0R1ηN1α0E23
A2001wr+1η0wr+1ηN1wr+1ηN=R2η0R2ηN1α0E24

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 I1, we use the SRM to compute estimate y˜ to y on I1. We make use of y˜ to compute initial values for the problem on I2. We repeat this procedure for the problem on I2 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.

1. Let 0L=IεεL where Iε0ε for some sufficiently small number ε>0.

2. Replace nonlinear Eq. (5) with linear Eq. (8) on Iε and compute solution yεx to Eq. (7) on Iε.

3. Set α0=yεε, set α0=yεε and let ε=x0<x1<<xM=L.

4. form=1,2,,Mdo

5.  Define Imxm1xm and make use of the SRM to solve

y=uxyy,xIm,yxm1=α0,yxm1=α0

for the estimate y˜ to y on Im.

6.  Set α0=y˜xm and set α0=y˜'xm.

7. end

### 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+1xy+3y5y3=0,0<x10y0=1,y0=0E25

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

analytical solution

yx=11+x2E26

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

1. Coefficients

u0=2,L0=0,H0=12,G0=104E27

for problem (8).

2. Initial values

yεα0=9.999999926×101andyεα0=1.264241093×104E28

for problem (11).

3. Entries

Ri1=wηiandRi2=3vr+15ηivr+13ηi,i=0,,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.5×106. Thus the numerical solution agrees with the analytical solution in the first 6 decimal places. 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

xMSRMSolution (26)Absolute error
10.707106780.707106783.1×109
20.447213630.447213603×108
30.316227820.316227775×108
40.242535700.242535637.1×108
50.196116230.196116149.4×108
60.164399110.164398991.2×107
70.141421510.141421361.5×107
80.124034920.124034731.9×107
90.110431750.110431532.3×107
100.099503990.099503722.7×107

### Table 1.

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

maxvr+1vrwr+1wrε,r=0,1,

for the iterative method, where ε=106 is the tolerance and w=max0iNwi 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+5xy+8ey+2ey/2=0,0<x10,y0=y0=0E29

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

yx=2ln1+x2E30

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

1. Coefficients

u0=24,L0=0,H0=16,G0=5×104E31

for problem (8).

2. Initial values

yεα0=3.846468396×108andyεα0=4.767657761×104E32

for problem (11).

3. Elements

Ri1=wηiandRi2=8evr+1ηi+2evr+1ηi/2,i=0,,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 010 is shown in Figure 2(b). We observe that the absolute error does not exceed 0.5×107. 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 N=60,M=10,rmax=20 and ε=106. We observed that the MSRM stopped at iteration r=5 for problem (29).

xMSRMSolution (34)Absolute error
11.386294371.386294369.9×109
23.218875833.218875821.6×109
34.605170184.605170191.4×109
45.666426695.666426691.4×109
56.516193086.516193081.4×1010
67.221835837.221835832.3×1010
77.824046017.824046013.4×1010
88.348774548.348774543.4×1010
98.813438498.813438491.2×109
109.230241039.230241031.9×109

### Table 2.

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

Example 3 We seek a numerical solution to

y+2xy6y4ylny=,0<x1,y0=1,y0=0,E33

where

yx=ex2E34

is the exact solution [25]. We restrict the problem to a relatively small interval 01 because the solution (34) grows rapidly.

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

1. Coefficients

u0=6,L0=0,H0=10andG0=2×104E35

for problem (8).

2. Initial values

α0=1.000000017andα0=2.593994191×1004E36

for problem (11).

3. Right hand sides

Ri1=wrηiandRi2=6vr+1ηi+4vr+1ηilnvr+1ηi,i=0,,NE37

for linear systems (19) and (20).

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 01 to within at least 7 decimal places. A plot of the absolute error at these grid points on 01 is shown in Figure 3(b). We observe that the absolute error in the MSRM increases as me move away from the singular point x=0, but the absolute error never exceeds 0.5×107. This degree of accuracy is achieved with N=60,M=10,rmax=20 and ε=106. Moreover, only 6 iterations were required to achieve this degree of accuracy.

xMSRMSolution (34)Absolute error
0.11.010050181.010050171.4×108
0.21.040810791.040810771.4×108
0.31.094174301.094174281.5×108
0.41.173510891.173510871.6×108
0.51.284025431.284025421.8×108
0.61.433329431.433329412×108
0.71.632316241.632316222.2×108
0.81.896480911.896480882.8×108
0.92.247908022.247907993.6×108
12.718281872.718281834.6×108

### Table 3.

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

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

## References

1. 1. Van Gorder R.A., Vajravelu K., Analytic and numerical solutions to the Lane–Emden equation, Physics Letters A, 372 (2008) 6060-6065.
2. 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. 3. Van Gorder R.A., Analytic and numerical solutions to the Lane-Emden equation, New Astronomy 16 (2011) 492-497.
4. 4. Momoniat E., Harley C., Approximate implicit solution of a Lane-Emden equation, New Astronomy 11 (2006) 520-526.
5. 5. Căruntu B., Bota C., Approximate polynomial solutions of the nonlinear Lane-Emden type equations arising in astrophysics using the squared remainder minimization method, Computer Physics Communications 184 (2013) 1643-1648.
6. 6. Răb M., Bounds for solutions of the equation pu+qu=hxuu, Arch.Math.2, Scripta Fac.Sci.Nat.UJEP Brunessis, X1, 1975, 79-84.
7. 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. 8. Van Gorder R.A., Relation between Lane-Emden solutions and radial solutions to the elliptic Heavenly equation on a disk, New Astronomy 37 (2015) 42-47.
9. 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.
10. 10. Wazwaz A.M., Adomian decomposition method for a reliable treatment of the Emden-Fowler equation, Applied Mathematics and Computation 161 (2005) 543-560.
11. 11. Wazwaz A.M., A reliable treatment of singular Emden-Fowler initial value problems and boundary value problems, Applied Mathematics and Computation 217 (2011) 10387-10395.
12. 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 34 (2015) 178-186.
13. 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. 14. Gümgüm S., Taylor wavelet solution of linear and nonlinear Lane-Emden equations, Applied Numerical Mathematics 158 (2020) 44-53.
15. 15. Parand K., Hashemi S., RBF-DQ method for solving non-linear differential equations of Lane-Emden type, Ain Shams Engineering Journal (2018) 9, 615-629.
16. 16. Iserles A., A first course in the numerical analysis of differential euqations, Cambridge, 2009.
17. 17. Motsa S.S., A new spectral relaxation method for similarity variable nonlinear boundary layer flow systems, Chemical Engineering Communications Volume 201, 2014 - Issue 2.
18. 18. Motsa S.S., Dlamini P.G., Khumalo M., Spectral relaxation method and spectral quasilinearization method for solving unsteady boundary layer flow problems, Advances in Mathematical Physics, vol. 2014, Article ID 341964, 12 pages, 2014.
19. 19. Shateyi S., Marewo G.T., Numerical analysis of unsteady MHD flow near a stagnation point of a two-dimensional porous body with heat and mass transfer, thermal radiation, and chemical reaction, Boundary Value Problems 2014, 2014:218.
20. 20. Gangadhar K., Kannan T., Sakthivel G., DasaradhaRamaiah K., Unsteady free convective boundary layer flow of a nanofluid past a stretching surface using a spectral relaxation method, International journal of ambient energy 2020, vol. 41, no. 6, 609-616. https://doi.org/10.1080/01430750.2018.1472648
21. 21. Motsa S.S., Animasaun I.L., Unsteady Boundary Layer Flow over a Vertical Surface due to Impulsive and Buoyancy in the Presence of Thermal-Diffusion and Diffusion-Thermo using Bivariate Spectral Relaxation Method, Journal of Applied Fluid Mechanics, Vol. 9, No. 5, pp. 2605-2619, 2016.
22. 22. Motsa S.S., Dlamini P.G., Khumalo M., Solving hyperchaotic systems using the spectral relaxation method, Abstract and Applied Analysis Volume 2012, Article ID 203461, 18 pages. https://doi.org/10.1155/2012/203461.
23. 23. Ramos J.I., Linearization techniques for singular initial-value problems of ordinary differential equations, Applied Mathematics and Computation 161 (2005) 525-542.
24. 24. Trefethen L.N., Spectral Methods in MATLAB, SIAM, 2000.
25. 25. Karimi Dizicheh A., Salahshour S., Ahmadian A., Balaeanu D., A novel algorithm based on the Legendre wavelets spectral technique for solving Lane-Emden equations, Applied Numerical Mathematics 153 (2020) 443-456. https://doi.org/10.1155/2014/341964.

Written By

Gerald Tendayi Marewo

Submitted: December 14th, 2020 Reviewed: February 13th, 2021 Published: April 27th, 2021