**3. 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*<sup>M</sup>* and T*<sup>M</sup>* methods, depending on whether the full sum in (6) is kept or a truncated sum is used.

**Theorem 3.** *Let N* ≥ 1 *and h* = *τ*/*N. For a given M* ≥ 1*, assume that the values of Xn, for n* = −*N* ... *MN, are computed using the exact scheme of Theorem 2. Define FM and TM schemes to compute successive values for any m* > *M by the expressions*

$$\mathcal{F}\_M := \quad X\_{n+1} = \mathbf{e}^{Al} \sum\_{k=0}^{m-1} \frac{B^k h^k}{k!} X\_{n-kN\prime} \quad (m-1)\pi \le nh < m\pi,\tag{8}$$

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

$$\mathcal{T}\_M := \quad \mathcal{X}\_{n+1} = \mathbf{e}^{A\mathbf{h}} \sum\_{k=0}^M \frac{B^k h^k}{k!} \mathcal{X}\_{n-kN\prime} \quad (m-1)\tau \le n\hbar < m\tau. \tag{9}$$

*Then, both numerical schemes,* <sup>F</sup>*<sup>M</sup> and* <sup>T</sup>*M, have local error O*(*hM*+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*(*tk*) <sup>−</sup> *Xk* <sup>=</sup> *<sup>O</sup>*(*hM*+1) for *<sup>k</sup>* <sup>≤</sup> *<sup>n</sup>*, which is the case for *nh* <sup>≤</sup> *<sup>M</sup>τ*. Then, for *<sup>m</sup>* <sup>≥</sup> *<sup>M</sup>* <sup>+</sup> <sup>1</sup> and (*m* − 1)*τ* ≤ *nh* < *mτ*, using (6), one gets

$$\begin{split} \|X(t\_{n+1}) - X\_{n+1}\| &\leq \|\mathbf{e}^{At}\| \sum\_{k=0}^{M} \frac{\|B\|^k h^k}{k!} \|X(t\_{n-kN}) - 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. \end{split} \tag{10}$$

Let *<sup>M</sup>* <sup>=</sup> *B*, and *<sup>M</sup>*1, *<sup>M</sup>*<sup>2</sup> <sup>&</sup>gt; 0 such that e*As* <sup>&</sup>lt; *<sup>M</sup>*<sup>1</sup> and *F*(*s*) <sup>&</sup>lt; *<sup>M</sup>*<sup>2</sup> for *<sup>s</sup>* <sup>∈</sup> [0, *<sup>h</sup>*]. Then, by the induction hypothesis, the first term in (10) is *O*(*hM*+1) and the second term is bounded by

$$\frac{M^m M\_1 M\_2}{(m-1)!} \int\_{t\_n - m\tau}^{t\_n - m\tau + h} (t\_n - m\tau + h - s)^{m-1} ds < \frac{M^m M\_1 M\_2}{(m-1)!} h^m \le O(h^{M+1}).$$

Similar arguments result in the same bounds holding for the scheme F*M*.

**Remark 2.** *The results of Theorem 3 also hold if the values for Xn in the first intervals are computed using any numerical method of order at least O*(*hM*+1)*, instead of using the exact scheme. Although both types of schemes,* F*<sup>M</sup> and* T*M, have the same order, more general dynamic consistency properties can be proved for the class of* F*<sup>M</sup> schemes, as shown in the next section.*

The error analysis of the methods provided by Theorem 3 is illustrated in the next two figures. Errors of numerical solutions for Example 1, computed using T*<sup>M</sup>* schemes of three different orders, are shown in Figure 2 (top). The corresponding errors relative to the expected order, i.e., errors divided by *hM*, are shown in Figure 2 (bottom), with results in agreement with the expected orders given by Theorem 3.

Figure 3 presents the errors in relation with the size of the mesh for numerical solutions of Example <sup>1</sup> computed using <sup>T</sup>3, the truncated method with *<sup>M</sup>* <sup>=</sup> 3. Errors overlap when divided by *<sup>h</sup>*3, clearly showing that the method is of order three, as established in Theorem 3.

**Figure 2.** Absolute errors (log-scale) of numerical solutions for Example 1, computed using three different T*<sup>M</sup>* schemes of increasing orders, with *h* = 0.1. (**a**,**b**) Absolute errors for the first and second component, respectively. (**c**,**d**) Errors divided by *hM*.

**Figure 3.** Errors (log-scale) of numerical solutions for Example 1, computed with the method T<sup>3</sup> using three different mesh sizes. (**a**,**b**) Absolute errors for the first and second component, respectively. (**c**,**d**) Errors divided by *h*3.

## **4. Dynamic Consistency Properties**

In this section we analyse the consistency between dynamic properties of the numerical solutions resulting from applying the F*<sup>M</sup>* and T*<sup>M</sup>* schemes defined in Theorem 3 and the continuous solutions of problem (3)–(4).
