# Enhanced Multistage Homotopy Perturbation Method: Approximate Solutions of Nonlinear Dynamic Systems

Authors: Daniel Olvera; Alex Elías-Zúñiga
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We introduce a new approach called the enhanced multistage homotopy perturbation method (EMHPM) that is based on the homotopy perturbation method (HPM) and the usage of time subintervals to find the approximate solution of differential equations with strong nonlinearities. We also study the convergence of our proposed EMHPM approach based on the value of the control parameterh by following the homotopy analysis method (HAM). At the end of the paper, we compare the derived EMHPM approximate solutions of some nonlinear physical systems with their corresponding numerical integration solutions obtained by using the classical fourth order Runge-Kutta method via the amplitude-time response curves.


---

## Body

## 1. Introduction

With the recent progress in nonlinear problems research, there has been an increasing interest in analytical techniques to solve the corresponding nonlinear equations. However, most of the current methods based their solution methodology on the assumption of small nonlinearities that limit its potential usage in the solution of physical and engineering applications.Recent applications of some asymptotic methods such as the variational iteration, the homotopy perturbation, the energy balance, the parameter-expansion, the variational approach, the improved amplitude frequency formulation, the max-min approach, the Hamiltonian approach, and the homotopy analysis, to name a few, have been used to obtain approximate solutions of highly nonlinear problems in which the traditional perturbation methods have some limitations [1, 2].To obtain approximate analytical solutions of strongly nonlinear differential equations by using the homotopy approach, the methods of homotopy perturbation, homotopy analysis, Adomian decomposition (ADM), and the variational iteration (VIM) are commonly used in the literature. The ADM is an iterative method which provides analytical approximate solutions in the form of an infinite power series for nonlinear equations without linearization [3]. The HPM developed by He [4–6] is not limited by the assumption of small nonlinearities. The HPM is an effective method that has several applications in science and engineering to find approximate solutions of nonlinear differential equations. The HPM couples the traditional perturbation method and homotopy theory used in topology and assumes a simple initial analytical solution that gradually approximates to the solution of the differential equation through a parameter that converges to unity. The HPM decomposes a complex problem into a series of simple problems that simplifies the derivation of the approximate solutions. Therefore, the usage of the homotopy theory has been effective in deriving the approximate solutions of different types of problems, including nonlinear second order differential equations [7, 8] or nonlinear fractional equations [9] that model the behavior of physical and engineering systems. Hojjati and Jafari compared the HPM and the ADM to obtain the distributions of stresses and displacements in a rotating annular elastic disk [10]. They have concluded that, although the numerical results are almost the same, the HPM is much easier, more convenient and efficient than the ADM and the VIM. Abbasbandy in [11] concluded that one advantage of the HPM when compared to the ADM is related to its capability of achieving the approximate solution of the quadratic Riccati differential equation by considering all the Taylor expansion series terms. These results were also confirmed by Pamuk in [12].The HPM is an asymptotic method with limited convergence away from the equation initial conditions [13, 14]. In order to overcome this weakness, the HPM must be applied by intervals. For instance, Abbasbandy used this approach in solving the Riccati equation [15]. A more formal approach was developed by Hashim and Chowdhury to solve a system of first order differential equations through multiple intervals. They used this approach, called the multistage homotopy perturbation method (MHPM), to solve nonlinear ordinary differential equations with numerical predictions that follow well the corresponding numerical integration solutions [16, 17]. Based on the MHPM, Gökdogan and Merdan derived an approximate solution for the Coullet nonlinear differential equation with numerical simulations that agree well with the Runge-Kutta numerical integration method [18]. Wang and Yu extended the usage of the MHPM to obtain the approximate numeric-analytic solutions not only of the chaotic fractional order Chen system, but also of the hyperchaotic fractional order Lorenz system with numerical simulations that exhibit good accuracy at long time spans [19].Based on the findings of these previous research works, the aim of this work focuses on developing a semi-numerical-analytic technique that generalizes the MHPM in an attempt to obtain accurate approximate solutions of nonlinear differential equations. To assess the accuracy of this new approach, its numerical predictions will be compared with respect to those of the MHPM techniques introduced in [20, 21]. At the end of the paper, we will derive the approximate solutions of some nonlinear differential equations and compare their numerical predictions with respect to their corresponding numerical integration solutions obtained by using the fourth order Runge-Kutta method.


## 2. Description of the Homotopy Perturbation Method

In order to illustrate the basic idea behind the HPM, let us consider the following nonlinear differential equation:(1)A(u)-f(r)=0,r∈Ω,
that has the following boundary condition:
(2)B(u,∂u∂n)=0,r∈Γ,
where A is a differential operator, B is a boundary operator, f(r) is a known analytic function, and Γ is the boundary in domain Ω.TheA operator can be divided into two parts: L and N, where L represents the linear terms while N is related to nonlinear terms. Equation (1) can therefore be written as follows:
(3)L(u)+N(u)-f(r)=0.
Applying the principles of topology, the following homotopy is established: v(r,p):Ω×[0,1]→ℜ, satisfying the following conditions:
(4)H(v,p)=L(v)-L(u0)+pL(u0)+p[N(v)-f(r)]=0,
where p∈[0,1] is an embedding parameter and u0 is the initial approximation to (1) that satisfies the boundary condition. It is evident from (4) that setting p from zero to unity makes v(r,p) change from u0(r) to u(r), that is, by letting
(5)H(v,0)=L(v)-L(u0)=0H(v,1)=A(v)-f(r)=0.
In topology, this is calleddeformation and, in the same context, L(v)-L(u0) and A(v)-f(r) are calledhomotopic. He in [4] assumed that the solution of (4) can be expressed with a series of p power terms of the form
(6)v=v0+pv1+p2v2+⋯.
Thus, by considering p=1, the approximate solution of (1) yields
(7)u=lim⁡p→1=v0+v1+v2+⋯.
In order to test the precision achieved with the election of the linear operator, we proceed to analyze the following nonlinear equation [20]:
(8)y˙-y+ety2=0y(0)=c,
where y is the independent variable, t is the running time, and c is a constant. First, the homotopy for (8) is established by following the standard homotopy methodology that requires the selection of an operator whose solution is unstable, such as L(y)=y˙-y, and therefore, N(y)=ety2. Thus, the homotopic representation of (8) is given as
(9)H(p,Y)=L(Y)-L(y0)+pL(y0)+pN(Y)=0,y0(t)=c.
By substituting the second order expansion Y=Y0+pY1+p2Y2 into (9) and by grouping the terms with the same p order, we obtain the following set of equations:
(10)p0:Y˙0-Y0=0Y0(0)=cp1:Y˙1-Y1=-c2e3tY1(0)=0p2:Y˙2-Y2=c3(e5t-e3t)Y2(0)=0
whose solutions are given by
(11)Y0=cetY1=12c2(et-e3t)Y2=c3(14et-12e3t+14e5t).
The approximate solution of (8), given by (11), corresponds to the second order solution in [20] when c=1.


### 2.1. The Enhanced Homotopy Perturbation Method

Hosein Nia et al. proposed to modify the HPM by introducing a stable linear operator. This modified method was called theenhanced homotopy perturbation method. By using the EHPM, (8) is rewritten by considering a stable linear operator of the form L(y)=y˙+y, while the nonlinear operator is defined as N(y)=ety2-2y. These operator definitions provide the following expression:
(12)y˙+y+ety2-2y=0y(0)=c.
In this case, the homotopic representation of (12) is similar to that given by (9). By following the previous procedure and by regrouping the terms with the same p order, we obtain that
(13)p0:Y˙0+Y0=0Y0(0)=cp1:Y˙1+Y1=(2c-c2)e-tY1(0)=0p2:Y˙2+Y2=(4c-6c2+2c3)te-tY2(0)=0.
Integration of (13) provides the following expressions:
(14)Y0(t)=ce-tY1(t)=(2c-c2)te-tY2(t)=(2c-3c2+c3)t2e-t.
Notice that the second order solution of (14) agrees well with the solutions obtained in [20] when c=1.


### 2.2. Representation of the HPM by Using Taylor Series

In order to improve the convergence of the HPM, Odibat proposed the expansion of the independent variable by using Taylor series [21]. Thus, (8) can be written as
(15)y˙-y+∑n=0∞tnn!y2=0y(0)=c.
By using the linear operator L(y0)=(d/dt)y0, the homotopic representation of (15) can be expressed as
(16)H(p,Y)=Y˙-y˙0+py˙0-pY+∑n=0∞pn+1tnn!Y2.
In order to get the first approximation of Y0=c, the term pn+1 has to be included in (16). Substituting the third order expansion of (15) into (16) and by regrouping the terms with the same p order, we get that
(17)p0:Y˙0-y˙0=0,Y0(0)=cp1:Y˙1-Y0+Y02=0,Y1(0)=0p2:Y˙2-Y1+2Y0Y1+tY02=0,Y2(0)=0p3:Y˙3-Y2+12t2Y02+2tY0Y1+2Y0Y2+Y12=0,Y3(0)=0p4:Y˙4-Y3+16t3Y02+t2Y0Y1+2tY0Y2+tY12+2Y0Y3+2Y1Y2=0,Y4(0)=0.
Thus, the fourth order solution of (8) can be obtained by integrating (17); this yields
(18)Y0=cY1=(c-c2)tY2=(c3-2c2+12c)t2Y3=(-c4+3c3-136c2+16c)t3Y4=(c5-4c4+296c3-53c2+124c)t4.
Figure 1 illustrates the approximate solutions of (8) obtained by using the EHPM and the HPM solutions derived, respectively, by He and Odibat. Notice from Figure 1 that the solution obtained by using Hosein Nia et al. approach tends to follow the Runge-Kutta numerical integration solution provided by the MATLAB functionode45, while the He and the Odibat solutions tend to diverge at increasing values of t. However if we solve (8) by using the MHPM proposed by Hashim and Chowdhury [16], the solutions given by (11), (14), and (18) do not converge to the exact numerical solution because of the selection of the form of the linear operator that does not guarantee the convergence of the MHPM solution. A similar conclusion can be drawn by using the Taylor series to expand the terms of the independent variable (Odibat). Details of the computation algorithms used to obtain the numerical predictions illustrated in Figure 2 are provided in Algorithm 1.Algorithm 1:The EMHPM MATLAB code.
f
unction revision_Hosein_paper
%% Inputs
tspan=
[
0,10
]
; % time span
x0=1; % initial condition
Delta_t=10; % Delta time (equal to tspan for no multistage)
points=100; % number of points to evaluate each Delta time
%% Application of the EMHPM algorithm to the proposed solutions
[
t1, z1
]
=Semhpm(@he_stable,tspan,x0,Delta_t,points); % EMHPM by He's solution
[
t2, z2
]
=Semhpm(@hn_unstable,tspan,x0,Delta_t,points); % EMHPM by Hosein Nia's solution
[
t3, z3
]
=Semhpm(@od_taylor,tspan,x0,Delta_t,points); % EMHPM by Odibat's solution
[
t4, z4
]
=ode45(@numeric,tspan,x0); % numerical ode45 solution
%% Display computer plots
plot(t1,z1,'--',t2,z2,':',t3,z3,'-.',t4,z4,'-');
legend('He -∖ity_2','Hosein Nia et al. -∖ity_2','Odibat -∖ity_4','Numerical - ode45');
ylabel('Displacement,∖itx');xlabel(Time, t'); axis([0 10 0 2])
%% Solution Definitions
function dydt=numeric(t,y) % numerical solution
dydt(1,1)=y(1)-exp(t).*y(1)∧2;
function y=he_stable(t,x0,T) % homotopy/multistage unstable, He's solution
y0=x0*exp(t);
y1=.5*x0
∧
2*(exp(t)-exp(3*t));
y2=.25*x0
∧
3*(exp(t)-2*exp(3*t)+exp(5*t));
y=y0+y1+y2;
function y=hn_unstable(t,x0,T) % homotopy/multistage stable, Hosein Nia's solution
y0=x0*exp(-t);
y1=(2*x0-x0
∧
2)*t.*exp(-t);
y2=(2*x0-3*x0
∧
2+x0
∧
3).*t.
∧
2.*exp(-t);
y=y0+y1+y2;
function y=od_taylor(t,c,T) % homotopy/multistage Taylor Expansion, Odibat's solution
y0=c;
y1=(c-c.
∧
2)*t;
y2=(c.
∧
3-2*c.
∧
2+c/2)*t.
∧
2;
y3=(-c.
∧
4+3*c.
∧
3-13*c.
∧
2/6+c/6)*t.
∧
3;
y4=(c.
∧
5-4*c.
∧
4+29*c.
∧
3/6-5*c.
∧
2/3+c/24)*t.
∧
4;
y=y0+y1+y2+y3+y4;
%% Enhanced Multistage Homotopy Perturbation Method
function
[
t,z
]
=Semhpm(nde,tspan,z0,Deltat,pnts)
tini=tspan(1); tfin=tspan(end); tstart=tini;
tini=tini-tstart; tfin=tfin-tstart;
if tini == tfin
error('The last entry in tspan must be different from the first entry.');
elseif abs(tini)> abs(tfin)
tspan=flipud(fliplr(tspan)); tini=tspan(1); tfin=tspan(end);
end
tdir=sign(tfin-tini);
if any( tdir*diff(tspan) <= 0 )
error('The entries in tspan must strictly increase or decrease.');
end
incT=Deltat/pnts;
z(1,:)=z0'; t(1)=tini;
iteT=2;
t(iteT)=tini+incT*tdir; % real time
while tdir*t(iteT)<tdir*(tfin+tdir*incT)
P_act=ceil(.99999*(t(iteT)-tini)*tdir/Deltat); % next subinterval
c=z((Deltat/incT)*(P_act-1)+1,:); % initial conditions
tsub(iteT-1)=t(iteT)-(P_act-1)*Deltat*tdir; % shifted time
temp=tsub(iteT-1);
z(iteT,:)=nde(temp,c,tstart+t(iteT)); % evaluate the solution
if tdir*t(iteT)>=tdir*tfin*0.99999
break % last subinterval
else
iteT=iteT+1;
t(iteT)=tini+(iteT-1)*incT*tdir; % next subinterval
end
end
t=t'+tstart;Figure 1
Comparison between theode45 Runge-Kutta algorithm with respect to some homotopy perturbation approximate solutions of (8).Numerical comparisons between the Runge-Kutta method (ode45) and the numerical predictions provided by the approximate solutions of (8) derived by using the MHPM approach. Here Δt=0.1.
(a)
(b)


## 2.1. The Enhanced Homotopy Perturbation Method

Hosein Nia et al. proposed to modify the HPM by introducing a stable linear operator. This modified method was called theenhanced homotopy perturbation method. By using the EHPM, (8) is rewritten by considering a stable linear operator of the form L(y)=y˙+y, while the nonlinear operator is defined as N(y)=ety2-2y. These operator definitions provide the following expression:
(12)y˙+y+ety2-2y=0y(0)=c.
In this case, the homotopic representation of (12) is similar to that given by (9). By following the previous procedure and by regrouping the terms with the same p order, we obtain that
(13)p0:Y˙0+Y0=0Y0(0)=cp1:Y˙1+Y1=(2c-c2)e-tY1(0)=0p2:Y˙2+Y2=(4c-6c2+2c3)te-tY2(0)=0.
Integration of (13) provides the following expressions:
(14)Y0(t)=ce-tY1(t)=(2c-c2)te-tY2(t)=(2c-3c2+c3)t2e-t.
Notice that the second order solution of (14) agrees well with the solutions obtained in [20] when c=1.


## 2.2. Representation of the HPM by Using Taylor Series

In order to improve the convergence of the HPM, Odibat proposed the expansion of the independent variable by using Taylor series [21]. Thus, (8) can be written as
(15)y˙-y+∑n=0∞tnn!y2=0y(0)=c.
By using the linear operator L(y0)=(d/dt)y0, the homotopic representation of (15) can be expressed as
(16)H(p,Y)=Y˙-y˙0+py˙0-pY+∑n=0∞pn+1tnn!Y2.
In order to get the first approximation of Y0=c, the term pn+1 has to be included in (16). Substituting the third order expansion of (15) into (16) and by regrouping the terms with the same p order, we get that
(17)p0:Y˙0-y˙0=0,Y0(0)=cp1:Y˙1-Y0+Y02=0,Y1(0)=0p2:Y˙2-Y1+2Y0Y1+tY02=0,Y2(0)=0p3:Y˙3-Y2+12t2Y02+2tY0Y1+2Y0Y2+Y12=0,Y3(0)=0p4:Y˙4-Y3+16t3Y02+t2Y0Y1+2tY0Y2+tY12+2Y0Y3+2Y1Y2=0,Y4(0)=0.
Thus, the fourth order solution of (8) can be obtained by integrating (17); this yields
(18)Y0=cY1=(c-c2)tY2=(c3-2c2+12c)t2Y3=(-c4+3c3-136c2+16c)t3Y4=(c5-4c4+296c3-53c2+124c)t4.
Figure 1 illustrates the approximate solutions of (8) obtained by using the EHPM and the HPM solutions derived, respectively, by He and Odibat. Notice from Figure 1 that the solution obtained by using Hosein Nia et al. approach tends to follow the Runge-Kutta numerical integration solution provided by the MATLAB functionode45, while the He and the Odibat solutions tend to diverge at increasing values of t. However if we solve (8) by using the MHPM proposed by Hashim and Chowdhury [16], the solutions given by (11), (14), and (18) do not converge to the exact numerical solution because of the selection of the form of the linear operator that does not guarantee the convergence of the MHPM solution. A similar conclusion can be drawn by using the Taylor series to expand the terms of the independent variable (Odibat). Details of the computation algorithms used to obtain the numerical predictions illustrated in Figure 2 are provided in Algorithm 1.Algorithm 1:The EMHPM MATLAB code.
f
unction revision_Hosein_paper
%% Inputs
tspan=
[
0,10
]
; % time span
x0=1; % initial condition
Delta_t=10; % Delta time (equal to tspan for no multistage)
points=100; % number of points to evaluate each Delta time
%% Application of the EMHPM algorithm to the proposed solutions
[
t1, z1
]
=Semhpm(@he_stable,tspan,x0,Delta_t,points); % EMHPM by He's solution
[
t2, z2
]
=Semhpm(@hn_unstable,tspan,x0,Delta_t,points); % EMHPM by Hosein Nia's solution
[
t3, z3
]
=Semhpm(@od_taylor,tspan,x0,Delta_t,points); % EMHPM by Odibat's solution
[
t4, z4
]
=ode45(@numeric,tspan,x0); % numerical ode45 solution
%% Display computer plots
plot(t1,z1,'--',t2,z2,':',t3,z3,'-.',t4,z4,'-');
legend('He -∖ity_2','Hosein Nia et al. -∖ity_2','Odibat -∖ity_4','Numerical - ode45');
ylabel('Displacement,∖itx');xlabel(Time, t'); axis([0 10 0 2])
%% Solution Definitions
function dydt=numeric(t,y) % numerical solution
dydt(1,1)=y(1)-exp(t).*y(1)∧2;
function y=he_stable(t,x0,T) % homotopy/multistage unstable, He's solution
y0=x0*exp(t);
y1=.5*x0
∧
2*(exp(t)-exp(3*t));
y2=.25*x0
∧
3*(exp(t)-2*exp(3*t)+exp(5*t));
y=y0+y1+y2;
function y=hn_unstable(t,x0,T) % homotopy/multistage stable, Hosein Nia's solution
y0=x0*exp(-t);
y1=(2*x0-x0
∧
2)*t.*exp(-t);
y2=(2*x0-3*x0
∧
2+x0
∧
3).*t.
∧
2.*exp(-t);
y=y0+y1+y2;
function y=od_taylor(t,c,T) % homotopy/multistage Taylor Expansion, Odibat's solution
y0=c;
y1=(c-c.
∧
2)*t;
y2=(c.
∧
3-2*c.
∧
2+c/2)*t.
∧
2;
y3=(-c.
∧
4+3*c.
∧
3-13*c.
∧
2/6+c/6)*t.
∧
3;
y4=(c.
∧
5-4*c.
∧
4+29*c.
∧
3/6-5*c.
∧
2/3+c/24)*t.
∧
4;
y=y0+y1+y2+y3+y4;
%% Enhanced Multistage Homotopy Perturbation Method
function
[
t,z
]
=Semhpm(nde,tspan,z0,Deltat,pnts)
tini=tspan(1); tfin=tspan(end); tstart=tini;
tini=tini-tstart; tfin=tfin-tstart;
if tini == tfin
error('The last entry in tspan must be different from the first entry.');
elseif abs(tini)> abs(tfin)
tspan=flipud(fliplr(tspan)); tini=tspan(1); tfin=tspan(end);
end
tdir=sign(tfin-tini);
if any( tdir*diff(tspan) <= 0 )
error('The entries in tspan must strictly increase or decrease.');
end
incT=Deltat/pnts;
z(1,:)=z0'; t(1)=tini;
iteT=2;
t(iteT)=tini+incT*tdir; % real time
while tdir*t(iteT)<tdir*(tfin+tdir*incT)
P_act=ceil(.99999*(t(iteT)-tini)*tdir/Deltat); % next subinterval
c=z((Deltat/incT)*(P_act-1)+1,:); % initial conditions
tsub(iteT-1)=t(iteT)-(P_act-1)*Deltat*tdir; % shifted time
temp=tsub(iteT-1);
z(iteT,:)=nde(temp,c,tstart+t(iteT)); % evaluate the solution
if tdir*t(iteT)>=tdir*tfin*0.99999
break % last subinterval
else
iteT=iteT+1;
t(iteT)=tini+(iteT-1)*incT*tdir; % next subinterval
end
end
t=t'+tstart;Figure 1
Comparison between theode45 Runge-Kutta algorithm with respect to some homotopy perturbation approximate solutions of (8).Numerical comparisons between the Runge-Kutta method (ode45) and the numerical predictions provided by the approximate solutions of (8) derived by using the MHPM approach. Here Δt=0.1.
(a)
(b)


## 3. Description of the Enhanced Multistage Homotopy Perturbation Method

The proposedenhanced multistage homotopy perturbation method (EMHPM) is an algorithm that utilizes the HPM solution in time subintervals, based on the following transformation: u(t)→ui(T), where ui(T) is the approximate solution of the ith subinterval that satisfies the initial solution ui(0)=ui-1(ti-1). The new time shifted variable T satisfies the following condition: 0≤T≤ti-ti-1. Also, we assume that the proposed solution can be given by ui0=ui-1(ti-1), where ti-1 represents the time at the end of the previous interval. This implies that the final value of the approximate solution in a given subinterval represents the initial condition of the next subinterval.In order to establish a homotopy algorithm that will allow us to obtain approximate solutions of nonlinear differential equations, we next assume the following.(1)
The linear operatorL(ui) can be chosen as L(ui)=(d/dT)ui, where the proposed initial solution ui0(T) is equal to the initial condition ui-1(ti-1); that is, ui0(T)=ui-1(ti-1). For simplicity in the notation, we write ui-1≡ui-1(ti-1).(2)
Since the homotopy is defined in theith subinterval, the following relation holds: T=t-ti-1, which satisfies the condition 0≤T≤(ti-ti-1), where t represents the current time. Therefore, the approximation of kth order can be obtained by integrating with respect to T, while the terms related to the independent variable t are assumed constant in the ith subinterval.Therefore, the approximate solution of order m for the differential equation can be written as
(19)ui(T,ui-1)=∑k=0mUik(T,ui-1).
Note that the solution ui(T,ui-1) is valid only in the ith subinterval (ti-1,ti]. In order to obtain the solution u(t) in the ith subinterval, we consider that the relation
(20)u(t)≈ui(t-ti-1)
holds. Therefore, the approximate solution of u at the time ti is given by
(21)ui(t-ti-1)|t=ti=ui+1(t-ti)|t=ti=ui+1(0)=ui.
In summary, each solution u(t) for a given subinterval (ta,tb) is in turn divided into j subintervals, which do not necessarily have to be equally spaced: [t0,t1],[t1,t2],…,[tj-1,tj]. Finally, the approximate solution for u(t) is obtained by coupling the solutions ui(t).


### 3.1. Solution Based on the EMHPM

In order to demonstrate the effectiveness of the proposed EMHPM, we next derive the approximate solution of (8). Notice that, in this case, the homotopic representation of (8) is given by (9), together with the linear operator L(yi)=(d/dT)yi and the nonlinear operator N(yi)=etyi2-yi. Utilizing the third order expansion, we can write the set of first order linear equations that result from regrouping terms that correspond to the same p order
(22)p0:Y˙i0=0Yi0(0)=cp1:Y˙i1=Yi0-etYi02Yi1(0)=0p2:Yi2=Yi1-2etYi0Yi1Yi2(0)=0p3:Y˙i3=Yi2-et(2Yi0Yi2+Yi12)Yi3(0)=0.
By solving (22), we get that
(23)Yi0=cYi1=cT-c2TetYi2=12cT2(2c2e2t-3cet+1)Yi3=-16cT3(6c3e3t-12c2e2t+7cet-1).
Note that the solutions given by (23) can be generalized in a recursive expression of the form
(24)Yik=Tk(Yi(k-1)-et∑n1=0k-1Yin1Yi(k-1-n1)),
where k is an integer number that is bigger than zero.Figure3 shows the comparison between the approximate solution of (8) obtained by using the EMHPM (utilizing the third approximation y3 and time subintervals of the same size Δt=0.1) and its Runge-Kutta numerical integration solution (ode45). We can notice from Figure 3 that now the third order EMHPM approximate solution converges to the numerical integration solution (ode45).Figure 3
Comparison of theode45 and the EMHPM numerical integration solutions.Before we examine the application of the EMHPM to derive approximate solutions of nonlinear differential equation, we will next address some convergence issues related to our proposed approach.


### 3.2. Convergence of the EMHPM

Recently Turkyilmazoglu introduced a convergence scheme for the homotopy series [22]. He showed that the value of the convergence control parameter h has strong influence on the accuracy of the solutions obtained from the homotopy analysis method (HAM). Also, Turkyilmazoglu showed that the constant ℏ-curves approach developed by Liao and Sherif [23] is equivalent to his convergence approach. To prove his findings, Turkyilmazoglu studied the convergence of the solution of several equations by applying the HAM and he discussed some of the convergence limitations of the HPM solutions based on the value of the control parameter h. By following Turkyilmazoglu results, we next study the convergence of our proposed EMHPM approach by obtaining the approximate solutions of some nonlinear differential equations and prove that our EMHPM approximate solutions converge to the Runge-Kutta numerical integration solutions and to the solutions derived by using the enhanced multistage homotopy analysis method (EMHAM).First, we focus on the solution of (8) by using our proposed EMHPM approach, to get that
(25)H(p,Y,h)=(1-p)ddT[Yi-y0]-ph[ddTYi+etYi2-Yi]=0,y0(t)=c.
Utilizing the third order expansion of (25) yields
(26)Yi0=cYi1=cTh(-1+cet)Yi2=ch(-1+cet)(T+Th-T2h2+cetT2h)Yi3=12ch(-1+cet)×(2T+4Th+2Th2+2(-1+2cet)T2h13T3h2+2(-1+2cet)T2h2+13(1-6cet+6c2e2t)T3h2).
Notice that if we take the value of h=-1 in (26) then this equation becomes similar to (23).Figure4 shows the absolute error between the numerical solution obtained from theode45 subroutine in MATLAB and the solution obtained from the EMHAM as a function of the convergent factor h. Notice from Figure 4 that the minimum absolute error occurs at h=-1. Also, at this value of the convergent control parameter h, the solutions provided by the EMHPM and the EMHAM have the same absolute error.Figure 4
Absolute error plotted against the control parameterh.We next consider the second order nonlinear ordinary differential equation used in [22] given as
(27)2u¨+u-u2=0,u(0)=0,u(∞)=1
which has an exact solution of the form
(28)u(x)=12(-1+3tanh[14(2x+4arctanh[13])]2)
with an exact homotopic solution (HAM) given as
(29)u(x,p)=(12x1+p-3+2p+3tanh[12(x1+p+2arctanh[3-2p3])]2)×(2p)-1.
We next derive the approximate solutions of (27) by using both the EMHPM and the EMHAM approaches and compute its corresponding absolute error. Notice that the approximate solution of (27) obtained by applying the EMHPM approach provides the following expressions:
(30)U1k=TkU2k-1,U2k=T2k(∑n1=0k=1U1n1U1k-n1-1-U1k-1),
where U1k is the amplitude and U2k is the first derivate of k-order. Figure 5 shows the influence of the subinterval size as well as that of the k-order solutions obtained by the EMHPM approach and its convergence to the exact solution. Notice that the minimum relative error occurs at the subinterval size of 0.1 with k=7. That implies that the convergence could be reaching by increasing the order solution and/or by decreasing the subinterval size.Figure 5
Influence of the subinterval sizeΔt and the k-order solution.Figure6 shows the predicted absolute error obtained by comparing the EMHAM solutions with respect to the exact solution, given by (28), as a function of the convergence factor h by considering the subinterval size of Δt=0.01 with k=7. We can see from Figure 6 that the absolute error has a minimum value at h=-1. Thus, our proposed EMHPM approach, which is a special case of the EMHAM at h=-1, overlaps the exact solution given by (28).Figure 6
Absolute error plotted against the control parameterh.To further evaluate the accuracy of our proposed EMHPM approach, in the next section, we will derive the approximate solution of some nonlinear differential equations that arise in several engineering applications.


## 3.1. Solution Based on the EMHPM

In order to demonstrate the effectiveness of the proposed EMHPM, we next derive the approximate solution of (8). Notice that, in this case, the homotopic representation of (8) is given by (9), together with the linear operator L(yi)=(d/dT)yi and the nonlinear operator N(yi)=etyi2-yi. Utilizing the third order expansion, we can write the set of first order linear equations that result from regrouping terms that correspond to the same p order
(22)p0:Y˙i0=0Yi0(0)=cp1:Y˙i1=Yi0-etYi02Yi1(0)=0p2:Yi2=Yi1-2etYi0Yi1Yi2(0)=0p3:Y˙i3=Yi2-et(2Yi0Yi2+Yi12)Yi3(0)=0.
By solving (22), we get that
(23)Yi0=cYi1=cT-c2TetYi2=12cT2(2c2e2t-3cet+1)Yi3=-16cT3(6c3e3t-12c2e2t+7cet-1).
Note that the solutions given by (23) can be generalized in a recursive expression of the form
(24)Yik=Tk(Yi(k-1)-et∑n1=0k-1Yin1Yi(k-1-n1)),
where k is an integer number that is bigger than zero.Figure3 shows the comparison between the approximate solution of (8) obtained by using the EMHPM (utilizing the third approximation y3 and time subintervals of the same size Δt=0.1) and its Runge-Kutta numerical integration solution (ode45). We can notice from Figure 3 that now the third order EMHPM approximate solution converges to the numerical integration solution (ode45).Figure 3
Comparison of theode45 and the EMHPM numerical integration solutions.Before we examine the application of the EMHPM to derive approximate solutions of nonlinear differential equation, we will next address some convergence issues related to our proposed approach.


## 3.2. Convergence of the EMHPM

Recently Turkyilmazoglu introduced a convergence scheme for the homotopy series [22]. He showed that the value of the convergence control parameter h has strong influence on the accuracy of the solutions obtained from the homotopy analysis method (HAM). Also, Turkyilmazoglu showed that the constant ℏ-curves approach developed by Liao and Sherif [23] is equivalent to his convergence approach. To prove his findings, Turkyilmazoglu studied the convergence of the solution of several equations by applying the HAM and he discussed some of the convergence limitations of the HPM solutions based on the value of the control parameter h. By following Turkyilmazoglu results, we next study the convergence of our proposed EMHPM approach by obtaining the approximate solutions of some nonlinear differential equations and prove that our EMHPM approximate solutions converge to the Runge-Kutta numerical integration solutions and to the solutions derived by using the enhanced multistage homotopy analysis method (EMHAM).First, we focus on the solution of (8) by using our proposed EMHPM approach, to get that
(25)H(p,Y,h)=(1-p)ddT[Yi-y0]-ph[ddTYi+etYi2-Yi]=0,y0(t)=c.
Utilizing the third order expansion of (25) yields
(26)Yi0=cYi1=cTh(-1+cet)Yi2=ch(-1+cet)(T+Th-T2h2+cetT2h)Yi3=12ch(-1+cet)×(2T+4Th+2Th2+2(-1+2cet)T2h13T3h2+2(-1+2cet)T2h2+13(1-6cet+6c2e2t)T3h2).
Notice that if we take the value of h=-1 in (26) then this equation becomes similar to (23).Figure4 shows the absolute error between the numerical solution obtained from theode45 subroutine in MATLAB and the solution obtained from the EMHAM as a function of the convergent factor h. Notice from Figure 4 that the minimum absolute error occurs at h=-1. Also, at this value of the convergent control parameter h, the solutions provided by the EMHPM and the EMHAM have the same absolute error.Figure 4
Absolute error plotted against the control parameterh.We next consider the second order nonlinear ordinary differential equation used in [22] given as
(27)2u¨+u-u2=0,u(0)=0,u(∞)=1
which has an exact solution of the form
(28)u(x)=12(-1+3tanh[14(2x+4arctanh[13])]2)
with an exact homotopic solution (HAM) given as
(29)u(x,p)=(12x1+p-3+2p+3tanh[12(x1+p+2arctanh[3-2p3])]2)×(2p)-1.
We next derive the approximate solutions of (27) by using both the EMHPM and the EMHAM approaches and compute its corresponding absolute error. Notice that the approximate solution of (27) obtained by applying the EMHPM approach provides the following expressions:
(30)U1k=TkU2k-1,U2k=T2k(∑n1=0k=1U1n1U1k-n1-1-U1k-1),
where U1k is the amplitude and U2k is the first derivate of k-order. Figure 5 shows the influence of the subinterval size as well as that of the k-order solutions obtained by the EMHPM approach and its convergence to the exact solution. Notice that the minimum relative error occurs at the subinterval size of 0.1 with k=7. That implies that the convergence could be reaching by increasing the order solution and/or by decreasing the subinterval size.Figure 5
Influence of the subinterval sizeΔt and the k-order solution.Figure6 shows the predicted absolute error obtained by comparing the EMHAM solutions with respect to the exact solution, given by (28), as a function of the convergence factor h by considering the subinterval size of Δt=0.01 with k=7. We can see from Figure 6 that the absolute error has a minimum value at h=-1. Thus, our proposed EMHPM approach, which is a special case of the EMHAM at h=-1, overlaps the exact solution given by (28).Figure 6
Absolute error plotted against the control parameterh.To further evaluate the accuracy of our proposed EMHPM approach, in the next section, we will derive the approximate solution of some nonlinear differential equations that arise in several engineering applications.


## 4. Approximate Solutions of Some Nonlinear Ordinary Differential Equations

### 4.1. The Helmholtz-Duffing Equation

In this section, we will explore the accuracy of the EMHPM approach when this is applied to obtain the approximate solution of the Helmholtz-Duffing equation that has quadratic and cubic nonlinearities. This equation is used to describe the nonlinear response of some materials in mechanical engineering applications [24]. The differential equation for the Helmholtz-Duffing oscillator is given as
(31)x¨+A1x+B1x2+εx3+D1=0,x(0)=x0,x˙(0)=x˙0=0,
where x is the system displacement, A1, B1, D1, and ε are physical parameters, and x0 and x˙0 are the system initial conditions. By introducing the following change of variables: z1i=xi, z2i=x˙i, we can write (31) in state space form:
(32)z˙1i=z2iz˙2i=-D1-A1z1i-B1z1i2-εz1i3.
We next write the homotopy representation of (32) as(33)[H(Z1i,p)H(Z2i,p)]=[Z˙1i+z˙10+pz˙10-pZ2iZ˙2i+z˙2(i-1)+pz˙2(i-1)+p(D1+A1Z1i+B1Z1i2+εZ1i3)]=0.In accordance with the EMHPM approach, the second order expansion Zli=Zli0+pZli1+p2Zli2+⋯ can be substituted into (33), with the initial approximation given as zli0=zli(0)=zl(i-1)(ti-1) for l=1,2. Then, by grouping the terms with the same p order, we obtain the following set of linear equations:
(34)p0:Z˙1i0=0;Z1i0(0)=c1=z1(i-1)(ti-1)Z˙2i0=0;Z2i0(0)=c2=z2(i-1)(ti-1)p1:Z˙1i1=Z2i0Z˙2i1=-D1-A1Z1i0-B1Z1i02-εZ1i03p2:Z˙1i2=Z2i1Z˙2i2=-A1Z1i1-2B1Z1i0Z1i1-3εZ1i02Z1i1.
By considering the initial conditions Z(1,2)ik(0)=0, for k=1,2,3,…, then we get that
(35)Z1i0=c1Z2i0=c2Z1i1=c2TZ2i1=-T(εc13+B1c12+A1c1+D1)Z1i2=-12T2(εc13+B1c12+A1c1+D1)Z2i2=-12T2c2(3εc12+2B1c1+A1).
Note that in (35) the terms Z(1,2)ik can be written in the following recursive form (for k>0):(36)[Z1ikZ2ik]=Tk[Z2i(k-1)-D1g(k)-A1Z1i(k-1)-B1∑n1=0k-1Z1in1Z1i(k-1-n1)-ε∑n2=0k-1∑n1=0n2Z1in1Z1i(n2-n1)Z1i(k-1-n2)].Here, the zeroth order approximation is equal to the initial condition. It is important to note that g(k) is a step function with value of unity for k=1 and zero otherwise.In order to evaluate the accuracy of the EMHPM approach, the derived approximate solution of (31), which is given by (36), is compared with respect to (a) the numerical integration solution and (b) the exact solution of (31) which was recently derived in [25]. We will next explore the applicability of the EMHPM by considering two different sets of system parameter values and use its computer algorithm listed in Algorithm 1.Case 1.
In this case, we suppose that the system parameter values are given asA1=-1.5, B1=5, ε=1, and D1=-0.75 with the following initial conditions: x(0)=1, x˙(0)=0. Table 1 shows the absolute error values computed from the EMHPM and theode45 numerical integration solutions compared with respect to those obtained from the exact solution of (31) at different subinterval size values and m-order solutions. Here the numerical solutions provided by theode45 subroutine were obtained by decreasing its error tolerance to guarantee its convergence to the exact solution of (31). Also Table 1 exhibits a comparison of the computational time performance of the EMHPM approach with respect to that of the numerical methodode45. As we can see from Table 1, the EMHPM approach increases its precision with higher order approximations which reduce the number of subintervals.
Figure7 shows that our derived EMHPM approximate solution closely follows the exact solution of (31), while the numerical solutionode45, with an error tolerance 1E-3, experiences a shift in its period value as time increases. It should be pointed out that the exact solution of (31) which is based on an elliptic integral of the first kind provides a constant period value through the asymmetrical behavior of the system [25]. It is important to mention that the values of m and Δt, in the EMHPM approach, should be carefully selected to guarantee fast convergence, as showed in Table 1.Table 1
Computation time and absolute error comparisons.
EMHPM
ode45
Δ
t
m-order solution
Absolute error
Computation time (ms)
Error tolerance
Absolute error
Computation time (ms)
0.05
7
0.1604
104
1
E
-
3
83.9929
59
0.10
9
0.1321
74
1
E
-
4
36.5743
72
0.15
16
0.1304
78
1
E
-
5
1.7213
88
0.20
22
0.1874
98
1
E
-
6
0.0499
107Figure 7
Solution of (31) with Δt=200 ms and m=22.Case 2.
In this case, we solve (31) by assuming the following parameter values: A1=-1, B1=-1, ε=0.01, and D1=0.25, with initial conditions x(0)=0.209, x˙(0)=0. As we can see from Figure 8, the numerical solutionode45, with an error tolerance equal to 1E-6, starts to deviate from the exact solution as time increases. However, our EMHPM approximate solution shows good agreement with the exact solution of (31).Figure 8
Solution of (31) with Δt=60 ms and m=27.Next, we will use our EMHPM approach to obtain the time-amplitude approximate solutions of an irrational nonlinear spring oscillator and of an elastomagnetic passive suspension system.


### 4.2. Irrational Nonlinear Spring Oscillator

In nonlinear dynamic systems, there are cases when the amplitude-frequency response curves have jumps due to the nonlinearities of the system. In the hard spring model with cubic nonlinearities, the amplitude-frequency response curves exhibit more than one possible dynamics response [26]. Here, we have found the approximate solution of the irrational nonlinear oscillator showed in Figure 9 by using the EMHPM approach. In this oscillator, two helicoidal springs are attached to both sides of a slider that moves in a perpendicular plane with respect to the fixed side of the springs. This arrangement results in an irrational nonlinear stiffness that depends on the slider position. The motion for this system is described by the following equation:
(37)y¨+βy˙+[1-λ1+y2]y=FLΩ2cos⁡(ωt),
where y is the nondimensional slider position relative to the initial undeformed spring length (y=x/L1), F is the force magnitude, Ω is the driving frequency ratio (driving frequency ω relative to the natural frequency ωn), and λ characterizes the spring amount of stretch.Figure 9
Nonlinear spring oscillator, right: elastomagnetic passive suspension.To find the approximate solution of (37) by using the EMHPM, we first write this (37) in state space form by introducing the change of variables z1i=yi, z2i=y˙i. This yields
(38)z˙1i=z2iz˙2i=FLΩ2cos⁡(Ωωnt)-βz2i-[1-λ1+z1i2]z1i.
We next find the homotopy representation of (38) and follow the EMHPM approach to obtain the first order approximate solution of (37):
(39)yi(T)=yi-1+y˙i-1Ty˙i(T)=y˙i-1+[FLΩ2cos⁡(Ωωnt)-βy˙i-1λ1+yi-12-(1-λ1+yi-12)yi-1]T.
Figure 10 shows both the EMHPM approximate solution (obtained by considering subintervals of the same size Δt=0.01) and the numerical integration solutionode45 for the system parameter values of β=0.06, FL=0.06, Ω=1.1, ωn=1.1 Hz, and λ=2. As we can see from Figure 10, both solutions agree well for the time span showed.Figure 10
Solution of the nonlinear spring oscillator withΔt=10 ms.


### 4.3. Elastomagnetic Passive Suspension System

The dynamic response of an elastomagnetic passive suspension system of one degree of freedom that is showed in Figure9 exhibits nonlinear behavior [27]. This type of suspension system has several engineering applications. In this system, the suspended mass is attached to the base by a linear spring device and is excited by a magnetic levitation force. The corresponding equation of motion of the system is given by
(40)mmξ+ζξ˙+k2(ξ-l0)+Fm+mmg=-mmy¨,
where ξ=x-y is the relative position of the suspended mass mm, ζ is the system damping coefficient, k2 is the spring stiffness, l0 is the unstretched length of the spring, g is the gravity, and Fm represents the magnetic levitation force which can be modeled by the following equation:
(41)Fm=-A2(ξ+B2)n,
where the parameters A2, B2, and n have the experimental values of A2=2.26×10-3 Nm3, B2=18.8 mm, and n=3. The base is assumed to be excited by an oscillatory signal given by the following equation: y=A3cos⁡(Ωt), where A3 is the maximum amplitude of motion. In order to obtain the approximate solution of (40) by the EMHPM approach, we first write (40) in state space form by introducing the following change of variables: z1i=ξi and z2i=ξ˙i. Therefore, (40) can be rewritten in equivalent forms as
(42)z˙1i=z2iz˙2i=A3Ω2cos⁡(Ωt)-g-1mm(ζz2i-Kz1i+A2(z1i+B2)-n+Kl0).
By considering the homotopy representation of (42) and by following the EMHPM approach, the first order approximate solution of (42) can be written as
(43)ξi=ξi-1+ξ˙i-1Tξ˙i=ξ˙i-1-(A3Ω2cos⁡(Ωt)+g+1mm1mm×[ζξ˙i-1+Kξi-1-A2(ξi-1+B2)-n-k2l0]T).
To assess the accuracy of our EMHPM approximate solution (43), we compare in Figure 11 the EMHPM approximate solution with respect to the numerical integration solutionode45, by considering the system parameter values of A3=0.7 mm, Ω=2π Hz, ξ(0)=0 m, ξ˙(0)=0 m/seg, mm=11.4 kg, ζ=100 Ns/m, K=2098 N/m, l0=19 mm, and g=9.81 m/s2. As we can see from Figure 11, once again, the proposed EMHPM approximate solution shows excellent agreement with the numerical integration solution for the time span showed in Figure 11.Figure 11
Solution of (40), Δt=1 ms.In all cases discussed here, we can use different system parameter values and still show that the derived semi-numerical-analytic technique that generalizes the MHPM provides good results when compared to the numerical integration solution of the corresponding equations of motion.


## 4.1. The Helmholtz-Duffing Equation

In this section, we will explore the accuracy of the EMHPM approach when this is applied to obtain the approximate solution of the Helmholtz-Duffing equation that has quadratic and cubic nonlinearities. This equation is used to describe the nonlinear response of some materials in mechanical engineering applications [24]. The differential equation for the Helmholtz-Duffing oscillator is given as
(31)x¨+A1x+B1x2+εx3+D1=0,x(0)=x0,x˙(0)=x˙0=0,
where x is the system displacement, A1, B1, D1, and ε are physical parameters, and x0 and x˙0 are the system initial conditions. By introducing the following change of variables: z1i=xi, z2i=x˙i, we can write (31) in state space form:
(32)z˙1i=z2iz˙2i=-D1-A1z1i-B1z1i2-εz1i3.
We next write the homotopy representation of (32) as(33)[H(Z1i,p)H(Z2i,p)]=[Z˙1i+z˙10+pz˙10-pZ2iZ˙2i+z˙2(i-1)+pz˙2(i-1)+p(D1+A1Z1i+B1Z1i2+εZ1i3)]=0.In accordance with the EMHPM approach, the second order expansion Zli=Zli0+pZli1+p2Zli2+⋯ can be substituted into (33), with the initial approximation given as zli0=zli(0)=zl(i-1)(ti-1) for l=1,2. Then, by grouping the terms with the same p order, we obtain the following set of linear equations:
(34)p0:Z˙1i0=0;Z1i0(0)=c1=z1(i-1)(ti-1)Z˙2i0=0;Z2i0(0)=c2=z2(i-1)(ti-1)p1:Z˙1i1=Z2i0Z˙2i1=-D1-A1Z1i0-B1Z1i02-εZ1i03p2:Z˙1i2=Z2i1Z˙2i2=-A1Z1i1-2B1Z1i0Z1i1-3εZ1i02Z1i1.
By considering the initial conditions Z(1,2)ik(0)=0, for k=1,2,3,…, then we get that
(35)Z1i0=c1Z2i0=c2Z1i1=c2TZ2i1=-T(εc13+B1c12+A1c1+D1)Z1i2=-12T2(εc13+B1c12+A1c1+D1)Z2i2=-12T2c2(3εc12+2B1c1+A1).
Note that in (35) the terms Z(1,2)ik can be written in the following recursive form (for k>0):(36)[Z1ikZ2ik]=Tk[Z2i(k-1)-D1g(k)-A1Z1i(k-1)-B1∑n1=0k-1Z1in1Z1i(k-1-n1)-ε∑n2=0k-1∑n1=0n2Z1in1Z1i(n2-n1)Z1i(k-1-n2)].Here, the zeroth order approximation is equal to the initial condition. It is important to note that g(k) is a step function with value of unity for k=1 and zero otherwise.In order to evaluate the accuracy of the EMHPM approach, the derived approximate solution of (31), which is given by (36), is compared with respect to (a) the numerical integration solution and (b) the exact solution of (31) which was recently derived in [25]. We will next explore the applicability of the EMHPM by considering two different sets of system parameter values and use its computer algorithm listed in Algorithm 1.Case 1.
In this case, we suppose that the system parameter values are given asA1=-1.5, B1=5, ε=1, and D1=-0.75 with the following initial conditions: x(0)=1, x˙(0)=0. Table 1 shows the absolute error values computed from the EMHPM and theode45 numerical integration solutions compared with respect to those obtained from the exact solution of (31) at different subinterval size values and m-order solutions. Here the numerical solutions provided by theode45 subroutine were obtained by decreasing its error tolerance to guarantee its convergence to the exact solution of (31). Also Table 1 exhibits a comparison of the computational time performance of the EMHPM approach with respect to that of the numerical methodode45. As we can see from Table 1, the EMHPM approach increases its precision with higher order approximations which reduce the number of subintervals.
Figure7 shows that our derived EMHPM approximate solution closely follows the exact solution of (31), while the numerical solutionode45, with an error tolerance 1E-3, experiences a shift in its period value as time increases. It should be pointed out that the exact solution of (31) which is based on an elliptic integral of the first kind provides a constant period value through the asymmetrical behavior of the system [25]. It is important to mention that the values of m and Δt, in the EMHPM approach, should be carefully selected to guarantee fast convergence, as showed in Table 1.Table 1
Computation time and absolute error comparisons.
EMHPM
ode45
Δ
t
m-order solution
Absolute error
Computation time (ms)
Error tolerance
Absolute error
Computation time (ms)
0.05
7
0.1604
104
1
E
-
3
83.9929
59
0.10
9
0.1321
74
1
E
-
4
36.5743
72
0.15
16
0.1304
78
1
E
-
5
1.7213
88
0.20
22
0.1874
98
1
E
-
6
0.0499
107Figure 7
Solution of (31) with Δt=200 ms and m=22.Case 2.
In this case, we solve (31) by assuming the following parameter values: A1=-1, B1=-1, ε=0.01, and D1=0.25, with initial conditions x(0)=0.209, x˙(0)=0. As we can see from Figure 8, the numerical solutionode45, with an error tolerance equal to 1E-6, starts to deviate from the exact solution as time increases. However, our EMHPM approximate solution shows good agreement with the exact solution of (31).Figure 8
Solution of (31) with Δt=60 ms and m=27.Next, we will use our EMHPM approach to obtain the time-amplitude approximate solutions of an irrational nonlinear spring oscillator and of an elastomagnetic passive suspension system.


## 4.2. Irrational Nonlinear Spring Oscillator

In nonlinear dynamic systems, there are cases when the amplitude-frequency response curves have jumps due to the nonlinearities of the system. In the hard spring model with cubic nonlinearities, the amplitude-frequency response curves exhibit more than one possible dynamics response [26]. Here, we have found the approximate solution of the irrational nonlinear oscillator showed in Figure 9 by using the EMHPM approach. In this oscillator, two helicoidal springs are attached to both sides of a slider that moves in a perpendicular plane with respect to the fixed side of the springs. This arrangement results in an irrational nonlinear stiffness that depends on the slider position. The motion for this system is described by the following equation:
(37)y¨+βy˙+[1-λ1+y2]y=FLΩ2cos⁡(ωt),
where y is the nondimensional slider position relative to the initial undeformed spring length (y=x/L1), F is the force magnitude, Ω is the driving frequency ratio (driving frequency ω relative to the natural frequency ωn), and λ characterizes the spring amount of stretch.Figure 9
Nonlinear spring oscillator, right: elastomagnetic passive suspension.To find the approximate solution of (37) by using the EMHPM, we first write this (37) in state space form by introducing the change of variables z1i=yi, z2i=y˙i. This yields
(38)z˙1i=z2iz˙2i=FLΩ2cos⁡(Ωωnt)-βz2i-[1-λ1+z1i2]z1i.
We next find the homotopy representation of (38) and follow the EMHPM approach to obtain the first order approximate solution of (37):
(39)yi(T)=yi-1+y˙i-1Ty˙i(T)=y˙i-1+[FLΩ2cos⁡(Ωωnt)-βy˙i-1λ1+yi-12-(1-λ1+yi-12)yi-1]T.
Figure 10 shows both the EMHPM approximate solution (obtained by considering subintervals of the same size Δt=0.01) and the numerical integration solutionode45 for the system parameter values of β=0.06, FL=0.06, Ω=1.1, ωn=1.1 Hz, and λ=2. As we can see from Figure 10, both solutions agree well for the time span showed.Figure 10
Solution of the nonlinear spring oscillator withΔt=10 ms.


## 4.3. Elastomagnetic Passive Suspension System

The dynamic response of an elastomagnetic passive suspension system of one degree of freedom that is showed in Figure9 exhibits nonlinear behavior [27]. This type of suspension system has several engineering applications. In this system, the suspended mass is attached to the base by a linear spring device and is excited by a magnetic levitation force. The corresponding equation of motion of the system is given by
(40)mmξ+ζξ˙+k2(ξ-l0)+Fm+mmg=-mmy¨,
where ξ=x-y is the relative position of the suspended mass mm, ζ is the system damping coefficient, k2 is the spring stiffness, l0 is the unstretched length of the spring, g is the gravity, and Fm represents the magnetic levitation force which can be modeled by the following equation:
(41)Fm=-A2(ξ+B2)n,
where the parameters A2, B2, and n have the experimental values of A2=2.26×10-3 Nm3, B2=18.8 mm, and n=3. The base is assumed to be excited by an oscillatory signal given by the following equation: y=A3cos⁡(Ωt), where A3 is the maximum amplitude of motion. In order to obtain the approximate solution of (40) by the EMHPM approach, we first write (40) in state space form by introducing the following change of variables: z1i=ξi and z2i=ξ˙i. Therefore, (40) can be rewritten in equivalent forms as
(42)z˙1i=z2iz˙2i=A3Ω2cos⁡(Ωt)-g-1mm(ζz2i-Kz1i+A2(z1i+B2)-n+Kl0).
By considering the homotopy representation of (42) and by following the EMHPM approach, the first order approximate solution of (42) can be written as
(43)ξi=ξi-1+ξ˙i-1Tξ˙i=ξ˙i-1-(A3Ω2cos⁡(Ωt)+g+1mm1mm×[ζξ˙i-1+Kξi-1-A2(ξi-1+B2)-n-k2l0]T).
To assess the accuracy of our EMHPM approximate solution (43), we compare in Figure 11 the EMHPM approximate solution with respect to the numerical integration solutionode45, by considering the system parameter values of A3=0.7 mm, Ω=2π Hz, ξ(0)=0 m, ξ˙(0)=0 m/seg, mm=11.4 kg, ζ=100 Ns/m, K=2098 N/m, l0=19 mm, and g=9.81 m/s2. As we can see from Figure 11, once again, the proposed EMHPM approximate solution shows excellent agreement with the numerical integration solution for the time span showed in Figure 11.Figure 11
Solution of (40), Δt=1 ms.In all cases discussed here, we can use different system parameter values and still show that the derived semi-numerical-analytic technique that generalizes the MHPM provides good results when compared to the numerical integration solution of the corresponding equations of motion.


## 5. Conclusions

In this work, we have modified the HPM and introduced a new numerical-analytic approach that is based on a linear operator defined asL(ui)=(d/dT)ui and time subintervals that do not need to be equally spaced to solve strongly nonlinear ordinary differential equations with variable coefficients. To assess the convergence of this EMHPM approach, we have used the convergence control parameter h value and derived the approximate solutions of some nonlinear differential equations to show that our EMHPM approach converges not only to the numerical integration solutions, computed from the classical Runge-Kutta fourth order algorithm (ode45), but also to the corresponding exact analytical solutions. Furthermore, we have proved that if the time interval and the m-order approximation of our EMHPM approach are properly chosen, then its computational time is less than the time spent by the standard numerical integration solutions provided by theode45 MATLAB function, as showed in Table 1. Future work includes the adaptation of EMHPM for the solution of delay and fractional nonlinear differential equations which are common in several physical and engineering applications. The results of this work will be presented in future publications.



---

*Source: 486509-2014-04-17.xml*