Adomian decomposition method

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. It is further extensible to stochastic systems by using the Ito integral. 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. 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 flexibly than direct Taylor series expansion.

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:

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

$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:

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

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

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

Replacing in the previous expression, we obtain:

$(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_{0}$ with some explicit expression on the right, and $y_{i},i=1,2,3,\cdots$, with some expression on the right containing terms of lower order than $i$. For instance:

${\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:

${\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:

${\frac {d^{3}u}{dx^{3}}}+{\frac {1}{2}}u{\frac {d^{2}u}{dx^{2}}}=0$

With the following conditions at the boundaries:

$u(0)=0$

$u^{\prime }(0)=0$

$u^{\prime }(x)\rightarrow 1,\;\;\;\;x\rightarrow \infty$

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

$Lu+Nu=0$

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

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

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

If:

$u=u^{0}+u^{1}+u^{2}+\cdots +u^{N}$

$u^{0}+u^{1}+u^{2}+\cdots +u^{N}=$

$\quad \alpha +\beta x+\gamma x^{2}/2-{\frac {1}{2}}L^{-1}(u^{0}+u^{1}+u^{2}+\cdots +u^{N}){\frac {d^{2}}{dx^{2}}}(u^{0}+u^{1}+u^{2}+\cdots +u^{N})$

and:

${\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:

$A_{n}={\frac {1}{n!}}{\frac {d^{n}}{d\lambda ^{n}}}f(u(\lambda ))\mid _{\lambda =0}$, where: ${\frac {d^{n}}{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:

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

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

$f(z)=\sum _{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 $\infty$ of this expression is $a_{L}/b_{M}$.

If we choose $b_{0}=1$, M linear equations for the b coefficients are obtained:

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

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

${\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:

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

Which when $\gamma =0.0408$ becomes:

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

With the limit:

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

Which is approximately equal to 1 (the 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. For instance, let us consider the following problem:

${\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:

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

$u(y=-y_{l})=g_{1}(x)$ and $u(y=y{l})=g_{2}(x)\qquad (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.

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

$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:

$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:

$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:

${\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}\;\;\;\;\;\;0<n<\infty \end{aligned}}$

where the couple ${a_{n}(y),b_{n}(y)}$ is the solution of the following system of equations:

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

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

$Nu=-b\partial _{x}u^{2}=-b\partial _{x}(u_{0}+u_{1}+u_{2}+u_{3}+...)(u_{0}+u_{1}+u_{2}+u_{3}+...)=$

$\quad -b\partial _{x}(u_{0}u_{0}+2u_{0}u_{1}+u_{1}u_{1}+2u_{0}u_{2}+...)=-b\partial _{x}\sum _{n=1}^{\infty }A(n-1)$,

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

Here $C(\nu ,n)$ are products (or sum of products) of &\nu; 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 $\sum _{n=0}^{\infty }A_{n}$ is equal to the sum of a generalized Taylor series about $u_{0}$ (Adomian, 1994).

For the example (1) the Adomian polynomials are:

${\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 (1989) 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 $\lim _{n\to \infty }\varphi ^{n}=u$.

With this method the solution can be found by systematically integrating along any of the two directions: in direction X we would use expression (3); in the alternative Y direction we would use the following expression:

$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}$:

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

If we call x-partial solution and y-partial solution to the two respective solutions so obtained, 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_{1},f_{2}\}$ or $\{g_{1},g_{2}\}$ 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_{1}$, $x=x_{2}$ must be consistent with those imposed at $y=y_{1}$ and $y=y_{2}$.

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

$u(x=0)=f_{1}(y)=0$

$u(x=x_{l})=f_{2}(y)=0$

By using Adomian's method and a symbolic processor (as Mathematica or Maple) it is easy to obtain the third order approximant to the solution. This approximant has an error lower than $5 \times 10^{-16}$ 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) (García-Olivares 2003-a).

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

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

and $g_{2}(x)=g_{1}(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 react with incredulity to these results. Is it not 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_{1},f_{2},g_{1},g_{2}$; however, given $f_{1},f_{2}$ it is always possible to find boundary functions $g_{1}^{*},g_{2}^{*}$ so close to $g_{1},g_{2}$ 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 tough visually not distinguishable. For instance, the solution with the boundary conditions:

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

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

$f_{1,2}(y)=0.00413683-0.00040048y-0.0813802y^{2}+0.0101279y^{3}+0.260417y^{4}-0.0694455y^{5}-0.277778y^{6}+0.15873y^{7}$

at x=0 and x=0.5, produce lateral functions with different sign convexity even tough 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 to obtain the boundary functions $g_{1}^{*},g_{2}^{*}$ that are compatible with any lateral set of conditions $\{f_{1},f_{2}\}$ 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_{1}$ (condition 1-a) with a N-order polynomial in y: $p_{1},p_{2}$ in such a way that: $f_{1}^{\prime }=f_{1}+p_{1},f_{2}^{\prime }=f_{2}+p_{2}$, where the norm of the two perturbation functions are smaller than the accuracy needed at the boundaries. These $p_{1},p_{2}$ depend on a set of polynomial coefficients $c_{i},i=1,\cdots ,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,\cdots ,N$. Finally, a boundary function $F(c_{1},c{2},\cdots ,c_{N})$ is defined as the sum of these four functions, and the distance between $F(c_{1},c{2},\cdots ,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_{1},c{2},\cdots ,c_{N})$ which has a global minimum for some combination of the parameters $c_{i},i=1,\cdots ,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).

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

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 to 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.