Efficient high order numerical methods for evolving the solution of an ordinary differential equation (ODE) are frequently desired, especially for the time-evolution of PDEs that are semi-discretized in space. Higher order methods are desirable as they enable more accurate solutions for larger step-sizes. In this paper we consider numerical solvers for ordinary differential equations (ODEs)
(where is sufficiently smooth). Any non-autonomous system of this form can be converted to an autonomous system (see ) so that without loss of generality, we restrict ourselves to the autonomous case:
The forward Euler method is the simplest numerical method for evolving this problem forward in time with time-steps :
where approximates the exact solution at time and we let . For this method, we have
where is called the local truncation error, and the approximation error, at any given time . At any such time we are interested in is the global error which is the difference between the numerical and exact solution. For the forward Euler method the global error is first order accurate:
Observe that the global error and the local truncation error are of the same order, . This is the common experience that one step methods and linear multi-step methods that are typically used produce solutions that have a global error that is of the same order as the local truncation error.
To generate methods with higher local truncation errors, we can generalize the forward Euler scheme by adding steps (e.g. linear multistep methods), stages (e.g. Runge–Kutta methods), or derivatives (e.g. Taylor series methods). The Dahlquist Equivalence Theorem  states that any zero-stable, consistent linear multistep method with truncation error will have global error , provided that the solution has at least smooth derivatives. This behavior is usually seen in other types of schemes, and is the general expectation: so much so that the order of a numerical method is typically defined solely by the order conditions derived by Taylor series analysis of the local truncation error. However, the Lax-Richtmeyer equivalence theorem (see e.g. , , ) states that if the numerical scheme is stable then its global error is at least of the same order as its local truncation error. Recent work [17, 29, 18, 28, 5] highlights 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.
Schemes that have global errors that are higher order than predicted by the local truncation errors were devised by Kulikov  and by Weiner and colleagues  following the quasi-consistency theory first introduced by Skeel in 1978 . These methods have truncation errors that are of order but have global errors of order . In  Ditkowski and Gottlieb derived sufficient conditions for a class of general linear methods (GLMs) 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.
, that allow the computation of the the precise form of the first surviving term in the global error (the vector multiplying
). In these works, this form was leveraged for error estimation. However, the conditions in[18, 28, 19], were overly restrictive and made it difficult to find methods. 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 they showed how the solution can be post-processed to obtain accuracy of order . In this work, they produced implicit method and explicit methods that have favorable stability properties. The coefficients of those methods can be downloaded from . In  Ditkowski, Gottlieb, and Grant derived sufficient conditions under which two-derivative GLMs of a certain form can produce solutions of order and can be post-processed to obtain solutions of order . The coefficients of those methods can be downloaded from .
In many problems, there is a component that is stiff and easy to invert (e.g. linear) and a component that is non-stiff but difficult to invert (e.g. nonlinear). Such examples include advection-diffusion PDEs, where the diffusion term requires a small time-step when used with an explicit time-stepping method but is linear, and the advection term is nonlinear but allows a relatively large time-step when used with an explicit time-stepping method. In such cases it is convenient to distinguish these two components and write the resulting autonomous ODE as
(where we assume, as before, the required smoothness on and ). Here, is a stiff but ’simple to invert’ operator while is non-stiff but complicated (or costly) to invert and typically non-linear.
Implicit-explicit (IMEX) methods were developed to treat the component explicitly and the component implicitly, thus alleviating the linear stability constraint on the stiff component while treating the difficult to invert component explicitly. IMEX methods were first introduced by Crouziex in 1980  for evolving parabolic equations. For time dependent PDEs (notably (convection-diffusion equations) Ascher, Ruuth, and Wetton  introduced IMEX multi-step methods and Ascher, Ruuth and Spiteri  IMEX Runge–Kutta schemes. Implicit methods are often particularly desirable when applied to a linear component, in which case the order conditions may simplify: such methods were considered by Calvo, Frutos, and Novo in . In 2003, Kennedy and Carpenter  derived IMEX Runge–Kutta methods based on singly diagonally implicit Runge–Kutta (SDIRK) methods. This work introduced sophisticated IMEX methods with good accuracy and stability properties, as well as high quality embedded methods for error control and other features that make these methods usable in complicated applications. Recently, there has been more interest in IMEX methods with multiple steps and stages, including , . Optimized stability regions for IMEX methods with multiple steps and stages have also been of recent interest, as in  and .
In recent work, Schneider and colleagues [23, 24] produced super-convergent IMEX Peer methods (a Peer method is a GLM where each stage is of the same order). These methods satisfy error inhibiting conditions similar to those in  and so produce order although their truncation error is of order , so we shall refer to them here as IMEX-EIS schemes. In this work, we extend the error inhibiting with post-processing theory in  and  to IMEX methods. We proceed to devise error inhibiting IMEX methods of up to five steps and order six that have truncation errors or order but which can be post-processed to attain a solution of order . We refer to these as IMEX-EIS+ methods. We test these methods on a number of numerical examples to demonstrate their enhanced accuracy properties. This paper is structured as follows: in Section 2 we present our notation and some information about our IMEX methods which have multiple steps and stages. In Section 3.1 we provide sufficient conditions under which these methods with local truncation error of order will produce solutions of order , which can be post-processed to order . We provide the construction of one such post-processor in Section 3.2. In Section 4 we present some of our new methods and show their linear stability properties, and in Section 5 we show how these methods perform in practice on some numerical test cases.
2 Multistep multi-stage additive methods
In this work we consider the class of additive general linear methods (GLMs) of the form
where is a vector of length that contains the numerical solution at times for :
The functions and are defined as the component-wise function evaluation on the vector :
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 a sufficiently accurate method.
Similarly, we define the projection of the exact solution of the ODE (3) onto the temporal grid by:
with 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
Note that the additive form (4) 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. A special case of additive methods is that the method is explicit in and implicit in : such methods are called implicit-explicit or IMEX methods. In this work we are primarily interested in IMEX schemes so we require that in the form (4) is strictly lower triangular, and that contains some diagonal or super-diagonal elements. However, in this section and in Section (3.1) we do not assume this is the case, and so our error inhibiting with post-processing theory is true for all additive methods of the form (4), which is a much more general class than IMEX methods.
2.1 Truncation errors
A method of the form (4) has an approximation error at at time
We denote the vector of abscissas by , and the vector of ones by . Any terms of the form are to be understood component-wise . Note that this notation matches the one used in .
Alternatively, observing that we can write this as:
This notation is the same as used in .
As this must be true for all and these become
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: this 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 an observation that will be useful in the remainder of the paper.
Observation 1 Given the smoothness of and , and the assumption that , we observe that
where and .
This is simply due to the fact that is smooth enough and is small enough that we can expand:
where the error vector is , and we use the notation . Each term can be expanded as
so that, using the assumption that we have
The proof for is identical.
We now describe the growth of the error:
Given a zero-stable method of the form (4) which satisfies the order conditions
and where the functions and are smooth, the evolution of the error can be described by:
using Observation 1 we have
The term final term is very small: so we can neglect it. We then have an equation for the error which is essentially linear because the terms , , and are all constant scalars, which depend only on the solution of the ODE.
Now move the terms to the left side:
Noting that and we have the desired result
3 Constructing IMEX-EIS+ methods
Schneider, Lang, and Hundsdorfer  showed how to construct IMEX methods of the form (4) which satisfy the order conditions to order but are error inhibiting and so produce a solution of order . In this section we show that under additional conditions we can express the exact form of the next term in error and define an associated post-processor that allow us to recover order from a scheme that would otherwise be only th-order accurate. We note that also our main interest is in IMEX methods, and we will focus on these later in the paper, the theory presented in this section is general enough for all additive methods.
3.1 IMEX error inhibiting schemes that can be post-processed
In this section we consider an additive method of the form (4) 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 to order , so that the method will give us a numerical solution that has error that is guaranteed of order .
In  it was shown that if the truncation error vector lives in the null-space of the operator then the order of the error is of order . In the following theorem we establish additional conditions on the coefficient matrices , , , , , which allow us to determine precisely what the leading term of this error will look like and therefore remove it by post-processing.
Consider a zero-stable additive general linear method of the form
where is a rank one matrix that satisfies the consistency condition , the coefficient matrices , , , , , satisfy the order conditions
for , defined by (10), and the error inhibiting condition equationparentequation
|is satisfied (so that numerical solution produced by this method will have error ). If the conditions|
are 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 and are smooth assures us that so that we can expand the first term
As we showed in Lemma 3 in , we can assume that for some reasonably large time interval .
The discrete Duhamel principle (given as Lemma 5.1.1 in ) states that given an iterative process of this form where is a linear operator, we have
To analyze the error (22), 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 (19a) that states that .
The result in this theorem requires a closer look at the terms in this product
where we applied (19a) and used the observation that and .
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 Truncation error
We use the facts that: (1) the form of means that ; (2) , ; and (3) for any integer , to obtain
The last terms disappear because of (19b),