Error Inhibiting Block One-Step Schemes for Ordinary Differential Equations

TL;DR

This paper introduces error inhibiting block one-step schemes for ODEs that control local error growth, achieving higher global accuracy than traditional methods, and demonstrates their effectiveness through constructed examples.

Contribution

It develops a methodology for constructing explicit error inhibiting block one-step methods that achieve higher global order accuracy by controlling local error growth.

Findings

Methods demonstrate higher order global error than local truncation error.

Constructed schemes show improved accuracy on test cases.

Theoretical framework for error inhibition is validated through numerical experiments.

Abstract

The commonly used one step methods and linear multi-step methods all have a global error that is of the same order as the local truncation error (as defined in \cite{gustafsson1995time,quarteroni2010numerical,AllenIsaacson,IsaacsonKeller,Sewell}). In fact, this is true of the entire class of general linear methods. In practice, this means that the order of the method is typically defined solely by the order conditions which are derived by studying the local truncation error. In this work, we investigate the interplay between the local truncation error and the global error, and develop a methodology which defines the construction of explicit {\em error inhibiting} block one-step methods (alternatively written as explicit general linear methods \cite{butcher1993a}). These {\em error inhibiting schemes} are constructed so that the accumulation of the local truncation error over time is controlled, which results in a global error that is one order higher than the local truncation error. In this work, we delineate how to carefully choose the coefficient matrices so that the growth of the local truncation error is inhibited. We then use this theoretical understanding to construct several methods that have higher order global error than local truncation error, and demonstrate their enhanced order of accuracy on test cases. These methods demonstrate that the error inhibiting concept is realizable. Future work will further develop new error inhibiting methods and will analyze the computational efficiency and linear stability properties of these methods.

Equations

$$u_t=F(t,u),t\ge 0$$

$$u(t=0)=u_0$$

$$v_{n+1}=v_n+\Delta tF(t_n,v_n).$$

$$v_{n+1}=\sum _{j=1}^s{a}_j{v}_{n+1-j}+\Delta t\sum _{j=0}^s{b}_{j}F(t_{n+1-j},v_{n+1-j}),$$

$$y_i$$

$$v_{n+1}$$

$$y_i$$

$$v_{n+1}^i$$

$$v_{n+1}=v_n+\Delta tF(t_n,v_n)+\frac{\Delta t^2}{2}{F}^{\prime }(t_n,v_n)+\frac{\Delta t^3}{3!}{F}^{\prime \prime }(v^n),$$

$$E_n=v_n-u(t_n),$$

$$\tau =\frac{u(t_{n+1})-u(t_n)-\Delta tF(t_n,u(t_n))}{\Delta t}\approx O(\Delta t).$$

$$u_t=f(t)u,t\ge 0$$

$$u(t=0)=u_0$$

$$v_{n+s}=\sum _{j=0}^{s-1}{a}_j{v}_{n+j}+\Delta t\sum _{j=0}^{s-1}{b}_{j}F(t_{n+j},v_{n+j})=\sum _{j=0}^{s-1}{a}_j{v}_{n+j}+\Delta t\sum _{j=0}^{s-1}{b}_{j}f(t_{n+j})v_{n+j}$$

$$U_n={\left(u(t_{n+s-1}),\dots ,u(t_n)\right)}^T$$

$$V_n={\left(v_{n+s-1},\dots ,v_n\right)}^T.$$

$$V_{n+1}=Q_n{V}_n$$

$$Q_{n}=\left(\begin{array}[]{cccc}a_{s-1}+{\Delta t}b_{s-1}f(t_{n+s-1})&\;\;a_{s-2}+{\Delta t}b_{s-2}f(t_{n+s-2})&\dots&\;\;a_{0}+{\Delta t}b_{0}f(t_{n})\\ I\\ &\ddots\\ &&I&0\end{array}\right).$$

$$V_n=S_{\Delta t}\left(t_n,t_{\nu }\right)V_{\nu },\nu \le n$$

$$S_{\Delta t}\left(t_n,t_{\nu }\right)=Q_{n-1}\dots Q_{\nu +1}{Q}_{\nu },S_{\Delta t}\left(t_n,t_n\right)=I.$$

$$\prod _{\mu =\nu }^{n-1}{Q}_{\mu }\equiv Q_{n-1}\dots Q_{\nu +1}{Q}_{\nu },\prod _{\mu =n}^{n-1}{Q}_{\mu }\equiv I.$$

$$S_{\Delta t}\left(t_n,t_{\nu }\right)=Q^{n-\nu }.$$

$${\Delta t}\,\mbox{\boldmath$\tau$}_{n}\,=\,U_{n+1}-Q_{n}U_{n}$$

$$\mbox{\boldmath$\tau$}_{n}\,=\,\left(\tau_{n},0,\ldots,0\right)^{T}.$$

$$E_n=U_n-V_n,$$

$$E_{n+1}\,=\,Q_{n}E_{n}\,+\,{\Delta t}\,\mbox{\boldmath$\tau$}_{n}\;.$$

$$E_{n}\,=\,S_{\Delta t}\left(t_{n},0\right)E_{0}\,+\,{\Delta t}\,\sum_{\nu=0}^{n-1}S_{\Delta t}\left(t_{n},t_{\nu+1}\right)\mbox{\boldmath$\tau$}_{\nu}\;,$$

$$E_{n}\,=\,\prod_{\mu=0}^{n-1}Q_{\mu}E_{0}\,+\,{\Delta t}\,\sum_{\nu=0}^{n-1}\left(\prod_{\mu=\nu+1}^{n-1}Q_{\mu}\right)\mbox{\boldmath$\tau$}_{\nu}\;.$$

$$\Vert S_{\Delta t}\left(t_n,t_{\nu }\right)\Vert \le K_s{e}^{{\alpha }_s\left(t_n-t_{\nu }\right)}$$

$$\left\|E_{n}\right\|\,\leq\,K_{s}\left[e^{\alpha_{s}t_{n}}\left\|E_{0}\right\|\,+\,\max_{0\leq\nu\leq n-1}\left\|\mbox{\boldmath$\tau$}_{\nu}\right\|\phi_{h}^{*}(\alpha_{s},\,t_{n})\right]\;.$$

$$\phi_{{\Delta t}}^{*}(\alpha_{s},\,t_{n})\,=\,{{\Delta t}}\,\sum_{\nu=0}^{n-1}e^{\alpha_{s}\left(t_{n}-t_{\nu+1}\right)}\,\approx\,\int_{0}^{t_{n}}e^{\alpha_{s}\left(t_{n}-\zeta\right)}d\zeta\,=\,\left\{\begin{array}[]{lll}\frac{e^{\alpha_{s}\,t_{n}}-1}{\alpha_{s}}&&\alpha_{s}\neq 0\\ t_{n}&&\alpha_{s}=0\end{array}\right.\;.$$

$${\mathbf{e}}_{n+s}=\Delta t\sum _{\nu =0}^{n-1}\left(1+O(\Delta t)\right){\tau }_{\nu }$$

Full text

††journal: Journal of Scientific Computing

1 Introduction

When solving an ordinary differential equation (ODE) of the form

[TABLE]

One can evolve the solution forward in time using the first order forward Euler method

[TABLE]

To obtain a more accurate solution, one can use methods with multiple steps:

[TABLE]

known as linear multistep methods [3]. Alternatively, one can use multiple stages, such as Runge–Kutta methods [3]:

[TABLE]

The class of general linear methods described in [2, 9] combines the use of multiple steps and stages, constructing methods of the form:

[TABLE]

The inclusion of multiple derivatives, such as Taylor series methods [3],

[TABLE]

is another possibility, and multiple stages and derivatives have also been developed and used successfully [17], [18], [11], [10], [4].

For time-dependent problems the global error, which is the difference between the numerical and exact solution at any given time $t_{n}=n\Delta t$:

[TABLE]

depends on the local truncation error which, roughly speaking, is the accuracy over one time step. In our case we define the local truncation error as the error of the method over one time-step, normalized by $\Delta t$. For example, the local truncation error for Euler’s method is (following [6, 13, 1, 8, 15])

[TABLE]

(To avoid confusion it is important to note that sometimes the truncation error is defined a little differently than we define it above and is not normalized by ${\Delta t}$).

A well known theoretical result, known as the Lax-Richtmeyer equivalence theorem (see e.g. [12], [6], [13]) states that if the numerical scheme is stable then the global error is at least of the same order as the local truncation error. In all the schemes for numerically solving ordinary differential equations (ODEs) that we are familiar with from the literature, the global errors are indeed of the same order as their local truncation errors111In the case where the truncation error is defined without the ${\Delta t}$ normalization the global error is one order lower than the truncation error.. This is common to other fields in numerical mathematics, such as for finite difference schemes for partial differential equations (PDEs), see e.g. [6, 13]. It was recently demonstrated, however, that finite difference schemes for PDEs can be constructed such that their convergence rates, or the order of their global errors, are higher than the order of the truncation errors [5]. In this work we adopt and adapt the ideas presented in [5] to show that it is possible to construct numerical methods for ODEs where the the global error is one order higher than the local truncation error. As we discuss below, these schemes achieve this higher order by inhibiting the lowest order term in the local error from accumulating over time, and so we name them Error Inhibiting Schemes.

Following an idea in [14], an interesting paper by Shampine and Watt in 1969 [16] describes a class of implicit one-step methods that obtain a block of $s$ new step values at each step. These methods take $s$ initial step values and generate the next $s$ step values, and so on, all in one step. These methods are in fact explicit block one-step methods, and can be written as general linear methods of the form (1) above. Inspired by this form, we construct explicit block one-step methods which are in the form (1), but where the matrix $\tilde{U}$ is an identity matrix, and the matrix $\tilde{A}$ is all zeros; these are known as Type 3 methods in [2]. The major feature of our methods is that in addition to satisfying the appropriate order conditions listed in [2], they have a special structure that mitigates the accumulation of the truncation error, so we obtain a global error that is one order higher than predicted by the order conditions in [2], which describe the local truncation error.

In Section 2 we motivate our approach by describing how typical multistep methods can be written and analyzed as block one-step methods: these methods obtain a block of $s$ new step values at each step. We show how this form allows us to precisely describe the growth of the error over the time-evolution. In Section 3 we then exploit this understanding to develop explicit error inhibiting block one-step methods that produce higher order global errors than possible for typical multistep methods. In Section 4 we present some methods developed according to the theory in Section 3 and we test these methods on several numerical examples to demonstrate that the order of convergence is indeed one order higher than the local truncation error. We also show that, in contrast to our error inhibiting Type 3 method, a typical Type 3 method developed by Butcher in [2] does not satisfy the critical condition for a method to be error inhibiting and therefore produces a global error that is of the same order as the local truncation error. Finally, we present our conclusions in Section 5, and suggest that further investigation of error inhibiting methods shall include the analysis of their linear stability properties, storage implications, and computational efficiency.

2 Motivation

In this section we present the analysis of explicit multistep methods in a block one-step form for a simple linear problem. In this familiar setting we define the local truncation error, the global error, and the solution operator that connects them. We also discuss the stability of a method of this form. We limit our analysis to the linear case so that we can clearly observe the process by which the solution operator interacts with the local truncation error, and results in a global error that is of the same order as the local truncation error. Although we are dealing for the moment with standard multistep methods, this will set the stage for the construction and analysis of error inhibiting block one-step methods.

In order to illustrate the main idea we start with a linear ordinary differential equation (ODE)

[TABLE]

where $f(t)<M\,,\;\forall t\geq 0$ and $f(t)$ is analytic.

An $s$-step explicit multistep method applied to (4) takes the form

[TABLE]

where the time domain is discretized by the sequence $t_{n}=n\,{\Delta t}$, and $v_{n}$ denotes the numerical approximation of $u(t_{n})$. The method (5) is defined by its coefficients $\{a_{j}\}_{j=0}^{s-1}$ and $\{b_{j}\}_{j=0}^{s-1}$, which are constant values.

Following [6] we rewrite the method (5) in its block form. To do this, we first introduce the exact solution vector

[TABLE]

and similarly, the numerical solution vector is

[TABLE]

Now (5) can be written in block form so that it looks like a one step scheme

[TABLE]

where

[TABLE]

From repeated applications of equation (8) we observe that the numerical solution vector $V_{n}$ at any time $t_{n}$ can be related to $V_{\nu}$ for any previous time $t_{\nu}$

[TABLE]

where $S_{\Delta t}$ is the discrete solution operator. This operator can be expressed explicitly by

[TABLE]

For simplicity we can define this by

[TABLE]

Note that if each matrix $Q_{\mu}$ is independent of $\mu$ (in other words, in the constant coefficient case where $f$ is independent of $t$), we simply have a product of matrices $Q$, and the discrete solution operator becomes

[TABLE]

The behavior of a method depends in large part on the accuracy of its solution operator. We begin by defining the local truncation error as the error of the method over one time-step, normalized by ${\Delta t}$:

Definition 1: The local truncation error \mbox{\boldmath\tau}_{n} is given by [6, 13, 1, 8, 15]

[TABLE]

Note that in the case of the standard multistep method, where $Q_{n}$ is given by the matrix (9), the truncation error has only one non-zero element:

[TABLE]

The error that we are most interested in is the difference between the exact error vector and the numerical error vector at time $t_{n}$,

[TABLE]

known as the global error. At the initial time, we have the error $E_{0}$ which is based on the starting values a method of this sort requires: the values $v_{j}$, $j=0,\ldots,s-1$ that are prescribed or somehow computed. Typically, $v_{0}$ is the initial condition defined in (1) and $v_{j}$, $j=1,\ldots,s-1$ are computed to sufficient accuracy using some other numerical scheme. Thus, the value of $E_{0}$ is assumed to be as small as needed.

The evolution of the global error (16) depends on the local truncation error defined by (14) and the discrete solution operator given in (8):

[TABLE]

Unraveling this equality all the way back to $E_{0}$ gives

[TABLE]

or, equivalently

[TABLE]

(This formula is obtained from the discrete version of Duhamel’s principle, see Lemma 5.1.1 in [6]).

It is clear from (18) that the behavior of the discrete solution operator $S_{\Delta t}(t_{n},t_{\nu+1})$ must be controlled for this error to converge. This property defines the stability of the method. Also here we use the stability definition presented in [6], namely:

Definition 2: The scheme (8) is called stable if there are constants $\alpha_{s}$ and $K_{s}$, independent of ${\Delta t}$, such that for all $0<{\Delta t}\leq{\Delta t}_{0}$

[TABLE]

If the scheme is stable, we can use (20) and (18) to bound the growth of the error:

[TABLE]

where

[TABLE]

Equation (21) means that stability implies convergence:222 For partial differential equations this result is known as one part of the celebrated Lax-Richtmeyer equivalence theorem. See e.g. [12], [6], [13]. if the scheme is stable than the global error is controlled by the local truncation error for any given final time. In the formula above it is clear that the global error must have order at least as high as the local truncation error, but the possibility of having a higher order global error is left open.

The first Dahlquist barrier [7, 3] states that any explicit $s$ step linear multistep method can be of order $p$ no higher than $s$. It is the common experience that methods have global error of the same order as the local truncation error. These two together greatly limit the accuracy of the methods we can derive.

Remark 1

In an Adams-Bashforth scheme the entry in the first row and first column in the term $S_{\Delta t}\left(t_{n},t_{\nu}\right)=\prod_{\mu=\nu}^{n-1}Q_{\mu}$ is equal to $1+O({\Delta t})$. Therefore the error, due to the accumulation of the contributions from the truncation errors, becomes:

[TABLE]

which is approximately the average value of $\tau_{\nu}$ over ${\nu=0,..,n-1}$. This suggests that we may need to look outside the family of linear multistep methods to attain a higher order global error.

The analysis in this section suggests that if the operator $Q_{n}$ is properly constructed, the growth of the global error described in Equation (19) may be controlled through the properties of the operator $Q_{n}$ and its relationship with the local truncation error \mbox{\boldmath\tau}_{n}. However, as implied by the example of the Adams-Bashforth scheme above, we need to construct methods where the operator $Q_{n}$ is not limited by the structure in this section. In the next section we present the construction of block one-step methods that are error inhibiting. The class of methods described by this block one-step structure is very broad: while all classical multistep methods can be written in this block form, not every such block one-step method can be written as a classical multistep method. Thus, we rely on the discussion in this section with one main change: the structure of the operator $Q_{n}$.

3 An Error Inhibiting Methodology

In Section 2 we rewrote explicit linear multistep methods in a block one-step form, and expressed the relationship between its local and global error. We observed that the growth of the local errors is driven by the behavior of the discrete solution operator $Q_{n}$, and in particular its interaction with the local truncation error. Using this insight we show in this section that it is possible to construct such explicit block one-step methods (which are also known as Type 3 DIMSIM methods in [2]) that inhibit the growth of the truncation error so that the global error (16) gains an order of accuracy over the local truncation error (14).

We begin in Section 3.1 by describing the construction and analysis of error inhibiting block one-step schemes for the case of linear constant coefficient equations. We then show that this approach yields methods that are also error inhibiting for variable coefficient linear equations in Section 3.2 and nonlinear equations in Section 3.3.

3.1 Error inhibiting schemes for linear constant coefficient equations

Given a linear ordinary differential equation with constant coefficients:

[TABLE]

where $f\in\mathbbm{R}$. We define a vector of length $s$ that contains the exact solution of (24) at times $\left(t_{n}+j\Delta t/s\right)$ for $j=0,\ldots,s-1$

[TABLE]

and the corresponding vector of numerical approximations

[TABLE]

Note that although we are assuming that the solution $u$ at any given time is a scalar, this entire discussion easily generalizes to the case where $u$ is a vector, with only some cumbersome notation needed. Thus without loss of generality we continue the discussion with scalar notation.

Remark 2

The notation above emphasizes that this scheme uses $s$ terms for generating the next $s$ terms, unlike the explicit linear multistep methods in the section above which use $s$ terms to generate one term. To match with the notation in Section 2 above, we can replace $\Delta t^{\prime}=s\Delta t$ thus defining this scheme on integer grid points.

We define the block one-step method for the linear constant coefficient problem (24)

[TABLE]

where

[TABLE]

and $A,B\in\mathbbm{R}^{s\times s}$. Unlike in the case of classical multistep methods, here we do not restrict the structure of the matrices $A$ and $B$. Thus, any multistep method of the form (5) can be written in this form (as we saw above), but not every method of the form (28) can be written as a multistep method. In fact, this methods is a general linear method of the DIMSIM form (1) with $\tilde{A}$ is all zeroes, $\tilde{U}$ is the identity matrix, $\tilde{V}=A$, and $\tilde{B}=B$. This particular formulation is, as we mentioned above, called a Type 3 DIMSIM in Butcher’s 1993 paper [2].

At any time $t_{n}$, we know that $\;u(t_{n}+{\Delta t})\,=\,u(t_{n})+O({\Delta t})$, so that for the numerical solution $V_{n}$ to converge to the analytic solution $U_{n}$ one of the eigenvalues of $Q$ must be equal to $1+O({\Delta t})$, and its eigenvector must have the form:

[TABLE]

The structure of the eigensystem of $A$, which is the leading part of $Q$, is critical to the stability of the scheme and the dynamics of the error.

Suppose $A$ is constructed such that:

${\bf rank}(A)=1$. 2. 2.

Its non-zero eigenvalue is equal to one and its corresponding eigenvector is $\left(1,\ldots,1\right)^{T}$ 3. 3.

$A$ can be diagonalized.

Property (2) assures that the method produces the exact solution for the case $f=0$. Now, since the term ${\Delta t}Bf$ is only an $O({\Delta t})$ perturbation to $A$, the matrix $Q$ will have one eigenvalue, $z_{1}=1+O({\Delta t})$ whose eigenvector has the form

[TABLE]

and the rest of the eigenvalues satisfy $z_{j}=O({\Delta t})$ for $j=2,\ldots,s$.

Since the $\|Q\|=1+O({\Delta t})$, we can conclude that there exist constants $K_{s}$ and $\alpha_{s}$ such that

[TABLE]

where $\alpha_{s}=\|B\|\,|f|$. Therefore, according to Definition 2, the scheme (27) is stable. By the same argument used above, we can show that the global error will have order that is no less than the order of the local truncation error.

We now turn to the task of investigating the truncation error, \mbox{\boldmath\tau}_{n}. The definition of the local truncation error in this case is still

[TABLE]

as defined in the previous section in Equation (14).

Remark 3

Since $Q=A+\Delta tBf$ and $u_{t}=fu$ the local truncation error can be written as

[TABLE]

Therefore \mbox{\boldmath\tau}_{n} does not explicitly depend on $f$. This observation is valid for the variable coefficients and the nonlinear case as well.

The definition of the error is

[TABLE]

as in Equation (16). The evolution of the error is still described by Equation (19)

[TABLE]

which in the linear constant coefficient case becomes

[TABLE]

The main difference between this case and the linear multistep method in Section 2 is that the structure of $Q$ is different, and that unlike (15), in this case all the entries in \mbox{\boldmath\tau}_{n} are typically non-zero.

Equation (32) indicates that there are several sources for the error at the time $t_{n}$:

The initial error $E_{0}$ which is the error in the initial condition $V_{0}$: This error is caused primarily by the numerical scheme used to compute the first $s-1$ elements in $V_{0}$. We assume these errors can be made arbitrary small. The initial value, which is the final element of $V_{0}$, is taken from the analytic initial condition and is considered to be accurate to machine precision. 2. 2.

The term {\Delta t}\,\mbox{\boldmath\tau}_{n-1}, which is the last term in the sum in the right hand side of (32): This term is clearly, by definition, of the size O({\Delta t})\|\mbox{\boldmath\tau}_{n-1}\|. 3. 3.

The summation

[TABLE]

which are all the rest of the terms in the sum in the right hand side of (18): This is the term we need to bound to control the growth of the truncation error.

The terms in the sum (33) are all comprised of the discrete solution operator $Q$ multiplying the local truncation error. This leads us to the major observation that is the key to constructing error inhibiting methods: if the local truncation error lives in the subspace of eigenvectors that correspond to the eigenvalues of $O({\Delta t})$, then the growth of the truncation error will be inhibited, and the global error will be one order higher than the local truncation error.

Recall that $Q$ has one dominant eigenvalue that has the form $1+O({\Delta t})$ and all the others are $O({\Delta t})$. Correspondingly, two subspaces can be defined

[TABLE]

where $\psi_{j}$ is the eigenvector associated with each eigenvalue $z_{j}$. As $\psi_{j}$ can be normalized, we assume that $\|\psi_{j}\|=O(1)$. It should be noted that while $\Psi_{1}$ and $\Psi_{1}^{c}$ are linearly independent, they are not orthogonal subspaces. Furthermore, since the matrix $A$ is diagonalizable by construction, its eigenvectors span $\mathbbm{R}^{s}$. Since \mbox{\boldmath\tau}_{\nu}\in\mathbbm{R}^{s}, it can be written as

[TABLE]

where $\gamma_{1}\psi_{1}\in\Psi_{1}$ and $\sum_{j=2}^{s}\gamma_{j}\psi_{j}\in\Psi_{1}^{c}$.

Of course, the truncation error \mbox{\boldmath\tau}_{\nu} is determined by the entries of $Q$. To ensure that the local truncation error is mostly in the space $\Psi_{1}^{c}$ of eigenvectors which correspond to the eigenvalues of size $O({\Delta t})$, we choose the entries of $Q$ (i.e. the entries of $A$ and $B$) such that $\gamma_{1}=O({\Delta t})$, which will mean that

[TABLE]

Using this, we can bound product of the discrete solution operator and the truncation error,

[TABLE]

where $z_{j}$ are the eigenvalues of $Q$. Therefore we have

[TABLE]

Whenever the condition (36) is satisfied, we can show that the sum (33) above is bounded:

[TABLE]

(Recall (22) for the definition of of $\phi_{{\Delta t}}^{*}(c,\,T)$.)

In the final equation, $T$ is the final time, and the term $\phi_{{\Delta t}}^{*}(c,\,T)$ is therefore a constant. Thus we have the bound

[TABLE]

Putting this all together into (32), we obtain

[TABLE]

Thus, if the coefficients of $A$ and $B$ are chosen so that we can control the size of \left\|Q\mbox{\boldmath\tau}_{\nu}\right\| in (36), we can obtain a scheme that inhibits the growth of the local truncation error, so that the global error is one order more accurate than its truncation error.

3.2 Linear variable-coefficient equations

In the previous section we showed how to construct an error inhibiting method by choosing the coefficients in $A$ and $B$ so that the local truncation error lives (mostly) in the space that is spanned by the eigenvectors corresponding to eigenvalues that are of $O({\Delta t})$. In this section we show that under the same criteria as above, these methods are also error inhibiting when applied to a variable coefficient linear ordinary differential equation:

[TABLE]

where $f(t)$ assumed to be analytic or as smooth as needed, and bounded. In this case the scheme is given by a time-dependent evolution operator $Q_{n}$ which may change each time-step:

[TABLE]

where

[TABLE]

and the matrices $A$ and $B$ are the same as described above for the constant coefficient scheme.

Since $f(t)$ is an analytic function, $Q_{n}$ can be written as

[TABLE]

We can also say then that

[TABLE]

Each $\tilde{Q}_{n}$ has the same structure as $Q$ in the constant coefficient case. In particular

[TABLE]

Furthermore, as was pointed out in Remark 3, since the local truncation error \mbox{\boldmath\tau}_{n} does not depend explicitly on $f(t)$ at any time $t_{n}$, we can write \mbox{\boldmath\tau}_{n} as a linear combination of the eigenvectors of $A$ that correspond to the zero eigenvalues. Thus, \mbox{\boldmath\tau}_{n} lives (mostly) in the space that is spanned by the eigenvectors of any matrix $\tilde{Q}_{n}$ corresponding to eigenvalues that are of $O({\Delta t})$. We can then follow the same analysis as in (35)–(3.1), to obtain the bound

[TABLE]

In this case, Equation (18) takes the modified form (for $n\geq 1$)

[TABLE]

The first term is negligible because we assume that the initial error can be made arbitrarily small, and the final term is clearly of order {\Delta t}\mbox{\boldmath\tau}_{n-1}. Using (45), (46) and the same analysis as in (35)–(38) we have

[TABLE]

Putting these all together we have

[TABLE]

This simple proof shows that even for the variable coefficient case, the schemes constructed as described above have a higher order error than would be expected from the truncation error. In the next subsection we extend this analysis to the general nonlinear case.

3.3 Nonlinear equations

Finally, we analyze the behavior of methods satisfying the assumptions in Section 3.1 when applied to nonlinear problems. Consider the nonlinear equation

[TABLE]

where $f(u,t)$ assumed to be analytic in $u$ and $t$. We now use the scheme

[TABLE]

where the matrices $A$ and $B$ are as constructed above for the constant coefficients problem.

As defined in (14), the exact solution to (3.3) and the truncation error are related by

[TABLE]

Note that by Taylor expansion

[TABLE]

where$f_{u}(u,t)=\partial f(u,t)/\partial u$ and $|r(v_{\nu}-u_{\nu})|\leq c_{1}|v_{\nu}-u_{\nu}|^{2}$. Subtracting (49) from (50) and assuming that $E_{n}=U_{n}-V_{n}\ll 1$ gives

[TABLE]

where $\|R(E_{n})\|\leq c_{1}\|E_{n}\|^{2}$. Equation (51) means that as long as O(E_{n}^{2})\ll O(\mbox{\boldmath\tau}_{n}), the equation for the error $E_{n}$ can be analyzed in essentially the same way as for the linear variable coefficient case, and the same estimates hold.

In order to evaluate the time interval in which O(E_{n}^{2})\ll O(\mbox{\boldmath\tau}_{n}) we note that although the term $R(E_{n})$ in (51) is not a non-homogeneous term but rather a function of $E_{n}$, we can still use the approach used in [6, Theorem 5.1.2]) to prove stability for a perturbed solution operator. As in [6, Theorem 5.1.2]), we use the discrete version of Duhamel’s principle to obtain

[TABLE]

where

[TABLE]

Taking the norm of (3.3) and using the triangle inequality we obtain

[TABLE]

As in the linear case we assume that the initial error, $E_{0}$ is arbitrary small, so the first term is negligible. If $\hat{Q}_{\nu+1}$ is constructed such that \|\hat{Q}_{\nu+1}\mbox{\boldmath\tau}_{\nu}\|={\Delta t}O(\mbox{\boldmath\tau}_{\nu}) then using the same analysis as in variable coefficient case the second term in (3.3) is less or equal to {\Delta t}c_{0}\phi_{h}^{*}(c,\,t_{n})\max_{\nu=0,...,n-1}\left\|\mbox{\boldmath\tau}_{\nu}\right\|. As to the third term, the same arguments can be used to show that it is bounded by

[TABLE]

so that (3.3) (with the substitution of (55) for the final term) can be re-arranged to obtain

[TABLE]

If $c_{1}\phi_{h}^{*}(c,\,t_{n})\left\|E_{n}\right\|<1/2$, we obtain

[TABLE]

This estimate holds as long as

[TABLE]

which is satisfied for all times $t_{n}$ such that ${\Delta t}\phi_{h}^{*}(c,\,t_{n})=O(1)$.

Therefore

[TABLE]

for the nonlinear case as well.

4 Some Error Inhibiting Explicit Schemes

In the previous section we define sufficient conditions for methods of the form

[TABLE]

where

[TABLE]

to be error inhibiting. These are

1. C1.

${\bf rank}(A)=1$. 2. C2.

Its non-zero eigenvalue is equal to $1$ and its corresponding eigenvector is

[TABLE] 3. C3.

$A$ can be diagonalized. 4. C4.

The matrices $A$ and $B$ are constructed such that when the local truncation error is multiplied by the discrete solution operator we have the bound:

[TABLE]

This is accomplished by requiring the local truncation error to live in the space of the eigenvectors of $A$ that correspond to the zero eigenvalues.

In this section we present several schemes which were constructed using the approach presented in the previous section. In Section 4.1, we present a block one-step method that evolves two steps ($v_{n}$ and $v_{n+\frac{1}{2}}$) to obtain the next two steps ($v_{n+1}$ and $v_{n+\frac{3}{2}}$). This method has truncation error (14) that is second order, while its global order (16) is third order. We demonstrate that the expected convergence rate is attained on several sample nonlinear problems. In this section we also show that a typical Type 3 DIMSIM method (derived in [2]) that satisfies the first three conditions above but not the fourth, has truncation error of order two, and its global error is of the same order. This demonstrates the importance of condition C4.

Next, in Section 4.2 we present a block one-step method that evolves three steps $v_{n}$, $v_{n+\frac{1}{3}}$ and$v_{n+\frac{2}{3}}$ to obtain $v_{n+1}$, $v_{n+\frac{4}{3}}$ and $v_{n+\frac{5}{3}}$. This method has truncation error (14) that is third order, while its global order (16) is fourth order, as we demonstrate on several sample problems. Finally, to show that the methods in each class are not unique, we present two other methods of this type and show that their global error is of one order higher than the local truncation error on a sample nonlinear system.

4.1 A third order error inhibiting method with $s=2$.

In this subsection we define an explicit block one-step with $s=2$ that satisfies the conditions C1 – C4 above. This method takes the values of the solution at the times $t_{n}$ and $t_{n+\frac{1}{2}}$ and obtains the solution at the time-level $t_{n+1}$ and $t_{n+\frac{3}{2}}$. The exact solution vector for this problem is

[TABLE]

and, similarly, the corresponding vector of numerical approximations is

[TABLE]

The scheme is given by:

[TABLE]

and has truncation error

[TABLE]

The matrix $A$ can be diagonalized as follows:

[TABLE]

Observe that the leading order of the truncation error (61) is in the space of the second eigenvector of $A$, the one that corresponds to the zero eigenvalue. Also, as was pointed out in Remark 3, \mbox{\boldmath\tau}_{n} depends only on this eigenvector of $A$ and a multiple that is not directly dependent on $f$ but only on the third derivative of the solution $u$. This underscores the analysis in Sections 3.2 and 3.3 that demonstrates that the error inhibiting property carries through for variable coefficient and nonlinear problems.

To study the behavior of the global error we use the fact shown in Section 3.3 that even for a nonlinear equation it is sufficient to analyze the matrix

[TABLE]

where $f$ is a constant. In this case:

[TABLE]

Recall that, neglecting the initial error $E_{0}$, we can say that the global error is (16)

[TABLE]

Putting together equations (61) and (72) we see that each term Q\,\mbox{\boldmath\tau}_{\nu} contributes to the error in two ways:

• •

The first contribution is due to the fact that \mbox{\boldmath\tau}_{\nu} is almost co-linear with the second eigenvector $\psi_{2}$. The order of this contribution is

[TABLE]

where the term $|z_{2}|$ is the second eigenvalue which is of order $O({\Delta t})$.

• •

The second contribution to the error comes from the component of \mbox{\boldmath\tau}_{\nu} that is a multiple $\gamma_{1}$ of the first eigenvector $\psi_{1}$,

[TABLE]

the term $\gamma_{1}$ is of $O({\Delta t})$ because \mbox{\boldmath\tau}_{\nu} lives mostly in the space of $\psi_{2}$.

While each of the terms in {\Delta t}Q\,\mbox{\boldmath\tau}_{\nu} has order O({\Delta t}^{2})\cdot O(\|\mbox{\boldmath\tau}_{\nu}\|)=O({\Delta t}^{4}), as the method is evolved forward, the errors accumulate over time, and sum of all contributions from all the times gives us a global error of order O({\Delta t})\cdot O(\|\mbox{\boldmath\tau}_{n}\|)=O({\Delta t}^{3}).

Example 1a: To demonstrate that this method indeed performs as designed we study its behavior on a nonlinear scalar equation of the form:

[TABLE]

We evolve the solution of this equation to time $T=1$ using the scheme (60). The initial steps are computed exactly. The plots of the errors and the truncation errors are presented in Figure 1(a). Both errors are shown for the first component, $v_{n+1/2}$ (denoted v(1) in the legend) and the second component, $v_{n}$ (denoted v(2) in the legend). Clearly, although the truncation error is only second order (denoted tr err v(1) and tr err v(2) in the legend), the global error is third order, as predicted by the theory.

Example 1b: It is important that the method will perform as designed on a nonlinear system as well. To demonstrate this, we solve the the van der Pol system

[TABLE]

using the same scheme (60). As this is a system, it is important that both components are examined. Thus, the vector of the numerical solution has two components for the time level $t_{n}$, denoted by v(2), and two components for the time level $t_{n+\frac{1}{2}}$, denoted by v(1). In Figure 1(b) the convergence plot of the components of $u^{(1)}$ and $u^{(2)}$ are presented. Once again, we see that the convergence rate is indeed third order.

Remark 4

It is important to note that not all Type 3 DIMSIM methods have the EIS property! The property that the local truncation error lives in the space spanned by the eigenvectors of $A$ that correspond to the zero eigenvalues is needed for the error inhibiting behavior to occur, and this property is not generally satisfied. To observe this, we study the DIMSIM scheme of types 3 presented by J. C. Butcher in [2].

Consider the scheme

[TABLE]

given in [2]. This scheme has truncation error

[TABLE]

The matrix $A$ can be diagonalized as follows:

[TABLE]

The truncation error \mbox{\boldmath\tau}_{n} can be written as a linear combination of the two eigenvectors of $A$ as follows:

[TABLE]

Unlike the EIS scheme (60), here the first term in this expansion is of the order of O(\mbox{\boldmath\tau}_{n})=O(\Delta t^{2}). Therefore a term of the order of \Delta tO(\mbox{\boldmath\tau}_{n})=O(\Delta t^{3}) is accumulated at each time step, so that the global error is second order.

We note that both this method (76) and our error inhibiting method (60) satisfy the order conditions in Theorem 3.1 of [2] only up to second order ($p=2$). However, as we see in Figure 2, when the method (76) is used to simulate the solution of the problems (4.1) and (4.1) we have second order accuracy, while the error inhibiting method (60) gave third order accuracy (Figure 1).

4.2 A fourth order error inhibiting method with $s=3$.

In this subsection we present an error inhibiting method with $s=3$ that takes the values of the solution at the times $t_{n}$, $t_{n+\frac{1}{3}}$, and $t_{n+\frac{2}{3}}$ and uses these three values to obtain the solution at the time-level $t_{n+1}$, $t_{n+\frac{4}{3}}$, and $t_{n+\frac{5}{3}}$. The exact solution vector is given by

[TABLE]

and the corresponding vector of numerical approximations is

[TABLE]

Consider the error inhibiting scheme

[TABLE]

which has a local truncation error of third order,

[TABLE]

However, it can be verified that for the linear case, the product

[TABLE]

Given the analysis in Section 3.3 above, this result will carry over to the nonlinear case, and thus this method will have a fourth order global error, despite the third order truncation error.

To demonstrate this result we revisit the two examples (4.1) and (4.1) in the previous subsection and use the scheme (90) to evolve them forward in time. The results, shown in Figure 3, are exactly as we expect: although the truncation errors (seen for the problem (4.1) in Figure 3(a)) are only third order, the errors are fourth order for both problems (4.1) and the van der Pol problem (4.1).

4.2.1 Other fourth order error inhibiting methods with $s=3$.

The methods above are not unique, in fact other methods can be derived using this approach. In this section we present two additional error inhibiting methods with $s=3$ that have local truncation error that is third order but demonstrate fourth order global error on a nonlinear system.

The first method is

[TABLE]

and has a local truncation error of third order,

[TABLE]

The second method is

[TABLE]

The truncation error is also third order

[TABLE]

Both these methods satisfy

[TABLE]

as well. As above, this property results in an error inhibiting mechanism that produced a global error of order four. This can be seen once again in Figure 4, using the nonlinear problem (4.1) above. The results of method (109) are on the left and of (127) are on the right.

5 Conclusions

While it is generally assumed that the global error will be of the order of the local truncation error, in this work we presented an approach to creating methods that have a global error of higher order than predicted by the local truncation error. To accomplish this, we used the block formulation of a method $V_{n+1}=Q_{n}V_{n}$ where the discrete solution operator $Q_{n}=A+{\Delta t}BF_{n}$ is comprised of matrices of coefficients $A$ and $B$, and the matrix operator $F_{n}.$

We show that if $A$ is a diagonalizable matrix of rank one, that has only one nonzero eigenvalue, $z_{1}=1$, that corresponds to the eigenvector of all ones, then the error inhibiting property will occur if the leading part of the local truncation error error for the linear constant coefficient case ($F_{n}=F=$ a constant) is spanned by the eigenvectors corresponding to the zero eigenvalues of $A$ (to the leading order). We show that a method that has these properties will have a global error that has higher order than the local error, on nonlinear problems.

After presenting the concept behind these methods we use the theoretical properties above to develop block one-step methods that are in the family of Type 3 DIMSIM methods presented in [2]. We demonstrate in numerical examples on nonlinear problems (including a nonlinear system) that these methods have global error that is one order higher than the local truncation errors. We also show that this is in contrast to another Type 3 DIMSIM method which has a matrix $A$ that satisfies the first three properties C1 – C3, but does not satisfy the error inhibiting property C4, that the local truncation error is in the space spanned by the eigenvectors of $A$ that correspond to the zero eigenvalues, and indeed does not give us a global error that is higher than the local truncation error on nonlinear test problems.

The major development in this work is the concept of an error inhibiting method and the new approach for developing methods that are constructed to control the growth of the local truncation error. While the newly developed methods presented in this work can be used in place of currently standard methods (particularly in place of type 3 DIMSIM methods) to obtain higher order accuracy, it is not yet known how they compare to other methods in terms of other important properties. In future work we intend to the study of the computational efficiency and storage requirements of these methods and the analysis of their linear stability regions. We expect that this will also lead to further development of error inhibiting methods that have other favorable properties.

Acknowledgements: *The authors wish to thank Professor John Butcher for a very helpful discussion, and in particular for his valuable advice on general linear methods, especially the Type 3 DIMSIM methods.

The work of Sigal Gottlieb was supported by AFOSR grant FA9550-15-1-0235.*