# Preface

This section discusses a completely different strategy to solve an initial value problem for a single first-order differential equation $$y' = f(x,y) , \quad y(x_0 )= y_0 .$$ In order to calculate an approximate value of the solution φ(tn+1) at the next mesh point t = tn+1, the values of the calculated solution at some previous mesh points are used. The numerical methods that use information at more than the last mesh point are referred to as multistep methods. This section presents two types of multistep methods: Adams methods and backward differentiation methods. Different levels of accuracy can be achieved with each type of methods.

# Variational Iteration Method

Recently, there has been great development of new powerful methods capable of handling linear and nonlinear equations. Among them, there are known the Adomian decomposition method (ADM), the variational iteration method, and the homotopy perturbation method. The variational iteration method (VIM, for short) was developed by Ji-Huan He in 1999--2007. This method is preferable over numerical methods as it is free from rounding off errors as it does not involve discretization and does not require large computer power or memory. The VIM gives rapidly convergent successive approximations of the exact solution if such a solution exists; otherwise, a few approximations can be used for numerical purposes. Unlike the ADM, where computational algorithms are normally used to deal with nonlinear terms, the VIM does not require the utilization of restrictive assumptions because it approaches linear and nonlinear problems directly in a like manner. Many researches in variety of scientific fields applied this method and showed the VIM has many merits and to be reliable for a variety of scientific application. In particular, Abassy et al., proposed the modified variational iteration method (MVIM) to avoid repeated computation of redundant terms.

To illustrate the basic concept of the He’s VIM, consider the following general nonlinear equation:

$L\left[ y(x) \right] = N\left[ y(x) \right] + g(x) ,$
where $$L = \texttt{D} = {\text d}/{\text d}x$$ is the linear differential operator (which in general can include other terms), N is the nonlinear operator, and g is the input (known) function. He has modified the general Lagrange multiplier method to an iteration method called as correction functional. The basic character of the method is to construct a correction functional for the aforementioned equation, which reads as
$y_{n+1} (x) = y_n (x) + \int_0^x \lambda (s) \left\{ L\left[ y_n (s) \right] - N\left[ \tilde{y}_n (s) \right] - g(s) \right\} {\text d}s ,$
where λ is a Lagrange multiplier. The subscript n denotes the n-th approximation, and $$\tilde{y}_n$$ is a restricted variation., i.e., $$\delta\tilde{y}_n = 0 .$$ We start our exposition with the linear case:
$y' (x) + p(x)\, y = g(x) , \qquad y(0) = \alpha .$
According to the VIM, the basic character of the method is to construct a correction functional for the equation, which reads
$y_{n+1} (x) = y_n (x) + \int_0^x \lambda (s) \left\{ y'_n (s) + p(s)\, \tilde{y}_n (s) - g(s)\right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Here λ(s,x) is so called a general Lagrange multiplier, which can be identified optimally via variational theory, and $$\tilde{y}_n$$ denotes a restricted variation, that is, $$\delta\tilde{y}_n =0.$$ Having λ obtained, an iteration formula should be used for the determination of the successive approximations yn+1(t) of the solution y(t). The initial approximation y0(t) can be arbitrary function; however, it is usually selected based on the initial conditions. Consequently, the solution is given by
$y(x) = \lim_{n\to\infty}\,y_n (x) .$

The main problem for obtaining highly accurate solutions with the VIM is the determination of the Lagrange multiplier, and there are known many forms of them. This parameter plays the important role in the speed of convergence of the VIM solution. Finding the Lagrange multiplier could be achieved by making the correction functional stationary. Taking variation with respect to the variable yn, noticing δyn(0)=0, we get

\begin{eqnarray*} \delta y_{n+1} (x) &=& \delta y_n (x) + \delta \int_0^x \lambda \left\{ y'_n (s) + \tilde{y}_n (s) - g(s) \right\} {\text d}s \\ &=& \delta u_n (x) + \left. \lambda\,\delta y_n \right\vert_{s=x} + \int_0^x \left[ -\lambda ' \delta y_n \right] {\text d}s =0 , \end{eqnarray*}
for all variations δyn. This yields the following equation
$\left. \lambda ' \right\vert_{s=x} =0 \qquad\mbox{and} \qquad \left. 1+ \lambda \right\vert_{s=x} =0 .$
Solving the above equations, we obtain λ = -1, and the correction functional gives the iteration formula
$y_{n+1} (x) = y_n (x) - \int_0^x \left\{ y'_n (s) + p(s)\,y_n (s) - g(s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Originally, Ji-Huan He did not consider the linear part as a restricted variation and applied it only to nonlinear terms. Furthermore, the lesser the application of restricted variations the faster the approximations converging to its exact solution. In this case, he obtained (using integration by parts) the following equations for λ determination:
$\begin{split} -\lambda ' (s) + p(s)\,\lambda (s) &=0 , \\ 1 + \left. \lambda (s) \right\vert_{s=t} &= 0. \end{split}$
Solving the above equations with Mathematica,
DSolve[{lambda'[s] == p[s]*lambda[s], lambda[x] == -1}, lambda[s], s]
we obtain another formula for the Lagrange multiplier
$\lambda (s,x) = - \exp \left\{ \int_x^s p(\tau )\,{\text d} \tau \right\} .$
In both approaches, there are repeated calculations in each step (similar to Picard's method); to cancel some of the repeated calculations, the iteration formula can be handled as follows.
$y_{n+1} = y_n - \int_0^x y'_n (s)\, {\text d}s - \int_0^x \left\{ p(s)\,y_n (s) - g(s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Hence,
$y_{n+1} = y_n (0) - \int_0^x \left\{ p(s)\,y_n (s) - g(s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
We can set $$y_n (0) = y_0 = \alpha ,$$ so
$y_{n+1} = y_0 - \int_0^x \left\{ p(s)\,y_n (s) - g(s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
To eliminate all repeated computations, let us rewrite the latter in the following iteration formula:
$y_{n+1} = y_0 - \int_0^x \left\{ p(s)\left( y_n (s) - y_{n-1} (s) \right) \right\} {\text d}s - \int_0^x \left\{ p(s)\,y_{n-1} (s) -g(s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Using the identity $$y_{n} = y_0 - \int_0^x \left\{ p(s)\,y_{n-1} (s) - g(s) \right\} {\text d}s ,$$ we obtain so called the modified variational iteration formula:
$y_{n+1} = y_n - \int_0^x \left\{ p(s)\left( y_n (s) - y_{n-1} (s) \right) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
We proceed further iterations based on $$y_{0} = y(0) = \alpha , \ y_{-1} =0 .$$ Now it is right time to show how it works in a series on examples.

Example: Consider the IVP for the linear differential equation

$y' -y = \sin x , \qquad y(0) =1 .$
Mathematica easily find its solution:
DSolve[{y'[x] == y[x] + Sin[x], y[0] == 1}, y[x], x]// Flatten
$y(x) = \frac{1}{2} \left( 3\,e^x - \cos x - \sin x \right) = 1 + x + x^2 + \frac{x^3}{3} + \frac{x^4}{24} + \frac{x^5}{5!} + \frac{x^6}{360} + \cdots .$
For solving this initial value problem by MVIM, we set $$y_{-1} =0, \ y_0 =1 .$$ Then we use the following iterative formula
$y_{n+1} = y_n + \int_0^x \left\{ y_n (s) - y_{n-1} (s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
The first term needs a special attention.
$y_{1} (x) = y_0 + \int_0^x \left\{ y_0 (s) - y_{-1} (s) + \sin s \right\} {\text d}s = 2 +x - \cos x .$
Other terms are calculated according to the above formula.
\begin{eqnarray*} y_{2} (x) &=& y_1 + \int_0^x \left[ y_1 (s) - y_0 (s) \right] {\text d}s = \frac{1}{2} \left( 2+x \right)^2 - \cos x - \sin x , \\ y_{3} (x) &=& y_2 + \int_0^x \left[ y_2 (s) - y_1 (s) \right] {\text d}s = \left( 1+x \right)^2 + \frac{x^3}{6} - \sin x = 1 + x + x^2 + \frac{x^3}{3} + \frac{x^4}{24} - \frac{x^5}{5!} + \cdots , \\ y_{4} (x) &=& y_3 + \int_0^x \left[ y_3 (s) - y_2 (s) \right] {\text d}s = 1 + x + x^2 + \frac{x^3}{3} + \frac{x^4}{24} , \end{eqnarray*}
and so on.
y1[x_] = 1 + Integrate[1 + Sin[s], {s, 0, x}];
y2[x_] = y1[x] + Integrate[y1[s] - 1, {s, 0, x}];
y3[x_] = Simplify[y2[x] + Integrate[y2[s] - y1[s], {s, 0, x}]];
y4[x_] = Simplify[y3[x] + Integrate[y3[s] - y2[s], {s, 0, x}]];
We plot last two approximations along with the true solution.
exact[x_] = (3*Exp[x]-Cos[x]-Sin[x])/2 ;
Plot[{Callout[y4[x], "four-term VIM approximation"], Callout[y3[x], "three-term VIM approximation"], Callout[exact[x], "exact"]}, {x, 0, 3.5}, PlotStyle->Thick]
■

Example: Consider the IVP for the generalized logistic equation

$\dot{u} = u \left( 2-u \right) \left( 4-u \right) , \qquad u(0) =1.$
Here dot represents differentiation with respect to time variable. First, find its solution:
exact[t_] = u[t]/.DSolve[{u'[t] == u[t]*(2 - u[t])*(4 - u[t]), u[0] == 1}, u[t], t]
$u(t) = 2 \,\frac{1 + 3\, e^{8t} - \sqrt{1 + 3\, e^{8t}}}{1 + 3\, e^{8t}} = 1 + 3t - \frac{3}{2}\, t^2 - \frac{17}{2}\, t^3 + \frac{125}{8}\, t^4 + \frac{721}{40}\, t^5 - \frac{23401}{240}\, t^6 + \cdots .$
In the hybrid model, u(t) is the portion of population of a certain characteristic, and t is the time measured in generations.

Since the Lagrange multiplier is λ = -1, we have the following recurrence:

$u_{n+1} = u_n - \int_0^x \left\{ u'_n (s) - u_n (s) \left( 2-u_n (s) \right) \left( 4-u_n (s) \right) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Since the given differential equation is homogeneous, we don't need to distinguish u1 from other terms. Settingu0 to be the initial value, we have
\begin{eqnarray*} u_{0} (t) &=& 1, \\ u_1 (t) &=& u_0 + \int_0^t \left\{ u_0 (s) \left( 2-u_0 (s) \right) \left( 4-u_0 (s) \right) \right\} {\text d}s = 1+ 3t , \\ u_2 (t) &=& u_1 - \int_0^t \left\{ u'_1 (s) - u_1 (s) \left( 2-u_1 (s) \right) \left( 4-u_1 (s) \right) \right\} {\text d}s = 1+ 3t - \frac{3}{2}\, t^2 - 9\,t^3 + \frac{27}{4} \, t^4 , , \\ u_3 (t) &=& u_2 - \int_0^t \left\{ u'_2 (s) - u_2 (s) \left( 2-u_2 (s) \right) \left( 4-u_2 (s) \right) \right\} {\text d}s = 1+ 3t - \frac{3}{2}\, t^2 - 9\,t^3 + \frac{27}{4} \, t^4 , , \\ u_4 (t) &=& u_3 - \int_0^t \left\{ u'_3 (s) - u_3 (s) \left( 2-u_3 (s) \right) \left( 4-u_3 (s) \right) \right\} {\text d}s = 1+ 3t - \frac{3}{2}\, t^2 - 9\,t^3 + \frac{27}{4} \, t^4 , , \end{eqnarray*}
and so on.
u1[t_] = 1 + 3*t;
u2[t_] = u1[t] - Integrate[D[u1[s], s] - u1[s]*(2 - u1[s])*(4 - u1[s]), {s, 0, t}];
u3[t_] = u2[t] - Integrate[D[u2[s], s] - u2[s]*(2 - u2[s])*(4 - u2[s]), {s, 0, t}];
u4[t_] = u3[t] - Integrate[D[u3[s], s] - u3[s]*(2 - u3[s])*(4 - u3[s]), {s, 0, t}]
We plot last two approximations along with the true solution.
Plot[{Callout[u4[t], "four-term VIM approximation"], Callout[u3[t], "three-term VIM approximation"], Callout[exact[t], "exact solution"]}, {t, 0, 1.4}, PlotStyle->Thick]
■

Example: Consider the IVP for the following Riccati equation

$y' = 3\,y - y^2 + 2\,x, \qquad y(0) =0.$
First, we write the variational iteration sequence:
$y_{n+1} (x) = y_n (x) + \int_0^x \lambda (x,s) \left\{ y'_n (s) - 3\,y_n (s) + \tilde{y}_n^2 - 2\,s \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Taking variation, we obtain
\begin{eqnarray*} \delta y_{n+1} (x) &=& \delta y_n (x) + \delta \int_0^x \lambda \left\{ y'_n (s) - 3\,y_n (s) + \tilde{y}_n^2 - 2\,s \right\} {\text d}s \\ &=& \delta y_n (x) + \delta \int_0^x \lambda \left\{ y'_n (s) - 3\,y_n (s) \right\} {\text d}s \\ &=& \delta y_n (x) + \left. \lambda\,\delta y_n (s) \right\vert_{s=0}^{s=t} - \int_0^t \left[ \lambda ' + 3\,\lambda \right] \delta y_n (s) \, {\text d}s . \end{eqnarray*}
This gives us the equations for determination of the Lagrange multiplier:
$\begin{split} \lambda ' (s) + 3\,\lambda (s) &=0 , \\ 1 + \left. \lambda (s) \right\vert_{s=t} &= 0. \end{split}$
So the multiplier can be identified as
$\lambda (x,s) = - e^{3(x-s)} .$
Then the iteration formula becomes
$y_{n+1} = y_n - \int_0^x e^{3(x-s)} \left\{ y'_n (s) - 3\,y_n (s) + y_n^2 - 2\,s \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Using previously defined subroutine seriesDSolve, we find series solution of the given problem:
seriesDSolve[ y'[x] == 2*x + 3*y[x] - (y[x])^2, y, {x, 0, 10}, {y[0] -> 0}]
$y = x^2 + x^3 + \frac{3}{4}\,x^4 + \frac{x^5}{4} - \frac{5}{24}\, x^6 - \frac{25}{56}\,x^7 - \frac{187}{448}\, x^8 - \frac{2551}{12096}\, x^9 + \frac{1217}{40320}\, x^{10} + \cdots .$
\begin{eqnarray*} y_1 (x) &=& x^2 - \int_0^x e^{3(x-s)} \left\{ 2s - 3\,s^2 + s^4 - 2\,s \right\} {\text d}s = x^2 + \frac{1}{81} \left( 27\,x^4 -45\,x^2 -30\,x -10 \right) + \frac{10}{81}\, e^{3x} , \\ y_2 (x) &=& y_1 - \int_0^x e^{3(x-s)} \left\{ y'_1 -3\,y_1 + y_1^2 -2\,s \right\} {\text d}s \approx x^2 - \frac{x^3}{3} - \frac{x^4}{4} - \frac{3}{20}\, x^5 - \frac{3}{40}\, x^6 + \frac{253}{840}\, x^7 + \cdots . \end{eqnarray*}
So we see that calculations become messy very shortly. We ask Mathematica for help.
y1[t_] = t^2 - Integrate[(s^4 - 3*s^2)*Exp[3*t - 3*s], {s, 0, t}];
y2[t_] = y1[t] - Integrate[(D[y1[s], s] - 3*y1[s] + (y1[s])^2 - 2*s)* Exp[3*t - 3*s], {s, 0, t}];
y3[t_] = y2[t] - Integrate[(D[y2[s], s] - 3*y2[s] + (y2[s])^2 - 2*s)* Exp[3*t - 3*s], {s, 0, t}];
Expand[y3[t]]
Series[%, {t, 0, 10}]
$y_3 (x) \approx x^2 + x^3 + \frac{3}{4}\,x^4 + \frac{x^5}{4} - \frac{5}{24}\, t^6 - \frac{25}{56}\, x^7 - \frac{187}{448}\, x^8 + \frac{2551}{12096} \, x^9 + \frac{1217}{40320}\, x^{10} + \cdots ,$
which gives the correct answer up to the power 10 (at least). Then we plot the third approximation y3 and the true solution.
sol = (y[x] /.
DSolve[{y'[x] == 2*x + 3*y[x] - (y[x])^2, y[0] == 0}, y[x], x])[[1]];
Plot[{Callout[sol, "exact"], Callout[y3[x], "three-term VIM approximation"]}, {x,0,1.4}, PlotStyle->Thick]

We demonstrate a modified variational iteration method and choose the Lagrange multiplier as λ = -1. This leads in our case to the recurrence:

$y_{n+1} (t) = y_n (t) + 3 \int_0^t \left( y_n (s) - y_{n-1} (s) \right) {\text d}s - \int_0^t \left( y_n^2 (s) - y_{n-1}^2 (s) \right) {\text d}s , \qquad n=1,2,\ldots ;$
where $$y_{-1} = 0, \ y_0 = 0$$ (from the initial condition). The first term should include the input function
$y_{1} (t) = y_0 + \int_0^t \left[ y'_0 (s) - 3 \, y_0 + y_0^2 - 2\,s \right] {\text d}s = t^2 .$
Therefore, we can obtain the following successive approximations:
\begin{eqnarray*} y_2 (t) &=& t^2 + 3\, \int_0^t \left[ y_1 (s) - y_{0} \right] {\text d}s - \int_0^t \left[ y_1^2 - y_{0}^2 \right] {\text d}s = t^2 + t^3 - \frac{t^5}{5} , \\ y_3 (t) &=& y_2 (t) + 3\, \int_0^t \left[ y_2 (s) - y_{1} (s) \right] {\text d}s - \int_0^t \left[ y_2^2 (s) - y_{1}^2 (s) \right] {\text d}s = t^2 + t^3 + \frac{3}{4}\,t^4 - \frac{t^5}{5} - \frac{13}{30}\, t^6 - \frac{t^7}{7} + \frac{t^8}{20} + \frac{2}{45}\,t^9 - \frac{t^{11}}{275} , \\ y_4 (t) &=& y_3 (t) + 3\, \int_0^t \left[ y_3 (s) - y_{2} (s) \right] {\text d}s - \int_0^t \left[ y_3^2 (s) - y_{2}^2 (s) \right] {\text d}s = t^2 + t^3 + \frac{3}{4}\,t^4 - \frac{t^5}{5} - \frac{13}{30}\, t^6 - \frac{5}{14}\, t^7 - \frac{11\,t^8}{80} + \frac{169}{2160}\,t^9 + \frac{61}{420} \, t^{10} + \cdots . \end{eqnarray*}
y2[t_] = t^2 + t^3 - t^5 /5
y3[t_] = y2[t] + Integrate[3*y2[s] - 3*s^2 - (y2[s])^2 + s^4, {s, 0, t}]
y4[t_] = y3[t] + Integrate[3*y3[s] - 3*y2[s] - (y3[s])^2 + (y2[s])^2, {s, 0, t}]
This yields
$y_{4} (t) = t^2 + t^3 + \frac{3}{4}\,t^4 + \frac{t^5}{4} - \frac{13}{30}\, t^6 - \frac{19}{35}\, t^7 - \frac{107}{560}\, t^8 + \frac{41}{432} \, t^9 + \frac{111}{700}\, t^{10} + \cdots$
Then we plot our outputs.
sol = (y[x] /.
DSolve[{y'[x] == 2*x + 3*y[x] - (y[x])^2, y[0] == 0}, y[x], x])[[1]];
Plot[{Callout[sol, "exact solution"], Callout[y3[x], "three-term MVIM approximation with $Lambda]=-1"], Callout[y4[x], "four-term MVIM approximation with \[Lambda]=-1"]}, {x, 0, 2.0}, PlotStyle -> Thick] ■ Example: Consider the initial value problem for the autonomous equation \[ y' = -\sqrt{1/4 + y^2} , \qquad y(0) =3.$
Its solution exists for all x and is given by
$y(x) = -\frac{1}{2}\,\sinh \left( x- \mbox{arcsinh}(6) \right) .$
This problem was considered in the previous sections where it was shown that Picard's iteration failed to provide the solution while the Adomian decomposition method was successful.

Using the initial approximation y0 = 3, we find the next VIM approximation:

$y_{1} (x) = u_0 - \int_0^x \sqrt{1/4 + y_0^2} {\text d}s = 3 - \frac{\sqrt{37}}{2} \,x$
because the Lagrange multiplier had been chosen to be λ=-1. We find all next approximants according to the MVIM formula:
$y_{n+1} (x) = u_n - \int_0^x \left[ \sqrt{1/4 + y_n^2} - \sqrt{1/4 + y_{n-1}^2} \right] {\text d}s , \qquad n=1,2,3,\ldots.$
Therefore,
\begin{eqnarray*} y_2 (x) &=& y_1 (x) - \int_0^x \left[ \sqrt{1/4 + y_1^2} - \sqrt{1/4 + y_{0}^2} \right] {\text d}s = 3 - \frac{\sqrt{37}}{2} \,x + \cdots , \\ y_3 (x) &=& y_2 (x) - \int_0^x \left[ \sqrt{1/4 + y_2^2} - \sqrt{1/4 + y_{1}^2} \right] {\text d}s . \end{eqnarray*}
y0[x_] = 3;
y1[x_] = 3 - Sqrt[37]*x/2;
y2[x_] = Simplify[ y1[x] - Integrate[Sqrt[1/4 + y1[s]] - Sqrt[1/4 + y0[s]], {s, 0, x}, Assumptions -> x > 0]];
y3[x_] = Simplify[ y1[x] - Integrate[Sqrt[1/4 + y2[s]] - Sqrt[1/4 + y1[s]], {s, 0, x}, Assumptions -> x > 0]]
and Mathematica fails to obtain the explicit formula. Therefore, we can obtain only the second approximation, and the VIM becomes unreliable.
exact[x_] = -Sinh[x - ArcSinh[6]]/2;
Plot[{exact[x], y2[x]}, {x, 0, 1.5}, PlotLegends -> {"exact solution", "second VIM approximation"}, PlotStyle -> Thick]
■

Example: Consider the IVP for the logistic equation

$u' = u^2 -2\,x\,u + 3\,x^2 , \qquad u(0) =2.$
While Mathematica is able to find its explicit solution
DSolve[{u'[x] == (u[x])^2 - 2*x*u[x] + 3*x^2, u[0] == 2}, u[x], x]
its representation is so messy that we prefer to deal with its approximation. First, we find its power series representation:
seriesDSolve[ u'[x] == (u[x])^2 - 2*x*u[x] + 3*x^2, u, {x, 0, 15}, {u[0] -> 2}]
$u(x) = 2 + 4\,x + 6\,x^2 + \frac{35}{3}\,x^3 + \frac{62}{3}\, x^4 + \frac{566}{15}\,x^5 + \frac{3112}{45}\,x^6 + \frac{1901}{15}\,x^7 + \frac{2089}{9}\,x^8 + \cdots .$
Using λ = -1, we have the variational iteration formula:
$u_{n+1} = u_n - \int_0^x \left[ u'_n (s) - u_n^2 (s) + 2\,s\,u_n (s) - 3\,s^2 \right] {\text d}s , \qquad n=0,1,2,\ldots .$
With $$u_0 =2, \ u_{-1} =0 ,$$ this allows us to find the next approximation:
$u_{1} = 2 + \int_0^x \left[ u_0^2 (s) - 2\,s\,u_0 (s) + 3\,s^2 \right] {\text d}s = 2 + 4\,x - 2\,x^2 + x^3 .$
Then we use the MVIM:
$u_{n+1} = u_n + \int_0^x \left[ u_n^2 (s) - u_{n-1}^2 (s) \right] {\text d}s - 2 \int_0^x s \left[ u_n (s) - u_{n-1} (s) \right] {\text d}s , \qquad n=1,2,3, \ldots .$
This yields
\begin{eqnarray*} u_2 (x) &=& u_1 (x) + \int_0^x \left[ u_1^2 (s) - u_{0}^2 \right] {\text d}s - \int_0^x s \left[ u_1 (s) - u_{0} (s) \right] {\text d}s = 2 + 4\,x + 6\,x^2 + \frac{7}{3}\, x^3 - \frac{5}{2}\, x^4 + \frac{11}{5}\, x^5 - \frac{2}{3}\,x^6 + \frac{x^7}{7} , \\ y_3 (t) &=& y_2 (t) + \int_0^x \left[ u_2^2 (s) - u_{1}^2 \right] {\text d}s - \int_0^x s \left[ u_2 (x) - u_{1} (s) \right] {\text d}s = 2 + 4\,x + 6\,x^2 + 13\,x^3 + \frac{77}[6}\,x^4 + \cdots , \\ y_4 (t) &=& y_3 (t) + \int_0^x \left[ u_3^2 (s) - u_{2}^2 \right] {\text d}s - \int_0^x s \left[ u_3 (x) - u_{2} (s) \right] {\text d}s = 2 + 4\,x + 6\,x^2 + 13\,x^3 + \frac{77}[6}\,x^4 + \frac{127}{15}\, x^5 + \frac{193}{60}\,x^6 + \cdots . \end{eqnarray*}
u0[x_] = 2;
u1[x_] = 2 + 4 x - 2 x^2 + x^3;
u2[x_] = u1[x] + Integrate[(u1[s])^2 - (2)^2 , {s, 0, x}] - Integrate[s*(u1[s] - 2), {s, 0, x}];
u3[x_] = u2[x] + Integrate[(u2[s])^2 - (u1[s])^2 , {s, 0, x}] - Integrate[s*(u2[s] - u1[s]), {s, 0, x}];
u4[x_] = u3[x] + Integrate[(u3[s])^2 - (u2[s])^2 , {s, 0, x}] - Integrate[s*(u3[s] - u2[s]), {s, 0, x}]
Finally, we plot the true solution along with its two VIM approximations.
sol = (y[x] /.
DSolve[{y'[x] == (y[x])^2 -2*x*y[x] + 3*x^2 , y[0] == 2}, y[x], x])[[1]];
Plot[{Callout[sol, "exact solution"], Callout[u3[x], "three-term MVIM approximation with $Lambda]=-1"], Callout[u4[x], "four-term MVIM approximation with \[Lambda]=-1"]}, {x, 0, 2.0}, PlotStyle -> Thick] Example: Consider the IVP for the variable coefficient Riccati equation \[ y' + e^x \,y + y^{-x}\, y^3 = \cosh x , \qquad y(0) = 1.$
Following the VIM, we obtain the recurrence:
$y_{n+1} = y_n - \int_0^x \left[ y'_n (s) + e^s y_n (s) + e^{-s} y_n^3 (s) - \cosh s \right] {\text d}s , \qquad n=0,1,2,\ldots .$
This in turn gives the successive approximations:
\begin{eqnarray*} y_0 (x) &=& 1 , \\ y_1 (x) &=& 1 - \int_0^x \left[ y'_0 (s) + e^s y_0 (s) + e^{-s} y_0^3 (s) - \cosh s \right] {\text d}s = 1- e^x + \cosh x , \\ y_2 &=& y_1 - \int_0^x \left[ y'_1 (s) + e^s y_1 (s) + e^{-s} y_1^3 (s) - \cosh s \right] {\text d}s = \frac{1}{32} \left( 67 + e^{-4x] + 8\,e^{-3} + 18\,e^{-2x} -32\,e^{-x} -40\,e^x + 10\,e^{2x} + 20\,x \right) \\ &=& 1 -x + 2\,x^2 - \frac{11}[6}\,x^3 + \frac{5}{3}\, x^4 - \frac{101}{120}\,x^5 + \frac{91}{180}\,x^6 - \cdots , \\ y_3 &=& y_2 - \int_0^x \left[ y'_2 (s) + e^s y_2 (s) + e^{-s} y_2^3 (s) - \cosh s \right] {\text d}s = 1 -x + 2\,x^2 - \frac{9}{2}\, x^3 + \frac{22}{3}\, x^4 - \frac{1397}{120}\,x^5 + \frac{967}{60}\,x^6 - \frac{103987}{5040}\, x^7 + \cdots . \end{eqnarray*}
y1[x_] = 1 - Integrate[Exp[s] + Exp[-s] - Cosh[s], {s, 0, x}]
y2[x_] = Simplify[ y1[x] - Integrate[ D[y1[s], s] + y1[s]*Exp[s] + Exp[-s]*(y1[s])^3 - Cosh[s], {s, 0, x}]]
y3[x_] = Simplify[ y2[x] - Integrate[ D[y2[s], s] + y2[s]*Exp[s] + Exp[-s]*(y2[s])^3 - Cosh[s], {s, 0, x}]]
Next we demonstrate the MVIM, which leads to the recurrence:
$u_{n+1} = u_n - \int_0^x \left[ e^s \left( u_n (s) - u_{n-1} (s) \right) + e^{-s} \left( u_n^3 (s) - u_{n-1}^3 (s) \right) \right] {\text d}s , \qquad n=1,2,\ldots .$
With $$u_{-1} =0 , \ u_0 =1 ,$$ we obtain
\begin{eqnarray*} u_1 (x) &=& 1 - \int_0^x \left[ e^s u_0 (s) + e^{-s} u_0^3 (s) - \cosh s \right] {\text d}s = 1- e^x + \cosh x =1-x - \frac{x^3}{6} - \frac{x^5}{120} - \cdots , \\ u_2 (x) &=& u_1 - \int_0^x \left[ e^s \left( u_1 (s) - u_0 (s) \right) + e^{-s} \left( u_1^3 (s) - u_0^3 (s) \right) \right] {\text d}s = y_2 (x) , \end{eqnarray*}
and so on. As we see from the above formulas, we obtain exactly the same sequence of functions, but we performed less calculations. Now we plot the third approximation along with the true solution.
ode1 = {y'[x] + Exp[x]*y[x] + Exp[-x]*(y[x])^3 == Cosh[x], y[0] == 1}
sol = NDSolve[ode1, y, {x, 0, 2}];
Plot[{Callout[y3[x], "three-term VIM approximation"], Callout[y[x] /. sol, "exact"]}, {x, 0, 1}, PlotStyle->Thick]
1. Abbasbandy S (2007), A new application of He’s variational iteration method for quadratic Riccati differential equation by using Adomian’s polynomials, Journal of Applied Mathematics and Computing, 207:59–63
2. Tamer A. Abassy, New Applications of the Modified Variational Iteration Method, Studies in Nonlinear Sciences, 2012, Vol. 3, No. 2, pp. 49--61, doi: 10.5829/idosi.sns.2012.3.2.246
3. Tamer A. Abassy, Magdy A. El-Tawil, and H. El-Zoheiry, Exact solutions of some nonlinear partial differential equations using the variational iteration method linked with Laplace transforms and the Padé technique, Computers and Mathematics with Applications, 2007, Vol. 54, pp. 940--954.
4. M.A. Fariborzi Araghi, S. Gholizadeh Dogaheh, and Z. Sayadi, The modified variational iteration method for solving linear and nonlinear ordinary differential equations, Australian Journal of Basic and Applied Science, 2011, Vol. 5, No. 10, pp. 406--416.
5. Inan Ates and Ahmet Yildirim, "Comparison between variational iteration method and homotopy perturbation method for linear and nonlinear partial differential equations with the nonhomogeneous initial conditions," Numerical Methods for Partial Differential Equations, 2010, Vol 26, Issue 6, pp. 1581--1593, doi: https://doi.org/10.1002/num.20511
6. Elzaki, T.M., Solution of Nonlinear Partial Differential Equations by New Laplace Variational Iteration Method, 2018, doi: 10.5772/intechopen.73291
7. Ganji D.D. and Sadighi A. (2006) Application of He’s methods to nonlinear coupled systems of reactions. International Journal of Nonlinear Sciences and Numerical Simulation, 2006, Vol, 7, No 4, pp. 411--418.
8. Ganji DD, Tari H, Jooybari MB (2007) Variational iteration method and homotopy perturbation method for evolution equations. Comput Math Appl 54:1018–1027
9. Ji-Huan He, "Variational iteration method---a kind of non-linear analytical technique: some examples," International Journal of Non-Linear Mechanics, 1999, Vol 34, 699--708.
10. Ji-Huan He, "Variational iteration method for autonomous ordinary differential systems," Applied Mathematics and Computation, 2000, Vol. 114, 115--123.
11. Ameina S. Nuseir and Mohammad Al-Towaiq, The modified variational iteration method for solving the Impenetrable Agar model problem, International Journal of Pure and Applied Mathematics, 2014, Vol. 96, No 4, pp. 445--456, doi: http://dx.doi.org/10.12732/ijpam.v96i4.3
12. S. S. Samaee, O. Yazdanpanah, and D. D. Ganji, "New approaches to identification of the Lagrange multiplier in the variational iteration method," Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2914, doi: 0.1007/s40430-014-0214-3
13. Tatari M, Dehghan M (2007), "He’s variational iteration method for computing a control parameter in a semi-linear inverse parabolic equation," Chaos, Solitons & Fractals, 2007, Vol. 33, Issue 2, pp. 671--677.
14. Tari H., Ganji D.D., and Babazadeh H. (2007) "The application of He’s variational iteration method to nonlinear equations arising in heat transfer," Physics Letters A, Vol. 363, No. 3, pp. 213--217, doi: 10.1016/j.physleta.2006.11.005
15. Tari H, Ganji D.D., and Rostamian M. (2007) "Approximate solution of K (2, 2), KdV and modified KdV equations by variational iteration method, homotopy perturbation method and homotopy analysis method," International Journal of Nonlinear Sciences and Numerical Simulation, 2007, Vol. 8, pp.203--210.
16. A.-M. Wazwaz (2009), "The variational iteration method for analytic treatment for linear and nonlinear ODEs," Applied Mathemathics and Computation, 2009, Vol. 212, Issue 1, pp. 120--143.
17. A.-M. Wazwaz (2006), "A study on linear and nonlinear Schrodinger equations by the variational iteration method," Chaos Solitons Fractals, doi: 10.1016/j.chaos.10.009
18. A.-M. Wazwaz (2007), The variational iteration method for the exact solution of Laplace equation, Physics Letters A, Vol. 363, No 4, pp. 260--262, doi: 10.1016/j.physleta.2006.11.014
19. A.-M. Wazwaz (2001) The numerical solution of fifth-order boundary value problems by the decomposition method. Journal of Applied Mathematics and Computing, 136:259–270
20. A.-M. Wazwaz (2007) The variational iteration method for solving linear and nonlinear systems of PDES, Computers and Mathematics with Applications, 2007, Vol. 54, pp. 895--902.