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 [13] 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 [1315].

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 .

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.

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).

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 .

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.

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/(N1); ktau=2; tspan=0,ktauT; % EMHPM Parameters
%% Solution
% Solution by dde23
tdde=linspace(0,ktautau,ktau(N1)+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(N1); ind1=2:2:ktau*(N1);
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(N1)/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=1E3+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=tinitstart; tfin=tfintstart; % 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(tfintini);
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_act1)+1,:); % set the corresponding initial condition
    tsub(iteT1)=t(iteT)-(P_act1)Deltattdir; % evaluate the solution at the shifted time
    temp=tsub(iteT1);
    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+(iteT1)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); % samples1 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:end1,:);z_aux; t=t(1:end1,:);t_aux; % joining the solutions
    itDlt=itDlt+1;
end
z=z(pntpntDelta+1:end,:);t=t(pntpntDelta+1:end,:);
end

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.

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 [2022]: 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: , ,  rads, ,  kg,  N/m2, and  N/m2. 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 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 .

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.