Preface
This section is devoted to a fascinating method for solving linear and nonlinear ordinary and partial differential equations invented by George Adomian (1922--1996). This method will be used in some other sections, so here we just give it a friendly introduction.
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340
Glossary
Adomian Decomposition Method
The Adomian decomposition method (ADM) is a systematic approximation method for solving nonlinear functional equations including ordinary and partial differential equations. It was named by Richard Bellman in honor of Adomian because it was developed from the 1970s --- 1990s by an American mathematician and aerospace engineer of Armenian descent George Adomian (1922-1996), chair of the Center for Applied Mathematics at the University of Georgia. The method is based on the assumption that the solution to a nonlinear functional equation can be represented by the infinite series \( y(x) = \sum_{n \ge 0} u_n (x) . \) The crucial aspect of the method is the decomposition of the problem under consideration into a sequence of recursively specified auxiliary problems that can be solved explicitly. This stage utilizes the employment of so called "Adomian polynomials" to represent the nonlinear portion of the equation as a convergent series with respect to these polynomials, without actual linearization of the system. These polynomials were named after G. Adomian by Randolph Rach in 1984 (A convenient computational form for the Adomian polynomials can be found in the Journal of Mathematical Analysis and Applications, Volume 102, Issue 2, September 1984, pages 415--419). There are several known algorithms for calculating Adomian polynomials; the following publications reflect efforts made in this area:
- Babolian, E. and Javadi, Sh., “New method for calculating Adomian polynomials”, Applied Mathematics and Computation, 2004, Volume 153, Issue 1, 25 May 2004, pages 253--259. https://doi.org/10.1016/S0096-3003(03)00629-5
- Biazar, J., Babolian, E., Kember, G., Nouri, A., Islam, R., An alternate algorithm for computing Adomian polynomials in special cases, Applied Mathematics and Computation, 2003, Volume 138, Issues 2–3, 20 June 2003, Pages 523--529; https://doi.org/10.1016/S0096-3003(02)00174-1
- Biazar, J. and Shafiof, S.M. A Simple Algorithm for Calculating Adomian Polynomials, International Journal of Contemporary Mathematical Sciences, 2007, Vol. 2, no. 20, pp. 975--982.
- Chang, S.H. and Chang, I.L., A new algorithm for calculating one-dimensional differential transform of nonlinear functions, Applied Mathematics and Computation, 2008, Vol. 195, Issue 2, pp. 799--808. https://doi.org/10.1016/j.amc.2007.05.026
- Chen, W., and Lu, Z., An algorithm for Adomian decomposition method, Applied Mathematics and Computation, 2004, Vol. 159, No. 1, (2004) pp. 221--235.
- Choi, H.-W., and Shin, J.G., Symbolic implementation of the algorithm for calculating Adomian polynomials, Applied Mathematics and Computation, 2003, Volume 146, Issue 1, 30 December 2003, Pages 257--271; https://doi.org/10.1016/S0096-3003(02)00541-6
- Daftardar-Gejji, V. and Jafari, H., An iterative method for solving non linear functional equations, Journal of Mathematical Analysis and Applications, 2006, vol. 316, Issue 2, pp. 753--763; https://doi.org/10.1016/j.jmaa.2005.05.009
- Jun-Sheng Duan, Recurrence triangle for Adomian polynomials, Applied Mathematics and Computation, 2010, Vol. 216, Issue 4, pp. 1235--1241. Elsevier. 15 April 2010; https://doi.org/10.1016/j.amc.2010.02.015
- Duan, J.-S., An efficient algorithm for the multivariable Adomian polynomials, Applied Mathematics and Computation, 2010, Vol. 217, Issue 6, pp. 2456--2467. https://doi.org/10.1016/j.amc.2010.07.046
- Jun-Sheng Duan, Convenient analytic recurrence algorithms for the Adomian polynomials, Applied Mathematics and Computation, 2011, Vol. 217, Issue 13, pp. 6337--6348. Elsevier. 1 March 2011. https://doi.org/10.1016/j.amc.2011.01.007
- Duan, J.-S., New recurrence algorithms for the nonclassic Adomian polynomials, Computers & Mathematics with Applications, 2011, Volume 62, Issue 8, pages 2961--2977. https://doi.org/10.1016/j.camwa.2011.07.074
- Elsaid, A., Adomian polynomials: a powerful tool for iterative methods of series solution of nonlinear equations, Journal of Applied Analysis and Computation, Volume 2, Number 4, November 2012, pp. 381--394.
- Fatoorehchi, H. and Abolghasemi, H. (2011b), “On calculation of Adomian polynomials by MATLAB”, Journal of Applied Computer Science and Mathematics (“Stefan cel Mare” University of Suceava), 2011, Vol. 11, No. 5, pp. 85--88.
- Jafari, H., Ghasempoor, S., and Khaliqu, C.M., A comparison between Adomian’s polynomials and He’s polynomials for nonlinear functional equations, Mathematical Problems in Engineering, Volume 2013, Article ID 943232, 4 pages.
- Kaliyappan, M. and Hariharan, S., Symbolic computation of Adomian polynomials based on Rach’s Rule, British Journal of Mathematics & Computer Science, Vol. 5, Issue 5, pp. 562--570, 2015, Article no.BJMCS.2015.
- Kataria, K.K. and Vellaisamy, Simple parametrization methods for generating Adomian polynomials, Applicable Analysis and Discrete mathematics, 2016, Vol. 10, N. 01, pp. 168--185; doi: 10.2298/AADM160123001K available online at http://pefmath.etf.rs
- Liao, S.J., On the homotopy analysis method for nonlinear problems, Applied Mathematics and Computation, 2004, Vol. 147, Issue 2, pp. 499--513; https://doi.org/10.1016/S0096-3003(02)00790-7
- Nigam, R. A New Formulation of Adomian Polynomials, International Journal of Mathematics and Scientific Computing, Vol. 5, No 2, 2015, pp. 92--97.
- Rach, R.C., “A new definition of the Adomian polynomials”, Kybernetes, Vol. 37 Issue: 7, 2008, pp.910--955, https://doi.org/10.1108/03684920810884342
- Wazwaz, A.M., A new algorithm for calculating Adomian polynomials for non-linear operators, Applied Mathematics and Computation, 111 (2000), pp. 53--69.
- Zaouagui I.N. and Badredine T., New Adomian’s Polynomials Formulas for the Non-linear and Nonautonomous Ordinary Differential Equations, Journal of Applied & Computational Mathematics, 2017, Vol 6(4), pp. doi: 10.4172/2168-9679.1000373
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.
R. Rach, Dr. George Adomian – distinguished scientist and mathematician, Kybernetes, Vol. 25, No. 9, 1996, 45-50.
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.
- 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.
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 a 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 the first finite number of components leads to an approximation to the true solution with any desired level of high precision.
If the given differential equation is homogeneous (without a driving term), then the ADM finite sum always provides the truncated series version to the true solution provided that the Adomian series converges. However, for inhomogeneous differential equations, the ADM finite sum approximation may include some noise terms that eventually are canceled out with subsequent iterations assuming that all iterations are complete. This phenomenon is well-known for Picard's iteration. 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 the 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 reviews on this topic can be found in the following articles: Ray's (New Approach for General Convergence of the Adomian Decomposition Method, World Applied Sciences Journal, 2014, Vol. 32, No. 11, pp. 2264--2268), Rach's (Bougoffa, L., Rach, R.C., El-Manouni, S., A convergence analysis of the Adomian decomposition method for an abstract Cauchy problem of a system of first-order nonlinear differential equations, International Journal of Computer Mathematics, 2013, Volume 90, Issue 2, pages 360--375), and Sunday's (Convergence analysis and implementation of Adomian decomposition method on second-order oscillatory problems, Asian Research Journal of Mathematics, 2017, Vol. 2, No. 5, pp. 1--12. Article no.ARJOM.32011) articles.
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 challenging 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 a demonstration of the Adomian decomposition method with the following initial value problem:
\[
y' = f(x,y) + g(x), \qquad y(x_0 ) = y_0 ,
\]
where f is the given (smooth) function, g is an input
(driving) term, y is the (unknown) output of the system, and
constants x0 and y0 are
prescribed. According to the ADM, we seek an approximation to the
solution of this initial value problem as an infinite series of
functions un to be determined recursively:
\[
y(x) = \sum_{n\ge 0} u_n = u_0 + u_1 + u_2 + \cdots ,
\]
while the nonlinear analytic function f(x,y) is assumed to be expanded into an Adomian series of polynomials
\[
f(x,y) = \sum_{n\ge 0} A_n (x) = A_0 (u_0 ) + A_1 (u_0 , u_1 ) + A_2 \left( u_0 , u_1 , u_2 \right) + \cdots ,
\]
where An are the (classical) Adomian’s polynomials defined recursively according to the formula:
\[
A_n \left( u_0 , u_1 , \ldots , u_n \right) = \frac{1}{n!} \, \frac{{\text d}^n}{{\text d} \lambda^n} \left[ f \left( x, \sum_{k\ge 0} u_k \lambda^k \right) \right]_{\lambda =0} , \qquad n=0,1,2,\ldots ;
\]
which is actually
the Faà
di Bruno's formula applied to the composition of two
functions: f and \( \sum_{k\ge 0} u_k . \)
These polynomials were first introduced by G. Adomian and R. Rach in their 1983 paper (Inversion of nonlinear stochastic operators,
Journal of Mathematical Analysis and Applications, Vol. 91, No. 1,
January 1983, pp. 39--46).
In the preceding formula, the parameter λ is a grouping parameter
that is used to identify a certain iteration term. The coefficient of
λn corresponds to the nth iteration term. One of the
ways to extract this term is to differentiate n times the
composition \( f \left( x, \sum_{k\ge 0} u_k
\lambda^k \right) \)
and set λ = 0 to eliminate the rest of the terms. Note that the previous
terms are eliminated by
differentiation. Fortunately, Mathematica has two dedicated
commands for extracting coefficients from a
polynomial: Coefficient that gives a particular coefficient
and CoefficientList that provides all coefficients of
the given polynomial (but not an arbitrary function).
The Adomian decomposition method assumes that the solution to the given initial value problem
\[
\frac{{\text d}y}{{\text d}x} = f(x,y), \qquad y(x_0 ) = y_0 ,
\]
admits the infinite series representation
\[
y(x) = u_0 + \sum_{n\ge 1} u_n .
\]
It is also assumed that the slope function can be expanded as an infinite series,
\[
f(x,y) = \sum_{n\ge 0} A_n \left( u_0 , u_1 , \ldots , u_n \right) ,
\]
where the An are Adomian's polynomials. To derive the expression for these polynomials, we consider the Taylor expansion of the slope function about u0:
\[
f(x,y) = f(u_0 ) + f' (u_0 ) \left( y - u_0 \right) + \frac{1}{2!}\, f''
(u_0 ) \left( y - u_0 \right)^2 + \frac{1}{3!}\, f^{(3)} (u_0 ) \left( y - u_0 \right)^3 + \cdots + \frac{1}{n!}\, f^{(n)} (u_0 ) \left( y - u_0 \right)^n + \cdots ,
\]
where we dropped writing the independent variable for simplicity because it does not affect the derivation. Substituting this expansion into the solution series and noting that
\[
y - u_0 = \left( u_0 + u_1 + u_2 + \cdots \right) - u_0 = u_1 + u_2 + u_3 + \cdots ,
\]
yields
\begin{align*}
f(y) &= f(u_0 ) + f' (u_0 ) \left( u_1 + u_2 + u_3 + \cdots \right)
\\
& \quad + \frac{1}{2!}\, f''
(u_0 ) \left( u_1 + u_2 + u_3 + \cdots \right)^2 + \frac{1}{3!}\, f''
(u_0 ) \left( u_1 + u_2 + u_3 + \cdots \right)^3 + \cdots .
\end{align*}
Expanding the first few terms of the Taylor expansion, we get
\begin{align*}
\left( u_1 + u_2 + u_3 + \cdots \right)^2 &= u_1^2 + 2u_1 u_2 + 2u_1 u_3 + u_2^2 + 2u_1 u_3 + \cdots ,
\\
\left( u_1 + u_2 + u_3 + \cdots \right)^3 &= u_1^3 + 3u_1^2 u_2 + 3u_1^2 u_3 + 3u_1 u_2^2 + 6u_1 u_2 u_3 + u_2^3 + \cdots .
\end{align*}
Substituting these expansions into the Taylor expansion of the slope function gives
\begin{align*}
f(y) &= f(u_0 ) + f' (u_0 ) \left( u_1 + u_2 + u_3 + \cdots \right)
\\
& \quad + \frac{1}{2!}\, f''
(u_0 ) \left( u_1^2 + 2u_1 u_2 + 2 u_1 u_3 + u_2^2 + 2 u_1 u_4 + 2 u_2 u_3 + \cdots \right)
\\
& \quad + \frac{1}{3!}\, f'''
(u_0 ) \left( u_1^3 + 3u_1^2 u_2 + 3 u_1^2 u_3 + 3 u_1 u_2^2 + 6u_1 u_2 u_3 + u_2^3 + \cdots \right) + \cdots .
\end{align*}
After rearranging the terms,
\begin{align*}
f(y) &= f(u_0 ) + f' (u_0 )\, u_1 + f' (u_0 ) \,u_2 + \frac{1}{2!}\, f''
(u_0 ) \, u_1^2 + f' (u_0 )\, u_3 + \frac{1}{2!}\, f''
(u_0 ) \, 2u_1 u_2
\\
& \quad + \frac{1}{3!}\, f'''
(u_0 ) \, u_1^3 + f' (u_0 )\, u_4 + \frac{1}{2!}\, f''
(u_0 ) \left( 2 u_1 u_3 + u_2^2 \right) + \frac{1}{4!}\, f^{(4)} (u_0 ) \, u_1^4 + \cdots .
\end{align*}
To calculate the first Adomian polynomial A0, consider any term where the sum of the subscripts of the ui add up 0. Next to get A1, consider each term where the subscripts of the ui add up 1. In a similar fashion, the general Adomian polynomial An includes any term where the sum of the subscripts of the ui add up to n.It is worth noting that the sum of the subscripts of the term uij is ij and not i.
▣ There are several known modifications of the Adomian polynomials, but we will use only one, known as accelerated Adomian polynomials (since there is no standard notation, we denote them by aAn) that were first introduced by George Adomian in his 1989 book (Nonlinear Stochastic Systems Theory and Applications to Physics, Kluwer Academic Publishers, Dordrecht):
\[
aA_n \left( u_0 , u_1 , \ldots , u_n \right) = f\left[
u_0 + u_1 + \cdots + u_n
\right] - f\left[ u_0 + u_1 + \cdots + u_{n-1}
\right] , \qquad J_0 = f\left[ u_0 \right] .
\]
In 2000, Khelifa and Cherruault established a bound for Adomian polynomials (Kybernetes, Vol. 29, No. 3, pp. 332--355):
\[
\left\vert A_n \right\vert \le \frac{(n+1)^n}{(n+1)!} \, M^{n+1} , \qquad \sup \left\vert f^{(k)} (x,y_0 ) \right\vert \le M .
\]
In general, since the Adomian method requires the existence of all derivatives (holomorphy) for the slope function, which is more restrictive than the Lipschitz condition required for Picard's procedure, it is expected that the ADM converges faster than its Picard counterpart.
The initial term u0 in the series expansion \( y(x) = \sum_{n \ge 0} u_n (x) \) is a solution of the truncated initial value problem:
\[
y' = g(x), \qquad y(x_0 ) = y_0 ,
\]
because it incorporates the influence of all inhomogeneous terms of the problem (that includes the initial condition and the forcing input function g). Once it is achieved, all other terms in the infinite series are solutions of the homogeneous initial value problems.
Therefore, u0 is usually not the same initial approximation as in Picard's iteration procedure that is just the initial value y0. It should be noted that when the Adomian method is applied to other problems (not necessarily initial value problems for first order differential equations as we consider in this section) one can choose another approximation to the nonhomgeneous differential equation: it is important to exclude, depending on x, the input term from the given slope function.
So in this section, we follow the standard ADM procedure applicable to differential equations of the form: \( y' = f(x,y) + g(x) . \)
Therefore, the initial term u0 could be chosen with some flexibility. Once u0 is determined, all other terms in Adomian's decomposition follow from the recursive initial value problems of Appell's type:
\[
\frac{{\text d}u_n}{{\text d}x} = A_{n-1} \left( u_0 , u_1 , \ldots , u_{n-1} \right) , \qquad u_n (x_0 ) = 0 , \qquad n=1,2,\ldots ,
\]
where the Adomian polynomials are determined recursively
\begin{eqnarray*}
A_0 &=& f(x, u_0 ) ,
\\
A_1 &=& u_1 f' (u_0 ) ,
\\
A_2 &=& u_2 f' (u_0 ) + \frac{1}{2!}\, u_1^2 f'' (u_0 ) ,
\\
A_3 &=& u_3 f' (u_0 ) + \frac{1}{3!}\, u_1^3 f''' (u_0 ) + \frac{2}{2!} \,u_1 u_2 f'' (u_0 ) ,
\\
A_4 &=& u_4 f' (u_0 ) + \left[ \frac{1}{2!}\, u_2^2 + u_1 u_3 \right] f'' (u_0 ) + \frac{1}{2!} \,u_1^2 u_2 f''' (u_0 ) + \frac{1}{4!} \, u_1^4 f'''' (u_0 ) ,
\\
A_5 &=& u_5 f' (u_0 ) + f'' (u_0 ) \left( u_1 u_4 + u_2 u_3 \right) + \frac{1}{2!}\,f''' (u_0 ) \left( u_1^2 u_3 + u_1 u_2^2 \right) + \frac{1}{3!}\, f^{(4)} (u_0 ) \, u_1^3 u_2 + \frac{1}{5!} \, f^{(5)} (u_0 ) \, u_1^5 .
\end{eqnarray*}
Determination of other Adomian polynomials requires more work.
\begin{eqnarray*}
A_6 &=& u_6 f' (u_0 ) + f'' (u_0 ) \left( u_1 u_5 + u_2 u_4 + u_3^2 \right) +
\frac{1}{3!}\,f''' (u_0 ) \left( 3\,u_1^2 u_4 + u_2^3 + 6\,u_1 u_2 u_3 \right)
\\
&\quad& + \frac{1}{4!}\,f^{(4)} (u_0 ) \left( 4\,u_1^3 + 6\,u_1^2 u_2^2 \right) + \frac{1}{4!}\,f^{(5)} (u_0 )\, u_1^4 u_2 + \frac{1}{6!} \,f^{(6)} (u_0 ) \, u_1^6 ,
\\
A_7 &=& u_7 f' (u_0 ) + f'' (u_0 ) \left( u_1 u_6 + u_2 u_5 + u_3 u_4 \right) +
\frac{1}{2}\,f''' (u_0 ) \left( u_1^2 u_5 + u_1 u_3^2 + u_3 u_2^2 + 2\,u_1 u_2 u_4 \right)
\\
& \quad& + \frac{1}{3!} \,f^{(4)} (u_0 ) \left( u_1^3 u_4 +3\, u_1^2 u_2 u_3 + u_1 u_2^3 \right) + \frac{1}{4!}\,f^{(5)} (u_0 ) \left( u_1^4 u_3 + 2\,u_1^3 u_2^2 \right) + \frac{1}{6!} \,f^{(6)} (u_0 ) \,u_1^5 u_2 + \frac{1}{7!} \,f^{(7)} (u_0 ) \,u_1^7 ,
\\
A_8 &=& u_8 f' (u_0 ) + f'' (u_0 ) \left( u_1 u_7 + u_2 u_6 + u_3 u_5 + u_4^2 /2 \right)
\\
& \quad& +
\frac{1}{2}\,f''' (u_0 ) \left( u_1^2 u_6 + u_2^2 u_4 + u_2 u_3^2 + 2\,u_1 u_2 u_5 + 2\, u_1 u_3 u_4 \right) + \frac{1}{3!} \,f^{(4)} (u_0 ) \left( u_1^3 u_5 +3\, u_1^2 u_2 u_4 + 3\, u_1 u_2^2 u_3 + 3\,u_1^2 u_3^2 /2 + u_2^4 /4 \right)
\\
&\quad& + \frac{1}{4!}\,f^{(5)} (u_0 ) \left( u_1^4 u_4 + 4\,u_1^3 u_2 u_3 + 2\, u_1^2 u_2^3 \right) + \frac{1}{6!} \,f^{(6)} (u_0 ) \left( u_1^5 u_3 + 15\, u_1^4 u_2^2 \right) + \frac{1}{6!} \,f^{(7)} (u_0 ) \,u_1^6 u_2 + \frac{1}{8!} \,f^{(8)} (u_0 ) \, u_1^8 ,
\\
\vdots &=& \vdots .
\end{eqnarray*}
Note: Adomian polynomials are not polynomials in the variable x. They are only polynomials in u1, … ,un, but not in u0. Additionally, components u1, … ,un might not be polynomials in the variable x.
Adomian's polynomials can be evaluated using the Duan's code (Applied Mathematical Modelling, 2013, 37, Issues 21-22, pp. 8687--8708) as follows:
Ado[M_]:= Module[{c,n,k,j,der},Table[c[n,k],{n,1,M},{k,1,n}];
A[0]=f[u[0]];
For[n = 1,n<=M,n++,c[n,1]=u[n];
For[k = 2,k<=n,k++,
c[n,k]=Expand[1/n*Sum[(j + 1)*u[j + 1]* c[n-1-j,k-1],{j,0,n-k}]]];
der = Table[D[f[u[0]],{u[0],k}],{k,1,M}];
A[n]=Take[der,n].Table[c[n,k],{k,1,n}]];
Table[A[n],{n,0,M}]]
A[0]=f[u[0]];
For[n = 1,n<=M,n++,c[n,1]=u[n];
For[k = 2,k<=n,k++,
c[n,k]=Expand[1/n*Sum[(j + 1)*u[j + 1]* c[n-1-j,k-1],{j,0,n-k}]]];
der = Table[D[f[u[0]],{u[0],k}],{k,1,M}];
A[n]=Take[der,n].Table[c[n,k],{k,1,n}]];
Table[A[n],{n,0,M}]]
Example: Let us consider the nonlinear (pendulum) operator: \( N\left[ \theta \right] = \sin\theta . \) Then the first few Adomian's polynomials for this nonlinear operator are
\begin{eqnarray*}
A_0 &=& \sin \theta_0 , \\
A_1 &=& \theta_1 \cos \theta_0 , \\
A_2 &=& - \frac{1}{2} \,\theta_1^2 \sin \theta_0 + \theta_2 \cos \theta_0 ,
\\
A_3 &=& - \frac{1}{6}\, \theta_1^3 \cos \theta_0 - \theta_1 \theta_2 \sin \theta_0 + \theta_3 \cos \theta_0 ,
\end{eqnarray*}
and so on. The modified polynomials are
\begin{eqnarray*}
J_1 &=& \sin \left( u_0 + u_1 \right) - \sin u_0 , \\
J_2 &=& \sin \left( u_0 + u_1 + u_2 \right) - \sin \left( u_0 + u_1 \right) ,
\end{eqnarray*}
and so on.
Let us consider the hyperbolic sine nonlinear function \( N\left[ x \right] = \sinh \left( x/2 \right) . \) The corresponding Adomian's polynomials are
\begin{eqnarray*}
A_1 &=& \frac{1}{2}\, x_1 \cosh \left( x_0 /2 \right) , \\
A_2 &=& \frac{1}{2}\, x_2 \cosh \left( \frac{x_0}{2} \right) + \frac{1}{2} \, x_1^2 \frac{1}{2^2} \,\sinh \left( \frac{x_0}{2} \right) , \\
A_3 &=& \frac{1}{3}\, x_3 \cosh \left( \frac{x_0}{2} \right) + \frac{1}{2^2} x_1 x_2 \sinh \left( \frac{x_0}{2} \right) + \frac{1}{6} \,x_1^3 \frac{1}{2^3} \,\cosh \left( \frac{x_0}{2} \right) \\
A_4 &=& \frac{1}{2}\, x_4 \cosh \left( \frac{x_0}{2} \right) + \left[ \frac{1}{2}\, x_2^2 + x_1 x_2 \right] \frac{1}{2^2} \sinh \left( \frac{x_0}{2} \right) + \frac{1}{2}\, x_1^2 x_2 \frac{1}{2^3} \,\cosh \left( \frac{x_0}{2} \right) + \frac{1}{24}\, \frac{1}{2^4} \,x_1^4 \sinh \left( \frac{x_0}{2} \right) ,
\end{eqnarray*}
and so on.
Our next example is the power function: \( N\left[ u \right] = u^{\gamma} , \) where γ is a real number. Then corresponding Adomian's polynomials are
\begin{eqnarray*}
A_0 &=& u_0^{\gamma} , \\
A_1 &=& \gamma \,u_0^{\gamma -1} u_1 , \\
A_2 &=& \gamma \,u_0^{\gamma -1} u_2 + \frac{1}{2}\,\gamma \left( \gamma -1 \right) u_0^{\gamma -2} u_1^2 , \\
A_3 &=& \gamma \,u_0^{\gamma -1} u_3 + \gamma \left( \gamma -1 \right) u_0^{\gamma -2} u_1 u_2 + \frac{1}{6}\, \gamma \left( \gamma -1 \right) \left( \gamma -2 \right) u_0^{\gamma -3} u_1^3 , \\
A_4 &=& \gamma \,u_0^{\gamma -1} u_4 + \gamma \left( \gamma -1 \right) u_0^{\gamma -2} \left[ \frac{1}{2}\, u_2^2 + u_1 u_3 \right] + \frac{1}{2} \, \gamma \left( \gamma -1 \right) \left( \gamma -2 \right) u_0^{\gamma -3} u_1^2 u_2 + \frac{1}{24} \, \gamma \left( \gamma -1 \right) \left( \gamma -2 \right) \left( \gamma -3 \right) u_0^{\gamma -4} u_1^4 .
\end{eqnarray*}
■ Example: Consider the simplest algebraic equation---the quadratic equation
\[
a\,x^2 + b\,x + c = 0, \qquad b\ne 0,
\]
where 𝑎, b, and c are some given constants. The case b = 0 can be easily included in our consideration by adding and subtracting the x term. The idea of isolating x is to transform the given problem to the form x = Φ(x), with some continuous function Φ, so that it appears as a fixed point problem. From the given quadratic equation, we obtain
\[
x = - \frac{c}{b} - \frac{a}{b}\, x^2 = \alpha + \gamma\,x^2 \qquad\mbox{with } \quad \alpha = - \frac{c}{b}, \ \gamma = - \frac{a}{b} .
\]
In general, one can add and subtract k x, k ≠ 0, and express x as
\[
x = -\frac{c}{k} + \frac{k-b}{k}\, x - \frac{a}{k}\, x^2 = \alpha + \beta\, x + \gamma \,x^2 \qquad (k \ne 0).
\]
In our case, we have k = b, but later we will consider other values of k.
It is helpful to work a numerical example and show all calculations for it. Hence, we consider the equation
\[
x^2 - x -6 =0,
\]
whose solutions are known to be x = 3 and x = -2. Let us see how ADM arrives at these values. We rewrite our quadratic equation in the form suitable for the fixed point theorem
\[
x = -6 + x^2 \qquad\Longrightarrow \qquad \alpha = -6, \ \beta = 0, \ \gamma = 1. \]
We decompose its solution into an infinite series
\[
x = x_0 + x_1 + x_2 + \cdots = \sum_{n\ge 0} x_n .
\]
The quadratic equation can rewritten in the operator form
\[
L \,x = N\left[ x \right] + g, \qquad L\,x = x, \quad N\left[ x \right] = \beta\,x + \gamma\,x^2 = x^2 , \quad g = \alpha = -6.
\]
Our next step is to decompose the nonlinear term x² as a series over Adomian's polynomials:
\[
N\left[ x \right] = \beta\,x + \gamma\,x^2 = \sum_{n\ge 0} A_n .
\]
The components xn of the decomposition solution
\( x= \sum_{n\ge 0} x_n \) are determined recursively
\begin{align*}
x_0 &= - \frac{c}{b} = -\frac{c}{k} = \alpha \qquad\Longrightarrow \qquad x_0 = -6,
\\
x_{n+1} &= A_n , \qquad n=0,1,2,\ldots .
\end{align*}
Using the explicit expression for N[x] = x², we find the first few Adomian polynomials:
\begin{align*}
A_0 &= x_0^2 = \frac{c^2}{b^2} \qquad \Longrightarrow \qquad A_0 = (-6)^2 = 36 = x_1,
\\
A_1 &= 2\,x_0 x_1 = -12\,x_1 \qquad \Longrightarrow \qquad A_1 = -12\,x_1 = -12 \cdot 36 = -432 = x_2 ,
\\
A_2 &= x^2_1 + 2\,x_0 x_2 \qquad \Longrightarrow \qquad A_2 = 3888 = x_3 ,
\\
A_3 &= 2\,x_1 x_2 + 2\,x_0 x_3 \qquad \Longrightarrow \qquad A_3 = -77760 = x_4 ,
\\
A_4 &= 2\,x_1 x_3 + 2\,x_0 x_4 + x_2^2 \qquad \Longrightarrow \qquad A_4 = 1399680 = x_5 ,
\\
A_5 &= 2\,x_1 x_4 + 2\,x_0 x_5 + 2\,x_2 x_3 \qquad \Longrightarrow \qquad A_5 = x_6 = ,
\end{align*}
and so on. In general for N[x] = βx + γx², we have x0 = α and
\begin{align*}
A_0 &= N(x_0 ) = N(\alpha ) = \beta\alpha + \gamma \alpha^2 = x_1,
\\
A_1 &= x_1 N' (x_0 ) = x_1 \left( \beta + 2\gamma\,\alpha \right) = \left( \beta\alpha + \gamma \alpha^2 \right) \left( \beta + 2\gamma\,\alpha \right) = x_2 ,
\\
A_2 &= x_2 N' (x_0 ) + x_1^2 N'' (x_0 )/2 = x_2 \left( \beta + 2\gamma\,\alpha \right) + \gamma x_1^2 = \left( \beta\alpha + \gamma \alpha^2 \right) \left( \beta + 2\gamma\,\alpha \right)^2 + \gamma \left( \beta\alpha + \gamma \alpha^2 \right)^2
\\
&= x_3 = \left( \beta\alpha + \gamma \alpha^2 \right) \left[ \left( \beta + 2\gamma\,\alpha \right)^2 + \gamma \left( \beta\alpha + \gamma \alpha^2 \right) \right] ,
\\
A_3 &= x_3 N' (\alpha ) + 2\gamma \,x_1 x_2 = \left( \beta\alpha + \gamma \alpha^2 \right) \left[ \left( \beta + 2\gamma\,\alpha \right)^2 + \gamma \left( \beta\alpha + \gamma \alpha^2 \right) \right] \left( \beta + 2\gamma\,\alpha \right) + 2 \gamma \left( \beta\alpha + \gamma \alpha^2 \right)^2 \left( \beta + 2\gamma\,\alpha \right)
\\
&= x_4 = \left( \beta\alpha + \gamma \alpha^2 \right) \left( \beta + 2\gamma \alpha \right) \left[ \left( \beta + 2\gamma \alpha \right) + 3\gamma \left( \beta\alpha + \gamma \alpha^2 \right) \right] ,
\\
A_4 &= x_4 N' (\alpha ) + N'' (\alpha ) \left[ x_2^2 + 2 x_1 x_3 \right] /2 = x_4 \left( \beta + 2\gamma \alpha \right) + \gamma \left[ x_2^2 + 2 x_1 x_3 \right]
\\
A_5 &= x_5 N' (\alpha ) + N'' (\alpha ) \left[ x_1 x_4 + x_2 x_3 \right]
\end{align*}
We immediately recognize that the series x = x0 + x1 + x2 + ··· is the Maclaurin series for the function
\[
\frac{-\beta - \beta \sqrt{1 - 4\alpha \gamma /\beta^2}}{2\gamma} =
\]
The series converges when \( \left\vert 4\alpha \gamma /\beta^2 \right\vert < 1 . \) So the series converges when
\( b^2 > |a\,c| . \) This condition is violated for our case and the Adomian series diverges. We can see this by finding derivatives of f(x) = x² - 5 at the roots:
\( f' (3) = 2\,(3) = 6 > 1 \) and \( f' (-2) = -4 , \) whose absolute value exceeds 1. Therefore, our fixed point equation x = f(x) is not a contractor (as required by the condition discovered by Cherruault for the ADM convergence).
To overcome this misfortune, we add and subtract 4x to the both sides of the given equation
\[
x^2 -x -6 + 4x = 4x \qquad \Longrightarrow \qquad x = f(x) = \frac{1}{4}\, x^2 + \frac{3}{4}\,x - \frac{3}{2} = N[x] - \frac{3}{2} .
\]
Now we choose the initial approximation to be x0 = -3/2, and expand the remaining function N[x] = x²/4 - 3x/4 into the Adomian series
\[
\frac{1}{4}\, x^2 + \frac{3}{4}\,x = \frac{x}{4} \left( x+3 \right) = \sum_{n\ge 0} A_n ,
\]
where N[x] = x(x + 3)/4. Substituting the series into the quadratic equation, we obtain
\[
x_{n+1} = A_n , \qquad n=0,1,2,\ldots .
\]
Now we calculate the Adomian polynomials using derivatives
\[
N(x) = \frac{x}{4} \left( x+3 \right) , \qquad N'(x) = \frac{2x+3}{4}, \qquad N''(x) = \frac{1}{2} .
\]
Therefore,
\begin{align*}
x_1 &= A_0 = N\left( -\frac{3}{2} \right) = - \frac{9}{16} = -0.5625,
\\
x_2 &= A_1 = x_1 N'\left( -\frac{3}{2} \right) = 0,
\\
x_3 &= A_2 = x_2 N'\left( -\frac{3}{2} \right) + x_1^2 N'' /2 = \frac{81}{1024} = 0.0791016 ,
\\
x_4 &= A_3 = x_3 N' \left( -\frac{3}{2} \right) + x_1 x_2 /2 = 0,
\\
x_5 &= A_4 = (x_2^2 + 2\,x_1 x_3)/4 = - \frac{729}{32 768} \approx -0.0222473 ,
\\
x_6 &= A_5 = 0,
\\
x_7 &= A_6 = \frac{19 683}{2 097 152} \approx 0.00938559 .
\end{align*}
f[x_] = x*(x + 3)/4
f1[x_] = (2*x + 3)/4
x0 = -3/2
x1 = f[x0]
x2 = x1*f1[x0]
x3 = x2*f1[x0] + x1^2/4
x4=x3*f1[x0]+x1*x2/2
x5=(x2^2 + 2*x1*x3)/4
x6=0
x7=(x1 x5 + x2 x4 + x3^2)/2
xx7=x0+x1+x3+x5+x7
f1[x_] = (2*x + 3)/4
x0 = -3/2
x1 = f[x0]
x2 = x1*f1[x0]
x3 = x2*f1[x0] + x1^2/4
x4=x3*f1[x0]+x1*x2/2
x5=(x2^2 + 2*x1*x3)/4
x6=0
x7=(x1 x5 + x2 x4 + x3^2)/2
xx7=x0+x1+x3+x5+x7
-(4186461/2097152)
Adding the first seven terms, we obtain the approximation of the root to be -1.99626. However, if we add one more term, we get a better approximation
-1.99978 to the actual root. Mathematica confirms
x9 = (x1 x7 + x3 x5)/2
x9 + xx7
x9 + xx7
-(33550737/16777216)
Now we find another root by choosing another modification:
\[
x^2 -x -6 - 4x = -4x \qquad \Longrightarrow \qquad x = f(x) = -\frac{1}{4}\, x^2 + \frac{5}{4}\,x + \frac{3}{2} = N[x] + \frac{3}{2} .
\]
Using N[x] = -x(x - 5)/4, we calculate the Adomian polynomials:
\begin{align*}
x_0 &= \frac{3}{2} , \\
x_1 &= A_0 = N\left( \frac{3}{2} \right) = \frac{21}{16} = 1.3125 ,
\\
x_2 &= A_1 = x_1 N'\left( \frac{3}{2} \right) = \frac{21}{32} \approx 0.65625 ,
\\
x_3 &= A_2 = x_2 N'\left( \frac{3}{2} \right) + x_1^2 N'' /2 = - \frac{105}{1024} \approx -0.102539 ,
\\
x_4 &= A_3 = x_3 N'\left( \frac{3}{2} \right) -x_1 x_2 /2 = - \frac{987}{2038} \approx -0.481934 .
\end{align*}
When we add the first seven components, we get 2.87075.
f[x_] = -x*(x - 5)/4
f1[x_] = (-2*x + 5)/4
x0 = 3/2
x1 = f[x0]
x2 = x1*f1[x0]
x3 = x2*f1[x0] - x1^2/4
x4 = x3*f1[x0] - x1*x2/2
x5 = x4*f1[x0] - (x2^2 + 2*x1*x3)/4
x6 = x5*f1[x0] - (x1*x4 + x2*x3)/2
x7 = x6*f1[x0] - (x1 x5 + x2 x4 + x3^2)/2
xx7 = x0 + x1 + x3 + x5 + x7
f1[x_] = (-2*x + 5)/4
x0 = 3/2
x1 = f[x0]
x2 = x1*f1[x0]
x3 = x2*f1[x0] - x1^2/4
x4 = x3*f1[x0] - x1*x2/2
x5 = x4*f1[x0] - (x2^2 + 2*x1*x3)/4
x6 = x5*f1[x0] - (x1*x4 + x2*x3)/2
x7 = x6*f1[x0] - (x1 x5 + x2 x4 + x3^2)/2
xx7 = x0 + x1 + x3 + x5 + x7
6020397/2097152
■ 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 this 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 the first few terms)
\[
y(x) = \sum_{n\ge 0} u_n (x) = u_0 (x) + u_1 (x) + u_2 (x) + \cdots ,
\]
where the first term u0 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 u0 differs from Picard's starting
term, which is 1. With u0, we can find A0, the zeroth Adomian polynomial:
\[
A_0 = f(x, u_0 ) = -2x \,u_0 = -2x - 3\,x^3 .
\]
f[x_, y_] = -2*x*y
f1[x_] = -2*x
u0[x_] = 1 + 3*x^2/2
A0[x_] = f[x, u0[x]]
f1[x_] = -2*x
u0[x_] = 1 + 3*x^2/2
A0[x_] = f[x, u0[x]]
-2 x (1 + (3 x^2)/2)
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_] := - 2*u*x
Ado[5]
Ado[5]
{-2 x 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 .
\]
This 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] /. %
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, \( \phi_2 = 1 + \frac{x^2}{2} - \frac{x^4}{4} , \) in the coefficient of the largest power of x. Next, we find u2 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*u[1][x], 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
u[2][x_] = y[x] /. %
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] /. %
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.
\]
According to Mathematica, its explicit solution is known to be (see Riccati section)
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}}]
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 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 a 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 u0 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_] = 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]
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 u1 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] /. %
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] /. %
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 φ2(x) with the true Maclaurin expansion, we see that φ2(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}} .
\]
The 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 A3:
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 the first four terms:
u4[x_] = 1/ 103524220800 (103524220800 x^4 + 27606458880 x^7 + 3524753232 x^10 + 273490560 x^13 + 13386165 x^16 + 382720 x^19)
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}]
error[x_] = y[x] - phi4[x]
Plot[{error[x]}, {x, 0, 0.7}, PlotStyle -> {Thick}]
 |
 |
|---|
| 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 |
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[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
\( \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} \)
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]
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 the 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
the power series approximation cannot detect the exact point no matter how
many terms we use.
Non standard ADM
\[
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 the 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] /. %
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
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] /. %
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} .
\]
The 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] /. %
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]]
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 a better approximation than the preceding algorithm.
Accelerated Adomian polynomials
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 formula above 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}}]
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 similar to the previous example:
\[
y(x) = \sum_{n\ge 0} u_n ,
\]
where the initial term u0 is exactly the same: \( u_0 = 1 + x^3 /3 . \)
Next, we find five Adomian polynomials using Mathematica
f[u_] := - u*u;
Ado[5]
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 the 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] /. %
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] /. %
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] /. %
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 ϕ2, 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]
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 ϕ3 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] /. %
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 accompanied by noise terms of powers greater than 4. The next term in the 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] /. %
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 the 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]
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]
| 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 the 5-term ADM approximation is less accurate than the Runge--Kutta algorithm of order four; Also, the fomer requires more labor 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
Series[exact[x], {x, 0, 12}]
\[
y(x) = \sum_{n\ge 0} u_n (x) ,
\]
where the initial approximation is taken to be u0 = 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]
Plot[{Callout[exact[x], "exact"], Callout[phi5[x], "ADM 5-approximation"]}, {x,0,3.0}, PlotStyle->Thick]
Some modifications of the ADM
\[
y' = f(x,y) + g(x) , \qquad y(x_0 ) = y_0
\]
as the sum of components \( y(x) = u_0 + \sum_{n\ge 1} u_n (x) , \) where u0 and un, n ≥ 1, are solutions of the following initial value problems
\[
u'_0 = g(x) , \qquad u_0 (x_0 ) = c_0 ,
\]
and
\[
u'_n = A_{n-1} \left( u_0 , u_1 , \ldots , u_{n-1} \right) , \qquad u_n (x_0 ) = c_n , \quad n=1,2,\ldots ,
\]
where { cn }n≥0 is an arbitrary sequence of numbers adding to y0:
\[
y_0 = c_0 + c_1 + c_2 + \cdots = \sum_{n\ge 0} c_n .
\]
The following examples demonstrate implementation of the latter.
Example: Consider the initial value problem
\[
y' + y^2 =0, \qquad y(0) = 2.
\]
Its solution is
\[
y(x) = \frac{2}{2x+1} , \qquad x > 1/2.
\]
DSolve[{y'[x] == -y[x]*y[x], y[0] == 2}, y[x], x] // Flatten
{y[x] -> 2/(1 + 2 x)}
We decompose the solution and the nonlinearity into series:
\[
y(x) = \sum_{n\ge 0} u_n (x) \qquad\mbox{and} \qquad y^2 = \sum_{n\ge 0} A_n \left( u_0 , u_1 , \ldots , u_n \right) ,
\]
where the Adomian polynomials are
\[
A_n = \sum_{k=0}^n u_k u_{n-k} , \qquad n=0,1,2,\ldots .
\]
Upon introducing a parameter c < 2, we decompose the initial condition as
\[
y(0) = 2 = \left( 2-c \right) + c .
\]
Then the first two components u0 and u1 in the Adomian series are solutions of the following initial value problems
\[
u'_0 = 0 , \qquad u_0 (0) = 2-c \qquad \Longrightarrow \qquad u_0 (x) = 2-c;
\]
and
\[
u'_1 = -u_0^2 , \qquad u_1 (0) = c \qquad \Longrightarrow \qquad u_1 (x) =c - \left( 2-c \right)^2 x .
\]
DSolve[{y'[x] == -(2-c)^2, y[0] == c}, y[x], x] // Flatten
{y[x] -> c - 4 x + 4 c x - c^2 x}
Other components are solutions of the IVPs with homogeneous initial conditions
\[
u'_{n+1} = -A_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_{n+1} (0) = 0 , \qquad n=1,2,\ldots .
\]
So
\begin{align*}
u_2 (x) &= \left( 2-c \right) \left[ -2cx + \left( 2-c \right)^2 x^2 \right] ,
\\
u_3 (x) &= -c^2 x + 3 \left( 2-c \right)^2 x^2 - \left( 2-c \right)^4 x^3
\\
u_4 (x) &= 6c^2 x^2 -3c^3 x^2 + \cdots .
\end{align*}
u0=2-c
u1[x_] = c - 4 x + 4 c x - c^2 x
DSolve[{y'[x] == -(2*u0*u1[x]), y[0] == 0}, y[x], x] // Flatten
u[2][x_] = y[x] /. %
DSolve[{y'[x] == -(2*u0*u[2][x] + u1[x]^2), y[0] == 0}, y[x], x] // Flatten
u[3][x_] = y[x] /. %
DSolve[{y'[x] == -(2*u0*u[3][x] + 2*u1[x]*u[2][x]), y[0] == 0}, y[x], x] // Flatten
u[4][x_] = y[x] /. %
DSolve[{y'[x] == -(2*u0*u[4][x] + 2*u1[x]*u[3][x] + u[2][x]^2), y[0] == 0}, y[x], x] // Flatten
u[5][x_] = y[x] /. %
phi5[x_,c_]= Simplify[u0 + u1[x]+u[2][x]+u[3][x]+u[4][x]+u[5][x]]
The fifth degree polynomial approximation becomes
u1[x_] = c - 4 x + 4 c x - c^2 x
DSolve[{y'[x] == -(2*u0*u1[x]), y[0] == 0}, y[x], x] // Flatten
u[2][x_] = y[x] /. %
DSolve[{y'[x] == -(2*u0*u[2][x] + u1[x]^2), y[0] == 0}, y[x], x] // Flatten
u[3][x_] = y[x] /. %
DSolve[{y'[x] == -(2*u0*u[3][x] + 2*u1[x]*u[2][x]), y[0] == 0}, y[x], x] // Flatten
u[4][x_] = y[x] /. %
DSolve[{y'[x] == -(2*u0*u[4][x] + 2*u1[x]*u[3][x] + u[2][x]^2), y[0] == 0}, y[x], x] // Flatten
u[5][x_] = y[x] /. %
phi5[x_,c_]= Simplify[u0 + u1[x]+u[2][x]+u[3][x]+u[4][x]+u[5][x]]
\[
\phi_5 (x,c) = u_0 (c) + u_1 (x,c) + u_2 (x,c) + u_3 (x,c) + u_4 (x,c) + u_5 (x,c) = \sum_{k=0}^5 u_k (x,c) .
\]
When c = 0, we get the classical Adomian series
\[
\phi (x) = y(x,0) = 2\sum_{n\ge 0} \left( -2x \right)^n ,
\]
which converges in the interval -1/2 < x < 1/2.
In order to see the effect of varying the value of c on the domain of convergence, we plot the curves of the explicit solution and fifth degree decomposition approximations for some values of c.
Plot[{2/(2*x + 1), phi5[x, 0], phi5[x, 1], phi5[x, 3/2]}, {x, -0.1,
1}, PlotStyle -> Thick,
PlotLabels -> {"true", "c=0", "c=1", "c=3/2"}]
Example: Consider the initial value problem with negative power nonlinearity
\[
\frac{{\text d}y}{{\text d}x} = \frac{2}{y} , \qquad y(0) = 3 .
\]
Its explicit solution is
\[
y(x) = \sqrt{9+4x} = 3 + \frac{2}{3}\, x - \frac{2}{3^3}\,x^2 + \frac{4}{243}\,x^3 - \frac{10}{2187}\, x^4 + \frac{28\,x^5}{19 683} - \frac{28\,x^6}{59 049} +\cdots .
\]
DSolve[{y'[x] == 2/y[x], y[0] == 3}, y[x], x] // Flatten
Series[Sqrt[9 + 4*x], {x, 0, 10}
The Maclaurin series above for the true solution converges when |x| < 9/4 = 2.25. The coefficients of this power series are expressed through the binomial coefficients
Series[Sqrt[9 + 4*x], {x, 0, 10}
\[
\sqrt{9+4x} = 3 \sum_{k\ge 0} \binom{1/2}{k} \left( \frac{4}{9}\,x \right)^k , \qquad \binom{1/2}{k} = \frac{\left( \frac{1}{2} \right)^{\underline{k}}}{k!} .
\]
By the ADM, we seek the solution to the given initial value problem in the form
\[
y(x) = u_0 (x) + \sum_{n\ge 1} u_n (x) ,
\]
where components un are determined recursively:
\begin{align*}
u'_0 (x) &= 0 , \qquad u_0 (0) = 3-c \qquad \Longrightarrow \qquad u_0 (x) = 3-c ,
\\
u'_1 (x) &= A_0 = \frac{2}{3-c} , \qquad u_1 (0) = c \qquad \Longrightarrow \qquad u_1 (x) =
\\
u'_{n+1} (x) &= A_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_{n+1} (0) =0, \qquad n=1,2,3,\ldots .
\end{align*}
u0[x_]= 3-c
DSolve[{y'[x] == 2/(3-c), y[0] == c}, y[x], x] // Flatten
u[1][x_] = (-3 c + c^2 - 2 x)/(-3 + c)
For determinantion of other components, we need the Adomian polynomials for the the nonlinearity 2/y:
\begin{align*}
A_1 &= -\frac{2\,u_1}{u_0^2} ,
\\
A_2 &= \frac{2\,u_1^2}{u_0^3} - \frac{2\,u_2}{u_0^2} ,
\\
A_3 &= \frac{4\,u_1 u_2}{u_0^3} - \frac{2\,u_1^3}{u_0^4} - \frac{2\,u_3}{u_0^2} ,
\\
A_4 &= \frac{2\,u_1^4}{u_0^5} - \frac{6\,u_1^2 u_2}{u_0^4} + \frac{2\,u_2^2 + 4\,u_1 u_3}{u_0^3} - \frac{2\,u_4}{u_0^2} ,
\\
A_5 &= \frac{8\,u_1^3 u_2}{u_0^5} - \frac{2\,u_1^5}{u_0^6} - \frac{6\,u_1 u_2^2 + 6\,u_1^2 u_3}{u_0^4} + \frac{4\,u_2 u_3 + 4\,u_1 u_4}{u_0^3} - \frac{2\,u_5}{u_0^2} .
\end{align*}
DSolve[{y'[x] == 2/(3-c), y[0] == c}, y[x], x] // Flatten
u[1][x_] = (-3 c + c^2 - 2 x)/(-3 + c)
f[y_] = 2/y
Ado[5]
We use Mathematica to find other components
Ado[5]
DSolve[{y'[x] == 2/(3 - c), y[0] == c}, y[x], x] // Flatten
u[1][x_] = y[x] /. %
u[1][x_] = y[x] /. %
(-3 c + c^2 - 2 x)/(-3 + c)
DSolve[{y'[x] == -2*u[1][x]/u0^2, y[0] == 0}, y[x], x] // Flatten
u[2][x_] = y[x] /. %
u[2][x_] = y[x] /. %
-((2 (-3 c x + c^2 x - x^2))/(-3 + c)^3)
DSolve[{y'[x] == 2*u[1][x]*u[1][x]/u0^3 - 2*u[2][x]/u0^2, y[0] == 0},
y[x], x] // Flattenn
u[3][x_] = y[x] /. %
u[3][x_] = y[x] /. %
-((2 (9 c^2 x - 6 c^3 x + c^4 x + 9 c x^2 - 3 c^2 x^2 + 2 x^3))/(-3 +
c)^5)
DSolve[{y'[x] ==
4*u[1][x]*u[2][x]/u0^3 - 2*u[1][x]^3 /u0^4 - 2*u[3][x]/u0^2,
y[0] == 0}, y[x], x] // Flatten
u[4][x_] = y[x] /. %
u[4][x_] = y[x] /. %
-((2 (-27 c^3 x + 27 c^4 x - 9 c^5 x + c^6 x - 54 c^2 x^2 +
36 c^3 x^2 - 6 c^4 x^2 - 30 c x^3 + 10 c^2 x^3 - 5 x^4))/(-3 +
c)^7)
DSolve[{y'[x] == 2*u[1][x]^4/u0^5 - 6*u[1][x]^2 *u[2][x] /u0^4 + 2*u[2][x]^2/u0^3 +4*u[1][x]*u[3][x]/u0^3 - 2*u[4][x]/u0^2,
y[0] == 0}, y[x], x] // Flatten
u[5][x_] = y[x] /. %
u[5][x_] = y[x] /. %
-((2 (-27 c^3 x + 27 c^4 x - 9 c^5 x + c^6 x - 54 c^2 x^2 +
36 c^3 x^2 - 6 c^4 x^2 - 30 c x^3 + 10 c^2 x^3 - 5 x^4))/(-3 +
c)^7)
We sum all these components to obtain fifth degree polynomial approximation
phi5[x_,c_]= u0 + u[1][x]+u[2][x]+ u[3][x]+u[4][x]+u[5][x]
Then we plot this function for different values of parameter c:
Plot[{Sqrt[9+4*x], phi5[x, 0], phi5[x, 1], phi5[x, 3/2]}, {x, -0,
1.5}, PlotStyle -> Thick,
PlotLabels -> {"true", "c=0", "c=1", "c=3/2"}]
- G. Adomian, K. Malakian, A comparison of the iterative method and Picard's successive approximations for deterministic and stochastic differential equations, Applied Mathematics and Computation, Vol. 8, No. 3, 187--204, 1981.
- N. Bellomo, D. Sarafyan, On Adomian's decomposition method and some comparisons with Picard's iterative scheme, Journal of Mathematical Analysis and Applications, Vol. 123, No. 2, 389--400, 1987.
- Duan, J.-S., Rach, R., Wang, Z., On the effective region of convergence of the decomposition series solution, Journal of Algorithms & Computational Technology, 2012, Vol. 7, No. 2, pp. 227--247.
- Ray, S.S., New Approach for General Convergence of the Adomian Decomposition Method, World Applied Sciences Journal, 2014, Vol. 32, No. 11, pp. 2264--2268; doi: 10.5829/idosi.wasj.2014.32.11.1317
Return to Mathematica page
Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)
Return to the Part 7 (Boundary Value Problems)