This section provides an illustration of application of the Adomian decomposition method (ADM for short) to second order differential equations. This method is one of the most effective methods for constructing analytic approximate solutions of nonlinear differential equations.

Return to computing page for the first course APMA0330
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
Return to Part IV of the course APMA0330

Adomian Decomposition Method

The Adomian decomposition method (ADM for short) and its modifications can be considered as an quantitative rather than merely qualitative analytic approximation method, which does not require accepting a priori assumptions in our cybernetics models that drastically alter the outcomes so that they do not faithfully replicate reality. It decomposes the solution into a series of components to be determined by recursion, and then decomposes the analytic nonlinearity into the series of the Adomian polynomials. The ADM is applicable in both nonlinear initial value problems and boundary value problems without nonphysical restrictive assumptions such as required by linearization, perturbation, ad hoc assumptions, guessing the initial term or a set of basis functions, and so forth, most of which would change the physical behavior of the problem. It will be very naive to assume that the Adomian decomposition scheme can solve any nonlinear problem. Nevertheless, the ADM is the first candidate for a trial.

Nonlinear differential equation models have been extensively applied in the majority of applications. However, very few nonlinear problems have simple, closed-form solutions. In most cases, solutions of nonlinear differential equations can be expressed only by means of special functions or infinite series. Before we start exposition of the Adomian decomposition method, it makes sense to look at some simple examples that show how solutions to the differential equations are sensitive to numerical values.

Consider an ordinary differential equation written in operator form:

\[ L\left[ y \right] = L_0 \left[ y \right] + N\left[ y \right] + g(x) , \]
where \( L\left[ y \right] = \texttt{D}^2 + a\,\texttt{D} + b\, \texttt{I} \) with \( \texttt{D} = {\text d}/{\text d}x \) being the derivative operator and \( \texttt{I} \) being the identity operator, is the linear differential operator acting on a function y(x), \( L_0 \left[ y \right] = p(x)\,\texttt{D} + q(x)\, \texttt{I}\) is a linear differential operator of order 1, N[y] represents an analytic nonlinear operator, and g(x) is a specified analytic input function (usually called the driving or forcing term).

As a rule, the linear operator L0 is included in the nonlinear one N[y] so we will drop it. The Adomian procedure does not react whether lower order terms as \( a\,\texttt{D} + b\, \texttt{I} \) are included in L or not. Including these terms depends whether you can obtain the explicit formula for the inverse L-1 or not.

The Adomian decomposition method assumes that the unknown solution can be represent by an infinite sum:

\[ y(x) = \sum_{n\ge 0} u_n (x) . \]
The crucial part of the method is decomposition of the nonlinear term by an infinite series of the Adomian polynomials:

\[ N \left[ y(x) \right] = \sum_{n\ge 0} A_n \left( u_0 , u_1 , \ldots , u_n \right) . \]
Representing the unknown function and nonlinear input through generating functions
\begin{eqnarray*} Y(x, \lambda ) &=& \sum_{n\ge 0} \lambda^n u_n (x) , \\ N\left( x, Y(x, \lambda ) \right) &=& \sum_{n\ge 0} \lambda^n A_n , \end{eqnarray*}
we find the Adomian's polynomials with the aid of the Faà di Bruno's formula:
\[ A_n (u) = \left. \frac{1}{n!}\,\frac{\partial^n}{\partial \lambda^n} \, N \left[ y(x,\lambda ) \right] \right\vert_{\lambda =0} , \qquad n=0,1,2,\ldots , \]
where λ is a grouping parameter. In computational practice, we truncate the decomposition series after n = M for some finite M because the higher order solution components un for n > M do not contribute to the accuracy of calculations.

Mathematica code for evaluating Adomian polynomials for the one-variable Adomian polynomials:

AP1[f_, M_]:=
Module[{F}, Subscript[F, 0] = f[Sum[Subscript[u, k] * s^k, {k, 0, M}] ];
For [i = 0, i <'= M, i++ , A[i] = Collect[Expand[1/i! * (Subscript[ F, i]/. s -> 0)], Derivative[_][f ][_]];
Subscript[F, i + 1] = D[Subscript[F, i], s]];
Table[A[i], {i, 0, M}]]

Another Mathematica code for evaluating Adomian polynomials for the one-variable Adomian polynomials:
AP2[f_, M_]:=
Module[{c, n, k, j, der}, Table[c[n, k], {n, 1, M}, {k, 1, n}];
der = Table [ D[f [Subscript[u, 0]], {Subscript[u, 0], k}], {k, 1, M}];
A[0] = f [Subscript[u, 0]];
For[n = 1, n <= M, n++ , c[n, 1] = Subscript[u, n];
For[k = 2, k <= n, k++, c[n, k] = Expand[1/n * Sum[(j + 1) * Subscript[u, j + 1] * c[n-1-j, k-1],{j, 0, n-k}]]];
A[n] = Take[der, n].Table[c[n, k], {k, 1, n}]];
Table[A[n], {n, 0, M}] ]

Once we know how to determine Adonian's polynomials, we proceed according to the main streamline. First, we find the initial term by solving nonhomogeneous initial value problem for the linear differential equation:

\[ L \left[ u_0 \right] = g(x), \qquad u_0 \left( x_0\right) = y\left( x_0 \right) = y_0 , \quad u'_0 \left( x_0 \right) = y'\left( x_0 \right) = v_0 . \]
When \( L = \texttt{D}^2 , \) its inverse is immediate
\[ u_0 (x) = L^{-1} \left[ g(x) \right] = y_0 + \left( x- x_0 \right) v_0 + \int_{x_0}^x \left( x- t \right) g(t) \,{\text d}t . \]
When L contains lower derivative terms, finding its inverse may require some nontrivial efforts. With u0 at hand, all other terms are obtained recursively by solving the same type of initial value problems:
\[ L \left[ u_{n+1} \right] = A_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_{n+1} \left( x_0\right) = 0, \quad u'_{n+1} \left( x_0 \right) = 0 , \qquad n=0,1,2,\ldots . \]
Now it is clear why it is crucial to have an explicit expression for L-1. Next, the Adomian iteration scheme goes smoothly, unless you cannot determine the Adomian polynomial at some stage, or you are satisfied with the obtained accuracy.

We demonstrate the Adomian decomposition method in a series of examples.

Example: Consider the initial value problem for undamped Duffing equation:

\[ y'' + k\, y + \varepsilon \,y^3 =0, \qquad y(0)=1 , \quad y' (0) =0. \]
We can rewrite the Duffing equation in the operator form:
\[ L\left[ y \right] \equiv \texttt{D}^2 y = N \left[ y \right] \equiv \varepsilon \,y^3 - k\, y. \]
Physically, the resilience of the oscillator is directly proportional to N[y]. When k ≥ 0 and ε < 0, there is one unique equilibrium point y = 0 where the nonlinear term vanishes. However, when k ≤ 0 and ε > 0, there are two equilibrium points: y = 0 and \( y = \sqrt{|k/\varepsilon |} . \) From the physical points of view, the sum of the kinetic and potential energy of the oscillator keeps the same, therefore it is clear that the oscillation motion is periodic, no matter k is positive or negative. Thus, from physical points of view, it is easy to know that y(t) is periodic, even if we do not directly solve the Duffing equation.    ■

Example: Consider the second-order nonlinear differential equation with the exponential input

\[ \frac{{\text d}^2 y}{{\text d}x^2} - \frac{2+x}{1+x} \, \frac{{\text d}y}{{\text d}x} + y^2 (x) = x^2 e^{2x}, \qquad y(0) =0, \quad y' (0) =1. \]
This IVP has the exact solution \( y = \phi (x) = x\,e^x . \)    ■

Example: Consider the second-order nonlinear differential equation with the exponential input

\[ y'' + x\,y' -2\,e^y = 2x\,\tan x, \qquad y(0) = y' (0) =0. \]
This IVP has the explicit solution \( y = \phi (x) = -2\,\ln | \cos x | , \) which is defined on the interval -π/2 < x < π/2.

We input in Mathematica notebook:

Init={{0,0}}; Alpha={{{0,x}}}; Be={−2}; F={E∧(v1)}; G={2x∗Tan[x]};
and then run
PSSOL[Init, Alpha, Be, F, G, 18]
to output the 21-term approximant of the solution,φ21(x). The Mathematica command PadeApproximant[φ21(x), {x, 0,{10, 10}}] outputs the diagonal Padé approximant[10/10]{φ21(x)}.    ■

Example: The system of partial differential equations, called KdV equations, can be reduced by appropriate substitution to the following ordinary differential equation

\[ v' = k\,v''' - v\,v' . \]


  1. Duan, J.-S., Rach, R., Baleanu, D., and Wazwaz, A.-M., “A review of the Adomian decomposition method and its applications to fractional differential equations”, Communications in Fractional Calculus, 2012, Vol. 3, No. 2, pp. 73--99.
  2. J.-S. Duan, J.-S., Rach, R., and Wazwaz, A.-M., A new modified Adomian decomposition method for higher-order nonlinear dynamical systems, CMES: Computer Modeling in Engineering & Sciences, 2013, Vol. 94, No. 1, pp. 77--118.
  3. Elsaid, A., Adomian polynomials: a powerful tool for iterative methods of series solution of nonlinear equations, Journal of Applied Analysis and Computation, 2012, Vol. 2, No 4, pp. 381--394.
  4. Rach, R., Wazwaz, A.-M., and Duan, J.-S., A reliable modification of the Adomian decomposition method for higher-order nonlinear differential equations, Kybernetes, 2013, Vol. 42, No 2, pp. 282--308, doi: 10.1108/03684921311310611


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)