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

Authors: D. Olvera; A. Elías-Zúñiga; L. N. López de Lacalle; C. A. Rodríguez
Journal: Abstract and Applied Analysis
 (2015)

---

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


---

## Body

## 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 B
=
0. 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 Section2 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:(1)
A
u
-
f
r
=
0
,
r
∈
Ω
with boundary conditions
(2)
B
u
,
∂
u
∂
n
=
0
,
r
∈
Γ
,
where A is a general differential operator, B is a boundary operator, f
r is a known analytic function, and Γ is the boundary of the domain Ω.The operatorA can generally be divided into two parts: L and N, where L involves the linear terms and N the nonlinear ones. Equation (1) therefore can be rewritten as follows:
(3)
L
u
+
N
u
-
f
r
=
0
.
By the homotopy perturbation technique, we construct a homotopy v
r
,
p
:
Ω
×
0,1
→
R that satisfies
(4)
H
v
,
p
=
L
v
-
L
u
0
+
p
L
u
0
+
p
N
v
-
f
r
=
0
,
where p
∈
0,1 is an embedding parameter and u
0 is an initial approximation of (1) which satisfies the boundary conditions (2). Thus, from (4), we have
(5)
H
v
,
0
=
L
v
-
L
u
0
=
0
,
H
v
,
1
=
A
v
-
f
r
=
0
.
The changing process of p from zero to unity is just that of v
r
,
p from u
0
r to u
r. In topology, this is called deformation, and L
v
-
L
u
0 and A
v
-
f
r are called homotopic.He in [12] uses the embedding parameter p as the small parameter and assumed that the solution of (4) can be written as a power series of p in the form
(6)
v
=
v
0
+
p
v
1
+
p
2
v
2
+
⋯
.
By setting p
=
1, He obtained the approximate solution of (1) as
(7)
u
=
lim
⁡
p
→
1
=
v
0
+
v
1
+
v
2
+
⋯
.
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Δ
t
→
0. 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:u
t
→
u
i
T, where u
i satisfies the initial condition u
i
(
0
)
=
u
i
-
1
(
t
i
-
1
), T is a shifted time scale used to determine the approximate solution in each subinterval, and u
i
T represents the approximate solution in the ith subinterval. In this case, the initial suggested solution in the ith subinterval is given by u
i
0
(
T
)
=
u
i
-
1
(
t
i
-
1
), where t
i
-
1 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 operatorL
u
i represents L
u
i
=
d
/
d
T
u
i, where the assumed approximate solution u
i
0
T is set equal to the initial condition u
i
-
1
t
i
-
1; that is, u
i
0
=
u
i
-
1
(
t
i
-
1
). To simplify the notation, we let u
i
-
1
≡
u
i
-
1
t
i
-
1.
(2)
The transformationT
=
t
-
t
i
-
1 on 0
<
T
≤
t
i
-
t
i
-
1 holds in the homotopy i-subinterval. Thus, higher order equations are integrated with respect to T, while the terms related to the independent variable t are assumed to remain constant.Therefore, we may conclude that the m order approximate solution, by applying the EMHPM, can be written as
(8)
u
i
T
,
u
i
-
1
=
∑
k
=
0
m
U
i
k
T
,
u
i
-
1
,
where the solution u
i
T
,
u
i
-
1 is valid only in the ith subinterval t
i
-
1
,
t
i. Hence, the solution u
t on the ith subinterval t
i
-
1
,
t
i can be written as
(9)
u
t
≈
u
i
t
-
t
i
-
1
with initial condition u
i
-
1
(
t
i
-
1
), and i
=
1,2
,
…
,
j. Thus, the approximate solution of u at the time t
i is given by
(10)
u
i
t
-
t
i
-
1
t
=
t
i
=
u
i
+
1
t
-
t
i
t
=
t
i
=
u
i
+
1
0
=
u
i
.
In summary, the solution u
t for an open-closed interval (t
0
,
t
1] is divided into j subintervals that, in general, are not equally spaced: [
t
0
,
t
1
]
,
[
t
1
,
t
2
]
,
…
,
[
t
j
-
1
,
t
j
]. Thus, the approximated solution of u
t for the span time interval is obtained by coupling the u
i
t 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(11)
x
˙
t
+
x
t
-
τ
=
0
with initial condition x
(
0
)
=
c. Here, the independent variable x is a scalar x
t
∈
R, the dot stands for differentiation with respect to time t, and τ is the time delay. To evaluate (11) on a
≤
t
≤
b, the term x
t
-
τ must represent a known function x
t on a
-
τ
≤
t
≤
a. For instance, if a
=
0, the solution of (11) can be obtained in the interval (
0
,
τ
] 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 ith interval i
-
1
τ
,
i
τ, i
=
2,3
,
…
,
j, where j is an integer number that can be chosen as 2
≤
j
≤
∞. With this approach, we can apply the HPM to find the solution of (11) by assuming that the previous delayed function is x
τ
0
T
=
c; thus the solution for the first interval is given by x
τ
1
T, valid on 0
,
τ. In terms of (4), we now construct the homotopy of (11):
(12)
H
X
τ
1
,
p
=
d
d
T
X
τ
1
+
p
x
τ
0
=
0
.
We next substitute the first order expansion X
τ
1
=
X
0
τ
1
+
p
X
1
τ
1 in (12) and balance the terms with identical power of p to obtain the following set of linear differential equations:
(13)
p
0
:
d
d
T
X
0
τ
1
=
0
X
0
τ
1
0
=
c
=
X
τ
0
τ
,
p
1
:
d
d
T
X
1
τ
1
=
-
X
τ
0
X
1
τ
1
0
=
0
.
Integration of (13) yields
(14)
X
0
τ
1
=
c
,
X
1
τ
1
=
-
c
T
.
Hence, the first order solution of (12) is given by
(15)
x
τ
1
T
=
c
-
c
T
.
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
(16)
x
τ
2
T
=
c
-
c
τ
-
c
T
+
1
2
c
T
2
,
x
τ
3
T
=
c
-
2
c
τ
+
1
2
c
τ
2
-
c
-
c
τ
T
+
1
2
c
T
2
-
1
6
c
T
3
.
Figure 1 shows the exact solution of (11) obtained by coupling at each interval the solution obtained by following HPM procedure for t
=
10
τ.Figure 1
Exact solution of (11) obtained by using the HPM and τ
=
1.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(17)
x
˙
+
x
t
-
τ
-
cos
⁡
π
t
x
2
=
0
,
τ
=
1
,
x
0
=
c
=
x
τ
0
τ
in which the solution x
τ
0
T
=
c
1 holds on (
-
τ
,
0
]. In order to find the solution x
τ
1 in the interval 0
,
τ, we assume that the homotopy representation of (17) can be given as
(18)
H
X
τ
1
,
p
=
d
d
T
X
τ
1
+
p
x
τ
0
-
cos
⁡
π
t
X
τ
1
2
=
0
.
Notice that the variable X depends on the time T for which 0
≤
T
≤
τ. If we now substitute the second order expansion X
τ
1
=
X
0
τ
1
+
p
X
1
τ
1
+
p
2
X
2
τ
1 in (18), and, after balancing the p terms, we get that
(19)
p
0
:
d
d
T
X
0
τ
1
=
0
,
X
0
T
=
0
=
c
1
=
X
τ
0
T
=
τ
,
p
1
:
d
d
T
X
1
τ
1
=
-
x
τ
0
+
cos
⁡
π
t
X
0
τ
1
2
=
0
,
X
1
0
=
0
,
p
2
:
d
d
T
X
2
τ
1
=
2
cos
⁡
π
t
X
0
τ
1
X
1
τ
1
,
X
2
0
=
0
,
p
3
:
d
d
T
X
3
τ
1
=
cos
⁡
π
t
2
X
0
τ
1
X
2
τ
1
+
X
1
τ
1
2
=
0
,
X
3
0
=
0
.
Equations (19) have the following solutions:
(20)
X
0
τ
1
=
c
1
,
X
1
τ
1
=
-
T
x
τ
0
-
c
1
2
cos
⁡
π
t
,
X
2
τ
1
=
-
c
1
T
2
cos
⁡
π
t
x
τ
0
-
c
1
2
cos
⁡
π
t
,
X
3
τ
1
=
1
3
T
3
cos
⁡
π
t
3
c
1
4
cos
⁡
2
π
t
-
4
c
1
2
x
τ
0
+
x
τ
0
2
.
Thus, the approximate solution of (17) by using the EMHPM is given by
(21)
x
τ
1
T
≈
X
0
τ
1
+
X
1
τ
1
+
X
2
τ
1
+
X
3
τ
1
.
In this case, the exact solution of x
τ
1
(
T
) is unknown. To obtain x
τ
2, we compute again the approximate solution of x
τ
1
(
T
) by applying our EMHPM and the value of the delayed time is assumed to remain constant in each subinterval. To determine x
τ
2, we next use the homotopy representation of (17) for the interval (
τ
,
2
τ
]:
(22)
H
X
τ
2
,
p
=
d
d
T
X
τ
2
+
p
x
τ
1
-
cos
⁡
π
t
X
τ
2
2
=
0
.
Substituting the second order expansion in (22), we get
(23)
X
0
τ
2
=
c
2
X
1
τ
2
=
-
T
x
τ
1
-
c
2
2
cos
⁡
π
t
X
2
τ
2
=
-
c
2
T
2
cos
⁡
π
t
x
τ
1
-
c
2
2
cos
⁡
π
t
X
3
τ
2
=
1
3
T
3
cos
⁡
π
t
3
c
2
4
cos
⁡
2
π
t
-
4
c
2
2
x
τ
1
+
x
τ
1
2
.
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:
(24)
H
X
τ
i
,
p
=
d
d
T
X
τ
i
+
p
X
τ
i
-
1
-
cos
⁡
π
t
X
τ
i
2
=
0
.
Then, by using our EMPHM approach, we have that
(25)
X
0
τ
i
=
c
X
1
τ
i
=
-
T
x
τ
i
-
1
-
c
2
cos
⁡
π
t
X
2
τ
i
=
-
c
T
2
cos
⁡
π
t
x
τ
i
-
1
-
c
2
cos
⁡
π
t
X
3
τ
i
=
1
3
T
3
cos
⁡
π
t
3
c
4
cos
⁡
2
π
t
-
4
c
2
x
τ
i
-
1
+
x
τ
i
-
1
2
.
Notice from (25) that the kth order approximate solution of (17) can be written as
(26)
X
0
τ
i
=
c
X
k
τ
i
=
T
k
-
x
τ
i
-
1
g
k
+
cos
⁡
π
t
∑
n
1
=
0
k
-
1
X
n
1
τ
i
X
k
-
1
-
n
1
τ
i
,
where k
>
0, g
k
=
1 when k
=
1 and zero otherwise.Figure2 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 x
τ
0
T
=
cos
⁡
π
T
+
1, x
τ
0
T
=
e
T
+
1, and a time subintervals Δ
t
=
0.01. 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.


## 4.1. Delay Differential Equations with Constant Coefficients

First, let us consider the simplest DDE of the form(11)
x
˙
t
+
x
t
-
τ
=
0
with initial condition x
(
0
)
=
c. Here, the independent variable x is a scalar x
t
∈
R, the dot stands for differentiation with respect to time t, and τ is the time delay. To evaluate (11) on a
≤
t
≤
b, the term x
t
-
τ must represent a known function x
t on a
-
τ
≤
t
≤
a. For instance, if a
=
0, the solution of (11) can be obtained in the interval (
0
,
τ
] 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 ith interval i
-
1
τ
,
i
τ, i
=
2,3
,
…
,
j, where j is an integer number that can be chosen as 2
≤
j
≤
∞. With this approach, we can apply the HPM to find the solution of (11) by assuming that the previous delayed function is x
τ
0
T
=
c; thus the solution for the first interval is given by x
τ
1
T, valid on 0
,
τ. In terms of (4), we now construct the homotopy of (11):
(12)
H
X
τ
1
,
p
=
d
d
T
X
τ
1
+
p
x
τ
0
=
0
.
We next substitute the first order expansion X
τ
1
=
X
0
τ
1
+
p
X
1
τ
1 in (12) and balance the terms with identical power of p to obtain the following set of linear differential equations:
(13)
p
0
:
d
d
T
X
0
τ
1
=
0
X
0
τ
1
0
=
c
=
X
τ
0
τ
,
p
1
:
d
d
T
X
1
τ
1
=
-
X
τ
0
X
1
τ
1
0
=
0
.
Integration of (13) yields
(14)
X
0
τ
1
=
c
,
X
1
τ
1
=
-
c
T
.
Hence, the first order solution of (12) is given by
(15)
x
τ
1
T
=
c
-
c
T
.
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
(16)
x
τ
2
T
=
c
-
c
τ
-
c
T
+
1
2
c
T
2
,
x
τ
3
T
=
c
-
2
c
τ
+
1
2
c
τ
2
-
c
-
c
τ
T
+
1
2
c
T
2
-
1
6
c
T
3
.
Figure 1 shows the exact solution of (11) obtained by coupling at each interval the solution obtained by following HPM procedure for t
=
10
τ.Figure 1
Exact solution of (11) obtained by using the HPM and τ
=
1.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(17)
x
˙
+
x
t
-
τ
-
cos
⁡
π
t
x
2
=
0
,
τ
=
1
,
x
0
=
c
=
x
τ
0
τ
in which the solution x
τ
0
T
=
c
1 holds on (
-
τ
,
0
]. In order to find the solution x
τ
1 in the interval 0
,
τ, we assume that the homotopy representation of (17) can be given as
(18)
H
X
τ
1
,
p
=
d
d
T
X
τ
1
+
p
x
τ
0
-
cos
⁡
π
t
X
τ
1
2
=
0
.
Notice that the variable X depends on the time T for which 0
≤
T
≤
τ. If we now substitute the second order expansion X
τ
1
=
X
0
τ
1
+
p
X
1
τ
1
+
p
2
X
2
τ
1 in (18), and, after balancing the p terms, we get that
(19)
p
0
:
d
d
T
X
0
τ
1
=
0
,
X
0
T
=
0
=
c
1
=
X
τ
0
T
=
τ
,
p
1
:
d
d
T
X
1
τ
1
=
-
x
τ
0
+
cos
⁡
π
t
X
0
τ
1
2
=
0
,
X
1
0
=
0
,
p
2
:
d
d
T
X
2
τ
1
=
2
cos
⁡
π
t
X
0
τ
1
X
1
τ
1
,
X
2
0
=
0
,
p
3
:
d
d
T
X
3
τ
1
=
cos
⁡
π
t
2
X
0
τ
1
X
2
τ
1
+
X
1
τ
1
2
=
0
,
X
3
0
=
0
.
Equations (19) have the following solutions:
(20)
X
0
τ
1
=
c
1
,
X
1
τ
1
=
-
T
x
τ
0
-
c
1
2
cos
⁡
π
t
,
X
2
τ
1
=
-
c
1
T
2
cos
⁡
π
t
x
τ
0
-
c
1
2
cos
⁡
π
t
,
X
3
τ
1
=
1
3
T
3
cos
⁡
π
t
3
c
1
4
cos
⁡
2
π
t
-
4
c
1
2
x
τ
0
+
x
τ
0
2
.
Thus, the approximate solution of (17) by using the EMHPM is given by
(21)
x
τ
1
T
≈
X
0
τ
1
+
X
1
τ
1
+
X
2
τ
1
+
X
3
τ
1
.
In this case, the exact solution of x
τ
1
(
T
) is unknown. To obtain x
τ
2, we compute again the approximate solution of x
τ
1
(
T
) by applying our EMHPM and the value of the delayed time is assumed to remain constant in each subinterval. To determine x
τ
2, we next use the homotopy representation of (17) for the interval (
τ
,
2
τ
]:
(22)
H
X
τ
2
,
p
=
d
d
T
X
τ
2
+
p
x
τ
1
-
cos
⁡
π
t
X
τ
2
2
=
0
.
Substituting the second order expansion in (22), we get
(23)
X
0
τ
2
=
c
2
X
1
τ
2
=
-
T
x
τ
1
-
c
2
2
cos
⁡
π
t
X
2
τ
2
=
-
c
2
T
2
cos
⁡
π
t
x
τ
1
-
c
2
2
cos
⁡
π
t
X
3
τ
2
=
1
3
T
3
cos
⁡
π
t
3
c
2
4
cos
⁡
2
π
t
-
4
c
2
2
x
τ
1
+
x
τ
1
2
.
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:
(24)
H
X
τ
i
,
p
=
d
d
T
X
τ
i
+
p
X
τ
i
-
1
-
cos
⁡
π
t
X
τ
i
2
=
0
.
Then, by using our EMPHM approach, we have that
(25)
X
0
τ
i
=
c
X
1
τ
i
=
-
T
x
τ
i
-
1
-
c
2
cos
⁡
π
t
X
2
τ
i
=
-
c
T
2
cos
⁡
π
t
x
τ
i
-
1
-
c
2
cos
⁡
π
t
X
3
τ
i
=
1
3
T
3
cos
⁡
π
t
3
c
4
cos
⁡
2
π
t
-
4
c
2
x
τ
i
-
1
+
x
τ
i
-
1
2
.
Notice from (25) that the kth order approximate solution of (17) can be written as
(26)
X
0
τ
i
=
c
X
k
τ
i
=
T
k
-
x
τ
i
-
1
g
k
+
cos
⁡
π
t
∑
n
1
=
0
k
-
1
X
n
1
τ
i
X
k
-
1
-
n
1
τ
i
,
where k
>
0, g
k
=
1 when k
=
1 and zero otherwise.Figure2 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 x
τ
0
T
=
cos
⁡
π
T
+
1, x
τ
0
T
=
e
T
+
1, and a time subintervals Δ
t
=
0.01. 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 ann-dimensional delay differential equation of the form
(27)
x
˙
t
=
A
t
x
+
B
t
x
t
-
τ
,
where A
t
+
τ
=
A
t, B
t
+
τ
=
B
t, x
t is the state vector, and τ is the time delay. By following our EMHPM procedure, we can write (27) in equivalent form as
(28)
x
˙
i
T
-
A
t
x
i
T
≈
B
t
x
i
τ
T
,
where x
i
T denotes the m order solution of (27) in the ith subinterval that satisfies the initial conditions x
i
0
=
x
i
-
1 and A
t and B
t represent the values of the periodic coefficients at the time t. In order to approximate the delayed term x
i
τ
T in (28), the period [
t
0
-
τ
,
t
0
] is discretized in N points equally spaced as shown in Figure 3. Here, we assume that the function x
i
τ
T in the delay subinterval [
t
i
-
N
,
t
i
-
N
+
1
] is approximated by a constant value
(29)
x
i
τ
T
=
x
i
-
N
+
1
T
≈
x
i
-
N
as shown in Figure 3. By following the homotopy perturbation technique, we can write the homotopy representation of (28) as
(30)
H
X
i
,
p
=
L
X
i
-
L
x
i
0
+
p
L
x
i
0
=
p
A
t
X
i
+
B
t
x
i
-
N
.
Substituting the m order expansion X
i
=
X
i
0
+
p
X
i
1
+
⋯
+
p
m
X
i
m in (30) and by assuming an initial approximation of the form x
i
0
=
x
i
-
1, we get, after applying the proposed EMHPM approach, the following set of first order linear delay differential equations:
(31)
p
0
:
d
d
T
X
i
0
+
d
d
T
x
i
-
1
=
0
,
X
i
0
=
x
i
-
1
,
p
1
:
d
d
T
X
i
1
=
A
t
X
i
0
+
B
t
x
i
-
N
,
X
i
1
0
=
0
,
p
2
:
d
d
T
X
i
2
=
A
t
X
i
1
,
X
i
2
0
=
0
,
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
⋮
p
m
:
d
d
T
X
i
m
=
A
t
X
i
m
-
1
,
X
i
m
0
=
0
.
By solving (31), we get
(32)
X
i
0
=
x
i
-
1
,
X
i
1
=
A
t
x
i
-
1
T
+
B
t
x
i
-
N
T
,
X
i
2
=
1
2
A
t
2
x
i
-
1
T
2
+
1
2
A
t
B
t
x
i
-
N
T
2
,
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
⋮
X
i
m
=
1
m
!
A
t
m
x
i
-
1
T
m
+
1
m
!
A
t
m
-
1
B
t
x
i
-
N
T
m
.
Equations (32) can be written as
(33)
X
i
k
=
T
k
A
t
X
i
k
-
1
+
g
k
B
t
x
i
-
N
,
k
=
1,2
,
3
,
…
,
where X
i
0
=
x
i
-
1 and g
k
=
1 for k
=
1 and g
k
=
0, otherwise. Thus, the solution of (27) is obtained by adding the X
i
k approximate solutions:
(34)
x
i
T
≈
∑
k
=
0
m
X
i
k
T
.
Notice, however, that solution (34) may be further improved by using a first order polynomial representation of x
i
τ
T as shown in Figure 4. Then, the function x
i
τ
T in the delay subinterval [
t
i
-
N
,
t
i
-
N
+
1
] takes the form
(35)
x
i
τ
T
=
x
i
-
N
+
1
T
≈
x
i
-
N
+
N
-
1
τ
x
i
-
N
+
1
-
x
i
-
N
T
.
Substituting (35) into (28) gives
(36)
x
˙
i
T
-
A
t
x
i
T
≈
B
t
x
i
-
N
-
N
-
1
τ
B
t
x
i
-
N
T
+
N
-
1
τ
B
t
x
i
-
N
+
1
T
.
We next assume that the homotopy representation of (36) is given as
(37)
H
X
i
,
p
=
L
X
i
-
L
x
i
0
+
p
L
x
i
0
-
p
A
X
i
+
B
x
i
-
N
-
(
N
-
1
)
τ
B
x
i
-
N
T
lllllll
;
lllll
+
N
-
1
τ
B
x
i
-
N
+
1
T
=
0
.
Substituting the m order expansion  X
i
T
=
X
i
0
T
+
p
X
i
1
T
+
⋯
p
m
X
i
m
T  in (37) and assuming that the initial approximation is given by x
i
0
=
x
i
-
1, we get
(38)
p
0
:
d
d
T
X
i
0
+
d
d
T
x
i
-
1
=
0
,
X
i
0
=
x
i
-
1
,
p
1
:
d
d
T
X
i
1
=
A
t
X
i
0
+
B
t
x
i
-
N
-
N
-
1
τ
B
t
x
i
-
N
T
llllllllllllllllllll
+
N
-
1
τ
B
t
x
i
-
N
+
1
T
,
X
i
1
0
=
0
,
p
2
:
d
d
T
X
i
2
=
A
X
i
1
,
X
i
2
0
=
0
,
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
⋮
p
m
:
d
d
T
X
i
m
=
A
X
i
m
-
1
,
X
i
m
0
=
0
.
By solving (38) and by following the EMHPM procedure, we get
(39)
X
i
0
=
x
i
-
1
,
X
i
1
=
A
t
x
i
-
1
T
+
B
t
x
i
-
N
T
-
1
2
N
-
1
τ
B
t
x
i
-
N
T
2
lllllllll
+
1
2
N
-
1
τ
B
t
x
i
-
N
+
1
T
2
,
X
i
2
=
1
2
A
t
2
x
i
-
1
T
2
+
1
2
A
t
B
t
x
i
-
N
T
2
lllllllll
-
1
6
N
-
1
τ
A
t
B
t
x
i
-
N
T
3
+
1
6
N
-
1
τ
A
t
B
t
x
i
-
N
+
1
T
3
,
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
⋮
X
i
m
=
1
m
!
A
t
m
x
i
-
1
T
m
+
1
m
!
A
t
m
-
1
B
t
x
i
-
N
T
m
llllllllll
-
1
m
+
1
!
N
-
1
τ
A
t
m
-
1
B
t
x
i
-
N
T
m
+
1
llllllllll
+
1
m
+
1
!
N
-
1
τ
A
t
m
-
1
B
t
x
i
-
N
+
1
T
m
+
1
.
Here, the recursive form of X
i
k
T is written as
(40)
X
i
k
=
X
i
k
a
+
X
i
k
b
k
=
1,2
,
3
,
…
,
where X
i
0
a
=
x
i
-
1
,
X
i
0
b
=
0 and
(41)
X
i
k
a
=
T
k
A
t
X
i
k
-
1
a
+
g
k
B
t
x
i
-
N
,
X
i
k
b
=
T
k
+
1
N
-
1
τ
T
-
B
t
x
i
-
N
+
B
t
x
i
-
N
+
1
A
t
X
i
k
-
1
b
llllllllllllllllllll
+
g
k
N
-
1
τ
T
-
B
t
x
i
-
N
+
B
t
x
i
-
N
+
1
.
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:(42)
x
¨
+
κ
x
˙
+
δ
+
ε
cos
⁡
2
π
t
T
x
=
b
x
t
-
τ
,
where κ, δ, ε, τ, and T 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:
(43)
x
¨
i
T
+
κ
x
˙
i
T
+
α
t
x
i
T
≈
b
x
i
-
N
+
1
T
,
where x
i
t denotes the m order solution of (43) in the ith subinterval that satisfies the following initial conditions: x
i
0
=
x
i
-
1 and x
˙
i
0
=
x
˙
i
-
1. The space state form representation of (43) is given by
(44)
x
˙
i
T
=
A
t
x
i
T
+
B
t
x
i
-
N
+
1
T
,
where
(45)
A
t
=
0
1
-
α
t
-
κ
,
B
t
=
0
0
b
t
0
and α
t
=
δ
+
ε
cos
⁡
t 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 κ
=
0.2, δ
=
3.0, ε
=
1, b
=
-
1, and T
=
τ
=
2
π. 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 2
T, by assuming that N
=
50 with the following initial values: x
-
50
T
=
x
-
49
T
=
⋯
x
0
T
=
0.001 and x
˙
-
50
T
=
x
˙
-
49
T
=
⋯
x
˙
0
T
=
0.Figure 5
Numerical solutions of the damped Mathieu equation with time delay by dde23, the zeroth EMHPM, and the first EMHPM withN
=
50 and m
=
4.It can be seen from Figure5 that in the interval [0, T] 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 [T, 2
T] 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 m and the discretized time intervals N 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. Them order solution of the EMPHM approach is chosen to guarantee the convergence of its approximate solution.
dde23[ms]
EMHPM
N
m
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
21Algorithm 1:MATLAB algorithm.
f
unction dde23_ddeEM_mathieu_paper
disp(‘Mathieu equation solution’)
% Solution by EMHPM with zeroth and first order solution
eps=1; kappa=0.2; T=2
*
pi; tau=T;  b=
−
1; delta=3; % Mathieu Parameters
pntDelay=1; N=50; m_ord=5; dt=tau/(N
−
1); ktau=2; tspan=
[
0,ktau
*
T
]
; % EMHPM Parameters
%% Solution
% Solution by dde23
tdde=linspace(0,ktau
*
tau,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)
-kapa
*
y(2)−
(dlt+eps
*
cos(2
*
pi
*
t/T))
*
y(1)+b
*
z(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+eps
*
cos(2
*
pi/T
*
tt);
for ik=1:m
Z(ik+1,1)=Z(ik,2)
*
t/ik;
Z(ik+1,2)=
−
kpa
*
Z(ik,2)
−
alf
*
Z(ik,1);
if ik==1, Z(ik+1,2)=Z(ik+1,2)+b
*
zm(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+eps
*
cos(2
*
pi/T
*
tt); Z=
[
c0(1),c0(2)]
;  Z_=
[
0,0
]
;
for ik=1:m
Z(ik+1,:)=
[
Z(ik,2)
*
t/ik,
−
kpa
*
Z(ik,2)
−
alf
*
Z(ik,1)
]
;
Z_(ik+1,:)=
[
Z_(ik,2)
*
t/ik,
−
kpa
*
Z_(ik,2)
−
alf
*
Z_(ik,1)
]
;
if ik==1,
Z(ik+1,2)=Z(ik+1,2)+b
*
zm(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+0
*
t 0+0
*
t
]
;
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+incT
*
tdir;
while tdir
*
t(iteT)<tdir*(tfin+tdir
*
incT)
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)
*
Deltat
*
tdir; % evaluate the solution at the shifted time
temp=tsub(iteT
−
1);
z(iteT,:)=nde(temp,c,tstart+t(iteT));
if tdir
*
t(iteT)>=tdir
*
tfin
*
0.99999 % repeat for the next sub-interval
break
else
iteT=iteT+1;  t(iteT)=tini+(iteT
−
1)
*
incT
*
tdir;
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:pnt
*
pntDelta)’*
incT; % evaluation of the initial solution
[
−
tau 0
]
z=history(t); % initial behavior in the inverval
[
−
tau 0
]
c=z(end,:); % initial condition
m=pnt
*
pntDelta; % samples between tau
itDlt=0;
while (t(end)+eps)<tfin;
tspan=
[
tini+itDlt
*
Deltat;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(pnt
*
pntDelta+1:end,:);t=t(pnt
*
pntDelta+1:end,:);
endFigure6 shows the relative error between our approximate EMPHM and the dde23 solution and its relationship with the order m and the discretized time intervals N. Notice that the relative error values coincide at values of N
≥
45. Also, we can see from Figure 6 that the computed relative error values for approximate solutions of order m
≥
5 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 ofm
= 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]:
(46)
x
¨
t
+
2
ζ
ω
n
x
˙
t
+
ω
n
2
x
t
=
-
a
p
K
s
t
m
m
x
t
-
x
t
-
τ
,
where ω
n is the angular natural frequency of the system, ζ is the damping ratio, a
p is the depth of cut, m
m is the modal mass of the tool, τ represents the time delay which is equal to the tooth passing period, and K
s
t is the specific cutting force coefficient which can be determined from
(47)
K
s
t
=
∑
j
=
1
z
n
g
ϕ
j
t
sin
ϕ
j
t
K
t
cos
⁡
ϕ
j
t
+
K
n
sin
ϕ
j
t
,
where z
n is the tool number of teeth, K
t and K
n are the tangential and the normal linear cutting force coefficients, respectively, ϕ
j
t is the angular position of the j-tooth defined as
(48)
ϕ
j
t
=
2
π
n
60
t
+
2
π
j
z
n
,
and n is the spindle speed in rpm [20]. The function g
φ
j
t is a switching function, which has a unity value when the j-tooth is cutting and zero otherwise:
(49)
g
ϕ
j
t
=
1
ϕ
st
<
ϕ
j
t
<
ϕ
ex
0
otherwise
.
Here, ϕ
st and ϕ
ex are the angles where the teeth enter and exit the workpiece. For upmilling, ϕ
st
=
0 and ϕ
ex
=
arccos
1
-
2
a
d, for downmilling, ϕ
st
=
arccos
2
a
d
-
1 and ϕ
ex
=
π, where a
d is the radial depth of cut ratio.By following the EMHPM procedure, we can write (46) in equivalent form as
(50)
x
¨
i
T
+
2
ζ
ω
n
x
˙
i
T
+
ω
n
2
x
i
T
≈
-
a
p
K
st
m
m
x
i
T
-
x
i
-
N
+
1
T
,
where x
i
T denotes the m order solution of (46) on the ith subinterval that satisfies the initial conditions x
i
0
=
x
i
-
1, x
˙
i
0
=
x
˙
i
-
1, and h
t
=
h
t and x
-
τ is given by (35). Introducing the transformation x
i
=
x
i
,
x
˙
i
T, (50) can be written as a system of first order linear delay differential equations of the form
(51)
x
˙
i
T
=
A
t
x
i
T
+
B
t
x
i
-
N
+
1
T
,
where
(52)
A
t
=
0
1
-
ω
n
2
-
w
m
m
K
st
-
2
ζ
ω
n
;
B
t
=
0
0
w
m
m
K
st
0
.
We next apply the EMHPM procedure to solve (46) by considering a downmilling operation with the following parameter values: z
n
=
2, a
d
=
0.1, ω
n
=
5793 rad⁄s, ζ
=
0.011, m
m
=
0.03993 kg, K
t
=
6
×
10
8 N/m2, and K
n
=
2
×
10
8 N/m2. As we can see from Figures 7 and 8 and for the depth of cut values of a
p
= 2 mm (stable) and a
p
= 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 a
p
=
2 mm, n
=
10000 rpm. Stable machine operation.Figure 8
EMHPM approximate solutions of (46) with parameter values of a
p
=
3 mm and n
=
10000 rpm. Unstable (chatter) machine operation.Figure9 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 N
=
75 since the average step size of the dde23 algorithm is around Δ
t
≈
τ
/
N. Note that, for m
=
7, 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 m
=
2, 7, and 10 and N
≥
20.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 (m = 7)
77
First order EMHPM (m = 7)
87Figure 9
Estimated relative error values between dde23 and the EMHPM approximate solutions. Here we have used the system parameter values ofa
p
=
2 mm and n
=
1000 rpm and values of m
= 2, 7, and 10.


## 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:(42)
x
¨
+
κ
x
˙
+
δ
+
ε
cos
⁡
2
π
t
T
x
=
b
x
t
-
τ
,
where κ, δ, ε, τ, and T 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:
(43)
x
¨
i
T
+
κ
x
˙
i
T
+
α
t
x
i
T
≈
b
x
i
-
N
+
1
T
,
where x
i
t denotes the m order solution of (43) in the ith subinterval that satisfies the following initial conditions: x
i
0
=
x
i
-
1 and x
˙
i
0
=
x
˙
i
-
1. The space state form representation of (43) is given by
(44)
x
˙
i
T
=
A
t
x
i
T
+
B
t
x
i
-
N
+
1
T
,
where
(45)
A
t
=
0
1
-
α
t
-
κ
,
B
t
=
0
0
b
t
0
and α
t
=
δ
+
ε
cos
⁡
t 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 κ
=
0.2, δ
=
3.0, ε
=
1, b
=
-
1, and T
=
τ
=
2
π. 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 2
T, by assuming that N
=
50 with the following initial values: x
-
50
T
=
x
-
49
T
=
⋯
x
0
T
=
0.001 and x
˙
-
50
T
=
x
˙
-
49
T
=
⋯
x
˙
0
T
=
0.Figure 5
Numerical solutions of the damped Mathieu equation with time delay by dde23, the zeroth EMHPM, and the first EMHPM withN
=
50 and m
=
4.It can be seen from Figure5 that in the interval [0, T] 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 [T, 2
T] 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 m and the discretized time intervals N 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. Them order solution of the EMPHM approach is chosen to guarantee the convergence of its approximate solution.
dde23[ms]
EMHPM
N
m
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
21Algorithm 1:MATLAB algorithm.
f
unction dde23_ddeEM_mathieu_paper
disp(‘Mathieu equation solution’)
% Solution by EMHPM with zeroth and first order solution
eps=1; kappa=0.2; T=2
*
pi; tau=T;  b=
−
1; delta=3; % Mathieu Parameters
pntDelay=1; N=50; m_ord=5; dt=tau/(N
−
1); ktau=2; tspan=
[
0,ktau
*
T
]
; % EMHPM Parameters
%% Solution
% Solution by dde23
tdde=linspace(0,ktau
*
tau,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)
-kapa
*
y(2)−
(dlt+eps
*
cos(2
*
pi
*
t/T))
*
y(1)+b
*
z(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+eps
*
cos(2
*
pi/T
*
tt);
for ik=1:m
Z(ik+1,1)=Z(ik,2)
*
t/ik;
Z(ik+1,2)=
−
kpa
*
Z(ik,2)
−
alf
*
Z(ik,1);
if ik==1, Z(ik+1,2)=Z(ik+1,2)+b
*
zm(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+eps
*
cos(2
*
pi/T
*
tt); Z=
[
c0(1),c0(2)]
;  Z_=
[
0,0
]
;
for ik=1:m
Z(ik+1,:)=
[
Z(ik,2)
*
t/ik,
−
kpa
*
Z(ik,2)
−
alf
*
Z(ik,1)
]
;
Z_(ik+1,:)=
[
Z_(ik,2)
*
t/ik,
−
kpa
*
Z_(ik,2)
−
alf
*
Z_(ik,1)
]
;
if ik==1,
Z(ik+1,2)=Z(ik+1,2)+b
*
zm(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+0
*
t 0+0
*
t
]
;
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+incT
*
tdir;
while tdir
*
t(iteT)<tdir*(tfin+tdir
*
incT)
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)
*
Deltat
*
tdir; % evaluate the solution at the shifted time
temp=tsub(iteT
−
1);
z(iteT,:)=nde(temp,c,tstart+t(iteT));
if tdir
*
t(iteT)>=tdir
*
tfin
*
0.99999 % repeat for the next sub-interval
break
else
iteT=iteT+1;  t(iteT)=tini+(iteT
−
1)
*
incT
*
tdir;
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:pnt
*
pntDelta)’*
incT; % evaluation of the initial solution
[
−
tau 0
]
z=history(t); % initial behavior in the inverval
[
−
tau 0
]
c=z(end,:); % initial condition
m=pnt
*
pntDelta; % samples between tau
itDlt=0;
while (t(end)+eps)<tfin;
tspan=
[
tini+itDlt
*
Deltat;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(pnt
*
pntDelta+1:end,:);t=t(pnt
*
pntDelta+1:end,:);
endFigure6 shows the relative error between our approximate EMPHM and the dde23 solution and its relationship with the order m and the discretized time intervals N. Notice that the relative error values coincide at values of N
≥
45. Also, we can see from Figure 6 that the computed relative error values for approximate solutions of order m
≥
5 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 ofm
= 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]:
(46)
x
¨
t
+
2
ζ
ω
n
x
˙
t
+
ω
n
2
x
t
=
-
a
p
K
s
t
m
m
x
t
-
x
t
-
τ
,
where ω
n is the angular natural frequency of the system, ζ is the damping ratio, a
p is the depth of cut, m
m is the modal mass of the tool, τ represents the time delay which is equal to the tooth passing period, and K
s
t is the specific cutting force coefficient which can be determined from
(47)
K
s
t
=
∑
j
=
1
z
n
g
ϕ
j
t
sin
ϕ
j
t
K
t
cos
⁡
ϕ
j
t
+
K
n
sin
ϕ
j
t
,
where z
n is the tool number of teeth, K
t and K
n are the tangential and the normal linear cutting force coefficients, respectively, ϕ
j
t is the angular position of the j-tooth defined as
(48)
ϕ
j
t
=
2
π
n
60
t
+
2
π
j
z
n
,
and n is the spindle speed in rpm [20]. The function g
φ
j
t is a switching function, which has a unity value when the j-tooth is cutting and zero otherwise:
(49)
g
ϕ
j
t
=
1
ϕ
st
<
ϕ
j
t
<
ϕ
ex
0
otherwise
.
Here, ϕ
st and ϕ
ex are the angles where the teeth enter and exit the workpiece. For upmilling, ϕ
st
=
0 and ϕ
ex
=
arccos
1
-
2
a
d, for downmilling, ϕ
st
=
arccos
2
a
d
-
1 and ϕ
ex
=
π, where a
d is the radial depth of cut ratio.By following the EMHPM procedure, we can write (46) in equivalent form as
(50)
x
¨
i
T
+
2
ζ
ω
n
x
˙
i
T
+
ω
n
2
x
i
T
≈
-
a
p
K
st
m
m
x
i
T
-
x
i
-
N
+
1
T
,
where x
i
T denotes the m order solution of (46) on the ith subinterval that satisfies the initial conditions x
i
0
=
x
i
-
1, x
˙
i
0
=
x
˙
i
-
1, and h
t
=
h
t and x
-
τ is given by (35). Introducing the transformation x
i
=
x
i
,
x
˙
i
T, (50) can be written as a system of first order linear delay differential equations of the form
(51)
x
˙
i
T
=
A
t
x
i
T
+
B
t
x
i
-
N
+
1
T
,
where
(52)
A
t
=
0
1
-
ω
n
2
-
w
m
m
K
st
-
2
ζ
ω
n
;
B
t
=
0
0
w
m
m
K
st
0
.
We next apply the EMHPM procedure to solve (46) by considering a downmilling operation with the following parameter values: z
n
=
2, a
d
=
0.1, ω
n
=
5793 rad⁄s, ζ
=
0.011, m
m
=
0.03993 kg, K
t
=
6
×
10
8 N/m2, and K
n
=
2
×
10
8 N/m2. As we can see from Figures 7 and 8 and for the depth of cut values of a
p
= 2 mm (stable) and a
p
= 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 a
p
=
2 mm, n
=
10000 rpm. Stable machine operation.Figure 8
EMHPM approximate solutions of (46) with parameter values of a
p
=
3 mm and n
=
10000 rpm. Unstable (chatter) machine operation.Figure9 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 N
=
75 since the average step size of the dde23 algorithm is around Δ
t
≈
τ
/
N. Note that, for m
=
7, 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 m
=
2, 7, and 10 and N
≥
20.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 (m = 7)
77
First order EMHPM (m = 7)
87Figure 9
Estimated relative error values between dde23 and the EMHPM approximate solutions. Here we have used the system parameter values ofa
p
=
2 mm and n
=
1000 rpm and values of m
= 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 ruleu
t
→
u
i
T, where u
i 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.



---

*Source: 382475-2015-03-02.xml*