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) thatwithFollowing the proof of GDN-stability (see Theorem 1), setting that , , then we haveNoting that and for , we conclude thatHence,with . This completes the proof.
5. Numerical Experiments
In this section, we aim to illustrate the theoretical results obtained in Sections 3 and 4. The following two exponential Runge–Kutta methods will be used. One-stage exponential Runge–Kutta method (ERK1): Two-stage exponential Runge–Kutta method (ERK2): with
In the above Butcher tableaus, are free parameters in and 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 of , we can simply choose as the identity matrix . Noting that the eigenvalues of are nonnegative, we find readily that is positive for all . In the numerical tests, we take . This leads to the exponential Euler method. In this case, the matrices and have the form . 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:with . This problem has analytical solution .
The ERK1 method with and ERK2 method with are applied to problem (75). The errors with respect to time step at final time 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: , . 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 in GDN-stability of the method is 1.
Example 2. Next, we consider a realistic modelThis system has been introduced in [23] as a delayed model for describing the standard bacterial growth in a chemostat. Here denotes the concentration of bacteria and is the concentration of substrate. In the numerical test, we take constants and . It is shown in [24] that if , the equilibrium is asymptotically stable. If , the equilibrium is unstable. In this case, moreover, there is a asymptotically stable equilibrium given by and . is called the “washout state” and is the “survival state.” (76) is advanced to by the ERK2 method with , , and . The phase portraits corresponding to different initial states are shown in Figure 2 for and in Figure 3 for . 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].
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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The author declares that there are no conflicts of interest.
Acknowledgments
This study was supported in part by the National Natural Science Foundation of China (no. 12001115).