# Analysis of Exponential Runge–Kutta Methods for Differential Equations with Time Delay

Authors: Rui Zhan
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Numerous mathematical models simulating the phenomenon in science and engineering use delay differential equations. In this paper, we focus on the semilinear delay differential equations, which include a wide range of mathematical models with time lags, such as reaction-diffusion equation with delay, model of bacteriophage predation on bacteria in a chemostat, and so on. This paper is concerned with the stability and convergence properties of exponential Runge–Kutta methods for semilinear delay differential equations. GDN-stability and D-convergence of exponential Runge–Kutta methods are investigated. These two concepts are generalizations of the classical AN-stability and B-convergence for ordinary differential equations to delay differential equations. Sufficient conditions for GDN-stability are given by a newly introduced concept of strong exponential algebraic stability. Further, with the aid of diagonal stability, we show that exponential Runge–Kutta methods are D-convergent. The D-convergent orders are also examined. Numerical experiments are presented to illustrate the theoretical results.


---

## Body

## 1. Introduction

Mathematical models are valuable tools for explaining the observed system or predicting future events. When simulating the phenomenon in chemistry, biological system, medicine, engineering control systems, and climate models, we need to consider the effects caused by the time delay. In these mathematical models, the change of the unknown quantity depends not only on the present state but also on the past state. The resulting equations are called by delay differential equations (DDEs). DDEs provide more realistic models by replacing the point-wise assumptions of traditional differential equations with distributed assumptions.For instance, the following DDE system is studied in [1] as a mathematical model for bacteriophage predation on bacteria in a chemostat:(1)dRdt=R0−Rω−εNAψARKA+R−εNBψBRKB+R,dNAdt=NAψARKA+R−αANAP−ωNA,dNBdt=NBψBRKB+R−αBNBP−ωNB,dPdt=βe−τωαANAt−τPt−τ+βe−τωαBNBt−τPt−τ−αANAP−αBNBP−ωP,where R,NA,NB,P denote the resource supporting bacterial growth, uninfected bacteria, phage infected bacteria, and phage, respectively. τ is the latent period of the phage. Note that the delay terms NAt−τ, NBt−τ, and Pt−τ increase the model’s suitability for simulating complicated real systems, but they also make it harder to analyze the equations..As another example, the following reaction-diffusion equations with time delay are investigated in [2]:(2)∂ux,t∂t=fux,t,ux,t−τ+A∂2ux,t∂x2+vx,t,where ux,t stands for the state, f is a nonlinear function, vx,t denotes the distributed input, and τ is the time delay. The authors considered the boundary input and distributed input in the model and presented a verified criterion for exponential input-to-state stability.Moreover, two delay differential equation models are established in [3] by considering the matrix game models between public goods suppliers and between the suppliers and the government. A two-delay HIV-1 virus model with delay-dependent parameters is considered in [4].Equations mentioned above can be recast into the following abstract semilinear form:(3)y′t=Lyt+ft,yt,yt−τ,where L is a linear operator. With regard to numerical methods for DDE (3), due to the similarity in form of DDEs and ordinary differential equations (ODEs), many numerical methods originally designed for ODEs have been attempted to apply to DDEs. The reader is referred to the monograph in [5] for a comprehensive review.On the other hand, exponential integrators have attracted a lot of attention in recent years. Along with the increase of computational efficiency, there has been much work on constructing and analyzing exponential integrators for semilinear stiff ODEs. Various exponential integrators have been studied such as exponential Runge–Kutta methods [6], exponential multistep methods [7], exponential Rosenbrock methods [8], parallel exponential Rosenbrock methods [9], exponential Taylor methods [10], and exponential general linear methods [11]. A good exposition of exponential integrators can be found in [12].However, as noted by pioneers in the field of numerical methods for DDEs, the stability properties and accuracy of numerical methods for ODEs may not be preserved if one naively applies these methods to DDEs. According to the observation in [5], above phenomena can be classified into order failure and stability failure. They concluded that integration of DDEs requires specifically designed methods other than the standard ODE methods, and the stability and convergence properties of the numerical methods should also be examined carefully.This paper is concerned with the stability and convergence property of exponential Runge–Kutta methods for semilinear stiff delay differential equations. There are various types of stability concepts considered in the literature. As the generalizations of the classical A-stability for ODEs to DDEs, P-stability and GP-stability of exponential Runge–Kutta methods of collocation type have been studied for linear autonomous DDE in [13]. For the linear nonautonomous test equation, we considered the PN-stability and GPN-stability of the second-order Magnus integrator [14]. These kinds of concepts are the generalizations of the AN-stability to DDEs. In this paper, we will focus on the GDN-stability which can be viewed as the generalization of the BN-stability to DDEs. Moreover, on the convergence analysis, D-convergence concept was proposed in [15] for the Runge–Kutta method. Note that the D-convergence concept is the generalization of the B-convergence. In this paper, we will follow the work in [16] to investigate the D-convergence property of exponential Runge–Kutta methods for DDEs.It is worth to mention that differential equations are widely used to describe practical problems in science and engineering. Abbas et al. [17] proposed a model for heat transfer and fluid flow over an inclined moving surface. The magnetohydrodynamic effects on heat and mass transfer phenomena in third-grade fluid were studied in [18]. Meanwhile, the artificial neural network models have also been introduced in many engineering applications. Shafiq et al. [19] developed this kind of model to predict the boundary layer flow of a single-walled carbon nanotube nanofluid. A comparative study of artificial neural network versus parametric method in COVID-19 data analysis was carried out in [20].The paper is outlined as follows. In Section2, the stability property of the exact solution to the problem we considered is investigated and exponential Runge–Kutta methods (ERKMs) are described. In Section 3, in order to build sufficient conditions for GDN-stability, we introduce the concept of strong exponential algebraic stability. We show that under this condition, the ERKMs are GDN-stable. Section 4 is devoted to the D-convergence of ERKMs. With the aid of strong exponential algebraic stability and diagonal stability, we prove that strongly exponentially algebraically stable and diagonally stable ERKMs with stage order p, together with a Lagrangian interpolation of order q, are D-convergent of order minp,q+1. Some numerical experiments are presented in Section 5 to illustrate the theoretical results.


## 2. Stability Property of the Problem and the Numerical Scheme

Consider the following semilinear delay differential equations:(4)y′t=Lyt+ft,yt,yt−τ,t∈t0,T,yt=ϕt,t∈t0−τ,t0,where constant τ is the delay term, L is a complex-valued N×N matrix, and f:t0,T×ℂN×ℂN⟶ℂN and ϕ:t0−τ,t0⟶ℂN are continuous functions.Let⋅,⋅ be the usual inner product on ℂN. The induced norm and logarithmic norm are denoted by ⋅ and μ⋅, respectively. In this paper, we assume that there are constants σ,γ such that(5)ℜft,y1,z−ft,y2,z,y1−y2≤σy1−y22,(6)ft,y,z1−ft,y,z2≤γz1−z2,(7)0≤γ≤−σ,and the logarithmic norm of matrix L, which is defined by(8)μL=supx≠0ℜLx,xx2,satisfies(9)μL≤0.In the following, problem (4) satisfying conditions (5)–(9) is denoted by Pμ,σ,γ. We further assume that the exact solution yt is sufficiently smooth in the means that(10)diytdti≤Mi,holds uniformly for t∈t0−τ,T.


### 2.1. Stability of Exact Solution

We begin with the stability analysis of problemPμ,σ,γ. The following stability result for ODEs will be used in the subsequent sections.Lemma 1.
(see [21]). Consider the following linear initial value problem:(11)y′t=atyt+ft,t≥t0,yt0=y0,where y,a,f:t0,+∞⟶ℂ. If ℜat<0 for t≥t0, then the exact solution yt fulfills(12)yt≤maxy0;maxt0≤x≤tfx−ℜax.
Now we state the contractive property of exact solution of problemPμ,σ,γ.Lemma 2.
(see [21]). Let yt,zt denote the exact solutions of the problem Pμ,σ,γ with different initial functions ϕt,ψt, respectively. Then, for the difference of yt and zt, we have(13)yt−zt≤maxt0−τ≤x≤t0ϕx−ψx,t≥t0.


### 2.2. Numerical Scheme

Noting that problem (4) is semilinear, we can express its exact solution by applying the variation-of-constants formula(14)yt=et−t0Lyt0+∫t0tet−ηLfη,yη,yη−τdη.Now we construct the exponential Runge–Kutta methods based on this formula. The linear part is solved exactly, which can avoid the stiffness from the linear part. Then, approximating the nonlinear partf in the integral yields the numerical methods. Specifically, a s−stage exponential Runge–Kutta method with coefficients aij,bi,ci can be given by the usual Butcher tableau form. Here aij,bi,i,j=1,2,…,s are linear combinations of φ− functions, which are defined as follows.Definition 1.
For any complex numbersz, let φ0z=ez. For all integers k≥1, define φkz as(15)φkz=∫01e1−szsk−1j−1!ds.
Note that we have the following recurrence relation:(16)φj+1z=φjz−φj0z.
Applyings−stage exponential Runge–Kutta methods aij,bi,ci with step size h to problem Pμ,σ,γ yields the scheme(17)yn+1=ehLyn+h∑i=1sbihLftni,yin,y¯in,(18)yin=ecihLyn+h∑j=1saijhLftnj,yjn,y¯jn,j=1,2,…,s,where tni=tn+cih and yn, yin, and y¯in are approximations of ytn, ytni, and ytni−τ, respectively. For convenience, aijhL and bihL are abbreviated as aij and bi.
We note that, in contrast to the classic Runge–Kutta methods which solve (4) directly, exponential Runge–Kutta methods are based on the formula (14). For stiff differential equation, an excessively small time step is needed for classic Runge–Kutta methods. However, by computing the exponential term exactly or with high-order stable algorithm, the stiffness of the equation is handled analytically. Hence, it is reasonable to expect that exponential Runge–Kutta methods will be stable and efficient for stiff equation.
The delay termy¯in can be given by the initial function if tn+cih−τ≤t0. Otherwise, we approximate y¯in by the interpolation based on previous inner stage values. Let integer m and η∈0,1 be such that τ=m−ηh. Recall that the Lagrange basis functions for uniform mesh have the form(19)ℒℓη=∏k=−rk≠ℓνη−kℓ−k,ℓ=−r,−r+1,…,ν,where r,ν are nonnegative integers. In this setting, we compute y¯in by(20)y¯in=φtn+cih−τ,tn+cih−τ≤t0,∑ℓ=−rνℒℓηyin−m+ℓ,tn+cih−τ>t0.
Note that we should takeν<m to ensure that no unknown values yin−m+ℓ are involved in the interpolation.


## 2.1. Stability of Exact Solution

We begin with the stability analysis of problemPμ,σ,γ. The following stability result for ODEs will be used in the subsequent sections.Lemma 1.
(see [21]). Consider the following linear initial value problem:(11)y′t=atyt+ft,t≥t0,yt0=y0,where y,a,f:t0,+∞⟶ℂ. If ℜat<0 for t≥t0, then the exact solution yt fulfills(12)yt≤maxy0;maxt0≤x≤tfx−ℜax.
Now we state the contractive property of exact solution of problemPμ,σ,γ.Lemma 2.
(see [21]). Let yt,zt denote the exact solutions of the problem Pμ,σ,γ with different initial functions ϕt,ψt, respectively. Then, for the difference of yt and zt, we have(13)yt−zt≤maxt0−τ≤x≤t0ϕx−ψx,t≥t0.


## 2.2. Numerical Scheme

Noting that problem (4) is semilinear, we can express its exact solution by applying the variation-of-constants formula(14)yt=et−t0Lyt0+∫t0tet−ηLfη,yη,yη−τdη.Now we construct the exponential Runge–Kutta methods based on this formula. The linear part is solved exactly, which can avoid the stiffness from the linear part. Then, approximating the nonlinear partf in the integral yields the numerical methods. Specifically, a s−stage exponential Runge–Kutta method with coefficients aij,bi,ci can be given by the usual Butcher tableau form. Here aij,bi,i,j=1,2,…,s are linear combinations of φ− functions, which are defined as follows.Definition 1.
For any complex numbersz, let φ0z=ez. For all integers k≥1, define φkz as(15)φkz=∫01e1−szsk−1j−1!ds.
Note that we have the following recurrence relation:(16)φj+1z=φjz−φj0z.
Applyings−stage exponential Runge–Kutta methods aij,bi,ci with step size h to problem Pμ,σ,γ yields the scheme(17)yn+1=ehLyn+h∑i=1sbihLftni,yin,y¯in,(18)yin=ecihLyn+h∑j=1saijhLftnj,yjn,y¯jn,j=1,2,…,s,where tni=tn+cih and yn, yin, and y¯in are approximations of ytn, ytni, and ytni−τ, respectively. For convenience, aijhL and bihL are abbreviated as aij and bi.
We note that, in contrast to the classic Runge–Kutta methods which solve (4) directly, exponential Runge–Kutta methods are based on the formula (14). For stiff differential equation, an excessively small time step is needed for classic Runge–Kutta methods. However, by computing the exponential term exactly or with high-order stable algorithm, the stiffness of the equation is handled analytically. Hence, it is reasonable to expect that exponential Runge–Kutta methods will be stable and efficient for stiff equation.
The delay termy¯in can be given by the initial function if tn+cih−τ≤t0. Otherwise, we approximate y¯in by the interpolation based on previous inner stage values. Let integer m and η∈0,1 be such that τ=m−ηh. Recall that the Lagrange basis functions for uniform mesh have the form(19)ℒℓη=∏k=−rk≠ℓνη−kℓ−k,ℓ=−r,−r+1,…,ν,where r,ν are nonnegative integers. In this setting, we compute y¯in by(20)y¯in=φtn+cih−τ,tn+cih−τ≤t0,∑ℓ=−rνℒℓηyin−m+ℓ,tn+cih−τ>t0.
Note that we should takeν<m to ensure that no unknown values yin−m+ℓ are involved in the interpolation.


## 3. GDN-Stability

In this section, we investigate the stability property of exponential Runge–Kutta methods.Definition 2.
A numerical method is called GDN-stable, if numerical solutionsyn,zn of the problem Pμ,σ,γ with different initial functions ϕt,ψt satisfy(21)yn−zn≤Cmaxt0−τ≤t≤t0ϕt−ψt,n≥0,where the constant C depends on the problem and the method but not on t nor n.Remark 1.
It is widely recognized that the numerical method should preserve certain property of the exact solution. We say a method is GRN-stable if the contractive property stated in Lemma2 holds. It is clear that GDN-stable is weaker than GRN-stable. On the one hand, it would seem difficult to verify the GRN-stability. On the other hand, we will show that even in the case of GDN-stability, some strong and nontrivial conditions are required.
In order to investigate the GDN-stability of exponential Runge–Kutta methods, we introduce the concept of strong exponential algebraic stability, which is the generalization of strong algebraic stability defined in [22] to the exponential integrators.Definition 3.
The exponential Runge–Kutta method ((17) and (18)) is called strongly exponentially algebraically stable, if for i=0,1,…,s, the matrices(22)Wi=AiA+ATAiT−aiaiT,are positive semidefinite and μaij≥0, j=1,…,s. Here c0=1 and(23)Ai=blockdiagai1Teci−c1hL,ai2Teci−c2hL,…,aisTeci−cshLfori=0,1,…,s;ai=ai1,ai2,…,aisTfori=0,1,…,s;a0j=bj,forj=1,…,s.
Now we point out that strong exponential algebraic stability implies the GDN-stability.Theorem 1.
Apply the exponential Runge–Kutta method to problemPμ,σ,γ. If the method is strongly exponentially algebraically stable, then the method is GDN-stable, i.e.,(24)yn−zn≤e1/2T−t0C0L02maxt0−τ≤t≤t0ϕt−ψt,where C0=max∑i=1sbie1−cihμL,∑j=1saijeci−cjhμL,i=1,2,…,s and L0=supη∈0,1∑j=−rvℓjη.Proof.
Letyn+1,yjn,y¯jn, zn+1,zjn,z¯jn be the numerical solutions of schemes (17) and (18) for problem (4) with differential initial functions ϕt,ψt, respectively. For the sake of simplicity, we use the following abbreviations:(25)un=yn−zn,uin=yin−zin,u¯in=y¯in−z¯in,Qin=hftni,yin,y¯in−ftni,zin,z¯in,Qi1n=hftni,yin,y¯in−ftni,zin,y¯in,Qi2n=hftni,zin,y¯in−ftni,zin,z¯in,Qin=Qi1n+Qi2n.
From the exponential Runge–Kutta schemes (17) and (18), we obtain the following error equations:(26)un+1=ehLun+∑i=1sbiQin,(27)uin=ecihLun+∑j=1saijQjn,i=1,2,…s.
We first derive a recurrence relation forun+1. Substituting (27) into (26) leads to(28)un+12=ehLun+∑i=1sbiQin,ehLun+∑i=1sbiQin=ehLun2+2∑i=1sℜbiQin,ehLun+∑i,j=1sbiQin,bjQjn=ehLun2+2∑i=1sℜbiQin,e1−cihLuin−∑j=1saijQjn+∑i,j=1sbiQin,bjQjn=ehLun2+2∑i=1sℜbiQin,e1−cihLuin−∑i,j=1sQin,biTe1−cihLaij+ajiTbjTe1−cjhLT−biTbjQjn≤ehL2un2+2∑i=1sℜQin,biTe1−cihLuin.
Recalling that problemPμ,σ,γ satisfies conditions (5)–(9), we obtain(29)un+12≤e2hμLun2+2∑i=1sℜQi1n+Qi2n,biTe1−cihLuin≤un2+2h∑i=1sbie1−cihμLσuin2−σuinu¯in≤un2+h∑i=1sbie1−cihμLσuin2−σu¯in2≤un2−hσ∑i=1sbie1−cihμLu¯in2.
In view of approximation (20) for delay term, we find that(30)u¯in2=∑ℓ=−rνℒℓηuin−m+ℓ2≤L02max−r≤ℓ≤νuin−m+ℓ2,with L0=supη∈0,1∑ℓ=−rvℒℓη. Substituting (30) into (29) yields a recurrence relation for the bound of un:(31)un+12≤1−hσCbL0maxun2,maxi,ℓ∈Guin−m+ℓ,where the constant Cb=∑i=1sbie1−cihμL and the index set G=i,ℓ|i=1,2,…,s;ℓ=−r,−r+1,…,ν.
Following the same line as above, using conditions (5)–(9) and bound (30) for u¯in, we can estimate the inner stage difference uin given by (27) as follows:(32)uin2=ecihLun+∑j=1saijQjn,ecihLun+∑j=1saijQjn=ecihLun2+2∑j=1sℜQjn,aijTeci−cjhLujn−∑j,k=1sQjn,aijTeci−cjhLajk+akjTaikTeci−ckhLT−aijTaikQkn≤ecihLun2+2∑j=1sℜQjn,aijTeci−cjhLujn≤un2−hσ∑j=1saijeci−cjhμLu¯jn2≤1−hσCiL02maxun2,maxj,ℓ∈Gujn−m+ℓ2,where Ci=∑j=1saijeci−cjhμL.
SettingC0=−σmaxi=1,…,sCb,CiL02, then we can rewrite bounds (31) and (32) into the same form(33)un+12≤1+hC0maxun2,maxj,ℓ∈Gujn−m+ℓ2,(34)uin2≤1+hC0maxun2,maxj,ℓ∈Gujn−m+ℓ2.
Finally, we can get the GDN-stability by the method of induction. LetΦ denote the difference in the initial function, i.e.,(35)Φ=maxt0−τ≤t≤t0ϕt−ψt.
We intend to show that(36)un2≤1+hC0nΦ2,(37)uin2≤1+hC0n+1Φ2.
These can be easily checked forn=0:(38)u02=y0−z02=yt0−zt02≤Φ2,ui02≤1+hC0maxu02,maxj,ℓ∈Guj−m+ℓ2≤1+hC0Φ2,and for n=1:(39)u12≤1+hC0maxu02,maxj,ℓ∈Guj−m+ℓ2≤1+hC0Φ2,ui12≤1+hC0maxu12,maxj,ℓ∈Guj1−m+ℓ2≤1+hC02Φ2.
Suppose the statement is true for alln≤k. Then, for n=k+1, we have(40)uk+12≤1+hC0maxuk2,maxj,ℓ∈Gujk−m+ℓ2≤1+hC0k+1Φ2,uik+12≤1+hC0maxuk+12,maxj,ℓ∈Gujk+1−m+ℓ2≤1+hC0k+2Φ2.
Therefore, bounds (36) and (37) hold for all n:(41)un2≤1+hC0nΦ2≤enhC0Φ2≤eT−t0C0Φ2.
It follows that(42)un≤e1/2T−t0C0Φ.


## 4. D-Convergence

This section is devoted to the D-convergence of exponential Runge–Kutta methods.Definition 4.
Apply an exponential Runge–Kutta method to problemPμ,σ,γ with yn=ϕtn for n≤0 and yin=ϕtni for n<0. The method is called to be D-convergent of order p, if there exists a constant H, for 0<h≤H, such that the global error satisfies(43)yn−ytn≤Ctnhp,where the function Ctn depends only on the bounds Mi, the delay τ, the parameters μL,σ,γ, and the method.
For convenience, let us denote(44)yn=y1ny2n⋮ysn,y¯n=y¯1ny¯2n⋮y¯sn,Yn=ytn1ytn2⋮ytns,Y¯n=ytn1−τytn2−τ⋮ytns−τ,A=aijs×s,b=b1b2⋮bs,E=ec1hLec2hL⋮ecshL,ftn,yn,y¯n=ftn1,y1n,y¯1nftn2,y2n,y¯2n⋮ftns,ysn,y¯sn,Ftn,Yn,Y¯n=ftn1,ytn1,ytn1−τftn1,ytn2,ytn2−τ⋮ftns,ytns,ytns−τ.
Moreover, we extend the inner product and norm onℂN to space ℂNs:(45)u,v=∑i=1Nui,vi,u=∑i=1sui2,ui,vi∈ℂN.
Using above notations, the exponential Runge–Kutta methods (17) and (18) can be written in the compact form(46)yn+1=ehLyn+hbTftn,yn,y¯n,yn=Εyn+hAftn,yn,y¯n.


### 4.1. Stage Order and Diagonal Stability

The D-convergent order is partly determined by the stage order, which is defined as follows.Definition 5.
The exponential Runge–Kutta methodA,b,c is called p−stage order, if the numerical solutions yn and yn have local truncation error of order p+1. That is, there exists a constant h0, such that for 0<h≤h0,(47)ytn+1=ehLytn+hbTFtn,Yn,Y¯n+CENhp+1,Yn=Eytn+hAFtn,Yn,Y¯n+CEREhp+1,where EJ=1,1,…,1JT for J=NorJ=RE, in which RE is the row number of E. Here constants h0 and C depend on the method and the derivative bounds Mi.
Using the variation-of-constants formula and the Taylor expansion, it is easy to obtain the following.Lemma 3.
An exponential Runge–Kutta method has stage orderp, if its coefficients A,b,c satisfy(48)∑i=1sbi⋅cik−1=φk⋅k−1!,k=1,2,…,p,∑j=1saij⋅cjk−1=φk⋅cikk−1!,i=1,2,…,s,k=1,2,…,p.
To build the D-convergent result, some stability properties of the method are required. Besides the strongly exponential algebraical stability studied in Section3, the following concept of diagonal stability is also needed.Definition 6.
The exponential Runge–Kutta methodA,b,c is said to be diagonally stable, if there exists a positive definite diagonal matrix D such that DA+ATD is positive definite.
Givenσ≤0, let Dσ be the set of Z such that Z=blockdiagZ1,…,Zs,μZi≤σ,Zi∈ℂN×N,i=1,2,…,s, and INs be the identity matrix of order Ns×Ns. Diagonal stability implies that the inversion of INs−hAZ is bounded.Lemma 4.
(see [22]). Assume that the exponential Runge–Kutta method is diagonally stable; then, there exist constants D1,D2, which depends only on the method, such that for any given h>0 and Z∈Dσ, there holds(49)INs−hAZ−1≤D1,haiTZINs−hAZ−1≤D2,i=0,1,2,…,s.


### 4.2. D-Convergence

In this section, we investigate the D-convergence of exponential Runge–Kutta methods. We first combine two stability properties to get the following result.Lemma 5.
If the exponential Runge–Kutta method is strongly exponentially algebraically stable and diagonally stable, then for anyh>0 and Z∈Dσ, we have(50)ecihL+haiTZINs−hAZ−1Ε≤1,i=0,1,2,…,s.Proof.
LetK be the abbreviation for INs−hAZ−1E and it can be rewritten as K=K1T,K2T,…,KsTT with Ki∈ℂN×N. We intend to prove that for any w∈ℂN, there holds(51)ecihL+haiTZKw,ecihL+haiTZKw≤w2.
Indeed, we have(52)ecihL+haiTZKw,ecihL+haiTZKw=ecihLw2+2ℜecihLw,haiTZKw+haiTZKw,haiTZKw=ecihLw2+2∑j=1sℜecihLw,haijZjKjw+∑j=1shaijhLZjKjw,∑j=1shaijZjKjw.
Next, we splitE as(53)E=INs−hAZINs−hAZ−1E=INs−hAZ−1E−hAZINs−hAZ−1E=K−hAZK.
In view of the definition ofE, we have the following identity:(54)ecihL=Ki−h∑j=1saijZjKj,i=1,…,s.
Substituting into (52) yields(55)ecihL+haiTZKw,ecihL+haiTZKw=ecihLw2+2∑j=1sℜeci−cjhLKj−∑l=1sajlZlKlw,haijZjKjw+∑j=1shaijZjKjw,∑j=1shaijZjKjw=ecihLw2+2∑j=1sℜaijTeci−cjhLKjw,hZjKjw−∑j,l=1saijTeci−cjhLajl+aljTailTeci−cjhLT−aijTailhZjKjw,hZlKlw.
Note that the last term in above equation is the component-wise form ofWi defined in (22). By virtue of the Definition 3 of strongly exponentially algebraically stable, this term is nonpositive. The penultimate term is also nonpositive due to the fact that μaij≥0 and hμZj≤hσ≤0s. Hence, we conclude that(56)ecihL+haiTZKw,ecihL+haiTZKw≤ecihLw2≤ecihμL2w2≤w2,where we used the assumption μL≤0. This ends the proof.
Now, we are ready to state our main result on the D-convergence of the exponential Runge–Kutta method.Theorem 2.
Let an exponential Runge–Kutta method of stage orderp is strongly exponentially algebraically stable and diagonally stable. Moreover, if the order of the interpolation for the delay term is q, then the exponential Runge–Kutta method is D-convergent of order minp,q+1.Proof.
For convenience, we denote(57)δn=yn−ytn,δin=yin−ytni,δn=δ1nT,…,δsnTT,δnτ,i=y¯in−ytni−τ.
Since the stage order of the method isp, by Definition 5, we have(58)δn+1=ehLδn+hbTZnδn+ftn,Yn,y¯n−Ftn,Yn,Y¯n+CENhp+1,(59)δn=Eδn+hAZnδn+ftn,Yn,y¯n−Ftn,Yn,Y¯n+CEREhp+1,with Zn=blockdiagZin,…,Zsn and(60)Zin=∫01∂f∂ytni,ytni+θδin,y¯indθ,where ∂f/∂y denotes the partial derivative of function f with respect to the second variable. Noting assumption (5), we find readily that μ∂f/∂y≤σ.
Rewriting (59) as(61)δn=INs−hAZn−1Eδn+hAftn,Yn,y¯n−Ftn,Yn,Y¯n+CEREhp+1,and inserting it into (58) yields(62)δn+1=ehL+hbTZnINs−hAzn−1Εδn+ChbTZnINs−hAZn−1ERΕhp+1+hbTZnINs−hAZn−1hAftn,Yn,y¯n−Ftn,Yn,Y¯n+hbTftn,Yn,y¯n−Ftn,Yn,Y¯n+CENhp+1.
Using the bounds given in Lemmas4 and 5, we have(63)δn+1≤δn+hD2Aftn,Yn,y¯n−Ftn,Yn,Y¯n+Chp+1D2ERE+EN+hbTftn,Yn,y¯n−Ftn,Yn,Y¯n≤δn+CD2s+1shp+1−hσD2∑i=1s∑j=1saij+∑i=1sbimax1≤i≤sy¯in−ytni−τ.
Note that the boundy¯in−ytni−τ in last term stands for the Lagrange interpolation error of the delay term (see (20)). It follows that(64)y¯in−ytni−τ≤∑ℓ=−rνℒℓηδin−m+ℓ+Mq+1q+1!ωηhq+1≤L0max−r≤ℓ≤νδin−m+ℓ+Mq+1q+1!ω0hq+1,with q=r+ν, ωη=∏ℓ=−rνη−ℓ, and ω0=maxη∈0,1ωη.
Therefore, we have from the above that(65)δn+1≤1+hγ0maxi,ℓ∈Gδn,δin−m+ℓ+γ1hq+2+Cγ2hp+1,with(66)γ0=−σL0D2∑i=1s∑j=1saij+∑i=1sbi,γ1=−σMq+1q+1!ω0D2∑i=1s∑j=1saij+∑i=1sbi,γ2=D2s+1s.
Following the same line as in the derivation of (65), for δin, we infer from (58) and (59) and the bound given in Lemmas 4 and 5 and (64) that(67)δin≤1+hγ0imaxi,ℓ∈Gδn,δin−m+ℓ+γ1ihq+2+Cγ2hp+1,with(68)γ0i=−σL0D2∑k=1s∑j=1sakj+∑j=1saij,γ1i=−σMq+1q+1!ω0D2∑k=1s∑j=1sakj+∑j=1saij.
Following the proof of GDN-stability (see Theorem1), setting that γ¯0=maxiγ0,γ0i, γ¯1=maxiγ1,γ1i,γ2, then we have(69)δn+1≤1+hγ¯0maxi,ℓ∈Gδn,δin−m+ℓ+Cγ¯1hminq+2,p+1,δin≤1+hγ¯0maxi,ℓ∈Gδn,δin−m+ℓ+Cγ¯1hminq+2,p+1.
Noting thatδ0=0 and δin=0 for n<0, we conclude that(70)δn≤∑j=0n−11+hγ¯0jCγ¯1hminq+2,p+1,δin≤∑j=0n1+hγ¯0jCγ¯1hminq+2,p+1.
Hence,(71)δn≤1+hγ¯0n−1hγ¯0Cγ¯1hminq+2,p+1≤Ctnhminq+1,p,with Ct=Cγ¯1expt−t01+hγ¯0−1/γ¯0. This completes the proof.


## 4.1. Stage Order and Diagonal Stability

The D-convergent order is partly determined by the stage order, which is defined as follows.Definition 5.
The exponential Runge–Kutta methodA,b,c is called p−stage order, if the numerical solutions yn and yn have local truncation error of order p+1. That is, there exists a constant h0, such that for 0<h≤h0,(47)ytn+1=ehLytn+hbTFtn,Yn,Y¯n+CENhp+1,Yn=Eytn+hAFtn,Yn,Y¯n+CEREhp+1,where EJ=1,1,…,1JT for J=NorJ=RE, in which RE is the row number of E. Here constants h0 and C depend on the method and the derivative bounds Mi.
Using the variation-of-constants formula and the Taylor expansion, it is easy to obtain the following.Lemma 3.
An exponential Runge–Kutta method has stage orderp, if its coefficients A,b,c satisfy(48)∑i=1sbi⋅cik−1=φk⋅k−1!,k=1,2,…,p,∑j=1saij⋅cjk−1=φk⋅cikk−1!,i=1,2,…,s,k=1,2,…,p.
To build the D-convergent result, some stability properties of the method are required. Besides the strongly exponential algebraical stability studied in Section3, the following concept of diagonal stability is also needed.Definition 6.
The exponential Runge–Kutta methodA,b,c is said to be diagonally stable, if there exists a positive definite diagonal matrix D such that DA+ATD is positive definite.
Givenσ≤0, let Dσ be the set of Z such that Z=blockdiagZ1,…,Zs,μZi≤σ,Zi∈ℂN×N,i=1,2,…,s, and INs be the identity matrix of order Ns×Ns. Diagonal stability implies that the inversion of INs−hAZ is bounded.Lemma 4.
(see [22]). Assume that the exponential Runge–Kutta method is diagonally stable; then, there exist constants D1,D2, which depends only on the method, such that for any given h>0 and Z∈Dσ, there holds(49)INs−hAZ−1≤D1,haiTZINs−hAZ−1≤D2,i=0,1,2,…,s.


## 4.2. D-Convergence

In this section, we investigate the D-convergence of exponential Runge–Kutta methods. We first combine two stability properties to get the following result.Lemma 5.
If the exponential Runge–Kutta method is strongly exponentially algebraically stable and diagonally stable, then for anyh>0 and Z∈Dσ, we have(50)ecihL+haiTZINs−hAZ−1Ε≤1,i=0,1,2,…,s.Proof.
LetK be the abbreviation for INs−hAZ−1E and it can be rewritten as K=K1T,K2T,…,KsTT with Ki∈ℂN×N. We intend to prove that for any w∈ℂN, there holds(51)ecihL+haiTZKw,ecihL+haiTZKw≤w2.
Indeed, we have(52)ecihL+haiTZKw,ecihL+haiTZKw=ecihLw2+2ℜecihLw,haiTZKw+haiTZKw,haiTZKw=ecihLw2+2∑j=1sℜecihLw,haijZjKjw+∑j=1shaijhLZjKjw,∑j=1shaijZjKjw.
Next, we splitE as(53)E=INs−hAZINs−hAZ−1E=INs−hAZ−1E−hAZINs−hAZ−1E=K−hAZK.
In view of the definition ofE, we have the following identity:(54)ecihL=Ki−h∑j=1saijZjKj,i=1,…,s.
Substituting into (52) yields(55)ecihL+haiTZKw,ecihL+haiTZKw=ecihLw2+2∑j=1sℜeci−cjhLKj−∑l=1sajlZlKlw,haijZjKjw+∑j=1shaijZjKjw,∑j=1shaijZjKjw=ecihLw2+2∑j=1sℜaijTeci−cjhLKjw,hZjKjw−∑j,l=1saijTeci−cjhLajl+aljTailTeci−cjhLT−aijTailhZjKjw,hZlKlw.
Note that the last term in above equation is the component-wise form ofWi defined in (22). By virtue of the Definition 3 of strongly exponentially algebraically stable, this term is nonpositive. The penultimate term is also nonpositive due to the fact that μaij≥0 and hμZj≤hσ≤0s. Hence, we conclude that(56)ecihL+haiTZKw,ecihL+haiTZKw≤ecihLw2≤ecihμL2w2≤w2,where we used the assumption μL≤0. This ends the proof.
Now, we are ready to state our main result on the D-convergence of the exponential Runge–Kutta method.Theorem 2.
Let an exponential Runge–Kutta method of stage orderp is strongly exponentially algebraically stable and diagonally stable. Moreover, if the order of the interpolation for the delay term is q, then the exponential Runge–Kutta method is D-convergent of order minp,q+1.Proof.
For convenience, we denote(57)δn=yn−ytn,δin=yin−ytni,δn=δ1nT,…,δsnTT,δnτ,i=y¯in−ytni−τ.
Since the stage order of the method isp, by Definition 5, we have(58)δn+1=ehLδn+hbTZnδn+ftn,Yn,y¯n−Ftn,Yn,Y¯n+CENhp+1,(59)δn=Eδn+hAZnδn+ftn,Yn,y¯n−Ftn,Yn,Y¯n+CEREhp+1,with Zn=blockdiagZin,…,Zsn and(60)Zin=∫01∂f∂ytni,ytni+θδin,y¯indθ,where ∂f/∂y denotes the partial derivative of function f with respect to the second variable. Noting assumption (5), we find readily that μ∂f/∂y≤σ.
Rewriting (59) as(61)δn=INs−hAZn−1Eδn+hAftn,Yn,y¯n−Ftn,Yn,Y¯n+CEREhp+1,and inserting it into (58) yields(62)δn+1=ehL+hbTZnINs−hAzn−1Εδn+ChbTZnINs−hAZn−1ERΕhp+1+hbTZnINs−hAZn−1hAftn,Yn,y¯n−Ftn,Yn,Y¯n+hbTftn,Yn,y¯n−Ftn,Yn,Y¯n+CENhp+1.
Using the bounds given in Lemmas4 and 5, we have(63)δn+1≤δn+hD2Aftn,Yn,y¯n−Ftn,Yn,Y¯n+Chp+1D2ERE+EN+hbTftn,Yn,y¯n−Ftn,Yn,Y¯n≤δn+CD2s+1shp+1−hσD2∑i=1s∑j=1saij+∑i=1sbimax1≤i≤sy¯in−ytni−τ.
Note that the boundy¯in−ytni−τ in last term stands for the Lagrange interpolation error of the delay term (see (20)). It follows that(64)y¯in−ytni−τ≤∑ℓ=−rνℒℓηδin−m+ℓ+Mq+1q+1!ωηhq+1≤L0max−r≤ℓ≤νδin−m+ℓ+Mq+1q+1!ω0hq+1,with q=r+ν, ωη=∏ℓ=−rνη−ℓ, and ω0=maxη∈0,1ωη.
Therefore, we have from the above that(65)δn+1≤1+hγ0maxi,ℓ∈Gδn,δin−m+ℓ+γ1hq+2+Cγ2hp+1,with(66)γ0=−σL0D2∑i=1s∑j=1saij+∑i=1sbi,γ1=−σMq+1q+1!ω0D2∑i=1s∑j=1saij+∑i=1sbi,γ2=D2s+1s.
Following the same line as in the derivation of (65), for δin, we infer from (58) and (59) and the bound given in Lemmas 4 and 5 and (64) that(67)δin≤1+hγ0imaxi,ℓ∈Gδn,δin−m+ℓ+γ1ihq+2+Cγ2hp+1,with(68)γ0i=−σL0D2∑k=1s∑j=1sakj+∑j=1saij,γ1i=−σMq+1q+1!ω0D2∑k=1s∑j=1sakj+∑j=1saij.
Following the proof of GDN-stability (see Theorem1), setting that γ¯0=maxiγ0,γ0i, γ¯1=maxiγ1,γ1i,γ2, then we have(69)δn+1≤1+hγ¯0maxi,ℓ∈Gδn,δin−m+ℓ+Cγ¯1hminq+2,p+1,δin≤1+hγ¯0maxi,ℓ∈Gδn,δin−m+ℓ+Cγ¯1hminq+2,p+1.
Noting thatδ0=0 and δin=0 for n<0, we conclude that(70)δn≤∑j=0n−11+hγ¯0jCγ¯1hminq+2,p+1,δin≤∑j=0n1+hγ¯0jCγ¯1hminq+2,p+1.
Hence,(71)δn≤1+hγ¯0n−1hγ¯0Cγ¯1hminq+2,p+1≤Ctnhminq+1,p,with Ct=Cγ¯1expt−t01+hγ¯0−1/γ¯0. This completes the proof.


## 5. Numerical Experiments

In this section, we aim to illustrate the theoretical results obtained in Sections3 and 4. The following two exponential Runge–Kutta methods will be used.One-stage exponential Runge–Kutta method (ERK1):(72)Two-stage exponential Runge–Kutta method (ERK2):(73)with(74)d1z=c1c1−c2φ1c1z,d2z=c2c1−c2φ1c2z,d3z=1c1−c2φ1z,d4z=c12c1−c2φ2c1z,d5z=c22c1−c2φ2c2z,d6z=1c1−c2φ2z.In the above Butcher tableaus,c1,c2 are free parameters in 0,1 and c1≠c2 in ERK2. Note that stage order condition (48) is involved.For the ERK1 method, in order to verify the diagonal stability by the positive definiteness ofDA+ATD, we can simply choose D as the identity matrix I. Noting that the eigenvalues of φ1c1z are nonnegative, we find readily that DA+ATD=c1φ1c1z+c1φ1c1zT is positive for all c1∈0,1. In the numerical tests, we take c1=1. This leads to the exponential Euler method. In this case, the matrices W0 and W1 have the form W0=W1=φ1zTφ1z. Then, the exponential Euler method is also strongly exponentially algebraically stable. With regard to the ERK2 method, it seems difficult for us to theoretically verify the strongly exponentially algebraical stability and diagonal stability. We will show the D-convergence and GDN-stability of the ERK2 method numerically. We remark that in all the tests, the ERK methods are implemented using MATLAB.Example 1.
We consider the following semilinear scalar delay differential equation:(75)y′t=−40yt+ytyt−1−2+gt,t>0,ϕt=cost,t∈−1,0,with gt=42cost−sint−costcost−1. This problem has analytical solution yt=cost.
The ERK1 method withc1=1 and ERK2 method with c1=0.4,c2=0.8 are applied to problem (75). The errors with respect to time step h at final time T=6 are recorded in Table 1. The desired first and second D-convergence orders are observed for ERK1 and ERK2 methods, respectively.
To investigate the GDN-stability, we compare the numerical solutions of (75) with two different initial functions: ϕ1t=cost, ϕ2t=cost+ϵ. Note that the difference of two initial functions is ϵ. The difference of two numerical solutions is plotted in Figure 1 as the function of time. This demonstrates that the difference of numerical solutions is bounded by the difference of initial functions. It further indicates that the constant C in GDN-stability of the method is 1.Table 1
Numerical errors and orders of ERK1c1=1 and ERK2 c1=0.4,c2=0.8 methods for Example 1.
hERK1ERK2ErrorOrderErrorOrder5/2009.4335E-3—3.1599E-06—5/4005.1735E-30.866669.2389E-071.77415/8002.7035E-30.936292.5157E-071.87685/16001.3811E-30.968996.5698E-081.93705/32006.9792E-40.984711.6789E-081.96835/64003.5080E-40.992414.2437E-091.9841Figure 1
GDN-stability of ERK2c1=0.4,c2=0.8,h=0.01 method for Example 1.Example 2.
Next, we consider a realistic model(76)S′t=1−St−mSta+Stxt,x′t=e−τmSt−τa+St−τxt−τ−xt.
This system has been introduced in [23] as a delayed model for describing the standard bacterial growth in a chemostat. Here xt≥0 denotes the concentration of bacteria and St≥0 is the concentration of substrate. In the numerical test, we take constants a=1 and m=2e. It is shown in [24] that if τ>1, the equilibrium 1,0 is asymptotically stable. If τ<1, the equilibrium 1,0 is unstable. In this case, moreover, there is a asymptotically stable equilibrium S¯,x¯ given by S¯=eτ/2e−eτ and x¯=e1−τ−1/e−1/2eτ. 1,0 is called the “washout state” and S¯,x¯ is the “survival state.” (76) is advanced to T=150 by the ERK2 method with c1=0.4, c2=0.8, and h=0.05. The phase portraits corresponding to different initial states are shown in Figure 2 for τ=1.1 and in Figure 3 for τ=0.9. Each phase curve starts with the initial state and gets closer to the predicted equilibrium state which has also been marked. It is clear that the numerical results perfectly fit into the theory given in [24].Figure 2
Phase portraits obtained by the ERK2 method withc1=0.4,c2=0.8,h=0.05 for Example 2 with τ=1.1.Figure 3
Phase portraits obtained by the ERK2 method withc1=0.4,c2=0.8,h=0.05 for Example 2 with τ=0.9.


## 6. Conclusions

In this paper, we investigate the numerical methods for semilinear delay differential equations arising from physics, biology, control theory, etc. A sufficient condition is provided to guarantee the stability of the equations. Exponential Runge–Kutta methods are constructed to solve this type of equations. The sufficient conditions for stability and convergence properties of the exponential Runge–Kutta method are derived and the corresponding methods are given.



---

*Source: 2693940-2022-10-14.xml*