The celebrated Lax-Richtmeyer equivalence theorem (, , ) states that if the numerical scheme is stable then its global error is at least of the same order as its local truncation error. Indeed, it is the common experience that the frequently used methods produce solutions that have global error of the same order as the normalized local truncation error. Indeed, the Dahlquist’s Equivalence Theorem  states this fact for linear multistep methods. In fact, this behavior is so expected for all one-step and multistep methods that we typically define the order of a stable numerical method solely by the order conditions derived by Taylor series analysis of the local truncation error. This relationship between the order of the local truncation error and global error is also seen in finite difference schemes for partial differential equations (PDEs) [12, 29]. Work over the past decade reminds us that, while for a stable scheme the global error must at least of the same order as its local truncation error, it may in fact be of higher order.
In this paper we consider numerical solvers for ordinary differential equations (ODEs) of the form
(where we assume is at least continuous with respect to and Lipschitz-continuous with respect to , and typically much smoother). Without loss of generality , we focus our attention on the autonomous ODE
The canonical numerical method for this problem is the forward Euler method
where approximates the exact solution at time and we let . The forward Euler method has local truncation error and its approximation error at any given time defined by
and it produces a global error which is first order accurate
To develop methods that have higher order, we can add steps and define a linear multistep method , use multiple stages as in a Runge–Kutta methods , or include higher derivatives as we do in Taylor series methods . We can also combine the approaches above: the combination of multiple steps and stages is commonly see as in the general linear methods described in [2, 14], and multiple derivatives and stages have also been used [4, 5, 10, 23, 24, 27, 30]. In all these cases the goal is to increase the order of the local truncation error and therefore of the global error. All these approaches typically produce methods that have global errors that are of the same order as their local truncation errors.
It is, however, possible in some cases to devise schemes that have global errors that are higher order than predicted by the local truncation errors. Kulikov (for Nordsiek methods)  and by Weiner and colleagues  (for explicit two-step peer methods) followed along the lines of the quasi-consistency theory first introduced by Skeel in 1978  to find methods that have solution of order although their truncation errors are of order . In  Ditkowski and Gottlieb derived sufficient conditions for a family of general linear methods (GLMs) of any number of steps and stages under which one can control the accumulation of the local truncation error over time evolution and produce methods that are one order higher than expected from the truncation error alone. A similar theory was developed for implicit-explicit methods in . In , , and 
, it was shown that under additional conditions on the method the precise form of the first surviving term in the global error (the vector multiplying
) can be computed explicitly and leveraged for error estimation. In Ditkowski, Gottlieb, and Grant showed that under less restrictive conditions the form of the first surviving term in the global error can be computed explicitly and exploited to enhance the accuracy to order , using post-processing.
In this paper, we extend the error inhibiting approach to methods that have multiple steps, multiple stages, and two derivatives. In Section 3.1 we show that the theory developed in  carries through easily to such methods, and by imposing the same EIS+ conditions as in  on on the two-derivative GLMs we consider we can not only produce solutions with accuracy one order higher than predicted by truncation error analysis but also precisely describe the leading coefficients of the error term. This allows us to post-process the final time solution and obtain a global error that is two orders higher than the local truncation error. The construction of the post-processor is exactly the same as in  and we review it in Section 3.2. We proceed to devise error inhibiting two-step GLMs this new EIS approach, present these methods in Section 4, and in Section 5 we test a selection of these methods on numerical examples to demonstrate their enhanced accuracy properties.
2 Two-derivative methods
Multiderivative Runge–Kutta methods, particularly two-derivative methods, were first considered in [25, 38, 36, 32, 33, 15, 16, 23, 27, 4], and later explored for use with partial differential equations (PDEs) [30, 37, 22, 28, 8]. Recent interest in exploring the strong-stability properties of multiderivative Runge–Kutta methods, (also known as multistage multiderivative methods) is evident in [5, 10, 24]. A 2015 work by Okuonghae and Ikhile  discussed two- and three-derivatives methods with multiple stages and steps. These two-derivative general linear methods (GLMs) are similar to those in this work. (Also cite the paper Zack refereed and look in there for references)
In this work we consider the class of two-derivative general linear methods (GLMs) of the form
where is a vector of length that contains the numerical solution at times for :
The function is defined as the component-wise function evaluation on the vector :
and, similarly, the function is
For convenience, we select so that the first element in the vector approximates the solution at time , and the abscissas are non-decreasing . To initialize these methods, we define the first element in the initial solution vector and the remaining elements are computed using some highly accurate method.
We define the corresponding projection of the exact solution of the ODE (1) onto the temporal grid:
where and are the component-wise function evaluation on the vector . We define the global error as the difference between the vectors of the exact and the numerical solutions at some time
Methods of the form (2) have a local truncation error at time that is related to the approximation error
where the truncation error vectors have the form equationparentequation
Here, is the vector of abscissas and is the vector of ones, and the terms are understood component-wise . For a method to have order it must satisfy the order conditions
Note that the form (2) includes both explicit and implicit schemes, as appears on both sides of the equation. However, if and are both strictly lower triangular the scheme is explicit. We are only interested in zero-stable methods. A sufficient condition for this is that the coefficient matrix is a rank one matrix that has row sum one (i.e. satisfies the consistency condition). For simplicity we assume this to be the case in this the remainder of this work.
In this subsection we make some observations that will be useful in the remainder of the paper. In all the following we assume sufficient smoothness of all the quantities such as etc. Furthermore, we assume that he order conditions hold for all , and that we are in the asymptotic regime so that the errors are small and we can say that .111We will verify that this is reasonable in Lemma 2.2
Observation 1 We observe that
This is simply due to the fact that is smooth enough that it can be expanded:
where the error vector is , and we use the notation . Each term can be expanded as
so that, exploiting the smoothness of we have
which, due to the fact that becomes
Observation 2 Given the smoothness of , we observe that
The proof of this observation is exactly the same as above, with the assumption that is smooth enough to be expanded.
We use these observations to develop an equation that describes the growth of the error:
Given a zero-stable method of the form (2) which satisfies the order conditions
and where the functions and are smooth, the evolution of the error can be described by:
and use this equation to subtract from
using Observations 1 and 2 we have
Note that and are all scalars. Now recall that since , and that terms of the form and are all bounded, so we have
Now move the terms to the left side:
Noting that we have the desired result
We now turn to verifying that the assumption we made in the first paragraph of this subsection, that we can assume that the error is small:
The proof of this Lemma is given in , Lemma 3.
3 Designing multi-derivative error inhibiting schemes that can be post-processed to order
In this section we show how to construct two-derivative GLMs of the form (2) that are error inhibiting and so produce a solution of order even though one would only expect order from their truncation errors. We further show that under additional conditions we can express the exact form of the error and define an associated post-processor that allow us to recover order from a scheme that would otherwise be only th order accurate.
3.1 Two-derivative error inhibiting schemes that have order
In this section we consider a multi-derivative GLM of the form (2) where we assume that is a rank one matrix that satisfies the consistency condition so that this scheme is zero-stable. Furthermore, we assume that the coefficient matrices , , , , are such that the order conditions are satisfied for all , so that the method will give us a numerical solution that has error that is guaranteed of order . In the following theorem we show that if the truncation error vector lives in the null-space of the operator the order of the error can be shown to be of order . Furthermore, the theorem shows that under additional conditions on the coefficient matrices , , , , , we can determine precisely what the leading term of this error will look like and therefore remove it by post-processing.
Given a zero-stable two-derivative general linear method of the form
where is a rank one matrix that satisfies the consistency condition , and the coefficient matrices , , , , satisfy the order conditions
if the error inhibiting condition equationparentequation
|is satisfied, then the numerical solution produced by this method will have error|
|Furthermore, if the conditions|
are also satisfied, then error vector will have the more precise form:
Using Lemma 1 we obtain the equation for the evolution of the error
the fact is smooth assures us that so that we can expand the first term to obtain
The discrete Duhamel principle (given as Lemma 5.1.1 in ) states that given an iterative process of the form
where is a linear operator, we have
To analyze the error, we separate it into four parts:
and discuss each part separately:
The method is initialized so that the numerical solution vector is accurate enough to ensure that the initial error is negligible and we can ignore the first term.
The final term in the summation is . Recall that
so that this term contributes to the final time error the term
The term is a product of the operator
and the approximation error
so we obtain
due to the condition (14a) that states that .
Finally we look at the rest of the sum and use the boundedness of the operator to observe
Clearly, the first term here is a constant that depends only on the final time, so it is the product we need to bound. Using the definition of the operators
and the approximation error
Using the fact that in our case and that we obtain
The second term is eliminated by (14c)
All these parts together give us the result
We note that generalizing this approach from our work in  to include a second derivative was quite simple. This is due to the fact that the second derivative appears with a term multiplying it, so that the operator is in fact exactly the same for a two-derivative method of the form (2) as for the methods considered in . For this reason, expanding this approach to higher derivatives will be simple: the only change would be in the definition of the order conditions (10).
3.2 Designing a post-processer to recover order
At every time-step the error vector has the particular form (15)
so that the leading order term can be removed by post-processing at the end of the computation and a final time solution of order can be obtained, as we showed in . We notice that the form of the error for multi-derivative method (2) is exactly the same as that of the error of the EIS methods presented in , the only difference being the definition of the truncation error vector. The theoretical discussion of this post-processor is presented in  and we refer the reader to the presentation there. We review one such construction below:
Select the number of computation steps that will be used, requiring that
where is the number of steps and the truncation-error order of the time-stepping method.
Define the time vector of all the temporal grid points in the last computation steps:
Correspondingly, concatenate the final