Approximate Solution of Lane-Emden Type Equations Using Variation of Parameters Method with an Auxiliary Parameter

Show more

1. Introduction

Many problems in mathematical physics and astrophysics can be modeled by the so-called Lane-Emden type equation defined in the form

${y}^{\u2033}\left(x\right)+\frac{\alpha}{x}{y}^{\prime}\left(x\right)+\mathcal{F}\left(x,y\right)=g\left(x\right),\text{\hspace{0.17em}}\alpha x\ge 0,$ (1)

With the initial conditions:

$y\left(0\right)=\mathfrak{a},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}^{\prime}\left(0\right)=0,$ (2)

where “ $\text{'}$ ” is the differentiation with respect to x, $\mathcal{F}\left(x,y\right)$ is a nonlinear function of x and y, $g\left(x\right)$ is the non-homogeneous term and $\alpha ,\mathfrak{a}$ are constants. The closed form solution of the Lane-Emden type Equation (1) is always enabled [1] in the neighborhood of the point $x=0$, which is called a singular point, for the above initial conditions.

Taking $\alpha =2,g\left(x\right)=0,\mathcal{F}\left(x,y\right)={y}^{k}$ and $\mathfrak{a}=1$ in Equations (1) and (2) respectively, we get

${y}^{\u2033}\left(x\right)+\frac{2}{x}{y}^{\prime}\left(x\right)+{y}^{k}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}x\ge 0.$ (3)

subject to initial conditions:

$y\left(0\right)=1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}^{\prime}\left(0\right)=0.$ (4)

Equations (3) and (4) are known as the classical Lane-Emden equation.

Choosing $\alpha =2,g\left(x\right)=0,\mathcal{F}\left(x,y\right)={\left({y}^{2}-C\right)}^{3/2}$ and $\mathfrak{a}=1$ in Equations (1) and (2) respectively, we get the so-called white-dwarf equation that expresses the gravitation potential of the degenerate white-dwarf stars. Isothermal gas spheres [1] are similarly modeled by

${y}^{\u2033}\left(x\right)+\frac{2}{x}{y}^{\prime}\left(x\right)+{\text{e}}^{y\left(x\right)}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}x\ge 0$ (5)

with the initial conditions:

$y\left(0\right)=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}^{\prime}\left(0\right)=0.$ (6)

The parameter k in Equation (3) has physical significance in the range $0\le k\le 5$ so Equation (3) has a closed form solutions for $k=0,1,5$ [2] but may be resorted to numerical solutions for other values of k but the singularity at $x=0$ is still a challenge to the numerical solutions. The semi-analytic solutions can be found by Adomian’s Decomposition Method (ADM) and perturbation techniques and mostly convergent in restricted regions so another approach such as Pade’s method was required to expand these regions [1] [3] [4] . Many semi-analytical methods and numerical methods like Homotopy analysis method [5] [6] [7] , Picard method [8] and spectral methods [9] [10] [10] were used to solve different types of differential equations. Many researches solved the Lane-Emden equation; like in [12] the authors applied Variational Iteration Method to Lane-Emden equation using some transformation. Singh et al. [13] suggested an effective analytic algorithm using Modified Homotopy Analysis Method (MHAM), at which convergence regions could be adjusted without using Pade’s technique. Also another algorithm was proposed in [14] , using an iterative method which is a hybrid of ADM and Variational Iteration Method (VIM). Lastly, Parand et al. [15] proposed an approximation algorithm using Hermite collocation method at which the solution is reduced to the solution of a system of algebraic equations.

In the present paper, we aim effectively to employ Variation of Parameters Method (VPM) coupling with an unknown auxiliary parameter h to solve the homogeneous as well as non-homogeneous Lane-Emden equation. VPM is free from calculation of the so-called Adomian’s polynomials and is widely applicable because of the reliability of it and the reduction in the size of the computational domain. A comparison with the standard VPM is made to show the reliability and efficiency of the applied algorithm which is simple and is accurately approximates the solution over a large domain.

2. Variation of Parameters Method with an Auxiliary Parameter

Variation of Parameters Method (VPM) was first proposed by Ma [16] [17] [18] and it has been thoroughly used by many researchers to solve a wide class of initial and boundary value problems [19] [20] [21] [22] . Ghaneai and Hosseini [23] inserted an auxiliary parameter in variational iteration algorithm to obtain solution of wave-like and heat-like equations and concluded that an auxiliary parameter provides a simple way to control and adjust the convergence region of approximate solution in a large domain. In [24] VIM coupled with an auxiliary parameter are used to predict the multiple solutions of nonlinear boundary value problems with great success. Also an auxiliary parameter is inserted in VPM to solve some equations arising in physics [25] [26] . To convey the basic step of VPM for solving Lane-Emden equation, consider the general nonlinear ordinary differential equation in operator form as follows:

$Ly\left(x\right)+Ry\left(x\right)+Ny\left(x\right)+\mathcal{H}\left(x\right)=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}x\ge 0,$ (9)

with the initial conditions:

$y\left(0\right)=\mathfrak{a},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}^{\prime}\left(0\right)=0,$ (10)

where

$Ly\left(x\right)=Ly\left(x\right)+h\left(Ly\left(x\right)+Ry\left(x\right)+Ny\left(x\right)+\mathcal{H}\left(x\right)\right),$ (11)

or

$Ly\left(x\right)=Ly\left(x\right)+hGy\left(x\right).$ (12)

According to the Variation of Parameters Method [21] , [25] , [26] and [27] , the solution of the homogeneous Equation (11):

$Ly\left(x\right)=0$ is given as

${y}_{H}\left(x\right)={c}_{1}+{c}_{2}x,$ (13)

For the particular solution, the constants ${c}_{1}$ and ${c}_{2}$ in Equation (13) are replaced by the functions $u\left(x\right)$ and $v\left(x\right)$ respectively, thus we have

${y}_{p}\left(x\right)=u\left(x\right){f}_{1}\left(x\right)+v\left(x\right){f}_{2}\left(x\right).$ (14)

where ${f}_{1}\left(x\right)=1$ and ${f}_{2}\left(x\right)=x$. The functions $u\left(x\right)$ and $v\left(x\right)$ can be determined as

$u\left(x\right)={\displaystyle {\int}_{0}^{x}\frac{{W}_{1}}{W}\text{d}\eta}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}v\left(x\right)={\displaystyle {\int}_{0}^{x}\frac{{W}_{1}}{W}\text{d}\eta},$ (15)

where W is the Wronskian of ${f}_{1}\left(x\right)$ and ${f}_{2}\left(x\right)$ i.e.

$W=\left|\begin{array}{cc}{f}_{1}& {f}_{2}\\ {{f}^{\prime}}_{1}& {{f}^{\prime}}_{2}\end{array}\right|=1,$ (16)

${W}_{1}=\left|\begin{array}{cc}0& {f}_{2}\\ Ly\left(x\right)+hGy\left(x\right)& {{f}^{\prime}}_{2}\end{array}\right|=-x\left(Ly\left(x\right)+hGy\left(x\right)\right),$ (17)

and

${W}_{2}=\left|\begin{array}{cc}{f}_{1}& 0\\ {{f}^{\prime}}_{1}& Ly\left(x\right)+hGy\left(x\right)\end{array}\right|=Ly\left(x\right)+hGy\left(x\right).$ (18)

Substituting the values of $W,{W}_{1}$ and ${W}_{2}$ in Equation (15), yields

$u\left(x\right)=-{\displaystyle {\int}_{0}^{x}\eta \left(Ly\left(\eta \right)+hGy\left(\eta \right)\right)\text{d}\eta},\text{\hspace{0.17em}}\text{\hspace{0.17em}}v\left(x\right)={\displaystyle {\int}_{0}^{x}\left(Ly\left(\eta \right)+hGy\left(\eta \right)\right)\text{d}\eta}.$ (19)

Hence, the general solution of Equation (12) is

$y\left(x\right)={y}_{H}\left(x\right)+{y}_{p}\left(x\right)={c}_{1}+{c}_{2}x+{\displaystyle {\int}_{0}^{x}\left(x-\eta \right)\left(Ly\left(\eta \right)+hGy\left(\eta \right)\right)\text{d}\eta}.$ (20)

Applying the initial conditions (10) in Equation (20) and solve for ${c}_{1}$ and ${c}_{2}$, we get ${c}_{1}=\mathfrak{a}$ and ${c}_{2}=0$.

Substituting the values of ${c}_{1}$ and ${c}_{2}$ in Equation (20) yields

$y\left(x\right)=\mathfrak{a}+{\displaystyle {\int}_{0}^{x}\left(x-\eta \right)\left(Ly\left(\eta \right)+h\left(Ly\left(x\right)+Ry\left(x\right)+Ny\left(x\right)+\mathcal{H}\left(x\right)\right)\right)\text{d}\eta}.$ (21)

which can be solved iteratively as

$\{\begin{array}{l}{y}_{o}\left(x\right)=a\\ {y}_{1}\left(x,h\right)={y}_{o}\left(x\right)+{\displaystyle {\int}_{0}^{x}\left(x-\eta \right)(L{y}_{0}\left(\eta \right)+h(L{y}_{0}\left(\eta \right)+R{y}_{0}\left(\eta \right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+N{y}_{0}\left(\eta \right)+H\left(\eta \right)))\text{d}\eta ,\\ {y}_{n+1}\left(x,h\right)={y}_{o}\left(x\right)+{\displaystyle {\int}_{0}^{x}\left(x-\eta \right)(L{y}_{n}\left(\eta ,h\right)+h(L{y}_{n}\left(\eta ,h\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+R{y}_{n}\left(\eta ,h\right)+N{y}_{n}\left(\eta ,h\right)+H\left(\eta \right))\text{d}\eta ,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}n\ge 1.\end{array}$ (22)

Consequently, an exact solution can be obtained when

$y\left(x,h\right)=\underset{n\to \infty}{\mathrm{lim}}{y}_{n}\left(x,h\right).$ (23)

The valid values of h can be determined by means of the so-called h-curves [24] [25] [26] . According to this h-curve, it can be easily determine the valid region of h, which corresponds to the line segment nearly parallel to the horizontal axis. Note that the Variation of Parameters Method [19] [20] [21] [22] provides the following iterative scheme for Equation (9):

$\{\begin{array}{l}{y}_{o}\left(x\right)=a,\\ {y}_{n+1}\left(x\right)={y}_{o}\left(x\right)+{\displaystyle {\int}_{0}^{x}\left(x-\eta \right)\left(-R{y}_{n}\left(\eta \right)-N{y}_{n}\left(\eta \right)-\mathcal{H}\left(\eta \right)\right)\text{d}\eta},n\ge 0.\end{array}$ (24)

3. Approximate Solution of the Lane-Emden Equation

Now we apply Variation of Parameters Method (VPM) coupled with an unknown auxiliary parameter developed in Section 2 for solving Equation (1). The recursive formulas (22) and (24) for this problem is obtained respectively as

$\begin{array}{l}{y}_{n+1}\left(x,h\right)={y}_{0}\left(x\right)+{\displaystyle {\int}_{0}^{x}\left(x-\eta \right)(\frac{{\partial}^{2}{y}_{n}\left(\eta ,h\right)}{\partial {\eta}^{2}}+h(\frac{{\partial}^{2}{y}_{n}\left(\eta ,h\right)}{\partial {\eta}^{2}}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}+\frac{\alpha}{\eta}\left(\frac{\partial {y}_{n}\left(\eta ,h\right)}{\partial \eta}\right)+f\left(\eta ,{y}_{n}\left(\eta ,h\right)\right)-g\left(\eta \right)))\text{d}\eta .\end{array}$

then

$\begin{array}{l}{y}_{n+1}\left(x,h\right)\\ ={y}_{0}\left(x\right)+{\displaystyle {\int}_{0}^{x}\left(x-\eta \right)\left(\left(-\frac{\alpha}{\eta}\left(\frac{\partial {y}_{n}\left(\eta ,h\right)}{\partial \eta}\right)-f\left(\eta ,{y}_{n}\left(\eta ,h\right)\right)+g\left(\eta \right)\right)\right)\text{d}\eta}.\end{array}$ (25)

Example 1

Consider the following form of Lane-Emden equation [1] [2] [3] :

${y}^{\u2033}\left(x\right)+\frac{2}{x}{y}^{\prime}\left(x\right)+{y}^{k}\left(x\right)=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}x\ge 0.$ (26)

where k is a constant. With the initial conditions:

$y\left(0\right)=1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}^{\prime}\left(0\right)=0.$ (27)

Note that the exact solutions for cases of $k=0,1,5$ are known, VPM developed in Section 2 is applied to solve Equations (26) and (27) for some cases of the constant k.

Case (1) for (k = 0):

The closed form solution is: $y\left(x\right)=1-\frac{{x}^{2}}{6}$. According to standard VPM, we have the following from iterative scheme (24):

$\{\begin{array}{l}{y}_{0}\left(x\right)=y\left(0\right)=1,\\ {y}_{n+1}\left(x,h\right)=1+{\displaystyle {\int}_{0}^{x}\left(x-\eta \right)\left(\left(-\frac{2}{\eta}\left(\frac{\partial {y}_{n}\left(\eta ,h\right)}{\partial \eta}\right)-1\right)\right)\text{d}\eta}.\end{array}$ (28)

The result obtained by standard VPM is not valid for large values of x. Now, using the iterative scheme (22) we have:

${y}_{0}\left(x\right)=y\left(0\right)=1,$

${y}_{1}\left(x,h\right)=1+{\displaystyle {\int}_{0}^{x}\left(x-\eta \right)\left(\frac{{\partial}^{2}{y}_{0}\left(\eta \right)}{\partial {\eta}^{2}}+h\left(\frac{{\partial}^{2}{y}_{0}\left(\eta \right)}{\partial {\eta}^{2}}+\frac{2}{\eta}\left(\frac{\partial {y}_{0}\left(\eta \right)}{\partial \eta}\right)+1\right)\right)\text{d}\eta}.$

Generally,

$\begin{array}{l}{y}_{n+1}\left(x,h\right)=1+{\displaystyle {\int}_{0}^{x}\left(x-\eta \right)(\frac{{\partial}^{2}{y}_{n}\left(\eta ,h\right)}{\partial {\eta}^{2}}+h(\frac{{\partial}^{2}{y}_{n}\left(\eta ,h\right)}{\partial {\eta}^{2}}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{2}{\eta}\left(\frac{\partial {y}_{n}\left(\eta ,h\right)}{\partial \eta}\right)+1))\text{d}\eta ,\text{\hspace{0.17em}}n\ge 1.\end{array}$ (29)

First few terms of solution are

$\begin{array}{l}{y}_{0}\left(x\right)=y\left(0\right)=1,\\ {y}_{1}\left(x,h\right)=1+h\frac{{x}^{2}}{2},\\ {y}_{2}\left(x,h\right)=1+\left(h+\frac{3{h}^{2}}{2}\right){x}^{2},\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\vdots \end{array}$

It can be seen in Figure 1 that the admissible range of h is $-0.6\le h\le -0.1$. The h curve of $y\left(x\right)$ for 20th order approximation when $x=1$, absolute error for15th order approximation by standard VPM and absolute error when $h=-0.3$ are shown in Figure 1.

(a)(b)(c)

Figure 1. For (k = 0) in Example 1 (a) the h curve for the 20th order approximation; (b) Absolute error by the standard VPM; (c) Absolute error when $h=-0.3$.

Case (2) for (k = 1):

The closed form solution is $y\left(x\right)=\frac{\mathrm{sin}x}{x}$. According to standard VPM, the absolute error of ${y}_{30}\left(x\right)$ is shown in Figure 2(a) which also shows that the result is not valid for large value of x. Now, using the iterative scheme (22) we have the first few terms of solution:

(a)(b)(c)(d)

Figure 2. For (k = 1) in Example 1 (a) Absolute error by the standard VPM; (b) Absolute error when
$h=-0.35$ ; (c) the h curve for the 30^{th} order approximation; (d) Series and exact solutions.

$\begin{array}{l}{y}_{0}\left(x\right)=1,\\ {y}_{1}\left(x,h\right)=1+h\frac{{x}^{2}}{2},\\ {y}_{2}\left(x,h\right)=1+\left(h+\frac{3{h}^{2}}{2}\right){x}^{2}+\frac{{h}^{2}{x}^{4}}{24},\end{array}$

Other terms of the solution can also be obtained in a similar way.

From Figure 2(c) the admissible range is $-0.59\le h\le -0.05$ and after plotting a number of h curves we get the best value of h which is equal to −0.35, therefore the approximate solution is given by

$y\left(x\right)=\underset{n\to \infty}{\mathrm{lim}}{y}_{n}\left(x,h\right)\approx \frac{\mathrm{sin}x}{x}.$ (30)

An approximate solution (not exact) was obtained, using VIM, by [12] and is valid only for $0\le x<1$. The solution is convergent in the interval $\subseteq \left[0,14\right]$ as shown in Figure 2(d). The accuracy of absolute error is remarkably improved for 30th order approximation, as illustrated in Figure 2(b).

Case (3) for (k = 2):

In this case we analyze the dependence of the convergence regions on the value of $\mathfrak{a}$, The h curve for the ninth-order approximation at $x=1$ is shown in Figure 3, evaluating $h=-0.4$. Figure 4 shows the solution at $\mathfrak{a}=1,4,5,10$ respectively, it is obvious that when the value of $\mathfrak{a}$ increases monotonically, the region of convergence monotonically decreases.

Case (4) for (k = 3, k = 4):

For these cases, a comparison is made through Table 1 and Table 2, respectively between the solution of Equation (26) and that obtained by [15] . It is shown that the obtained solution is accurate. h-curves at $x=1$ are shown in Figure 5(a) and Figure 5(b), respectively. The obtained results by the used method are compared with standard VPM through Table 3 and Table 4 respectively. It is clear from tables that the standard VPM is not efficient for large values of x and insertion of h to the method is very effective step in the solution.

Case (5) for (k = 5):

The closed form solution is $y\left(x\right)={\left(1+\frac{{x}^{2}}{3}\right)}^{-1/2}$. Adopting the similar

Figure 3. The h curve for (k = 2) in Example 1 for the 9th order approximation.

(a)(b)(c)(d)

Figure 4. For (k = 2) in Example 1 (a) solution with $\mathfrak{a}=1$ ; (b) solution with $\mathfrak{a}=4$ ; (c) solution with $\mathfrak{a}=5$ ; (d) solution with $\mathfrak{a}=10$.

(a)(b)

Figure 5. F or Example 1 (a) the h curve for the seventh-order approximation for (k = 3); (b) the h curve for the fifth-order approximation for (k = 4).

Table 1. Exact value, ${y}_{7}\left(x\right)$, the absolute error ${E}_{7}$ and the error obtained by [15] for $k=3$.

Table 2. Exact value, ${y}_{5}\left(x\right)$, the absolute error ${E}_{5}$ and the error obtained by [15] for $k=4$.

Table 3. Solutions obtained by the standard VPM and VPM with an auxiliary parameter for $k=3$.

Table 4. Solutions obtained by the standard VPM and VPM with an auxiliary parameter for $k=4$.

procedure as in the previous examples, Figure 6(a) shows the absolute error of ${y}_{5}\left(x\right)$ which confirm that the obtained result by standard VPM is not valid for large value of x. The first few terms of solution are

$\begin{array}{l}{y}_{0}\left(x\right)=1,\\ {y}_{1}\left(x,h\right)=1+h\frac{{x}^{2}}{2},\\ {y}_{2}\left(x,h\right)=1+\left(h+\frac{3{h}^{2}}{2}\right){x}^{2}+\frac{5{h}^{2}{x}^{4}}{24}+\frac{{h}^{3}{x}^{6}}{12}+\frac{5{h}^{4}{x}^{8}}{224}+\frac{{h}^{5}{x}^{10}}{288}+\frac{{h}^{6}{x}^{12}}{4224},\end{array}$

Other terms of the solution can also be obtained in a similar way, h can be chosen in the range
$-0.28\le h\le -0.42$, from Figure 6(b), the h curve for
$x=1$, the value of h equal to

Example 2

For the Isothermal gas spheres equation:

${y}^{\u2033}+\frac{2}{x}{y}^{\prime}+{\text{e}}^{y}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}x\ge 0,$ (31)

with the initial conditions ( $y\left(0\right)=0,{y}^{\prime}\left(0\right)=0$ ), we can use Taylor series expansion of ${\text{e}}^{y}$

${\text{e}}^{y}\cong 1+y+\frac{{y}^{2}}{2}+\frac{{y}^{3}}{6}+\frac{{y}^{4}}{24}.$

(a)(b)(c)(d)

Figure 6. For (k = 5) in Example 1 (a) Absolute error by the standard VPM; (b) The h curve; (c) Absolute error when $h=-0.4$ ; (d) Series and exact solutions.

Table 5 shows the results obtained with standard VPM and VPM with h. Now applying the iterative scheme (22) we have:

${y}_{0}\left(x\right)=y\left(0\right)=0,$ (32)

${y}_{1}\left(x,h\right)=\frac{h{x}^{2}}{2},$ (33)

${y}_{2}\left(x,h\right)=\left(h+\frac{3{h}^{2}}{2}\right){x}^{2}+\frac{{h}^{2}{x}^{4}}{24}+\frac{{h}^{3}{x}^{6}}{240}+\frac{{h}^{4}{x}^{8}}{2688}+\frac{{h}^{5}{x}^{10}}{34560}$ (34)

Other terms of the solution can also be obtained in a similar way to get:

$y\left(x\right)\cong -\frac{1}{6}{x}^{2}+\frac{1}{120}{x}^{4}-\frac{1}{1890}{x}^{6}+\cdots $

(with h equal to −0.395) which is the same as the solution obtained by [4] using ADM, and [13] using MHAM. Figure 7 shows the h curve for $x=1$. Table 5 shows the comparison of $y\left(x\right)$ obtained here and those obtained by [4] and [15] . The comparison between the solutions obtained by the standard VPM and VPM with an auxiliary parameter for isothermal gas sphere equation is tabulated in Table 6.

Figure 7. The h curve for isothermal gas sphere equation in Example 2.

Table 5. Comparison of $y\left(x\right)$ between present method and series solution given by [4] and [15] for isothermal gas sphere equation in Example 2..

Table 6. Comparison between the solutions obtained by the standard VPM and VPM with an auxiliary parameter for isothermal gas sphere equation in Example 2.

Example 3

Consider the following problem [13] :

${y}^{\u2033}\left(x\right)+\frac{2}{x}{y}^{\prime}\left(x\right)-2\left(2{x}^{2}+3\right)y\left(x\right)=0,$ (35)

With the initial conditions

$y\left(0\right)=1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}^{\prime}\left(0\right)=0.$ (36)

having $y\left(x\right)={\text{e}}^{{x}^{2}}$ as exact solution. Now, using the iterative scheme (22) we have

$\begin{array}{l}{y}_{0}\left(x\right)=y\left(0\right)=1,\\ {y}_{1}\left(x,h\right)=1-3h{x}^{2}-h\frac{{x}^{4}}{3},\\ {y}_{2}\left(x,h\right)=1-6h{x}^{2}-9h{x}^{2}-2h\frac{{x}^{4}}{3}+17{h}^{2}\frac{{x}^{4}}{18}+7{h}^{2}\frac{{x}^{6}}{5}+{h}^{2}\frac{{x}^{8}}{42}\end{array}$

Hence, $y\left(x\right)\cong {\text{e}}^{{x}^{2}}$ which is the exact solution of Equation (45) with $h=-0.4$ as the range of h is $-0.2\le h\le -0.6$ as shown in Figure 8(a). The exact solution and series solution for ( $h=-0.4$ ) absolute error by standard VPM and by VPM with h are shown in Figures 8(b)-(e), respectively. Figure 8(c) shows that VPM with an auxiliary parameter gives better approximation than MHAM [13] . The obtained solution, in comparison to that obtained by [15] , is very accurate. Table 7 shows the comparison of $y\left(x\right)$ obtained here and those obtained by [15] .

In the next examples, the nonhomogeneous Lane-Emden equations are considered.

Example 4

${y}^{\u2033}\left(x\right)+\frac{2}{x}{y}^{\prime}\left(x\right)+y\left(x\right)=6+12x+{x}^{2}+{x}^{3}$ (37)

With the initial conditions

$y\left(0\right)=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}^{\prime}\left(0\right)=0,$ (38)

Which has the following exact solution: $y\left(x\right)={x}^{2}+{x}^{3}$. Now using the iterative scheme (22) we have:

(a)(b)(c)(d)(e)

Figure 8. For Example 3 (a) The h curve; (b) Exact solution; (c) Series solution with $h=-0.4$ ; (d) Absolute error by standard VPM; (e) Absolute error by VPM with $h=-0.4$.

Table 7. Comparison between the solutions obtained by the present method and series solution given by Parand et al. [15] for Example 3.

$\begin{array}{l}{y}_{0}\left(x\right)=y\left(0\right)=0,\\ {y}_{1}\left(x,h\right)=-3h{x}^{2}-2h{x}^{3}-h\frac{{x}^{4}}{12}-h\frac{{x}^{5}}{20},\\ {y}_{2}\left(x,h\right)=-6h{x}^{2}-9{h}^{2}{x}^{2}-4h{x}^{3}-4{h}^{2}{x}^{3}-h\frac{{x}^{4}}{6}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}-7{h}^{2}\frac{{x}^{4}}{18}-h\frac{{x}^{5}}{10}-7{h}^{2}\frac{{x}^{5}}{40}-{h}^{2}\frac{{x}^{6}}{360}-{h}^{2}\frac{{x}^{7}}{840},\end{array}$

After number of iterations we get an approximate solution for h in the range of $-0.2\le h\le -0.5$ at which $y\left(x\right)\cong {x}^{2}+{x}^{3}$ by choosing an auxiliary parameter $h=-0.4$. Figure 9(a) shows the h curve at $x=1$. Exact solution and series solution are shown in Figure 9(b) and Figure 9(c) respectively. Like all previous examples, the present method is more accurate than standard VPM especially for large value of x. The absolute error by standard VPM is shown in Figure 9(d) and the absolute error by VPM with $h=-0.4$ is shown in Figure 9(e) which shows that the solution obtained here is more accurate in comparison to that obtained with MHAM [13] and also to that obtained by [15] . Table 8 shows the comparison of $y\left(x\right)$ obtained after 20 approximations and those obtained by [15] .

Example 5

Consider the equation:

${y}^{\u2033}\left(x\right)+\frac{8}{x}{y}^{\prime}\left(x\right)+xy\left(x\right)={x}^{5}-{x}^{4}+44{x}^{2}-30x,$ (39)

with the initial conditions

$y\left(0\right)=0,\text{\hspace{0.17em}}{y}^{\prime}\left(0\right)=0.$ (40)

Having $y\left(x\right)={x}^{4}-{x}^{3}$ as exact solution, using the iterative scheme (18) we have:

(a)(b)(c)(d)(e)

Figure 9. For Example 4 (a) the h curve; (b) Exact solution; (c) series solution when $h=-0.4$ ; (d) absolute error by standard VPM; (e) absolute error by VPM when $h=-0.4$.

$\begin{array}{l}{y}_{0}\left(x\right)=y\left(0\right)=0,\\ {y}_{1}\left(x,h\right)=5h{x}^{3}-11h\frac{{x}^{4}}{3}+h\frac{{x}^{6}}{30}-h\frac{{x}^{7}}{42},\\ {y}_{2}\left(x,h\right)=10h{x}^{3}+25{h}^{2}{x}^{3}-22h\frac{{x}^{4}}{3}-121{h}^{2}\frac{{x}^{4}}{9}+h\frac{{x}^{6}}{15}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+19{h}^{2}\frac{{x}^{6}}{75}-h\frac{{x}^{7}}{21}-{h}^{2}\frac{{x}^{7}}{7}+\frac{{h}^{2}{x}^{9}}{2160}-\frac{{h}^{2}{x}^{10}}{3780},\end{array}$

The exact solution ( ${x}^{4}-{x}^{3}$ ) can be obtained when n approaches to infinity by choosing a suitable value for h in the range $-0.3\le h\le -0.12$ which can be obtained from the h curve at $x=1$ shown in Figure 10(a), exact and series solutions with $h=-0.2$ are shown in Figure 10(b).

(a)(b)

Figure 10. For Example 5 (a) The h curve; (b) Exact and series solution.

Table 8. Comparison of $y\left(x\right)$ given by present method at ( $h=-0.4$ ) and series solution given by [15] for Example 4.

4. Conclusion

In this article, Variation of Parameter Method is used coupled with an auxiliary parameter for the Lane-Emden equation. Moreover, the iterative scheme (22) reduces to standard VPM for $h=-1$. The numerical results obtained in this research are very effective, completely reliable and powerful in obtaining analytical solutions of many problems. The approach used in the solution provides a high degree of accuracy only with a few iterations and reduces the size of calculations without using any restrictive assumptions.

References

[1] Davis, H.T. (1962) Introduction to Nonlinear Differential and Integral Equation, Dover Publications, New York.

[2] Chandrasekhar, S. (1967) Introduction to the Study of Stellar Structure, Dover Publications, New York.

[3] Shawagfeh, N.T. (1993) Nonperturbative Approximate Solution for Lane-Emden Equation. Journal of Mathematical Physics, 34, 4364-4369.

https://doi.org/10.1063/1.530005

[4] Wazwaz, A.M. (2001) A New Algorithm for Solving Differential Equations of Lane-Emden Type. Applied Mathematics and Computation, 118, 287-310.

https://doi.org/10.1016/S0096-3003(99)00223-4

[5] Semary, M.S. and Hassan, H.N. (2018) The Homotopy Analysis Method for q-Difference Equations. Ain Shams Engineering Journal, 9, 415-421.

https://doi.org/10.1016/j.asej.2016.02.005

[6] Semary, M.S. and Hassan, H.N. (2012) Series Solution for Continuous Population Models for Single and Interacting Species by the Homotopy Analysis Method. Communications in Numerical Analysis, 2012, Article ID: cna-00106.

https://doi.org/10.5899/2012/cna-00106

[7] Semary, M.S. and Hassan, H.N. (2016) An Effective Approach for Solving MHD Viscous Flow Due to a Shrinking Sheet. Applied Mathematics and Information Sciences, 10, 1425-1432.

https://doi.org/10.18576/amis/100421

[8] Semary, M.S., Hassan, H.N. and Radwan, A.G. (2017) Single and Dual Solutions of Fractional Order Differential Equations Based on Controlled Picard’s Method with Simpson Rule. Journal of the Association of Arab Universities for Basic and Applied Sciences, 24, 247-253.

https://doi.org/10.1016/j.jaubas.2017.06.001

[9] Hassan, H.N. (2016) An Accurate Numerical Solution for the Modified Equal Width Wave Equation Using the Fourier Pseudo-Spectral Method. Journal of Applied Mathematics and Physics, 4, 1054-1067.

https://doi.org/10.4236/jamp.2016.46110

[10] Hassan, H.N. (2017) An Efficient Numerical Method for the Modified Regularized Long Wave Equation Using Fourier Spectral Method. Journal of the Association of Arab Universities for Basic and Applied Sciences, 24, 198-205.

https://doi.org/10.1016/j.jaubas.2016.10.002

[11] Hassan, H.N. and Saleh, H.K. (2010) The Solution of the Regularized Long Wave Equation Using the Fourier Leap-Frog Method. Zeitschrift für Naturforschung A, 65, 268-276.

https://doi.org/10.1515/zna-2010-0402

[12] Dehghan, M. and Shakeri, F. (2008) Approximate Solution of a Differential Equation Arising in Astrophysics Using Variational Iteration Method. New Astronomy, 13, 53-59.

https://doi.org/10.1016/j.newast.2007.06.012

[13] Singh, O.P., Pandey, R.K. and Singh, V.K. (2009) An Analytic Algorithm of Lane-Emden Type Equations Arising in Astrophysics Using Modified Homotopy Analysis Method. Computer Physics Communications, 180, 1116-1124.

https://doi.org/10.1016/j.cpc.2009.01.012

[14] Baranwal, V.K., Pandey, R.K., Tripathi, M.P. and Singh, O.P. (2012) An Analytic Algorithm of Lane-Emden-Type Equations Arising in Astrophysics—A Hybrid Approach. Journal of Theoretical and Applied Physics, 6, 2251-7235.

[15] Parand, K., Dehghan, M., Rezaei, A.R. and Ghaderi, S.M. (2010) An Approximation Algorithm for the Solution of the Nonlinear Lane-Emden Type Equations Arising in Astrophysics Using Hermite Functions Collocation Method. Computer Physics Communications, 181, 1096-1108.

https://doi.org/10.1016/j.cpc.2010.02.018

[16] Ma, W.-X. and You, Y.-C. (2005) Solving the Korteweg-de Vries Equation by Its Bilinear Form: Wronskian Solutions. Transactions of the American Mathematical Society, 357, 1753-1778.

https://doi.org/10.1090/S0002-9947-04-03726-2

[17] Ma, W.-X. and You, Y.-C. (2004) Rational Solutions of the Toda Lattice Equation in Casoratian Form. Chaos, Solitons and Fractals, 22, 395-406.

https://doi.org/10.1016/j.chaos.2004.02.011

[18] Ma, W.-X., Wu, H.-Y. and He, J.-S. (2007) Partial Differential Equations Possessing Frobenius Integrable Decompositions. Physics Letters A, 364, 29-32.

https://doi.org/10.1016/j.physleta.2006.11.048

[19] Noor, M.A., Mohyud-Din, S.T. and Waheed, A. (2008) Variation of Parameters Method for Solving Fifth-Order Boundary Value Problems. Applied Mathematics & Information Sciences, 2, 135-141.

[20] Mohyud-Din, S.T., Noor, M.A. and Waheed, A. (2009) Variation of Parameter Method for Solving Sixth-Order Boundary Value Problems. Communications of the Korean Mathematical Society, 24, 605-615.

[21] Mohyud-Din, S.T., Noor, M.A. and Waheed, A. (2010) Variation of Parameter Method for Initial and Boundary Value Problems. World Applied Sciences Journal, 11, 622-639.

[22] Mohyud-Din, S.T., Noor, M.A. and Noor, K.I. (2009) Modified Variation of Parameters Method for Second-Order Integro-Differential Equations and Coupled Systems. World Applied Sciences Journal, 6, 1139-1146.

[23] Ghaneai, H. and Hosseini, M.M. (2015) Variational Iteration Method with an Auxiliary Parameter for Solving Wave-Like and Heat-Like Equations in Large Domains. Computers & Mathematics with Applications, 69, 363-373.

https://doi.org/10.1016/j.camwa.2014.11.007

[24] Semary, M.S. and Hassan, H.N. (2015) A New Approach for a Class of Nonlinear Boundary Value Problems with Multiple Solutions. Journal of the Association of Arab Universities for Basic and Applied Sciences, 17, 27-35.

https://doi.org/10.1016/j.jaubas.2013.11.002

[25] Sikandar, W., Khan, U., Ahmed, N. and Mohyud-Din, S.T. (2018) Variation of Parameters Method with an Auxiliary Parameter for Initial Value Problems. Ain Shams Engineering Journal, 9, 1959-1963.

https://doi.org/10.1016/j.asej.2016.09.014

[26] Sikandar, W., Khan, U., Ahmed, N. and Mohyud-Din, S.T. (2017) Optimal Solutions for Homogeneous and Non-Homogeneous Equations Arising in Physics. Results in Physics, 7, 216-224.

https://doi.org/10.1016/j.rinp.2016.12.018

[27] Ramos, J.I. (2008) On the Variational Iteration Method and Other Iterative Techniques for Nonlinear Differential Equations. Applied Mathematics and Computation, 199, 39-69.

https://doi.org/10.1016/j.amc.2007.09.024