Exact and Nonstandard Finite Difference Schemes for Coupled Linear Delay Differential Systems

: In recent works, exact and nonstandard ﬁnite difference schemes for scalar ﬁrst order linear delay differential equations have been proposed. The aim of the present work is to extend these previous results to systems of coupled delay differential equations X (cid:48) ( t ) = AX ( t ) + BX ( t − τ ) , where X is a vector, and A and B are commuting real matrices, in general not simultaneously diagonalizable. Based on a constructive expression for the exact solution of the vector equation, an exact scheme is obtained, and different nonstandard numerical schemes of increasing order are proposed. Dynamic consistency properties of the new nonstandard schemes are illustrated with numerical examples, and proved for a class of methods.


Introduction
Due to the presence of time lags in the dynamics of most real systems, delay differential equations (DDE) have become basic instruments in the mathematical modelling of a wide range of problems in science and engineering, such as in population biology, physiology, epidemiology, economics, and control problems (see, e.g., [1][2][3][4][5], and references therein), and special methods have been developed to compute numerical solutions for DDE [6]. In the case of differential problems without delay, exact schemes have been defined for different particular problems, and the use of nonstandard finite difference (NSFD) numerical schemes has gained increasing interest in the last years [7][8][9]. The NSFD numerical schemes can be competitive in terms of accuracy while providing dynamically consistent solutions, i.e., they can provide numerical discrete solutions that inherit the structural properties defining the dynamical behaviour of the original continuous equation [10]. The possibility of defining NSFD schemes that reproduce the qualitative behaviour of the continuous solutions has made them specially attractive for population and epidemiology models (e.g., [11][12][13][14][15]), and they have also been proposed for some problems with delay [16][17][18][19][20][21]. However, for DDE models the construction of exact schemes, and consequently of NSFD methods derived from them, has not been much developed.
In [22], a NSFD method was proposed for the scalar first order linear delay problem where α, β ∈ R, τ > 0, and f : [−τ, 0] → R is a continuous function. The method of [22] was exact in the initial time interval 0 ≤ t ≤ τ, and then switched to a NSFD method of second order at most.
More recently [23], an exact scheme for problem (1)- (2) was constructed, valid in the whole domain of definition, and a family of increasing order NSFD schemes was defined. The NSFD methods presented in [22,23] were shown to be consistent with different dynamical properties of the continuous problem (1)- (2). In the present work, we consider the coupled linear delay system satisfying the initial condition where X(t) and F(t) are d-dimensional real vector functions, and A and B are d × d commuting real matrices, in general not simultaneously diagonalizable. The usefulness of nonstandard schemes for scalar linear delay problems and their possible advantages over alternative numerical methods have been discussed in [22,23]. Particularly, the family of schemes proposed in [23] allows the computation of numerical solutions for scalar linear delay problems with the required degree of accuracy and with comparatively low computational complexity. Moreover, the numerical approximations obtained with these nonstandard schemes reproduce dynamical properties of the exact continuous solutions, such as asymptotic stability, positivity, and oscillation behaviour.
The aim, and main contribution, of the present work is to make available, for a wide class of coupled linear delay differential systems, NSFD methods that possess analogous advantages to those in the scalar setting, exhibiting similar properties in terms of accuracy and dynamic consistency. It is to be remarked that for a class of the new NSFD schemes proposed in this work, the F M schemes as defined in Theorem 3, it is rigorously proved that they preserve delay dependent stability. This is a property that usual alternative methods, such as natural Runge-Kutta methods, do not possess, and that is challenging to prove for numerical methods for linear delay systems [6] (p. 356).
There are two main difficulties when dealing with problem (3)-(4), compared with the corresponding scalar problem (1)- (2). Firstly, the obtention of an exact constructive solution that would allow deriving an exact scheme. Secondly, once the new NSFD schemes are defined, the process of proving dynamical properties, which is much more complex than in the scalar case. To overcome these difficulties, the key point is to assume commutativity of the coupled coefficient matrices, a property also considered in other problems involving delay systems [24]. With this assumption, a compact expression for the exact solution of problem (3)-(4), analogous to the scalar case, can be obtained. Also, for commuting matrices, a common Schur basis exists and both matrix coefficients in (3), A and B, can be simultaneously reduced to triangular form, which facilitates analyzing the dynamical properties of the new proposed NSFD schemes. This paper is structured as follows. In the next section, based on a constructive expression for the exact solution of the initial value vector problem (3)-(4), an exact scheme that is valid in the whole domain of definition is obtained. In Section 3, a family of new nonstandard schemes of increasing order of accuracy are proposed. Next, in Section 4, dynamic consistency properties of the new nonstandard schemes are illustrated with numerical examples and proved for a class of methods. In the final section, the results are summarized and discussed.
. For m > 1, it is also immediate to check that X (t) = AX(t) + BX(t − τ), and that the expressions given by (5) for two consecutive intervals agree at the connecting points t = mτ. Thus, X(t) is continuous for t > −τ, with continuous derivative for t > 0, and satisfies (3)-(4).
From the exact solution given in Theorem 1, an exact numerical difference scheme can be obtained, in a similar way as done in [23] for the scalar case, as shown in the next theorem.

Theorem 2.
Let h > 0 such that Nh = τ, for some integer N ≥ 1. Writing t n ≡ nh and X n ≡ X(t n ), for n ≥ −N, the numerical solution given by X n = F(t n ), for −N ≤ n ≤ 0, and, for (m − 1)τ ≤ nh < mτ and m ≥ 1 by defines an exact numerical scheme for problem (3)-(4).
, corresponding to the three terms in expression (5). Then, expanding the binomial terms and rearranging and renaming indices, one has In a similar way, one gets The expressions given in Theorems 1 and 2 are also valid when A = 0, i.e., for the particular case of the pure delay problem If A and B are diagonal, or in the case where they are simultaneously diagonalizable after the usual change of variables, problem (3)-(4) consists of d independent scalar problems, and the expressions given by Theorems 1 and 2 for each component of X(t) coincide with those given in [23] for the corresponding scalar problems. Example 1. Figure 1 presents a numerical example of application of the results of this section, showing the continuous solution given by Theorem 1 (lines) and the exact numerical solution of Theorem 2 with N = 5 (points), for the problem (3)-(4) with parameters τ = 1 and (a) (b)

Nonstandard Finite Difference Methods of Increasing Orders
The exact numerical solution given by Theorem 2 has the drawback of the integral term in (6), as an exact expression could be obtained for certain initial functions F(t), but in general a numerical approximation would be needed. A class of methods could be derived by approximating this integral term, either by using some numerical integration algorithm or by approximating the initial function with some family of functions that allowed the explicit computation of the integral. Instead, as proposed in [23] for the scalar problem, a family of nonstandard methods of increasing orders can be derived by computing the exact solution in the first M intervals and then discarding the integral term, as shown in the next theorem. We define two classes of methods of order M, F M and T M methods, depending on whether the full sum in (6) is kept or a truncated sum is used.
Then, both numerical schemes, F M and T M , have local error O(h M+1 ) and order M.
Proof. Let be any vector norm and a compatible norm for matrices, and consider the scheme T M . Assume that X(t k ) − X k = O(h M+1 ) for k ≤ n, which is the case for nh ≤ Mτ. Then, for m ≥ M + 1 and (m − 1)τ ≤ nh < mτ, using (6), one gets Then, by the induction hypothesis, the first term in (10) is O(h M+1 ) and the second term is bounded by Similar arguments result in the same bounds holding for the scheme F M .

Dynamic Consistency Properties
In this section we analyse the consistency between dynamic properties of the numerical solutions resulting from applying the F M and T M schemes defined in Theorem 3 and the continuous solutions of problem (3)-(4).

Asymptotic Stability
First we will show that the F M schemes defined in Theorem 3 preserve delay-dependent stability, i.e., that they are τ(0)-stable [28].
It is well known that for the trivial solution of (3)-(4) to be asymptotically stable it is necessary and sufficient that all the roots λ i of the characteristic equation where I is the d × d identity matrix, have negative real parts, (λ i ) < 0. This condition, involving a transcendental equation with an infinite number of roots, is difficult to verify in general. However, when A and B commute, there is a common Schur basis for them, and they can be simultaneously reduced to triangular form, with elements in the diagonal corresponding to the eigenvalues of each matrix [29]. Thus, in this case, condition (11) is equivalent to where (α i , β i ) are pairs of eigenvalues of A and B, as they appear in the i diagonal position in the common triangular form. Hence, writing (α, β) for any of these pairs, it follows that if the trivial solution of (3)-(4) is asymptotically stable then implies (λ) < 0. Consider now the difference equations system (8) defining the F M scheme. For any n such that (m − 1)τ ≤ nh = nτ/N < mτ, the integer part of n/N is [n/N] = m − 1. Thus, we can write (8) in the form of a Volterra difference system of convolution type, by setting B j = 0, the d-dimensional zero matrix, when j = kN, and when j = kN, for integer k. Thus, using the Z-transform method, it holds that the system (14) is asymptotically stable if all roots of the characteristic equation Now we have the basis to prove our next theorem.

Theorem 4 (τ(0)-stability). Consider problem (3)-(4) and the F M schemes defined in Theorem 3. If the trivial solution of (3)-(4)
is asymptotically stable then the numerical solutions computed using F M schemes are also asymptotically stable.

Remark 3.
For the class of T M schemes, a general and unconditional result similar to Theorem 4 is not to be expected, as shown by considering the simple case where A = 0, M = 1, and N = 1, so that the T 1 scheme reduces to X n+1 = X n + BhX n−1 .

Delay Independent Stability
Our next theorem shows that the class of T M schemes do preserve absolute or delay independent stability, i.e., that they are P-stable [6] (p. 296). This is also trivially the case for F M schemes, as P-stability is a weaker condition than τ(0)-stability.
Theorem 5 (P-stability). Consider problem (3)-(4) and the T M schemes defined in Theorem 3. If the trivial solution of (3)-(4) is asymptotically stable for all values of τ, then the numerical solutions computed using T M schemes are also asymptotically stable.
The solution of the difference system (9) defining the T M scheme is asymptotically stable if all roots of the characteristic equation are inside the unit disc. A nonzero z is a root of (22) if for a pair (α, β) it holds that Thus, if condition (21) hold and we assume that there is a root with |z| ≥ 1, we would get a contradiction, since, from (23), The stability analysis provided by Theorems 4 and 5 assures that, for a fixed delay, the region of asymptotic stability for (3)-(4) is contained in the region of asymptotic stability of F M schemes, while for T M schemes it can only be assured that the region of asymptotic stability of (3)-(4) for all delays is contained in the corresponding region for the numerical solution. However, T M schemes usually perform much better than can be guaranteed, as shown in the next example.
Example 2. Figure 4 shows the numerical solutions computed with the method T 2 , with N = 5, for the pure delay problem (7) with parameters and two different values of delay, τ = 3 and τ = 3.3. Matrix B has real eigenvalues, λ 1 = −0.37 and λ 2 = −0.5. Hence, in this case the trivial solution of (7) is asymptotically stable if all the eigenvalues β of B satisfy |β|τ < π/2 [6] (p. 289), i.e., for τ < π. As shown in Figure 4, the numerical solutions present the correct behaviour, even for values of τ close to the limit of stability. For both components, the solution approach zero as t increases for τ = 3, inside the region of stability (Figure 4, top), while they diverge for τ = 3.3, outside the region of stability (Figure 4, bottom).

Oscillation and Positivity
Our next theorem shows that F M schemes also preserve the oscillation properties of exact solutions for problem (3)-(4).
We recall that a solution of (3) is said to oscillate if every component of the solution has arbitrary large zeros; otherwise it is called non-oscillatory [33] (Definition 5.0.1). It is known that every solution of the delay differential system (3) oscillates if and only if the characteristic Equation (11) has no real roots [33] (Theorem 5.1.1).
We will use the result of the following lemma, whose proof is similar to that of Theorem 7.1.1 in [33]. Proof of Theorem 6. Using common triangular decompositions of A and B, if every solution of (3)-(4) oscillates then, for any pair of ordered eigenvalues (α, β), Equation (13), or equivalently Equation (19), has no real roots. If we assume that there is a non-oscillatory solution of (8) we get a contradiction, since, from Lemma 1, there would be a positive z satisfying and writing z = exp(λh), we would get Equation (19) with λ a real root.

Remark 4.
For the class of T M schemes, a general result similar to Theorem 6 seems difficult, although particular cases could be dealt with, as shown in our next proposition.

Proposition 1.
If every solution of the pure delay problem (7) oscillates, then the numerical solutions computed using the T 1 scheme also oscillate.
Proof. For the pure delay problem (7), an equivalent condition for every solution to oscillate is that B has no real eigenvalues in the interval [−1/eτ, +∞) [33] (Theorem 5.2.2). The characteristic equation for the system of difference equations (9) defining the T 1 scheme, i.e., Equation (22) with A = 0 and and every solution oscillates if (26) has no positive roots [33] (Theorem 7.1.1). But z is a root of (26) if for an eigenvalue β of B it holds that Thus, if every solution of (7) oscillates, so that any possible real eigenvalue β of B satisfies βτ < −1/e, and we assume that there is a positive root of (27), we get a contradiction. From (27), if z is positive, then β is real and z N (z − 1) = βτ/N < 0. Hence, it follows that z < 1 and But for 0 < z < 1, the maximum value of Nz N (1 − z) is attained when z = N/(N + 1), so that Example 3. Figure 5 shows the numerical solution computed with the method T 2 , with N = 10, for the first component of the pure delay problem (7) with parameters τ = 1 and B and F(t) as in Example 2. In this case, every solution oscillates if all the eigenvalues β of B satisfy |β|τ > 1/e ≈ 0.3679. As shown in Figure 4, the numerical solutions preserve the correct behaviour, even for a value of τ very close to the limit of oscillation.

Positivity
Conditions for the solution of a DDE system to preserve positivity, in the sense that for any component-wise positive initial function F(t) the solution always remains positive, are necessarily very restrictive.
Consider the pure delay problem (7). If B = (b ij ) > 0 element-wise, i.e., b ij > 0, i, j = 1 . . . d, then it is clear from the expression of the exact solution given in (5) that for any component-wise positive initial function F(t) all components of the solution X(t) remain positive for all t > 0. In this case, it is also clear from the expressions of F M and T M schemes given in Theorem 3 that the numerical solutions computed with both methods also remain positive for all t > 0.
If B is only non-negative, i.e., B ≥ 0 element-wise, then the exact as well as the numerical solutions remain non-negative for any non-negative initial function and all t > 0. The condition of all elements of B being non-negative is also necessary to preserve positivity, for if there is an element of B, say b 1r , negative, then it is possible to find an initial function, component-wise positive, for which some component of X(t) becomes negative, already in the first interval 0 < t < τ. To see this, take F(t) with components F r (t) = t 2 and F j (t) = δt 2 , j = r, and choose δ such that Taking into account that, from (5), for 0 < t < τ one gets X(t) = BG(t), where the components of G(t) are G r (t) = h(t) and G j (t) = δh(t), j = r, with h(t) = ((t − τ) 3 − (−τ) 3 )/3, it follows that the first component of X(t) becomes negative, For the general linear problem (3)-(4), if B > 0 and also A > 0 element-wise, then it is also immediate that positivity is preserved both in the exact solution and in the numerical solutions computed using the F M and T M schemes. For B ≥ 0, non-negativity of the solutions, both exact and numerical, is preserved if A is Metzler, i.e., with non-negative off-diagonal elements, as then exp(At) is non-negative for any t > 0.

Conclusions
Despite the growing interest in NSFD methods, including their application to some problems with delay, the scheme presented in Theorem 2 is possibly the first example of an exact scheme for a system of delay differential equations, generalising to systems of linear DDE with commuting matrix coefficients the results presented in [23] for scalar linear DDE problems.
The families of F M and T M schemes defined in Theorem 3 allow the computation of numerical solutions for problem (3)-(4) with high accuracy and low computational costs for extended time intervals, showing good dynamic consistency properties. In particular, F M schemes have been proved to preserve delay-dependent asymptotic stability of the continuous solution, i.e., they are τ(0)-stable difference methods, while T M schemes have been proved to preserve delay-independent asymptotic stability, i.e., they are P-stable methods. Also, F M schemes preserve the oscillation behaviour of the exact solution, which has also been proved for the T 1 scheme when applied to the pure delay problem (7). Both types of scheme also provide numerical solutions that remain positive, or non-negative, when the original problem satisfy conditions assuring the corresponding property.
Several problems and lines of research are open from the results presented in this work. Proving dynamic consistency properties similar to those of F M schemes for some particular T M schemes, either in general or when applied to some type of problems or under certain conditions, could deserve further attention, as T M schemes offer the same accuracy than F M schemes with reduced computational needs. Applying the new schemes to low order systems, e.g., with coefficients being 2 × 2 or 3 × 3 matrices, might allow to express the systems of difference equations defining the schemes in the more usual form of a NSFD method, with derivatives for each component being approximated by the corresponding increments divided by a scalar function ϕ(h) = h + O(h 2 ), as has been done for some examples of systems without delay [34,35]. This could also be the case when considering problems where the matrix coefficients A and B posses some special structure.
Author Contributions: Conceptualization, M.Á.C. and F.R.; methodology and formal analysis, all authors; writing-original draft preparation, F.R.; writing-review and editing, all authors.
Funding: This research was funded by Ministerio de Economía y Competitividad grant number CGL2017-89804-R.