1 Introduction
Efficient high order numerical methods for evolving the solution of an ordinary differential equation (ODE) are frequently desired, especially for the timeevolution of PDEs that are semidiscretized in space. Higher order methods are desirable as they enable more accurate solutions for larger stepsizes. In this paper we consider numerical solvers for ordinary differential equations (ODEs)
(1)  
(where is sufficiently smooth). Any nonautonomous system of this form can be converted to an autonomous system (see [15]) so that without loss of generality, we restrict ourselves to the autonomous case:
(2)  
The forward Euler method is the simplest numerical method for evolving this problem forward in time with timesteps :
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 multistep 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 [27] states that any zerostable, 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 LaxRichtmeyer equivalence theorem (see e.g. [21], [12], [22]) 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 [17] and by Weiner and colleagues [29] following the quasiconsistency theory first introduced by Skeel in 1978 [25]. These methods have truncation errors that are of order but have global errors of order . In [5] 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.
Additional conditions were derived in [18, 28, 19]
, 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 [6] 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 postprocessed 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 [8]. In [7] Ditkowski, Gottlieb, and Grant derived sufficient conditions under which twoderivative GLMs of a certain form can produce solutions of order and can be postprocessed to obtain solutions of order . The coefficients of those methods can be downloaded from [9].In many problems, there is a component that is stiff and easy to invert (e.g. linear) and a component that is nonstiff but difficult to invert (e.g. nonlinear). Such examples include advectiondiffusion PDEs, where the diffusion term requires a small timestep when used with an explicit timestepping method but is linear, and the advection term is nonlinear but allows a relatively large timestep when used with an explicit timestepping method. In such cases it is convenient to distinguish these two components and write the resulting autonomous ODE as
(3)  
(where we assume, as before, the required smoothness on and ). Here, is a stiff but ’simple to invert’ operator while is nonstiff but complicated (or costly) to invert and typically nonlinear.
Implicitexplicit (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 [4] for evolving parabolic equations. For time dependent PDEs (notably (convectiondiffusion equations) Ascher, Ruuth, and Wetton [1] introduced IMEX multistep methods and Ascher, Ruuth and Spiteri [2] 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 [3]. In 2003, Kennedy and Carpenter [16] 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 [14], [26]. Optimized stability regions for IMEX methods with multiple steps and stages have also been of recent interest, as in [30] and [20].
In recent work, Schneider and colleagues [23, 24] produced superconvergent 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 [5] and so produce order although their truncation error is of order , so we shall refer to them here as IMEXEIS schemes. In this work, we extend the error inhibiting with postprocessing theory in [6] and [7] 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 postprocessed to attain a solution of order . We refer to these as IMEXEIS+ 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 postprocessed to order . We provide the construction of one such postprocessor 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 multistage additive methods
In this work we consider the class of additive general linear methods (GLMs) of the form
(4) 
where is a vector of length that contains the numerical solution at times for :
(5) 
The functions and are defined as the componentwise function evaluation on the vector :
(6) 
and
(7) 
For convenience, we select so that the first element in the vector approximates the solution at time , and the abscissas are nondecreasing . 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:
(8) 
with and are the componentwise 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
(9) 
Remark
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 implicitexplicit 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 superdiagonal elements. However, in this section and in Section (3.1) we do not assume this is the case, and so our error inhibiting with postprocessing 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
(10) 
where
(11) 
where equationparentequation
(12a)  
(12b)  
(12c) 
We denote the vector of abscissas by , and the vector of ones by . Any terms of the form are to be understood componentwise . Note that this notation matches the one used in [20].
Alternatively, observing that we can write this as:
where equationparentequation
(13a)  
(13b)  
(13c) 
This notation is the same as used in [26].
Regardless of whether we use the notation (12) or (13), for a method to have order truncation error of order it must satisfy the order conditions
(14) 
As this must be true for all and these become
or, equivalently,
We are only interested in zerostable 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.
2.2 Preliminaries
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
(15) 
and
(16) 
where and .
Proof
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:
Lemma
Given a zerostable 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:
Proof
Subtracting (10) from (4), we obtain
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 IMEXEIS+ methods
Schneider, Lang, and Hundsdorfer [23] 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 postprocessor that allow us to recover order from a scheme that would otherwise be only thorder 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 postprocessed
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 zerostable. 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 [23] it was shown that if the truncation error vector lives in the nullspace 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 postprocessing.
Theorem 1
Consider a zerostable additive general linear method of the form
(18) 
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
(19a)  
is satisfied (so that numerical solution produced by this method will have error ). If the conditions  
(19b)  
(19c)  
(19d) 
are satisfied, then error vector will have the more precise form:
(20) 
Proof
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 [6], we can assume that for some reasonably large time interval .
The discrete Duhamel principle (given as Lemma 5.1.1 in [12]) states that given an iterative process of this form where is a linear operator, we have
(22) 
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 have
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),
Comments
There are no comments yet.