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 parameter h 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.

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 [46] 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: that has the following boundary condition: where is a differential operator, is a boundary operator, is a known analytic function, and Γ is the boundary in domain Ω.

The operator can be divided into two parts: and , where represents the linear terms while is related to nonlinear terms. Equation (1) can therefore be written as follows: Applying the principles of topology, the following homotopy is established: , satisfying the following conditions: where is an embedding parameter and is the initial approximation to (1) that satisfies the boundary condition. It is evident from (4) that setting from zero to unity makes change from to , that is, by letting In topology, this is called deformation and, in the same context, and are called homotopic. He in [4] assumed that the solution of (4) can be expressed with a series of power terms of the form Thus, by considering , the approximate solution of (1) yields In order to test the precision achieved with the election of the linear operator, we proceed to analyze the following nonlinear equation [20]: where is the independent variable, is the running time, and 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 , and therefore, . Thus, the homotopic representation of (8) is given as By substituting the second order expansion into (9) and by grouping the terms with the same order, we obtain the following set of equations: whose solutions are given by The approximate solution of (8), given by (11), corresponds to the second order solution in [20] when .

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 the enhanced homotopy perturbation method. By using the EHPM, (8) is rewritten by considering a stable linear operator of the form , while the nonlinear operator is defined as . These operator definitions provide the following expression: 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 order, we obtain that Integration of (13) provides the following expressions: Notice that the second order solution of (14) agrees well with the solutions obtained in [20] when .

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 By using the linear operator , the homotopic representation of (15) can be expressed as In order to get the first approximation of , the term has to be included in (16). Substituting the third order expansion of (15) into (16) and by regrouping the terms with the same order, we get that Thus, the fourth order solution of (8) can be obtained by integrating (17); this yields 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 function ode45, while the He and the Odibat solutions tend to diverge at increasing values of . 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.

function 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*x02*(exp(t)-exp(3*t));
y2=.25*x03*(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-x02)*t.*exp(-t);
y2=(2*x0-3*x02+x03).*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;

3. Description of the Enhanced Multistage Homotopy Perturbation Method

The proposed enhanced multistage homotopy perturbation method (EMHPM) is an algorithm that utilizes the HPM solution in time subintervals, based on the following transformation: , where is the approximate solution of the th subinterval that satisfies the initial solution . The new time shifted variable satisfies the following condition: . Also, we assume that the proposed solution can be given by , where 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 operator can be chosen as , where the proposed initial solution is equal to the initial condition ; that is, . For simplicity in the notation, we write .(2)Since the homotopy is defined in the th subinterval, the following relation holds: , which satisfies the condition , where represents the current time. Therefore, the approximation of th order can be obtained by integrating with respect to , while the terms related to the independent variable are assumed constant in the th subinterval.Therefore, the approximate solution of order for the differential equation can be written as Note that the solution is valid only in the th subinterval . In order to obtain the solution in the th subinterval, we consider that the relation holds. Therefore, the approximate solution of at the time is given by In summary, each solution for a given subinterval is in turn divided into subintervals, which do not necessarily have to be equally spaced: . Finally, the approximate solution for is obtained by coupling the solutions .

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 and the nonlinear operator . 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 order By solving (22), we get that Note that the solutions given by (23) can be generalized in a recursive expression of the form where is an integer number that is bigger than zero.

Figure 3 shows the comparison between the approximate solution of (8) obtained by using the EMHPM (utilizing the third approximation and time subintervals of the same size ) 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).

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 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 . 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 Utilizing the third order expansion of (25) yields Notice that if we take the value of in (26) then this equation becomes similar to (23).

Figure 4 shows the absolute error between the numerical solution obtained from the ode45 subroutine in MATLAB and the solution obtained from the EMHAM as a function of the convergent factor . Notice from Figure 4 that the minimum absolute error occurs at . Also, at this value of the convergent control parameter , the solutions provided by the EMHPM and the EMHAM have the same absolute error.

We next consider the second order nonlinear ordinary differential equation used in [22] given as which has an exact solution of the form with an exact homotopic solution (HAM) given as 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: where is the amplitude and is the first derivate of -order. Figure 5 shows the influence of the subinterval size as well as that of the -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 . That implies that the convergence could be reaching by increasing the order solution and/or by decreasing the subinterval size.

Figure 6 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 by considering the subinterval size of with . We can see from Figure 6 that the absolute error has a minimum value at . Thus, our proposed EMHPM approach, which is a special case of the EMHAM at , overlaps the exact solution given by (28).

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 where is the system displacement, , , , and are physical parameters, and and are the system initial conditions. By introducing the following change of variables: , , we can write (31) in state space form: We next write the homotopy representation of (32) asIn accordance with the EMHPM approach, the second order expansion can be substituted into (33), with the initial approximation given as for . Then, by grouping the terms with the same order, we obtain the following set of linear equations: By considering the initial conditions , for , then we get that Note that in (35) the terms can be written in the following recursive form (for ):Here, the zeroth order approximation is equal to the initial condition. It is important to note that is a step function with value of unity for 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 as , , , and with the following initial conditions: , . Table 1 shows the absolute error values computed from the EMHPM and the ode45 numerical integration solutions compared with respect to those obtained from the exact solution of (31) at different subinterval size values and -order solutions. Here the numerical solutions provided by the ode45 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 method ode45. As we can see from Table 1, the EMHPM approach increases its precision with higher order approximations which reduce the number of subintervals.
Figure 7 shows that our derived EMHPM approximate solution closely follows the exact solution of (31), while the numerical solution ode45, 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 and , in the EMHPM approach, should be carefully selected to guarantee fast convergence, as showed in Table 1.

Case 2. In this case, we solve (31) by assuming the following parameter values: , , , and , with initial conditions , . As we can see from Figure 8, the numerical solution ode45, 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).

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: where is the nondimensional slider position relative to the initial undeformed spring length , is the force magnitude, is the driving frequency ratio (driving frequency relative to the natural frequency ), and characterizes the spring amount of stretch.

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 , . This yields We next find the homotopy representation of (38) and follow the EMHPM approach to obtain the first order approximate solution of (37): Figure 10 shows both the EMHPM approximate solution (obtained by considering subintervals of the same size ) and the numerical integration solution ode45 for the system parameter values of , , ,  Hz, and . As we can see from Figure 10, both solutions agree well for the time span showed.

4.3. Elastomagnetic Passive Suspension System

The dynamic response of an elastomagnetic passive suspension system of one degree of freedom that is showed in Figure 9 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 where is the relative position of the suspended mass , is the system damping coefficient, is the spring stiffness, is the unstretched length of the spring, is the gravity, and represents the magnetic levitation force which can be modeled by the following equation: where the parameters , , and have the experimental values of  Nm3,  mm, and . The base is assumed to be excited by an oscillatory signal given by the following equation: , where 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: and . Therefore, (40) can be rewritten in equivalent forms as By considering the homotopy representation of (42) and by following the EMHPM approach, the first order approximate solution of (42) can be written as 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 solution ode45, by considering the system parameter values of  mm,  Hz,  m,  m/seg,  kg,  Ns/m,  N/m,  mm, and  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.

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 as 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 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 -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 the ode45 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.

Conflict of Interests

The authors declare that they have no conflict of interests with any mentioned entities in the paper.

Acknowledgments

This work was funded by Tecnológico de Monterrey, Campus Monterrey, through the Research Chair in Nanomaterials for Medical Devices and the Research Chair in Intelligent Machines. Additional support was provided from Consejo Nacional de Ciencia y Tecnología (Conacyt), México.