Abstract and Applied Analysis
Volume 2015 (2015), Article ID 382475, 12 pages
http://dx.doi.org/10.1155/2015/382475

Research Article

Approximate Solutions of Delay Differential Equations with Constant and Variable Coefficients by the Enhanced Multistage Homotopy Perturbation Method

D. Olvera,¹ A. Elías-Zúñiga,¹ L. N. López de Lacalle,² and C. A. Rodríguez¹

¹Center for Innovation in Design and Technology, Tecnológico de Monterrey, Campus Monterrey, E. Garza Sada 2501, 64849 Monterrey, NL, Mexico
²Department of Mechanical Engineering, University of the Basque Country, Alameda de Urquijo s/n, Bilbao, 48013 Bizkaia, Spain

Received 26 December 2013; Revised 21 July 2014; Accepted 16 August 2014

Academic Editor: Zhichun Yang

Copyright © 2015 D. Olvera et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We expand the application of the enhanced multistage homotopy perturbation method (EMHPM) to solve delay differential equations (DDEs) with constant and variable coefficients. This EMHPM is based on a sequence of subintervals that provide approximate solutions that require less CPU time than those computed from the dde23 MATLAB numerical integration algorithm solutions. To address the accuracy of our proposed approach, we examine the solutions of several DDEs having constant and variable coefficients, finding predictions with a good match relative to the corresponding numerical integration solutions.

1. Introduction

Delayed differential equations (DDEs) are used to describe many physical phenomena of interest in biology, medicine, chemistry, physics, engineering, and economics, among others. Since the introduction of the first delayed models, many publications have appeared as summarizing theorems and homotopy methods of solution that deal with the stability properties of delayed systems (see [1–3] and references cited there in).

For instance, Shakeri and Dehghan introduced an approach to find the solution of delay differential equations by means of the homotopy perturbation technique (HPM) with results that agree well with exact solutions [1]. Wu in [2] used the homotopy analysis method to obtain the approximate solution of a strong nonlinear ENSO delayed oscillator model that provides good agreement when compared to its exact solution under the condition of . Alomari and coworkers in [3] developed an algorithm to obtain approximate analytical solutions for DDEs by using the homotopy analysis method (HAM) and the modified homotopy analysis method (MHAM). They used their derived method to obtain the approximate solution of various linear and nonlinear DDEs with numerical predictions that agree well with the numerical integration solutions, and they also proved that their derived solutions converge to the exact ones. By applying the homotopy perturbation method (HPM), Biazar and Behzad found approximate solutions of neutral differential equations with proportional delays which describe well their corresponding numerical integration solutions [4]. Recently, Anakira and co-workers in [5] extended the applicability of the so called optimal homotopy asymptotic method (OHAM) that does not depend on small or large parameters, to find the approximate analytic solution of DDEs. They used their proposed approach to compare the derived approximate solutions of several DDEs with their exact analytical solutions with predictions that compare well with the exact ones.

On the other hand, Insperger and Stépán in [6] used the semidiscretization method to determine the stability lobes of DDEs that model the dynamics of cutting machine operations. Based on the properties of the Chebyshev polynomials, Butcher and coworkers in [7] developed a methodology to obtain the stability lobes of milling machine operations and they proved that this technique is faster than that of the full and the semidiscretization methods since these solution techniques approximate the original DDEs by a series of ODEs [8].

Here in this paper, we develop a generalized procedure to solve linear and nonlinear DDEs by introducing some modifications to the multistage homotopy perturbation method (MHPM) derived by Hashim and Chowdhury to obtain approximate solutions of ordinary differential equations [9]. The proposed enhanced multistage homotopy perturbation method (EMHPM) is based on a sequence of subintervals that allow us to find more accurate approximated solutions under a numerical-analytical procedure that requires less CPU time when compared to the numerical integration solutions provided by the MATLAB dde23 algorithm written by Shampine and Thompson in [10]. The EMHPM is based on a homotopy function that could be divided into a linear operator and a nonlinear operator to satisfy its assumed initial solution. This split of the homotopy function allows us to modify the nonlinear operator to guarantee, by using the enhanced homotopy perturbation method, the stability of the proposed approximate solutions of nonlinear differential equations [11].

To clarify our proposed method, we briefly review in Section 2 some basic concepts of the homotopy perturbation method, and, then in Section 3, we introduce the EMHPM to solve DDEs. The difference between the HPM and the EMHPM is discussed in Section 4 by addressing the approximate solutions of a nonlinear delayed differential equation with variable coefficients. Finally, the general solution of two DDEs that describe the dynamics of two engineering problems, by using the EMHPM, is discussed in Section 5.

2. Homotopy Perturbation Method

The homotopy perturbation method (HPM) is a coupling of the traditional perturbation method and homotopy in topology which eliminates the limitation of the small parameter assumed in the perturbation methods [12]. Under this approach, a nonlinear problem can be transformed into an infinite number of simple problems without the restriction of having small nonlinear parameter values. This homotopy perturbation method takes the main advantages of traditional perturbation methods together with homotopy analysis [13–15].

To illustrate the basic ideas of the HPM, let us consider the following nonlinear differential equation:

with boundary conditions

where is a general differential operator, is a boundary operator, is a known analytic function, and is the boundary of the domain .

The operator can generally be divided into two parts: and , where involves the linear terms and the nonlinear ones. Equation (1) therefore can be rewritten as follows:

By the homotopy perturbation technique, we construct a homotopy that satisfies

where is an embedding parameter and is an initial approximation of (1) which satisfies the boundary conditions (2). Thus, from (4), we have

The changing process of from zero to unity is just that of from to . In topology, this is called deformation, and and are called homotopic.

He in [12] uses the embedding parameter as the small parameter and assumed that the solution of (4) can be written as a power series of in the form

By setting , He obtained the approximate solution of (1) as

Then, this method was applied to obtain the approximate solution of some nonlinear ordinary differential equations valid not only for small, but also for large nonlinear parameter values.

We next will introduce an approach based on homotopy methods, to obtain the solution of DDEs with constant and variable coefficients.

3. The EMHPM Methodology to Solve DDEs

The HPM is an asymptotic method that depends on the auxiliary linear operator form and the initial guess of the initial conditions. Therefore, the convergence of the approximate solution cannot be guaranteed in some cases [16]. Hashim and Chowdhury showed in [9] that the solutions obtained by the standard HPM were not valid for large time span unless more terms are calculated. Thus, they proposed a multistage homotopy perturbation method (MHPM) which treated the HPM algorithm in a sequence of subintervals in an attempt to improve the accuracy of the approximate solutions of linear and nonlinear ordinary differential equations (ODEs).

However, when the MHPM is applied to obtain the approximate solutions of ODEs which contain coefficients as a function of time, this method cannot provide accurate solutions when . In this work, we introduce some modifications to the MHPM and focus on the derivation of approximate solutions of DDEs equations with variable coefficient terms. This new approach is based on the enhanced multistage homotopy perturbation method (EMHPM) introduced in [17] to obtain the solution of nonlinear ordinary differential equations.

The EMHPM is an algorithm which approximates the HPM solution by subintervals, utilizing the following transformation rule: , where satisfies the initial condition , is a shifted time scale used to determine the approximate solution in each subinterval, and represents the approximate solution in the th subinterval. In this case, the initial suggested solution in the th subinterval is given by , where represents the time at the end of the previous subinterval (i.e., the value of the approximate solution at the end of the previous subinterval represents the initial conditions of the next subinterval under consideration).

To apply the homotopy technique to solve delay differential equations, we also assume the following.

(1)The linear operator represents , where the assumed approximate solution is set equal to the initial condition ; that is, . To simplify the notation, we let .(2)The transformation on holds in the homotopy -subinterval. Thus, higher order equations are integrated with respect to , while the terms related to the independent variable are assumed to remain constant.

Therefore, we may conclude that the order approximate solution, by applying the EMHPM, can be written as

where the solution is valid only in the th subinterval . Hence, the solution on the th subinterval can be written as

with initial condition , and . Thus, the approximate solution of at the time is given by

In summary, the solution for an open-closed interval (] is divided into subintervals that, in general, are not equally spaced: . Thus, the approximated solution of for the span time interval is obtained by coupling the solutions.

4. Approximate Solutions of Some DDEs by Applying the EMHPM

In this section we focus on the solution of DDEs with constant and variable coefficients and examine the applicability of the EMHPM to find the corresponding approximate solutions.

4.1. Delay Differential Equations with Constant Coefficients

First, let us consider the simplest DDE of the form

with initial condition . Here, the independent variable is a scalar , the dot stands for differentiation with respect to time , and is the time delay. To evaluate (11) on , the term must represent a known function on . For instance, if , the solution of (11) can be obtained in the interval by assuming an initial function that satisfies the initial condition. By using this solution, it becomes possible to obtain the solution of (11) in the next th interval , , where is an integer number that can be chosen as . With this approach, we can apply the HPM to find the solution of (11) by assuming that the previous delayed function is ; thus the solution for the first interval is given by , valid on . In terms of (4), we now construct the homotopy of (11):

We next substitute the first order expansion in (12) and balance the terms with identical power of to obtain the following set of linear differential equations:

Integration of (13) yields

Hence, the first order solution of (12) is given by

Notice that (15) represents the exact solution of (11) on the first interval. By following the same procedure, it is easy to show that the exact solution of (11), for the second and third intervals, is given, respectively, as

Figure 1 shows the exact solution of (11) obtained by coupling at each interval the solution obtained by following HPM procedure for .

Figure 1: Exact solution of (11) obtained by using the HPM and .

It is easy to show that the solution of (11) by the EMHPM coincides with the solution obtained by using the HPM since (11) is a delay differential equation with constant coefficients.

4.2. Delay Differential Equations with Variable Coefficients

We next show how the EMHPM approach can be applied to obtain the approximate solution of nonlinear delay differential equation with variable coefficients. In this case, we obtain the approximate solutions of a DDE of the form

in which the solution holds on . In order to find the solution in the interval , we assume that the homotopy representation of (17) can be given as

Notice that the variable depends on the time for which . If we now substitute the second order expansion in (18), and, after balancing the terms, we get that

Equations (19) have the following solutions:

Thus, the approximate solution of (17) by using the EMHPM is given by

In this case, the exact solution of is unknown. To obtain , we compute again the approximate solution of by applying our EMHPM and the value of the delayed time is assumed to remain constant in each subinterval. To determine , we next use the homotopy representation of (17) for the interval :

Substituting the second order expansion in (22), we get

Note that (20) and (23) provide approximate solutions to (17) but evaluated at different interval time delays. To find the third order approximate solution of (17), we can use a homotopy of the form:

Then, by using our EMPHM approach, we have that

Notice from (25) that the th order approximate solution of (17) can be written as

where , when and zero otherwise.

Figure 2 shows the approximate solution of (17) obtained by using the EMHPM approach compared to its numerical integration solution by using the dde23 MATLAB subroutine program. This case assumes two different initial solutions of the form , , and a time subintervals . We can see from Figure 2 that both simulations agree well for the time span showed.

Figure 2: EMHPM and dde23 solution of (17).

To further assess the applicability of our proposed EMHPM approach to high order delay differential equations, we will next describe a methodology to obtain the approximate solutions of well-known high order delay differential equations by generalizing our EMHPM approach.

5. Generalized Solution of Linear DDEs by the EMHPM Approach

Let us consider an -dimensional delay differential equation of the form

where , , is the state vector, and is the time delay. By following our EMHPM procedure, we can write (27) in equivalent form as

where denotes the order solution of (27) in the th subinterval that satisfies the initial conditions and and represent the values of the periodic coefficients at the time . In order to approximate the delayed term in (28), the period is discretized in points equally spaced as shown in Figure 3. Here, we assume that the function in the delay subinterval is approximated by a constant value

as shown in Figure 3. By following the homotopy perturbation technique, we can write the homotopy representation of (28) as

Substituting the order expansion in (30) and by assuming an initial approximation of the form , we get, after applying the proposed EMHPM approach, the following set of first order linear delay differential equations:

By solving (31), we get

Equations (32) can be written as

where and for and , otherwise. Thus, the solution of (27) is obtained by adding the approximate solutions:

Notice, however, that solution (34) may be further improved by using a first order polynomial representation of as shown in Figure 4. Then, the function in the delay subinterval takes the form

Substituting (35) into (28) gives

We next assume that the homotopy representation of (36) is given as

Substituting the order expansion    in (37) and assuming that the initial approximation is given by , we get

By solving (38) and by following the EMHPM procedure, we get

Here, the recursive form of is written as

where and

Thus, the approximate solution of (27) by the EMHPM can be obtained by substituting (40) into (34).

Figure 3: Schematic of the zeroth order polynomial used to fit the approximate EMHPM solution.

Figure 4: Schematic EMHPM solution using first polynomial to approximate delay subinterval.

In the next section, we will apply our EMHPM procedure to obtain the solution of two second order delay differential equations: (a) the damped Mathieu equation with time delay, and (b) the well-known delay differential equation that describes the dynamics in one degree-of-freedom milling machine operations.

5.1. Solution of the Damped Mathieu Equation with Time Delay

In order to assess the accuracy of our EMHPM approach, we first obtain the solution of the damped Mathieu differential equation with time delay that combines the effect of parametric excitation and damping. This equation is described by the following equation:

where , , , , and are system parameters whose value depends on the physics of the system. The approximate solution of (42) obtained by using the semidiscretization method is widely discussed in [18, 19]. Here, we focus our attention on applying the EMHPM to find the approximate solution of (42) and we also assess the accuracy of the derived solution by comparing it with the corresponding numerical integration solution of (42).

By following the EMHPM procedure, we first write (42) in the following equivalent form:

where denotes the order solution of (43) in the th subinterval that satisfies the following initial conditions: and . The space state form representation of (43) is given by

where

and is a time periodic term. The EMHPM approximate solution of (44) is illustrated in Figure 5 where we have assumed an unstable system behavior for which , , , , and . See [20]. As we can see from Figure 5, our approximate EMHPM solution to (42) is compared with its numerical integration solution obtained from dde23 MATLAB algorithm for the time interval of , by assuming that with the following initial values: and .

Figure 5: Numerical solutions of the damped Mathieu equation with time delay by dde23, the zeroth EMHPM, and the first EMHPM with and .

It can be seen from Figure 5 that in the interval [0, ] both the zeroth and the first order solutions are the same since the delay subintervals are constant. See Figure 3. However, in the next interval [, ] it is clear that the first order EMPHM solution provides a better approximation on the delay subinterval. The computation total time to calculate the solutions in the MATLAB code is listed in Table 1. The order and the discretized time intervals in the EMPHM approach are chosen to guarantee the convergence of our approximate solution to the exact one. To provide a full understanding of how the solution is computed by the EMHPM approach, we attached in Algorithm 1 the corresponding MATLAB code.

Table 1: Computer time needed to solve the damped Mathieu equation with time delay. The order solution of the EMPHM approach is chosen to guarantee the convergence of its approximate solution.

dde23 ms
EMHPM

Zeroth ms
First ms

19
15
5
4
5

20
4
6
6

40
3
12
12

60
2
17
18

60
5
18
19

60
10
20
21

function dde23_ddeEM_mathieu_paper

disp(‘Mathieu equation solution’)

% Solution by EMHPM with zeroth and first order solution

eps=1; kappa=0.2; T=2pi; tau=T; b=−1; delta=3; % Mathieu Parameters

pntDelay=1; N=50; m_ord=5; dt=tau/(N−1); ktau=2; tspan=0,ktauT; % EMHPM Parameters

%% Solution

% Solution by dde23

tdde=linspace(0,ktautau,ktau(N−1)+1);

dde=@(t,y,z) mathieu_dde(t,y,z,kappa,delta,eps,b,tau);

te_dde=tic;sol = dde23(dde,tau,@history,tspan); toc(te_dde); xdde=deval(sol,tdde);

% Solution by zeroth EMHPM

dde_emhpm_fun=@(t,c0,tt,zm,zn,tau,N) mathieu_zeroth(t,c0,tt,zm,zn,tau,N,m_ord,

kappa,delta,eps,b,T);

tic;t0,z0=ddeEMHPM(dde_emhpm_fun,tspan,@history,tau,dt,pntDelay); toc

% Solution by First EMHPM

dde_emhpm_fun=@(t,c0,tt,zm,zn,tau,N) mathieu_first(t,c0,tt,zm,zn,tau,N,m_ord,kappa,delta,eps,b,T);

tic;t1,z1=ddeEMHPM(dde_emhpm_fun,tspan,@history,tau,dt,pntDelay); toc

% Plot results

ind0=1:2:ktau(N−1); ind1=2:2:ktau*(N−1);

Parent1=figure (1);

axes1 = axes(‘Parent’,Parent1,‘FontSize’,12,‘FontName’,‘Times New Roman’);

box(axes1,‘on’); hold(axes1,‘all’);

plot(tdde,xdde(1,:),‘Parent’,axes1,‘LineWidth’,2,‘Color’,0.502 0.502 0.502,‘DisplayName’,‘Numerical

dd23’);

plot(t0(ind0),z0(ind0,1),‘MarkerSize’,5,‘Marker’,‘o’,‘LineStyle’, ‘none’,‘DisplayName’,‘Zeroth

EMHPM’,‘Color’,0 0 0);

plot(t1(ind1),z1(ind1,1),‘MarkerSize’,7,‘Marker’,‘x’,‘LineStyle’,‘none’,‘DisplayName’,‘First

EMHPM’,‘Color’,0 0 0);

xlabel(‘∖itt’,‘FontSize’,12,‘FontName’,‘Times New Roman’);

ylabel(‘∖itx’,‘FontSize’,12,‘FontName’,‘Times New Roman’);

end

%% Mathieu definitions

function dydt = mathieu_dde(t,y,z,kapa,dlt,eps,b,T)

dydt = y(2)

    -kapay(2)−(dlt+epscos(2pit/T))y(1)+bz(1);

end

function Z = mathieu_zeroth(t,c0,tt,zm,zn,tau,N,m,kpa,dlt,eps,b,T)

Z=c0(1),c0(2); alf=dlt+epscos(2pi/Ttt);

for ik=1:m

    Z(ik+1,1)=Z(ik,2)t/ik;

    Z(ik+1,2)=−kpaZ(ik,2)−alfZ(ik,1);

    if ik==1, Z(ik+1,2)=Z(ik+1,2)+bzm(1); end

    Z(ik+1,2)=Z(ik+1,2)t/ik;

end

Z=sum(Z);

end

function Z = mathieu_first(t,c0,tt,zm,zn,tau,N,m,kpa,dlt,eps,b,T)

alf=dlt+epscos(2pi/Ttt); Z=c0(1),c0(2); Z_=0,0;

for ik=1:m

    Z(ik+1,:)=Z(ik,2)t/ik, −kpaZ(ik,2)−alfZ(ik,1);

    Z_(ik+1,:)=Z_(ik,2)t/ik, −kpaZ_(ik,2)−alfZ_(ik,1);

    if ik==1,

        Z(ik+1,2)=Z(ik+1,2)+bzm(1);

        Z_(ik+1,2)=Z_(ik+1,2)+b(N−1)/tau(zn(1)−zm(1))t;

    end

    Z(ik+1,2)=Z(ik+1,2)t/ik;

    Z_(ik+1,2)=Z_(ik+1,2)t/(ik+1);

end

Z=sum(Z)+sum(Z_);

end

function out=history(t)

out=1E−3+0t 0+0t;

end

%% EMHPM algorithm for ODE solutions

function t,z= odeEMHPM(nde,tspan,z0,Deltat,pnts)

% ode solver by Enhanced Multistage Homotopy Perturbation Method

tini=tspan(1); tfin=tspan(end); tstart=tini;

tini=tini−tstart; tfin=tfin−tstart; % shifted time set to zero

% Handle errors

if tini == tfin

    error(‘The ending and starting time values must be different.’);

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(‘tspan entries must be strictly sorted.’);

end

incT=Deltat/pnts;

if incT<=0

    error (‘Increasing must be greater than zero.’)

end

z(1,:)=z0′; t(1)=tini; iteT=2; % set initial values

t(iteT)=tini+incTtdir;

while tdirt(iteT)<tdir*(tfin+tdirincT)

    P_act=ceil(.99999(t(iteT)−tini)tdir/Deltat); % count sub-intervals

    c=z((Deltat/incT)(P_act−1)+1,:); % set the corresponding initial condition

    tsub(iteT−1)=t(iteT)-(P_act−1)Deltattdir; % evaluate the solution at the shifted time

    temp=tsub(iteT−1);

    z(iteT,:)=nde(temp,c,tstart+t(iteT));

    if tdirt(iteT)>=tdirtfin0.99999 % repeat for the next sub-interval

        break

    else

        iteT=iteT+1; t(iteT)=tini+(iteT−1)incTtdir;

    end

end

t=t′+tstart;

end

%% EMHPM algorithm for DDE solutions

function t z=ddeEMHPM(dde,tspan,history,tau,Deltat,pntDelta)

% dde solver by Enhanced Multistage Homotopy Perturbation Method

tini=tspan(1); tfin=tspan(end); % span where the solution is founded

if tau<Deltat

    error(‘Subtinterval must not be greater than tau’);

end

if mod(tau,Deltat)~=0

    pastDlt=Deltat;

    Deltat=tau/round(tau/Deltat);

    warning(‘Subtinterval was modified from %0.5g to %0.5g.’,pastDlt,Deltat);

end

pnt=round(tau/Deltat); % samples−1 −tau 0

incT=Deltat/pntDelta; % step for set resolution

t=−tau+tini+(0:pntpntDelta)’incT; % evaluation of the initial solution −tau 0

z=history(t); % initial behavior in the inverval −tau 0

c=z(end,:); % initial condition

m=pntpntDelta; % samples between tau

itDlt=0;

while (t(end)+eps)<tfin;

    tspan=tini+itDltDeltat;tini+(itDlt+1)Deltat; % preparing the span for next Deltat

    zm=z(length(z)−m,:)’; zn=z(1+length(z)−m,:)’; % previous tau solution

    emhpm_fun=@(time,c,T) dde(time,c,T,zm,zn,tau,m); % application of the odeEMHPM

    t_aux z_aux=odeEMHPM(emhpm_fun,tspan,z(end,:),Deltat,pntDelta);

    z=z(1:end−1,:);z_aux; t=t(1:end−1,:);t_aux; % joining the solutions

    itDlt=itDlt+1;

end

z=z(pntpntDelta+1:end,:);t=t(pntpntDelta+1:end,:);

end

Algorithm 1: MATLAB algorithm.

Figure 6 shows the relative error between our approximate EMPHM and the dde23 solution and its relationship with the order and the discretized time intervals . Notice that the relative error values coincide at values of . Also, we can see from Figure 6 that the computed relative error values for approximate solutions of order remain unchanged.

Figure 6: Estimated relative error values between the numerical solution dde23 and the EMHPM approximate solutions. Here we use for the EMHPM the values of 2, 5, and 10.

5.2. A Practical Application: Cutting Operation on Milling Machine

We next use our EMHPM procedure to obtain the solution of the single degree-of-freedom milling operation. We use the simplified form based on [20–22]:

where is the angular natural frequency of the system, is the damping ratio, is the depth of cut, is the modal mass of the tool, represents the time delay which is equal to the tooth passing period, and is the specific cutting force coefficient which can be determined from

where is the tool number of teeth, and are the tangential and the normal linear cutting force coefficients, respectively, is the angular position of the -tooth defined as

and is the spindle speed in rpm [20]. The function is a switching function, which has a unity value when the -tooth is cutting and zero otherwise:

Here, and are the angles where the teeth enter and exit the workpiece. For upmilling, and , for downmilling, and , where is the radial depth of cut ratio.

By following the EMHPM procedure, we can write (46) in equivalent form as

where denotes the order solution of (46) on the th subinterval that satisfies the initial conditions , , and and is given by (35). Introducing the transformation , (50) can be written as a system of first order linear delay differential equations of the form

where

We next apply the EMHPM procedure to solve (46) by considering a downmilling operation with the following parameter values: , ,  rad⁄s, ,  kg,  N/m², and  N/m². As we can see from Figures 7 and 8 and for the depth of cut values of 2 mm (stable) and 3 mm (unstable), our EMHPM approximate solutions follow closely the numerical integration solutions of (46) obtained by using the dde23 algorithm.

Figure 7: EMHPM approximate solutions of (46) with parameter values of  mm,  rpm. Stable machine operation.

Figure 8: EMHPM approximate solutions of (46) with parameter values of  mm and  rpm. Unstable (chatter) machine operation.

Figure 9 shows the relative error between the EMPHM and the dde23 numerical solution, while Table 2 shows the CPU time needed for each solution. Here we use since the average step size of the dde23 algorithm is around . Note that, for , the zeroth order EMHPM approximate solution has the fastest CPU time. We can see from Figure 9 that the value of the relative error becomes basically the same for , 7, and 10 and .

Table 2: CPU time comparison among the approximate solutions of (46) by using dde23, the zeroth order, and the first order EMHPM solutions.

Solution method
Time [ms]

dde23
139

Zeroth order EMHPM ( = 7)
77

First order EMHPM ( = 7)
87

Figure 9: Estimated relative error values between dde23 and the EMHPM approximate solutions. Here we have used the system parameter values of  mm and  rpm and values of 2, 7, and 10.

6. Conclusions

We have developed a new algorithm based on the homotopy perturbation method to solve delay differential equations. The proposed EMHPM approach is based on a sequence of subintervals that approximate the solution of delayed differential equations by using the transformation rule , where satisfies the initial conditions. We have shown that our proposed EMHPM approach can be applied to obtain the approximate solution of a delay differential equation not only with constant, but also with variable coefficients with theoretical predictions that follow well the numerical integration solutions. To further assess the validity of this new approach, we have compared the approximate solutions of two delayed differential equations with respect to their corresponding numerical integration solutions obtained from the MATLAB dde23 algorithm. The test cases were (a) the damped Mathieu differential equation with time delay and (b) the governing equation of motion of downmilling operations. We have found that the EMHPM closely follows the numerical integration solutions of the corresponding equations and that these require less CPU time and have smaller relative errors.

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 Group in Nanomaterials for Medical Devices and the Research Group in Advanced Manufacturing. Additional support was provided from Consejo Nacional de Ciencia y Tecnología (Conacyt), Mexico.

References

F. Shakeri and M. Dehghan, “Solution of delay differential equations via a homotopy perturbation method,” Mathematical and Computer Modelling, vol. 48, no. 3-4, pp. 486–498, 2008.

Z.-K. Wu, “Solution of the enso delayed oscillator with homotopy analysis method,” Journal of Hydrodynamics, vol. 21, no. 1, pp. 131–135, 2009.

A. K. Alomari, M. S. Noorani, and R. Nazar, “Solution of delay differential equation by means of homotopy analysis method,” Acta Applicandae Mathematicae, vol. 108, no. 2, pp. 395–412, 2009.

J. Biazar and G. Behzad, “The homotopy perturbation method for solving neutral functional-differential equations with proportional delays,” Journal of King Saud University—Science, vol. 24, no. 1, pp. 33–37, 2012.

N. R. Anakira, A. K. Alomari, and I. Hashim, “Optimal homotopy asymptotic method for solving delay differential equations,” Mathematical Problems in Engineering, vol. 2013, Article ID 498902, 11 pages, 2013.

T. Insperger and G. Stépán, Semi-Discretization for Time-Delay Systems: Stability and Engineering Applications, Springer, New York, NY, USA, 2011.

E. A. Butcher, P. Nindujarla, and E. Bueler, “Stability of up- and down-milling using chebyshev collocation method,” in Proceedings of the ASME International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, vol. 6, part A–C, pp. 841–850, 2005.

T. Insperger, “Full-discretization and semi-discretization for milling stability prediction: some comments,” International Journal of Machine Tools and Manufacture, vol. 50, no. 7, pp. 658–662, 2010.

I. Hashim and M. S. Chowdhury, “Adaptation of homotopy-perturbation method for numeric-analytic solution of system of ODEs,” Physics Letters A, vol. 372, no. 4, pp. 470–481, 2008.

L. F. Shampine and S. Thompson, “Solving DDEs in MATLAB,” Applied Numerical Mathematics, vol. 37, no. 4, pp. 441–458, 2001.

S. H. H. Nia, A. N. Ranjbar, D. D. Ganji, H. Soltani, and J. Ghasemi, “Maintaining the stability of nonlinear differential equations by the enhancement of HPM,” Physics Letters A, vol. 372, no. 16, pp. 2855–2861, 2008.

J.-H. He, “Homotopy perturbation technique,” Computer Methods in Applied Mechanics and Engineering, vol. 178, no. 3-4, pp. 257–262, 1999.

J.-H. He, “Homotopy perturbation method: a new nonlinear analytical technique,” Applied Mathematics and Computation, vol. 135, no. 1, pp. 73–79, 2003.

J. H. He, “Comparison of homotopy perturbation method and homotopy analysis method,” Applied Mathematics and Computation, vol. 156, no. 2, pp. 527–539, 2004.

J. H. He, “New interpretation of homotopy perturbation method,” International Journal of Modern Physics B, vol. 20, no. 18, pp. 2561–2568, 2006.

S. Liao, “Comparison between the homotopy analysis method and homotopy perturbation method,” Applied Mathematics and Computation, vol. 169, no. 2, pp. 1186–1194, 2005.

D. Olvera and A. Elías-Zúñiga, “Enhanced multistage homotopy perturbation method: approximate solutions of nonlinear dynamic systems,” Abstract and Applied Analysis, vol. 2014, Article ID 486509, 12 pages, 2014.

T. Insperger and G. Stépán, “Stability of the damped Mathieu equation with time delay,” Journal of Dynamic Systems, Measurement, and Control, vol. 125, no. 2, pp. 166–171, 2003.

T. Insperger and G. Stépán, “Stability chart for the delayed Mathieu equation,” The Royal Society of London. Proceedings Series A: Mathematical, Physical and Engineering Sciences, vol. 458, no. 2024, pp. 1989–1998, 2002.

T. Insperger and G. Stépán, “Updated semi-discretization method for periodic delay-differential equations with discrete delay,” International Journal for Numerical Methods in Engineering, vol. 61, no. 1, pp. 117–141, 2004.

T. Insperger, B. P. Mann, G. Stépán, and P. V. Bayly, “Stability of up-milling and down-milling, part 1: alternative analytical methods,” International Journal of Machine Tools and Manufacture, vol. 43, no. 1, pp. 25–34, 2003.

F. I. Compeán, D. Olvera, F. J. Campa, L. N. L. de Lacalle, A. Elías-Zúñiga, and C. A. Rodríguez, “Characterization and stability analysis of a multivariable milling tool by the enhanced multistage homotopy perturbation method,” International Journal of Machine Tools and Manufacture, vol. 57, pp. 27–33, 2012.