The Adomian decomposition method (ADM) is a semi-analytical method for solving ordinary and partial nonlinear differential equations. The method was developed from the 1970s to the 1990s by George Adomian, chair of the Center for Applied Mathematics at the University of Georgia.^[1] It is further extensible to stochastic systems by using the Ito integral.^[2] The aim of this method is towards a unified theory for the solution of partial differential equations (PDE); an aim which has been superseded by the more general theory of the homotopy analysis method.^[3] The crucial aspect of the method is employment of the "Adomian polynomials" which allow for solution convergence of the nonlinear portion of the equation, without simply linearizing the system. These polynomials mathematically generalize to a Maclaurin series about an arbitrary external parameter; which gives the solution method more flexibility than direct Taylor series expansion.^[4]

Ordinary differential equations

Adomian method is well suited to solve Cauchy problems, an important class of problems which include initial conditions problems.

Application to a first order nonlinear system

An example of initial condition problem for an ordinary differential equation is the following:

${\displaystyle y^{\prime }(t)+y^{2}(t)=-1,}$

${\displaystyle y(0)=0.}$

To solve the problem, the highest degree differential operator (written here as L) is put on the left side, in the following way:

${\displaystyle Ly=-1-y^{2},}$

with L = d/dt and ${\displaystyle L^{-1}=\int _{0}^{t}()}$. Now the solution is assumed to be an infinite series of contributions:

${\displaystyle y=y_{0}+y_{1}+y_{2}+y_{3}+\cdots .}$

Replacing in the previous expression, we obtain:

${\displaystyle (y_{0}+y_{1}+y_{2}+y_{3}+\cdots )=y(0)+L^{-1}[-1-(y_{0}+y_{1}+y_{2}+y_{3}+\cdots )^{2}].}$

Now we identify y₀ with some explicit expression on the right, and y_i, i = 1, 2, 3, ..., with some expression on the right containing terms of lower order than i. For instance:

${\displaystyle {\begin{aligned}&y_{0}&=&\ y(0)+L^{-1}(-1)&=&-t\\&y_{1}&=&-L^{-1}(y_{0}^{2})=-L^{-1}(t^{2})&=&-t^{3}/3\\&y_{2}&=&-L^{-1}(2y_{0}y_{1})&=&-2t^{5}/15\\&y_{3}&=&-L^{-1}(y_{1}^{2}+2y_{0}y_{2})&=&-17t^{7}/315.\end{aligned))}$

In this way, any contribution can be explicitly calculated at any order. If we settle for the four first terms, the approximant is the following:

${\displaystyle {\begin{aligned}y&=y_{0}+y_{1}+y_{2}+y_{3}+\cdots \\&=-\left[t+{\frac {1}{3))t^{3}+{\frac {2}{15))t^{5}+{\frac {17}{315))t^{7}+\cdots \right]\end{aligned))}$

Application to Blasius equation

A second example, with more complex boundary conditions is the Blasius equation for a flow in a boundary layer:

${\displaystyle {\frac {\mathrm {d} ^{3}u}{\mathrm {d} x^{3))}+{\frac {1}{2))u{\frac {\mathrm {d} ^{2}u}{\mathrm {d} x^{2))}=0}$

With the following conditions at the boundaries:

${\displaystyle {\begin{aligned}u(0)&=0\\u^{\prime }(0)&=0\\u^{\prime }(x)&\to 1,\qquad x\to \infty \end{aligned))}$

Linear and non-linear operators are now called ${\displaystyle L={\frac {\mathrm {d} ^{3)){\mathrm {d} x^{3))))$ and ${\displaystyle N={\frac {1}{2))u{\frac {\mathrm {d} ^{2)){\mathrm {d} x^{2))))$, respectively. Then, the expression becomes:

${\displaystyle Lu+Nu=0}$

and the solution may be expressed, in this case, in the following simple way:

${\displaystyle u=\alpha +\beta x+\gamma x^{2}/2-L^{-1}Nu}$

where: ${\displaystyle L^{-1}\xi (x)=\int dx\int \mathrm {d} x\int \mathrm {d} x\;\;\xi (x)}$ If:

${\displaystyle {\begin{aligned}u&=u^{0}+u^{1}+u^{2}+\cdots +u^{N}\\&=\alpha +\beta x+\gamma x^{2}/2-{\frac {1}{2))L^{-1}(u^{0}+u^{1}+u^{2}+\cdots +u^{N}){\frac {\mathrm {d} ^{2)){\mathrm {d} x^{2))}(u^{0}+u^{1}+u^{2}+\cdots +u^{N})\end{aligned))}$

and:

${\displaystyle {\begin{aligned}u^{0}&={}\alpha +\beta x+\gamma x^{2}/2\\u^{1}&=-{\frac {1}{2))L^{-1}(u^{0}u^{0''})&=&-L^{-1}A_{0}\\u^{2}&=-{\frac {1}{2))L^{-1}(u^{1}u^{0''}+u^{0}u^{1''})&=&-L^{-1}A_{1}\\u^{3}&=-{\frac {1}{2))L^{-1}(u^{2}u^{0''}+u^{1}u^{1''}+u^{0}u^{2''})&=&-L^{-1}A_{2}\\&\cdots \end{aligned))}$

Adomian’s polynomials to linearize the non-linear term can be obtained systematically by using the following rule:

${\displaystyle A_{n}={\frac {1}{n!)){\frac {\mathrm {d} ^{n)){\mathrm {d} \lambda ^{n))}f(u(\lambda ))\mid _{\lambda =0},}$

where: ${\displaystyle {\frac {\mathrm {d} ^{n)){\mathrm {d} \lambda ^{n))}u(\lambda )\mid _{\lambda =0}=n!u_{n))$

Boundary conditions must be applied, in general, at the end of each approximation. In this case, the integration constants must be grouped into three final independent constants. However, in our example, the three constants appear grouped from the beginning in the form shown in the formal solution above. After applying the two first boundary conditions we obtain the so-called Blasius series:

${\displaystyle u={\frac {\gamma }{2))x^{2}-{\frac {\gamma ^{2)){2))\left({\frac {x^{5)){5!))\right)+{\frac {11\gamma ^{3)){4))\left({\frac {x^{8)){8!))\right)-{\frac {375\gamma ^{4)){8))\left({\frac {x^{11)){11!))\right)+\cdots }$

To obtain γ we have to apply boundary conditions at ∞, which may be done by writing the series as a Padé approximant:

${\displaystyle f(z)=\sum _{n=0}^{L+M}c_{n}z^{n}={\frac {a_{0}+a_{1}z+\cdots +a_{L}z^{L)){b_{0}+b_{1}z+\cdots +b_{M}z^{M))))$

where L = M. The limit at ${\displaystyle \infty }$ of this expression is a_L/b_M.

If we choose b₀ = 1, M linear equations for the b coefficients are obtained:

${\displaystyle \left[{\begin{array}{cccc}c_{L-M+1}&c_{L-M+2}&\cdots &c_{L}\\c_{L-M+2}&c_{L-M+3}&\cdots &c_{L+1}\\\vdots &\vdots &&\vdots \\c_{L}&c_{L+1}&\cdots &c_{L+M-1}\end{array))\right]\left[{\begin{array}{c}b_{M}\\b_{M-1}\\\vdots \\b_{1}\end{array))\right]=-\left[{\begin{array}{c}c_{L+1}\\c_{L+2}\\\vdots \\c_{L+M}\end{array))\right]}$

Then, we obtain the a coefficients by means of the following sequence:

${\displaystyle {\begin{aligned}a_{0}&=c_{0}\\a_{1}&=c_{1}+b_{1}c_{0}\\a_{2}&=c_{2}+b_{1}c_{1}+b_{2}c_{0}\\&\cdots \\a_{L}&=c_{L}+\sum _{i=1}^{\min(L,m)}b_{i}c_{L-i}.\end{aligned))}$

In our example:

${\displaystyle u'(x)=\gamma x-{\frac {\gamma ^{2)){2))\left({\frac {x^{4)){4!))\right)+{\frac {11\gamma ^{3)){4))\left({\frac {x^{7)){7!))\right)-{\frac {375\gamma ^{4)){8))\left({\frac {x^{10)){10!))\right)}$

Which when γ = 0.0408 becomes:

${\displaystyle u'(x)={\frac {0.0204+0.0379\,z-0.0059\,z^{2}-0.00004575\,z^{3}+6.357\cdot 10^{-6}z^{4}-1.291\cdot 10^{-6}z^{5)){1-0.1429\,z-0.0000232\,z^{2}+0.0008375\,z^{3}-0.0001558\,z^{4}-1.2849\cdot 10^{-6}z^{5))},}$

with the limit:

${\displaystyle \lim _{x\to \infty }u'(x)=1.004.}$

Which is approximately equal to 1 (from boundary condition (3)) with an accuracy of 4/1000.

Partial differential equations

Application to a rectangular system with nonlinearity

One of the most frequent problems in physical sciences is to obtain the solution of a (linear or nonlinear) partial differential equation which satisfies a set of functional values on a rectangular boundary. An example is the following problem:

${\displaystyle {\frac {\partial ^{2}u}{\partial x^{2))}+{\frac {\partial ^{2}u}{\partial y^{2))}-b{\frac {\partial u^{2)){\partial x))=\rho (x,y)\qquad (1)}$

with the following boundary conditions defined on a rectangle:

${\displaystyle u(x=0)=f_{1}(y)\quad {\text{and))\quad u(x=x_{l})=f_{2}(y)\qquad {\text{(1-a)))}$

${\displaystyle u(y=-y_{l})=g_{1}(x)\quad {\text{and))\quad u(y=y_{l})=g_{2}(x)\qquad {\text{(1-b)))}$

This kind of partial differential equation appears frequently coupled with others in science and engineering. For instance, in the incompressible fluid flow problem, the Navier–Stokes equations must be solved in parallel with a Poisson equation for the pressure.

Decomposition of the system

Let us use the following notation for the problem (1):

${\displaystyle L_{x}u+L_{y}u+Nu=\rho (x,y)\qquad (2)}$

where L_x, L_y are double derivate operators and N is a non-linear operator.

The formal solution of (2) is:

${\displaystyle u=a(y)+b(y)x+L_{x}^{-1}\rho (x,y)-L_{x}^{-1}L_{y}u-L_{x}^{-1}Nu\qquad (3)}$

Expanding now u as a set of contributions to the solution we have:

${\displaystyle u=u_{0}+u_{1}+u_{2}+u_{3}+\cdots }$

By substitution in (3) and making a one-to-one correspondence between the contributions on the left side and the terms on the right side we obtain the following iterative scheme:

${\displaystyle {\begin{aligned}u_{0}&=a_{0}(y)+b_{0}(y)x+L_{x}^{-1}\rho (x,y)\\u_{1}&=a_{1}(y)+b_{1}(y)x-L_{x}^{-1}L_{y}u_{0}+b\int dxA_{0}\\&\cdots \\u_{n}&=a_{n}(y)+b_{n}(y)x-L_{x}^{-1}L_{y}u_{n-1}+b\int dxA_{n-1}\quad 0<n<\infty \end{aligned))}$

where the couple {a_n(y), b_n(y)} is the solution of the following system of equations:

${\displaystyle {\begin{aligned}\varphi ^{n}(x=0)&=f_{1}(y)\\\varphi ^{n}(x=x_{l})&=f_{2}(y),\end{aligned))}$

here ${\displaystyle \varphi ^{n}\equiv \sum _{i=0}^{n}u_{i))$ is the nth-order approximant to the solution and N u has been consistently expanded in Adomian polynomials:

${\displaystyle {\begin{aligned}Nu&=-b\partial _{x}u^{2}=-b\partial _{x}(u_{0}+u_{1}+u_{2}+u_{3}+\cdots )(u_{0}+u_{1}+u_{2}+u_{3}+\cdots )\\&=-b\partial _{x}(u_{0}u_{0}+2u_{0}u_{1}+u_{1}u_{1}+2u_{0}u_{2}+\cdots )\\&=-b\partial _{x}\sum _{n=1}^{\infty }A(n-1),\end{aligned))}$

where ${\displaystyle A_{n}=\sum _{\nu =1}^{n}C(\nu ,n)f^{(\nu )}(u_{0})}$ and f(u) = u² in the example (1).

Here C(ν, n) are products (or sum of products) of ν components of u whose subscripts sum up to n, divided by the factorial of the number of repeated subscripts. It is only a thumb-rule to order systematically the decomposition to be sure that all the combinations appearing are utilized sooner or later.

The ${\displaystyle \sum _{n=0}^{\infty }A_{n))$ is equal to the sum of a generalized Taylor series about u₀.^[1]

For the example (1) the Adomian polynomials are:

${\displaystyle {\begin{aligned}A_{0}&=u_{0}^{2}\\A_{1}&=2u_{0}u_{1}\\A_{2}&=u_{1}^{2}+2u_{0}u_{2}\\A_{3}&=2u_{1}u_{2}+2u_{0}u_{3}\\&\cdots \end{aligned))}$

Other possible choices are also possible for the expression of A_n.

Series solutions

Cherruault established that the series terms obtained by Adomian's method approach zero as 1/(mn)! if m is the order of the highest linear differential operator and that ${\displaystyle \lim _{n\to \infty }\varphi ^{n}=u}$.^[5] With this method the solution can be found by systematically integrating along any of the two directions: in the x-direction we would use expression (3); in the alternative y-direction we would use the following expression:

${\displaystyle u=c(x)+d(x)y+L_{y}^{-1}\rho (x,y)-L_{y}^{-1}L_{x}u-L_{y}^{-1}Nu}$

where: c(x), d(x) is obtained from the boundary conditions at y = - y_l and y = y_l:

${\displaystyle {\begin{aligned}u(y=-y_{l})&=g_{1}(x)\\u(y=y_{l})&=g_{2}(x)\end{aligned))}$

If we call the two respective solutions x-partial solution and y-partial solution, one of the most interesting consequences of the method is that the x-partial solution uses only the two boundary conditions (1-a) and the y-partial solution uses only the conditions (1-b).

Thus, one of the two sets of boundary functions {f₁, f₂} or {g₁, g₂} is redundant, and this implies that a partial differential equation with boundary conditions on a rectangle cannot have arbitrary boundary conditions on the borders, since the conditions at x = x₁, x = x₂ must be consistent with those imposed at y = y₁ and y = y₂.

An example to clarify this point is the solution of the Poisson problem with the following boundary conditions:

${\displaystyle {\begin{aligned}u(x=0)&=f_{1}(y)=0\\u(x=x_{l})&=f_{2}(y)=0\end{aligned))}$

By using Adomian's method and a symbolic processor (such as Mathematica or Maple) it is easy to obtain the third order approximant to the solution. This approximant has an error lower than 5×10⁻¹⁶ in any point, as it can be proved by substitution in the initial problem and by displaying the absolute value of the residual obtained as a function of (x, y).^[6]

The solution at y = -0.25 and y = 0.25 is given by specific functions that in this case are:

${\displaystyle g_{1}(x)=0.0520833\,x-0.347222\,x^{3}+9.25186\times 10^{-17}x^{4}+0.833333\,x^{5}-0.555556\,x^{6))$

and g₂(x) = g₁(x) respectively.

If a (double) integration is now performed in the y-direction using these two boundary functions the same solution will be obtained, which satisfy u(x=0, y) = 0 and u(x=0.5, y) = 0 and cannot satisfy any other condition on these borders.

Some people are surprised by these results; it seems strange that not all initial-boundary conditions must be explicitly used to solve a differential system. However, it is a well established fact that any elliptic equation has one and only one solution for any functional conditions in the four sides of a rectangle provided there is no discontinuity on the edges. The cause of the misconception is that scientists and engineers normally think in a boundary condition in terms of weak convergence in a Hilbert space (the distance to the boundary function is small enough to practical purposes). In contrast, Cauchy problems impose a point-to-point convergence to a given boundary function and to all its derivatives (and this is a quite strong condition!). For the first ones, a function satisfies a boundary condition when the area (or another functional distance) between it and the true function imposed in the boundary is so small as desired; for the second ones, however, the function must tend to the true function imposed in any and every point of the interval.

The commented Poisson problem does not have a solution for any functional boundary conditions f₁, f₂, g₁, g₂; however, given f₁, f₂ it is always possible to find boundary functions g₁^*, g₂^* so close to g₁, g₂ as desired (in the weak convergence meaning) for which the problem has solution. This property makes it possible to solve Poisson's and many other problems with arbitrary boundary conditions but never for analytic functions exactly specified on the boundaries. The reader can convince himself (herself) of the high sensitivity of PDE solutions to small changes in the boundary conditions by solving this problem integrating along the x-direction, with boundary functions slightly different even though visually not distinguishable. For instance, the solution with the boundary conditions:

${\displaystyle f_{1,2}(y)=0.00413682-0.0813801\,y^{2}+0.260416\,y^{4}-0.277778\,y^{6))$

at x = 0 and x = 0.5, and the solution with the boundary conditions:

${\displaystyle {\begin{aligned}f_{1,2}(y)=0.00413683&-0.00040048\,y-0.0813802\,y^{2}+0.0101279\,y^{3}+0.260417\,y^{4}\\&-0.0694455\,y^{5}-0.277778\,y^{6}+0.15873\,y^{7}+\cdots \end{aligned))}$

at x = 0 and x = 0.5, produce lateral functions with different sign convexity even though both functions are visually not distinguishable.

Solutions of elliptic problems and other partial differential equations are highly sensitive to small changes in the boundary function imposed when only two sides are used. And this sensitivity is not easily compatible with models that are supposed to represent real systems, which are described by means of measurements containing experimental errors and are normally expressed as initial-boundary value problems in a Hilbert space.

Improvements to the decomposition method

At least three methods have been reported ^{[6] [7] [8]} to obtain the boundary functions g₁^*, g₂^* that are compatible with any lateral set of conditions {f₁, f₂} imposed. This makes it possible to find the analytical solution of any PDE boundary problem on a closed rectangle with the required accuracy, so allowing to solve a wide range of problems that the standard Adomian's method was not able to address.

The first one perturbs the two boundary functions imposed at x = 0 and x = x₁ (condition 1-a) with a Nth-order polynomial in y: p₁, p₂ in such a way that: f₁' = f₁ + p₁, f₂' = f₂ + p₂, where the norm of the two perturbation functions are smaller than the accuracy needed at the boundaries. These p₁, p₂ depend on a set of polynomial coefficients c_i, i = 1, ..., N. Then, the Adomian method is applied and functions are obtained at the four boundaries which depend on the set of c_i, i = 1, ..., N. Finally, a boundary function F(c₁, c₂, ..., c_N) is defined as the sum of these four functions, and the distance between F(c₁, c₂, ..., c_N) and the real boundary functions ((1-a) and (1-b)) is minimized. The problem has been reduced, in this way, to the global minimization of the function F(c₁, c₂, ..., c_N) which has a global minimum for some combination of the parameters c_i, i = 1, ..., N. This minimum may be found by means of a genetic algorithm or by using some other optimization method, as the one proposed by Cherruault (1999).^[9]

A second method to obtain analytic approximants of initial-boundary problems is to combine Adomian decomposition with spectral methods.^[7]

Finally, the third method proposed by García-Olivares is based on imposing analytic solutions at the four boundaries, but modifying the original differential operator in such a way that it is different from the original one only in a narrow region close to the boundaries, and it forces the solution to satisfy exactly analytic conditions at the four boundaries.^[8]

Integral Equations

The Adomian decomposition method may also be applied to linear and nonlinear integral equations to obtain solutions.^[10] This corresponds to the fact that many differential equation can be converted into integral equations.^[10]

Adomian Decomposition Method

The Adomian decomposition method for nonhomogenous Fredholm integral equation of the second kind goes as follows:^[10]

Given an integral equation of the form:

${\displaystyle u(x)=f(x)+\lambda \int _{a}^{b}K(x,t)u(t)dt}$

We assume we may express the solution in series form:

${\displaystyle u(x)=\sum _{n=0}^{\infty }u_{n}(x)}$

Plugging the series form into the integral equation then yields:

${\displaystyle \sum _{n=0}^{\infty }u_{n}(x)=f(x)+\lambda \int _{a}^{b}K(x,t)(\sum _{n=0}^{\infty }u_{n}(t))dt}$

Assuming that the sum converges absolutely to ${\displaystyle u(x)}$ we may integerchange the sum and integral as follows:

${\displaystyle \sum _{n=0}^{\infty }u_{n}(x)=f(x)+\lambda \int _{a}^{b}\sum _{n=0}^{\infty }K(x,t)u_{n}(t)dt}$

${\displaystyle \sum _{n=0}^{\infty }u_{n}(x)=f(x)+\lambda \sum _{n=0}^{\infty }\int _{a}^{b}K(x,t)u_{n}(t)dt}$

Expanding the sum on both sides yields:

${\displaystyle u_{0}(x)+u_{1}(x)+u_{2}(x)+...=f(x)+\lambda \int _{a}^{b}K(x,t)u_{0}(t)dt+\lambda \int _{a}^{b}K(x,t)u_{1}(t)dt+\lambda \int _{a}^{b}K(x,t)u_{2}(t)dt+...}$

Hence we may associate each ${\displaystyle u_{i}(x)}$ in the following recurrent manner:

${\displaystyle u_{0}(x)=f(x)}$

${\displaystyle u_{i}(x)=\lambda \int _{a}^{b}K(x,t)u_{i-1}dt,\,\,\,\,\,\,\,\,i\geq 1}$

which gives us the solution ${\displaystyle u(x)}$ in the solution form above.

Example

Given the Fredholm integral equation:

${\displaystyle u(x)=\cos(x)+2x+\int _{0}^{\pi }xt\cdot u(t)dt}$

Since ${\displaystyle f(x)=\cos(x)+2x}$, we can set:

${\displaystyle u_{0}(x)=\cos(x)+2x}$

${\displaystyle u_{1}(x)=\int _{0}^{\pi }xt\cdot u_{0}dt=\int _{0}^{\pi }xt\cdot (cosx+2x)dt=(-2+{\frac {2\pi ^{3)){3)))x}$

${\displaystyle u_{2}(x)=\int _{0}^{\pi }xt\cdot u_{1}dt=\int _{0}^{\pi }xt\cdot (-2+{\frac {2\pi ^{3)){3)))t\,dt=({\frac {-2\pi ^{3)){3))+{\frac {2\pi ^{6)){9)))x}$

...

Hence the solution ${\displaystyle u(x)}$ may be written as:

${\displaystyle u(x)=\cos(x)+2x+(-2+{\frac {2\pi ^{3)){3)))x+({\frac {-2\pi ^{3)){3))+{\frac {2\pi ^{6)){9)))x+...}$

Since this is a telescoping series, we can see that every terms after ${\displaystyle cos(x)}$ cancels and may be regarded as "noise",^[10] Thus, ${\displaystyle u(x)}$ becomes:

${\displaystyle u(x)=\cos(x)}$

Gallery