Two-derivative error inhibiting schemes with post-processing

12/09/2019 ∙ by Adi Ditkowski, et al. ∙ Oak Ridge National Laboratory Tel Aviv University University of Massachusetts Dartmouth 0

High order methods are often desired for the evolution of ordinary differential equations, in particular those arising from the semi-discretization of partial differential equations. In prior work in we investigated the interplay between the local truncation error and the global error to construct error inhibiting general linear methods (GLMs) that control the accumulation of the local truncation error over time. Furthermore we defined sufficient conditions that allow us to post-process the final solution and obtain a solution that is two orders of accuracy higher than expected from truncation error analysis alone. In this work we extend this theory to the class of two-derivative GLMs. We define sufficient conditions that control the growth of the error so that the solution is one order higher than expected from truncation error analysis, and furthermore define the construction of a simple post-processor that will extract an additional order of accuracy. Using these conditions as constraints, we develop an optimization code that enables us to find explicit two-derivative methods up to eighth order that have favorable stability regions, explicit strong stability preserving methods up to seventh order, and A-stable implicit methods up to fifth order. We numerically verify the order of convergence of a selection of these methods, and the total variation diminishing performance of some of the SSP methods. We confirm that the methods found perform as predicted by the theory developed herein.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The celebrated Lax-Richtmeyer equivalence theorem ([21], [12], [29]) 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 [35] 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 [14], 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 [3], use multiple stages as in a Runge–Kutta methods [3], or include higher derivatives as we do in Taylor series methods [3]. 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) [18] and by Weiner and colleagues [40] (for explicit two-step peer methods) followed along the lines of the quasi-consistency theory first introduced by Skeel in 1978 [31] to find methods that have solution of order although their truncation errors are of order . In [6] 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 [34]. In [19], [39], and [20]

, 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

[7] 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 [7] carries through easily to such methods, and by imposing the same EIS+ conditions as in [7] 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 [7] 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 [26] 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)

2.1 Essentials

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.

2.2 Preliminaries

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


where .


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


where .


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:



Recall that

and use this equation to subtract from


using Observations 1 and 2 we have

so that

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:


Given a zero-stable scheme (2), if the order conditions are satisfied for , and and are bounded, then there is some time interval [0,T] such that the error (given by (13) ) satisfies

The proof of this Lemma is given in [7], 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.

Theorem 1

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 [12]) 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 .

    Looking more closely at the terms that remain in this product and applying (14a) again

    where we used the observation that . The product is then

    we now use the two error inhibiting conditions (14b) and (14c) to obtain

  • 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

    we have

    Using the fact that in our case and that we obtain

    The first term and third terms disappear because of (14a) and (14b)

    The second term is eliminated by (14c)

    So that

All these parts together give us the result



We note that generalizing this approach from our work in [7] 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 [7]. 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 [7]. 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 [7], the only difference being the definition of the truncation error vector. The theoretical discussion of this post-processor is presented in [7] and we refer the reader to the presentation there. We review one such construction below:

  1. 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.

  2. Define the time vector of all the temporal grid points in the last computation steps:

  3. Correspondingly, concatenate the final