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.

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:where 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 , , and 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]:where stands for the state, is a nonlinear function, 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:where 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 Section 2, 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 , together with a Lagrangian interpolation of order , are D-convergent of order . 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:where constant is the delay term, is a complex-valued matrix, and and are continuous functions.

Let be the usual inner product on . The induced norm and logarithmic norm are denoted by and , respectively. In this paper, we assume that there are constants such thatand the logarithmic norm of matrix , which is defined bysatisfies

In the following, problem (4) satisfying conditions (5)–(9) is denoted by . We further assume that the exact solution is sufficiently smooth in the means thatholds uniformly for .

2.1. Stability of Exact Solution

We begin with the stability analysis of problem . The following stability result for ODEs will be used in the subsequent sections.

Lemma 1. (see [21]). Consider the following linear initial value problem:where . If for , then the exact solution fulfillsNow we state the contractive property of exact solution of problem .

Lemma 2. (see [21]). Let denote the exact solutions of the problem with different initial functions , respectively. Then, for the difference of and , we have

2.2. Numerical Scheme

Noting that problem (4) is semilinear, we can express its exact solution by applying the variation-of-constants formula

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 part in the integral yields the numerical methods. Specifically, a stage exponential Runge–Kutta method with coefficients can be given by the usual Butcher tableau form. Here are linear combinations of functions, which are defined as follows.

Definition 1. For any complex numbers , let . For all integers , define asNote that we have the following recurrence relation:Applying stage exponential Runge–Kutta methods with step size to problem yields the schemewhere and , , and are approximations of , , and , respectively. For convenience, and are abbreviated as and .
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 term can be given by the initial function if . Otherwise, we approximate by the interpolation based on previous inner stage values. Let integer and be such that . Recall that the Lagrange basis functions for uniform mesh have the formwhere are nonnegative integers. In this setting, we compute byNote that we should take to ensure that no unknown values 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 solutions of the problem with different initial functions satisfywhere the constant depends on the problem and the method but not on nor .

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 Lemma 2 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 , the matricesare positive semidefinite and , . Here andNow we point out that strong exponential algebraic stability implies the GDN-stability.

Theorem 1. Apply the exponential Runge–Kutta method to problem . If the method is strongly exponentially algebraically stable, then the method is GDN-stable, i.e.,where and .

Proof. Let , be the numerical solutions of schemes (17) and (18) for problem (4) with differential initial functions , respectively. For the sake of simplicity, we use the following abbreviations:From the exponential Runge–Kutta schemes (17) and (18), we obtain the following error equations:We first derive a recurrence relation for . Substituting (27) into (26) leads toRecalling that problem satisfies conditions (5)–(9), we obtainIn view of approximation (20) for delay term, we find thatwith . Substituting (30) into (29) yields a recurrence relation for the bound of :where the constant and the index set .
Following the same line as above, using conditions (5)–(9) and bound (30) for , we can estimate the inner stage difference given by (27) as follows:where .
Setting , then we can rewrite bounds (31) and (32) into the same formFinally, we can get the GDN-stability by the method of induction. Let denote the difference in the initial function, i.e.,We intend to show thatThese can be easily checked for :and for :Suppose the statement is true for all . Then, for , we haveTherefore, bounds (36) and (37) hold for all :It follows that

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 problem with for and for . The method is called to be D-convergent of order , if there exists a constant , for , such that the global error satisfieswhere the function depends only on the bounds , the delay , the parameters , and the method.
For convenience, let us denoteMoreover, we extend the inner product and norm on to space :Using above notations, the exponential Runge–Kutta methods (17) and (18) can be written in the compact form

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 method is called stage order, if the numerical solutions and have local truncation error of order . That is, there exists a constant , such that for ,where for , in which is the row number of . Here constants and depend on the method and the derivative bounds .
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 order , if its coefficients satisfyTo build the D-convergent result, some stability properties of the method are required. Besides the strongly exponential algebraical stability studied in Section 3, the following concept of diagonal stability is also needed.

Definition 6. The exponential Runge–Kutta method is said to be diagonally stable, if there exists a positive definite diagonal matrix such that is positive definite.
Given , let be the set of such that , and be the identity matrix of order . Diagonal stability implies that the inversion of is bounded.

Lemma 4. (see [22]). Assume that the exponential Runge–Kutta method is diagonally stable; then, there exist constants , which depends only on the method, such that for any given and , there holds

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 any and , we have

Proof. Let be the abbreviation for and it can be rewritten as with . We intend to prove that for any , there holdsIndeed, we haveNext, we split asIn view of the definition of , we have the following identity:Substituting into (52) yieldsNote that the last term in above equation is the component-wise form of 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 and . Hence, we conclude thatwhere we used the assumption . 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 order is strongly exponentially algebraically stable and diagonally stable. Moreover, if the order of the interpolation for the delay term is , then the exponential Runge–Kutta method is D-convergent of order .

Proof. For convenience, we denoteSince the stage order of the method is , by Definition 5, we havewith andwhere denotes the partial derivative of function with respect to the second variable. Noting assumption (5), we find readily that .
Rewriting (59) asand inserting it into (58) yieldsUsing the bounds given in Lemmas 4 and 5, we haveNote that the bound in last term stands for the Lagrange interpolation error of the delay term (see (20)). It follows thatwith , , and .
Therefore, we have from the above thatwithFollowing the same line as in the derivation of (65), for , we infer from (58) and (59) and the bound given in Lemmas 4 and 5 and (64) thatwith