Preface

This is a tutorial made solely for the purpose of education and it was designed for students taking Applied Math 0330. It is primarily for students who have very little experience or have never used Mathematica before and would like to learn more of the basics for this computer algebra system. As a friendly reminder, don't forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

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 APMA0340
Return to the main page for the course APMA0330
Return to Part III of the course APMA0330

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 an American mathematician and aerospace engineer of Armenian descent George Adomian (1922-1996), chair of the Center for Applied Mathematics at the University of Georgia. The crucial aspect of the method is the employment of the "Adomian polynomials" to represent the nonlinear portion of the equation as a convergent series with respect to these polynomials, without actual linearization of the system. These polynomials mathematically generalize Maclaurin series about an arbitrary external parameter, which gives the solution method more flexibility than direct Taylor series expansion.

Let me introduce some scientists (in alphabetic order) who contributed and greatly improved the Adomian Decomposition Method to make it available to mathematical and engineering community.

Dr. George Adomian (1922--1996) was an American mathematician, theoretical physicist, and electrical engineer of Armenian descent. He received his Ph.D. degree from UCLA. He first proposed and considerably developed the Adomian Decomposition Method (ADM) for solving nonlinear differential equations, ordinary, and partial. He was a professor of mathematics at the University of Georgia, the founder and chief scientist of General Analytics Corporation, a winner of the 1989 Richard Bellman Prize for outstanding contributions to nonlinear stochastic analysis, and a 1988 National Academy of Sciences Scholar. He is the author of eight books and over three hundred journal papers. Dr. Adomian was the Radar Officer (naval rank of Lieutenant, the equivalent of a Captain in the Army by direct commission) aboard the USS Antietam (CV-36), and served in the Pacific Theater of Operations during the Second World War. He also attended the secret radar school at MIT which trained the first radar officers for the U.S. Navy.

Dr. Yves Cherruault (1937--2010) was a French mathematician. He received his Ph.D. degree from the University of Paris. He was a professor of mathematics at the University Pierre et Marie Curie, Paris, and the Director of MÉDIMAT (Laboratory of Mathematics Applied to Biomedicine). Professor Cherruault is one of the founders of the field of Biomathematics and an author of seven books and over two hundred journal papers. Dr. Cherruault developed some of the convergence theorems of the ADM.

Jun-Sheng Duan.
Dr. Jun-Sheng Duan is a Chinese applied mathematician and computer scientist born in Inner Mongolia, China, in 1965. He received his Ph.D. degree from Shandong University, China. He has made extensive contributions to solutions of differential equations in mathematical physics and engineering using ADM and the Modified Decomposition Method (MDM) in collaboration with Drs. Rach and Wazwaz. He has been a professor of mathematics at the Inner Mongolia Polytechnic University, the Tianjin University of Commerce, and currently at the School of Sciences, Shanghai Institute of Technology, China. He is the author of more than eighty journal papers. Dr. Duan is an outstanding player of wei qi (the Chinese “game of go”), and ping-pong (table tennis). He enjoys photography and traveling.










Randolph Rach.
Dr. Randolph Rach is a retired Army veteran and formerly a Senior Engineer at Microwave Laboratories Inc. His experience includes research and development in microwave electronics and traveling-wave tube technology, and his research interests span nonlinear system analysis, nonlinear ordinary and partial differential equations, nonlinear integral equations and nonlinear boundary value problems. He published more than one hundred and thirty papers in applied mathematics, and he is an early contributor to the ADM. Dr. Rach’s fundamental theorems established the basis for the early development of ADM. He was also the first to propose the MDM. His past and current research continue to advance ADM/MDM work.












Sergio E. Serrano.

Dr. Sergio E. Serrano was born in Santander (Colombia) in 1953. He is a professor of environmental engineering, hydrologic science, applied mathematics, and philosophy at Temple University in Philadelphia. He received his Ph.D. degree at the University of Waterloo (Canada). He has more than one hundred research publications in international science, engineering, and mathematical journals. He is also the author of eight books in environmental engineering, statistics, philosophy, and psychology. Dr. Serrano has been awarded four times with nationally-competitive research grants by the National Science Foundation, Washington, DC. Using adaptations and modifications of the ADM, he has developed hundreds of practical engineering models of flood wave propagation, contaminant transport, and groundwater flow in irregular geometries. He published influential textbooks on applications of differential equations in CAS Maple and Hydrology:

Dr. Serrano has a passion for hiking. In 1973, he explored the rain forest of the Darien Sierra on foot from Acandi (Colombia) to Yaviza (Panama). He plays the recorder (Renaissance flute). He enjoys alchemy, archaeology, home wine making, herbal medicine, and cooking. He believes that the joy of meaningful living and meaningful relationships can be found in the simplicity of everyday life. He lives in Philadelphia with his wife of thirty years and his daughter.




Abdul-Majid Wazwaz.
Dr. Abdul-Majid Wazwaz is a Professor of Mathematics at Saint Xavier University in Chicago. He received his Ph.D. from the University of Illinois at Chicago. He was the author and co-author of more than four hundred and fifty papers in applied mathematics and mathematical physics. He is the author of five books on the subjects of discrete mathematics, integral equations, and partial differential equations. He has contributed extensively to theoretical advances in solitary waves theory, the ADM, and other computational methods. He is a member of the editorial board of the journals Nonlinear Dynamics (Springer) and Physica Scripta (IOP). For these three years in a row, 2014, 2015, and 2016, Thomas Reuters granted him three different badges for being a “Highly Cited Researcher.”

We show how the method works in the series of examples. Since the Adomian method has its roots in quantum mechanics and operator methods, we need to reconsider our old friend---the derivative operator. It is convenient to set a special symbol for the derivative operator, which is naturally denoted by \( \texttt{D} = {\text d} / {\text d} x . \) Its inverse is studied in calculus, and we denote it as

\[ \texttt{D}^{-1} y (x) = y_0 + \int_0^x y(s)\,{\text d} s , \]
where the initial point of integration, s = 0 was chosen for simplicity (any point from the domain of y(s) will work) and \( y(0) = y_0 . \) Of course, the derivative operator is defined on the space of smooth functions and its inverse acts in the space of integrable (for instance, in Riemann sense) functions. It should be noted that such defined operator \( \texttt{D}^{-1} \) is only right inverse to the derivative operator: \( \texttt{D}\,\texttt{D}^{-1} = \texttt{I} , \) which is the identity operator. Note that application of \( \texttt{D}^{-1} \) to the derivative yields
\[ \texttt{D}^{-1} \texttt{D}\,y = \texttt{D}^{-1} y' (x) = y' (0) + \int_0^x y'(s)\,{\text d} s = y' (0) + y(x) - y(0) . \]

Example: Let us take a constant function \( y(x) \equiv 1 , \) then

\[ \texttt{D}^{-1} y (x) = \texttt{D}^{-1} \,1 = C + x , \]
where C is a constant of integration. The derivative of y(x) is identically zero, so
\[ \texttt{D}^{-1} y' (x) = \texttt{D}^{-1} \texttt{D}\,1 = \texttt{D}^{-1} 0 = 0 \ne 1 = \texttt{D} \, \texttt{D}^{-1} 1 = \texttt{D} \,(C+x ) . \]
So we see that \( \texttt{D} \) and \( \texttt{D}^{-1} \) do not commute. This means that regular indefinite integral with arbitrary constant is only left inverse of the derivative operator. ■

Example: We start with a simple initial value problem for linear equation:

\[ \frac{{\text d}y}{{\text d}x} + a\,y =0 , \qquad y(0) = y_0 , \]
where a is a constant. Introducing the derivative operator \( \texttt{D} = {\text d} / {\text d} x , \) we rewrite the linear differential equation as
\[ \texttt{D}\,y + a\,y =0 , \qquad \Longrightarrow \qquad y = - \texttt{D}^{-1} ay = -a\, \texttt{D}^{-1} y , \]
where
\[ \texttt{D}^{-1} y (x) = y_0 + \int_0^x y(s)\,{\text d} s \]
Assuming that the solution is the infinite sum
\[ y(x) = y_0 + y_1 + y_2 + \cdots = \sum_{k\ge 0} y_k , \]
we substitute this series into either the differential equation
\[ \texttt{D}\,\sum_{k\ge 0} y_k + a\,\sum_{k\ge 0} y_k = \sum_{k\ge 1} \texttt{D}\,y_k + a\,\sum_{k\ge 0} y_k = 0 , \]
or the formula
\[ \sum_{k\ge 0} y_k = -a\, \texttt{D}^{-1} \sum_{k\ge 0} y_k = -a \, \sum_{k\ge 0} \texttt{D}^{-1} y_k . \]
Comparing like terms, we get the recurrence relation
\begin{align*} \texttt{D} y_0 &= 0 \qquad \mbox{subject to}\quad y_0 (0) = y_0 , \\ \texttt{D} y_{k+1} + a\, y_k &=0 \qquad \mbox{or} \qquad y_{k+1} = -a\, \texttt{D}^{-1} y_k , \quad k=0,1,2,\ldots , \end{align*}
subject to \( y_{k+1} (0) =0 . \) Each of these recurrences is actually a linear difference-differential equation, which is easy to solve:
\begin{align*} y_1 &= -a\, \texttt{D}^{-1} y_0 = -ax\,y_0 , \qquad \mbox{with } y_1 (0) = 0 , \\ y_2 &= -a\, \texttt{D}^{-1} y_1 = a^2 y_0 \,\texttt{D}^{-1} x = a^2 \,\frac{x^2}{2} \, y_0 , \qquad \mbox{with } y_2 (0) = 0 , \\ y_3 &= -a\, \texttt{D}^{-1} y_2 = - a^3 \,\frac{x^3}{3!} \, y_0 , \qquad \mbox{with } y_3 (0) = 0 , \end{align*}
and so on. Adding all these terms, we get the solution
\[ y(x) = y_0 \left( 1- ax + \frac{a^2x^2}{2!} - \cdots \right) = y_0 e^{-ax} . \]
Therefore, we have the solution expressed in closed form: \( y(x) = y_0 e^{-ax} . \) ■

Example: Our next example is also related to the linear differential equation, but now we consider variable coefficient nonhomogeneous equation:

\[ \texttt{D}\, y + 2x\,y =x, \qquad y(1) =3 . \]
Expanding the unknown solution y as \( y(x) = \sum_{k\ge 0} y_k , \) and substituting into the differential equation, we obtain
\[ \texttt{D}\left( y_0 + y_1 + y_2 + \cdots \right) = x - 2x \left( y_0 + y_1 + y_2 + \cdots \right) . \]
Clearly y0 in the above equation is a solution of the initial value problem:
\[ \texttt{D}\, y =x, \qquad y(1) =3 . \]
This is a calculus problem and Mathematica easily solves it:
DSolve[{y'[x] == x, y[1] == 3}, y[x], x]// Flatten
Out[2]= {y[x] -> 1/2 (5 + x^2)}
Therefore, \( y_0 (x) = \frac{5}{2} + \frac{x^2}{2} . \) Any other subsequent term is derived from the recurrence:
\[ \texttt{D}\, y_n =-2x\, y_{n-1} \qquad y_n (1) =0 ; \qquad n=1,2,3,\ldots . \]
We solve this recurrence using Mathematica:
y[0, x_] = f[x] /. First @ DSolve[{f'[x] == x, f[1] == 3}, f, x]
Table[y[n, x], {n, 0, 6}]
Do[y[n, x_] =
f[x] /. First @ DSolve[{f'[x] == -2*x*y[n - 1, x], f[1] == 0}, f, x], {n, 1, 5}]
This gives first five approximations:
\[ y_1 = \frac{1}{4} \left( 11- 1x^2 - x^4 \right) , \quad y_2 = \frac{1}{12} \left( 17 - 33x^2 +15 x^4 + x^6 \right) , \quad y_3 = \frac{1}{48} \left( 23 - 68x^2 + 66 x^4 -20x^6 -x^8 \right) , \quad y_4 = \frac{1}{240} \left( 29 - 115 x^2 + 170 x^4 -110 x^6 + 25 x^8 + x^{10} \right) , \]
and \( y_5 = \frac{1}{1440} \left( 35 - 174 x^2 + 345 x^4 -340 x^6 + 165 x^8 - 30 x^{10} - x^{12} \right) . \) Adding these terms, we get 5-term approximation to the actual solution \( y(x) = \frac{1}{2} + \frac{5}{2}\,e^{1-x^2} : \)
\[ \phi_5 (x) = y_0 + y_1 + y_2 + y_3 + y_4 + y_5 = \frac{1}{1440} \left( 10499 -9744 x^2 + 4785 x^4 - 1480 x^6 + 285 x^8 - 24 x^{10} - x^{12} \right) . \]
We plot the actual solution in blue and approximation in orange:
z[x_] := 1/2 + (5/2)*Exp[1 - x^2]
phi[x_] := Simplify[Sum[y[n, x], {n, 0, 5}]]
Plot[{z[x], phi[x]}, {x, 1, 2}, PlotStyle -> {{Thick, Blue}, {Thick, Orange}}]
From the above graph, we see that the ADM solution \( \phi_5 (x) \) is almost indistinguishable from the exact within the interval [1, 1.5]. However, for \( x \ge 1.6 , \) the ADM solution begins to diverge from the exact one. This is typical: the convergence rate is a function of the window of observation. If number of terms in ADM approximation is taken more, the approximate solution improves at high values x. From here we conclude that as x increases, the number of ADM terms needed to obtain an accurate solution increases. ■

Now we turn our attention to initial value problems for nonlinear equations such as

\[ \frac{{\text d}y}{{\text d}t} + a(t)\,y(t) + b(t) \, N[y] = f(t), \qquad y(0) = y_0 , \]
where we have chosen the initial point to be the origin for simplicity. Here a(t), b(t), and f(t) are known and bounded functions in variable t, and N[y] is a nonlinear (analytic) operator acting on dependent variable (function) y(t). As before, we define the derivative operator \( \texttt{D} = {\text d}/ {\text d}t \) with respect to independent variable t and rewrite the given initial value problem in operator form:
\[ \texttt{D} \,y + a(t)\,y(t) + b(t) \, N[y] = f(t), \qquad y(0) = y_0 , \]
The main idea of the Adomian decomposition method is to represent the nonlinear term as the sum of so called Adomian polynomials:
\[ N[y] = \sum_{n\ge 0} A_n (y_0 , y_1 , \ldots , y_n) , \]
where An are the Adomian polynomials specifically generated for each nonlinear operator according to the formula:
\[ A_n (y_0 , y_1 , \ldots , y_n) = \left. \frac{1}{n!} \, \frac{{\text d}^n}{{\text d} \lambda^n} \, N \left[ \sum_{k=0}^n y_k \lambda^k \right] \right\vert_{\lambda =0} , \qquad n=0,1,2,3,\ldots . \]
In particular,
\begin{align*} A_0 &= N[y_0 ] , \\ A_1 &= y_1 \,\frac{{\text d}\, N[y_0 ]}{{\text d} y_0} , \\ A_2 &= y_2 \,\frac{{\text d}\, N[y_0 ]}{{\text d} y_0} + \frac{y_1^2}{2!} \,\frac{{\text d}^2 \, N[y_0 ]}{{\text d} y_0^2} , \\ A_3 &= y_3 \,\frac{{\text d}\, N[y_0 ]}{{\text d} y_0} + y_1 y_2 \,\frac{{\text d}^2 \, N[y_0 ]}{{\text d} y_0^2} + \frac{y_1^3}{3!} \,\frac{{\text d}^3 \, N[y_0 ]}{{\text d} y_0^3} , \\ A_4 &= y_4 \,\frac{{\text d}\, N[y_0 ]}{{\text d} y_0} + \left( y_1 y_3 + \frac{y_2^2}{2!} \right) \frac{{\text d}^2 \, N[y_0 ]}{{\text d} y_0^2} + \frac{y_1^2 y_2}{2!} \,\frac{{\text d}^3 \, N[y_0 ]}{{\text d} y_0^3} + \frac{y_1^4}{4!} \,\frac{{\text d}^4 \, N[y_0 ]}{{\text d} y_0^4} , \\ A_5 &= y_5 \,\frac{{\text d}\, N[y_0 ]}{{\text d} y_0} + \left( y_1 y_4 + y_2 y_3 \right) \frac{{\text d}^2 \, N[y_0 ]}{{\text d} y_0^2} + \frac{1}{2} \left( y_1 y_2^2 + y_1^2 y_3 \right) \frac{{\text d}^3 \, N[y_0 ]}{{\text d} y_0^3} + \frac{y_1^3 y_2}{3!} \,\frac{{\text d}^4 \, N[y_0 ]}{{\text d} y_0^4} + \frac{y_1^5}{5!} \,\frac{{\text d}^5 \, N[y_0 ]}{{\text d} y_0^5} , \\ A_6 &= y_6 \,\frac{{\text d}\, N[y_0 ]}{{\text d} y_0} + \left( y_2 y_4 + \frac{y_3^2}{2!} + y_1 y_5 \right) \frac{{\text d}^2 \, N[y_0 ]}{{\text d} y_0^2} + \left( \frac{y_1^2 y_4}{2!} + \frac{y_2^3}{3!} + y_1 y_2 y_3 \right) \frac{{\text d}^3 \, N[y_0 ]}{{\text d} y_0^3} + \cdots . \end{align*}
For the approximations of the solution \( y = u_0 + u_1 + u_2 + \cdots \) and the nonlinear operator \( N[y] = A_0 + A_1 + A_2 + \cdots \) correspond ordinary generating functions
\[ Y( t, \lambda ) = \sum_{n\ge 0} u_n (t) \,\lambda^n \qquad\mbox{and} \qquad N[Y(t,\lambda )] = \sum_{n\ge 0} A_n \left( u_0 (t) , u_1 (t) , \ldots , u_n (t) \right) \lambda^n . \]
where λ is a parameter. The decomposition method consists in identifying the un's by means of the formulae
\begin{align*} \texttt{D}\, u_0 &= f(t ), \qquad u_0 (0) = y_0 , \\ \texttt{D}\, u_n &= -a(t) \,u_{n-1} - b(t)\,A_{n-1} \left( u_0 , u_1 , \ldots , u_{n-1} \right) , \qquad u_n (0) =0 , \quad n=1,2,\ldots . \end{align*}

Example: We apply Adomian decomposition method to find approximate solution to the initial value problem for logistic equation:

\[ \dot{P} = 0.02\,P(1-P) , \qquad P(0) = P_0 = 10. \]
This problem has a unique solution
\[ P(t) = \frac{10\, e^{t/50}}{10\, e^{t/50} -9} , \qquad t\ge 0, \]
as Mathematica confirms
DSolve[{P'[t] == 2/100*P[t]*(1 - P[t]), P[0] == 10}, P[t], t]
We expand the nonlinear term N[P] = P2 into the series:
\[ N[P] = P^2 = \sum_{n\ge 0} A_n (P_0 , P_1 , \ldots P_n ) , \]
where first five polynomials are not hard to determine manually:
\begin{align*} A_0 &= y_0^2 = 10^2 = 100, \\ A_1 &= y_1 \,2\, y_0 , \\ A_2 &= y_2 \,2\, y_0 + y_1^2 , \\ A_3 &= y_3 \,2\, y_0 + 2\,y_1 y_2 , \\ A_4 &= y_4 \,2\, y_0 + 2\, y_1 y_3 + y_2^2 , \\ A_5 &= y_5 \,2\, y_0 + 2\,y_1 y_4 + 2\,y_2 y_3 . \end{align*}
Substituting expansion \( P = u_0 + u_1 + u_2 + \cdots \) into the given logistic equation, we get
\[ \texttt{D} \left( u_0 + u_1 + u_2 + \cdots \right) = 0.02 \left( u_0 + u_1 + u_2 + \cdots \right) - 0.02 \left( A_0 + A_1 + A_2 + \cdots \right) . \]
For the first term, we have the initial value problem
\[ \texttt{D} \, u_0 = 0, \qquad u_0 (0) = 10 . \]
Solving the initial value problem for the above linear equation we get the first term:
\[ u_0 (t) = 10 . \]
For next terms, we have
\[ \texttt{D} \, u_n \equiv \dot{u}_n = 0.02 \, u_{n-1} - 0.02\,A_{n-1} , \qquad u_n (0) = 0 , \quad n=1,2,3,\ldots ; \]
which is a linear differential equation with respect to un. Solving these equations recursively, we obtain
\begin{align*} u_1 &= -\frac{9}{5}\,t , \qquad &A_1 = -36\,t , \\ u_2 &= \frac{171}{500} \,t^2 , \qquad &A_2 = \frac{252}{25}\,t^2 , \\ u_3 &= -\frac{1623}{25000} \,t^3 , \qquad &A_3 = -\frac{1581}{625} \, t^3 , \\ u_4 &= \frac{61617}{5000000} \,t^4 , \qquad &A_4 = \frac{74643}{125000} \, t^4 , \\ u_5 &= -\frac{2924103}{1250000000} \, t^5 , \qquad &A_5 = - \frac{4236099}{31250000} \, t^5 . \end{align*}
u[0, t] := 10
A[0, t] := 100
u[1, t_] = f[t] /. First@ DSolve[{f'[t] == (1/50)*u[0,t] - (1/50)*100, f[0] == 0}, f, t]
A[1, t] = -2*u[0,t]*u[1,t]
u[2, t_] = f[t] /. First@ DSolve[{f'[t] == (1/50)*u[1, t] - (1/50)*A[1, t], f[0] == 0}, f, t]
A[2,t_]=20*u[2,t] + (u[1,t])^2
u[3, t_] = f[t] /. First@ DSolve[{f'[t] == (1/50)*u[2, t] - (1/50)*A[2, t], f[0] == 0}, f, t]
A[3,t_]= 20*u[3,t] + 2*u[1,t]*u[2,t]
u[4, t_] = f[t] /. First@ DSolve[{f'[t] == (1/50)*u[3, t] - (1/50)*A[3, t], f[0] == 0}, f, t]
A[4,t_]= 20*u[4,t] + (u[2,t])^2 + 2*u[1,t]*u[3,t]
u[5, t_] = f[t] /. First@ DSolve[{f'[t] == (1/50)*u[4, t] - (1/50)*A[4, t], f[0] == 0}, f, t]
Then we plot the exact solution in blue and ADM approximation in orange:
z[t_] = 10*Exp[t/50]/(10*Exp[t/50] - 9)
phi[t_] := 10 - (9 t)/5 + (171 t^2)/500 - (1623 t^3)/25000 + (61617 t^4)/ 5000000 - (2924103 t^5)/1250000000
Plot[{z[t], phi[t]}, {t, 0, 0.5}, PlotStyle -> {{Thick, Blue}, {Thick, Orange}}]
As we see, the ADM approximation with six terms (green) is more accurate to the actual solution (in blue) than the ADM approximation with five terms (in red).

Example: Consider the initial value problem for the integro-differential equation

\[ \dot{y} = 2.3\, y -0.01\,y^2 -0.1\,y\, \int_0^t y(\tau )\,{\text d}\tau , \qquad y(0) =50. \]
Using ADM method, we seek approximate solution in the form
\[ y = u_0 + u_1 + u_2 + \cdots , \]
where \( u_0 =50 \) and nonlinear operator N[y] is expanded into Adomian polynomials:
\[ N[y] \equiv -0.01\,y^2 -0.1\,y\, \int_0^t y(\tau )\,{\text d}\tau = A_0 + A_1 + A_3 + \cdots . \]
Since
\[ \frac{{\text d}\, N[y ]}{{\text d} y} = -0.02\,y - 0.1\,y^2 - 0.1\,\int_0^t y(\tau )\,{\text d}\tau , \quad \frac{{\text d}^2\, N[y ]}{{\text d} y^2} = -0.02 -0.3\,y , \quad \frac{{\text d}^3\, N[y ]}{{\text d} y^3} = -0.3 , \quad \frac{{\text d}^k\, N[y ]}{{\text d} y^k} = 0, \quad k= 4,5 , \ldots , \]
we have
\begin{align*} A_0 &= -0.01\,50^2 - 0.1\,50\, \int_0^t 50\,{\text d}\tau = -25 -250\,t , \\ A_1 &= -u_1 \left( 1+ 5t \right) - 5\,\int_0^t u_1 (\tau ) \,{\text d}\tau , \\ A_2 &= -u_2 \left( 1+ 5t \right) - \frac{u_1^2}{100} - \frac{u_1}{10} \, \int_0^t u_1 (\tau ) \,{\text d}\tau -5 \int_0^t u_2 (\tau ) \,{\text d}\tau , \\ A_3 &= - u_3 \left( 1+ 50t \right) -\frac{1}{50}\,u_1 u_2 - \frac{u_2}{10} \, \int_0^t u_1 (\tau ) \,{\text d}\tau - \frac{u_1}{10} \, \int_0^t u_2 (\tau ) \,{\text d}\tau - 50 \int_0^t u_3 (\tau ) \,{\text d}\tau , \\ A_4 &= -u_4 \left( 1+ 5t \right) -\frac{1}{100}\,u_2^2 - \frac{1}{50}\,u_1 u_3 - \frac{u_3}{10} \, \int_0^t u_1 (\tau ) \,{\text d}\tau - \frac{u_2}{10} \, \int_0^t u_2 (\tau ) \,{\text d}\tau - \frac{u_1}{10} \, \int_0^t u_3 (\tau ) \,{\text d}\tau -5 \int_0^t u_4 (\tau ) \,{\text d}\tau , \\ A_5 &= - u_5 \left( 1+ 50t \right) -\frac{1}{50}\,u_2 u_3 -\frac{1}{50}\,u_1 u_4 - \frac{u_4}{10} \, \int_0^t u_1 (\tau ) \,{\text d}\tau - \frac{u_3}{10} \, \int_0^t u_2 (\tau ) \,{\text d}\tau - \frac{u_2}{10} \, \int_0^t u_3 (\tau ) \,{\text d}\tau - \frac{u_1}{10} \, \int_0^t u_4 (\tau ) \,{\text d}\tau -50 \int_0^t u_5 (\tau ) \,{\text d}\tau . \end{align*}
Now we are ready to construct 6-term ADM approximation
\[ {\texttt D} u_n = 2.3\,u_{n-1} + A_{n-1}, \quad u_n (0) =0, \qquad n=1,2,3,\ldots . \]
Using Mathematica, we perform every step explicitly (without automatic looping).
u[0, t_] := 50
A[0, t_] = -1/100*50^2 - 1/10*50*Integrate[50, {x, 0, t}]
u[1, t_] = f[t] /. First@ DSolve[{f'[t] == 23*5 + A[0, t], f[0] == 0}, f, t]
Out[20]= -25 - 250 t
Out[21]= -5 (-18 t + 25 t^2)
For the second stage, we get
A[1, t_] = -Simplify[u[1, t]*(1 + 5*t) + 5*Integrate[u[1,x], {x, 0, t}]]
u[2, t_] = f[t] /. First@ DSolve[{f'[t] == 23*u[1,t]/10 + A[1, t], f[0] == 0}, f, t]
Out[22]= -(10/3) t (27 + 165 t - 250 t^2)
Out[23]= 1/6 (351 t^2 - 1675 t^3 + 1250 t^4)
A[2, t_] = -Simplify[u[2, t]*(1 + 5*t) + (u[1,t])^2 /100 + 5*Integrate[u[2,x], {x, 0, t}] + u[1,t]*Integrate[u[1,x], {x, 0, t}]/10]
u[3, t_] = f[t] /. First@ DSolve[{f'[t] == 23*u[2,t]/10 + A[2, t], f[0] == 0}, f, t]
Out[24]= -(1/24) t^2 (3348 + 6980 t - 55625 t^2 + 42500 t^3)
Out[25]= 1/720 (-1188 t^3 - 167925 t^4 + 402750 t^5 - 212500 t^6)
Then on the third stage, we have
A[3, t_] = -Simplify[u[3, t]*(1 + 50*t) + u[1, t]*u[2, t]/50 + u[2,t]*Integrate[u[1,x], {x, 0, t}]/10 + u[1,t]*Integrate[u[2,x], {x, 0, t}]/10 + 50 *Integrate[u[3,x], {x, 0, t}] ]
u[4, t_] = f[t] /. First@ DSolve[{f'[t] == 23*u[3, t]/10 + A[3, t], f[0] == 0}, f, t]
Out[26]= (t^3 (-522396 + 2753625 t + 74256000 t^2 - 177218125 t^3 + 92000000 t^4))/5040
Out[27]= (-1353807 t^4 + 100065 t^5 + 134567125 t^6 - 258056250 t^7 + 115000000 t^8)/50400
Finally, we obtain
A[4, t_] = -Simplify[u[4, t]*(1 + 5*t) + (u[2,t])^2 /100 + u[3,t]*u[1,t]/50 + u[3,t]*Integrate[u[1,x], {x, 0, t}]/10 + u[2,t]*Integrate[u[2,x], {x, 0, t}]/10 + u[1,t]*Integrate[u[3,x], {x, 0, t}]/10 + 5*Integrate[u[4,x], {x, 0, t}] ]
u[5, t_] = f[t] /. First@ DSolve[{f'[t] == 23*u[4, t]/10 + A[4, t], f[0] == 0}, f, t]
Out[28]= -((t^4 (7967484 - 1448901972 t + 5233879350 t^2 + 26405797500 t^3 - 60180328125 t^4 + 27762500000 t^5))/1814400)
Out[29]= (1/45360000)(-600313518 t^5 + 6071613975 t^6 + 21100995000 t^7 - 149290171875 t^8 + 193617578125 t^9 - 69406250000 t^10)
Now we add these five terms and plot two ADM approximations
phi5[t_] = u[0, t]
Do[phi5[t] = phi5[t] + u[j, t], {j, 1, 5}]
A[5, t_] = -Simplify[ u[5, t]*(1 + 50*t) + u[2, t]*u[3, t]/50 + u[4, t]*u[1, t]/50 + u[4, t]*Integrate[u[1, x], {x, 0, t}]/10 + u[3, t]*Integrate[u[2, x], {x, 0, t}]/10 + u[2, t]*Integrate[u[3, x], {x, 0, t}]/10 + u[1, t]*Integrate[u[4, x], {x, 0, t}]/10 + 50*Integrate[u[5, x], {x, 0, t}]]
u[6, t_] = f[t] /. First@ DSolve[{f'[t] == 23*u[5, t]/10 + A[5, t], f[0] == 0}, f, t]
phi6[t_] = phi5[t]+u[6,t]
Plot[Tooltip@{phi5[t], phi6[t]}, {t, 0, 0.8}, PlotStyle -> {Red, Green}, PlotLegends -> {"ADM(5)", "ADM(6)"}]

 

 

 

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)