# Euler-Maclaurin Method for Linear Differential Equations with Piecewise Constant Arguments with One Delay: Stability and Oscillations

Authors: Qi Wang; Jiechang Wen; Shenshan Qiu
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This paper focuses on the stability and oscillations of Euler-Maclaurin method for linear differential equations with piecewise constant argumentsu′(t)=au(t)+bu([t]). The necessary and sufficient conditions under which the numerical stability region contains the analytical stability region are given. Furthermore, the conditions of oscillation for the Euler-Maclaurin method are obtained. We prove that the Euler-Maclaurin method preserves the oscillations of the analytic solution. Moreover, the relationships between stability and oscillations are discussed for analytic solution and numerical solution, respectively. Finally, some numerical experiments for verifying the theoretical analysis are also provided.


---

## Body

## 1. Introduction

In the present paper we will consider the following differential equations with piecewise constant arguments (EPCA):(1)u′(t)=au(t)+bu([t]),t≥0,u(0)=u0,wherea,b,u0∈R and [·] denotes the greatest integer function. Moreover, we assume that b<0. The general form of this type  of  equations is
(2)u′(t)=f(t,u(t),u(β(t))),t≥0,u(0)=u0,where the argumentβ(t) has intervals of constancy.It is well known that studies of EPCA were motivated by the fact that they represent a hybrid of continuous and discrete dynamical systems and combine the properties of both the differential and difference equations. The local asymptotic stabilization of a class of instable autonomous nonlinear difference equations was studied in [1]. Convergence and stability of adaptation algorithm for a class of nominally first-order hybrid time-invariant linear systems were considered in [2], and stability for hybrid difference-differential systems was analyzed in [3]. The first work devoted to EPCA is the paper of Shah and Wiener [4] in 1983. In 1984 Cooke and Wiener [5] studied these equations with delay. The investigation of EPCA is an important subject because they play an important role in numerous applications [6–8]. In the literature, there are many papers dealing with the properties of EPCA, including the existence and uniqueness of solutions [9], the stability [10], the oscillations [11], and the existence of integral manifolds [12]. Moreover, the concept of EPCA has been generalized in [13–15] by introducing arbitrary piecewise constant functions as arguments. Significant parts of pioneer results for EPCA can be found in [16]. For more details of EPCA, the reader can see [17–22] and the references cited therein.In recent ten years, much work has been done in the numerical solution of EPCA. For instance, in [23, 24], the stability of numerical solution in Runge-Kutta methods for (1) and the mixed type EPCA was studied, respectively. In [25], the numerical stability of θ-methods was considered for EPCA of advanced type. For EPCA with variable coefficients, Song and Liu [26] studied the asymptotic stability of the analytic and numerical solutions. The oscillations of θ-methods and Runge-Kutta methods were investigated in [27, 28] for EPCA of retarded type, respectively. However, as far as we know, very few results were obtained on combining the stability with the oscillations of the numerical solution in the same article except for [29]. Different from [29], the novel idea of our paper is that we will study both stability and oscillations of the numerical solution in the Euler-Maclaurin method for (1), and their relationships are analyzed quantitatively.Euler-Maclaurin method is an important tool of numerical analysis which was discovered independently and almost simultaneously by Euler and Maclaurin in the first half of 18th century. Rota [30] called Euler-Maclaurin “one of the most remarkable formulas of mathematics.” After that, it shows us how to trade a finite sum for an integral. It works much like Taylor’s formula. The equation involves an infinite series that may be truncated at any point, leaving an error term that can be bounded.It is well known that the trapezoidal rule can be derived from the Euler-Maclaurin formula. See, for example, Munro’s paper [31]. In [32], the author indicated how the Newton-Cotes quadrature formulas and various other quadrature formulas can be developed from special cases of periodic Euler-Maclaurin formula. With respect to the numerical analysis of delay differential equations (DDEs), few results on the Euler-Maclaurin method were obtained except for [33]. In our present paper, we will apply this method to (1) and investigate numerical stability and oscillations.The rest of this paper is organized as follows. In Section2, some necessary notations and results are given. In Section 3, we propose Euler-Maclaurin method to solve (1), and the asymptotical stability of numerical solution is discussed. In Section 4, we investigate the oscillations of the numerical solution and prove that the oscillations of the analytic solution are inherited by the Euler-Maclaurin method. In Section 5, we establish some relationships between stability and oscillations, and several numerical experiments are shown in Section 6. Finally, concluding remarks are given in Section 7.


## 2. Preliminaries

In this section, we will focus our attention on some preliminary results which will be used in the further analysis.Proposition 1 (see [16]).
A solution of (1) on [0,∞) is a function u(t) which satisfies the following conditions:(i)
u(t)  is continuous on  [0,∞);(ii)
the derivative u′(t) exists at each point t∈[0,∞), with the possible exception of the points [t]∈[0,∞), where one-sided derivatives exist;(iii)
equation (1) is satisfied on each interval  [n,n+1)⊂[0,∞) with integral end points.Theorem 2 (see [16]).
Equation (1) has a unique solution on [0,∞):(3)u(t)=m0({t})b0[t]u0,
where {t} is the fractional part of t({t}=t-[t]) and
(4)b0=m0(1),m0(t)=eat+(eat-1)a-1b.Theorem 3 (see [16]).
The solution of (1) is asymptotically stable (u(t)→0 as t→∞) for all u0, if and only if
(5)-a(ea+1)ea-1<b<-a,a≠0,-2<b<0,a=0.Definition 4.
A non-trivial solution of (1) is said to be oscillatory if there exists a sequence {tk}k=1∞ such that  tk→∞ as k→∞ and u(tk)u(tk-1)<0; otherwise it is called non-oscillatory. We say (1) is oscillatory if all the non-trivial solutions of (1) are oscillatory; we say (1) is non-oscillatory if all the non-trivial solutions of (1) are non-oscillatory.Theorem 5 (see [16]).
Every solution of (1) is non-oscillatory if and only if
(6)b>aea1-ea.Remark 6.
By Definition4, we can see that if a solution u(t) of (1) is non-oscillatory and continuous, then it must be eventually positive or eventually negative. That is, there exists a T0∈R such that u(t)>0 or u(t)<0 for t≥T0.


## 3. Numerical Stability

### 3.1. Bernoulli’s Number and Bernoulli’s Polynomial

It is well known that(7)zez-1=∑j=0∞Bjj!zj,|z|<2π,zexzez-1=∑j=0∞Bj(x)j!zj,|z|<2π,whereBj and Bj(x), j=0,1,…,  are called Bernoulli’s number and Bernoulli’s polynomial, respectively.Proposition 7 (see [33]).
The Bernoulli’s numberBj satisfies the following properties:(i)
B0=1,B1=-1/2,(ii)
B2j=2(-1)j+1(2j)!∑k=1∞(2kπ)-2j,B2j+1=0,j≥1,(iii)
(B2(2i-1)/(2(2i-1))!)+(B2(2i)/(2(2i))!)>0,
i=1,2,…,(iv)
(B2(2i+1)/(2(2i+1))!)+(B2(2i)/(2(2i))!)<0,i=1,2,….Proposition 8 (see [33]).
The Bernoulli’s polynomialBj(x) satisfies the following properties:(i)
B0(x)=1,B1(x)=x-1/2, Bk(x)=∑j=0kCkjBjxk-j,(ii)
B2k+1(1)=B2k+1(0)=B2k+1=0,(iii)
B2k(1)=B2k(0)=B2k,(iv)
Bk(x)=(1/(k+1))Bk+1′(x),k=1,2,….


### 3.2. The Euler-Maclaurin Method and Discretization

Letq>0 be an integer, and assume that the function f(t) is at least 2q+2  times continuously differentiable on [c,d]. Further assume that h evenly divides c and d; then Atkinson’s version of the Euler-Maclaurin formula [34] is as follows:
(8)∫cdf(t)dt=h∑n=cdf(n)-h2(f(c)+f(d))∫cdf(t)dt=-∑i=1qh2iB2i(2i)!(D2i-1f(d)-D2i-1f(c)),whereBj denotes the jth Bernoulli number and D denotes the differentiation operator.Leth be a given stepsize; m>1 is a given integer and satisfies h=1/m. The grid points ti are defined by ti=ih(i=0,1,2,…); let i=km+l(l=0,1,…,m-1). Applying formula (8) to (1), we have
(9)ui+1=ui+ha2(ui+1+ui)-∑j=1nB2j(ha)2j(2j)!(ui+1-ui)+hbukm.Lemma 9 (see [35]).
Assume thatf(x) has 2n+3rd continuous derivative on the interval [ti,ti+1]; then we have
(10)|∫titi+1f(t)dt-h2(f(ti+1)+f(ti))+∑j=1nB2jh2j(2j)!(f(2j-1)(ti+1)-f(2j-1)(ti))+∑j=1nB2jh2j(2j)!(f(2j-1)(ti+1)-f(2j-1)(ti))|=O(h2n+3).According to Theorem  2.2 in [33], we obtain the following theorem for convergence.Theorem 10.
For any givenn∈N, the Euler-Maclaurin method is of order 2n+2.
Equation (9) can be written as
(11)ukm+l+1=α(x)ukm+l+ba(α(x)-1)ukm,whena≠0,ukm+l+1=ukm+l+hbukm,whena=0,
where(12)x=ha,α(x)=1+xΦ(x),(13)Φ(x)=1-x2+∑j=1nB2jx2j(2j)!.It is not difficult to see that (11) is equivalent to
(14)u(k+1)m=(α(x)m+ba(α(x)m-1))ukm,whena≠0,u(k+1)m=(1+b)ukm,whena=0,(15)ukm+l=(α(x)l+ba(α(x)l-1))ukm,whena≠0,ukm+l=(1+hbl)ukm,whena=0,where0≤l≤m-1.


### 3.3. Stability Analysis

Definition 11.
Process (9) for (1) is called asymptotically stable at (a,b) for all sufficiently small h if and only if there exists a constant h0, such that for any given u0, relation (15) defines ui that satisfies ui→0 as i→∞ whenever 0<h<h0.Definition 12.
The set, denoted byS, of all points (a,b) at which process (9) for (1) is asymptotically stable, for all sufficiently small h, is called the asymptotical stability region.Lemma 13 (see [33]).
If|x|≤1, then φ(x)≥1/2 for x>0 and φ(x)≥1 for x≤0.By relation (15), we can easily prove the following lemma.Lemma 14.
For anyk and 0≤l≤m-1, there exists a constant C>0 independent of k and l such that
(16)|ukm+l|≤C|ukm|.So we have the following corollary.Corollary 15.
u
i
→
0 as i→∞ if and only if ukm→0 as k→∞.According to Definition11, Corollary 15, and relation (14), the following theorem is obtained.Theorem 16.
The numerical solution of (1) is asymptotically stable if and only if
(17)-a(α(x)m+1)α(x)m-1<b<-a,whena≠0,-2<b<0,whena=0.For convenience, we divideS into two parts S=S*∪S0, where
(18)S*={(a,b)∈S:a≠0},S0={(a,b)∈S∖S*:a=0},therefore, by Definition12 and Theorem 16, we have
(19)S*={(a,b)∈R2:-a(α(x)m+1)α(x)m-1<b<-a,a≠0},S0={(0,b)∈R2:-2<b<0}.Let us introduce the following sets of all points(a,b)∈R2 which satisfy condition (5) by H; in a similar way, we denote H=H*∪H0, where
(20)H*={(a,b)∈R2:-a(ea+1)ea-1<b<-a,a≠0},H0={(0,b)∈R2:-2<b<0}.We will study which conditions lead toH*⊆S* and H0⊆S0.In the following, we assumeM>|a|, which implies that |x|<1 for the stepsize h=1/m with m≥M. The following lemma will be useful to determine the stability and oscillatory conditions.Lemma 17 (see [33]).
If|x|≤1, then
(21)Φ(x)≤xex-1
whenn is even and
(22)Φ(x)≥xex-1
whenn is odd.Then the first main result for stability is as follows.Theorem 18.
For the Euler-Maclaurin method,H*⊆S* if and only if n is odd.Proof.
In view of (19) and (20), we know that H*⊆S* if and only if
(23)-a(α(x)m+1)α(x)m-1≤-a(ea+1)ea-1,
that is,(24)aα(x)m-1≥aea-1,
which is equivalent to(25)α(x)≤ex,a>0,α(x)≥ex,a<0.
In view ofα(x)=1+(x/Φ(x)), thus (25) leads to
(26)Φ(x)≥xex-1.
As a consequence of Lemma17, the proof is complete.Theorem 19.
For all the Euler-Maclaurin method, we haveH0=S0.


## 3.1. Bernoulli’s Number and Bernoulli’s Polynomial

It is well known that(7)zez-1=∑j=0∞Bjj!zj,|z|<2π,zexzez-1=∑j=0∞Bj(x)j!zj,|z|<2π,whereBj and Bj(x), j=0,1,…,  are called Bernoulli’s number and Bernoulli’s polynomial, respectively.Proposition 7 (see [33]).
The Bernoulli’s numberBj satisfies the following properties:(i)
B0=1,B1=-1/2,(ii)
B2j=2(-1)j+1(2j)!∑k=1∞(2kπ)-2j,B2j+1=0,j≥1,(iii)
(B2(2i-1)/(2(2i-1))!)+(B2(2i)/(2(2i))!)>0,
i=1,2,…,(iv)
(B2(2i+1)/(2(2i+1))!)+(B2(2i)/(2(2i))!)<0,i=1,2,….Proposition 8 (see [33]).
The Bernoulli’s polynomialBj(x) satisfies the following properties:(i)
B0(x)=1,B1(x)=x-1/2, Bk(x)=∑j=0kCkjBjxk-j,(ii)
B2k+1(1)=B2k+1(0)=B2k+1=0,(iii)
B2k(1)=B2k(0)=B2k,(iv)
Bk(x)=(1/(k+1))Bk+1′(x),k=1,2,….


## 3.2. The Euler-Maclaurin Method and Discretization

Letq>0 be an integer, and assume that the function f(t) is at least 2q+2  times continuously differentiable on [c,d]. Further assume that h evenly divides c and d; then Atkinson’s version of the Euler-Maclaurin formula [34] is as follows:
(8)∫cdf(t)dt=h∑n=cdf(n)-h2(f(c)+f(d))∫cdf(t)dt=-∑i=1qh2iB2i(2i)!(D2i-1f(d)-D2i-1f(c)),whereBj denotes the jth Bernoulli number and D denotes the differentiation operator.Leth be a given stepsize; m>1 is a given integer and satisfies h=1/m. The grid points ti are defined by ti=ih(i=0,1,2,…); let i=km+l(l=0,1,…,m-1). Applying formula (8) to (1), we have
(9)ui+1=ui+ha2(ui+1+ui)-∑j=1nB2j(ha)2j(2j)!(ui+1-ui)+hbukm.Lemma 9 (see [35]).
Assume thatf(x) has 2n+3rd continuous derivative on the interval [ti,ti+1]; then we have
(10)|∫titi+1f(t)dt-h2(f(ti+1)+f(ti))+∑j=1nB2jh2j(2j)!(f(2j-1)(ti+1)-f(2j-1)(ti))+∑j=1nB2jh2j(2j)!(f(2j-1)(ti+1)-f(2j-1)(ti))|=O(h2n+3).According to Theorem  2.2 in [33], we obtain the following theorem for convergence.Theorem 10.
For any givenn∈N, the Euler-Maclaurin method is of order 2n+2.
Equation (9) can be written as
(11)ukm+l+1=α(x)ukm+l+ba(α(x)-1)ukm,whena≠0,ukm+l+1=ukm+l+hbukm,whena=0,
where(12)x=ha,α(x)=1+xΦ(x),(13)Φ(x)=1-x2+∑j=1nB2jx2j(2j)!.It is not difficult to see that (11) is equivalent to
(14)u(k+1)m=(α(x)m+ba(α(x)m-1))ukm,whena≠0,u(k+1)m=(1+b)ukm,whena=0,(15)ukm+l=(α(x)l+ba(α(x)l-1))ukm,whena≠0,ukm+l=(1+hbl)ukm,whena=0,where0≤l≤m-1.


## 3.3. Stability Analysis

Definition 11.
Process (9) for (1) is called asymptotically stable at (a,b) for all sufficiently small h if and only if there exists a constant h0, such that for any given u0, relation (15) defines ui that satisfies ui→0 as i→∞ whenever 0<h<h0.Definition 12.
The set, denoted byS, of all points (a,b) at which process (9) for (1) is asymptotically stable, for all sufficiently small h, is called the asymptotical stability region.Lemma 13 (see [33]).
If|x|≤1, then φ(x)≥1/2 for x>0 and φ(x)≥1 for x≤0.By relation (15), we can easily prove the following lemma.Lemma 14.
For anyk and 0≤l≤m-1, there exists a constant C>0 independent of k and l such that
(16)|ukm+l|≤C|ukm|.So we have the following corollary.Corollary 15.
u
i
→
0 as i→∞ if and only if ukm→0 as k→∞.According to Definition11, Corollary 15, and relation (14), the following theorem is obtained.Theorem 16.
The numerical solution of (1) is asymptotically stable if and only if
(17)-a(α(x)m+1)α(x)m-1<b<-a,whena≠0,-2<b<0,whena=0.For convenience, we divideS into two parts S=S*∪S0, where
(18)S*={(a,b)∈S:a≠0},S0={(a,b)∈S∖S*:a=0},therefore, by Definition12 and Theorem 16, we have
(19)S*={(a,b)∈R2:-a(α(x)m+1)α(x)m-1<b<-a,a≠0},S0={(0,b)∈R2:-2<b<0}.Let us introduce the following sets of all points(a,b)∈R2 which satisfy condition (5) by H; in a similar way, we denote H=H*∪H0, where
(20)H*={(a,b)∈R2:-a(ea+1)ea-1<b<-a,a≠0},H0={(0,b)∈R2:-2<b<0}.We will study which conditions lead toH*⊆S* and H0⊆S0.In the following, we assumeM>|a|, which implies that |x|<1 for the stepsize h=1/m with m≥M. The following lemma will be useful to determine the stability and oscillatory conditions.Lemma 17 (see [33]).
If|x|≤1, then
(21)Φ(x)≤xex-1
whenn is even and
(22)Φ(x)≥xex-1
whenn is odd.Then the first main result for stability is as follows.Theorem 18.
For the Euler-Maclaurin method,H*⊆S* if and only if n is odd.Proof.
In view of (19) and (20), we know that H*⊆S* if and only if
(23)-a(α(x)m+1)α(x)m-1≤-a(ea+1)ea-1,
that is,(24)aα(x)m-1≥aea-1,
which is equivalent to(25)α(x)≤ex,a>0,α(x)≥ex,a<0.
In view ofα(x)=1+(x/Φ(x)), thus (25) leads to
(26)Φ(x)≥xex-1.
As a consequence of Lemma17, the proof is complete.Theorem 19.
For all the Euler-Maclaurin method, we haveH0=S0.


## 4. Oscillations Analysis

Similar to Definition4, we can easily get the definitions of oscillations and nonoscillations of numerical solution; for the sake of simplicity, we omit them. Then the relationships of the nonoscillations between the integer nodes and any nodes are as follows.Theorem 20.
{
u
i
} and {ukm} are given by (15) and (14), respectively; then {ui} is non-oscillatory if and only if {ukm} is non-oscillatory.Proof.
The necessity is easy to prove, so we only consider the sufficiency. If{ukm} is non-oscillatory, without loss of generality, we can assume that {ukm} is an eventually negative solution of (14); that is, there exists a K0∈R such that ukm<0 for k>K0. We will prove ukm+l<0 for all k>K0+1 and l=0,1,…,m-1.
Ifa<0, then 0<α(x)<1 and α(x)m≤α(x)l; hence
(27)ukm+l=(α(x)l+ba(α(x)l-1))ukmukm+l≤(α(x)m+ba(α(x)m-1))ukmukm+l=ukm+m<0.
Ifa>0, then α(x)>1 and α(x)-m≤α(x)-l; thus
(28)α(x)-lukm+l=(1+ba(1-α(x)-l))ukmα(x)-lukm+l≤(1+ba(1-α(x)-m))ukmα(x)-lukm+l=α(x)-mukm+m<0.
Thereforeukm+l<0. This completes the proof.From Theorem20, we get the following corollary.Corollary 21.
{
u
i
} and {ukm} are given by (15) and (14), respectively; then {ui} is oscillatory if and only if {ukm} is oscillatory.Theorem 22.
Equation (14) is oscillatory if and only if
(29)b<aα(x)m1-α(x)m,whena≠0,b<-1,whena=0.Proof.
Equation (14) is oscillatory if and only if the corresponding characteristic equation has no positive roots, that is,
(30)α(x)m+ba(α(x)m-1)<0,whena≠0,1+b<0,whena=0,
which is equivalent to(31)b<aα(x)m1-α(x)m,whena≠0,(32)b<-1,whena=0,
thus (29) is obtained. The proof is finished.Let(33)Π=aea1-ea,Π(m)=aα(x)m1-α(x)m,we have the following lemma.Lemma 23.
(I)   Π(m)→Π as h→0;
(II)   Π≥Π(m) if either of the following conditions is satisfied:
(i)   a>0 and ex≥α(x),
(ii)   a<0 and ex≤α(x);
(
III
)
Π
<
Π
(
m
) if either of the following conditions is satisfied:
(i)  a>0 and ex<α(x),
(ii)  a<0 and ex>α(x).Proof.
(II)   Ifa>0 and ex≥α(x), then ea≥α(x)m, which is equivalent to
(34)11-ea≥11-α(x)m,that is,
(35)ea1-ea≥α(x)m1-α(x)m,
so we haveΠ≥Π(m). Similarly, we can prove the other cases.Definition 24.
We say the Euler-Maclaurin method preserves the oscillations of (1) if (1) oscillates, then there is an  h0>0, such that (15) oscillates for h<h0. Similarly, we say the Euler-Maclaurin method preserves the nonoscillations of (1) if (1) nonoscillates then there is an h0>0, such that (15) non-oscillates for h<h0.By Theorems5, 20, and 22 and Corollary 21, we have the following theorem.Theorem 25.
(i)  The Euler-Maclaurin method preserves the oscillations of (1) if and only if  Π<Π(m).
(ii)  The Euler-Maclaurin method preserves the non-oscillations of (1) if and only if  Π≥Π(m).Applying Theorem25 and Lemmas 17 and 23, we can obtain the following main theorems for oscillations.Theorem 26.
The Euler-Maclaurin method preserves the oscillations of (1) if n is even.Theorem 27.
The Euler-Maclaurin method preserves the non-oscillations of (1) if n is odd.


## 5. Relationships between Stability and Oscillations

Let(36)Π~=-a(ea+1)ea-1,Π~(m)=-a(α(x)m+1)α(x)m-1.By Theorems3, 5, 16, and 22, we can get the following other main theorems.Theorem 28.
Whena≠0, the analytic solution of (1) is(i)
oscillatory and unstable ifb∈(-∞,Π~),(ii)
oscillatory and asymptotically stable ifb∈(Π~,Π),(iii)
non-oscillatory and asymptotically stable ifb∈(Π,-a),(iv)
non-oscillatory and unstable ifb∈(-a,0).Theorem 29.
Whena≠0, the numerical solution of (1) is(i)
oscillatory and unstable ifb∈(-∞,Π~(m)),(ii)
oscillatory and asymptotically stable ifb∈(Π~(m),Π(m)),(iii)
non-oscillatory and asymptotically stable ifb∈(Π(m),-a),(iv)
non-oscillatory and unstable ifb∈(-a,0).Theorem 30.
Whena=0, the analytic solution and the numerical solution of (1) both are(i)
oscillatory and unstable ifb∈(-∞,-2),(ii)
oscillatory and asymptotically stable ifb∈(-2,-1),(iii)
non-oscillatory and asymptotically stable ifb∈(-1,0).Remark 31.
In Theorem30, we can see that the solution of (1) is not non-oscillatory and unstable at the same time. Otherwise, if the solution of (1) is both non-oscillatory and unstable, we can get b∈(0,+∞); it is inconsistent with b<0.


## 6. Numerical Experiments

In this section, we will use the following four equations to demonstrate the main theorems(37)u′(t)=u(t)-1.5u([t]),u(0)=1,(38)u′(t)=-2u(t)-u([t]),u(0)=1,(39)u′(t)=3u(t)-3.2u([t]),u(0)=1,(40)u′(t)=0.8u(t)-0.1u([t]),u(0)=1.For (37) and (38), it is not difficult to see that (1,-1.5),(-2,-1)∈H by straightforward computing. Let the stepsize h=1/m; we shall use the Euler-Maclaurin method with n=2, the θ-methods with θ=0.5, and the 2-Radau IA method to get the numerical solution at t=10. On the other hand, the theoretical solutions are u(10)≈3.0750×10-9 and u(10)≈5.3398×10-6 for (37) and (38), respectively. In Tables 1 and 2, we have listed the absolute errors (AE), the relative errors (RE) at t=10, and the ratio of the errors of the case m=20 over that of m=40. We can see from these tables that the Euler-Maclaurin method with n=2 is of order 6; that is, the method preserves its order of convergence. Furthermore, the errors of the θ-methods and the 2-Radau IA method are both larger than that of the Euler-Maclaurin method. Therefore, compared with the two mentioned numerical methods, the Euler-Maclaurin method has higher accuracy.Table 1
Errors and ratio for (37).
Euler-Maclaurin method
θ-methods
2-Radau IA method
AE
RE
AE
RE
AE
RE
m
=
2
1.5555
E
-
13
5.0586
E
-
5
2.7882
E
-
9
9.0670
E
-
1
6.5320
E
-
10
2.1240
E
-
1
m
=
3
1.3547
E
-
14
4.4054
E
-
6
1.8944
E
-
9
6.1610
E
-
1
1.7228
E
-
10
5.6000
E
-
2
m
=
5
6.2942
E
-
16
2.0469
E
-
7
8.6318
E
-
10
2.8070
E
-
1
3.5051
E
-
11
1.1400
E
-
2
m
=
10
9.8172
E
-
18
3.1926
E
-
9
2.3894
E
-
10
7.7700
E
-
2
4.2380
E
-
12
1.4000
E
-
3
m
=
20
1.5337
E
-
19
4.9876
E
-
11
6.1287
E
-
11
1.9900
E
-
2
5.2217
E
-
13
1.6981
E
-
4
m
=
40
2.7140
E
-
21
8.8259
E
-
13
1.5420
E
-
11
5.0000
E
-
3
6.4824
E
-
14
2.1081
E
-
5
Ratio
56.5107
56.5110
3.9745
3.9800
8.0552
8.0551Table 2
Errors and ratio for (38).
Euler-Maclaurin method
θ-methods
2-Radau IA method
AE
RE
AE
RE
AE
RE
m
=
2
2.5585
E
-
9
4.7914
E
-
4
1.1595
E
-
5
2.1715
E
+
0
8.9866
E
-
7
1.6830
E
-
1
m
=
3
2.1747
E
-
10
4.0727
E
-
5
3.5420
E
-
6
6.6330
E
-
1
2.6576
E
-
7
4.9800
E
-
2
m
=
5
9.9800
E
-
12
1.8690
E
-
6
1.0695
E
-
6
2.0030
E
-
1
5.9321
E
-
8
1.1100
E
-
2
m
=
10
1.5485
E
-
13
2.9000
E
-
8
2.4905
E
-
7
4.6600
E
-
2
7.7177
E
-
9
1.4000
E
-
3
m
=
20
2.4153
E
-
15
4.5232
E
-
10
6.1184
E
-
8
1.1500
E
-
2
9.8799
E
-
10
1.8502
E
-
4
m
=
40
3.7816
E
-
17
7.0819
E
-
12
1.5230
E
-
8
2.9000
E
-
3
1.2508
E
-
10
2.3424
E
-
5
Ratio
63.8698
63.8699
4.0173
3.9655
7.8989
7.8987In Figures1 and 3, we draw the numerical solution with m=20 for (37) and (38), respectively. It is easy to see that the numerical solution is asymptotically stable. In Figure 2, we draw the stability region for (37); it can be seen that the analytical stability region is contained in the numerical stability region, which is in accordance with Theorem 18.Figure 1
The numerical solution of (37) with n=2 and m=20.The stability region of (37): (a) the analytical stability region, (b) the numerical stability region with n=3 and m=20.
(a)
(b)Figure 3
The numerical solution of (38) with n=3 and m=20.For (39) and (40), the analytic solution of (39) is oscillatory, and the analytic solution of (40) is non-oscillatory according to Theorem 5. From Figures 4 and 5, we can see that the numerical solution of (39) is oscillatory and the numerical solution of (40) is non-oscillatory, which are in accordance with Theorems 26 and 27.The analytic solution and the numerical solution of (39) with n=2 and m=40.
(a)
(b)The analytic solution and the numerical solution of (40) with n=3 and m=40.
(a)
(b)In Figure4, let m=40; it is easy to calculate that Π≈-3.1572, Π~≈-3.3144, Π(m)≈-3.1572, and Π~(m)≈-3.3144. Obviously, b=-3.2∈(Π~,Π) and b=-3.2∈(Π~(m),Π(m)). Therefore, the analytic solution and the numerical solution of (39) are both oscillatory and asymptotically stable according to Theorems 28 and 29, which are in agreement with Figure 4. For (37), (38), and (40), we can verify them in the same way (see Figures 1, 3, and 5).From Figures1–5, we can see that the numerical results coincide with our theoretical analysis.


## 7. Conclusions

In this paper, the Euler-Maclaurin method is presented for solving an important EPCA. Asymptotical stability and oscillations of numerical solution are considered, respectively. Furthermore, the relationships between stability and oscillations are also revealed. Results of numerical examples have shown that the Euler-Maclaurin method has good convergence and stability, so it is applicable for the study of EPCA. That is, the Euler-Maclaurin method is an alternative numerical method for DDEs except for the classicalθ-methods and Runge-Kutta methods. We will consider the nonlinear and multidimensional problems in the future work.



---

*Source: 232484-2013-05-09.xml*