Preface

This section is devoted to applications of the modified decomposition method (MDM for short) to the second order differential equations.

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

Glossary

Modified Decomposition Method 2

The modified decomposition method can be extended for higher order differential equations. We consider the initial values problems for the following second order differential equation \begin{equation} y'' + a(x)\, y + b(x)\, y' = g(x) + f(y) , \qquad y\left( x_0 \right) = \alpha , \quad y'\left( x_0 \right) = \beta , \label{Eadm.8} \end{equation} where 𝑎(x), b(x), and g(x) are holomorphic functions (in applications, these functiosn are usually polynomails) in some neighborhood of the initial point x0, but f(y) is a holomorphoc function in some open domain including points α and β. It is assumed that we know Taylor series expansions of the former functions: \begin{align*} a(x) = \sum_{n\ge 0} a_n \left( x - x_0 \right)^n , \\ b(x) = \sum_{n\ge 0} b_n \left( x - x_0 \right)^n , \\ g(x) = \sum_{n\ge 0} g_n \left( x - x_0 \right)^n . \end{align*} The solution to the initial value problem \eqref{Eadm.8} is asuumed to be represented by a convergent power series \begin{equation} y(x) = \sum_{n\ge 0} u_n (x) = c_0 + c_1 \left( x - x_0 \right) + \sum_{n\ge 2} c_n \left( x - x_0 \right)^n , \qquad c_0 = \alpha , \quad c_1 = \beta . \label{Eadm.1} \end{equation} For application of the modified decomposition method, we also need to know the expansion of the nonlinearity f(y) into the Adomian--Rach series: \[ f(y) = \sum_{n\ge 0} A_n \left( u_0 , u_1 , \ldots , u_n \right) = \sum_{n\ge 0} A_n \left( c_0 , c_1 , \ldots , c_n \right) \left( x - x_0 \right)^n , \] where the Adomian polynomials depend on the components of the solution series \eqref{Eadm.1}. Substituting each power series into the given differential equation and using the convolution formula, we obtain \begin{align*} &\sum_{n\ge 0} \left( n+2 \right) \left( n+1 \right) c_{n+2} \left( x - x_0 \right)^n + \sum_{n\ge 0} \left( a \ast y \right)_n \left( x - x_0 \right)^n + \sum_{n\ge 0} \left( b \ast y' \right)_n \left( x - x_0 \right)^n \\ &= \sum_{n\ge 0} g_n \left( x - x_0 \right)^n + \sum_{n\ge 0} A_n \left( x - x_0 \right)^n \end{align*} where \[ \left( a \ast y \right)_n = \sum_{k=0}^n a_{n-k} c_k \] is the n-component of the convolution of two series (as a result of multiplication of two series). Matching the coefficients of these series, we derive the full-history recurrence \begin{equation} \left( n+2 \right) \left( n+1 \right) c_{n+2} + \sum_{k=0}^n a_{n-k} c_k + \sum_{k=0}^n b_{n-k} c_{k+1} \left( k+1 \right) = g_n + A_n , \qquad n=0,1,2,\ldots . \label{Eadm.10} \end{equation} The initial two terms must be chosen in accordance to the given initial conditions \[ c_0 = \alpha , \qquad c_1 = \beta . \] We demonstrate this technique in the following examples.

Example: The Gross–Pitaevskii equation (GPE, named after Eugene P. Gross and Lev Petrovich Pitaevskii) describes the ground state of a quantum system of identical bosons using the Hartree–Fock approximation and the pseudopotential interaction model. Its simplified stationary version is read as

\[ - \frac{{\text d}^2 y}{{\text d}t^2} + \left( 1 - \frac{3}{\cosh^2 t} \right) y(t) + y^3 (t) = 0, \qquad y(0) =1, \quad \dot{y}(0) =0. \]
The given initial value problem admits the explicit solution y(t) = sech(t).

We first compute the Adomian polynomails for noninear part f(y) = y³:

\begin{align*} A_0 &= u_0^3 , \\ A_1 &= 3\,u_0^2 u_1 , \\ A_2 &= 3\,u_0 u_1^2 + 3\,u_0^2 u_2 , \\ A_3 &= 3\,u_3 u_0^2 + 6\, u_0 u_1 u_2 + u_1^3 . \end{align*}
The general formular is alos avaiulable as a double convolution:
\[ A_k = \sum_{i=0}^k \sum_{j=0}^{k-i} u_i u_j u_{k-i-j} , \qquad k=0,1,2,\ldots . \]
The solution is assumed to be decomposed as
\[ y(t) = u_0 (t) + \sum_{k\ge 1} u_k (t) = 1 + \sum_{k\ge 1} u_k (t) , \]
where the initial term is u0 = 1. Other components are determined recursively:
\begin{align*} u_{n+1} (t) &= \int_0^t {\text d}\tau \int_0^{\tau} {\text d}s \left( \left[ 1 - \frac{3}{\cosh^2 s} \right] u_n (s) + A_n \left( u_0 , u_1 , \ldots u_n \right) \right) \\ &= \int_0^t {\text d}\tau \left( t- \tau \right) \left( \left[ 1 - \frac{3}{\cosh^2 \tau} \right] u_n (\tau ) + A_n \left( \tau \right) \right) . \end{align*}
   ■

Example: Consider the initial value problem for the nonlinear differential equation

\[ \ddot{y} - e^{-4t} y^3 = 3\, e^{2t} , \qquad y(0) = 0, \quad \dot{y}(0) = 2. \]
The above problem has the exact solution y(t) = e2t. With the initial term
\[ u_{n+1} (t) = 3\, e^{2t} -1-2t - \int_0^t {\text d}\tau \int_0^{\tau} {\text d}s e^{-4s}\.y^3 (s) , \qquad n=0,1,2,\ldots . \]
We expand the nonlinear term into the Adomian series
\[ u_{n+1} (t) = 3\, e^{2t} -1-2t - \int_0^t {\text d}\tau \int_0^{\tau} {\text d}s e^{-4s}\.A_n (s) , \qquad n=0,1,2,\ldots . \]
   ■

Example: Consider the initial value problem for the Lane--Emden equation

\[ \frac{{\text d}^2 y}{{\text d}t^2} + \frac{2}{t}\,\frac{{\text d}y}{{\text d}t} + y^5 =0 , \qquad y(0) = 1. \]
This second order differential equation is named after the American astrophysicist Jonathan Homer Lane (1819--1880) and the Swiss astrophysicist and meteorologist Jacob Robert Emden (1862--1940). It is a so called the singular differential equation having a critical point at t = 0. We don't know a uniqueness and/or existence theorem for the formulated above initial value problem. The solutions of the Lane--Emden equation that are finite at the origin have necessarily dy/dt = 0 at t = 0. The Lane--Emden operator admits a factorization:
\[ t^{-2} \frac{\text d}{{\text d}t} \left( t^2 \frac{{\text d}y}{{\text d}t} \right) + y^n =0 \qquad\mbox{or} \qquad t^{-1} \frac{{\text d}^2}{{\text d}t^2} \left( t\,y \right) . \]
The Lane--Emden equation is known to have an explicit solution for n = 0, 1, 5. Since for his equation we don't have uniqueness and existence theorems, there could be simgular solutions. In paticular, an explicit solution for n = 5 is
\[ y = \psi (t) = \left( 1 + \frac{t^2}{3} \right)^{-1/2} = 1 - \frac{t^2}{6} + \frac{t^4}{24} - \frac{5\,t^6}{432} + \frac{35}{10 368}\,t^8 - \frac{7\, t^{10}}{6912} + \cdots . \]
The Adomian polynomials for the quintic nonlinearity y5 are \begin{align*} A_0 &= u_0^5 , \\ A_1 &= 5\,u_0^4 u_1 , \\ A_2 &= 10\,u_0^3 u_1^2 + 5\,u_0^4 u_2 , \\ A_3 &= 10\,u_0^2 u_1^3 + 20\, u_0^3 u_1 u_2 + 5\,u_0^4 u_3 , \\ A_4 &= 5\,u_0 u_1^4 + 30\,u_0^2 u_1^2 u_2 + 10\, u_0^3 u_2^2 + 20\, u_0^3 u_1 u_3 + 5\,u_0^4 u_4 , \\ A_5 &= u_1^5 + 20\,u_0 u_1^3 u_2 + 30\,u_0^2 u_1 u_2^2 + 30\, u_0^2 u_1^2 u_3 + 20\,u_0^3 u_2 u_3 + 20\,u_0^3 u_1 u_4 + 5\,u_0^4 u_5 . \end{align*} The components of the decomposition solution \( \displaystyle y(t) = 1 + \sum_{n\ge 1} u_n \) are determined by the recursive scheme consisting of initial value problems:
\[ L \left[ u_{n+1} \right] = - A_n , \qquad u_{n+1} (0) = u'_{b+1} (0) = 0 , \qquad n=0,1,2,, \ldots . \]
where
\[ L \left[ y(t) \right] = t^{-2} \frac{\text d}{{\text d}t} \left( t^2 \frac{{\text d}y}{{\text d}t} \right) \qquad \Longrightarrow \qquad L^{-1} \left[ f \right] = y(0) + \int_0^t t^{-2} {\text d}t \int_0^t s^2 f(s)\,{\text d}s = y(0) + \int_0^t \frac{s}{t} \left( t-s \right) f(s) \,{\text d}s . \]
Computations show that \begin{align*} u_1 (t) &= - L \left[ 1 \right]^{-1} = - \frac{t^2}{6} , \\ u_2 (t) &= - L \left[ - \frac{5}{6}\, t^2 \right]^{-1} = \frac{t^4}{24} , \\ u_3 (t) &= - L \left[ 10\,u_1^2 + 5\,u_2 \right]^{-1} = - \frac{5\,t^6}{432} , \\ u_4 (t) &= - L \left[ 10\,u_1^3 + 20\,u_1 u_2 + 5\,u_3 \right]^{-1} = \frac{35\,t^8}{10 368} , \\ u_5 (t) &= - L \left[ 5\,u_1^4 + 30\,u_1^2 u_2 + 10\,u_2^2 + 20\,u_1 u_3 + 5\,u_4 \right]^{-1} = - \frac{7\, t^{10}}{6912} . \end{align*}
Integrate[s/t*(t - s)*s^2/6, {s, 0, t}]*5
t^4/24
Integrate[s/t*(t - s)*(10*s^4/6^2 + 5*s^4 /24), {s,0,t}]
(5 t^6)/432
Integrate[s/t*(t - s)*(10*s^6/6^3 + 20*s^2/6*s^4/24+ 5*5*s^6 /432), {s,0,t}]
(35 t^8)/10368
Finally, we plot 10-degree approximation
\[ \phi_{5} (t) = 1 + u_1 + u_2 + u_3 + u_4 + u_5 = 1 - \frac{t^2}{6} + \frac{t^4}{24} - \frac{5\,t^6}{432} + \frac{35}{10 368}\,t^8 - \frac{7\, t^{10}}{6912} \]
along with the explicit solution
     
phi5[t_]= 1- t^2/6 + t^4/24 -(5 t^6)/432 + (35 t^8)/10368 - 7* t^(10)/6912;
y[t_]= (1+t^2 /3)^(-1/2);
Plot[{phi5[t], y[t]}, {t, -3, 3}, PlotRange -> {-1.2, 1.1}, PlotStyle -> Thick, PlotLabels -> {"MDM-approximation", "true solution"}]
       MDM approximation and the true solution.            Mathematica code

   ■

Example: In astrophysics, the Emden–Chandrasekhar equation is a dimensionless form of the Poisson equation for the density distribution of a spherically symmetric isothermal gas sphere subjected to its own gravitational force, named after Robert Emden (1862--1940) and Subrahmanyan Chandrasekhar (1910--1995). The equation was first introduced by Robert Emden in 1907. The equation reads

\[ \frac{1}{\xi^2} \,\frac{\text d}{{\text d}\xi} \left(\xi^2 \frac{{\text d}\psi}{{\text d}\xi} \right) = e^{-\psi} , \]
where ξ is the dimensionless radius and ψ is the related to the density of the gas sphere as \( \displaystyle \rho = \rho_c e^{-\psi} , \) where ρ is the density of the gas at the centre. This differential equation is usually subject to the homogeneous initial conditions. Moreover, the corresponding initial value problem has a singular solution
\[ e^{-\psi_s} = 2x^2 \qquad\mbox{or} \qquad \psi_s = -2\,\ln x - \ln 2. \]
If ψ(ξ) is a solution to Emden–Chandrasekhar equation, then ψ(Cξ) - 2 ln(C) is also a solution of the equation, where C is an arbitrary constant. The solutions of the Emden–Chandrasekhar equation which are finite at the origin have necessarily dψ/dξ = 0 at ξ = 0.

The Lane--Emden operator admits a factorization:

\[ L\left[ \psi \right] = \xi^{-2} \frac{\text d}{{\text d}\xi} \left( \xi^2 \frac{{\text d}\psi}{{\text d}\xi} \right) \qquad \Longrightarrow \qquad L^{-1} \left[ f \right] = \int_0^{\xi} t^{-2} {\text d}t \int_0^t s^2 f(s)\,{\text d}s = \int_0^{\xi} \frac{s}{\xi} \left( \xi -s \right) f(s) \,{\text d}s . \]
We decompose the solution and exponential nonlinarity:
\[ \psi (\xi ) = \sum_{n\ge 1} u_n \qquad (u_0 =0) \]
and
\[ e^{-\psi} = \sum_{n\ge 0} A_n . \]
The Adomian polynomials are \begin{align*} A_0 &= e^0 = 1 , \\ A_1 &= - u_1 , \\ A_2 &= \frac{1}{2}\, u_1^2 - u_2 , \\ A_3 &= -\frac{1}{6}\,u_1^3 + u_1 u_2 - u_3 , \\ A_4 &= \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 , \\ A_5 &= -\frac{1}{120}\, 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 . \end{align*}
f[y_] = Exp[-y]
Ado[5]
where
Ado[M_] := Module[{c, n, k, j, der}, Table[c[n, k], {n, 1, M}, {k, 1, n}]; A[0] = f[u[0]];
For[n = 1, n <= M, n++, c[n, 1] = u[n];
For[k = 2, k <= n, k++, c[n, k] = Expand[1/n* Sum[(j + 1)*u[j + 1]*c[n - 1 - j, k - 1], {j, 0, n - k}]]];
der = Table[D[f[u[0]], {u[0], k}], {k, 1, M}];
A[n] = Take[der, n].Table[c[n, k], {k, 1, n}]]; Table[A[n], {n, 0, M}]]
Since the initial component u0 = 0, the other terms are solutions of the initial value problems:
\[ u_{n+1} = L^{-1} \left[ A_n \right] , \qquad n=0,1,2,\ldots . \]
Hence, \begin{align*} u_1 &= L^{-1} \left[ 1 \right] = \frac{\xi^2}{6} , \\ u_2 &= L^{-1} \left[ -u_1 \right] = - \frac{\xi^4}{120} , \\ u_3 &= L^{-1} \left[ \frac{1}{2}\, u_1^2 - u_2 \right] = \frac{\xi^6}{1890} , \\ u_4 &= L^{-1} \left[ - \frac{1}{6}\, u_1^2 - u_2 \right] = - \frac{37}{1632960} \, \xi^8 \\ u_5 &= L^{-1} \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] = \frac{431}{898128000}\, \xi^{10} . \end{align*}
u[1][t_] = Integrate[s/t*(t - s)*1, {s, 0, t}]
t^2/6
u[2][t_] = -Integrate[s/t*(t - s)*s^2 /6, {s, 0, t}]
-(t^4/120)
u[3][t_] = Integrate[s/t*(t - s)*(u[1][s]^2 /2 - u[2][s]), {s, 0, t}]
(t^6/1890)
u[4][t_] = Integrate[s/t*(t - s)*(-u[1][s]^3 /6 + u[1][s]* u[2][s] - u[3][s]), {s, 0, t}]
-((37 t^8)/1632960)
u[5][t_] = Integrate[s/t*(t - s)*(-u[1][s]^4 /24 - u[1][s]^2 * u[2][s]/2 + u[2][s]^2 /2 + u[1][s]*u[3][s] -u[4][s]), {s, 0, t}]
(431 t^10)/898128000
This allows us to find the five-term approximation
\[ \psi_5 (\xi ) = \frac{\xi^2}{6} - \frac{\xi^4}{120} + \frac{\xi^6}{1890} - \frac{37}{1632960} \, \xi^8 + \frac{431}{898128000}\, \xi^{10} . \]
phi5[x_]= x^2 /6 - x^4/120 + x^6 /1890 -37*x^8 /1632960 + 431*x^(10) /898128000;
Plot[phi5[x], {x,0,2}, PlotStyle->Thick]
Fifth-term MDM approximation to the solution
   ■

Example: Let us consider the pendulum equation

      \[ \ddot{\theta} + \omega^2 \sin \theta = 0, \qquad \theta (0) = \theta_0 = \frac{\pi}{3} , \quad \dot{\theta} (0) = 0. \]
a = {Graphics[{Arrowheads[{0.0, 0.07}], Arrow[{{0, -0.5}, {0, 0.5}}]}]};
dash = ParametricPlot[{2.02*Cos[t], 2.02*Sin[t]}, {t, -2.5, -0.5}, PlotStyle -> Dashed, Axes -> False]; vline = Graphics[Line[{{0, 0}, {0, -2.1}}]];
angle = ParametricPlot[#[[1]]*{Cos[\[Theta]], Sin[\[Theta]]}, {\[Theta], #[[ 2]], #[[3]]}, Axes -> False, PlotStyle -> #[[4]]] /. Line[x_] :> Sequence[Arrowheads[{-0.08, 0.08}], Arrow[x]] & /@ {{-1.5, 90 Degree, 130 Degree, Cyan}};
line = Graphics[{Thickness[0.01], Line[{{0, 0}, {1.25, -1.45}}]}]; disk = Graphics[{Gray, Disk[{1.33, -1.53}, 0.1]}];
sin = Graphics[Text[ Style["A sin(\[Omega]t)", FontSize -> 16, Black], {-0.46, 0.0}]];
theta = Graphics[Text[ Style["\[Theta](t)", FontSize -> 16, Black], {0.4, -1.0}]]
; m = Graphics[Text[ Style["m", FontSize -> 16, Black], {1.55, -1.53}]];
Show[dash, a, vline, angle, line, disk, m, theta, PlotRange -> All]
       Pendulum.            Mathematica code

Here θ is the angle between the rod of length ℓ and the downward vertical direction with the counterclockwise direction taken as positive. Also, ω² = g/ℓ, with g being the acceleration due to gravity; as usual, derivatives with respect to time variable t are denoted by dots: \( \displaystyle \dot{\theta} = {\text d}\theta / {\text d}t \quad\mbox{and}\quad \ddot{\theta} = {\text d}^2 \theta / {\text d}t^2 . \)

The exact solution can be expressed in terms of a Jacobi elliptic function \[ \theta (t) = 2\,\arcsin \left\{ \sin\frac{\theta_0}{2} \, \mbox{sn} \left[ K\left( \sin^2 \frac{\theta_0}{2} \right) - \omega_0 t ; \sin^2 \frac{\theta_0}{2} \right] \right\} , \] where K(m) is the complete elliptic integral of the first kind and sn(x, k) is the elliptic sine function. Here θ0 = π/3 is the initial displacement. We seek the solution of the given initial value problem in the form of convergent power series \[ \theta (t) = c_0 + \sum_{n\ge 1} c_n t^n , \qquad \mbox{with} \quad c_0 = \theta_0 = \frac{\pi}{3} , \quad c_1 = 0 . \] Since the initial velocity is homogeneous, all Adomian polynomials of odd indices are zeroes. The Adomian polynomials of even indices in terms of the decomposition coefficients $c_n$ for the sinusoidal nonlinearity $f(y) = \sin y$ are \begin{align*} A_0 &= f\left( c_0 \right) = \sin \theta_0 = \sin \frac{\pi}{3} = \frac{\sqrt{3}}{2} , %\\ %A_1 &= c_1 \cos \theta_0 = c_1 \frac{1}{2} = 0, \\ A_2 &= c_2 \cos \theta_0 - \frac{c_1^2}{2}\,\sin \theta_0 = \frac{c_2}{2} , %\\ %A_3 &= c_3 \cos \theta_0 - c_1 c_2 \sin \theta_0 - \frac{c_1^3}{6}\, \cos %\theta_0 = \frac{c_3}{2} , \\ A_4 &= c_4 \cos \theta_0 - \frac{c_2^2}{2}\, \sin \theta_0 = \frac{c_4}{2} - \frac{c_2^2}{2}\,\frac{\sqrt{3}}{2} , \\ A_6 &= c_6 \cos \theta_0 - c_2 c_4 \sin\theta_0 , \\ A_8 &= \left( c_8 - \frac{1}{2}\, c_2^2 c_4 \right) \cos \theta_0 + \left( \frac{c_2^4}{4} - c_2 c_6 - \frac{c_4^2}{2} \right) \sin\theta_0 , \end{align*} and so on. Substituting these series into the pendulum equation, we get the recurrence: \[ \left( n+2 \right) \left( n+1 \right) c_{n+2} + \omega^2 A_n = 0 , \qquad n=0,1,2,\ldots . \] Solving it, we obtain the values of coefficients: \begin{align*} c_2 &= - \frac{\omega^2}{2}\, A_0 = - \frac{\omega^2}{2}\,\frac{\sqrt{3}}{2} \\ c_3 &= - \frac{\omega^2}{6}\, A_1 = 0, \\ c_4 &= - \frac{\omega^2}{12}\, A_2 = - \frac{\omega^2}{12}\, \frac{c_2}{2} = \frac{\omega^4}{24} \, ,\frac{\sqrt{3}}{2} , \\ c_5 &= - \frac{\omega^2}{20}\, A_3 = 0, \\ c_6 &= - \frac{\omega^2}{30}\, A_4 = - \frac{\omega^2}{30}\left( \frac{c_4}{2} - \frac{c_2^2}{2}\,\frac{\sqrt{3}}{2} \right) , \end{align*} and so on. This gives us a polynomial approximation \begin{align*} \theta (t) &= \theta_0 - \frac{\left( \omega\,t \right)^2}{2!}\,\sin\theta_0 + \frac{\left( \omega\,t \right)^4}{4!}\, \sin\theta_0 \cos\theta_0 - \frac{\left( \omega\,t \right)^6}{6!}\left( \sin\theta_0 \cos^2 \theta_0 - 3\,\sin^3 \theta_0 \right) \\ &\quad + \frac{\left( \omega\,t \right)^8}{8!}\left( \frac{19}{8}\,\sin (4\,\theta_0 ) - \frac{17}{4}\,\,\sin (2\,\theta_0 ) \right) - \frac{\left( \omega\,t \right)^{10}}{10!}\left( 410\,\sin (\theta_0 ) - \frac{1929}{8}\,\sin (3\theta_0 ) + \frac{503}{8}\,\sin (5\theta_0 ) \right) + \cdots . \end{align*}

     
s = NDSolve[{theta''[t] + theta[t] == 0, theta[0] == Pi/3, theta'[0] == 0}, theta, {t, 0, 30}];
pendulum = Plot[Evaluate[theta[t] /. s], {t, 0, 5}, PlotStyle -> {Thickness[0.01]}];
t0 = Pi/3; c2 = -omega^2*Sin[t0]/2;
A2 = c2*Cos[t0];
c4 = -omega^2 *A2/12;
A4 = c4*Cos[t0] - (c2)^2 *Sin[t0]/2;
c6 = -omega^2 *A4/30;
A6 = Simplify[c6*Cos[t0] - c2*c4*Sin[t0]];
c8 = -omega^2*A6/8/7;
p6[t_] = (t0 + c2*t^2 + c4*t^4 + c6*t^6) /. omega -> 1
p8[t_] = (t0 + c2*t^2 + c4*t^4 + c6*t^6 + c8*t^8) /. omega -> 1
pp6 = Plot[p6[t], {t, 0, 30}, PlotRange -> {{0, 5}, {-1.1, 1.1}}, PlotStyle -> {Purple, Thick}];
pp8 = Plot[p8[t], {t, 0, 30}, PlotRange -> {{0, 5}, {-1.1, 1.1}}, PlotStyle -> {Orange, Thick}];
txt = Graphics[{Black, Text[Style["True solution", Bold, 18], {1.3, -1}]}];
txt6 = Graphics[{Black, Text[Style["degree 6", Italic, 16], {3.9, 1}]}];
txt8 = Graphics[{Black, Text[Style["degree 8", Italic, 16], {4.5, -0.88}]}];
Show[pendulum, pp6, pp8, txt, txt6, txt8]
       True solution and two Adomian polynomial approximations.            Mathematica code

Here is Mathematica code from Duan & rach (2011) article:
f[u] := Sin[u]; \[Gamma]=25; C1=9.99; k=C1/2/Sqrt[\[Gamma]]; y=2*ArcSin[k*JacobiSN[Sqrt[\[Gamma]]*t, kˆ2]]; AP[n] := Module[{Z, der, i, j, k, m, s}, A[0]=f[v[0]]; Table[Z[i,j], i,1,n, {j,1,i}]; der=Table[Derivative[k][f][v[0]], {k,1,n}]; For[m=1, m<=n, m++, Z[m,1]=v[m]; For[k=2, k<=m, k++, If[Head[Z[m-1,k-1]]===Plus, Z[m,k]=Map[#*v[1]/(Exponent[#,v[1]]+1)&,Z[m-1,k-1]], Z[m,k]=Map[#*v[1]/(Exponent[#,v[1]]+1)&,Z[m-1,k-1],0]]]; For[k=2,k<=Floor[m/2], k++, Z[m,k]=Z[m,k]+(Z[m-k,k]/.v[s_] -> v[s+1])]; A[m]=Table[Z[m,k],k,1,m].Take[der,m];];]; AP[37]; phi[t_,t0_,C0_,C1_,n_]:=Module[{m},a[0]=C0; a[1]=C1; Do[a[m+2]= Expand[-\[Gamma]/(m+1)/(m+2)*(A[m]/.v->a)] ,m,0,n-3]; Sum[a[m]*(t-t0)ˆm, {m,0,n-1}]]; f1=Prepend[Table[phi[t,0,0,C1,10*2ˆk], {k,0,2}], y]; F1=Show[Plot[f1,{t,-1,1}, AxesLabel->{t}, PlotRange->{-4,4}]]; soln[t0_,T_,C0_,C1_,h_,n_]:=Module[{k, m, a}, N1=Floor[(T-t0)/h]; t[0]=t0; u[0]=C0; u1[0]=C1; For[k=0, k<=N1–1, k++, t[k+1]=t[k]+h; u[k+1]=phi[t[k+1], t[k], u[k], u1[k], n+1]; u1[k+1]=D[phi[t, t[k], u[k],u1[k], n+1], t]/.t->t[k+1];]; f2=Table[{t[k], u[k]}, {k, 0, N1}]; ListPlot[f2]]; ef=Plot[y, {t, 0, 10}, PlotRange->{-4, 4}, AxesLabel->{t}]; {F1, Show[ef,soln[0,10,0,9.99,0.1,9]],Show[ef,soln[0,10,0,9.99,0.2,20]]}
   ■

Example: We consider a version of Troesch’s problem

\[ \ddot{y} = n\,\sinh (n\,y) , \qquad y(0) =0, \quad \dot{y}(0) = \sqrt{2} . \]
Numerical solutions of this problem may show instability when the parameter n exceeds 5.    ■

 

  1. Duan, J.-S., Rach, R., New higher-order numerical one-step methods based on the Adomian and the modified decomposition methods, Applied Mathematics and Computation, 2011, Vol. 218, Issue 6, pp. 2810--2828. https://doi.org/10.1016/j.amc.2011.08.024
  2. Roberts, S.M., Shipman, J.S., On the closed form solution of Troesch’s problem, Journal of Computational Physics, 1976, Vol. 21,, Issue 3, pp. 291--304.

 

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)