# Preface

The homotopy pertubation method was proposed by the professor of Soochow University (China) Ji-Huan He in 1999. He studied Construction Engineering in the middle 1980s in Xi’an University of Architecture & Technology, China, and received his master's degree of mechanical engineering in 1990 from Shanghai University, China, with the thesis "Reliability Analysis of Pneumatical Cylinders." Subsequently He worked as an engineer at two manufactories for about three years. He defended in 1997 a Ph.D. degree at Shanghai University.

Upon graduation, Ji-Huan He focused on analytical methods for nonlinear equations, and suggested some new approximate analytical methods, e.g., the variational iteration method, the homotopy perturbation method, and the parameter-expansion method, which are now widely used to solve various nonlinear equations. In 2002, Ji-Huan He moved to Donghua University conducting research on nanotechnology. In paricular, He participated in projects aimed to develope some new devices for producing nanofibers, such as vibration-electrospinning and magneto-electrospinning.

Introduction to Linear Algebra

# Variational Iteration Method

In the previous sections (see part III and part IIV), we discussed the basic terms related to variational iteration method. Now we present its extension for higher order differential equations and systems of ODEs.

Consider the differential equation

$L \left[ y \right] = N\left[ y \right] + g(x)$
subject to some auxiliary (initial and/or boundary) conditions. Here L is a linear differential operator, N is a nonlinear operator, and g(x) is the source inhomogeneous term. The variational iteration method admits the use of a correction functional in the form
$y_{n+1} (x) = y_n (x) + \int_0^x \lambda (s) \left\{ L \left[ y \right] - N\left[ \tilde{y} \right] - g(x) \right\} {\text d}s ,$
where λ is a general Lagrange multiplier, and $$\tilde{y}$$ has a restricted variation, that is, $$\delta\tilde{y} =0 .$$ Having λ found, an iteration formula should be used for determination of the successive approximations yn of the solution y(x):
$y(x) = \lim_{n\to\infty} \, y_n (x) .$

When L is the n-th order differential operator, $$L = \texttt{D}^n = {\text d}^n / {\text d}x^n ,$$ the Lagrange multiplier could be chosen as

$\lambda = \frac{(-1)^n}{(n-1)!} \left( s-x \right)^{n-1} .$
Although the zeroth approximation y0 can be chosen arbitrary, it is preferable to select it in the form:
$y_0 (x) = y(0) + y' (0)\,x + \frac{1}{2!}\, y'' (0) \,x^2 + \cdots + \frac{1}{n-1}\, y^{(n-1)} (0)\, x^{n-1} ,$
where n is the order of the differential operator L.

Generally, for the initial value problem for n-th order differential equation

$y^{(n)} (x) = f\left( y, y', \ldots , y^{(n-1)} \right) + g(x) , \qquad y(0) = \alpha_0 , \ y' (0) = \alpha_1 , \ \ldots , \ y^{(n-1)} (0) = \alpha_{n-1} ,$
the iteration formula takes the form
$y_{n+1} (x) = y_n (x) + \frac{(-1)^n}{(n-1)!} \, \int_0^x \left( s-x \right)^{n-1} \left\{ y^{(n)} (s) - f \left( y, y', \ldots , y^{(n-1)} \right) - g(s) \right\} {\text d}s , \qquad n=0,1,2, \ldots .$

Example: Consider the second order linear ODE with constant coefficients and extend our analysis to this equation that given by

$y'' + p\,y' + q\,y = g(x) , \qquad y(0) = \alpha , \quad y' (0) = \beta ,$
where p and q are some given constants. According to the VIM, the basic character of the method is to construct a correction functional for the equation:
$y_{n+1} (x) = y_n (x) + \int_0^x \lambda (s) \left\{ y''_n (s) + p\, \tilde{y}'_n (s) + q\,\tilde{y}_n (s) - g(s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Calculating variation with respect to yn yields the following stationary conditions
$\left. \lambda '' \right\vert_{s=x} =0 , \qquad \left. 1 - \lambda ' \right\vert_{s=x} =0 , \qquad \left. \lambda \right\vert_{s=x} =0 .$
The Lagrange multiplier, therefore, can be identified as λ = s-x. Substituting this value of the Lagrange multiplier into the functional gives the iteration formula
$y_{n+1} (x) = y_n (x) + \int_0^x (s-x) \left\{ y''_n (s) + p\, y'_n (s) + q\,y_n (s) - g(s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
In order to eliminate repeated calculations, we rewrite the above formula as
$y_{n+1} (x) = y_n (x) + \int_0^x (s-x) \,y''_n (s) {\text d}s + \int_0^x (s-x) \left\{ p\, y'_n (s) + q\,y_n (s) - g(s) \right\} {\text d}s .$
Integrating the integral containing the second derivative by parts gives
$\int_0^x (s-x) \,y''_n (s) {\text d}s = -y_n (x) + y_n (0) + x\, y'_n (0) .$
This allows us to rewrite the iteration formula as
$y_{n+1} (x) = y_n (0) + x\,y'_n (0) + \int_0^x (s-x) \left\{ p\, y'_n (s) + q\,y_n (s) - g(s) \right\} {\text d}s .$
We can set $$y_0 = y_n (0) = x\,y'_n (0) = \alpha + x\,\beta ,$$ so
\begin{eqnarray*} y_{n+1} (x) &=& y_0 + \int_0^x (s-x) \left\{ p\, y'_n (s) + q\,y_n (s) - g(s) \right\} {\text d}s \\ y_0 + \int_0^x (s-x) \, p \left[ y'_n (s) - y'_{n-1} (s) \right] {\text d}s + \int_0^x (s-x) \, q \left[ y_n (s) - y_{n-1} (s) \right] {\text d}s + \int_0^x (s-x) \left\{ p\, y'_{n-1} (s) + q\,y_{n-1} (s) - g(s) \right\} {\text d}s . \end{eqnarray*}
However, we know that
$y_{n} (x) = y_0 + \int_0^x (s-x) \left\{ p\, y'_{n-1} (s) + q\,y_{n-1} (s) - g(s) \right\} {\text d}s .$
Therefore, we simplify the iteration formula to
$y_{n+1} (x) = y_n (x) + p\,\int_0^x (s-x) \left\{ y'_n (s) - y'_{n-1} (s) \right\} {\text d}s + q\,\int_0^x (s-x) \left\{ y_n (s) - y_{n-1} (s) \right\} {\text d}s , \qquad n=0, 1,2,3,\ldots .$
As usual, we set $$y_{-1} = 0.$$    ■

Example: Consider the IVP for harmonic oscillator

$\ddot{y} + \omega^2 y = f(t) \qquad\mbox{with}\quad f(t) = A\,\sin \omega t +B\,\cos t ,$
subject to the initial conditions
$y(0) =1, \qquad \dot{y} (0) =-1.$
Its correction functional can be written as follows
$y_{n+1} (t) = y_n (t) + \int_0^t \lambda \left\{ y''_n (s) + \omega^2 \,y_n (s) - f(s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Making the above correction functional stationary, and noticing that δy(0) = 0, we get
\begin{eqnarray*} \delta y_{n+1} (t) &=& \delta y_n (t) + \delta \int_0^t \lambda \left\{ y''_n (s) + \omega^2 \,y_n (s) - f(s) \right\} {\text d}s \\ &=& \delta y_n (t) + \left. \lambda (s)\,\delta y'_n (s) \right\vert_{s=t} - \left. \lambda ' (s)\, \delta y_n (s) \right\vert_{s=t} + \int_0^t \left[ \lambda '' + \omega^2 \lambda \right] \delta y_n (s) \,{\text d}s . \end{eqnarray*}
From the right hand side, it follows
$\begin{split} \lambda '' + \omega^2 \lambda &= 0 , \\ \left. \lambda (s) \right\vert_{s=t} &= 0 , \\ \left. 1 - \lambda ' (s)\right\vert_{s=t} &= 0 . \end{split}$
The Lagrange multiplier can be readily identified as
$\lambda (t, s) = \frac{1}{\omega}\,\sin \omega (s-t) .$
With this in hands, we obtain the following iteration formula
$y_{n+1} (t) = y_n (t) + \frac{1}{\omega} \int_0^t \sin \omega (s-t) \left\{ y''_n (s) + \omega^2 \,y_n (s) - f(s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Now we consider the homogeneous equation corresponding to the given driven equation $$\ddot{y} + \omega^2 y =0$$ subject to the given initial conditions; its solution, called the complementary function, can be chosen as the initial approximation
$y_{0} (t) = C_1 \cos \omega t + C_2 \sin \omega t = \cos \omega t - \frac{\sin \omega t}{\omega} ,$
where C1 and C2 are such constants to satisfy the given initial conditions, y(0) = 1 and y'(0) = -1. Then the iteration formula gives
\begin{eqnarray*} y_1 (t) &=& y_0 (t) - \frac{1}{\omega} \int_0^t \sin \omega (s-t) \left[ A\,\sin \omega s + B\,\sin \omega s \right] {\text d}s \\ &=& C_1 \cos \omega t + C_2 \sin \omega t - \frac{A}{2} \left[ t\,\cos \omega t - \frac{\sin \omega t}{\omega} \right] + \frac{B\,\omega}{\omega^2 -1} \left[ \cos \omega t - \cos t \right] , \end{eqnarray*}
which is the general solution of the given differential equation (harmonic oscillator).

However, is we apply restricted variations to the correction function, then its exact solution can be arrived at only by successive iterations. Considering the corresponding homogeneous differential equation $$\ddot{y} + \omega^2 y =0 ,$$ we rewrite the correction functional as follows:

$y_{n+1} (t) = y_n (t) + \int_0^t \lambda \left\{ y''_n (s) + \omega^2 \,\tilde{y}_n (s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
Here $$\tilde{y}_n (s)$$ is considered a restricted variation, then the stationary conditions ($$\delta\tilde{y}_n =0$$ ) of the above correction functional can be expressed as follows
$\begin{split} \lambda '' (s) &=0 , \\ \left. \lambda (s) \right\vert_{s=t} &= 0, \\ \left. 1 - \lambda ' (s) \right\vert_{s=t} &= 0 . \end{split}$
Therefore, the Lagrange multiplier becomes $$\lambda = s-t .$$ This leads to the following iteration formula
$y_{n+1} (t) = y_n (t) + \int_0^t (s-t) \left\{ y''_n (s) + \omega^2 \,y_n (s) - f(s) \right\} {\text d}s , \qquad n=0,1,2,\ldots .$
It we choose the initial iteration as $$y_0 = 1 -t$$ in order to satisfy the initial conditions, we have
\begin{eqnarray*} y_1 (t) &=& y_0 (t) + \omega^2 \int_0^t \left( s-t \right) y_0 (s)\, {\text d}s = 1-t + \omega^2 \left( -\frac{t^2}{2} + \frac{t^3}{6} \right) , \\ y_2 (t) &=& y_1 (t) + \int_0^t \left( s-t \right) \left( y''_1 (s) + \omega^2 y_1 (s) \right) {\text d}s = 1-t + \omega^2 \left( -\frac{t^2}{2} + \frac{t^3}{6} \right) + \omega^4 \left( \frac{t^4}{24} - \frac{t^5}{120} \right) , \\ y_3 (t) &=& y_2 (t) + \int_0^t \left( s-t \right) \left( y''_2 (s) + \omega^2 y_2 (s) \right) {\text d}s = 1-t + \omega^2 \left( - \frac{t^2}{2!} + \frac{t^3}{3!} \right) + \omega^4 \left( \frac{t^4}{4!} - \frac{t^5}{5!} \right) + \omega^6 \left( - \frac{t^6}{6!} + \frac{t^7}{7!} \right) , \\ y_n (t) &=& \sum_{k=0}^n \frac{(-1)^n}{(2k)!}\, \omega^{2k} t^{2k} - \frac{1}{\omega} \,\sum_{k=0}^n (-1)^k \frac{\omega^{2k+1} t^{2k+1}}{(2k+1)!} . \end{eqnarray*}
Thus, we obtain
$\lim_{n\to\infty}\, y_n (t) = \cos\omega t - \frac{1}{\omega}\,\sin \omega t ,$
which is the exact solution of the corresponding homogeneous equation.    ■

Example: We examine the Thomas--Fermi equation

$y'' = y^{3/2} x^{-1/2} \qquad y(0) =1, \quad \lim_{x\to\infty} \,y(x) =0.$
This equation was considered by Abdul-Majid Wazwaz in 2001 to model the effective nuclear charge in heavy atoms. The above nonlinear differential equation is named after the British physicist and applied mathematician Llewellyn Thomas (1903--1992) and the 1938 Nobel prize winner and the "architect of the atomic bomb" Enrico Fermi (1901--1954).

To overcome the difficulty of the fractional exponent of y, we use the transformation y(x) = 1 + u(x) and get

$u'' = \left( 1+ u \right)^{3/2} x^{-1/2} \qquad u(0) =0, \quad u' (0) = B.$
A direct application of the VIM to the given problem is problematic because the product of a polynomial and the slope function cannot be obtained in closed form. Therefore, we use slope approximation:
$\left( 1+ u \right)^{3/2} \approx 1 + \frac{3}{2}\,u + \frac{3}{8}\, u^2 - \frac{1}{16}\, u^3 .$
The initial approximation can be taken as u0 = B x. Then
$u_1 (x) = u_0 + \int_0^x (s-x) \left[ 1 + \frac{3}{2}\,u_0 + \frac{3}{8}\, u^2_0 - \frac{1}{16}\, u^3_0 \right] {\text d} s = B\,x - \frac{x^2}{2} - \frac{B}{4} x^3 - \frac{B^2}{32}\,x^4 + \frac{B^3}{320}\, x^5 .$
u0[x_] = B*x;
u1[x_] = u0[x] + Integrate[(s - x)*(1 + 3*u0[s]/2 + 3*u0[s]^2 /8 - (u0[s])^3/16), {s, 0, x}]
For next iterations, we use VIM:
\begin{eqnarray*} u_2 &=& u_1 + \int_0^x (s-x) \left[ 1 + \frac{3}{2}\,u_1 + \frac{3}{8}\, u^2_1 - \frac{1}{16}\, u^3_1 \right] {\text d} s = B\, x - \frac{B}{2}\, x^2 + \frac{x^4}{16} - \frac{B^2}{16}\, x^4 + \frac{3\,B}{80}\, x^5 + \frac{B^3}{160}\, x^5 + \cdots , \\ u_3 &=& u_2 + \int_0^x (s-x) \left[ 1 + \frac{3}{2}\,u_2 + \frac{3}{8}\, u^2_2 - \frac{1}{16}\, u^3_2 \right] {\text d} s = B\, x - \frac{B}{2}\, x^2 - \frac{3\,B}{4}\,x^3 + \frac{x^4}{16} - \frac{B^2}{16}\, x^4 - \frac{1}{32} \left( B^2 -4 \right) x^4 + \frac{3\,B}{80}\, x^5 + \frac{B^3}{160}\, x^5 + \cdots . \end{eqnarray*}
u2[x_] = u1[x] + Integrate[(s - x)*(1 + 3*u1[s]/2 + 3*u1[s]^2 /8 - (u1[s])^3/16), {s, 0, x}]
u3[x_] = u2[x] + Integrate[(s - x)*(1 + 3*u2[s]/2 + 3*u2[s]^2 /8 - (u2[s])^3/16), {s, 0, x}]
To determine B, we take the limit as x goes to infinity of the Padé approximants.    ■

Example: Consider the IVP for the pendulum equation

$\ddot{\theta} + \omega^2 \sin \theta =0 , \qquad \theta (0) = \theta_0 , \quad \dot{\theta} (0) =0.$

Example: Consider the IVP for Duffing equation    ▣

Example: Consider the Lane--Emden differential equation of index m:

$\theta'' (x) + \frac{2}{x}\,\theta' (x) + \theta^m (x) =0, \qquad \theta (0) =1, \quad \theta' (0) =0 .$
It is assumed that m ≠ 1 or 0, otherwise the equation becomes linear. This equation is a dimensionless form of Poisson's equation for the gravitational potential of a Newtonian self-gravitating, spherically symmetric, polytropic fluid. This equation was used to model the thermal behaviour of a spherical cloud of gas acting under the mutual attraction of its molecules and subject to the classical lows of thermodynamics. It is named after astrophysicists Jonathan Homer Lane (1819--1880) from the US and Robert Emden (1862--1940) from Switzerland.

Because of singularity at x = 0, we use the transformation

$u(x) = x\,\theta (x) \qquad \Longrightarrow \qquad \theta (x) = u(x) / x.$
Substituting into the Lane--Emden equation, we obtain
$u'' (x) + x^{1-m} u^m =0 , \qquad u(0) =0, \quad u' (0) =1.$
Using the MVIM, we derive the following iterative sequence
$u_{n+1} (x) = u_n (x) + \int_0^x \left( s-x \right) \left[ s^{1-m} \left( u_n^m (s) - u_{n-1}^m (s) \right) \right] {\text d} s , \qquad n=0,1,2,\ldots.$
Upon choosing u-1 = 0 and u0 = x, we find
\begin{eqnarray*} u_1 &=& u_0 (x) + \int_0^x (s-x)*s\,{\text d}s = x - \frac{x^3}{6} , \\ u_2 &=& u_1 + \int_0^x (s-x)*s^{1-m} \left[ u_1^m - u_0^m \right] {\text d}s = x - \frac{x^3}{6} + \frac{m}{120}\,x^5 , \\ u_3 &=& u_2 + \int_0^x (s-x)*s^{1-m} \left[ u_2^m - u_1^m \right] {\text d}s = x - \frac{x^3}{6} + \frac{m}{120}\,x^5 - \frac{m(8m-5)}{3\, 7!}\, x^7 , \\ u_4 &=& u_3 + \int_0^x (s-x)*s^{1-m} \left[ u_3^m - u_2^m \right] {\text d}s = x - \frac{x^3}{6} + \frac{m}{120}\,x^5 - \frac{m(8m-5)}{3\, 7!}\, x^7 + \frac{m \left( 70-183m+122 m^2 \right)}{9\, 9!}\, x^9 . \end{eqnarray*}
u0[x_] = x;
u1[x_] = u0[x] + Integrate[(s - x)*s, {s, 0, x}];
u2[x_] = u1[x] + Integrate[(s - x)*s^(1 - m) *(u1[s]^m - u0[s]^m), {s, 0, x}, Assumptions -> x > 0]
u3[x_] = u2[x] + Integrate[(s - x)*s^(1 - m) *(u2[s]^m - u1[s]^m), {s, 0, x}, Assumptions -> x > 0]

Example: Consider the third order differential equation

$u''' + \frac{3}{2+x}\left( u' \right)^2 = e^u , \qquad u(0) = 1, \quad u' (0) =-2, \quad u'' (0) =2.$
Using the Lagrange multiplier as $$\lambda = - \frac{1}{2} \left( s-x \right)^2 ,$$ and the initial approximation
$u_0 = u(0) + x\, u' (0) + \frac{x^2}{2}\, y'' (0) = 1 -2x + x^2 ,$
we set the iteration formula
$u_{n+1} (x) = u_n (x) - \frac{1}{2} \int_0^x \left( s-t \right)^2 \left[ u'''_n (s) - \frac{3}{2+s} \left( u'_n \right)^2 - e^{u_n}] \right] {\text d} s , \qquad n=0,1,2,\ldots .$
Calculations show
\begin{eqnarray*} u_0 (x) &=& 1 -2x + x^2 = (1-x)^2 , \\ u_1 (x) &=& u_0 (x) + \frac{1}{2} \int_0^x (s-x)^2 \left[ \frac{3}{2+s} \left( 2x-2 \right)^2 - e^{u_0 (s)}] \right] {\text d} s = \\ u_2 (x) &=& u_1 (x) - \frac{1}{2} \int_0^x (s-x)^2 \left[ u'''_1 (s) - \frac{3}{2+s} \left( u'_1 \right)^2 - e^{u_1 (s)} \right] {\text d} s = \\ u_3 (x) &=& u_2 (x) - \frac{1}{2} \int_0^x (s-x)^2 \left[ u'''_2(s) - \frac{3}{2+s} \left( u'_2 \right)^2 - e^{u_2 (s)} \right] {\text d} s = \end{eqnarray*}
u0[x_] = (1 - x)^2
u1[x_] = Simplify[ u0[x] + 1/2* Integrate[(s - x)^2 *(3/(2 + s) *(2*s - 2)^2 - Exp[u0[s]]) , {s, 0, x}]]
ConditionalExpression[ 1/8 (E (2 - 4 x) + 2 E^(-1 + x)^2 (-1 + x) + 8 (-1 + x)^2 + 4 x (-216 - 162 x - 16 x^2 + x^3) - Sqrt[\[Pi]] (1 - 4 x + 2 x^2) Erfi[1] + Sqrt[\[Pi]] (1 - 4 x + 2 x^2) Erfi[1 - x] - 432 (2 + x)^2 Log[2] + 432 (2 + x)^2 Log[2 + x]), Re[x] > -2 || x \[NotElement] Reals]

Example: Consider the system of differential equations ■

1. Abassy, Tamer A. New Applications of the Modified Variational Iteration Method, Studies in Nonlinear Sciences, 2012, Vol. 3, No. 2, doi: 10.5829/idosi.sns.2012.3.2.246
2. Chen, Y.M., Meng, G., and Liu, J.K., Variational iteration method for conservative oscillators with complicated nonlinearities, Mathematical and Computational Applications, 2010, Vol. 15, No 5, pp. 802--809.
3. Fernandez, Francisco M., On some approximate methods for nonlinear models
4. Ganji D.D., Jannatabadi M. and Mohseni E. (2006) Application of He’s variational iteration method to nonlinear Jaulent--Miodek equations and comparing it with ADM, Journal of Applied Mathematics and Computing, 2007, Vol. 207, pp. 35--45.
5. Ganji, D.D., Karimpour, S., and Ganji, S.S., Approximate analytical solution to nonlinear oscillations of non-natural systems using He's energy balance method, Progreass in Electromagnetics Research, 2008, Vol. 5, pp. 43--54.
6. Ganji, D.D., Nourollahi, M., Rostamian, M., A comparison of variational iteration method with Adomian’s Decomposition Method in some highly nonlinear equations. International Journal of Science and Technology, 2007, 2(2):179–188
7. Ghotbi, Abdoul R., Barari, A., Omidvar, M., and Ganji, D.D., Application of variational iteration method to coupled system of nonlinear partial differential physical equations, African Journal of Mathematics and Computer Science Research, 2009, Vol. 2, No. 4, pp. 063--068,
8. Gorji, M., Ganji, D.D., Soleimani, S., 2007. New application of He’s homotopy perturbation method. International Journal of Nonlinear Science and Numerical Simulation, 8(3):319–328.
9. He, Ji-Huan, Selected papers published by Ji-Huan He (email: hejihuan@suda.edu.cn) can be found on the websites:
https://works.bepress.com/ji_huan_he/
10. He, J.H., Homotopy perturbation technique, Computer Methods in Applied Mechanics and Engineering, 1999, Vol. 178, 257--262.
11. Kaya, D. and El-Sayed, S.M., A numerical method for solving Jaulent--Miodek equation, Phys. Lett. A 318 (2003) 345–353.
12. Merdan, M., He's variational iteration method for solving modelling the pollution of a system of lakes (in Turkish)
13. Rafei, M., Daniali, H., Ganji, D.D., Variational iteration method for solving the epidemic model and the prey and predator problem, Applied Mathematics and Computation, 2007, Vol 186, Issue 2, pp. 1701--1709.
14. Yildirim Ahmet and Kelleci Alev, Numerical Simulation of the Jaulent-miodek Equation by He’s Homotopy Perturbation Method, World Applied Sciences Journal, 2009, Vol. 7 (Special Issue for Applied Math), pp. 84--89.
15. Wazwaz, A.-M., The modified decomposition method and Pade approximants for solving the Thomas--Fermi equation, Applied Mathematics and Computation, 1999, Vol. 105, Issue , pp. 11--19.
16. Wazwaz, A.-M., The variational iteration method for solving two forms of Blasius equation on a half-infinite domain, Applied Mathematics and Computation, 2007, Vol. 188, Issue 1, pp. 485--491.