Adomian Modification Methods for the Solution of Chebyshev’s Differential Equations

Mariam Al Mazmumy, Aishah Alsulami, Huda Bakodah, Nawal Alzaid^*

Department of Mathematics, College of Science, University of Jeddah, Jeddah, Saudi Arabia.

DOI: 10.4236/am.2023.148032                

Abstract

The current study examines the important class of Chebyshev’s differential equations via the application of the efficient Adomian Decomposition Method (ADM) and its modifications. We have proved the effectiveness of the employed methods by acquiring exact analytical solutions for the governing equations in most cases; while minimal noisy error terms have been observed in a particular method modification. Above all, the presented approaches have rightly affirmed the exactitude of the available literature. More to the point, the application of this methodology could be extended to examine various forms of high-order differential equations, as approximate exact solutions are rapidly attained with less computation stress.

Keywords

ADM, Modifications Methods, Chebyshev’s Differential Equations, IVPs, Series Solutions

1. Introduction

Differential equations are used to model problems arising in physical and social sciences, to mention a few. Recently, the development of computational procedures to tackle functional equations stumbled upon in different areas of science and technology has broadened the inquisitiveness of so many researchers. For instance, let’s consider the case of Chebyshev’s differential equation that reads

$\left(1-x^2\right){u^{\prime }}^{\prime }-x{u}^{\prime }+n^{2}u=0,-1<x<1,$ (1)

where n is a non-negative integer. This equation is undoubtedly a special type of second-order Ordinary Differential Equation (ODE) that arises in various science and engineering problems [1] . This equation is further typified with the possession of orthogonal functions that are known to play vital parts in numerical methods and approximation theories.

Further, it will also be good to remember that in recent times, research activities regarding second-order ODEs with initial data have drawn the nosiness of different researchers. Various methods have been proposed in the past and present literature towards devising promising unified techniques to tackle a verity of ODEs; one could easily find the Adomian Decomposition Method (ADM). [2] [3] , and its related modifications and extensions to attract so many minds in this regards, see the following references [4] - [10] and the references therewith to explore numerous scientific models in the presence of the method.

However, as the literature is kind of vacant with regards to the application of ADM on the solution of Chebyshev’s differential equation, the current study attempts to fill this huge gap by demonstrating the advantages and the efficiency of ADM and its modifications in attaining optimal approximate solutions to the Chebyshev’s differential equation; and by extension, to the general differential and integral equations. Additionally, we organize the paper in the following manner: Section 2 gives the outline of the classical Adomian method; while its modifications to be utilized in this study are presented in Section 3. Section 4 makes consideration to certain numerical applications, and lastly, we give certain concluding remarks in Section 5.

2. Standard ADM Procedure

The present section gives a generalized derivation procedure for tackling nonlinear Initial-Value Problems (IVPs) based on the ADM. To do so, let us consider the following differential equation

$G\left(u\left(x\right)\right)=g\left(x\right)\mathrm{,}$ (2)

with G representing a generalized ordinary (or partial) differential operator, and $g\left(x\right)$ as a source term. This operator being general, it can equally be expressed to involve both linear and nonlinear operators. Thus, we decompose the operator further and rewrite the above equation as follows

$Lu+Ru+Nu=g,$ (3)

where L is the highest linear operator that is invertible, with $R<L$ ; while N is specifically the nonlinear operator. More so, we rewrite the latter equation as follows

$Lu=g-Ru-Nu,$

such that applying the inverse linear operator $L^{-1}$ to both sides of the above equation yields

$u=\varphi \left(x\right)+L^{-1}g-L^{-1}Ru-L^{-1}Nu.$ (4)

where $\varphi \left(x\right)$ is the function emanating from the prescribed initial data.

Further, the iterative procedure by the name ADM decomposes the solution $u\left(x\right)$ using an infinite series of the following form

$u\left(x\right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{u}_n\left(x\right),$ (5)

while the nonlinear component $Nu$ is equally decomposed using the following infinite series

$N\left(u\right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{A}_n,$ (6)

where A_n’s are polynomials devised by Adomian, and recursively determined using the following scheme

$A_n=\frac{1}{n!}\frac{{\text{d}}^n}{\text{d}{\lambda }^n}{\left[N\left({\displaystyle \underset{j=0}{\overset{n}{\sum }}}{\lambda }^j{u}_j\right)\right]}_{\lambda =0},n=0,1,2,\cdots$ (7)

Therefore, upon substituting Equations (5) and (6) into Equation (4), one gets

where $\varphi \left(x\right)$ is the function emanating from the prescribed initial data.

Furthermore, the ADM procedure swiftly reveals the generalized recursive solution for the problem from the above equation as follows

$\{\begin{array}{l}{u}_0=\varphi \left(x\right)+L^{-1}g\left(x\right),\\ u_{k+1}=-L^{-1}R{u}_k-L^{-1}{A}_k,k\ge 0,\end{array}$ (8)

where A_k’s are the Adomian polynomials computed from Equation (7). Expressing few of these terms, we get

$A_0=N\left(u_0\right),$

$A_1=\frac{\text{d}N\left(u_0\right)}{\text{d}{u}_0}{u}_1,$

$A_2=\frac{\text{d}N\left(u_0\right)}{\text{d}{u}_0}{u}_2+\frac{1}{2}\frac{{\text{d}}^{2}N\left(u_0\right)}{\text{d}{u}_0^2}{u}_1^2,$

$A_3=\frac{\text{d}N\left(u_0\right)}{\text{d}{u}_0}{u}_3+\frac{{\text{d}}^{2}N\left(u_0\right)}{\text{d}{u}_0^2}{u}_1{u}_2+\frac{1}{3!}\frac{{\text{d}}^{3}N\left(u_0\right)}{\text{d}{u}_0^3}{u}_1^3,$

$⋮$

Remarkable, it is obvious that the Adomian polynomials A_n’s depend on the solution components u_n. For instance, A₀ relies merely on u₀; A₁ relies merely on u₀ and u₁; A₂ relies merely on u₀, u₁ and u₂, and so on.

Finally, a realistic solution is obtained by considering the following m-term approximations as

${\Psi }_n={\displaystyle \underset{j=0}{\overset{n-1}{\sum }}}{u}_j,$ (9)

where

$u\left(x\right)=\underset{n\to \infty }{\lim }{\Psi }_n\left(x\right)={\displaystyle \underset{j=0}{\overset{\infty }{\sum }}}{u}_j\left(x\right).$ (10)

3. Adomian Modification Methods for Chebyshev’s Equations

Here, we give some important modification algorithms based on the standard Adomian method for the solution of the class of Chebyshev’s differential equations under consideration.

3.1. Standard ADM

Consider the Chebyshev’s differential equation as follows

$\left(1-x^2\right)u^{″}-x{u}^{\prime }+n^{2}u=0,-1<x<1.$ (11)

Now, based on the standard ADM procedure, the Chebyshev’s equation given above is rewritten as follows

$u^{″}=x^2{u}^{″}+x{u}^{\prime }-n^{2}u,-1<x<1.$ (12)

So, we consider the right-hand side of the equation as a normal non-homogeneous term, where the differential operator L is defined by $\frac{{\text{d}}^2}{\text{d}{x}^2}$ . More

so, we consider the inverse operator $L^{-1}$ as a two-fold integral operator defined by

$L^{-1}\left(\mathrm{.}\right)={\displaystyle {\int }_0^x}{\displaystyle {\int }_0^x}\left(\mathrm{.}\right)\text{d}x\text{d}x\mathrm{.}$

Thus, upon applying the inverse operator $L^{-1}$ to both sides of Equation (12), one gets

$u=\varphi \left(x\right)+L^{-1}\left[x^2{u}^{″}+x{u}^{\prime }-n^{2}u\right]\mathrm{,}$

such that

$L\varphi \left(x\right)=0.$

Thus, based on the ADM procedure, the solution $u\left(x\right)$ is introduced through an infinite summation of components $u_n\left(x\right)$ earlier discussed in the methodology. Hence, the recursive solution of the equation is obtained as follows

$\{\begin{array}{l}{u}_0\left(x\right)=\varphi \left(x\right)=c_0+c_{1}x,\\ u_{k+1}\left(x\right)=L^{-1}\left[x^2{u^{″}}_k+x{u^{\prime }}_k-n^2{u}_k\right],k\ge 0,\end{array}$ (13)

where the overall solution $u\left(x\right)$ follows immediately by summing the above components as follows

$u\left(x\right)=\underset{n\to \infty }{\lim }{\Psi }_n\left(x\right)={\displaystyle \underset{j=0}{\overset{\infty }{\sum }}}{u}_j\left(x\right).$ (14)

3.2. Algorithm 1

This algorithm is aimed at tackling ODEs with non-constant coefficients, just like the Chebyshev’s equation which is a singular equation. The idea is more general, and could be extended to other higher-order ODEs that arise in mathematical physics [11] . To this end, let us consider the Chebyshev’s equation given in Equation (11). We then further divide the equation by $\left(1-x^2\right)$ , thereby expressing ${\left(1-x^2\right)}^{-1}$ as follows

$\frac{1}{1-x^2}={\displaystyle \underset{m=0}{\overset{\infty }{\sum }}}{x}^{2m},$

such that expressing the equation in an operator form as $Lu+Ru=0$ suitable for the decomposition method; where $L=\frac{{\text{d}}^2}{\text{d}{x}^2}$ , and

$R=-x{\displaystyle \underset{m=0}{\overset{\infty }{\sum }}}{x}^{2m}\frac{\text{d}}{\text{d}x}+n^2{\displaystyle \underset{m=0}{\overset{\infty }{\sum }}}{x}^{2m}.$ (15)

More so, applying the inverse operator $L^{-1}$ , which is an indefinite two-fold, to both sides of the governing equation, after being written in operator notation yields

$u=\varphi -L^{-1}Ru,\varphi =c_0+c_{1}x.$ (16)

So, on using $u={\displaystyle {\sum }_{m=0}^{\infty }}{u}_m$ , we rewrite the above equation as follows

${\displaystyle \underset{m=0}{\overset{\infty }{\sum }}}{u}_m=c_0+c_{1}x-L^{-1}R{\displaystyle \underset{m=0}{\overset{\infty }{\sum }}}{u}_m,$ (17)

where we can conveniently identify the recursive solution as follows

$\{\begin{array}{l}{u}_0=\varphi \left(x\right)=c_0+c_{1}x,\\ u_{k+1}\left(x\right)=-L^{-1}R{u}_k,k\ge 0,\end{array}$ (18)

where all the components are recursively determined. Thus, on making use of ${\Psi }_n={\displaystyle {\sum }_{m=0}^{n-1}}{u}_m$ , one gets n-term approximation suitable for numerical simulation.

Furthermore, one could equally decompose the term $\varphi$ resulting from the initial conditions as follows

such that the recursive relation above becomes

$\{\begin{array}{l}{u}_0\left(x\right)=c_{0,0}+x{c}_{1,0},\\ u_{k+1}\left(x\right)=c_{0,\left(k+1\right)}+x{c}_{1,\left(k+1\right)}-L^{-1}R{u}_k,k\ge 0.\end{array}$ (19)

This last development would no doubt ease the computational integration process; having further sub-decomposed the first component, as asserted in [11] with the initial data $c_0=u\left(0\right)$ and $c_1=u^{\prime }\left(0\right)$ . Additionally, one would equally apply this procedure to boundary-value situations such that a constant of integration is added to satisfy the conditions for ${\varphi }_n$ ; nonlinear ODEs could in the same fashion be solved iteratively using the same procedure.

3.3. Algorithm 2

Here, we consider yet another algorithm based on the modification by Dita and Grama [12] to ODEs. To present the method, let us make consideration to the following generalized second-order linear ODE

$L\left(x\mathrm{,}D\right)u\left(x\right)-R\left(x\mathrm{,}u\mathrm{,}Du\mathrm{,}{D}^{2}u\right)=\mathrm{0,}$ (20)

where $L\left(x\mathrm{,}D\right)$ is the principal part, with $D=\frac{\text{d}}{\text{d}x}$ , R denotes the reminder that

is not included in $L\left(x\mathrm{,}D\right)$ ; while the negative sign “−” in the equation is considered for handiness. Additionally, the principal part has the following format

$L\left(x\mathrm{,}D\right)=h\left(x\right)Dp\left(x\right)D\mathrm{,}$ (21)

where both $h\left(x\right)$ and $p\left(x\right)$ are considered to be smooth functions of x. Remember also that, the biggest target here is to suitably devise a way of making an optimal choice for R and L in such a way that the resultant pseudo-Volterra integral equation would be solved so easily.

Now, since tackling Equation (20) is complicated, we, therefore, resort to re-expressing the equation as a usual nonhomogeneous equation and thereafter make use of a modification of the famous variation of parameters method to obtain its corresponding pseudo-Volterra integral equation. Equally, we can constructively determine an inverse to Equation (21) and be utilized. Thus, without much delay, a proper inversion for Equation (20) is chosen in the following form

$L^{-1}\left(x,D\right)u\left(x\right)={\displaystyle {\int }_0^x}\frac{\text{d}t}{p\left(t\right)}{\displaystyle {\int }_0^t}\frac{\text{d}y}{h\left(y\right)}u\left(y\right),$ (22)

where one could easily note that

$\left(L{L}^{-1}\right)u\left(x\right)=I\cdot u\left(x\right)\mathrm{,}\text{but}\left(L^{-1}L\right)u\left(x\right)\ne I\cdot u\left(x\right)\mathrm{,}I\equiv \text{identityoperator}\mathrm{.}$

Apparently, it is inferred from the second relation that the inverse operator $L^{-1}$ is truly not an inverse operator. However, upon considering initial conditions, say $u\left(0\right)$ and $u^{\prime }\left(0\right)$ , it is so. In fact,

$\begin{array}{c}\left(L^{-1}L\right)u\left(x\right)={\displaystyle {\int }_0^x}\frac{\text{d}t}{p\left(t\right)}{\displaystyle {\int }_0^t}\frac{\text{d}y}{h\left(y\right)}h\left(y\right)\frac{\text{d}}{\text{d}y}p\left(y\right)\frac{\text{d}u\left(y\right)}{\text{d}y}\\ ={\displaystyle {\int }_0^x}\frac{\text{d}t}{p\left(t\right)}{\left[p\left(y\right)\frac{\text{d}u}{\text{d}y}\right]}_0^t={\displaystyle {\int }_0^x}\frac{\text{d}u\left(t\right)}{\text{d}t}-p\left(0\right)u^{\prime }\left(0\right){\displaystyle {\int }_0^x}\frac{\text{d}t}{p\left(t\right)}\\ =u\left(x\right)-u\left(0\right)-p\left(0\right)u^{\prime }\left(0\right){\displaystyle {\int }_0^x}\frac{\text{d}t}{p\left(t\right)}.\end{array}$

Thus, we obtain the following Volterra integral equation from Equation (22)

$u\left(x\right)=u\left(0\right)+p\left(0\right)u^{\prime }\left(0\right){\displaystyle {\int }_0^x}\frac{\text{d}t}{p\left(t\right)}+{\displaystyle {\int }_0^x}\frac{\text{d}t}{p\left(t\right)}{\displaystyle {\int }_0^t}\frac{\text{d}y}{h\left(y\right)}R\left(y,Du\right),$ (23)

where its solution is sought to be of the following infinite series form

$u\left(x\right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{u}_n\left(x\right),$

such that the recursive relation takes the following form

(24)

What’s more, the components $u_k\left(x\right)$ in the above scheme are determined recursively via the application of the Picard successive approximation, that is,

$u_k=L^{-1}\left(x,D\right)u_{k-1}.$

3.4. Algorithm 3

Let us consider a more general second-order linear ODE of the form [12]

$u^{″}+a\left(x\right)u^{\prime }+b\left(x\right)u=h\left(x\right)\mathrm{.}$ (25)

Thus, assuming a solution $\phi \left(x\right)\ne 0$ to the associated homogeneous equation, we further make use of a modification of variation of parameters method to acquire the solution Equation (25) as follows

$u\left(x\right)=C_1\phi \left(x\right)+C_2\phi \left(x\right){\displaystyle \int }\frac{\text{d}x}{E\left(x\right){\phi }^2\left(x\right)}+\phi \left(x\right){\displaystyle \int }\frac{\text{d}x}{E\left(x\right){\phi }^2\left(x\right)}{\displaystyle \int E}\left(x\right)\phi \left(x\right)h\left(x\right)\text{d}x\mathrm{,}$ (26)

where $C_j$ , for $j=1,2$ are constants, and $E\left(x\right)={\text{e}}^{{\displaystyle \int a\left(x\right)\text{d}x}}$ .

Moreover, the inverse operator $L^{-1}$ in the above equation is an indefinite integral, which has to be transformed to a definite one based on the solution we are aiming for. So, considering Equation (26) as the starting point for the derivation of this modification method, let us make use of Equation (25) to transform the governing linear homogeneous equation into a pseudo-Volterra integral equation. By this, we separate a suitable part in the left-hand side of Equation (25) in such a way that the resulting solution is obtainable. Of course, the separation process takes into concern the behaviour of the solution in the neighborhood of the point where the solution is sought, in order to increase the possibilities of obtaining the full series in explicit form. Therefore, the easiest separation is through shifting the term $b\left(x\right)u\left(x\right)$ on the right-hand side of the equation, which warrants the solution of the left-side alone to be $u\left(x\right)=1$ .

Besides, Equation (26) denotes the starting point for the acquisition of the generalized decomposition approach. In fact, once could recover the Adomian’s form by assuming $a\left(x\right)=b\left(x\right)=0$ in Equation (25). Consequently, if $\phi \left(x\right)=1$ and $E\left(x\right)=1$ , a regular solution in the neighborhood of $x=0$ from Equation (26) takes the following form

$u\left(x\right)=C_1+C_{2}x+{\displaystyle {\int }_0^x\text{d}x}{\displaystyle {\int }_0^{x}h}\left(x\right)\text{d}x\mathrm{,}$ (27)

of which by identification

$C_1=u\left(0\right),C_2=u^{\prime }\left(0\right).$ (28)

4. Applications

The current section demonstrates the application of the approaches given earlier to study different variants of Chebyshev’s differential equations. These approaches which are based on the ADM, are reliable as can rightly be seen in what follows.

Example 4.1. Let us consider an IVP featuring Chebyshev’s differential equation as follows

1) Case 1: $n=1$

$\begin{array}{l}\left(1-x^2\right)u^{″}-x{u}^{\prime }+u=0,-1<x<1,\\ u\left(0\right)=0,u^{\prime }\left(0\right)=1.\end{array}$ (29)

Standard ADM

Therefore, we re-express the governing differential equation using operator notation as follows

$Lu=x^2{u}^{″}+x{u}^{\prime }-u\mathrm{,}$ (30)

and further apply the inverse operator $L^{-1}$ to the both sides of the latter equation to obtain

$u=x+L^{-1}\left[x^2{u}^{″}+x{u}^{\prime }-u\right]\mathrm{.}$

Hence, the application of the standard ADM gives the general recursive relation for the model as follows

$\begin{array}{l}{u}_0\left(x\right)=x,\\ u_{k+1}\left(x\right)=L^{-1}\left[x^2{u^{″}}_k+x{u^{\prime }}_k-u_k\right]=0,k\ge 0,\end{array}$

which gives upon summing the above iterates the following exact solution

$u\left(x\right)=x.$ (31)

In fact, this is a well-known exact analytical solution for the Chebyshev’s differential equation when $n=1$ .

Algorithm 1

Also, based on this algorithm, we rewrite the governing model using operator notation as follows

$Lu=\frac{1}{1-x^2}\left[x{u}^{\prime }-u\right],$ (32)

such that after applying the inverse operator $L^{-1}$ to both sides of the above equation gives

$u=x+L^{-1}\frac{1}{1-x^2}\left[x{u}^{\prime }-u\right]\mathrm{.}$

Thus, making use of the Taylor’s expansion on $\frac{1}{1-x^2}$ to obtain

$\frac{1}{1-x^2}={\displaystyle \underset{m=0}{\overset{3}{\sum }}}{x}^{2m}=1+x^2+x^4+x^6,$

we then proceed to write the resulting recursive relation of the governing model as follows

$\begin{array}{l}{u}_0=x,\\ u_{k+1}\left(x\right)=L^{-1}\left(1+x^2+x^4+x^6\right)\left[x{u^{\prime }}_k-u_k\right]=0,k\ge 0,\end{array}$

which obviously sums to similar exact solution as in above as follows

$u\left(x\right)=x.$

2) Case 2: $n=2$

$\begin{array}{l}\left(1-x^2\right)u^{″}-x{u}^{\prime }+4u=0,-1<x<1,\\ u\left(0\right)=-1,u^{\prime }\left(0\right)=0.\end{array}$ (33)

Standard ADM

In an operator notation, the governing differential equation is re-expressed as follows

$Lu=\left[x^2{u}^{″}+x{u}^{\prime }-4u\right]\mathrm{,}$ (34)

such that after applying the inverse operator $L^{-1}$ to both sides of the above equation yields

$u=-1+L^{-1}\left[x^2{u}^{″}+x{u}^{\prime }-4u\right].$

In the same fashion as in the preceding problem, the following recursive relation is thus obtained

$\begin{array}{l}{u}_0\left(x\right)=-1,\\ u_1\left(x\right)=L^{-1}\left[x^2{u^{″}}_0+x{u^{\prime }}_0-4u_0\right]=2x^2,\\ u_2\left(x\right)=L^{-1}\left[x^2{u^{″}}_1+x{u^{\prime }}_1-4u_1\right]=0,\\ u_{k+1}\left(x\right)=L^{-1}\left[x^2{u^{″}}_k+x{u^{\prime }}_k-4u_k\right]=0,k\ge 2,\end{array}$

which leads to the following exact solution

$u\left(x\right)=-1+2x^2.$ (35)

Notably, the exact solution obtained above is a well-known exact analytical solution for the Chebyshev’s differential equation when. In the same fashion, various exact solutions of the Chebyshev’s differential equations could be obtained upon choosing different values of n.

Algorithm 1

This approach requires the model to be written in the following operator notational form

$Lu=\frac{1}{1-x^2}\left[x{u}^{\prime }-4u\right]\mathrm{.}$ (36)

More, with the application of the inverse operator $L^{-1}$ on both sides of the above equation, one gets

$u=-1+L^{-1}\frac{1}{1-x^2}\left[x{u}^{\prime }-4u\right]\mathrm{,}$

which subsequently leads to the following final recursive relation

$\begin{array}{l}{u}_0=-1,\\ u_1\left(x\right)=L^{-1}\left(1+x^2+x^4+x^6\right)\left[x{u^{\prime }}_0-4u_0\right]=2x^2+\frac{1}{3}{x}^4+\frac{2}{15}{x}^6+\frac{1}{14}{x}^8,\\ u_2\left(x\right)=L^{-1}\left(1+x^2+x^4+x^6\right)\left[x{u^{\prime }}_1-4u_1\right]=-\frac{1}{3}{x}^4-\frac{2}{15}{x}^6-\frac{1}{15}{x}^8-\frac{181}{4725}{x}^{10}+\cdots \\ u_3\left(x\right)=L^{-1}\left(1+x^2+x^4+x^6\right)\left[x{u^{\prime }}_2-4u_2\right]=-\frac{1}{210}{x}^8-\frac{4}{675}{x}^{10}+\cdots \\ u_4\left(x\right)=L^{-1}\left(1+x^2+x^4+x^6\right)\left[x{u^{\prime }}_3-4u_3\right]=-\frac{1}{4725}{x}^{10}+\cdots \\ ⋮\end{array}$

and leading to the following series solution

$u\left(x\right)={\displaystyle \underset{k=0}{\overset{6}{\sum }}}{u}_k\left(x\right)=-1+2x^2-\frac{2}{45}{x}^{10}+\cdots ,$ (37)

where the exact solution of Chebyshev’s equation for $n=2$ is $T_2\left(x\right)=-1+2x^2$ . Thus, we report the absolute error difference between the exact and the obtained approximate solution of $u\left(x\right)$ in Table 1.

Table 1. Absolute errors of Chebyshev’s equation for $n=2$ .

Remark: In the same way, we compute the cases when $n=3,n=4$ , and $n=5$ and reported the results in Table 2 using the following initial conditions:

$\begin{array}{l}n=3:u\left(0\right)=0,u^{\prime }\left(0\right)=-3;\\ n=4:u\left(0\right)=1,u^{\prime }\left(0\right)=0;\\ n=5:u\left(0\right)=0,u^{\prime }\left(0\right)=5.\end{array}$

Remarkably, the application of the present algorithm on the solution of Chebyshev’s equation reveals that the error increases as the value of n increases.

Table 2. Absolute errors of Chebyshev’s equation for $n=3$ , $n=4$ , and $n=5$ .

Example 4.2. Let us consider the following form of Chebyshev’s differential equation

$\sqrt{1-x^2}\frac{\text{d}}{\text{d}x}\left(\sqrt{1-x^2}{u}^{\prime }\right)=-n^{2}u\mathrm{.}$ (38)

Algorithm 2 will be considered for the solution of this form of Chebyshev’s equation, with

$h\left(x\right)=p\left(x\right)=\sqrt{1-x^2}.$

1) Case 1: $n=1$

Consider the following IVP

$\begin{array}{l}\sqrt{1-x^2}\frac{\text{d}}{\text{d}x}\left(\sqrt{1-x^2}{u}^{\prime }\right)=-u,\\ u\left(0\right)=0,u^{\prime }\left(0\right)=1.\end{array}$ (39)

Decomposing the solution $u\left(x\right)$ as usual, and further make use of the form of Equation (23) by setting the Taylor’s expansion of $\frac{1}{\sqrt{1-x^2}}$ as

$\frac{1}{\sqrt{1-x^2}}=1+\frac{1}{2}{x}^2+\frac{3}{8}{x}^4+\frac{5}{16}{x}^6\mathrm{.}$

We, therefore, get via the algorithm the following

$u_0=u\left(0\right)+p\left(0\right)u^{\prime }\left(0\right){\displaystyle {\int }_0^x}\frac{\text{d}t}{p\left(t\right)}={\displaystyle {\int }_0^x}\left(1+\frac{1}{2}{t}^2+\frac{3}{8}{t}^4+\frac{5}{16}{t}^6\right)\text{d}t$

$\begin{array}{c}{u}_{k+1}=L^{-1}\left[-u_k\left(x\right)\right],k\ge 0\\ ={\displaystyle {\int }_0^x}\left(1+\frac{1}{2}{x}^2+\frac{3}{8}{x}^4+\frac{5}{16}{x}^6\right){\displaystyle {\int }_0^x}\left(1+\frac{1}{2}{x}^2+\frac{3}{8}{x}^4+\frac{5}{16}{x}^6\right)\left[-u_k\left(x\right)\right]\text{d}x\text{d}x.\end{array}$

That is,

$\begin{array}{l}{u}_0=x+\frac{1}{6}{x}^3+\frac{3}{40}{x}^5+\frac{5}{112}{x}^7\\ u_1=L^{-1}\left[-u_0\left(x\right)\right]=-\frac{1}{6}{x}^3-\frac{1}{12}{x}^5-\frac{37}{720}{x}^7-\frac{3229}{90720}{x}^9-\frac{253}{22400}{x}^{11}+\cdots \\ u_2=L^{-1}\left[-u_1\left(x\right)\right]=\frac{1}{120}{x}^5+\frac{1}{144}{x}^7+\frac{47}{8640}{x}^9+\frac{1571}{362880}{x}^{11}+\cdots \\ u_3=L^{-1}\left[-u_2\left(x\right)\right]=-\frac{1}{5040}{x}^7-\frac{1}{4320}{x}^9-\frac{19}{86400}{x}^{11}+\cdots \\ u_4=L^{-1}\left[-u_3\left(x\right)\right]=\frac{1}{362880}{x}^9+\frac{1}{241920}{x}^{11}+\cdots \\ ⋮\end{array}$

and further leads to the following series solution

$u\left(x\right)={\displaystyle \underset{k=0}{\overset{6}{\sum }}}{u}_k\left(x\right)=x-\frac{35}{1152}{x}^9-\frac{91}{12672}{x}^{11}-\frac{1561}{658944}{x}^{13}+\cdots$ (40)

where the exact solution of Chebyshev’s equation for $n=1$ is $T_1\left(x\right)=x$ as found in the above approach; the remaining terms could surely be noisy terms because of the method. Thus, we report the absolute error difference between the exact and the obtained approximate solution of $u\left(x\right)$ in Table 3.

Table 3. Absolute errors of Chebyshev’s equation for $n=1$ .

2) Case 2: $n=2$

Consider the following IVP

$\begin{array}{l}\sqrt{1-x^2}\frac{\text{d}}{\text{d}x}\left(\sqrt{1-x^2}{u}^{\prime }\right)=-4u,\\ u\left(0\right)=-1,u^{\prime }\left(0\right)=0.\end{array}$ (41)

In same manner, we make use of the Taylor’s expansion on $\frac{1}{\sqrt{1-x^2}}$ , and eventually via the algorithm arrive at the following recursive relation

$\begin{array}{c}{u}_{k+1}=L^{-1}\left[-4u_k\left(x\right)\right],k\ge 0\\ ={\displaystyle {\int }_0^x}\left(1+\frac{1}{2}{x}^2+\frac{3}{8}{x}^4+\frac{5}{16}{x}^6\right){\displaystyle {\int }_0^x}\left(1+\frac{1}{2}{x}^2+\frac{3}{8}{x}^4+\frac{5}{16}{x}^6\right)\left[-4u_k\left(x\right)\right]\text{d}x\text{d}x,\end{array}$

such that

$\begin{array}{l}{u}_0=-1\\ u_1=L^{-1}\left[-4u_0\left(x\right)\right]=2x^2+\frac{2}{3}{x}^4+\frac{16}{45}{x}^6+\frac{8}{35}{x}^8+\frac{689}{16800}{x}^{10}+\cdots \\ u_2=L^{-1}\left[-4u_1\left(x\right)\right]=-\frac{2}{3}{x}^4-\frac{4}{9}{x}^6-\frac{14}{45}{x}^8-\frac{656}{2835}{x}^{10}+\cdots \\ u_3=L^{-1}\left[-4u_2\left(x\right)\right]=\frac{4}{45}{x}^6+\frac{4}{45}{x}^8+\frac{52}{675}{x}^{10}+\cdots \end{array}$

$\begin{array}{l}{u}_4=L^{-1}\left[-4u_3\left(x\right)\right]=-\frac{2}{315}{x}^8-\frac{8}{945}{x}^{10}+\cdots \\ u_5=L^{-1}\left[-4u_4\left(x\right)\right]=\frac{4}{14175}{x}^{10}+\cdots \\ ⋮\end{array}$

which further leads to the following series solution

$u\left(x\right)={\displaystyle \underset{k=0}{\overset{6}{\sum }}}{u}_k\left(x\right)=-1+2x^2-\frac{35}{288}{x}^{10}+\cdots$ (42)

where the exact solution of Chebyshev’s equation for $n=2$ is $T_2\left(x\right)=-1+2x^2$ as found in the previous approach. Thus, we report the absolute error difference between the exact and the obtained approximate solution of $u\left(x\right)$ in Table 4.

Table 4. Absolute errors of Chebyshev’s equation for $n=2$ .

Remark: In the same way, we compute the cases when $n=3$ , $n=4$ , and $n=5$ as reported in Table 5 using the following initial conditions:

$\begin{array}{l}n=3:u\left(0\right)=0,u^{\prime }\left(0\right)=-3;\\ n=4:u\left(0\right)=1,u^{\prime }\left(0\right)=0;\\ n=5:u\left(0\right)=0,u^{\prime }\left(0\right)=5.\end{array}$

Table 5. Absolute errors of Chebyshev’s equation for $n=3$ , $n=4$ , and $n=5$ .

Example 4.3. Let us consider the Chebyshev’s differential equation expressed in the following form

$u^{″}+\frac{1}{x}{u}^{\prime }-\frac{n^2}{x^2}u=\frac{1}{x^2}{u}^{″}\mathrm{.}$ (43)

Algorithm 3 will be considered for the solution of this form of Chebyshev’s equation. What’s more, we consider the right-hand side of the equation as a normal nonhomogeneous term, and further seek a solution of the corresponding homogeneous part on the left-hand side. Let us try to establish a more general procedure based on the present algorithm to tackle the governing equation. Firstly, we assume ${\phi }_{+}=x^n$ and ${\phi }_{-}=x^{-n}$ to be two solutions of the equation.

Now, with the first solution ${\phi }_{+}=x^n$ , and let $C_1=1$ and $C_2=0$ . We then seek a regular solution at $x=0$ by considering the inverse operator $L^{-1}$ as a two-fold integration from 0 to x as in Equation (26). Thus, the procedure results in the following the Volterra integral equation

$u\left(x\right)=x^n+x^n{\displaystyle {\int }_0^x}{x}^{-2n-1}{\displaystyle {\int }_0^x}{x}^{n+1}\left[\frac{1}{x^2}{u}^{″}\right]\text{d}x\text{d}x.$ (44)

More so, solution of the above equation emerges by considering the zero-order approximation the nonhomogeneous term as follows

$u_0=x^n,$

and

$u_{k+1}=x^n{\displaystyle {\int }_0^x}{x}^{-2n-1}{\displaystyle {\int }_0^x}{x}^{n+1}\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x,k\ge 0.$

We, therefore, find the components as

$\begin{array}{l}{u}_1=x^n{\displaystyle {\int }_0^x}{x}^{-2n-1}{\displaystyle {\int }_0^x}{x}^{n+1}\left[\frac{1}{x^2}{u^{″}}_0\right]\text{d}x\text{d}x=-\frac{1}{4}n{x}^{n-2},\\ u_2=x^n{\displaystyle {\int }_0^x}{x}^{-2n-1}{\displaystyle {\int }_0^x}{x}^{n+1}\left[\frac{1}{x^2}{u^{″}}_1\right]\text{d}x\text{d}x=\frac{n\left(n-3\right)}{32}{x}^{n-4},\\ u_3=x^n{\displaystyle {\int }_0^x}{x}^{-2n-1}{\displaystyle {\int }_0^x}{x}^{n+1}\left[\frac{1}{x^2}{u^{″}}_2\right]\text{d}x\text{d}x=-\frac{n\left(n-4\right)\left(n-5\right)}{384}{x}^{n-6},\\ ⋮\end{array}$

In this way, we get

$u\left(x\right)=x^n-\frac{1}{4}{x}^{n-2}+\frac{n\left(n-3\right)}{32}{x}^{n-4}-\frac{n\left(n-4\right)\left(n-5\right)}{384}{x}^{n-6}+\cdots$ (45)

Also, from the second solution ${\phi }_{-}=x^{-n}$ , we let $C_1=0$ and $C_2=n\left(2^n\right)$ . Then, without lost of generality, we get the following integral equation

$u\left(x\right)=n\left(2^n\right)x^{-n}{\displaystyle {\int }_0^x}{x}^{2n-1}\text{d}x+x^{-n}{\displaystyle {\int }_0^x}{x}^{2n-1}{\displaystyle {\int }_0^x}{x}^{-n+1}\left[\frac{1}{x^2}{u}^{″}\right]\text{d}x\text{d}x.$ (46)

Therefore, iterative solution of the above equation is obtained recursively as in the preceding case as fallows

$u_0=n\left(2^n\right)x^{-n}{\displaystyle {\int }_0^x}{x}^{2n-1}\text{d}x=2^{n-1}{x}^n,$

and

$u_{k+1}=x^{-n}{\displaystyle {\int }_0^x}{x}^{2n-1}{\displaystyle {\int }_0^x}{x}^{-n+1}\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x,k\ge 0,$

where the components revealed explicitly are as follows

$\begin{array}{l}{u}_1=x^{-n}{\displaystyle {\int }_0^x}{x}^{2n-1}{\displaystyle {\int }_0^x}{x}^{-n+1}\left[\frac{1}{x^2}{u^{″}}_0\right]\text{d}x\text{d}x=-\frac{n\left(2^n\right)}{8}{x}^{n-2},\\ u_2=x^{-n}{\displaystyle {\int }_0^x}{x}^{2n-1}{\displaystyle {\int }_0^x}{x}^{-n+1}\left[\frac{1}{x^2}{u^{″}}_1\right]\text{d}x\text{d}x=\frac{n\left(n-3\right)\left(2^n\right)}{64}{x}^{n-4},\\ u_3=x^{-n}{\displaystyle {\int }_0^x}{x}^{2n-1}{\displaystyle {\int }_0^x}{x}^{-n+1}\left[\frac{1}{x^2}{u^{″}}_2\right]\text{d}x\text{d}x=-\frac{n\left(n-4\right)\left(n-5\right)\left(2^n\right)}{768}{x}^{n-6},\\ ⋮\end{array}$

which leads to the series solution

$u\left(x\right)=2^{n-1}{x}^n-\frac{n\left(2^n\right)}{8}{x}^{n-2}+\frac{n\left(n-3\right)\left(2^n\right)}{64}{x}^{n-4}-\frac{n\left(n-4\right)\left(n-5\right)\left(2^n\right)}{768}{x}^{n-6}+\cdots$

Thus, in what follows, we have made consideration to the cases of the governing model for some specific values of n.

1) Case 1: $n=1$

Let us consider the solutions to be of the following forms ${\phi }_{\pm }=x^{\pm 1}$ .

Now, for ${\phi }_{+}=x$ . We prescribe $C_1=1$ and $C_2=0$ . Therefore, without loss of generality, the following Volterra integral equation is acquired via the process

$u\left(x\right)=x+x{\displaystyle {\int }_0^x}{x}^{-3}{\displaystyle {\int }_0^x}{x}^2\left[\frac{1}{x^2}{u}^{″}\right]\text{d}x\text{d}x.$ (47)

which leads to the following iterative solution

$u_0=x,$

and

$u_{k+1}=x{\displaystyle {\int }_0^x}{x}^{-3}{\displaystyle {\int }_0^x}{x}^2\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x=0,k\ge 0,$

In fact, the series summation gives

$u\left(x\right)=x\mathrm{,}$ (48)

which is the exact analytical solution of the Chebyshev equation for $n=1$ .

In a similar way, we consider the solution ${\phi }_{-}=x^{-1}$ , with $C_1=0$ , and $C_2=2$ . Via the same process, the following integral equation is obtained

$u\left(x\right)=2x^{-1}{\displaystyle {\int }_0^x}x\text{d}x+x^{-1}{\displaystyle {\int }_0^x}x{\displaystyle {\int }_0^x}\left[\frac{1}{x^2}{u}^{″}\right]\text{d}x\text{d}x,$ (49)

that solves to the following scheme

$u_0=2x^{-1}{\displaystyle {\int }_0^x}x\text{d}x=x,$

and

$u_{k+1}=x^{-1}{\displaystyle {\int }_0^x}x{\displaystyle {\int }_0^x}\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x=0,k\ge 0,$

which equally gives the same exact solution

$u\left(x\right)=x.$

2) Case 2: $n=2$

Here, the solutions are considered as follows ${\phi }_{\pm }=x^{\pm 2}$ .

Now, for ${\phi }_{+}=x^2$ , we prescribe $C_1=1$ and $C_2=0$ . Therefore, without loss of generality, the following Volterra integral equation is acquired via the process

$u\left(x\right)=x^2+x^2{\displaystyle {\int }_0^x}{x}^{-5}{\displaystyle {\int }_0^x}{x}^3\left[\frac{1}{x^2}{u}^{″}\right]\text{d}x\text{d}x\mathrm{,}$ (50)

which leads to the following iterative solution

$u_0=x^2,$

and

$u_{k+1}=x^2{\displaystyle {\int }_0^x}{x}^{-5}{\displaystyle {\int }_0^x}{x}^3\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x,k\ge 0.$

More explicitly, we get

$\begin{array}{l}{u}_0=x^2,\\ u_1=x^2{\displaystyle {\int }_0^x}{x}^{-5}{\displaystyle {\int }_0^x}{x}^3\left[\frac{1}{x^2}{u^{″}}_0\right]\text{d}x\text{d}x=-\frac{1}{2},\\ u_{k+1}=x^2{\displaystyle {\int }_0^x}{x}^{-5}{\displaystyle {\int }_0^x}{x}^3\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x=0,k\ge 1.\end{array}$

In this way we get

$u\left(x\right)=x^2-\frac{1}{2},$ (51)

which is the exact solution of Chebyshev equation for $n=2$ .

In a similar way, we consider the solution ${\phi }_{-}=x^{-2}$ , with $C_1=0$ , and $C_2=8$ . Via the same process, the following integral equation is obtained

$u\left(x\right)=8x^{-2}{\displaystyle {\int }_0^x}{x}^3\text{d}x+x^{-2}{\displaystyle {\int }_0^x}{x}^3{\displaystyle {\int }_0^x}{x}^{-1}\left[\frac{1}{x^2}{u}^{″}\right]\text{d}x\text{d}x.$ (52)

which leads to the following iterative solution

$u_0=8x^{-2}{\displaystyle {\int }_0^x}{x}^3\text{d}x=2x^2,$

and

$u_{k+1}=x^{-2}{\displaystyle {\int }_0^x}{x}^3{\displaystyle {\int }_0^x}{x}^{-1}\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x,k\ge 0.$

More explicitly, we get

$\begin{array}{l}{u}_0=8x^{-2}{\displaystyle {\int }_0^x}{x}^3\text{d}x=2x^2,\\ u_1=x^{-2}{\displaystyle {\int }_0^x}{x}^3{\displaystyle {\int }_0^x}{x}^{-1}\left[\frac{1}{x^2}{u^{″}}_0\right]\text{d}x\text{d}x=-1,\\ u_{k+1}=x^{-2}{\displaystyle {\int }_0^x}{x}^3{\displaystyle {\int }_0^x}{x}^{-1}\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x=0,k\ge 1.\end{array}$

Thus, upon summing the above iterates, it yields

$u\left(x\right)=2x^2-1,$

which is the exact solution of Chebyshev equation for $n=2$ .

3) Case 3: $n=3$

We consider the following forms of solutions are ${\phi }_{\pm }=x^{\pm 3}$ .

Then, for ${\phi }_{+}=x^3$ , we prescribe $C_1=1$ and $C_2=0$ . Also, without loss of generality, the following Volterra integral equation is acquired via the process

$u\left(x\right)=x^3+x^3{\displaystyle {\int }_0^x}{x}^{-7}{\displaystyle {\int }_0^x}{x}^4\left[\frac{1}{x^2}{u}^{″}\right]\text{d}x\text{d}x.$ (53)

which leads to the following iterative solution

$u_0=x^3,$

and

$u_{k+1}=x^3{\displaystyle {\int }_0^x}{x}^{-7}{\displaystyle {\int }_0^x}{x}^4\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x,k\ge 0.$

That is,

$\begin{array}{l}{u}_0=x^3,\\ u_1=x^3{\displaystyle {\int }_0^x}{x}^{-7}{\displaystyle {\int }_0^x}{x}^4\left[\frac{1}{x^2}{u^{″}}_0\right]\text{d}x\text{d}x=-\frac{3}{4}x,\\ u_{k+1}=x^3{\displaystyle {\int }_0^x}{x}^{-7}{\displaystyle {\int }_0^x}{x}^4\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x=0,k\ge 1.\end{array}$

Also, upon summing the above iterates, we get

$u\left(x\right)=x^3-\frac{3}{4}x,$ (54)

which is the exact analytical solution of Chebyshev equation for $n=3$ .

In a similar pattern, we consider the other solution ${\phi }_{-}=x^{-3}$ , with $C_1=0$ , and $C_2=24$ . Also, without further delay, we get following integral equation

$u\left(x\right)=24x^{-3}{\displaystyle {\int }_0^x}{x}^5\text{d}x+x^{-3}{\displaystyle {\int }_0^x}{x}^5{\displaystyle {\int }_0^x}{x}^{-2}\left[\frac{1}{x^2}{u}^{″}\right]\text{d}x\text{d}x,$ (55)

that solves iteratively to

$u_0=24x^{-3}{\displaystyle {\int }_0^x}{x}^5\text{d}x=4x^3,$

and

$u_{k+1}=x^{-3}{\displaystyle {\int }_0^x}{x}^5{\displaystyle {\int }_0^x}{x}^{-2}\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x,k\ge 0.$

That is

$\begin{array}{l}{u}_0=24x^{-3}{\displaystyle {\int }_0^x}{x}^5\text{d}x=4x^3,\\ u_1=x^{-3}{\displaystyle {\int }_0^x}{x}^5{\displaystyle {\int }_0^x}{x}^{-2}\left[\frac{1}{x^2}{u^{″}}_0\right]\text{d}x\text{d}x=-3x,\\ u_{k+1}=x^{-3}{\displaystyle {\int }_0^x}{x}^5{\displaystyle {\int }_0^x}{x}^{-2}\left[\frac{1}{x^2}{u^{″}}_k\right]\text{d}x\text{d}x=0,k\ge 1.\end{array}$

Lastly, summation of the respective components yields

$u\left(x\right)=4x^3-3x\mathrm{,}$

which is the exact analytical solution for the Chebyshev equation when $n=3$ . In general one could obtain the exact solution of Chebyshev’s equation for any value of n following the same procedure.

5. Conclusion

As a closing remark, the Adomian decomposition method and its modifications were successfully employed in the current study to examine the important class of Chebyshev’s differential equations. These methods have been shown to reveal exact analytical solutions of the governing models in most cases; while minimal noisy error terms have been observed in a particular modification. Above all, these modification methods have rightly affirmed the exactness of the available results in the literature. Besides, the application of the presented sets of modification methods could be extended to examine various singular and nonsingular higher-order differential equations, as approximate exact solutions are rapidly attained with less computation stress.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Bell, W.W. (2004) Special Functions for Scientists and Engineers. Courier Corporation, London, UK.
[2]	Adomian, G. (1990) A Review of the Decomposition Method and Some Recent Results for Nonlinear Equations. Mathematical and Computer Modelling, 13, 17-43. https://doi.org/10.1016/0895-7177(90)90125-7
[3]	Adomian, G. (1994) Solving Frontier Problems of Physics: The Decomposition Method. Kluwer Academic Publishers, Dordrecht. https://doi.org/10.1007/978-94-015-8289-6
[4]	Wazwaz, A.-M. (1999) A Reliable Modification of Adomian Decomposition Method. Applied Mathematics and Computation, 102, 77-86. https://doi.org/10.1016/S0096-3003(98)10024-3
[5]	Wazwaz, A.-M. (2002) A New Method for Solving Singular Initial Value Problems in the Second-Order Ordinary Differential Equations. Applied Mathematics and Computation, 128, 45-57. https://doi.org/10.1016/S0096-3003(01)00021-2
[6]	Almazmumy, M., Hendi, F., Bakodah, H. and Alzumi, H. (2012) Recent Modifications of Adomian Decomposition Method for Initial Value Problem in Ordinary Differential Equations. American Journal of Computational Mathematics, 2, 228-234. https://doi.org/10.4236/ajcm.2012.23030
[7]	Bakodah, H.O., Banaja, M.A., Alrigi, B.A., Ebaid, A. and Rach, R. (2019) An Efficient Modification of the Decomposition Method with a Convergence Parameter for Solving Korteweg de Vries Equations. Journal of King Saud University-Science, 31, 1424-1430. https://doi.org/10.1016/j.jksus.2018.11.010
[8]	Alsulami, A.A., AL-Mazmumy, M., Bakodah, H.O. and Alzaid, N. (2022) A Method for the Solution of Coupled System of Emden-Fowler-Type Equations. Symmetry, 14, Article No. 843. https://doi.org/10.3390/sym14050843
[9]	Al-Mazmumy, M., Alsulami, A.A., Bakodah, H.O. and Alzaid, N. (2022) Modified Adomian Method through Efficient Inverse Integral Operators to Solve Nonlinear Initial-Value Problems for Ordinary Differential Equations. Axioms, 11, Article No. 698. https://doi.org/10.3390/axioms11120698
[10]	Al-Mazmumy, M., Alsulami, A., Bakodah, H. and Alzaid, N. (2022) Modified Adomian Method for the Generalized Inhomogeneous Lane-Emden-Type Equations. Nonlinear Analysis and Differential Equations, 10, 15-35. https://doi.org/10.12988/nade.2022.91142
[11]	Adomian, G. (1992) Differential Coefficients with Singular Coefficients. Applied Mathematics and Computation, 47, 179-184. https://doi.org/10.1016/0096-3003(92)90045-3
[12]	Dita, P. and Grama, N. (1997) On Adomian’s Decomposition Method for Solving Differential Equations. https://arxiv.org/abs/solv-int/9705008