Preface
This section provides an illustration of application of the Adomian decomposition method (ADM for short) to second order singular differential equations.
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 V of the course APMA0330
Glossary
Lane--Emden equation
\[
\frac{1}{\xi^2}\,\frac{\text d}{{\text d}\xi} \left( \xi^2
\frac{{\text d} \theta}{{\text d}\xi} \right) + \theta^n = 0 \quad \mbox{or}
\quad\frac{{\text d}^2 \theta}{{\text d}\xi^2} + \frac{2}{\xi}\,
\frac{{\text d} \theta}{{\text d}\xi} + \theta^n = 0 .
.
\]
This equation is named after astrophysicists Jonathan Homer Lane (1819--1880) and Robert Emden (1862--1940). The Lane--Emden equation has been used to model several phenomena in mathematical physics and astrophysics such as thermal explosions, stellar structure, the thermal behavior of a spherical cloud of gas, isothermal gas spheres, and thermionic currents. Writing the above equation in the operator form
\[
L \left[ \xi , \texttt{D} \right] = \frac{1}{\xi^2}\,\frac{\text d}{{\text d}\xi} \left( \xi^2
\frac{{\text d} \theta}{{\text d}\xi} \right) , \qquad\mbox{where }
\texttt{D} = \frac{\text d}{{\text d}\xi} ,
\]
we can invert it:
\[
L^{-1} \left[ f(\xi ) \right] = \theta_0 + \left( \frac{1}{\xi_0} - \frac{1}{\xi} \right) \left( \xi^2 \frac{{\text d} \theta}{{\text d}\xi} \right)_{\xi = \xi_0} + \int_{\xi_0}^{\xi} \left( t - \frac{t^2}{\xi} \right) f(t) \,{\text d}t ,
\]
where ξ0 is the initial point and θ0 =
θ(ξ0) is the initial value.
Factorization of singular differential operators
\[
y'' + \frac{2n}{x}\, y' + \frac{n(n-1)}{x^2}\, y = g(x) + f(x,y).
\]
We will discuss more general equations later. Let us consider
the initial conditions
\[
y\left( x_0 \right) = a, \qquad y' \left( x_0 \right) = b ,
\]
where x0 is an arbitrary real number other than zero because
the origin is a singular point for this differential equation. Without any loss of generality, we can assume that x0 is positive. The driving function g(x) is assumed to be known and f(x, y)
is the nonlinearity term. In preparation to apply the Adomian decomposition method, it is convenient to introduce the differential operator
\[
L\left[ x, \texttt{D} \right] = \texttt{D}^2 + \frac{2n}{x}\, \texttt{D} +
\frac{n(n-1)}{x^2}\,\texttt{I} = x^{-n} \texttt{D}^2 \left[ x^n \left( \cdot \right) \right] ,
\]
where \( \texttt{D} = {\text d}/{\text d}x \) is the
derivative operator and \( \texttt{I} \) is the identity operator. Since the inverse operator is known to be
\[
L^{-1} \left[ u(x) \right] = \left[ (1-n) \left( \frac{x_0}{x} \right)^n + n
\left( \frac{x_0}{x} \right)^{n-1} \right] y(x_0 ) +
\left( \frac{x_0}{x} \right)^n \left( x - x_0 \right) y' \left( x_0 \right)
+ x^{-n} \int_{x_0}^x \left( x -t \right) t^n u(t) \,{\text d}t ,
\]
we are ready to apply the ADM to solve the given initial value problem.
Suppose that we want to incorporate a constant term to our differential
operator. This task is accomplished with the differential operator:
(D[(Sin[w*x])^2 *D[x/Sin[w*x]*f[x], x], x])/x/Sin[w*x]
Simplify[%]
Simplify[%]
w^2 f[x] + (2 Derivative[1][f][x])/x + (f^\[Prime]\[Prime])[x]
\[
L \left[ y(x) \right] = \frac{1}{x\,\sin \omega x}\,
\frac{\text d}{{\text d}x}\, \sin^2 \omega x \,
\frac{\text d}{{\text d}x} \left( \frac{x}{\sin \omega x} \,y(x) \right) =
\frac{{\text d}^2 y}{{\text d}x^2} + \frac{2}{x}\,
\frac{{\text d} y}{{\text d}x} - \omega^2 y
\]
Its inverse is
\[
u(x) = L^{-1} \left[ f(x) \right] = \frac{x_0 \sin \omega x}{x\,\sin \omega x_0} \, u(x_0 ) + \frac{\sin^2 \omega x_0}{\omega} \left( \cot \omega x - \cot \omega x_0 \right) \left( \frac{\text d}{{\text d}x} \left( \frac{x}{\sin \omega x} \,u(x) \right) \right)_{x=x_0} +
\frac{\sin \omega x}{\omega\,x} \int_{x_0}^x \left( \cot \omega x - \cot \omega t \right) t\,\sin \omega t\, f(t)\,{\text d}t .
\]
Next, we consider the differential operator
\[
L \left[ x, \texttt{D} \right] = x^{-m-n} \, \texttt{D}\, x^m \texttt{D}\, x^n =\texttt{D}^2 + \frac{m+ 2n}{x}\, \texttt{D} + \frac{n(m+n-1)}{x^2}\,\texttt{I} .
\]
Simplify[D[x^m *D[x^n *f[x], x], x]/x^(m + n)]
(n (-1 + m + n) f[x] +
x ((m + 2 n) Derivative[1][f][x] + x (f^\[Prime]\[Prime])[x]))/x^2
This allows us to determine its inverse:
\[
L^{-1} \left[ u(x) \right] = u \left( x_0 \right) \left( n\,\frac{x_0^{m+n-1}}{x^{n}} \,\frac{x^{1-m} - x_0^{1-m}}{1-m} + \frac{x_0^n}{x^n} \right) +
u' \left( x_0 \right) \frac{x_0^{m+n}}{x^n} \, \frac{x^{1-m} - x_0^{1-m}}{1-m} + \frac{1}{1-m} \, x^{-n} \int_{x_0}^x t^{m+n}\, f(t) \left( x^{1-m} - t^{1-m} \right) {\text d}t .
\]
Example: Consider the second-order nonlinear differential equation with the product of the derivative and a sinusoidal nonlinearity of the solution
\[
\frac{{\text d}^2 y}{{\text d}x^2} - \frac{2}{x}\,
\frac{{\text d}y}{{\text d}x} + \pi^2 y(x) =0, \qquad y(1) = 0, \quad
y' (1) =-\pi .
\]
Applying the inverse operator formular to zero function, we obtain
\[
y(x) = \frac{\sin \pi x}{x} .
\]
■ Example: Consider the initial value problem
\[
y'' +
\]
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 Dufing
equation.
■
Factorization of Euler's operators
\[
\frac{{\text d}^2 y}{{\text d} x^2} + \frac{p}{x}\,\frac{{\text d}y}{{\text d}x} + \frac{q}{x^2}\, y(x) = g(x) ,
\]
which we rewrite in the operator form
\[
L \left[ x, \texttt{D} \right] y = \left[ \texttt{D}^2 + \frac{p}{x}\, \texttt{D} + \frac{q}{x^2}\, \texttt{I}
\right] y(x) = g(x) ,
\]
where p, q are some real numbers, g(x) is a given
input function,
\( \texttt{D} = {\text d}/{\text d}x \) is the
derivative operator and \( \texttt{I} \) is the
identity operator. One of the ways to determine a solution f the inhomogeneous
differential equation, is its factorization. If the following system of
algebraic equations
\[
\begin{split}
p &= m +2n , \\
q &= n \left( m+n-1 \right) ,
\end{split}
\]
has two distinct real solutions (i.e., \( (p-1)^2 - 4q > 0 \) ), then the operator L admits factorization:
\[
L \left[ x, \texttt{D} \right] = x^{-m-n} \texttt{D} \left[ x^m \texttt{D}
\,x^n \left( \cdot \right) \right] .
\]
Example: Consider the following Euler's operators: ■
Another type of factorization is possible for the differential operator:
\[
L \left[ x, \texttt{D} \right] = \left( \texttt{D} + \frac{n}{x}\,
\texttt{I} \right)^2 = \texttt{D}^2 + \frac{2n}{x}\,\texttt{D} +
\frac{n(n-1)}{x^2} \,\texttt{I}
\]
which was considereed previously. The Euler operator admits such factorization,
\[
\texttt{D}^2 + \frac{p}{x}\,\texttt{D} + \frac{q}{x^2}\,\texttt{I} =
\left( \texttt{D} + \frac{n}{x}\,
\texttt{I} \right)^2
\]
if and only if there exist a real nuber n such that
\[
\begin{split}
p &= 2n , \\
q &= n \left( n-1 \right) .
\end{split}
\]
This may happen only when
\[
4q = p^2 -2p .
\]
Adomian Decomposition Metod for Second Order ODEs
The Adomian decomposition method and its modifications can be considered as an 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 permits us to solve both nonlinear initial value problems and boundary value problems without unphysical 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.
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 term).
The linear operator is usually included in the nonlinear one N[y] so we will drop it. The Adomian method 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.
Very few nonlinear problems have simple, closed-form solutions. In most cases, solutions of nonlinear differential equations can be expressed only by means of an infinite series. Application of the Adomian decomposition method is based on the assumption that the unknown solution is represented as an infinite sum:
\[
y(x) = \sum_{n\ge 0} u_n (x) .
\]
The crusual 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 unknown function and nonlinear input through generating functions
\begin{eqnarray*}
y(x, \lambda ) &=& \sum_{n\ge 0} \lambda^n u_n (x) ,
\\
N\left( y(x, \lambda ) , \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:
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}]]
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}] ]
S the required solution is represented by the series
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}] ]
\[
y(x) = u_0 (x) + \sum_{n\ge 1} u_n (x) ,
\]
where the initial term is dtermined by
\[
u_0 (x) = L^{-1} \left[ g(x) \right] ,
\]
and all other terms correspond to the initial initial conditions and homogeneous equation. So we apply
\[
L^{-1} \left[ u(x) \right] =
\left( \frac{x_0}{x} \right)^n \left( x - x_0 \right) y' \left( x_0 \right)
+ x^{-n} \int_{x_0}^x \left( x -t \right) t^n u(t) \,{\text d}t ,
\]
with y(x0) = 0 and
y'(x0) = 0; to the recurrence
\[
y_{n+1} = L^{-1} \left[ A_n (x) \right] \qquad y_{n+1} = L^{-1} \left[ J_n (x) \right]
\]
depending whether the classical Adomian decomposition will be used or its
generalization.
Example: Consider the initial value problem
\[
y'' + \frac{2}{x}\, y' = y^2 - x^4 +6 , \qquad y(1) = 1, \quad y' (1) =2.
\]
The exact solution is y(x) = x².
We rewrite the given differential equation in operator form:
\[
y'' + \frac{2}{x}\, y' = x^{-1} \frac{{\text d}^2}{{\text d} x^2} \left( x\, y(x) \right) .
\]
The nonlinear operator and the driven term are
\[
N\left[ y \right] = y^2 \qquad\mbox{and} \qquad g(x) = 6 - x^4 ,
\]
respectively. The required solution is decomposed as
\[
y(x) = u_0 (x) + \sum_{n\ge 1} u_n (x) ,
\]
where the first term becomes
\[
u_0 (x) = L^{-1} \left[ 6 - x^4 \right] =
\frac{1-n + n\,x}{x^{n}} + 2\,\frac{x-1}{x^n}
+ x^{-n} \int_{1}^x \left( x -t \right) t^n (6-t^4 ) \,{\text d}t =
x^{-n} \left[ \frac{6}{1+n} -1-n - \frac{1}{6+n} +2x + nx - \frac{6x}{1+n} + \frac{x}{5+n} + \frac{6}{2+3n+n^2}\,x^{2+n} - \frac{x^{6=n}}{(5+n)(6+n)} \right] .
\]
Simplify[1 - n + n*x + 2*x - 2 +
Integrate[(x - t)*t^n *(6 - t^4), {t, 1, x}]]
-1 + 6/(2 + n) - 1/(6 + n) +
n (-1 + x) + (2 - 6/(1 + n) + 1/(5 + n)) x + (6 x^(2 + n))/(
2 + 3 n + n^2) - x^(6 + n)/((5 + n) (6 + n))
■ Example: Consider the initial value problem
\[
y''' + \frac{2}{x}\, y'' = y^2 - x^4 +6 , \qquad y(1) = 1, \quad y' (1) =2.
\]
The exact solution is y(x) = x².
■ Example:
Example: Move to chapter 4, application4 Consider the Lane--Emden type equation with an exponential nonlinearity of the solution with respect to the variable base
\[
\frac{{\text d}^2 y}{{\text d}x^2} + \frac{2}{x}\, \frac{{\text d}y}{{\text d}x} + (x+1)^{y(x)}, \qquad y(0) =1, \quad y' (0) =0.
\]
■ Example: The system of partial differential equations, called KdV equations, cna be reduced by appropriate substitution to the following ordinary differential equation
\[
v' = k\,v''' - v\,v' .
\]
■
ADM for singular differential equations
- 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.
- 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.
- Hasan, Y.Q., A new development to the Adomian decomposition for solving singularIVPs of Lane-Emden Type, United States of America Research Journal (USARJ), 2014, Vol. 2, No.3, 2014, ISSN 2332-2160
- 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)