Exact numerical solutions and high order nonstandard difference schemes for a second order delay differential equation

In this work, we obtain an expression, given in terms of some special functions, for the exact numerical solution of the second order delay differential equation x′′(t)=ax(t)+bx(t−τ),$$ {x}^{\prime \prime }(t)= ax(t)+ bx\left(t-\tau \right), $$ with general initial condition x(t)=f(t)$$ x(t)=f(t) $$ for −τ≤t≤0$$ -\tau \le t\le 0 $$ . Based on this exact solution, we propose a family of nonstandard finite difference methods that combine high order accuracy and efficient computational properties. We also show that the numerical solutions obtained with the new schemes are consistent with asymptotic stability properties of the exact solutions.


INTRODUCTION
Exact numerical schemes for delay models, and nonstandard finite difference (NSFD) methods derived from them, have been previously proposed for some first order delay differential equations (DDE) and systems.Garba et al. [1] considered the first order scalar initial value DDE problem x ′ (t) = ax(t) + bx(t − ), t > 0, (1) proposing a numerical method that was exact in the first interval 0 ≤ t ≤ , and then switched to a NSFD scheme of second order at best.García et al. [2] constructed an exact numerical solution for problem (1)-( 2) valid in the whole domain of definition, and proposed a derived family of increasing order NSFD methods.These authors generalized their work in García et al. [3] to coupled linear delay systems, with X(t) and F(t) being d-dimensional real vector functions, and A and B being d × d real matrices, in general not simultaneously diagonalisable, but under the condition that they commuted.It was also shown in previous works [2,3] that the proposed NSFD methods were consistent with dynamic properties of the continuous problems.Castro et al. [4] considered problem (3)-( 4) in the general case of non-commuting matrix coefficients A and B. They provided an expression for the exact numerical solution of (3)-(4) involving infinite series, which were truncated to construct NSFD methods of any desire order.Although it was suggested by numerical examples that the proposed methods preserved the delay-dependent stability of the corresponding continuous problems, that is, that they seemed to be (0)-stable [5], the proof of this property was not attempted there.
It was suggested in Castro et al. [4] that the expressions given there could be used to obtain closed form exact numerical solutions for particular problems, depending on the structure of the specific coefficient matrices, and in particular that they could provide the basis to tackle higher order scalar problems by reducing them to first order vector problems, as the resulting coefficient matrices would be in general non-commuting.In the present work, we show that this strategy can be successfully applied to the second order initial-value DDE problem x ′′ (t) = ax(t) + bx(t − ), t > 0, ( 5) where a, b ∈ R,  > 0, and  ∶ [−, 0] → R is a continuous function.Different types of representations for the solution of DDE linear systems have been proposed, including solutions by definite integrals and series expansions [6], convolutions with particular solutions or fundamental systems obtained through the method of steps or by Laplace transform [3,4,[7][8][9], or infinite series involving matrix Lambert function [10].Representations of the exact solution of problem ( 5)- (6) in the case of pure delay, that is, with a = 0, have been given in terms of delayed sine and cosine functions [11], and of delayed exponentials [12].An expression for the exact solution of ( 5)- (6) was presented in Rodríguez et al. [13], as a preliminary problem resulting from applying the method of separation of variables to an unidimensional wave equation with delay.In that context, the condition a < 0 was assumed, writing a = − 2 , but, as shown in Jornet [14], where it was applied to solve a more general multidimensional wave equation, it is straightforward to adapt the expression given in Rodríguez et al. [13] to the case a > 0. Notwithstanding, these expressions do not seem to be suitable for constructing an exact numerical scheme, and so providing the basis for an efficient computational method.
The structure of this paper is as follows.In the next section, we transform problem ( 5)-( 6) into a vector problem of the form (3)-( 4), and use the results in Castro et al. [4] as the basis to derive an exact numerical solution in the form of a perturbed difference system.Next, in Section 3, a family of NSFD methods, of any desired order of accuracy, are proposed, which are proved to be dynamically consistent with the asymptotic delay-dependent stability of the original continuous problem in Section 4. In the final section, the results are summarized and discussed.

EXACT NUMERICAL SOLUTIONS
In the next Lemma, we recall a result from Theorem 1 in Castro et al. [4] providing an expression for the exact solution of the general vector problem (5)- (6).

, and assume A and I
and define the matrix constants K m r,p , ∀m ≥ 1, by Let X(t) = F(t) for t ∈ [−, 0] and for m ≥ 1 and (m − 1) ≤ t < t + h ≤ m.Then, X(t) is a well-defined function satisfying (3) and (4).
We also recall that, from Lemma 2 in Castro et al. [4], the infinite sums in ( 9) are absolutely convergent, and that the matrices defined in (7) satisfy From Lemma 1, considering a regular mesh of amplitude h = ∕N, for some integer N ≥ 1, an exact numerical scheme can be immediately derived [4,Theorem 2].Writing t n ≡ nh and X n ≡ X(t n ), for n ≥ −N, an exact numerical solution is given by X n = F(t n ), for −N ≤ n ≤ 0, and, for (m − 1) ≤ nh < m and m ≥ 1, by The expressions given in ( 9) and ( 11) are not quite satisfactory, as they include infinite sums that, in general, cannot be reduced to finite expressions.However, for some particular problems, it can provide the basis to obtain finitely expressed explicit solutions, as will be shown next for problem (5)- (6).
To this end, write the second order scalar problem ( 5)- (6) in the form of the first order vector problem (3)-(4), by writing Now, using the relations where 0 as a matrix denotes the null matrix, explicit expressions for the constant matrices K m r,p defined in (8) can be obtained, as given in our next Lemma.(12), corresponding to vector form for the scalar second order delay equation (5), the constant matrices K m r,p defined in (8), with r ≥ p, can be expressed as follows, depending on being p = 0,

Lemma 2. For equation (3) with coefficients A and B as given in
Proof.The expressions for K m r,0 given in ( 14) are immediate, since K m r,0 = A r and A 2 = aI.The rest of the expressions can be proved by induction on the different indices involved.
Thus, for m = 1, since K 1 r,1 = K 1 r−1,0 B, using the expressions in (14) for r − 1 one has as given in (16).Now consider m ≥ 2 and p = 1.From ( 8) and ( 13) it is immediate that = aB + bA, as given in (15).Assume now that the expressions in (15) are valid for p = 1 and up to a certain k ≥ 1.Then, for k + 1 one has, since and hence, (15).
We assume now, for m ≥ 2, that the expressions in (15) are valid up to a certain p, with 1 ≤ p < m − 1, and prove that they also hold for p + 1.
Since K m r,p+1 = BK m r−1,p + AK m r−1,p+1 , taking into account that BK m r−1,p = 0 for r ≤ 2p − 1, and K m r,p+1 = 0 for r ≤ p, it follows that K m r,p+1 = 0 for r ≤ 2p.Then, for r = 2p + 1 one has and hence, (15).Thus, assuming that the expressions in (15) are valid for p + 1 up to a certain k ≥ p + 1, one gets and hence, as given in (15).Finally, consider the case p = m, with m ≥ 2. For r ≤ 2(m − 1) one has K m+1 r,m = 0, so that and for r = 2k one gets which are the expressions for m + 1 in (16).□ Using the expressions in Lemma 2 for the constant matrices K m r,p , we can obtain closed forms for the infinite series in (11).Although this task can be carried out for general values of the coefficients in (5), in what follows, we will consider the case a < 0, which is a necessary condition for the solution of ( 5) to be asymptotically stable [15].To make this case explicit, and simplify some results, we will consider equation ( 5) with coefficients a = − 2 and b = , with ,  ∈ R.
We note that, although the matrix constants K m r,p given in Lemma 1 might in general depend on m, the expressions obtained in Lemma 2 for our particular problem do not depend on the super-index m, so that it will be dropped henceforth.
r! K r,p , integrating by parts and noting that F(t n−mN ) = X n−mN , we write (11) in the form where , and the summation is assumed to be empty for m = 1.In the next Lemma, we give expressions for the matrices G p (h) in terms of Bessel and hypergeometric functions (see, e.g., previous works [16][17][18], for general properties of these special functions).

Lemma 3. The matrices G p (h) in (17) can be expressed in the form
for 1 ≤ p < m, and where J  (h) are Bessel functions of the first kind and 1 F 2 ) is a generalized hypergeometric function.
Proof.Since the series defining G p (h) are absolutely convergent, they can be rearranged by summing separately the terms multiplying different matrices, so that, for 1 , where, taking into account that K r,p = 0 for r ≤ 2(p − 1), and Now, using the duplication formula for the Gamma function [17, p. 138], so that and, similarly, and taking into account the series expansions of the Bessel function of the first kind [17, p. 217], , one gets, with some rearrangements and simplifications, and and taking into account the definition of generalized hypergeometric functions [17, p. 404], so that and, using again the duplication formula for the Gamma function, for (2m)! and (2k)!, ) .

□
From Lemma 3, and the structure of the matrices I, A, B, and AB, it is immediate to write down the elements of the matrices G p (h), as given in the next corollary.
Remark 1.The coefficient of X n in ( 17), e Ah , can be easily computed, since and we note that this expression is recovered taking p = 0 in (18), h cos(h), and h sin(h) [17, p. 228].Hence, the term corresponding to X n can be grouped with the summation in (17) by writing ∑ m−1 p=0 G p (h)X n−pN .Next, we will show that the coefficients of X n−mN in (17) cancel each other, and so this term can be dropped.

Lemma 4. The matrices G m (h) and Q
) , it follows that ) , and we only need to prove that m g i2 (h) =  − 2 m q i2 (h), for i = 1, 2. In fact, as will be shown next, it holds that m g ′ 22 (h) = m g 12 (h) and m q ′ 22 (h) = m q 12 (h), so that checking the previous relation between m g 22 (h) and m q 22 (h) will suffice to prove (20).
From the expressions for the elements of G m (h) given in Corollary 1, using a differentiation formula for the generalized hypergeometric function [17, p. 405, 16.3.4],and the relation between functions J  and 0 F 1 [17, p. 228, 10.16.9], one gets, writing z = −(h∕2) 2 , The corresponding relation for matrices Q m (h) can be easily checked for m = 1, since, from their definition in (7), and it is immediate that 1 q ′ 22 (h) = 1 q 12 (h).For m > 1, from the recursive definition in (7), and taking into account that e A(h−s) B = ( 0  cos ( (h − s)) ) , it follows that and it is also immediate that m q ′ 22 (h) = m q 12 (h).
To complete the proof, we only need to show that m g 22 (h) =  − 2 m q 22 (h).For m = 1, we can directly compute the respective values.From the value of 1 q 22 (h) shown in (21), one gets which is the same that can be obtained for 1 g 22 (h), either from the expression in Corollary 1, or particularizing for m = 1 the series corresponding to m g 22 (h), Next, we will show that both sequences of functions, { m g 22 (h)} m≥1 and { m q 22 (h)} m≥1 , satisfy the sequence of initial value problems, 1 q 22 (h), which will complete the proof.
The initial conditions are obvious for both sequences.The series in ( 22) can be differentiated term-by-term, and one gets,for m > 1, □ From ( 17) and the previous Lemmas, the main result of this section follows immediately.Theorem 1.An exact numerical solution of ( 5)-( 6), with a = − 2 and b = , in the regular mesh t n = nh, with h = ∕N, for n ≥ −N,is given by X n (t) = F(t n ), for −N ≤ n ≤ 0, and Remark 2. We note that the results of Lemma 4 provide explicit expressions for the otherwise recursively defined matrices ) .
Example 1. Figure 1 presents a numerical example of application of the results of this section, showing the continuous solution given by Theorem 6 in Rodríguez et al. [13] (lines) and the exact numerical solution of Theorem 1 with N = 10 (points), for the problem ( 5)-( 6) with coefficients a = − 2 = −4, b =  = 0.5, delay  = 1 and initial function  (t) = (t + 1) 2 .It is also shown the first derivative of the solution (Fig. 1, (b)), with continuous values obtained by analytically differentiating the expression in Theorem 6 in Rodríguez et al. [13] and with numerical values given by the first component of X n in the exact numerical solution of Theorem 1.

HIGH ORDER NONSTANDARD FINITE DIFFERENCE SCHEMES
Although the exact numerical solution given by Theorem 1 avoids the infinite sums present in the solution of the general problem (3)-( 4) given in Theorem 2 in Castro et al. [4], recalled in (11), it still includes an integral term, R m (h), which in general needs to be numerically approximated, except for particular initial functions F(t) that could allow an exact computation.Therefore, in our next Theorem, we propose two families of nonstandard methods of as high precision as needed, by the same approach previously developed in García et al. [3] for coupled systems with commuting coefficients, computing in the first M intervals the exact solution, and then discarding the integral term in the following intervals, either computing the full sum in (24) or truncating it up to the M term.
Theorem 2. Let N ≥ 1 and h = ∕N.Fix M ≥ 1, and assume that the values of X n , for n = −N … MN, are computed using the exact scheme of Theorem 1.Then, for m > M, the numerical schemes  M and  M defined by have local error O(h 2M+1 ) and order 2M.
In the proof of Theorem 2, we will use the bounds given in the next Lemma.In what follows, || || denotes the maximum vector norm and the corresponding compatible norm for matrices.Lemma 5.The matrices G p (h) and R m (h) in (24) satisfy where Proof.From the expressions for the elements of G p (h) given in Corollary 1, and using the bound for Bessel functions [17, p. 227] is obtained.Similarly, from the expressions for the integrand in R m (h) given in Remark 2, one has Therefore, using a formula for the derivative of the Bessel functions [17, p. 222], and the same bounds as above, one gets and also and thus, ( 27) is obtained.
) for k ≤ n, which is the case for nh ≤ M, we will show that the scheme  M maintains this bound, and hence, a fortiori, also the scheme  M .
For m ≥ M + 1 and (m − 1) ≤ nh < m, from ( 24) and ( 25), one gets where the second summation is assumed to be empty for M = m−1.The first summation is O(h 2M+1 ) by the induction hypothesis.From the bounds for ||G p (h)|| obtained in Lemma 5, the second one, when it exists, is also O(h 2M+1 ), since the first term in the sum is the one corresponding to M + 1, with ||G p (h)|| = O(h 2(M+1)−1 ), and the rest of the terms are of lower order.Finally, also from Lemma 5, A schematic pseudocode for computing the numerical solutions given in Theorem 2 is shown in the Appendix A. Next figures illustrate the error analysis of the schemes defined in Theorem 2. Figure 2 (top) shows the errors of numerical solutions for Example 1, computed using  M schemes (left) and  M schemes (right), with M = 3, for three different mesh sizes.The respective errors relative to the expected order, that is, errors divided by h 6 , are shown in Figure 2 (bottom), with graphics overlapping, in agreement with the expected order given by Theorem 2.
Figure 3 presents the errors for numerical solutions of Example 1 computed using  M and  M schemes of three different orders.

DELAY-DEPENDENT STABILITY
In this section, we analyze the consistency of the numerical solutions resulting from applying the  M and  M schemes defined in Theorem 2 with respect to the delay-dependent asymptotic stability properties of the corresponding continuous solutions of problem ( 5)-( 6), which are recalled in the next Lemma.Lemma 6.The zero solution of ( 5) is asymptotically stable if and only if a < 0, b > 0, and there is a nonnegative k ∈ Z such that 2k < √ −a < (2k + 1), and b 2 < min ( −a 2 − (2k) 2 , (2k + 1) 2  2 + a 2 ) .When the zero solution of ( 5) is asymptotically stable for some  > 0, there are a finite number of stability switches between instability and stability as the delay increases, until a value of delay beyond which the solution remains unstable.
Proof.The results in this Lemma are particular cases of the results for the more general second order retarded and neutral delay differential equations analyzed in Cahlon and Schmidt [15] and Kuang [19,.□ Our next Theorem shows that the numerical solutions of problem ( 5)-( 6) obtained with the  M schemes defined in Theorem 2 preserve the stability properties of the continuous solution.
Theorem 3. Consider problem ( 5)- (6), with a = − 2 and b = , and the  M schemes defined in Theorem 2. The numerical solutions computed using  M schemes are asymptotically stable if and only if the zero solution of ( 5) is asymptotically stable.
Proof.The trivial solution of ( 5) is asymptotically stable if and only if all the roots  i of the characteristic equation have negative real part.
The system of difference equations that defines the  M scheme in (25) can be expressed in the form of a Volterra system of convolution type, by writing Hence, using the Z transform method, a necessary and sufficient condition for the asymptotic stability of ( 29) is [20,Theorem 5.21] where B(z) is the Z transform of the sequence {B r }, and, from the expressions in Corollary 1 and Remark 1, we can proceed to compute the elements of the matrix B(z) = ( bi ).
Using the Lommel's expansion for Bessel functions [18, p. 140] one has for b11 = b22 , since p g 11 (h) = p g 22 (h), with z = h and t ) . ( The other elements in B(z) can be immediately obtained, by noting that p g ′ 21 (h) = p g 11 (h) and p g ′ 22 (h) = p g 12 (h), and that the corresponding series can be term-by-term differentiated and integrated.In effect, using suitable formulae for the derivative of Bessel functions [17, p. 222], with the relation and with the relation, As a consequence, we deduce that and hence, P(z) = z 2 − 2 cos . Therefore, ) N, one has z N = e  for any root of P(z), and, since Nh = , or, equivalently,  2 +  2 − e − = 0, that is, (28).Hence, a root z of P(z) satisfies |z| ≥ 1 if and only if there is a root  of (28) with ℜ() ≥ 0. This completes the proof.□ In the proof of Theorem 3 it has been essential the possibility of expressing  M schemes in the form of Volterra systems of convolution type, which is not the case for  M schemes.Notwithstanding, although we cannot guarantee that  M schemes also preserve delay-dependent stability, numerical examples show that they seem to behave rather well in this aspect.
Figure 4 shows the long-term behavior of the numerical solutions computed with method  3 , with N = 5, for the same problem in Example 1 with increasing delay values.Eq. ( 5) with coefficients a = − 2 and b =  as in Example 1 represents a demanding test of the delay-dependent preserving properties of numerical methods, the solution becoming stable for  > 0 and then presenting seven stability switches, from stability to instability at  i = (2i + 1)∕ √  2 + , i = 0, … , 3, and from instability to stability at  ′ i = 2i∕ √  2 − , i = 1, … , 3. Thus, rounding to three decimal places, it is stable at  ∈ (0, 1.481) ∪ (3.358, 4.443) ∪ (6.717, 7.404) ∪ (10.075, 10.367), and unstable otherwise.As shown in the set of graphs in Figure 4, where the maximum absolute value of the numerical solution in each interval of  amplitude is represented for a sequence of delay values very close to the stability switches, the numerical solutions provided by the 3 scheme faithfully reflect the delay-dependent behavior of the continuous solutions.

CONCLUSIONS
Models based on DDEs, and specifically on second order DDEs, appear in numerous scientific and technical problems (see, e.g., Kolmanovskii and Myshkis [21] and references therein), and NSFD methods [22] have been proposed for different problems with delay (e.g., previous studies [23][24][25][26]).However, the construction of exact numerical schemes for DDEs, and NSFD methods derived from them, has been limited, with the exact scheme defined in Theorem 1 being the first, to our knowledge, for a second order DDE.
The  M and  M schemes defined in Theorem 2 provide numerical solutions for problem ( 5)-( 6) with as high accuracy as needed, similarly to the methods previously developed for first order scalar delay equations and systems with commuting coefficients [2,3], but with order O(h 2M ) instead of the O(h M ) order achieved in the previous problems.Although higher order schemes possess the obvious advantage of their greater precision, they also have some drawbacks, as a scheme of order 2M requires computing the exact solution for the mesh points in the first M intervals, as given in Theorem 1, which can be demanding for large M.
As shown in Theorem 3,  M schemes faithfully preserve the delay-dependent asymptotic stability of the exact solutions, a property difficult to prove and not common in numerical methods for DDEs [27].Although this property has not been proven for  M schemes, it has been shown by numerical examples that they may also reproduce consistently the dynamic behavior of the exact solutions, with the same accuracy as  M schemes and with some reduced computational effort.
The coefficients of the exact numerical solution in Theorem 1, presented in Corollary 1, are expressed in terms of Bessel and hypergeometric functions.It is suggestive that representations of solutions for problem (5)-( 6) connected to these special functions could also be obtained.

FIGURE 2
FIGURE 2 Absolute errors (log-scale) of the numerical solution for Example 1, computed using  M (A) and  M (B) schemes, with M = 3, for three different mesh sizes.(C, D): Errors divided by h 6 , that is, relative to the expected order from Theorem 2. [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 3
FIGURE 3 Errors (log-scale) of numerical solutions for Example 1, computed with three different  M (A) and  M (B) schemes, with h = 0.1.[Colour figure can be viewed at wileyonlinelibrary.com] FIGURE 4Maximum absolute values, in each interval of  amplitude, of the numerical solutions computed with method  3 for Example 1 with increasing values of delays, from (A) to (O).Unstable and stable behaviors (up and down trends, respectively) are consistent with the delay-dependent stability of the continuous solution.