Preface


This section discusses applications of the Adomian decomposition method (ADM for short) to first order nonlinear differential equations with singularities (application of the ADM for regular nonlinear differential equations was done previously in section). The ADM calculates the solutions of nonlinear equations as infinite series in which each term can be determined in principle. The crucial part of the ADM is the decomposition of the nonliearity as the series with respect to generalized polynomials called Adomian’s polynomials. The series-solution is convergent towards an accurate solution, and in some cases (at least for homogeneous equations) it gives the exact solution. To surpress noise terms that appear in Adomian's solutions for inhomogeneous equations, Abdul-Majid Wazwaz (A reliable modification of Adomian decomposition method, Applied Mathematics and Computation, 1999, Vol. 102 No. 1, pp. 77-86.) proposed an ADM modification based to division of the input into two parts that are incorporated in the first two terms of ADM series solution. Since the calculations of the Adomian’s polynomials are based on the Faà di Bruno's formula, the number of terms in this formula grows exponentially; so this task cannot be accomplished on an electrical computer, generally speaking. Therefore, a truncated version of the ADM is usually used for approximations. There are known four kinds of Adomian polynomials (see: Duan, J.-S., Rach, R., Băleanu, 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.) and that certain circumstances could benefit from a different choice; however for most problems the classical Adomian polynomials are satisfactory. A special class of them constitute the accelerated Adomian polynomials that were first introduced by George Adomian in his 1989 book, and were detailed and intensively used in 2008 Randolph Rach article.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
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

ADM for first order ODEs


As a typical example, we consider the general singular differential equation
\[ \frac{{\text d} y}{{\text d}x} + \frac{p}{x}\, y(x) = g(x) + f(x,y) , \qquad y(x_0 ) = y_0 , \]
where p and x0 are some positive numbers, while y0 is an arbitrary constant, but it is kept fixed for the given initial value problem. Here g(x) is a forcing (or driving) function (input) and f(x, y) is a nonlinear operator. It is convenient to rewrite the given singular differential equation in operator form. Here we have two options to do this: either to isolate the derivative operator \( \texttt{D} = {\text d}/{\text d}x \) or to combine it together with the linear term. So we set two linear differential operators:
\[ L\left[ y \right] = \frac{{\text d} y}{{\text d}x} + \frac{p}{x}\, y = \left( \texttt{D} + \frac{p}{x}\,\texttt{I} \right) y\qquad\mbox{and} \qquad L_0 \left[ y \right] = \frac{{\text d} y}{{\text d}x} = \texttt{D}\, y , \]
where \( \texttt{I} \) is the identity operator. Correspondingly, we represent the given differential equation in two forms:
\[ L\left[ y \right] = g(x) + f(x,y) \qquad\mbox{or} \qquad L_0 \left[ y \right] = g(x) + f_0 (x,y) , \]
where \( f_0 (x,y) = f(x,y) - p\,y(x)/x . \) The inverse to the unbounded differential operator L was found previously:
\[ L^{-1} \left[ u \right] = y(t_0 )\,x^{-p}x_0^p + x^{-p} \int_{x_0}^x t^p u(t)\,{\text d}t . \]
Since \( L_0 = \texttt{D} \) is just the familiar derivative operator, its (left) inverse is known to be
\[ L^{-1}_0 \left[ u \right] = y(t_0 ) + \int_{x_0}^x u(t)\,{\text d}t . \]
These two inverse operators involve y0 = y(x0), which plays the role of arbitrary constant. Since our objective is education, we demonstrate three different approaches how to solve the given initial value problem based on the Adomian decomposition method. In the main streamline of the ADM, the solution y is represented as the infinite sum of series
\[ y(x) = u_0 (x) + \sum_{n\ge 1} u_n (x) , \]
and the nonlinear function is decomposed as follows:
\[ f(x,y) = \sum_{n\ge 0} A_n (x) , \]
where the n-th Adomian polynomial of \( u_0 , u_a , \ldots , u_n \) are calculated by 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( \sum_{k\ge 0} \lambda^k u_k \right) \right]_{\lambda =0} , \qquad n=0,1,2,\ldots . \]
These polynomials could be be obtained without the above Faà di Bruno's formula recursively according to Rach's rule:
\[ A_0 = f(u_0 ) , \qquad A_n = \sum_{k=1}^n C_n^k f^{(k)} (u_0 ) , \quad \mbox{for } n \ge 1, \]
where the coefficients Ckn are defined recursively following Jun-Sheng Duan algorithm as
\[ \begin{split} C_n^1 &= u_n , \quad n\ge 1, \\ C_n^k &= \frac{1}{n} \,\sum_{j=0}^{n-k} \left( j+1 \right) u_{j+1} C_{n-1-j}^{k-1} , \quad 2 \le k \le n . \end{split} \]
Upon substitution of the Adomian decomposition series for the solution y(x) and the series of Adomian polynomials tailored to the nonlinearity term f(x, y) and applying the inverse operator L-1, we obtain
\[ u_0 + \sum_{n\ge 1} u_n (x) = L^{-1} \left[ g(x) + \sum_{n\ge 0} A_n \left( u_0 , u_1 , \ldots , u_n \right) \right] . \]
In order to apply the fixed point theorem, we need the homogeneous initial conditions as well as the homogeneous equation. This is achieved by solving the initial value problem for u0:
\[ L \left[ u_0 \right] = g(x) , \qquad u_0 (x_0 ) = y_0 . \]
For other terms in the infinite series, we have the homogeneous conditions and find all term recursively:
\[ L \left[ u_{n+1} \right] = A_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_{n+1} (x_0 ) = 0, \qquad n=0,1,2,\ldots . \]
The generalized decomposition method follows the same procedure, but instead of Adomian polynomials it calls for another polynomials---called the accelerated Adomian polynomails
\[ \begin{split} J_0 &= A_0 = f(u_0 ) , \\ J_n &= f\left( u_0 + u_1 + \cdots , u_n \right) - f\left( u_0 + u_1 + \cdots , u_{n-1} \right) , \qquad n=1,2,\ldots . \end{split} \]
Then the solution's terms \( y(x) = \sum_{n\ge 0} u_n (x) \) are determined recursively
\[ L \left[ u_{n+1} \right] = J_n \left( u_0 + u_1 + \cdots , u_n \right) , \qquad u_{n+1} (x_0 ) =0 , \qquad n=0,1,2,\ldots . \]
Therefore, we consider three options:

Option 1: using ADM with the two-term operator: With u0(x) to be determined separately as the solution to the initial value problem

\[ L \left[ u_0 (x) \right] = g(x) , \qquad u_0 (x_0 ) = y_0 , \]
we find all other terms recursively
\[ u_{n+1} (x) = L^{-1} \left[ A_n \left( u_0 + u_1 + \cdots , u_n \right) \right] = x^{-p} \int_{x_0}^x t^p A_n (t)\, {\text d}t , \qquad n=0,1,2,\ldots . \]

Option 2: using accelerated Adomian polynomials for the two-term operator: The initial term is the same as in Option 1:

\[ u_0 (x) = L^{-1} \left[ g(x) \right] = y(x_0 ) \left( \frac{x_0}{x} \right)^p + x^{-p} \int_{x_0}^x t^p g(t) \,{\text d}t . \]
The unknown function representing the solution is assumed to be represented by an infinite decomposition series
\[ y(x) = u_0 + \sum_{n\ge 1} u_n (x) \]
where the solution components un, n > 1, are ro be determined recursively by solving the initial value problems:
\[ L \left[ u_{n+1} \right] = J_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_n (x_0 ) = 0, \qquad m=0,1,\ldots . \]
We list the formulas of the first several generalized polynomials for convenient reference as
\begin{align*} J_0 \left( u_0 \right) &= f(x, u_0 ) , \\ J_1 \left( u_0 , u_1 \right) &= f\left( x, u_0 + u_1 \right) - f\left( x, u_0 \right) , \\ J_2 \left( u_0 , u_1 , u_2 \right) &= f\left( x, u_0 + u_1 + u_2 \right) - f\left( x, u_0 + u_1 \right) , \\ J_3 \left( u_0 , u_1 , u_2 , u_3 \right) &= f\left( x, u_0 + u_1 + u_2 + u_3 \right) - f\left( x, u_0 + u_1 + u_2 \right) . \end{align*}

Option 3: using ADM with the isolated derivative: The initial term is expressed explicitly as

\[ u_0 (x) = y(x_0 ) + \int_{x_0}^x g(t) \, {\text d}t . \]
All other terms are determined recursively:
\[ u_{n+1} (x) = \int_{x_0}^x \left[ A_n (t) - \frac{p}{x}\, u_n (t) \right] \, {\text d}t , \qquad n=0,1,2,\ldots . \]
In the above formula, the Adomian polynomials Am are evaluated according to the same formula as in Option 1.

Example: Consider the initial value problem for nonlinear differential equation with the exponential slope:

\[ \dot{y} \overset{\mbox{def}}{=} \frac{{\text d} y}{{\text d}t} = \left( 2t-1 \right) e^{y-2t}, \qquad y(1) = 2 . \]
This problem has an explicit solution (a very rare case):
\[ y(t) = 2t - \ln t . \]
Of course, with Mathematica, it is easy to find its Taylor's series centered at t = 1:
Series[2*t -Log[t], {t, 1, 8}]
\[ y(t) = 2 + \left( t-1 \right) + \frac{1}{2} \left( t-1 \right)^2 - \frac{1}{3} \left( t-1 \right)^3 + \frac{1}{4} \left( t-1 \right)^4 -\frac{1}{5} \left( t-1 \right)^5 + \frac{1}{6} \left( t-1 \right)^6 -\frac{1}{7} \left( t-1 \right)^7 + \cdots . \]

Option 1: using classical ADM We seek the solution of the given initial value problem in the form of infinite series

\[ y(t) = u_0 (t) + \sum_{k\ge 1} u_k (t) , \]
where u0(t) is the solution of the initial value problem:
\[ \dot{u} = 0, \qquad u(1) = 2 \qquad \Longrightarrow \qquad u_0 = 2 . \]
All other terms are solutions of the homogeneose initial value problem
\[ \dot{u}_{m+1} = A_{m}, \qquad u(0) = 0, \]
where Am(u0, … , um) are the classical Adomian polynomials that correspond to the given slope function \( f(y) = \left( 2t-1 \right) e^{y-2t} . \) The first one is
\[ A_0 \left( u_0 \right) = f(2) = \left( 2t-1 \right) e^{2-2t} , \]
and all other polynomials are avaluated according to the formula:
\[ A_m \left( u_0 , \ldots , u_m \right) = \frac{1}{m!} \, \frac{{\text d}^m}{{\text d} \lambda^m} \left[ f \left( t, \sum_{k\ge 0} u_k \lambda^k \right) \right]_{\lambda =0} , \qquad m=0,1,2,\ldots . \]
For the given slope function, we have
\begin{eqnarray*} A_1 &=& u_1 f' (u_0 ) = u_1 \left( 2t-1 \right) e^{2-2t} = \left( 2t-1 \right) e^{2-2t} \left[1 - t\, e^{2-2t} \right] , \\ A_2 &=& u_2 f' (u_0 ) + \frac{1}{2!}\, u_1^2 f'' (u_0 ) = \left( 2t-1 \right) e^{2-2t} \left[ \frac{1}{2}\,u_1^2 + u_2 \right] = \left( 2t-1 \right) e^{2-2t} e^{-4t} \left( e^{2t} - e^2 t \right)^2 , \\ 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 ) = \left( 2t-1 \right) e^{2-2t} \left[ \frac{1}{6}\, u_1^3 + u_1 u_2 + u_3 \right] = \left( 2t-1 \right) e^{2-8t} \left[ e^{2t} + e^2 t \right]^3 , \\ 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^2 f'''' (u_0 ) = \left( 2t-1 \right) e^{2-2t} \left[ \frac{1}{24}\, u_1^4 + \frac{1}{2}\,u_1^2 u_2 + \frac{1}{2}\,u_2^2 + u_1 u_3 + u_4 \right] , \\ 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 \\ &=& \frac{1}{120}\left( 2t-1 \right) e^{2-2t} \left[ u_1^5 + \frac{1}{6}\,u_1^3 u_2 + \frac{1}{2}\,u_1 u_2^2 + \frac{1}{2}\,u_1^2 u_3 + u_2 u_3 + u_1 u_4 + u_5 \right] , \\ 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) + \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 \\ &=& \frac{1}{120}\left( 2t-1 \right) e^{2-2t} \left[ \frac{1}{720} \,u_1^6 + \frac{1}{24}\,u_1^4 u_2 = \frac{1}{4}\, U_1^2 U_2^2 + \frac{1}{6}\,U_1^3 U_3 + \frac{1}{6}\, u_2^3 + u_1 u_2 u_3 + \frac{1}{2}\,u_1^2 u_4 + \frac{1}{2}\, u_3^2 + u_2 u_4 + u_1 u_5 \right] , \\ 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) + \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) + \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) + \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*}
Now we are ready to calculate the terms in the infinite series for the solution:
\begin{eqnarray*} u_1 (t) &=& L^{-1} \left[ A_0 \right] = \int_1^t \left[ A_0 (s) \right] {\text d}s = \int_1^t \left( 2s-1 \right) e^{2-2s} {\text d}s = 1-t\,e^{2-2t} , \\ u_2 (t) &=& L^{-1} \left[ A_1 \right] = \int_1^t \left[ u_1 \left( 2s-1 \right) e^{2-2s} \right] {\text d}s = \frac{1}{2} -t\, e^{2-2t} + + \frac{t^2}{2}\, e^{4-4t} = \frac{1}{2}\, e^{-4t} \left( e^{2t} - e^2 t \right)^2 , \\ s_2 (t) &=& u_0 + u_1 + u_2 = \frac{3}{2} - e^{2-2t} + \frac{t^2}{2}\, e^{4-4t} , \\ u_3 (t) &=& L^{-1} \left[ A_2 \right] = \int_1^t A_2 (s)\,{\text d}s = \frac{1}{3}\, e^{-6t} \left( e^{2t} - e^2 t \right)^3 . \end{eqnarray*}
Mathematica confirms
Integrate[(2*s - 1)*Exp[2 - 2*s], {s, 1, t}]
1 - E^(2 - 2 t) t
Simplify[(1 - E^(2 - 2 t) t)*(2*t - 1)*Exp[2 - 2*t]]
E^(2 - 4 t) (-1 + 2 t) (E^(2 t) - E^2 t)
Integrate[(2*s - 1)*Exp[2 - 2*s]*(1 - s*Exp[2 - 2*s]), {s, 1, t}]
1/2 E^(-4 t) (E^(2 t) - E^2 t)^2
u0[t_] = (2*t - 1)*Exp[2 - 2*t];
u1[t_] = 1 - t*Exp[2 - 2*t];
u2[t_] = 1/2 - E^(2 - 2 t) t + 1/2 E^(4 - 4 t) t^2;
s2[t_] = Simplify[u0[t] + u1[t] + u2[t]]
1/2 (3 - 2 E^(2 - 2 t) + E^(4 - 4 t) t^2)
Now we find next terms:
A3[t_] = (2*t - 1)*Exp[2 - 2*t]*(u1[t]^3 /6 + u1[t]*u2[t] + u3[t]);
u4[t_] = Simplify[Integrate[A3[s], {s, 1, t}]]
1/4 E^(-8 t) (E^(2 t) - E^2 t)^4
\[ u_4 (t) = \frac{1}{4}\, e^{-8t} \left( e^{2t} - e^2 t \right)^4 . \]
Now we see the partern:
\[ u_n (t) = \frac{1}{n}\, e^{-2nt} \left( e^{2t} - e^2 t \right)^n , \qquad n= 1,2, \ldots . \]
Adding all terms, we obtain the required solution:
\[ y(t) = u_0 (t) + \sum_{n\ge 1} u_n (t) = 2 + \sum_{n\ge 1} \frac{1}{n}\, e^{-2nt} \left( e^{2t} - e^2 t \right)^n = 2 + \sum_{n\ge 1} \frac{1}{n}\, z^n , \]
where \( z = e^{-2t} \left( e^{2t} - e^2 t \right) = 1 - e^{2-2t} < 1. \) Using Taylor's series expansion
\[ -\ln \left( 1 - z \right) = \sum_{k\ge 1} \frac{z^k}{k} , \]
we obtain the explicit formula for the solution
\[ y(t) = 2 - \ln \left( 1 - z \right) = 2 - \ln \left( t\, e^{2-2t} \right) = 2 - \ln t - \ln \left( e^{2-2t} \right) = 2t - \ln t . \]

Option 2: using accelerated Adomian polynomials: The first two initial terms u0 = 2, \( u_1 = \int_1^t A_0 (s)\,{\text d}s = 1 - t\, e^{2-2t} \) and the Adomian polynomail A0 remains the same. The accelerated Adomian polynomails are

\[ J_m = \left( 2t-1 \right) e^{2-2t} \left[ e^{u_1 + \cdots + u_m} - e^{u_1 + \cdots + u_{m-1}} \right] = \left( 2t-1 \right) e^{2-2t} e^{u_1 + \cdots + u_{m-1}} \left[ e^{u_m} -1 \right] , \qquad m = 1,2, \ldots . \]
So we get some of them
\begin{align*} J_1 &= \left( 2t-1 \right) e^{2-2t} \left[ e^{u_1} - 1 \right] = \\ J_2 &= \left( 2t-1 \right) e^{2-2t} \left[ e^{u_1 + u_2} - e^{u_1} \right] = \end{align*}

Option 3: using MDM We assume that the required solution is represented by infinite Taylor's series

\[ y(t) = \sum_{n\ge 0} c_n \left( t-1 \right)^k = 2 + \sum_{n\ge 1} c_n \left( t-1 \right)^k , \]
and ther slope function can be decomposed in a similar series
\[ f(t,y) = \left( 2t -1 \right) e^{y -2t} = \sum_{n\ge 0} A_n \left( c_0 , c_1, \ldots , c_n \right)\left( t-1 \right)^k , \]
where An are Adomian's polynomails.    ■

Example: Consider another initial value problem for nonlinear differential equation with the exponential slope:

\[ \dot{y} \overset{\mbox{def}}{=} \frac{{\text d} y}{{\text d}t} = e^{2t}\, e^{-2y}, \qquad y(0) = \frac{1}{2} . \]
Note that we use the Newton's dot notation for the derivative with respect to time variable t. This problem has an explicit solution (a very rare case):
\[ y(t) = t + \frac{1}{2}\, \ln \left( 1 + \left( e-1 \right) e^{-2t} \right) . \]
Of course, with Mathematica, it is easy to find its Maclaurin series:
Series[t + (1/2)*Log[1 + (Exp[1] - 1)*Exp[-2*t]], {t, 0, 8}]
\[ y(t) = \frac{1}{2} + \frac{t}{e} + \frac{e-1}{e^2}\,t^2 + \frac{\left( e-2 \right)\left( e-1 \right)}{3\, e^3}\,t^3 + \frac{\left( e-1 \right)\left( 6-6e + e^2 \right)}{3\,e^4}\,t^4 + \frac{2 \left( e - 1 \right) \left( e^3 -14 e^2 + 36e -24 \right)}{15\,e^6}\,t^5 + \cdots . \]

Option 1: using classical ADM

Option 2: using accelerated Adomian polynomials:

Option 3: using MDM    ■

 

Inversion of singular differential operator


Let us consider the following singular differential operator
\[ L\left[ x, \texttt{D} \right] y(x) = \left( \texttt{D} + \frac{p}{x}\,\texttt{I} \right) y(x) = y' (x) + \frac{p}{x}\,y(x) , \]
where p is a real constant, \( \texttt{D} ≝ {\text d}/{\text d}x \) is the derivative operator, and \( \texttt{I} \) is the identity operator. The kernel (= set of solutions L y = 0; in other words, the kernel of a linear operator includes all functions that are annihilated by L) of the linear unbounded differential operator L is an one-dimensional space. To find its inverse L-1, we need to solve the differential equation
\[ L\left[ x, \texttt{D} \right] y(x) = w(x) \qquad \mbox{or}\qquad y' (x) + \frac{p}{x}\,y(x) = w(x). \]
Since this differential equation involves the derivative operator \( \texttt{D} ≝ {\text d}/{\text d}x , \) the above equation has infinite many solutions depending on arbitrary constant. Such family of solutions is called the general solution. From mathematical point of view, L-1 is not an operator because it assigns infinite many outputs to every input. To single out one of them, we consider the differential operator \( L \left[ x, \texttt{D} \right] \) on a set of smooth functions \( C^1 (a,b) \) on some open interval \( (a,b) \) that have a specific value at one of its point: y(x0) = y0. Then \( L \left[ x, \texttt{D} \right] : C^1 \to C^0 , \) and it has a unique inverse. Recall that we denote by C¹ the set of functions having continuous derivatives in a specified interval.

We start with a simple case when the differential operator has the only one singlar point at the origin:

\[ L\left[ x, \texttt{D} \right] = \texttt{D} + \frac{p}{x}\,\texttt{I} \qquad \mbox{or}\qquad L\left[ x, \texttt{D} \right] = x^{p-1} \,\texttt{D} x^p \,\texttt{I} . \]
Multiplying both side of the differential equation \( y; + p\,y/x = w \) by an integrating factor (xp in our case), we obtain an exact differential equation
\[ x^p \left( y' (x) + \frac{p}{x}\,y(x) \right) = \frac{\text d}{{\text d}x} \left( x^p y(x) \right) = x^p w(x) . \]
Its integration yields
\[ x^p y(x) = y(x_0 )\, x_0^p + \int_{x_0}^x t^p w(t)\,{\text d}t \qquad \Longrightarrow \qquad L^{-1} \left[ x, \texttt{D} \right] w(x) = y(x_0 ) \left( \frac{x_0}{x} \right)^p + x^{-p} \int_{x_0}^x t^p w(t)\,{\text d}t , \]
where x0 is an arbitrary number other than zero because the origin is a singular point for this differential equation. Here y0 plays the role of arbitrary constant to single out one inverse operator. So the singular point x0 = 0 divides the real axis ℝ into two parts, and the solution of the above differential equation \( y' +p\,y/x = f \) exists in either of these two parts. In future, we assume that x is positive and the initial condition is also specified at a positive point x0 > 0.

The situation can be extended for arbitrary function p(x); so we consider the first order differential operator:

\[ L\left[ x, \texttt{D} \right] y(x) = \left( \texttt{D} + p(x)\,\texttt{I} \right) y(x) = y' (x) + p(x)\,y(x) = e^{-\int p(x)\,{\text d}x} \, \frac{\text d}{{\text d}x} \left( e^{\int p(x)\,{\text d}x} \,y(x) \right) . \]
Its inverse is
\[ L^{-1} \left[ x, \texttt{D} \right] w(x) = y(x_0 ) \,e^{-\int p(x)\,{\text d}x} + e^{-\int p(x)\,{\text d}x} \int_{x_0}^x {\text d}t\,e^{\int p(t)\,{\text d}t} \, w(t) , \]

 

Singular differential equations


In this section we consider only differential equations with regular singular points of the form
\[ \frac{{\text d} y}{{\text d}x} + a(x) \, y(x) = f(x,y) , \]
where f is a holomorphic function in some two-dimensional domain and function 𝑎(x) has finite number of poles of first order. This means that there are finite number of m points xk, k = 1, 2, ⃛, m such that \( \left( x - x_k \right) a(x) \) is holomorphic.

Example: Consider the following singular differential equation:

\[ \frac{{\text d} y}{{\text d}x} + \frac{p}{x}\, y(x) = x^2 + y^2 (x), \qquad y(x_0 ) = y_0 . \]
Here x0 is any number other than zero and p is a real positive number. We rewrite the given differential equation in operator form:
\[ L\left[ y \right] = f(y) + g(x) , \]
where
\[ L\left[ y \right] = \frac{{\text d} y}{{\text d}x} + \frac{p}{x}\, y(x) , \qquad g(x) = x^2 , \quad f(y) = y^2 (x) . \]
The inverse of the unbounded differential operator L was found previously:
\[ L^{-1} \left[ u \right] = y(t_0 )\,x^{-p}x_0^p + x^{-p} \int_{x_0}^x t^p u(t)\,{\text d}t . \]
Application of the bounded operator to the given differential equation yields the functional equation for which the fixed point theorem can be applied:
\[ y(x) = L^{-1} \left[ g(x) \right] + L^{-1} \left[ f(y) \right] . \]
The Adomian decomposition method (ADM) states that the dependent variable y(x) and the nonlinear term f(y) should be written as the following infinite series
\begin{align*} y(x) &= \sum_{n\ge 0} u_k (x) , \\ f(y(x)) &= \sum_{n\ge 0} A_k \left( u_1 , u_2 , \ldots , u_k \right) , \end{align*}
where Ak is the k-th Adomian polynomial.

We demonstrate the Adomian decomposition method for solving the formulated initial value problem by utilizing three options.

Option 1: using ADM with two-term operator

Option 2: using GDM with two-term operator

Option 3: using ADM with isolated derivative

We rewrite the given differential equation by using only derivative operator:

\[ L_0 \left[ y \right] = f_0 (y) + g(x) , \]
where
\[ L_0 \left[ y \right] = \frac{{\text d} y}{{\text d}x} , \qquad g(x) = x^2 , \quad f_0 (y) = y^2 (x) - \frac{p}{x}\, y(x) . \]
The required solution should be written as the following infinite series
\[ y(x) = u_0 (x) + \sum_{n\ge 1} u_n (x) , \]
where u0 is the unique solution of the following initial value problem:
\[ u' (x) = g(x) , \qquad u(x_0 ) = y_0 . \]
So
\[ u_0 (x) = y(x_0 ) + \int_{x_0}^x g(x)\,{\text d}t . \]
To determine u1, we need to find the first Adomian polynomial:
\[ A_0 (x) = f_0 (u_0 ) = u_0^2 (x) - \frac{p}{x}\, u_0 (x) . \]
Then the first term u1 is the solution to the initial value problem
\[ u' = A_0 (x) = f_0 (u_0 ) = u_0^2 (x) - \frac{p}{x}\, u_0 (x) , \quad u(x_0 ) = 0 . \]
So it will be
\[ u_1 (x) = L^{-1}_0 \left[ A_0 (x) \right] = \int_{x_0}^x A_0 (t)\,{\text d}t = f_0 (u_0 ) \left( x- x_0 \right) \]
   ■
  1. Daftardar-Gejji, V. and Jafari, H., An iterative method for solving nonlinear functional equations, Journal of Mathematical Analysis and Applications, 2006, Vol. 316, No. 2, pp. 753--763. https://doi.org/10.1016/j.jmaa.2005.05.009
  2. Duan, J.-S. and Rach, R., The degenerate form of the Adomian polynomials in the power series method for nonlinear ordinary differential equations, Journal of Mathematics and System Science, Volume 5, Pages 411--428, doi: 10.17265/2159-5291/2015.10.003
  3. Jiao, Y.-C., Dang, C., Yamamoto, Y., An extension of the decomposition method for solving nonlinear equations and its convergence, Computers and Mathematics with Applications, 2008, Vol. 55, pp. 760--775

 

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)