*4.1. Asymptotic Stability*

First we will show that the F*<sup>M</sup>* 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 *λ<sup>i</sup>* of the characteristic equation

$$\det(\lambda I - A - \mathbf{e}^{-\lambda \tau} B) = 0,\tag{11}$$

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

$$\prod\_{i=1}^{d} (\lambda - \alpha\_i - \mathbf{e}^{-\lambda \tau} \beta\_i) = 0,\tag{12}$$

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

$$
\lambda - \mathfrak{a} - \mathfrak{e}^{-\lambda \tau} \mathfrak{E} = 0 \tag{13}
$$

implies "(*λ*) < 0.

Consider now the difference equations system (8) defining the F*<sup>M</sup>* 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,

$$X\_{n+1} = \sum\_{j=0}^{n} B\_j X\_{n-j\prime} \tag{14}$$

by setting *Bj* = 0, the *d*-dimensional zero matrix, when *j* = *kN*, and

$$B\_j = \mathbf{e}^{Ah} \frac{\mathbf{B}^{j/N} h^{j/N}}{(j/N)!} \tag{15}$$

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

$$\det(zI - \mathcal{B}(z)) = 0,\tag{16}$$

satisfy <sup>|</sup>*z*<sup>|</sup> <sup>&</sup>lt; 1 [30] (Theorem 5.21), where *<sup>B</sup>*˜(*z*) is the Z transform of *<sup>B</sup>*. In this case,

$$\mathcal{B}(z) = \sum\_{j=0}^{\infty} B\_j z^{-j} = \mathbf{e}^{Ah} \sum\_{k=0}^{\infty} \frac{B^k h^k}{k!} z^{-kN} = \mathbf{e}^{Ah} \mathbf{e}^{Bh/z^N}. \tag{17}$$

Now we have the basis to prove our next theorem.

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

**Proof.** From the common triangular decompositions of *A* and *B*, it follows that every root of (16) must satisfy, for some pair of ordered eigenvalues (*α*, *β*),

$$z - \mathbf{e}^{ah} \mathbf{e}^{\S h / z^N} = 0 \Longrightarrow \ln(z) - ah - \beta h / z^N = 0. \tag{18}$$

Writing ln(*z*) = *λτ*/*N*, so that *<sup>z</sup>*−*<sup>N</sup>* <sup>=</sup> exp(−*λτ*), one gets from (18)

$$
\lambda \pi / N - ah - \beta h \exp(-\lambda \pi) = 0,\tag{19}
$$

which is equivalent to (13), since *h* = *τ*/*N*. Hence, if the trivial solution of (3)–(4) is asymptotically stable it must hold that "(*λ*) < 0, and therefore |*z*| = exp("(*λτ*/*N*)) < 1.

**Remark 3.** *For the class of* T*<sup>M</sup> 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<sup>1</sup> *scheme reduces to*

$$X\_{n+1} = X\_n + BhX\_{n-1}.\tag{20}$$

*If B has a real eigenvalue β, the trivial solution of the pure delay problem (7) is asymptotically stable if* |*β*| < *π*/2*, while the asymptotic stability of (20) requires the more stringent condition* |*β*| < 1 *[31].*

Delay Independent Stability

Our next theorem shows that the class of T*<sup>M</sup>* 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*<sup>M</sup>* schemes, as *P*-stability is a weaker condition than *τ*(0)-stability.

**Theorem 5** (*P*-stability)**.** *Consider problem (3)–(4) and the* T*<sup>M</sup> 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*<sup>M</sup> schemes are also asymptotically stable.*

**Proof.** Using the common triangular forms of *A* and *B*, and considering a pair of ordered eigenvalues (*α*, *β*), a necessary condition for the trivial solution of (3)–(4) to be delay-independent asymptotically stable is [31,32]

$$
\partial \mathcal{R}(\kappa) + |\beta| < 0. \tag{21}
$$

The solution of the difference system (9) defining the T*<sup>M</sup>* scheme is asymptotically stable if all roots of the characteristic equation

$$\det\left(z^{MN+1}I - \mathbf{e}^{Ah} \sum\_{k=0}^{M} \frac{B^k h^k}{k!} z^{(M-k)N}\right) = 0\tag{22}$$

are inside the unit disc. A nonzero *z* is a root of (22) if for a pair (*α*, *β*) it holds that

$$z - \mathbf{e}^{ah} \sum\_{k=0}^{M} \frac{\beta^k h^k}{k!} z^{-kN} = 0. \tag{23}$$

Thus, if condition (21) hold and we assume that there is a root with |*z*| ≥ 1, we would get a contradiction, since, from (23),

$$|z| \le \mathbf{e}^{\Re(ah)} \sum\_{k=0}^{M} \frac{|\beta|^k h^k}{k!} |z|^{-kN} < \mathbf{e}^{(\Re(a) + |\beta|)h} < 1. \tag{24}$$

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*<sup>M</sup>* schemes, while for T*<sup>M</sup>* 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*<sup>M</sup>* 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* T2*, with N* = 5*, for the pure delay problem (7) with parameters*

$$B = \begin{pmatrix} -0.435 & 0.0325 \\ 0.13 & -0.435 \end{pmatrix}, \quad F(t) = \begin{pmatrix} \cos(\pi t) \\ (t+1)^2 \end{pmatrix}.$$

*and two different values of delay, τ* = 3 *and τ* = 3.3*.*

*Matrix B has real eigenvalues, λ*<sup>1</sup> = −0.37 *and λ*<sup>2</sup> = −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).*

**Figure 4.** Numerical solutions computed with the method T<sup>2</sup> for Example 2 with two different values of delay, showing stable and unstable behaviours. (**a**,**b**) First and second component, respectively, with delay *τ* = 3. (**c**,**d**) First and second component, respectively, with delay *τ* = 3.3.

## *4.2. Oscillation and Positivity*

Our next theorem shows that F*<sup>M</sup>* 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).

**Theorem 6** (Oscillation)**.** *If every solution of (3)–(4) oscillates, then the numerical solutions computed using* F*<sup>M</sup> schemes also oscillate.*

We will use the result of the following lemma, whose proof is similar to that of Theorem 7.1.1 in [33].

**Lemma 1.** *Consider the linear system of difference Equation (8) defining the* F*<sup>M</sup> scheme. Every solution of (8) oscillates if and only if the characteristic Equation (16) has no positive roots.*

**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

$$z - \mathbf{e}^{ah} \mathbf{e}^{\otimes h/z^N} = 0,\tag{25}$$

and writing *z* = exp(*λh*), we would get Equation (19) with *λ* a real root.

**Remark 4.** *For the class of* T*<sup>M</sup> 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<sup>1</sup> *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<sup>1</sup> scheme, i.e., Equation (22) with *A* = 0 and *M* = 1, reads

$$\det\left( (z^{N+1} - z^N)I - Bh \right) = 0,\tag{26}$$

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

$$z^{N+1} - z^N = \beta \hbar. \tag{27}$$

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 *<sup>β</sup>* is real and *<sup>z</sup>N*(*<sup>z</sup>* <sup>−</sup> <sup>1</sup>) = *βτ*/*<sup>N</sup>* <sup>&</sup>lt; 0. Hence, it follows that *<sup>z</sup>* <sup>&</sup>lt; 1 and

$$Nz^N(1-z) = -\beta \tau > \mathbf{e}^{-1}.$$

But for 0 <sup>&</sup>lt; *<sup>z</sup>* <sup>&</sup>lt; 1, the maximum value of *NzN*(<sup>1</sup> <sup>−</sup> *<sup>z</sup>*) is attained when *<sup>z</sup>* <sup>=</sup> *<sup>N</sup>*/(*<sup>N</sup>* <sup>+</sup> <sup>1</sup>), so that

$$Nz^N(1-z) \le \left(\frac{N}{N+1}\right)^{N+1} < \mathbf{e}^{-1}.$$

**Example 3.** *Figure 5 shows the numerical solution computed with the method* T2*, 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.*

**Figure 5.** Numerical solution for the first component of Example 1 and zoom-views in different intervals. (**a**) *t* ∈ [0, 50]. (**b**) *t* ∈ [13, 22]. (**c**) *t* ∈ [20, 25]. (**d**) *t* ∈ [46, 50].

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* = (*bij*) > 0 element-wise, i.e., *bij* > 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*<sup>M</sup>* and T*<sup>M</sup>* 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*1*r*, 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 *Fr*(*t*) = *t* <sup>2</sup> and *Fj*(*t*) = *δt* 2, *<sup>j</sup>* <sup>=</sup> *<sup>r</sup>*, and choose *<sup>δ</sup>* such that

$$0 < \delta < -b\_{1r} / |\sum\_{j \neq r} b\_{1j}|.$$

Taking into account that, from (5), for 0 < *t* < *τ* one gets *X*(*t*) = *BG*(*t*), where the components of *<sup>G</sup>*(*t*) are *Gr*(*t*) = *<sup>h</sup>*(*t*) and *Gj*(*t*) = *<sup>δ</sup>h*(*t*), *<sup>j</sup>* <sup>=</sup> *<sup>r</sup>*, with *<sup>h</sup>*(*t*) = ((*<sup>t</sup>* <sup>−</sup> *<sup>τ</sup>*)<sup>3</sup> <sup>−</sup> (−*τ*)3)/3, it follows that the first component of *X*(*t*) becomes negative,

$$X\_1(t) = b\_{1r}h(t) + \delta \sum\_{j \neq r} b\_{1j}h(t) \\ < \left(b\_{1r} - b\_{1r} \frac{\sum\_{j \neq r} b\_{1j}}{|\sum\_{j \neq r} b\_{1j}|}\right) h(t) \\ < 0,$$

since *h*(*t*) > 0 for *t* ∈ (0, *τ*).

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

## **5. 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*<sup>M</sup>* and T*<sup>M</sup>* 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*<sup>M</sup>* schemes have been proved to preserve delay-dependent asymptotic stability of the continuous solution, i.e., they are *τ*(0)-stable difference methods, while T*<sup>M</sup>* schemes have been proved to preserve delay-independent asymptotic stability, i.e., they are *P*-stable methods. Also, F*<sup>M</sup>* schemes preserve the oscillation behaviour of the exact solution, which has also been proved for the T<sup>1</sup> 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*<sup>M</sup>* schemes for some particular T*<sup>M</sup>* schemes, either in general or when applied to some type of problems or under certain conditions, could deserve further attention, as T*<sup>M</sup>* schemes offer the same accuracy than F*<sup>M</sup>* 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.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
