Adomian iterations

Complete List of Publications using the Adomian Decomposition Method

Let me introduce some scientists (in alphabetic order) who contributed to and greatly improved the Adomian Decomposition Method to make it available to the mathematical and engineering community.

Dr. George Adomian (1922--1996) was an American mathematician, theoretical physicist, and electrical engineer of Armenian descent. He received his Ph.D. degree from UCLA. He first proposed and considerably developed the Adomian Decomposition Method (ADM) for solving nonlinear differential equations, both ordinary, and partial, deterministic and stochastic, also integral equations, algebraic and transcendental (functional), and matrix equations. He was a Distinguished Professor (academic rank), the David C. Barrow Professor of Mathematics (Chair), and the Director of the Center for Applied Mathematics at the University of Georgia, the founder and Chief Scientist of General Analytics Corporation, a winner of the 1989 Richard Bellman Prize for outstanding contributions to nonlinear stochastic analysis, and a 1988 National Academy of Sciences Scholar. He is the author of eight books and over three hundred journal papers. Dr. Adomian was the Radar Officer (naval rank of Lieutenant, the equivalent of a Captain in the Army by direct commission) aboard the USS Antietam (CV-36), and served in the Pacific Theater of Operations during the Second World War. He also attended the secret radar school at MIT that trained the first radar officers for the U.S. Navy. His obituary was prepared by Randolph Rach and published in
R. Rach, Dr. George Adomian – distinguished scientist and mathematician, Kybernetes, Vol. 25, No. 9, 1996, 45-50.

Dr. Yves Cherruault (1937--2010) was a French mathematician. He received his Ph.D. degree from the University of Paris. He was a professor of mathematics at the University Pierre et Marie Curie, Paris, and the Director of MÉDIMAT (Laboratory of Mathematics Applied to Biomedicine). Professor Cherruault is one of the founders of the field of Biomathematics and an author of seven books and over two hundred journal papers. Dr. Cherruault developed some of the convergence theorems of the ADM.

Dr. Jun-Sheng Duan is a Chinese applied mathematician and computer scientist born in Inner Mongolia, China, in 1965. He received his Ph.D. degree from Shandong University, China. He has made extensive contributions to solutions of differential equations in mathematical physics and engineering using ADM and the Modified Decomposition Method (MDM) in collaboration with Drs. Rach and Wazwaz. He has been a professor of mathematics at the Inner Mongolia Polytechnic University, the Tianjin University of Commerce, and currently at the School of Sciences, Shanghai Institute of Technology, China. He is the author of more than eighty journal papers. Dr. Duan is an outstanding player of wei qi (the Chinese “game of go”), and ping-pong (table tennis). He enjoys photography and traveling.

Dr. Randolph Rach is a retired Army veteran and former Senior Engineer at Microwave Laboratories Inc. His experience includes research and development in microwave electronics and traveling-wave tube technology, and his research interests span nonlinear system analysis, nonlinear ordinary and partial differential equations, nonlinear integral equations and nonlinear boundary value problems. He published more than one hundred and thirty papers in applied mathematics, and he is an early contributor to the ADM. Dr. Rach’s fundamental theorems established the basis for the early development of ADM. He was also the first to propose the modified decomposition method (MDM for short). His past and current research continue to advance ADM/MDM work. He prepared the comprehensive bibliography on ADM:
R. Rach, A bibliography of the theory and applications of the Adomian decomposition method, 1961-2011, Kybernetes, Volume 41, Nos. 7/8, 2012, 1087--1148.

Dr. Sergio E. Serrano was born in Santander (Colombia) in 1953. He is a professor of environmental engineering, hydrologic science, applied mathematics, and philosophy at Temple University in Philadelphia. He received his Ph.D. degree at the University of Waterloo (Canada). He has more than one hundred research publications in international science, engineering, and mathematical journals. He is also the author of eight books in environmental engineering, statistics, philosophy, and psychology. Dr. Serrano has been awarded four times with nationally-competitive research grants by the National Science Foundation, Washington, DC. Using adaptations and modifications of the ADM, he has developed hundreds of practical engineering models of flood wave propagation, contaminant transport, and groundwater flow in irregular geometries. He published influential textbooks on applications of differential equations in CAS Maple and Hydrology:

• Sergio E. Serrano, Differential Equations, HydroScience Inc. Amber, Pennsylvania, 2016, ISBN 97809888865211
• Sergio E. Serrano, Hydrology for Engineers, Geologists, and Environmental Professionals. An Integrated Treatment of Surface, Subsurface, and Contaminant Hydrology. HydroScience Inc., Ambler, PA.

Dr. Serrano has a passion for hiking. In 1973, he explored the rain forest of the Darien Sierra on foot from Acandi (Colombia) to Yaviza (Panama). He plays the recorder (Renaissance flute). He enjoys alchemy, archaeology, home wine making, herbal medicine, and cooking. He believes that the joy of meaningful living and meaningful relationships can be found in the simplicity of everyday life. He lives in Philadelphia with his wife of thirty years and his daughter.

Dr. Abdul-Majid Wazwaz is a Professor of Mathematics at Saint Xavier University in Chicago. He received his Ph.D. from the University of Illinois at Chicago. He was the author and co-author of more than four hundred and fifty papers in applied mathematics and mathematical physics. He is the author of five books on the subjects of discrete mathematics, integral equations, and partial differential equations. He has contributed extensively to theoretical advances in solitary waves theory, the ADM, and other computational methods. He is a member of the editorial board of the journals Nonlinear Dynamics (Springer) and Physica Scripta (IOP). For three years in a row, 2014, 2015, and 2016, Thomas Reuters granted him three different badges for being a “Highly Cited Researcher.”

The Adomian decomposition method (ADM for short) has led to several modifications on the method made by various researchers in an attempt to improve the accuracy or expand the application of the original approach. Obviously, this tutorial cannot cover and explain all available improvements for the method. Instead, we show how the Adomian's decomposition method works in a series of examples. The basic spirit of the decomposition method consists of three steps. The first one includes a decomposition of the true solution to the initial value problem (or boundary value problem) into the infinite sum of the solution's components to auxiliary problems. The second step asks to find the solution's components recursively based on previously obtained components, so it involves solution of full-history recurrence. This stage requires solving very basic differential equations with explicitly determined right-hand side expressions that are obtained by decomposition of nonlinear terms into the so called Adomian's polynomials. Finally, summation of first finite number of components leads to the approximation to the true solution with any desired level of high precision.

If the given differential equation is homogeneous (without driving terms), then the ADM finite sum always provides the truncated series version to the true solution. However, for some special classes of inhomogeneous differential equations, the ADM finite sum approximation may include some noise terms that eventually are canceled out with next iterations subject that all iterations are complete. This phenomenon is well-known for Pirad's ietration. However, unwanted terms in the ADM applications were first discovered by G. Adomian and R. Rach in 1986 (The noisy convergence phenomena in decomposition method solutions, Journal of Computational and Applied Mathematics, Volume 15, Issue 3, July 1986, pages 379--381). Although it is possible to reduce these unwanted noise terms, as shown by Abdul-Majid Wazwaz and other researchers, the issue remains. However, these noise terms usually do not effect the approximations because they are rapidly damped out numerically.

The first proof of convergence of the ADM was given by Cherruault in 1989, who used fixed point theorems for abstract functional equations. Since conditions of these theorems are too restrictive for most physical and engineering applications to be verified in practice, many other articles on convergence of the ADM were published. In spite of the variety of publications on convergence, computational complexity, improvements, and applications of the ADM, no precise criterion of convergence was formulated in the literature, at least in the context of initial-value problems for ordinary and partial differential equations. When slope function is a holomorphic function, then the formalism of the Cauchy--Kowalevskaya theorem guarantees that solutions of initial-value problems for systems of ordinary differential equations exist and are analytic for small time intervals. Recent review on this topic can be found in Ray's article.

Let us compare the ADM and Picard's iteration scheme. The latter can be successfully applied only for differential equations with polynomial slope functions. The ADM extends Picard's requirements for analytically defined nonlinearities. The advantage of Adomian’s method over the Picard scheme is the ease of computation of successive terms. At each iterative step, the Adomian decomposition method actually requires solving the same (very simple) initial value problem with homogeneous initial conditions. Since the ADM presents the solution as (infinite) series, it avoids repetitions that slow down Picard's procedure. Although it is a changeable problem to explicitly determine the radius of convergence for the solution obtained on the basis of the ADM, it is usually much larger than that of Picard's prediction.

These advantages of the ADM over Picard's iteration scheme are abated by preprocessing: one should evaluate Adomian's polynomials that are used at every iteration step. The labor involved in such evaluation grows exponentially in general. This drawback becomes negligible in an explicit one step numerical implementation based on the ADM when only a few first terms are taken into account. It should be noted that the amount of work for the ADM preprocessing is larger than, say, the well known Runge--Kutta or cubic spline algorithms where a user enters only the slope functions and the initial conditions. Usually, the labor required by ADM preprocessing is compensated by a larger step size in the computations than in standard one step algorithms.

We start demonstration of the Adomian decomposition method with the following initial value problem:

Khelifa and Cherruault proved in 2000 (Kybernetes, Vol. 29, No. 3, pp. 332--355) a bound for Adomian polynomials:

The initial term u₀ in the series expansion \( y(x) = \sum_{n \ge 0} u_n (x) \) is a solution of the truncated initial value problem:

Adomian's polynomials can be evaluated using the following Duan's code (Applied Mathematical Modelling, 2013, 37, Issues 21-22, pp. 8687--8708):

Example: Consider the initial value problem for the linear equation

\[ y' = x \left( 3 - 2y \right) = f(x,y) + g(x), \qquad y(0) =1, \]

where \( g(x) = 3x \) and \( f(x,y) = -2x\,y . \) In this example, we will find the first 5 terms of the Adomian approximation for the above IVP. We will then compare our approximation with the true solution found with Mathematica

DSolve[{y'[x] == x*(3 - 2*y[x]), y[0] == 1}, y[x], x]// Flatten

{y[x] -> 1/2 E^-x^2 (-1 + 3 E^x^2)}

\[ y= \phi (x) = \frac{3}{2} - \frac{1}{2}\, e^{-x^2} = 1 + \frac{x^2}{2} - \frac{x^4}{4} + \frac{x^6}{12} - \frac{x^8}{2\cdot 4!} + \frac{x^{10}}{2\cdot 5!} + \cdots , \]

in order to see how good our approximation is. According to Adomian, we seek its solution as an infinite sum (which we truncate to keep a few first terms)

\[ y(x) = \sum_{n\ge 0} u_n (x) = u_0 (x) + u_1 (x) + u_2 (x) + \cdots , \]

where the first term u₀ is the solution to the following initial value problem:

\[ u' = g(x) = 3x , \qquad u(0) =1 , \]

because we need to start with the nonhomogeneous part of the original initial value problem. Since this is a calculus problem, its solution is \( u_0 (x) = C + \frac{3\,x^2}{2} , \) where C is a constant of integration which happens to be the initial condition of the IVP. Hence, C=1. This u₀ differs from Picard's starting term, which is 1. With u₀, we can find A₀, the zeroth Adomian polynomial:

\[ A_0 = f(x, u_0 ) = x \left( 3 -2\,u_0 \right) = x - 3\,x^3 . \]

All Adomian polynomials, except the first one, are computed according to a single formula (which is valid only for linear equations)

\[ A_n (x) = -2x\,u_n , \qquad n=1,2,\ldots . \]

Mathematica confirms:

f[u_] := x*(3 - 2*u)
Ado[5]

{x (3 - 2 u[0]), -2 x u[1], -2 x u[2], -2 x u[3], -2 x u[4], -2 x u[ 5]}

Also, we can achieve this with another Mathematica command (which can only be applied for polynomial input):

CoefficientList[-2 x*(u0 + u1*s + u2*s^2 + u3*s^3 + u4*s^4 + u5*s^5), s, 5]

{ -2 u0 x, -2 u1 x, -2 u2 x, -2 u3 x, -2 u4 x}

However, we don't need Adomian polynomials in our simple linear differential equation. We just substitute the infinite series form \( y(x) = \sum_{n\ge 0} u_n (x) \) into the given equation:

\[ \frac{\text d}{{\text d}x} \left[ u_0 + u_1 + u_2 + \cdots \right] = 3x - 2x \left[ u_0 + u_1 + u_2 + \cdots \right] . \]

This leads to a sequence of initial value problems:

\[ \frac{\text d}{{\text d}x} \, u_{n+1} = - 2x \, u_n , \quad u_{n+1} (0) = 0, \qquad n=0,1,2,\ldots . \]

The above differential equation is similar to the recursive equation for Appell polynomials. We ask Mathematica to find the solution to the above IVP.

DSolve[{y'[x] == -2*x*(1 +3* x^2 /2), y[0] == 0}, y[x], x]// Flatten
u[1][x_] = y[x] /. %

{y[x] -> 1/4 (-4 x^2 - 3 x^4)}

Hence, \( u_1 = - x^2 - 3\,x^4 /4 \) and we get the first order approximation

\[ y_1 = u_0 + u_1 = 1 + \frac{3\,x^2}{2} - x^2 - \frac{3\,x^4}{4} = 1 + \frac{1}{2}\,x^2 - \frac{3\,x^4}{4} , \]

phi1[x_] = Simplify[1 + 3*x^2/2 + u[1][x]]

1 + x^2 /2 - 3 x^4 /4

which differs from Picard's one, \( \phi_2 = 1 + \frac{x^2}{2} - \frac{x^4}{4} \) in the coefficient of the largest power of x. Next, we find u₂ by solving the initial value problem

\[ y' = A_1 = -2x \, u_1 = 2\,x^3 + 3\,x^5 /2, \quad y(0) =0 . \]

So we ask Mathematica again to do this job for us:

DSolve[{y'[x] == 2*x^3 + 3* x^5 /2, y[0] == 0}, y[x], x] //Flatten

{y[x] -> 1/4 (2 x^4 + x^6)}
u[2][x_] = y[x] /. %

Adding to the previous sum, we get the next approximation

phi2[x_] = Simplify[phi1[x] + u[2][x]]

1/4 (4 + 2 x^2 - x^4 + x^6)

Similarly,

DSolve[{y'[x] == -2*x*u[2][x], y[0] == 0}, y[x], x] //Flatten
u[3][x_] = y[x] /. %

{y[x] -> 1/48 (-8 x^6 - 3 x^8)}

So we get

phi3[x_] = Simplify[phi2[x] + u[2][x]]

1/48 (48 + 24 x^2 - 12 x^4 + 4 x^6 - 3 x^8)

\[ u_0 = 1 + \frac{3\,x^2}{2} , \quad u_1 = -x^2 - \frac{3}{4}\, x^4 , \quad u_2 = \frac{1}{2} \,x^4 + \frac{1}{4}\, x^6 , \quad u_3 = -\frac{1}{6}\,x^6 - \frac{1}{16}\, x^8 , \quad u_4 = \frac{1}{24} \,x^8 + \frac{1}{80}\,x^{10} , \quad u_5 = -\frac{1}{120} \,x^{10} - \frac{1}{480}\, x^{12} . \]

Adding these functions, we obtain the 5-term Adomian approximation

\[ \phi_5 (x) = u_0 + u_1 + u_2 +u_3 +u_4 +u_5 = 1+ \frac{x^2}{2} - \frac{x^4}{4} + \frac{1}{12}\,x^6 - \frac{1}{48}\, x^8 + \frac{1}{240}\, x^{10} - \frac{1}{480}\, x^{12} . \]

Now we plot the exact solution along with our 5-th order Adomian approximation.

Plot[{Callout[3/2-Exp[-x^2]/2, "exact"], Callout[phi5[x], "ADM approximation"]}, {x,0,2.0}, PlotStyle->Thick]

   ■

Example: Consider the initial value problem for the Riccati equation

\[ y' = x^2 + y^2 , \qquad y(0) =1. \]

Its explicit solution is known to be (see Riccati section) according to Mathematica

riccati = DSolve[{y'[x] == x^2 + (y[x])^2, y[0] == 1}, y[x], x]; phi[x_] = y[x] /. riccati

\[ y(x) = \frac{x^2 J_{5/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{3}{4} \right) - x^2 J_{1/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{1}{4} \right) + J_{1/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{3}{4} \right) - x^2 J_{3/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{3}{4} \right) }{x\, J_{1/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{1}{4} \right) -2\,J_{-1/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{3}{4} \right) } , \]

DSolve[{y'[x] == x^2 + (y[x])^2, y[0] == 1}, y[x], x]

where J_ν(z) is the Bessel function of the first kind and \( \Gamma (z) = \int_0^{\infty} e^{-t} t^{z-1}\,{\text d}t \) is the gamma function of Euler.

This function blows up when the denominator is zero. This point can be determined approximately from the following graph.

y[x_] = (-x^2 BesselJ[-(3/4), x^2/2] Gamma[1/4] + x^2 BesselJ[-(5/4), x^2/2] Gamma[3/4] + BesselJ[-(1/4), x^2/2] Gamma[3/4] - x^2 BesselJ[3/4, x^2/2] Gamma[3/ 4])/(x (BesselJ[1/4, x^2/2] Gamma[1/4] - 2 BesselJ[-(1/4), x^2/2] Gamma[3/4]))
deno[x_] = (x (BesselJ[1/4, x^2/2] Gamma[1/4] - 2 BesselJ[-(1/4), x^2/2] Gamma[3/4]))
Plot[{y[x], deno[x]}, {x, 0, 3}, PlotStyle -> {{Thick, Blue}, {Thick, Red}}]

By plotting the denominator of Ricatti solution, we find approximately the point where its graph (in red) crosses the horizontal axis. This point indicates that the position where the solution to the Riccati equation blows up. It is impossible to determine this point from the initial value problem only. Fortunately, Mathematica can help with its determination.

FindRoot[deno[x], {x, 1}]

{x -> 0.969811}

Therefore, we know now that the solution (which is unique) to the given initial value problem exists only within the interval [0, 0.9698).

Our next step is to determine a power series representation for the solution. We can find its 10-term Maclaurin series expansion with one line Mathematica command:

AsymptoticDSolveValue[{y'[x] == x^2 + (y[x])^2, y[0] == 1}, y[x], {x, 0, 10}]

1 + x + x^2 + (4 x^3)/3 + (7 x^4)/6 + (6 x^5)/5 + (37 x^6)/30 + ( 404 x^7)/315 + (369 x^8)/280 + (428 x^9)/315 + (1961 x^10)/1400

\[ y(x) = 1 + x + x^2 + \frac {4} {3}\, x^3 + \frac {7} {6}\, x^4 + \frac{6}{5}\, x^5 + \frac{37}{30}\, x^6 + \frac{404}{315}\, x^7 + \cdots . \]

At this point, although we know many terms in its Maclaurin series, it is impossible to find its radius of convergence because we don't know how coefficients grow.

According to G. Adomian, we seek the solution as an infinite sum

\[ y(x) = \sum_{n\ge 0} u_n , \]

where u₀ is the unique solution of the following IVP:

\[ y' = x^2 , \qquad y(0) =1; \]

which includes only the nonhomogeneous components of the original IVP. Using Mathematica, we obtain

DSolve[{y'[x] == x^2 , y[0] == 1}, y[x], x]
u[0][x_] = y[x] /. %

u[0][x_]:= 1/3 (3 + x^3)

So

\[ u_0 = 1 + \frac{1}{3}\, x^3 . \]

Since the problem is nonlinear, we need Adomian's polynomials to expand the nonlinear term (y²):

f[u_] := u*u
Ado[5]

{u[0]^2, 2 u[0] u[1], u[1]^2 + 2 u[0] u[2], 2 u[1] u[2] + 2 u[0] u[3], 2 (u[2]^2/2 + u[1] u[3]) + 2 u[0] u[4], 2 (u[2] u[3] + u[1] u[4]) + 2 u[0] u[5]}

We can calculate the Adomian polynomials manually upon substituting the series \( y= \sum_n u_n \) into the differential equation:

\[ \frac{{\text d}}{{\text d}x} \left[ u_0 + u_1 + u_2 + \cdots \right] = x^2 + \left[ u_0 + u_1 + u_2 + \cdots \right]^2 = x^2 + A_0 + A_1 + A_2 + \cdots . \]

In our case of the standard Riccati equation, we have

\begin{eqnarray*} A_0 &=& u_0^2 = \frac{1}{9} \left( 3 + x^3 \right)^2 , \\ A_1 &=& 2\,u_1 \,u_0 = \frac{2}{3} \left( 3 + x^3 \right) u_1 , \\ A_2 &=& 2\,u_2 \,u_0 + u_1^2 = \frac{2}{3} \left( 3 + x^3 \right) u_2 + u_1^2 , \\ A_3 &=& 2\,u_3 \,u_0 + 2\,u_1 u_2 = \frac{2}{3} \left( 3 + x^3 \right) u_3 + 2\,u_1 u_2 , \\ A_4 &=& 2\,u_4 \,u_0 + u_2^2 + 2\, u_1 u_3 = \frac{2}{3} \left( 3 + x^3 \right) u_4 + u_2^2 + 2\, u_1 u_3 , \\ A_5 &=& 2\,u_5 \,u_0 + 2 \left( u_1 u_4 + u_2 u_3 \right) , \end{eqnarray*}

and so on. Then we calculate components in the Adomian series solution \( y= \sum_n u_n \) by solving a sequence of initial value problems:

\[ \frac{\text d}{{\text d}x}\, u_{n+1} = A_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_{n+1} (0) =0 . \]

Since

\[ A_0 (x) = \left( 1 + \frac{1}{3}\, x^3 \right)^2 , \]

we find u₁ by solving the initial value problem

\[ u'_1 = A_0 (x) = \left( 1 + \frac{1}{3}\, x^3 \right)^2 , \qquad u_1 (0) =0 . \]

A[0][x_] = Expand[(1 + x^3 /3!)^2]

1 + (2 x^3)/3 + x^6/9

DSolve[{y'[x] == A[0][x], y[0] == 0}, y[x], x] // Flatten
u[1][x_] = y[x] /. %

1/126 (126 x + 21 x^4 + 2 x^7)

Next iteration:

A[1][x_] := Simplify[2*u[0][x]*u[1][x]]

2 x + x^4 + x^7/7 + (2 x^10)/189

DSolve[{y'[x] == A[1][x], y[0] == 0}, y[x], x] // Flatten
u[2][x_] := y[x] /. %

(83160 x^2 + 16632 x^5 + 1485 x^8 + 80 x^11)/83160

phi2[x_] := Simplify[u[0][x] + u[1][x] + u[2][x]]

1 + x + x^2 + x^3/3 + x^4/6 + x^5/5 + x^7/63 + x^8/56 + (2 x^11)/2079

Comparing the two-term Adomian approximation φ₂(x) with the true Maclaurin expansion, we see that φ₂(x) has correct first three terms and the noise terms:

\[ \phi_2 (x) = 1 + x + x^2 + \underbrace{\frac{x^3}{3} + \frac{x^4}{6} + \frac{x^5}{5} + \frac{x^7}{63} + \frac{x^8}{56} + \frac{2\,x^{11}}{2079}}_{\mbox{noise terms}} . \]

Next Adomian polynomial is found with the aid of Mathematica

A2[x_] = Simplify[2*u0[x]*u2[x] + (u1[x])^2]

3 x^2 + (7 x^5)/5 + (8 x^8)/35 + (53 x^11)/2772 + (13 x^14)/14553

DSolve[{y'[x] == A2[x], y[0] == 0}, y[x], x] // Flatten

(3492720 x^3 + 814968 x^6 + 88704 x^9 + 5565 x^12 + 208 x^15)/3492720

u[3][x_] = (3492720 x^3 + 814968 x^6 + 88704 x^9 + 5565 x^12 + 208 x^15)/3492720

phi3[x_] = Simplify[phi2[x] + u3[x]]

1 + x + x^2 + (4 x^3)/3 + x^4/6 + x^5/5 + ( 7 x^6)/30 + x^7/63 + x^8/56 + (8 x^9)/315 + (2 x^11)/2079 + ( 53 x^12)/33264 + (13 x^15)/218295

The three term approximation \( \phi_3 (x) = u_0 + u_1 + u_2 + u_3 \) is

\begin{align*} \phi_3 (x) &= u_0 + u_1 + u_2 + u_3 = 1 + x + x^2 + \frac{4}{3}\, x^3 \\ &\quad + \underbrace{\frac{x^4}{6} + \frac{x^5}{5} + \frac{7}{30}\, x^6 + \frac{x^7}{63} + \frac{1}{56}\, x^8 + \frac{8}{315}\, x^9 + \frac{2}{2079}\, x^{11} + \frac{53}{33264}\, x^{12} + \frac{13}{218295}\,x^{15}}_{\mbox{noise terms}} . \end{align*}

For the next term in Adomian's series, we need to find A₃:

A3[x_] = Simplify[2*u0[x]*u3[x] + 2*u1[x]*u2[x]]

4 x^3 + (28 x^6)/15 + (143 x^9)/420 + (17 x^12)/495 + ( 3613 x^15)/1746360 + (46 x^18)/654885

DSolve[{y'[x] == A3[x], y[0] == 0}, y[x], x] // Flatten
u4[x_] = 1/ 103524220800 (103524220800 x^4 + 27606458880 x^7 + 3524753232 x^10 + 273490560 x^13 + 13386165 x^16 + 382720 x^19)

Then we calculate the corresponding sum of first four terms:

phi4[x_] = Simplify[phi3[x] + u4[x]]

1 + x + x^2 + (4 x^3)/3 + (7 x^4)/6 + x^5/5 + (7 x^6)/30 + ( 89 x^7)/315 + x^8/56 + (8 x^9)/315 + (143 x^10)/4200 + ( 2 x^11)/2079 + (53 x^12)/33264 + (17 x^13)/6435 + ( 13 x^15)/218295 + (3613 x^16)/27941760 + (46 x^19)/12442815

\begin{align*} \phi_4 (x) &= \phi_3 + u_4 = 1 + x + x^2 + \frac{4}{3}\, x^3 + \frac{7}{6}\, x^4 \\ & \quad + \underbrace{\frac{x^5}{5} + \frac{7}{30}\, x^6 + \frac{89}{315}\, x^7 + \frac{x^8}{56} + \frac{8}{315}\, x^9 + \frac{143}{4200}\, x^{10} + \frac{2}{2079} \, x^{11} + \frac{53}{33264} \, x^{12} + \frac{17}{6435}\, x^{13} + \frac{13}{218295}\, x^{15} + \frac{3613}{27941760} \, x^{16} + \frac{46}{12442815} \,x^{19}}_{\mbox{noise terms}} . \end{align*}

This series reveals that the first terms up to the power 4 are correct. We next compare the last Adomian approximation with the true solution by calculating the error function.

Plot[{Callout[y[x], "exact solution"], Callout[phi4[x], "ADM 4-approximation"]}, {x, 0, 0.7}, PlotStyle -> Thick]
error[x_] = y[x] - phi4[x]
Plot[{error[x]}, {x, 0, 0.7}, PlotStyle -> {Thick}]

We compare our ADM approximation with true solution and approximation obtained with the aid of classical Runge--Kutta of order 4 (abbreviated by RK4) with step size h = 0.1.

Point	True value	ADM-4 approximation	RK4 approximation
x = 0.1	1.11146	1.11145	1.11146
x = 0.2	1.25302	1.25262	1.25302
x = 0.3	1.43967	1.43617	1.43967

Expantion of the true solution into the power series can be done directly:

Clear[seriesDSolve];
seriesDSolve[ode_, y_, iter : {x_, x0_, n_}, ics_: {}] :=
Module[{dorder, coeff},
dorder = Max[0, Cases[ode, Derivative[m_][y][x] :> m, Infinity]];
(coeff =
Fold[Flatten[{#1,
Solve[#2 == 0 /. #1,
Union@Cases[#2 /. #1, Derivative[m_][y][x0] /; m >= dorder,
Infinity]]}] &, ics,
Table[SeriesCoefficient[
ode /. Equal -> Subtract, {x, x0, k}], {k, 0, Max[n - dorder, 0]}]];
Series[y[x], iter] /. coeff) /; dorder > 0]

seriesDSolve[y'[x] == x^2 + (y[x])^2, y, {x, 0, 15}, {y[0] -> 1} ]

\( 1 + x + x^2 + \frac{4\,x^3}{3} + \frac{7\,x^4}{6} + \frac{6\,x^5}{5} + \frac{37\,x^6}{30} + \frac{404\,x^7}{315} + \frac{369\,x^8}{315} + \frac{1961\,x^{10}}{1400} + \)
\( \frac{75092\,x^{11}}{51975} + \frac{1238759\,x^{12}}{831600} + \frac{9884\,x^{13}}{6435} + \frac{17121817\,x^{14}}{10810800} + \frac{115860954\,x^{15}}{70945875} + O[x]^{16} \)

Then we plot the true solution and its Adomian approximation with N=2 terms

sol = (y[x] /.
DSolve[{y'[x] == x^2 + (y[x])^2, y[0] == 1}, y[x], x])[[1]]
Plot[{Callout[sol, "exact"], Callout[phi2[x], "ADM 2-approximation"], Callout[phi3[x], "ADM 3-approximation"]}, {x,0,1.5}, PlotStyle->Thick]

As we see from this graph, we can get a good approximation within a small interval [0,0.2]. To get a better approximation, however, we need to to add more terms. The solution has a singular point near x = 1, which we find by equating the denominator to zero in explicit representation of the true solution:

FindRoot[(x (BesselJ[1/4, x^2/2] Gamma[1/4] - 2 BesselJ[-(1/4), x^2/2] Gamma[3/4])) == 0, {x, 1}]

{x -> 0.969811}

Therefore, the true solution blows up around the point 0.969811, but power series approximation cannot detect the exact point no matter how many terms we use.

 

Not standard ADM

---

We show how the ADM works when we split the initial condition between several few first Adomian components. So we choose the initial approximation as the solution of the following initial value problem:

\[ u' = x^2 , \qquad u(0) = 1/3 . \]

Its solution is

DSolve[{y'[x] == x^2 , y[0] == 1/3}, y[x], x]

\[ u_0 (x) = \frac{1}{3} \left( 1 + x^3 \right) . \]

Since all Adomian polynomails were determined previously, we find next two components in the Adomian series by solving the following initial value problems:

DSolve[{y'[x] == (1+x^3)^2 /9, y[0] == 1/3}, y[x], x] // Flatten
u[1][x_] = y[x] /. %

1/126 (42 + 14 x + 7 x^4 + 2 x^7)

u[0][x_]= (1+x^3)/3
A[1][x_] := Simplify[2*u[0][x]*u[1][x]]
DSolve[{y'[x] == A[1][x], y[0] == 1/3}, y[x], x] // Flatten
u[2][x_] = y[x] /. %

Then two-term Adomian approximation becomes

phi2[x_] = Simplify[u[0][x] + u[1][x] + u[2][x]]

1 + x/3 + x^2/27 + x^3/3 + x^4/9 + x^5/45 + x^7/63 + x^8/168 + ( 2 x^11)/2079

\[ \phi_2 (x) = u_0 + u_1 + u_2 = 1 + \frac{x}{3} + \frac{x^2}{27} + \frac{x^3}{3} + \frac{x^4}{9} + \frac{x^5}{45} + \frac{x^7}{63} + \frac{x^8}{168} + \frac{2\,x^{11}}{2079} . \]

Next componet will be

A[2][x_] = Simplify[2*u[0][x]*u[2][x] + (u[1][x])^2]
DSolve[{y'[x] == A[2][x], y[0] == 0}, y[x], x] // Flatten
u[3][x_] = y[x] /. %

(1/3492720)(1164240 x + 388080 x^2 + 43120 x^3 + 194040 x^4 + 155232 x^5 + 30184 x^6 + 20790 x^8 + 9856 x^9 + 1855 x^12 + 208 x^15)

phi3[x_] = Simplify[phi2[x] + u[3][x]]

1 + (2 x)/3 + (4 x^2)/27 + (28 x^3)/81 + x^4/6 + x^5/15 + ( 7 x^6)/810 + x^7/63 + x^8/84 + (8 x^9)/2835 + (2 x^11)/2079 + ( 53 x^12)/99792 + (13 x^15)/218295

A[3][x_] = Simplify[2*u[0][x]*u[3][x] + 2*u[1][x]*u[2][x]]
DSolve[{y'[x] == A[3][x], y[0] == 0}, y[x], x] // Flatten
u[4][x_] = y[x] /. %
phi4[x_] = Simplify[phi3[x] + u[4][x]]

1 + (8 x)/9 + (10 x^2)/27 + (32 x^3)/81 + (83 x^4)/486 + ( 2 x^5)/15 + (14 x^6)/405 + (163 x^7)/8505 + x^8/56 + (8 x^9)/945 + ( 143 x^10)/113400 + (2 x^11)/2079 + (53 x^12)/49896 + ( 17 x^13)/57915 + (13 x^15)/218295 + (3613 x^16)/83825280 + ( 46 x^19)/12442815

\begin{align*} \phi_4 (x) &= u_0 + u_1 + u_2 + u_3 + u_4 \\ &= 1 + \frac{8\,x}{9} + \frac{10\,x^2}{27} + \frac{32\,x^3}{81} + \frac{83\,x^4}{486} + \frac{2\,x^5}{15} \\ & \quad + \frac{14\,x^6}{405} + \frac{163\,x^7}{8505} + \frac{x^8}{56} + \frac{8\,x^9}{945} + \cdots . \end{align*}

So we see that the standard ADM gives better approximation than the above algorithm.

 

Accelerated Adomian polynomials

---

We perform the ADM approximation based on the accelerated Adomian polynomials. The initial Adomian polynomails is the same:

u[0][x_] = (1 + x/3)^3 DSolve[{y'[x] == (1 + x^3/3)^2, y[0] == 0}, y[x], x] // Flatten u[1][x_] = y[x] /. %

1/126 (126 x + 21 x^4 + 2 x^7)}

J[1][x_]=(u[0][x]+u[1][x])^2 - (u[0][x])^2 DSolve[{y'[x] == J[1][x], y[0] == 0}, y[x], x] // Flatten u[2][x_] = y[x] /. %

(1/5239080)(5239080 x^2 + 1746360 x^3 + 1047816 x^5 + 291060 x^6 + 93555 x^8 + 34650 x^9 + 5040 x^11 + 2310 x^12 + 88 x^15)

And the third component becomes

J[2][x_] = (u[0][x] + u[1][x] + u[2][x])^2 - (u[0][x] + u[1][x])^2 DSolve[{y'[x] == J[1][x], y[0] == 0}, y[x], x] // Flatten u[3][x_] = y[x] /. %

(1/5239080)(5239080 x^2 + 1746360 x^3 + 1047816 x^5 + 291060 x^6 + 93555 x^8 + 34650 x^9 + 5040 x^11 + 2310 x^12 + 88 x^15)

phi3[x_] = Simplify[u[0][x]+u[1][x]+u[2][x] + u[3][x]]

1 + x + 2 x^2 + x^3 + x^4/6 + (2 x^5)/5 + x^6/9 + x^7/63 + x^8/28 + ( 5 x^9)/378 + (4 x^11)/2079 + x^12/1134 + (2 x^15)/59535

J[3][x_] = (u[0][x] + u[1][x] + u[2][x] + u[3][x])^2 - (u[0][x] + u[1][x] + u[2][x])^2 DSolve[{y'[x] == J[1][x], y[0] == 0}, y[x], x] // Flatten u[4][x_] = y[x] /. %

phi4[x_] = Simplify[phi3[x] + u[4][x]]

1 + x + 3 x^2 + (4 x^3)/3 + x^4/6 + (3 x^5)/5 + x^6/6 + x^7/63 + ( 3 x^8)/56 + (5 x^9)/252 + (2 x^11)/693 + x^12/756 + x^15/19845

   ■

Example: Consider the initial value problem for the Riccati equation

\[ y' = x^2 - y^2 , \qquad y(0) =1. \]

We can use Mathematica to find its solution.

DSolve[{y'[x] == x^2 - (y[x])^2, y[0] == 1}, y[x], x]

{{y -> Function[{x}, ((1/2 + I/ 2) ((1 + I) x^2 BesselJ[-(3/4), (I x^2)/2] Gamma[1/4] + I Sqrt[2] x^2 BesselJ[-(5/4), (I x^2)/2] Gamma[3/4] + Sqrt[2] BesselJ[-(1/4), (I x^2)/2] Gamma[3/4] - I Sqrt[2] x^2 BesselJ[3/4, (I x^2)/2] Gamma[3/4]))/(x (BesselJ[1/4, ( I x^2)/2] Gamma[1/4] + (1 + I) Sqrt[2] BesselJ[-(1/4), (I x^2)/2] Gamma[3/4]))]}}

This solution is not yet useful because it contains the imaginary unit. We need to transfer it to a real value expression.

\[ y(x) = x\,\frac{I_{-3/4} \left( \frac{x^2}{2} \right) \pi \left( \Gamma^2 \left( \frac{3}{4} \right) \sqrt{2} + \pi \right) - 2\,K_{3/4} \left( \frac{x^2}{2} \right) \Gamma^2 \left( \frac{3}{4} \right)}{\pi \left( \Gamma^2 \left( \frac{3}{4} \right) \sqrt{2} + \pi \right) I_{1/4} \left( \frac{x^2}{2} \right) + 2\,K_{1/4} \left( \frac{x^2}{2} \right) \Gamma^2 \left( \frac{3}{4} \right)} \]

where I_ν(z) is the modified Bessel function of the first kind and K_ν(z) is the modified Bessel function of the second kind (alternatively, it is also called after William Thomson, 1st Baron Kelvin); here \( \Gamma (\nu ) = \int_0^{\infty} x^{\nu -1} e^{-x}\,{\text d}x \) is the gamma function.

The denominator in the above formula is never zero, as its graph shows. Therefore, the solution to the given IVP exists for all positive x.

deno[x_] = Pi*(Sqrt[2]*(Gamma[3/4])^2 + Pi)* BesselI[1/4, x^2 /2] + 2*(Gamma[3/4])^2 * BesselK[1/4, x^2 /2]
Plot[deno[x], {x, 0, 1.5}, PlotStyle -> Thick, PlotRange -> {{0.01, 1.4}, {14, 25}}]

The Maclaurin power series expansion of the solution can be found with Mathematica:

AsymptoticDSolveValue[{y'[x] == x^2 - (y[x])^2, y[0] == 1}, y[x], {x, 0, 10}]

1 - x + x^2 - (2 x^3)/3 + (5 x^4)/6 - (4 x^5)/5 + (23 x^6)/30 - ( 236 x^7)/315 + (201 x^8)/280 - (218 x^9)/315 + (2803 x^10)/4200

\[ y(x) = 1-x + x^2 - \frac{2}{3}\,x^3 + \frac{5}{6}\,x^4 - \frac{4}{5}\, x^5 + \frac{23}{30}\,x^6 - \frac{236}{315}\, x^7 + \frac{201}{280}\, x^8 - \frac{218}{315}\,x^9 + \frac{2803}{4200}\,x^{10} - \frac{33412}{51975}\,x^{11} + \cdots . \]

According to G. Adomian, we seek the solution in a way similarly to the previous example:

\[ y(x) = \sum_{n\ge 0} u_n , \]

where the initial term u₀ is exactly the same: \( u_0 = 1 + x^3 /3 . \) Next, we find five Adomian polynomials using Mathematica

f[u_] := - u*u;
Ado[5]

{x^2 - u[0]^2, -2 u[0] u[1], -u[1]^2 - 2 u[0] u[2], -2 u[1] u[2] - 2 u[0] u[3], -2 (u[2]^2/2 + u[1] u[3]) - 2 u[0] u[4], -2 (u[2] u[3] + u[1] u[4]) - 2 u[0] u[5]}

\begin{align*} A_0 &= - u_0^2 , \\ A_1 &= -2\,u_0 u_1 , \\ A_2 &= -u_1^2 -2\,u_0 u_2 , \\ A_3 &= -2\, u_1 u_2 - 2\,u_0 u_3 , \\ A_4 &= - u_2^2 - 2\,u_1 u_3 -2\,u_0 u_4 , \\ A_5 &= -2\,u_2 u_3 - 2\, u_1 u_4 -2\,u_0 u_5 . \end{align*}

We proceed in exactly the same way as in the previous example. Now we are ready to determine the first five terms in Adomian series by solving the recursive initial value problems:

\[ \frac{{\text d} u_n}{{\text d}x} = A_{n-1} \left( u_0 , u_1 , \ldots , u_{n-1} \right) , \qquad u_n (0) =0 , \qquad n=1,2,\ldots . \]

u[0][x_] = 1 + x^3 /3 ;
DSolve[{y'[x] == - (u[0][x])^2, y[0] == 0}, y[x], x] // Flatten ;
u[1][x_] = y[x] /. %

{1/126 (-126 x + 42 x^3 - 21 x^4 - 2 x^7)}

\[ u_1 (x) = \frac{1}{126} \left( -126 x + 42 x^3 - 21 x^4 - 2 x^7 \right) . \]

The second term is

DSolve[{y'[x] == -2*u[0][x]*u[1][x], y[0] == 0}, y[x], x] // Flatten;
u[2][x_] = y[x] /. %

(83160 x^2 - 13860 x^4 + 16632 x^5 - 2640 x^7 + 1485 x^8 + 80 x^11)/83160

So

\[ u_2 (x) = \frac{1}{83160} \left( 83160 x^2 - 13860 x^4 + 16632 x^5 - 2640 x^7 + 1485 x^8 + 80 x^{11} \right) . \]

Now we add the first two terms:

phi2[x_] = Simplify[u[0][x] + u[1][x] + u[2][x]]

1 - x + x^2 + (2 x^3)/3 - x^4/3 + x^5/5 - x^7/21 + x^8/56 + ( 2 x^11)/2079

Next approximation:

DSolve[{y'[x] == -2*u[0][x]*u[2][x] - u[1][x]*u[1][x], y[0] == 0}, y[x], x] // Flatten;
u[3][x_] = y[x] /. %

(1/3492720)(-3492720 x^3 + 698544 x^5 - 814968 x^6 - 55440 x^7 + 124740 x^8 - 88704 x^9 + 10080 x^11 - 5565 x^12 - 208 x^15)

Adding to ϕ₂, we obtain

phi3[x_] = Simplify[phi2[x] + u[3][x]]

1 - x + x^2 - (2 x^3)/3 - x^4/6 + x^5/5 - ( 7 x^6)/30 - x^7/63 + x^8/56 - (8 x^9)/315 + (2 x^11)/2079 - ( 53 x^12)/33264 - (13 x^15)/218295

Hence,

\begin{align*} \phi_3 (x) &= u_0 (x) + u_1 (x) + u_2 (x) +u_3 (x) = 1 - x + x^2 - \frac{2\,x^3}{3} \\ & \quad - \underbrace{\frac{x^4}{6} + \frac{x^5}{5} - \frac{7\, x^6}{30} - \frac{x^7}{63} + \frac{x^8}{56} - \frac{8\, x^9}{315} + \frac{2\, x^{11}}{2079} - \frac{53\, x^{12}}{33264} - \frac{13\, x^{15}}{218295}}_{\mbox{noise terms}} . \end{align*}

Then we plot the true solution along with Adomian's approximation:

exact[x_] = Simplify[y[x] /. DSolve[{y'[x] == x*x - y[x]*y[x], y[0] == 1}, y[x], x]];
Plot[{Callout[sol, "exact"], Callout[phi2[x], "ADM 2-approximation"], Callout[phi3[x], "ADM 3-approximation"]}, {x,0,2.0}, PlotStyle->Thick]

As we see, the Adomian approximation with three terms ϕ₃ matches the true solution up to power 3, and the rest are noise terms that contain x powers greater than 3. We eliminate some noise terms with the next iteration.

DSolve[{y'[x] == -2*u[1][x]*u[2][x] - 2*u[0][x]*u[3][x], y[0] == 0}, y[x], x] // Flatten;
u[4][x_] = y[x] /. %

(1/103524220800)(103524220800 x^4 + 27606458880 x^7 + 3524753232 x^10 + 273490560 x^13 + 13386165 x^16 + 382720 x^19)

The next Adomian approximation becomes

phi4[x_] = Simplify[phi3[x] + u[4][x]]

1 - x + x^2 - (2 x^3)/3 + (5 x^4)/6 + x^5/5 - (7 x^6)/30 + ( 79 x^7)/315 + x^8/56 - (8 x^9)/315 + (143 x^10)/4200 + ( 2 x^11)/2079 - (53 x^12)/33264 + (17 x^13)/6435 - ( 13 x^15)/218295 + (3613 x^16)/27941760 + (46 x^19)/12442815

which gives the correct approximation up to the power 4 that are accompanied by noise terms of powers greater than 4. The next term in Adomian series

DSolve[{y'[x] == -2*u[1][x]*u[2][x] - 2*u[0][x]*u[3][x], y[0] == 0}, y[x], x] // Flatten;;
u[5][x_] = y[x] /. %

(1/18700822293753600)(-18700822293753600 x^5 - 5610246688126080 x^8 - 819274119535872 x^11 - 74948248743984 x^14 - 4598143925760 x^17 - 185499373059 x^20 - 4293552640 x^23)

The fifth ADM approximation becomes

phi5[x_] = Simplify[phi4[x] + u[5][x]]

1 - x + x^2 - (2 x^3)/3 + (5 x^4)/6 - (4 x^5)/5 - (7 x^6)/30 + ( 79 x^7)/315 - (79 x^8)/280 - (8 x^9)/315 + (143 x^10)/4200 - ( 2227 x^11)/51975 - (53 x^12)/33264 + (17 x^13)/6435 - ( 43327 x^14)/10810800 - (13 x^15)/218295 + (3613 x^16)/27941760 - ( 1318 x^17)/5360355 + (46 x^19)/12442815 - (7523 x^20)/758419200 - ( 15178 x^23)/66108676095

We compare our ADM approximations with the true solution and the approximation obtained with the aid of classical Runge--Kutta of order 4 (abbreviated by RK4) with step size h = 0.1.

numerator[x_] = x*((Sqrt[2]*(Gamma[3/4])^2 + Pi)*BesselI[-3/4, x^2/2]*Pi - 2*(Gamma[3/4])^2*BesselK[3/4, x^2/2])
deno[x_] = Pi*(Sqrt[2]*(Gamma[3/4])^2 + Pi)*BesselI[1/4, x^2/2] + 2*(Gamma[3/4])^2*BesselK[1/4, x^2/2]
phi[x_] = numerator[x]/deno[x]
Clear[y]; h = 0.1;
f[x_, y_] := x^2 - y^2
y[0] = 1;
Do[k1 = h f[h n, y[n]];
k2 = h f[h (n + .5), y[n] + k1/2];
k3 = h f[h (n + .5), y[n] + k2/2];
k4 = h f[h (n + 1), y[n] + k3];
y[n + 1] = y[n] + 1/6*(k1 + 2 k2 + 2 k3 + k4), {n, 0, 4}]
y[1]

0.835786

Point	True value	ADM-4 approximation	ADM-5 approximation	RK4 approximation
x = 0.1	0.909409	0.909418	0.909408	0.90941
x = 0.2	0.835785	0.836052	0.835732	0.835732
x = 0.3	0.777238	0.779122	0.776672	0.777238

So we see that with 5-term ADM approximation is less accurate that the Runge--Kutta algorithm of order four; Also, the fomer requires more labor work compared to the RK-4.    ■

Example: Consider the initial value problem for the autonomous equation

\[ y' = -\sqrt{1/4 + y^2} , \qquad y(0) =3. \]

Its solution exists for all x and is given by

\[ y(x) = -\frac{1}{2}\,\sinh \left( x- \mbox{arcsinh}(6) \right) . \]

Its series representation is

\[ y(x) = 3 - \frac{\sqrt{37}}{2}\,x + \frac{3\,x^2}{2} - \frac{\sqrt{37}}{12}\,x^3 + \frac{x^4}{8} - \frac{\sqrt{37}}{240}\,x^5 + \frac{x^6}{240} - \frac{\sqrt{37}}{10080}\,x^7 + \frac{x^8}{13440} - \frac{\sqrt{37}}{725760}\,x^9 + \frac{x^{10}}{1209600} - \frac{\sqrt{37}}{79833600}\,x^{11} + \cdots , \]

because

exact[x_] = -Sinh[x - ArcSinh[6]]/2 ;
Series[exact[x], {x, 0, 12}]

According to the ADM, we assume that the solution can be represented by a convergent series

\[ y(x) = \sum_{n\ge 0} u_n (x) , \]

where the initial approximation is taken to be u₀ = 3 because the given differential equation is autonomous. Substituting the infinite series into the given differential equation, we obtain

\[ \texttt{D} \, \sum_{n\ge 0} u_n (x) = \sum_{n\ge 0} \texttt{D} \, u_n (x) = \sum_{n\ge 0} A_n , \qquad\mbox{with} \quad \texttt{D} = \frac{\text d}{{\text d}x} , \]

where the Adomian's polynomials are

\begin{eqnarray*} A_0 &=& - \sqrt{\frac{1}{4} + u_0^2} = - \sqrt{\frac{1}{4} + 9} = -\frac{\sqrt{37}}{2} , \\ A_1 &=& - \frac{u_0 u_1}{\sqrt{1/4 + u_0^2}} = -\frac{6}{\sqrt{37}}\,u_1 , \\ A_2 &=& \frac{1}{2} \left[ \frac{u_0^2}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{1}{\sqrt{1/4 + u_0^2}} \right] u_1^2 - \frac{u_0 u_2}{\sqrt{1/4 + u_0^2}} = -\frac{6}{\sqrt{37}}\,u_2 - \frac{1}{37\sqrt{37}}\,u_1^2 , \\ A_3 &=& \frac{1}{2} \left[ \frac{u_0}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{u_0^3}{\left( 1/4 + u_0^2 \right)^{5/2}} \right] u_1^3 + \left[ \frac{u_0^2}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{1}{\sqrt{1/4 + u_0^2}} \right] u_1 u_2 - \frac{u_0 u_3}{\sqrt{1/4 + u_0^2}} \\ &=& -\frac{6}{\sqrt{37}}\,u_3 - \frac{2}{37\sqrt{37}}\,u_1 u_2 + \frac{12}{1369\sqrt{37}}\, u_1^3 , \\ A_4 &=& \frac{1}{8} \left[ \frac{5\,u_0^4}{\left( 1/4 + u_0^2 \right)^{7/2}} - \frac{6\, u_0^2}{\left( 1/4 + u_0^2 \right)^{5/2}} + \frac{1}{\left( 1/4 + u_0^2 \right)^{3/2}}\right] u_1^4 + \frac{3}{2} \left[ \frac{u_0}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{u_0^3}{\left( 1/4 + u_0^2 \right)^{5/2}} \right] u_1^2 u_2 + \frac{1}{2} \left[ \frac{u_0^2}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{1}{\left( 1/4 + u_0^2 \right)^{1/2}} \right] \left( u_2^2 + 2\, u_1 u_3 \right) - \frac{u_0 u_4}{\left( 1/4 + u_0^2 \right)^{1/2}} \\ &=& -\frac{6}{\sqrt{37}}\,u_4 - \frac{1}{37\sqrt{37}}\left[ u_2^2 + 2\,u_1 u_3 \right] + \frac{36}{1369\sqrt{37}}\, u_1^2 u_2 - \frac{143}{50653\sqrt{37}}\, u_1^4 , \\ A_5 &=& \frac{1}{120} \left[ - \frac{105\,u_0^5}{\left( 1/4 + u_0^2 \right)^{9/2}} + \frac{150\, u_0^3}{\left( 1/4 + u_0^2 \right)^{7/2}} - \frac{45\, u_0}{\left( 1/4 + u_0^2 \right)^{5/2}} \right] u_1^5 + \frac{1}{6} \left[ \frac{15\, u_0^4}{\left( 1/4 + u_0^2 \right)^{7/2}} - \frac{18\, u_0^2}{\left( 1/4 + u_0^2 \right)^{5/2}} + \frac{3}{\left( 1/4 + u_0^2 \right)^{3/2}}\right] u_1^3 u_2 + \frac{3}{2} \left[ \frac{u_0}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{u_0^3}{\left( 1/4 + u_0^2 \right)^{5/2}} \right] \left( u_1 u_2^2 + u_1^2 u_3 \right) + \left[ \frac{u_0^2}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{1}{\left( 1/4 + u_0^2 \right)^{1/2}} \right] \left( u_2 u_3 + u_1 u_4 \right) - \frac{u_0u_5}{\left( 1/4 + u_0^2 \right)^{1/2}} \\ &=& - \frac{6}{\sqrt{37}}\,u_5 - \frac{2}{37\sqrt{37}} \left[ u_2 u_3 + u_1 u_4 \right] + \frac{36}{1369\sqrt{37}} \left[ u_1 u_2^2 + u_1^2 u_3 \right] - \frac{572}{50653\sqrt{37}}\, u_1^3 u_2 + \frac{1692}{1874161\sqrt{37}} \, u_1^5 , \end{eqnarray*}

and so on. Then we get a sequence of initial value problems:

\[ \frac{\text d}{{\text d}x} \, u_{n+1} = A_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_{n+1} (0) =0, \quad n=0,1,2,\ldots . \]

Consequently, we have

\[ \frac{\text d}{{\text d}x} \, u_{1} = -\frac{\sqrt{37}}{2} \qquad \Longrightarrow \qquad u_1 (x) = -\frac{\sqrt{37}}{2}\, x , \]

\[ \frac{\text d}{{\text d}x} \, u_{2} = \frac{6}{\sqrt{37}} \,\frac{\sqrt{37}}{2}\,x \qquad \Longrightarrow \qquad u_2 (x) = \frac{3}{2}\, x^2 , \]

Now we proceed similarly.

\begin{eqnarray*} \frac{\text d}{{\text d}x} \, u_{3} &=& A_2 \qquad \Longrightarrow \qquad u_3 (x) = - \frac{\sqrt{37}}{12}\, x^3 , \\ \frac{\text d}{{\text d}x} \, u_{4} &=& A_3 \qquad \Longrightarrow \qquad u_4 (x) = \frac{x^4}{8} , \\ \frac{\text d}{{\text d}x} \, u_{5} &=& A_4 \qquad \Longrightarrow \qquad u_5 (x) = - \frac{\sqrt{37}}{240} \, x^5 , \end{eqnarray*}

Adding these functions, we obtain the 5-term Adomian approximation

\[ \phi_5 (x) = u_0 + u_1 + u_2 + u_3 + u_4 + u_5 = 3 - \frac{\sqrt{37}}{2}\,x + \frac{3\,x^2}{2} - \frac{\sqrt{37}}{12}\,x^3 + \frac{x^4}{8} - \frac{\sqrt{37}}{240}\,x^5 . \]

Finally we plot this approximation along with the true solution.

phi5[x_] = 3 - Sqrt[37]*x/2 + 3*x^2 /2 - Sqrt[37]*x^3 /12 + x^4 /8 - Sqrt[37]*x^5 /240 ;
Plot[{Callout[exact[x], "exact"], Callout[phi5[x], "ADM 5-approximation"]}, {x,0,3.0}, PlotStyle->Thick]

   ■

(a*b)/c+13*d
ur code
another line

(a*b)/c+13*d
ur code
another line