# Two-derivative error inhibiting schemes with post-processing

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.

## Authors

• 5 publications
• 9 publications
• 8 publications
• ### IMEX error inhibiting schemes with post-processing

High order implicit-explicit (IMEX) methods are often desired when evolv...
12/19/2019 ∙ by Adi Ditkowski, et al. ∙ 0

• ### Explicit and implicit error inhibiting schemes with post-processing

Efficient high order numerical methods for evolving the solution of an o...
10/07/2019 ∙ by Adi Ditkowski, et al. ∙ 0

• ### High order unconditionally strong stability preserving multi-derivative implicit and IMEX Runge–Kutta methods with asymptotic preserving properties

In this work we present a class of high order unconditionally strong sta...
02/23/2021 ∙ by Sigal Gottlieb, et al. ∙ 0

• ### Exponentially fitted two-derivative DIRK methods for oscillatory differential equations

In this work, we construct and derive a new class of exponentially fitte...
12/27/2020 ∙ by Julius O. Ehigie, et al. ∙ 0

• ### Discrete Adjoint Implicit Peer Methods in Optimal Control

It is well known that in the first-discretize-then-optimize approach in ...
02/27/2020 ∙ by Jens Lang, et al. ∙ 0

• ### Predictable Feature Analysis

Every organism in an environment, whether biological, robotic or virtual...
11/11/2013 ∙ by Stefan Richthofer, et al. ∙ 0

• ### A development of Lagrange interpolation, Part I: Theory

In this work, we introduce the new class of functions which can use to s...
04/27/2019 ∙ by Mehdi Delkhosh, et al. ∙ 0

##### 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

 ut=F(t,u),t≥0 u(t0)=u0.

(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

 ut=F(u),t≥0 (1) u(t0)=u0.

The canonical numerical method for this problem is the forward Euler method

 vn+1=vn+ΔtF(vn),

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

 τn=ΔtLTEn=u(tn−1)+ΔtF(u(tn−1))−u(tn)≈O(Δt2)

and it produces a global error which is first order accurate

 en=vn−u(tn)≈O(Δt).

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

 Vn+1=DVn+ΔtAF(Vn)+ΔtRF(Vn+1)+Δt2^A˙F(Vn)+Δt2^R˙F(Vn+1), (2)

where is a vector of length that contains the numerical solution at times for :

 Vn=(v(tn+c1Δt)),v(tn+c2Δt),…,v(tn+csΔt))T. (3)

The function is defined as the component-wise function evaluation on the vector :

 F(Vn)=(F(v(tn+c1Δt)),F(v(tn+c2Δt)),…,F(v(tn+csΔt)))T, (4)

and, similarly, the function is

 ˙F(Vn)=(dFdt(v(tn+c1Δt)),dFdt(v(tn+c2Δt)),…,dFdt(v(tn+csΔt)))T. (5)

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:

 Un=(u(tn+c1Δt),u(tn+c2Δt),…,u(tn+csΔt))T, (6)

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

 En=Vn−Un. (7)

Methods of the form (2) have a local truncation error at time that is related to the approximation error

 ΔtLTEn=τn=[DUn−1+ΔtAF(Un−1)+ΔtRF(Un)]−Un (8)

where

 τn=∞∑j=0τnjΔtj=∞∑j=0τjΔtjdjudtj∣∣∣t=tn (9)

where the truncation error vectors have the form equationparentequation

 τ0 = (D−I)1 (10a) τj = 1(j−1)!(1jD(c−1)j+A(c−1)j−1+(j−1)^A(c−1)j−2 + Rcj−1+(j−1)^Rcj−2−1jcj)for j=1,2, ....

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

 τj=0j≤p.

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

 F(Un+En)=F(Un)+FnyEn+O(Δt)En, (11)

where .

###### Proof

This is simply due to the fact that is smooth enough that it can be expanded:

 F(Un+En) =F(Un)+⎛⎜ ⎜ ⎜ ⎜ ⎜⎝Fy(u(tn+c1Δt))en+c1Fy(u(tn+c2Δt))en+c2⋮Fy(u(tn+csΔt))en+cs⎞⎟ ⎟ ⎟ ⎟ ⎟⎠+O(∥En∥2), =F(Un)+⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝Fn+c1y0⋯00Fn+c2y⋯0⋮⋮⋱⋮00⋯Fn+csy⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠En+O(∥En∥2),

where the error vector is , and we use the notation . Each term can be expanded as

 Fn+cjy=Fy(u(tn+cjΔt))=Fy(u(tn))+cjΔtFyy(u(tn))+O(Δt2)

so that, exploiting the smoothness of we have

 F(Un+En) =F(Un)+(FnyI+O(Δt))En+O(∥En∥2),

which, due to the fact that becomes

 F(Un+En) =F(Un)+FnyEn+O(Δt)En.

Observation 2 Given the smoothness of , we observe that

 ˙F(Un+En)=˙F(Un)+˙FnyEn+O(Δt)En, (12)

where .

###### Proof

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:

###### Lemma

Given a zero-stable method of the form (2) which satisfies the order conditions

 τj=0forj=0,...,p

and where the functions and are smooth, the evolution of the error can be described by:

 (I−ΔtRFny+O(Δt2))En+1 = (D+ΔtAFny+O(Δt2))En+τn+1. (13)

###### Proof

Recall that

 τn+1=DUn+ΔtAF(Un)+Δt2^A˙F(Un)+ΔtRF(Un+1)+Δt2^R˙F(Un+1)−Un+1.

and use this equation to subtract from

 Vn+1=DVn+ΔtAF(Vn)+Δt2^A˙F(Vn)+ΔtRF(Vn+1)+Δt2^R˙F(Vn+1)

obtaining

 Vn+1−Un+1 = DVn−DUn+ΔtA(F(Vn)−F(Un)) +Δt2^A(˙F(Vn)−˙F(Un))+ΔtR(F(Vn+1)−F(Un+1)) +Δt2^R(˙F(Vn+1)−˙F(Un+1))+τn+1

using Observations 1 and 2 we have

 En+1 = D(Un+En)−DUn+ΔtA(F(Un+En)−F(Un)) +Δt2^A(˙F(Un+En)−˙F(Un))+ΔtR(F(Un+1+En+1)−F(Un+1)) +Δt2^R(˙F(Un+1+En+1)−˙F(Un+1))+τn+1 = DEn+Δt(AFny+O(Δt))En+Δt2(˙Fny^A+O(Δt))En +Δt(RFn+1y+O(Δt))En+1+Δt2(˙Fn+1y^R+O(Δt))En+1+τn+1,

so that

 En+1 = DEn+ΔtAFnyEn+ΔtRFn+1yEn+1+Δt2^A˙FnyEn+Δt2^R˙Fn+1yEn+1 +O(Δt2)En+O(Δt2)En+1+τn+1.

Note that and are all scalars. Now recall that since , and that terms of the form and are all bounded, so we have

 En+1 = DEn+ΔtAFnyEn+ΔtRFn+1yEn+1+O(Δt2)En+O(Δt2)En+1.

Now move the terms to the left side:

 (I−ΔtRFn+1y+O(Δt2))En+1 = DEn+ΔtAFnyEn+O(Δt2)+τn+1.

Noting that we have the desired result

 (I−ΔtRFny+O(Δt2))En+1=(D+ΔtAFny+O(Δt2))En+τn+1.

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:

###### Lemma

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

 ∥En∥≤O(Δtp)<<1.

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 p+2

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 p+1 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

 Vn+1=DVn+ΔtAF(Vn)+ΔtRF(Vn+1)+Δt2^A˙F(Vn)+Δt2^R˙F(Vn+1)

where is a rank one matrix that satisfies the consistency condition , and the coefficient matrices , , , , satisfy the order conditions

 τj=0forj=1,...,p,

where

 τj = 1(j−1)!(1jD(c−1)j+A(c−1)j−1+(j−1)^A(c−1)j−2 + Rcj−1+(j−1)^Rcj−2−1jcj)for j=1,2, ...,

if the error inhibiting condition equationparentequation

 Dτp+1=0 (14a) is satisfied, then the numerical solution produced by this method will have error En=O(Δtp+1). Furthermore, if the conditions Dτp+2=0 (14b) D(A+R)τp+1=0 (14c)

are also satisfied, then error vector will have the more precise form:

 En=Δtp+1τnp+1+O(Δtp+2). (15)

###### Proof

Using Lemma 1 we obtain the equation for the evolution of the error

 En+1 = (I−ΔtRFny+O(Δt2))−1[(D+ΔtAFny+O(Δt2))En+τn+1] (16)

the fact is smooth assures us that so that we can expand the first term to obtain

 En+1 = (I+ΔtRFny+O(Δt2))[(D+ΔtAFny+O(Δt2))En+τn+1] = (D+ΔtFny(RD+A)+O(Δt2))En +Δt(Δtpτn+1p+1+Δtp+1(τn+1p+2+FnyRτn+1p+1)+O(Δt2)) = QnEn+ΔtTne.

The discrete Duhamel principle (given as Lemma 5.1.1 in [12]) states that given an iterative process of the form

 En+1=QnEn+ΔtTne

where is a linear operator, we have

 En=n−1∏μ=0QμE0+Δtn−1∑ν=0(n−1∏μ=ν+1Qμ)Tνe. (17)

To analyze the error, we separate it into four parts:

 En=n−1∏μ=0QμE0I+ΔtTn−1eII+ΔtQn−1Tn−2eIII+Δtn−3∑ν=0(n−1∏μ=ν+1Qμ)TνeIV

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

 Tn−1e=(Δtpτnp+1+Δtp+1(τnp+2+Fn−1yRτnp+1)+O(Δtp+2))

so that this term contributes to the final time error the term

 ΔtTn−1e=Δtp+1τnp+1+O(Δtp+2).
• The term is a product of the operator

 Qn−1=(D+ΔtFn−1y(RD+A)+O(Δt2))

and the approximation error

 Tn−2e=(Δtpτn−1p+1+Δtp+1(τn−1p+2+Fn−2yRτn−1p+1)+O(Δtp+2)),

so we obtain

 Qn−1Tn−2e=ΔtpDτn−1p+1+O(Δtp+1)=O(Δtp+1),

due to the condition (14a) that states that .

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

 D(τn−1p+2+Fn−2yRτn−1p+1)+Fn−1y(RD+A)τn−1p+1+O(Δt) =D(τn−1p+2+Fn−2yRτn−1p+1)+Fn−1yAτn−1p+1+O(Δt) =D(τn−1p+2+Fn−1yRτn−1p+1)+Fn−1yAτn−1p+1+O(Δt) =Dτn−1p+2+Fn−1y(DR+A)τn−1p+1+O(Δt)

where we used the observation that . The product is then

 Qn−1Tn−2e=Δtp+1[Dτn−1p+2+Fn−1y(DR+A)τn−1p+1]+O(Δtp+2),

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

 Qn−1Tn−2e=O(Δtp+2).
• Finally we look at the rest of the sum and use the boundedness of the operator to observe

 ∥∥ ∥∥Δtn−3∑ν=0(n−1∏μ=ν+1Qμ)Tνe∥∥ ∥∥ = ∥∥ ∥∥Δtn−3∑ν=0(n−1∏μ=ν+3^Qμ)(^Qν+2^Qν+1Tνe)∥∥ ∥∥ ≤ Δtn−3∑ν=0∥∥ ∥∥n−1∏μ=ν+3Qμ∥∥ ∥∥∥∥Qν+2Qν+1Tνe∥∥ ≤ Δtn−3∑ν=0(1+cΔt)n−ν−3∥∥Qν+2Qν+1Tνe∥∥ = exp(ctn)−1cmaxν=0,...,n−3∥∥Qν+2Qν+1Tνe∥∥.

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

 Qν+2Qν+1=[D2+Δt(Fν+2y(RD+A)D+Fν+1yD(RD+A))+O(Δt2)]

and the approximation error

 Tνe=(Δtpτν+1p+1+Δtp+1(τν+1p+2+FνyRτν+1p+1)+O(Δtp+2))

we have

 Qν+2Qν+1Tνe = Δtp(D+ΔtFν+2y(RD+A)+ΔtFν+1yDR)Dτν+1p+1 + Δtp+1(Fν+1yD2R+Fν+1yDA)τν+1p+1+Δtp+1D2τν+1p+2+O(Δtp+2).

Using the fact that in our case and that we obtain

 Qν+2Qν+1Tνe = Δtp(D+ΔtFν+2y(RD+A)+ΔtFν+1yDR)Dτν+1p+1 + Δtp+1Fν+1yD(R+A)τν+1p+1+Δtp+1Dτν+1p+2+O(Δtp+2).

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

 Dτp+1=0andDτp+2=0.

The second term is eliminated by (14c)

 D(R+A)τp+1=0.

So that

 Qν+2Qν+1Tνe=O(Δtp+2).

All these parts together give us the result

 En = n−1∏μ=0QμE0I+ΔtTn−1eII+ΔtQn−1Tn−2eIII+Δtn−1∑ν=0(n−3∏μ=ν+1Qμ)TνeIV (18) = 0+Δtp+1τnp+1II+O(Δtp+2)III+O(Δtp+2)IV = Δtp+1τnp+1+O(Δtp+2).

###### Remark

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 p+2 order

At every time-step the error vector has the particular form (15)

 En=Δtp+1τnp+1+O(Δtp+2),

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

 ms≥p+3

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:

 ~t=(tn−m+1+c1Δt,...,tn−m+1+csΔt,...,tn+c1Δt,...,tn+csΔt)T
3. Correspondingly, concatenate the final