**1. 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–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–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–15]), and they have also been proposed for some problems with delay [16–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

$$\mathbf{x}'(t) \quad = \quad \mathbf{a}\mathbf{x}(t) + \beta \mathbf{x}(t-\tau), \qquad t > 0,\tag{1}$$

$$\mathbf{x}(t) \quad = \quad f(t), \qquad \qquad -\pi \le t \le 0,\tag{2}$$

where *<sup>α</sup>*, *<sup>β</sup>* <sup>∈</sup> <sup>R</sup>, *<sup>τ</sup>* <sup>&</sup>gt; 0, and *<sup>f</sup>* : [−*τ*, 0] <sup>→</sup> <sup>R</sup> 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

$$X'(t) = AX(t) + BX(t - \pi), \qquad t > 0,\tag{3}$$

satisfying the initial condition

$$X(t) = F(t), \qquad -\tau \le t \le 0,\tag{4}$$

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*<sup>M</sup>* 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.

#### **2. Exact Numerical Scheme**

In our next theorem we present an explicit expression for the solution of problem (3)–(4), derived by using the method of steps [25] (pp. 45–47) and an integral convolution [26] (p. 67), in a similar way as was done in [27] for the scalar problem (1)–(2).

*Mathematics* **2019**, *7*, 1038

**Theorem 1.** *The exact solution of (3)–(4) is given by X*(*t*) = *F*(*t*)*, for* −*τ* ≤ *t* ≤ 0*, and, for* (*m* − 1)*τ* < *t* ≤ *mτ and m* ≥ 1*,*

$$X(t) = \sum\_{k=0}^{m-1} \frac{B^k (t - k\tau)^k}{k!} \mathbf{e}^{A(t - k\tau)} F(0) + \sum\_{k=0}^{m-2} \frac{B^{k+1}}{k!} \int\_{-\tau}^0 (t - (k+1)\tau - s)^k \mathbf{e}^{A(t - (k+1)\tau - s)} F(s) ds$$

$$+ \frac{B^m}{(m-1)!} \int\_{-\tau}^{t - m\tau} (t - m\tau - s)^{m-1} \mathbf{e}^{A(t - m\tau - s)} F(s) ds, \quad \text{(5)}$$

*where the second summation is assumed to be empty for m* = 1*.*

**Proof.** For *m* = 1, one has *X*(*t*) = e*AtF*(0) + *B t*−*τ* <sup>−</sup>*<sup>τ</sup>* <sup>e</sup>*A*(*t*−*τ*−*s*)*F*(*s*)*ds*, so that *<sup>X</sup>*(0) = *<sup>F</sup>*(0) and *<sup>X</sup>* (*t*) = *AX*(*t*) + *BF*(*t* − *τ*) = *AX*(*t*) + *BX*(*t* − *τ*). 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 tn* ≡ *nh and Xn* ≡ *X*(*tn*)*, for n* ≥ −*N, the numerical solution given by Xn* = *F*(*tn*)*, for* −*N* ≤ *n* ≤ 0*, and, for* (*m* − 1)*τ* ≤ *nh* < *mτ and m* ≥ 1 *by*

$$X\_{n+1} = \mathbf{e}^{Ah} \sum\_{k=0}^{m-1} \frac{B^k h^k}{k!} X\_{n-kN} + \frac{B^m}{(m-1)!} \int\_{t\_n - m\tau}^{t\_n - m\tau + h} (t\_n - m\tau + h - s)^{m-1} \mathbf{e}^{A(t\_n - m\tau + h - s)} F(s) ds, \tag{6}$$

*defines an exact numerical scheme for problem (3)–(4).*

**Proof.** Write *X*(*t*) = *E*1(*t*) + *E*2(*t*) + *E*3(*t*), corresponding to the three terms in expression (5). Then, expanding the binomial terms and rearranging and renaming indices, one has

$$\begin{split} E\_1(t\_{n+1}) &= \quad E\_1(t\_n + h) = \sum\_{k=0}^{m-1} \frac{B^k(t\_n - k\tau + h)^k}{k!} \mathbf{e}^{A(t\_n - k\tau + h)} F(0) \\ &= \quad \mathbf{e}^{At} \sum\_{k=0}^{m-1} \sum\_{r=0}^k \frac{B^r h^r}{r!} \frac{B^{k-r}(t\_n - r\tau - (k-r)\tau)^{k-r}}{(k-r)!} \mathbf{e}^{A(t\_n - r\tau - (k-r)\tau)} F(0) \\ &= \quad \mathbf{e}^{At} \sum\_{k=0}^{m-1} \frac{B^k h^k}{k!} \sum\_{r=0}^{m-1-k} \frac{B^r (t\_n - k\tau - r\tau)^r}{r!} \mathbf{e}^{A(t\_n - k\tau - r\tau)} F(0) = \mathbf{e}^{At} \sum\_{k=0}^{m-1} \frac{B^k h^k}{k!} E\_1(t\_n - kN). \end{split}$$

In a similar way, one gets

$$E\_2(t\_{n+1}) = \mathbf{e}^{Al} \sum\_{k=0}^{m-2} \frac{B^k h^k}{k!} E\_2(t\_n - kN) = \mathbf{e}^{Al} \sum\_{k=0}^{m-1} \frac{B^k h^k}{k!} E\_2(t\_n - kN),$$

since *E*2(*tn* − (*m* − 1)*N*) = 0 for (*m* − 1)*τ* ≤ *tn* < *mτ*. Also,

$$E\_3(t\_{n+1}) = \mathbf{e}^{t\_{\mathrm{fl}}} \sum\_{k=0}^{m-1} \frac{B^k h^k}{k!} E\_3(t\_n - kN) + \frac{B^m}{(m-1)!} \int\_{t\_n - m\tau}^{t\_n - m\tau + h} (t\_n - m\tau + h - s)^{m-1} \mathbf{e}^{A(t\_n - m\tau + h - s)} F(s) ds, \quad t = 0, 1, \ldots, m$$

so that expression (6) is recovered.

**Remark 1.** *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*

$$X'(t) = BX(t - \tau), \ t > 0, \qquad X(t) = F(t), \ -\tau \le t \le 0. \tag{7}$$

*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 = \begin{pmatrix} -3/2 & 1\\ -2 & 3/2 \end{pmatrix}, \quad B = \begin{pmatrix} 5/4 & -1\\ 2 & -7/4 \end{pmatrix}, \quad F(t) = \begin{pmatrix} 2(t+1) \\ (t+1)^2 \end{pmatrix}.$$

**Figure 1.** Exact solutions (lines) and numerical solutions provided by the exact scheme (points) for the two components of Example 1. (**a**) First component, *X*1(*t*). (**b**) Second component, *X*2(*t*).
