Exponential Multistep Methods for Stiff Delay Differential Equations

: Stiff delay differential equations are frequently utilized in practice, but their numerical simulations are difﬁcult due to the complicated interaction between the stiff and delay terms. At the moment, only a few low-order algorithms offer acceptable convergent and stable features. Exponential integrators are a type of efﬁcient numerical approach for stiff problems that can eliminate the inﬂuence of stiffness on the scheme by precisely dealing with the stiff term. This study is concerned with two exponential multistep methods of Adams type for stiff delay differential equations. For semilinear delay differential equations, applying the linear multistep method directly to the integral form of the equation yields the exponential multistep method. It is shown that the proposed k -step method is stifﬂy convergent of order k . On the other hand, we can follow the strategy of the Rosenbrock method to linearize the equation along the numerical solution in each step. The so-called exponential Rosenbrock multistep method is constructed by applying the exponential multistep method to the transformed form of the semilinear delay differential equation. This method can be easily extended to nonlinear delay differential equations. The main contribution of this study is that we show that the k -step exponential Rosenbrock multistep method is stifﬂy convergent of order k + 1 within the framework of a strongly continuous semigroup on Banach space. As a result, the methods developed in this study may be utilized to solve abstract stiff delay differential equations and can be served as time matching methods for delay partial differential equations. Numerical experiments are presented to demonstrate the theoretical results.


Introduction
Delay differential equations (DDEs) have received substantial attention in recent decades. They are used to model many systems in science and engineering. In these systems, the change of the unknown quantity depends not only on the present state but also on the past state. A typical DDE takes the form of y (t) = F(t, y(t), y(t − τ)), where y(t − τ) stands for the past effects, and τ is the time lag. Such equations arise in various fields such as biology [1], control theory [2], population dynamics [3], and others [4][5][6][7].
To solve delay differential equations numerically, the intuitive approach is to use the mature method for ordinary differential equations (ODEs). Many numerical methods originally designed for ODEs have been applied to DDEs. The reader is referred to the monograph by Bellen and Zennaro [8] for a comprehensive review.
In the last few years, exponential integrators have been getting a lot of attention in the field of numerical solutions of ODEs. Exponential integrators eliminate the influence of stiffness on the scheme by handling the stiff term precisely. The first two exponential integrators of order two and three based on the Adams-Moulton method have been presented in [9]. Along with improving computational efficiency, considerable effort has been directed at developing and analyzing exponential integrators for semilinear stiff ODEs. Various exponential integrators have been investigated, such as exponential Runge-Kutta methods [10,11], exponential multistep methods [12,13], exponential Rosenbrock methods [14][15][16], exponential Taylor methods [17], and exponential general linear methods [18]. A good exposition of exponential integrators can be found in [19].
However, a large number of numerical findings demonstrate that, when a method designed for ODEs is directly applied to the DDEs, the method's accuracy and stability properties cannot be maintained. The majority of current theoretical results on the convergence and stability of the classic methods have the following shortcomings. To preserve the stiff convergent order of the original numerical methods for ODEs, the addition of delay terms necessitates that the method fulfills more stiff convergence order conditions. Moreover, the method's stability conditions get more complicated. More and stronger stability conditions are required to ensure that the scheme meets certain stability requirements. Integration of DDEs requires specifically designed methods other than the stand ODE methods. Furthermore, the stability and convergence properties of numerical methods should be thoroughly investigated.
Analysis of exponential integrators to DDEs has been addressed in the literature. In [20], D-convergence and conditional GDN-the stability of implicit exponential Runge-Kutta methods have been studied for semilinear DDEs. For exponential Runge-Kutta methods, the sufficient conditions of D-convergence and conditional GDN-stability given in [20] imply that the method must be implicit. In [21,22], the stiff convergence and conditional DN-stability of explicit exponential Runge-Kutta methods have been investigated for stiff DDEs. In view of the excellent properties of exponential integrators, the motivation of this work is to develop efficient exponential integrators for DDEs. We aim to construct exponential integrators that have fewer restrictions.
This study is concerned with two exponential multistep methods of Adams type for stiff delay differential equations. For semilinear delay differential equations, applying the linear multistep method directly to the integral form of the equation yields the exponential multistep method. It is shown that the proposed k-step method is stiffly convergent of order k. On the other hand, we can follow the strategy of the Rosenbrock method to linearize the equation along the numerical solution in each step. The so-called exponential Rosenbrock multistep method is constructed by applying the exponential multistep method to the transformed form of the semilinear delay differential equation. This method can be easily extended to nonlinear delay differential equations. The main contribution of this work is that we show that the k-step exponential Rosenbrock multistep method is stiffly convergent of order k + 1 because of the favorable features of the transformed form. The novelty of this study is that the presented methods have fewer restriction conditions than classical linear multistep methods. Moreover, the convergence analysis is carried out within the framework of a strongly continuous semigroup on Banach space. As a result, the methods developed in this study may be utilized to solve abstract stiff delay differential equations and can be served as time matching methods for delay partial differential equations.
The paper is outlined as follows: in Section 2, exponential multistep methods are constructed for semilinear DDEs. It is shown that the proposed k-step method is stiffly convergent of order k. Section 3 is devoted to the exponential Rosenbrock multistep methods. Based on the interpolation with a Hermite node at the beginning of each step, the k-step exponential Rosenbrock multistep method is shown to be stiffly convergent of order k + 1. Finally, numerical experiments are shown to demonstrate the theoretical findings.

Exponential Multistep Methods
In this section, we consider the following semilinear delay differential equation: This form may arise from the space discretization of some semilinear delay parabolic equations [15]. In this case, y ∈ R d is the unknown function, A is a given matrix in R d×d , and g is a nonlinear continuous function. The methods presented in this paper are also applicable to abstract delay differential equations in certain Banach space X with norm · . In this instance, y ∈ X and A is a linear operator (can be unbounded) from X to X.

Construction of Exponential Multistep Methods
We divide the interval [0, T] by time step size h and denote t n = nh; then, for a given integer s, applying the variation-of-constants formula to Equation (1) yields y(t n+s ) = e shA y(t n ) + sh 0 e (sh−θ)A g(t n + θ, y(t n + θ), y(t n + θ − τ)dθ. ( This is the starting point of the exponential multistep methods. Let y n denote the approximation of exact solution y(t) at time t = t n . Similarly, y nτ is the numerical solution at t = t n − τ. Exponential multistep methods are constructed by treating the integral part in (2) numerically.
Here, we set G i = g(t i , y i , y iτ ) for simplicity. According to the Newton backward interpolation formula, the polynomial p n can be expressed as follows: Here, ∇ j G n denotes the jth backward difference, defined recursively by Note that the delay terms y nτ are involved in the interpolation data G n . For t n − τ ≤ 0, we can simply set y nτ by the initial function, i.e., y nτ = φ(t nτ ). When t n − τ > 0, we approximate it by the Lagrange interpolation polynomial q nτ through points (t n−m−v , y n−m−v ), (t n−m−v+1 , y n−m−v+1 ), . . . , (t n−m+r , y n−m+r ).
Here, integer m satisfies τ = mh − δh, 0 ≤ δ < 1. r ≤ m is required to ensure that no unknown values y k are used. The Lagrange interpolation polynomial q nτ can be expressed as follows: We now have all of the pieces necessary to construct the exponential multistep methods. This paper focuses on multistep methods of Adams type; hence, we set s = 1 in (2). The numerical method is obtained by replacing the nonlinearity g in (2) by the interpolation polynomial p n . Therefore, we can write y n+1 = e hA y n + h 1 0 e (1−θ)hA p n (t n + θh)dθ; then, we have the k-step exponential multistep method for semilinear DDEs (1) The weights β j are defined as Remark 1. We remark that the proposed method can be viewed as a small perturbation of the exponential Euler method, which is given by The weights β j given in (5) can be expressed as the linear combination of ϕ-functions, which are defined as follows: Note that ϕ-functions enjoy the following recurrence relation: By introducing the ϕ-functions, the weights of the k-step exponential multistep methods up to k = 7 are given by Remark 2. It is obvious that efficiently computing the ϕ-functions is critical for the implementation of exponential multistep methods. Several efforts have been made on this line; see [23][24][25][26]. For large-scale problems, it is more efficient to compute the product of a matrix exponential function with a vector directly without the explicit form of the matrix exponential.

Remark 3.
For multistep methods, the starting values y 1 , y 2 , . . . , y k−1 are unknown. Normally, they are obtained by one-step methods. Here, we take the method proposed in [12]. The details are listed in Appendix A.

Stiff Convergence of Exponential Multistep Methods
In this subsection, we investigate the stiff convergence of the k-step exponential multistep method (4). Stiff convergence refers to the error bounds being of form Ch p , where the constant C is independent of the stiffness and the step size h. The analysis is carried out under the framework of semigroup. The following assumptions on the semilinear DDEs (1) are suggested. Assumption 1. Let X be a Banach space with norm · . The solution y(t) of Equation (1) is sufficiently smooth. The linear operator A : X → X is the infinitesimal generator of a strongly continuous semigroup e tA on X. The nonlinearity g : [0, T] × X × X → X is Frechet differentiable in a strip along the exact solution y(t). All occurring derivatives are supposed to be uniformly bounded.
Under this assumption, there exists a Lipschitz constant L such that holds in a neighborhood of the exact solution. Moreover, the following parabolic smoothing property holds: e tA + Ae tA ≤ C.
For the ϕ−functions defined by (6), we have It follows that the weights β j (hA) are bounded operators. The main result is the following.
then the error bound holds uniformly for n such that 0 ≤ nh ≤ T. The constant C depends on T, but not on A , n, and h.
Proof. Let e n = y n − y(t n ) denote the error of the method at time t n . It is easy to see that e 0 = 0. The first step is to build the error equation for e n . Note that the exact solution y(t n ) of (1) has the similar form as the numerical solution y n . Indeed, according to the variation-of-constants formula, the solution of (1) has the form with f (t) = g(t, y(t), y(t − τ)) for simplicity. Denote p n as the interpolation polynomial through the exact data (t n−k+1 , f (t n−k+1 )), (t n−k , f (t n−k )), . . . , (t n , f (t n )).
It easily follows that and the interpolation error is Here (10), we obtain Here, the defect δ n+1 is Subtracting (4) from (12) gives the error recursion Solving this recursion with e 0 = 0, we find that We will now show that the error e n indeed satisfies the bound (9). First, for the defect δ n+1 , we infer from (8) and (11) that Next, we have to estimate the term p j (t j + θh) − p j (t j + θh) . Note that functions p j and p j contain the information of the delay terms. To this end, let q nτ denote the interpolation polynomial through the exact history data Then, there exists a polynomial q(δ) of degree r + v + 1 such that Now, we are in the position to estimate the error p j (t j + θh) − p j (t j + θh) .
Combining the above estimation with (13) yields Here, the bound (8) is involved. Using the discrete Gronwall's lemma, we have The assertion is obtained by setting r + v = k − 1.

Exponential Rosenbrock Multistep Methods
In the last section, we designed the exponential multistep methods for semilinear DDEs and investigated their stiff convergence. However, in some instances, the numerical results indicate that this treatment may cause certain issues. Especially when the initial state is far from the equilibrium point, the scheme must take a very tiny time step to make sure it stays stable. This reduces the method's computational efficiency.
In this section, we follow the strategy of the Rosenbrock method [14] to linearize the equation along the numerical solution in each step. The so-called exponential Rosenbrock multistep method is constructed by applying the exponential multistep method to the transformed form of the delay differential equation. This method can be easily extended to nonlinear delay differential equations. Furthermore, we can show that the k-step exponential Rosenbrock multistep method is stiffly convergent of order k + 1 because of the favorable features of the transformed form.

Construction of Exponential Rosenbrock Multistep Methods
To advance the semilinear DDE (1) from t n to t n+1 , we perform the following transformation based on the state (t n , y n , y nτ ). Let d, J, and J τ be the Jacobians of the right-hand side of Equation (1) with respect to the first term t, the second term y, and the third term y τ , respectively, i.e., We evaluate at state (t n , y n , y nτ ) and denote by d n = d(t n , y n , y nτ ), J n = J(t n , y n , y nτ ), J nτ = J τ (t n , y n , y nτ ).
Then, we rewrite the semilinear DDE (1) as Here, the reminder r n is given by It enjoys ∂ ∂t r n (t n , y n , y nτ ) = 0, ∂ ∂y r n (t n , y n , y nτ ) = 0, ∂ ∂y τ r n (t n , y n , y nτ ) = 0.
As will be indicated in Section 3.2, these properties are critical in deriving the convergence order of the exponential Rosenbrock multistep method.

Remark 4.
It appears that the semilinear DDE has been changed into a more complicated form. Nevertheless, the Lipschitz constant of the remainder r n in (14) is smaller than that of the nonlinear term g in (1). The numerical results presented in Section 4 indicate that the exponential multistep method based on (14) is more accurate than that based on the original form (1). Indeed, we will show that the k-step exponential Rosenbrock multistep method is stiffly convergent of order k + 1. The order of this method is higher than that of the exponential multistep method.
We now construct the exponential Rosenbrock multistep method for semilinear DDE (1) based on the form (14).
and evaluate it at t n + θh − τ. That is, For the reminder term r n , we approximate it by the interpolation polynomial p n through points (t n−k+1 , r n (t n−k+1 , y n−k+1 , y n−k+1,τ )), . . . , (t n , r n (t n , y n , y nτ )).
Noting the property (15), for the given data, we have r n (t n ) = 0. Hence, we further require that p n (t n ) = 0. In other words, we employ an incomplete Hermite interpolation to approximate r n . This is the source for the k + 1 order convergent of a k-step exponential Rosenbrock multistep method. Denoting R n i = r n (t i , y i , y iτ ) for simplicity, the Hermite interpolation polynomial p n reads Replacing the delay term y(t n + θh − τ) and the reminder term r n in (16) by polynomials η n and p n , respectively, gives y n+1 = e hJ n y n + hϕ 1 (hJ n )t n d n + h 2 ϕ 2 (hJ n )d n + h k ∑ j=1 β j (hJ n )J nτ ∇ j y nτ Here, ϕ j (z) and β j (z) are defined by (6) and (7). This is the exponential Rosenbrock multistep method for semilinear DDE (1). The procedure for determining the stating values y 1 , y 2 , . . . , y k−1 are sketched in Appendix B.

Remark 5.
The exponential multistep method can be easily extended to a nonlinear delay differential equation The exponential Rosenbrock multistep method is based on the transformed form y (t) = ∂F ∂y y(t) + ∂F ∂t t + ∂F ∂y τ y(t − τ) + r n (t, y(t), y(t − τ)), 0 ≤ t ≤ T, with the corresponding reminder r n . All of the derivatives take values at state (t n , y n , y nτ ).

Stiff Convergence of Exponential Rosenbrock Multistep Methods
This subsection is devoted to investigating the stiff convergence of the exponential Rosenbrock multistep methods for semilinear DDEs (1).
Under Assumption (1), these Jacobians satisfy the Lipschitz condition: Note that the nonlinear operator J can be viewed as a perturbation of g. We infer from Assumption (1) and the perturbation result in [27] that J generates a strongly continuous semigroup e tJ . Hence, there exist constants C, ω, such that e tJ n ≤ Ce ωt , J n e tJ n ≤ C hold uniformly in a neighborhood of the exact solution. It further implies that ϕ j (hJ n ) and β j (hJ n ) are bounded operators. Indeed, we have We are now in a position to state our main result on the convergence of the exponential Rosenbrock multistep methods. In contrast to the case of exponential multistep methods, the Jacobians vary with each step. Hence, the proof is more involved.

Theorem 2.
Apply the k−step exponential Rosenbrock multistep method (17) to the semilinear DDE (1) satisfying Assumption 1. If the Lagrange interpolation (3) satisfies r + v = k and the starting values satisfy then the error bound holds uniformly for n such that 0 ≤ nh ≤ T. The constant C depends on T, but not on A , n and h.
In this case, we have the interpolation error Replacing the nonlinearity term f n (t) and the delay term y(t n + θh − τ) in (19) by p n and η n respectively yields y(t n+1 ) = e hJ n y(t n ) + hϕ 1 (hJ n )t n d n + h 2 ϕ 2 (hJ n )d n This means that the exact solution y(t n+1 ) can be expressed in a similar form as the numerical solution. Here, the defects ρ n+1 and σ n+1 are Thus, we have the error recursion by taking the difference between (17) and (20): Solving this recursion and using e 0 = 0 gives Under Assumption 1, it follows from (18) that defects ρ n+1 and σ n+1 are bounded by with To bound the terms η j (t j + θh − τ) − η j (t j + θh − τ) and p j (t j + θh) − p j (t j + θh) in error Equation (21), we proceed in the same manner as for the term p j (t j + θh) − p j (t j + θh) in Theorem 1. In fact, we have Here, we have taken r + v = k.
Gathering the above bounds with (22), we infer from (18) that e n ≤ C max j=1,2,...,k−1 The result follows now from the application of the discrete Gronwall's lemma. This completes the proof.

Numerical Experiments
In this section, we present some experiments to demonstrate the theoretical results.
Example 1. We consider the following semilinear delay reaction diffusion equation on interval where D is the diffusion coefficient, σ, a, b, and c are positive constants. The added term f 1 (x, t) is specified so that the exact solution is u( The space derivative is approximated by the standard second order central difference method with the mesh size of 0.01. For time discretization, the exponential multistep method (4) and the exponential Rosenbrock multistep method (17) are employed. The relative L 2 errors of the proposed methods are displayed in Figures 1 and 2. We conclude that the k-step exponential multistep method gives a kth order of accuracy, while the k-step exponential Rosenbrock multistep method provides an accuracy of order k + 1. Moreover, the error of the exponential Rosenbrock multistep method is smaller than that of the exponential multistep method when the step size is the same.  Error k=1 k=2 k=3 k=4 slope=2 slope=3 slope=4 slope=5 Figure 2. Errors of k-step exponential Rosenbrock multistep method for problem (23). k = 1, 2, 3, 4.
To demonstrate the superiority of the proposed exponential methods over the classical methods, we measure the computing times of different methods. Adams linear multistep methods (see, e.g., [28]), exponential multistep methods, and exponential Rosenbrock multistep methods are compared in terms of computational times. The CPU times in seconds when all of the methods achieve a common error tolerance 10 −8 are recorded in Table 1. The results reported here show that exponential methods are more efficient than linear multistep methods. Furthermore, the advantages of high order methods over low order methods are proven clearly. Example 2. Next, we consider the following reaction diffusion system with delay We solve this problem on [0, 1] × [−1, 15] with initial conditions The numerical solutions are obtained by the methods used in Example 1. The relative L 2 errors are presented in Figures 3 and 4. Clearly, all the methods achieve their theoretical convergence order. It also confirms that the exponential Rosenbrock multistep methods are more accurate than the exponential multistep methods. Error k=1 k=2 k=3 k=4 slope=1 slope=2 slope=3 slope=4 Figure 3. Errors of k-step exponential multistep method for problem (24). k = 1, 2, 3, 4.  . Errors of k-step exponential Rosenbrock multistep method for problem (24). k = 1, 2, 3, 4.

Conclusions
In this paper, we propose the Adams type exponential multistep methods for semilinear stiff delay differential equations. The exponential multistep methods are obtained by applying the linear multistep method directly to the integral form of the equation, while the exponential Rosenbrock multistep methods are constructed by applying the exponential multistep methods to a transformed form of the equation. We show that the k-step exponential multistep method is stiffly convergent of order k and the stiff convergence order of the k-step exponential Rosenbrock multistep method is k + 1. We demonstrate that the exponential Rosenbrock multistep methods are superior to the exponential multistep methods in terms of accuracy. In summary, for semilinear stiff delay differential equations, this work provided efficient numerical methods with fewer restriction conditions.
Finally, we point out that, while the proposed exponential Rosenbrock multistep methods can be easily applied to nonlinear DDEs, the convergence analysis is restricted to semilinear stiff DDEs. The convergence property of the methods for nonlinear DDEs should be investigated in future research. Moreover, fractional [29] and stochastic [30] differential equations are also widely used in a variety of fields of applied mathematics. The exponential integrators for these problems need to be validated in the future.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
To determine the starting values y 1 , y 2 , . . . , y k−1 for a k-step exponential multistep method, we approximate g in the integral form (2) by the interpolation polynomial p 0 through the points (t 0 , G 0 ), (t 1 , G 1 ), . . . , (t k−1 , G k−1 ). Then, we have Here, j G n denotes the jth forward difference, which defined recursively by 0 G n = G n , j G n = j−1 G n+1 − j−1 G n , j ≥ 1.
As pointed out in [12], under appropriate assumptions, the nonlinear system (A1) has a unique solution. The solution can be obtained by fixed point iteration.