A Strong Stability Preserving Analysis for Explicit Multistage Two‑Derivative Time‑Stepping Schemes Based on Taylor Series Conditions

Zachary Grant·Sigal Gottlieb·David C.Seal

Abstract

Keywords Strong stability preserving·Taylor series·Hyperbolic conservation laws·Two derivative Runge Kutta

1 Introduction

The solution to a hyperbolic conservation law

may develop sharp gradients or discontinuities,which results in significant challenges to the numerical simulation of such problems.To ensure that it can handle the presence of a discontinuity,a spatial discretization is carefully designed to satisfy some nonlinear stability properties,often in a non-inner-product sense,e.g.,total variation diminishing,maximum norm preserving,or positivity-preserving properties.The development of high-order spatial discretizations that can handle discontinuities is a major research area [4,6,15,28,30,35,46–49,54].

When the partial differential equation (1) is semi-discretized we obtain the ordinary differential equation (ODE)

whereuis a vector of approximations toU.This formulation is often referred to as a method of lines (MOL) formulation,and has the advantage of decoupling the spatial and time discretizations.The spatial discretizations designed to handle discontinuities ensure that when the semi-discretized equation (2) is evolved using a forward Euler method

whereunis a discrete approximation toUat timetn,and the numerical solution satisfies the desired strong stability property,such as total variation stability or positivity.If the desired nonlinear stability property such as a norm,semi-norm,or convex functional,is represented by ‖⋅‖ and the spatial discretization satisfies the monotonicity property

under the time-step restriction

In practice,in place of the first-order time discretization (3),we typically require a higher order time integrator,that preserves the strong stability property

perhaps under a modified time-step restriction.For this purpose,time discretizations with good linear stability properties or even with nonlinear inner-product stability properties are not sufficient.Strong stability preserving (SSP) time discretizations were developed to address this need.SSP multistep and Runge-Kutta methods satisfy the strong stability property (6) for any functionF,any initial condition,and any convex functional ‖⋅‖ under some time-step restriction,provided only that (4) is satisfied.

Recently,there has been interest in exploring the SSP properties of multi-derivative Runge-Kutta methods,also known as multistage multi-derivative methods.Multi-derivative Runge-Kutta (MDRK) methods were first considered in [3,21,22,31,33,34,40,41,45,51],and later explored for use with partial differential equations (PDEs) [9,29,36,39,50].These methods have a form similar to Runge-Kutta methods but use an additional derivative=utt≈Uttto allow for higher order.The SSP properties of these methods were discussed in [5,32].In [5],a method is defined as SSP-SD ifit satisfies the strong stability property (6) for any functionF,any initial condition,and any convex functional ‖⋅‖ under some time-step restriction,provided that (4) is satisfied for any Δt≤ΔtFEand an additional condition of the form

is satisfied for any Δt≤ΔtFE.These conditions allow us to find a wide variety of time discretizations,called SSP-SD time-stepping methods,but they limit the type of spatial discretization that can be used in this context.

In this paper,we present a different approach to the SSP analysis,which is more along the lines of the idea in [32].For this analysis,we use,as before,the forward Euler base condition (4),but add to it a Taylor series condition of the form

that holds for any Δt≤KΔtFE.Compared to those studied in [5],this pair of base conditions allows for more flexibility in the choice of spatial discretizations (such as the methods that satisfy a Taylor series condition in [7,10,37,39]),at the cost of more limited variety of time discretizations.We call the methods that preserve the strong stability properties of these base conditions strong stability preserving Taylor series (SSP-TS) methods.The goal of this paper is to study a class of methods that is suitable for use with existing spatial discretizations,and present families of such SSP-TS methods that are optimized for the relationship between the forward Euler time step ΔtFEand the Taylor series time stepKΔtFE.

In the following subsections we describe SSP Runge-Kutta time discretizations and present explicit multistage two-derivative methods.We then motivate the need for methods that preserve the nonlinear stability properties of the forward Euler and Taylor series base conditions.In Sect.2 we formulate the SSP optimization problem for finding explicit twoderivative methods which can be written as the convex combination of forward Euler and Taylor series steps with the largest allowable time step,which we will later use to find optimized methods.In Sect.2.1 we explore the relationship between SSP-SD methods and SSP-TS methods.In Sect.2.2 we prove that there are order barriers associated with explicit two-derivative methods that preserve the properties of forward Euler and Taylor series steps with a positive time step.In Sect.3 we present the SSP coefficients of the optimized methods we obtain.The methods themselves can be downloaded from our github repository [14].In Sect.4 we demonstrate how these methods perform on specially selected test cases,and in Sect.5 we present our conclusions.

1.1 SSP Methods

It is well known [13,42] that some multistep and Runge-Kutta methods can be decomposed into convex combinations of forward Euler steps,so that any convex functional property satisfied by (4) will be preserved by these higher order time discretizations.If we rewrite thes-stage explicit Runge-Kutta method in the Shu-Osher form [43],it is clear that if all the coefficientsandare non-negative,and providedis zero only ifits correspondingis zero,then each stage can be written as a convex combination of forward Euler steps of the form (3),and be bounded by

If a method can be decomposed into such a convex combination of (3),with a positive value of C>0 then the method is called strong stability preserving (SSP),and the value C is called the SSP coefficient.SSP methods guarantee the strong stability properties of any spatial discretization,provided,only,that these properties are satisfied when using the forward Euler method.The convex combination approach guarantees that the intermediate stages in a Runge-Kutta method satisfy the desired strong stability property as well.The convex combination approach clearly provides a sufficient condition for preservation of strong stability.Moreover,it has also been shown that this condition is necessary [11,12,16,17].

Second- and third-order explicit Runge-Kutta methods [43] and later fourth-order methods [23,44] were found that admit such a convex combination decomposition with C>0 .However,it has been proven that explicit Runge-Kutta methods with positive SSP coefficient cannot be more than fourth-order accurate [27,38].

The time-step restriction (8) is comprised of two distinct factors: (1) the term ΔtFEthat is a property of the spatial discretization,and (2) the SSP coefficient C that is a property of the time discretization.Research on SSP time-stepping methods for hyperbolic PDEs has primarily focused on finding high-order time discretizations with the largest allowable time step Δt≤CΔtFEby maximizing the SSP coefficient C of the method.

High-order methods can also be obtained by adding more steps (e.g.,linear multistep methods) or more derivatives (Taylor series methods).Multistep methods that are SSP have been found [13],and explicit multistep SSP methods exist of very high orderp>4 ,but have severely restricted SSP coefficients [13].These approaches can be combined with Runge-Kutta methods to obtain methods with multiple steps,and stages.Explicit multistep multistage methods that are SSP and have orderp>4 have been developed as well [1,24].

1.2 Explicit Multistage Two‑Derivative Methods

Another way to obtain higher order methods is to use higher derivatives combined with the Runge-Kutta approach.An explicit multistage two-derivative time integrator is given by:

wherey(1)=un.

The coefficients can be put into matrix-vector form,where

As in our prior work [5],we focus on using explicit multistage two-derivative methods as time integrators for evolving hyperbolic PDEs.For our purposes,the operatorFis obtained by a spatial discretization of the termUt=−f(U)xto obtain the systemut=F(u) .Instead of computing the second-derivative termdirectly from the definition of the spatial discretizationF,we approximateby employing the Cauchy-Kovalevskaya procedure which uses the PDE (1) to replace the time derivatives by the spatial derivatives,and discretize these in space.

If the termF(u) is computed using a conservative spatial discretizationDxapplied to the flux:

then we approximate the second derivative

where a (potentially different) spatial differentiation operatoris used.Although these two approaches are different,the differences between them are of high order in space,so that in practice,as long as the spatial errors are smaller than the temporal errors,we see the correct order of accuracy in time,as shown in [5].

1.3 Motivation for the New Base Conditions for SSP Analysis

In [5] we considered explicit multistage two-derivative methods and developed sufficient conditions for a type of strong stability preservation for these methods.We showed that explicit SSPSD methods within this class can break this well known order barrier for explicit Runge-Kutta methods.In that work we considered two-derivative methods that preserve the strong stability property satisfied by a functionFunder a convex functional ‖⋅‖ ,provided that the conditions:

and

because it is more general.If the forward Euler (12) and second-derivative (13) conditions are both satisfied,then the Taylor series condition (14) will be satisfied as well.Thus,a spatial discretization that satisfies (12) and (13) will also satisfy (14),so that the SSP-SD concept in [5] allows for the most general time discretizations.Furthermore,some methods ofinterest and importance in the literature cannot be written using a Taylor series decomposition,most notably the unique two-stage fourth-order method1Note that here we use to indicate that these methods are designed for the exact time derivative of F.However,in practice we use the approximation as explained above.

which appears commonly in the literature on this subject [9,29,36].For these reasons,it made sense to first consider the SSP-SD property which relies on the pair of base conditions (12) and (13).

However,as we will see in the example below,there are spatial discretizations for which the second-derivative condition (13) is not satisfied but the forward Euler condition (12) and the Taylor series condition (14) are both satisfied.In such cases,the SSP-SD methods derived in [5] may not preserve the desired strong stability properties.The existence of such spatial discretizations is the main motivation for the current work,in which we re-examine the strong stability properties of the explicit two-derivative multistage method (9) using the base conditions (12) and (14).Methods that preserve the strong stability properties of (12) and (14) are called,herein,SSP-TS methods.The SSP-TS approach increases our flexibility in the choice of spatial discretization over the SSP-SD approach.Of course,this enhanced flexibility in the choice of spatial discretization is expected to result in limitations on the time discretization (e.g.,the two-stage fourth-order method is SSP-SD but not SSP-TS).

To illustrate the need for time discretizations that preserve the strong stability properties of spatial discretizations that satisfy (12) and (14),but not (13),consider the one-way wave equation

(heref(U)=U) whereFis defined by the first-order upwind method

When solving the PDE,we compute the operatorby simply applying the differentiation operator twice (note that(U)=1)

We note that when computed this way,the spatial discretizationFcoupled with the forward Euler satisfies the total variation diminishing (TVD) condition:

while the Taylor series term usingFandsatisfies the TVD

In other words,these spatial discretizations satisfy the conditions (12) and (14) withK=1 ,in the total variation semi-norm.However,(13) is not satisfied,so the methods derived in [5] cannot be used.Our goal in the current work is to develop time discretizations that will preserve the desired strong stability properties (e.g.,the total variation diminishing property) when using spatial discretizations such as the upwind approximation (17) that satisfy (12) and (14) but not (13).

Remark 1This simple first-order motivating example is chosen because these spatial discretizations are provably TVD and allow us to see clearly why the Taylor series base condition (14) is needed.In practice,we use higher order spatial discretizations such as WENO that do not have a theoretical guarantee of TVD,but perform well in practice.Such methods are considered in Examples 2 and 4 in the numerical tests,and provide us with similar results.

In this work we develop explicit two-derivative multistage SSP-TS methods of the form (9) that preserve the convex functional properties of forward Euler and Taylor series terms.When the spatial discretizationsFandthat satisfy (12) and (14) are coupled with such a timestepping method,the strong stability condition

will be preserved,perhaps under a different time-step condition

If a method can be decomposed in such a way,withwe say that it is SSP-TS.In the next section,we define an optimization problem that will allow us to find SSP-TS methods of the form (9) with the largest possible SSP coefficient

2 SSP Explicit Two‑Derivative Runge‑Kutta Methods

We consider the system of ODEs

resulting from a semi-discretization of the hyperbolic conservation law (1) such thatFsatisfies the forward Euler (first derivative) condition (12)for the desired stability property indicated by the convex functional ‖⋅‖.

The methods we are interested in also require an appropriate approximation to the second derivative in time

We assume in this work thatFandsatisfy an additional condition of the form (14)

in the same convex functional ‖⋅‖ ,whereKis a scaling factor that compares the stability condition of the Taylor series term to that of the forward Euler term.

We wish to show that given conditions (12) and (14),the multi-derivative method (9) satisfies the desired monotonicity condition under a given time step.This is easier if we re-write the method (9) in an equivalent matrix-vector form

We can now easily establish sufficient conditions for an explicit method of the form (22) to be SSP:

Theorem 1Given spatial discretizationsFand(12)and(14),an explicit two-derivative multistage method of the form(22)preserves the strong stability property‖un+1‖≤‖un‖under the time-step restrictionΔt≤rΔtFEifit satisfies the conditions

for somer>0 .In the above conditions,the inequalities are understood component-wise.

ProofWe begin with the method

where

If the elements ofP,Q,andare all non-negative,and ifthenis a convex combination of strongly stable terms

and so is also strongly stable under the time-step restrictions Δt≤rΔtFEandIn such cases,the optimal time step is given by the minimum of the two.In the cases we encounter here,this minimum occurs when these two values are set equal,so we requireConditions (23a)–(23c) now ensure thatP≥0 ,Q≥0 ,and≥0 componentwise forand so the method preserves the strong stability condition ‖un+1‖≤‖un‖ under the time-step restriction Δt≤rΔtFE.Note that if this fact holds for a given value ofr>0 then it also holds for all smaller positive values.

Definition 1A method that satisfies the conditions in Theorem 1 for valuesr∈(0,rmax] is called a Strong Stability Preserving Taylor Series (SSP-TS) method with an associated SSP coefficient

Remark 2Theorem 1 gives us the conditions for the method (22) to be SSP-TS for any time stepWe note,however,that while the corresponding conditions for Runge-Kutta methods have been shown to be necessary as well as sufficient,for the multiderivative methods we only show that these conditions are sufficient.This is a consequence of the fact that we define this notion of SSP based on the conditions (12) and (14),but if a spatial discretization also satisfies a different condition (for example,(13)) many other methods of the form (22) also give strong stability preserving results.Notable among these is the two-derivative two-stage fourth-order method (15) which is SSP-SD but not SSPTS.This means that solutions of (15) can be shown to satisfy the strong stability property ‖un+1‖≤‖un‖ for positive time steps,for the appropriate spatial discretizations,even though the conditions in Theorem 1 are not satisfied.

This result allows us to formulate the search for optimal SSP-TS methods as an optimization problem,as in [5,13,23,25].

Find the coefficient matrices S and

that maximize the value of=maxr

such that the relevant order conditions (summarized in Appendix 1)

and the SSP conditions (23a)–(23c) are all satisfied.

However,before we present the optimal methods in Sect.3,we present the theoretical results on the allowable order of multistage multi-derivative SSP-TS methods.

2.1 SSP Results for Explicit Two‑Derivative Runge‑Kutta Methods

In this paper,we consider explicit SSP-TS two-derivative multistage methods that can be decomposed into a convex combination of (12) and (14),and thus preserve their strong stability properties.In our previous work [5] we studied SSP-SD methods of the form (9) that can be written as convex combinations of (12) and (13).The following lemma explains the relationship between these two notions of strong stability.

Lemma 1Any explicit method of the form(9)that can be written as a convex combination of the forward Euler formula(12)and the Taylor series formula(14)can also be written as a convex combination of the forward Euler formula(12)and the second-derivative formula(13).

ProofWe can easily see that any Taylor series step can be rewritten as a convex combination of the forward Euler formula (12) and the second-derivative formula (13):

for any 01 .Clearly then,if a method can be decomposed into a convex combination of (12) and (14),and in turn (14) can be decomposed into a convex combination of (12) and (13),then the method itself can be written as a convex combination of (12) and (13).

This result recognizes that the SSP-TS methods we study in this paper are a subset of the SSP-SD methods in [5].This allows us to use results about SSP-SD methods when studying the properties of SSP-TS methods.

The following lemma establishes the Shu-Osher form of an SSP-SD method of the form (9).This form allows us to directly observe the convex combination of steps of the form (12) and (13),and thus easily identify the SSP coefficient CSD.

Lemma 2If an explicit method of the form(9) written in theShu-Osher form

has the properties that

(i)all the coefficients are non-negative,

(ii)=0whenever=0,

(iii)=0whenever=0,

then this method preserves the strong stability properties of(12)and(13) (i.e.,is SSP-SD)forΔt≤CSDΔtFEwith

ProofFor each stage we have

assuming only that for each one of theseThe result immediately follows from the fact that for eachiwe havefor consistency.

When a method is written in the block Butcher form (22),we can decompose it into a canonical Shu-Osher form,

This allows us to define an SSP-SD method directly from the Butcher coefficients.

Definition 2Given spatial discretizationsFandthat satisfy (12) and (13),an explicit two-derivative multistage method of the form (22) is called a Strong Stability Preserving Second Derivative (SSP-SD) method with and associated SSP coefficientifit satisfies the conditions

for allr=(0,rmax] andIn the above conditions,the inequalities are understood component-wise.

The relationship between the coefficients in (9) and (24) allows us to conclude that the matricesSandmust contain only non-negative coefficients.

Lemma 3If an explicit method of the form(9)can be converted to the Shu-Osher form(24)with all non-negative coefficientsfor alli,j,then the coefficientsmust be all non-negative as well.

ProofThe transformation between (9) and (24) is given byand,recursively,

Clearly then,

and

From there we proceed recursively: given that all≥0 and≥0 for alli,j,and thatakj≥0 and≥0 for all 1 ≤j

Now given≥0 and≥0 for alli,jand≥0 and≥0 for all 1 ≤j

We wish to study only those methods for which the Butcher form (9) is unique.To do so,we follow Higueras [18] in extending the reducibility definition of Dahlquist and Jeltsch [19].Other notions of reducibility exist,but for our purposes it is sufficient to define irreducibility as follows:

Definition 3A two-derivative multistage method of the form (9) isDJ-reducibleif there exist setsT1andT2such thatT1≠∅ ,T1∩T2=∅ ,T1∪T2=[1,2,…,s] ,and

We say a method isirreducibleifit is not DJ-reducible.

Lemma 4An irreducible explicit SSP-SD method of the form(22)must satisfy the(component-wise)condition

Proof An SSP-SD method in the block form (22),must satisfy conditions (25a)–(25c) forThe non-negativity of (25b) and (25c) requires their sum to be non-negative as well,

Note that these matrices commute,so we have

Recalling the definition of the matrixS,we have

Now,we can expand the inverse as

Because the positivity must hold for arbitrarily smallr<<1 and<<1 we can stop our expansion after the linear term,and require

which is

Now we are ready to address the proof.Assume thatj=Jis the largest value for which we havebJ==0,

Clearly,then we have

Since the method is explicit,the matricesAandare lower triangular (i.e.,aiJ==0 fori≤J),so this condition becomes

By the assumption above,we have (bi+)>0 fori>J,andr>0,0 .Clearly,then,for (27) to hold we must require

which,together withbJ==0 ,makes the method DJ-reducible.Thus we have a contradiction.

We note that this same result,in the context of additive Runge-Kutta methods,is due to Higueras [18].

2.2 Order Barriers

Explicit SSP Runge-Kutta methods with>0 are known to have an order barrier of four,while the implicit methods have a barrier of six [13].This follows from the fact that the orderpofirreducible methods with non-negative coefficients depends on the stage orderqsuch that

For explicit Runge-Kutta methods the first stage is a forward Euler step,soq=1 and thusp≤4 ,whereas for implicit Runge-Kutta methods the first stage is at most of order two,so thatq=2 and thusp≤6.

For two-derivative multistage SSP-TS methods,we find that similar results hold.A stage order ofq=2 is possible for explicit two-derivative methods (unlike explicit Runge-Kutta methods) because the first stage can be second order,i.e.,a Taylor series method.However,since the first stage can be no greater than second order we have a bound on the stage orderq≤2 ,which results in an order barrier ofp≤6 for these methods.In the following results we establish these order barriers.

Lemma 5Given an irreducible SSP-TS method of the form(9),ifbj=0,then the corresponding=0.

ProofIn any SSP-TS method the appearance of a second-derivative termcan only happen as part of a Taylor series term.This tells us thatmust be accompanied by the correspondingF,meaning that whenever we have a non-zeroorterm,the correspondingaijorbjterm must be non-zero.

Lemma 6Any irreducible explicit SSP-TS method of the form(9)must satisfy the(component-wise)condition

ProofAny irreducible method (9) that can be written as a convex combination of (12) and (14) can also be written as a convex combination of (12) and (13),according to Lemma 1.Applying Lemma 4 we obtain the conditionb+0 ,component-wise.Now,Lemma 5 tells us that if any componentbj=0 then its corresponding=0 ,so thatbj+0 for eachjimplies thatbj >0 for eachj.

Theorem 2Any irreducible explicit SSP-TS method of the form(9)with orderp≥5must satisfy the stage orderq=2condition

where the termc2is a component-wise squaring.

ProofA method of orderp≥5 must satisfy the 17 order conditions presented in the Appendix 1.Three of those necessary conditions are2In this work we use ⊙ to denote component-wise multiplication.

From this,we find that the following linear combination of these equations gives

(once again,the squaring here is component-wise).Given the strict component-wise positivity of the vectorbaccording to Lemma 6 and the non-negativity of,this condition becomes

Theorem 3Any irreducible explicit SSP-TS method of the form(9)cannot have orderp=7.

ProofThis proofis similar to the proof of Theorem 2.The complete list of additional order conditions for seventh order is lengthy and beyond the scope of this work.However,only three of these conditions are needed for this proof.These are

Combining these three equations we have

From this we see that any seventh order method of the form (9) which admits a decomposition of a convex combination of (12) and (14),must satisfy the stage orderq=3 condition

However,as noted above,the first stage of the explicit two-derivative multistage method (9) has the form

which can be at most of second order.This means that the stage order of explicit twoderivative multistage methods can be at mostq=2 ,and so the=0 condition cannot be satisfied.Thus,the result of the theorem follows.

Note that the order barriers do not hold for SSP-SD methods,because SSP-SD methods do not require that all components of the vectorbmust be strictly positive.

3 Optimized SSP Taylor Series Methods

In Sect.2 we formulated the search for optimal SSP two-derivative methods as:

Find the coefficient matricesSandthat maximize the value of=maxr

such that the relevant order conditions and the SSP conditions (23a)–(23b) are all satisfied.

To accomplish this,we develop and use a matlab optimization code [14] (similar to Ketcheson’s code [26]) for finding optimal two-derivative multistage methods that preserve the SSP properties (12) and (14).The SSP coefficients of the optimized SSP explicit multistage two-derivative methods of order up top=6 (for different values ofK) are presented in this section.

We considered three types of methods:

• (M1) Methods that have the general form (9) with no simplifications.

• (M2) Methods that are constrained to satisfy the stage order two (q=2 ) requirement (28),

• (M3) Methods that satisfy the stage order two (q=2 ) (28) requirement and require only(un) ,so they have only one second-derivative evaluation.This is equivalent to requiring that all values inand,except those on the first column of the matrix and the first element of the vector,be zero.

We refer to the methods by type,number of stages,order of accuracy,and value ofK.For example,an SSP-TS method of type (M1) withs=5 andp=4 ,optimized for the value ofK=1.5 would be referred to as SSP-TS M1(5,4,1.5) or as SSP-TS M1(5,4,1.5).For comparison,we refer to methods from [5] that are SSP in the sense that they preserve the properties of the spatial discretization coupled with (12) and (13) as SSP-SD MDRK(s,p,K) methods.

3.1 Fourth‑Order Methods

Using the optimization approach described above,we find fourth-order methods withs=3,4,5 stages for a range ofK=0.1,…,2.0 .In Fig.1 we show the SSP coefficients of methods of SSP-TS methods of type (M1) and (M2) withs=3,4,5 (in blue,red,green) plotted against the value ofK.The open stars indicate methods of type (M1) while the filled circles are methods of type (M2).Filled stars are (M1) markers overlaid with (M2) markers indicating close if not equal SSP coefficients.

Fig.1 The SSP-TS coefficient CTS (on the y-axis) of fourth-order SSP-TS M1 and M2 methods with s=3,4,5 stages plotted against the value of K (on the x-axis).The open stars indicate methods of type (M1) while the filled circles are methods of type (M2).Filled stars are (M1) markers overlaid with (M2) markers indicating close if not equal SSP coefficients

Three Stage MethodsThree-stage SSP-TS methods with fourth-order accuracy exist,and all these have stage order two (q=2 ),so they are all of type (M2).Figure 1 shows the SSP coefficients of these methods in blue.The (M3) methods have an SSP coefficient

For the case whereK≥1 we obtain the following optimal (M3) scheme with an SSP coefficient=1:

WhenK≤1 we have to modify the coefficients accordingly to obtain the maximal value of CTSas defined above.Here we provide the non-zero coefficients for this family of M3(3,4,K) as a function ofK:

In Table 1 we compare the SSP coefficient of three-stage fourth-order SSP-TS methods of type (M2) and (M3) for a selection of values ofK.Clearly,the (M3) methods have a much smaller SSP coefficient than the (M2) methods.However,a better measure of efficiency is the effective SSP coefficient computed by normalizing for the number of function evaluations required,which is 2sfor the (M2) methods,ands+1 for the (M3) methods.If we consider the effective SSP coefficient,we find that while the (M2) methods are more efficient for the larger values ofK,for smaller values ofKthe (M3) methods are more efficient.

Four Stage SSP-TS MethodsWhile four-stage fourth-order explicit SSP Runge-Kutta methods do not exist,four-stage fourth-order SSP-TS explicit two-derivative Runge-Kutta methods do.Four-stage fourth-order methods do not necessarily satisfy the stage order two (q=2 ) condition.These methods have a more nuanced behavior: for very smallK <0.2 ,the optimized SSP methods have stage orderq=1 .For 0.2

Table1 SSP-TS coefficients of three-stage fourth-order SSP-TS methods

Table2 SSP-TS coefficients of four-stage fourth-order SSP-TS methods

As seen in Table 2,the methods with the special structure (M3) have smaller SSP coefficients.But when we look at the effective SSP-TS coefficient we notice that,once again,for smallerKthey are more efficient.Table 2 shows that the (M3) methods are more efficient whenK≤1.5 ,and remain competitive for larger values ofK.

It is interesting to consider the limiting case,SSP-TS M2(4,4,∞ ),in which the Taylor series formula is unconditionally stable (i.e.,K=∞ ).This provides us with an upper bound of the SSP coefficient for this class of methods by ignoring any time-step constraint coming from condition (14).A four-stage fourth-order method that is optimal forK=∞ is

This method has an SSP coefficient=4 ,with an effective SSP coefficientThis method also has stage orderq=2 .This method is not intended to be useful in the SSP context but gives us an idea of the limiting behavior: i.e.,what the best possible value ofcould be if the Taylor series condition had no constraint (K=∞) .We observe in Table 2 that the SSP coefficient of the M2(4,4,K) method is within 10% of this limitingfor values ofK=2.

Table3 SSP-TS coefficients of five-stage fourth-order methods

Five-Stage MethodsThe optimized five-stage fourth-order methods have stage orderq=2 for the values of 0.5 ≤K≤7 ,and otherwise have stage orderq=1 .The SSP coefficients of these methods are shown in the green line in Fig.1,and the SSP and effective SSP coefficients for all three types of methods are compared in Table 3.We observe that these methods have higher effective SSP coefficients than the corresponding four-stage methods.

3.2 Fifth‑Order SSP‑TS Methods

While fifth-order explicit SSP Runge-Kutta methods do not exist,the addition of a second derivative which satisfies the Taylor Series condition allows us to find explicit SSPTS methods of fifth order.For fifth order,we have the result (in Sect.2.2 above) that all methods must satisfy the stage orderq=2 condition,so we consider only (M2) and (M3) methods.In Fig.2 we show the SSP-TS coefficients of M2(s,5,K) methods fors=4,5,6.

Four-Stage MethodsFour-stage fifth-order methods exist,and their SSP-TS coefficients are shown in blue in Fig.2.We were unable to find M3(4,5,K) methods,possibly due to the paucity of available coefficients for this form.

Five-Stage MethodsThe SSP coefficient of the five-stage M2 methods can be seen in red in Fig.2.We observe that the SSP coefficient of the M2(5,5,K) methods plateaus with respect toK.As shown in Table 4,methods with the form (M3) have a significantly smaller SSP coefficient than that of (M2).However,the effective SSP coefficient is more informative here,and we see that the (M3) methods are more efficient for small values ofK≤0.5 ,but not for larger values.

Fig.2 SSP-TS coefficients (on the y-axis) for M2(4,5,K) and M2(5,5,K) and M2(6,5,K),where K is on the x-axis

Table4 SSP-TS coefficients and effective SSP-TS coefficients of fifth-order methods

Six-Stage MethodsThe SSP coefficient of the six-stage M2 methods can be seen in green in Fig.2.In Table 4 we compare the SSP coefficients and effective SSP coefficients of (M2) and (M3) methods.As in the case above,the methods with the form (M3) have a significantly smaller SSP coefficient than that of (M2),and the SSP coefficient of the (M3) methods plateaus with respect toK.However,the effective SSP coefficient shows that the (M3) methods are more efficient for small values ofK≤0.7 ,but not for larger values.

3.3 Sixth‑Order SSP‑TS Methods

As shown in Sect.2.2,the highest order of accuracy this class of methods can obtain isp=6 ,and these methods must satisfy (28).We find sixth-order methods withs=5,6,7 stages of type (M2).As to methods with the special structure M3,we are unable to find methods withs≤p,but we find M3(7,6,K) methods and M3(8,6,K) methods.In the first six rows of Table 5 we compare the SSP-TS and effective SSP-TS coefficients of the (M2) methods withs=5,6,7 stages.In the last four rows of Table 5 we compare the SSP coefficients and effective SSP coefficients for sixth-order methods withs=7,8 stages.Figure 3 shows the SSP-TS coefficients of the optimized (M3) methods for seven and eight stages,which clearly plateau with respect toK(as can be seen in the tables as well).For the sixth-order methods,it is clear that M3(8,6,K) methods are most efficient for all values ofK.

Table5 SSP-TS coefficients and effective SSP-TS coefficients of sixth-order SSP-TS methods

Fig.3 SSP-TS coefficients (on the y-axis) for M3(7,6,K) and M3(8,6,K),where K is on the x-axis

3.4 Comparison with Existing Methods

First,we wish to compare the methods in this work to those in our prior work [5].If a spatial discretization satisfies the forward Euler condition (12) and the second-derivative condition (13) it will also satisfy the Taylor series condition (14),with

In this case,it is preferable to use the SSP-SD MDRK methods in [5].However,in the case that the second-derivative condition (13) is not satisfied for any value of>0 ,or if the Taylor series condition is independently satisfied with a largerKthan would be established from the two conditions,i.e.,then it may be preferable to use one of the SSP-TS methods derived in this work.

Next,we wish to compare the methods in this work to those in [32],which was the first paper to consider an SSP property based on the forward Euler and Taylor series base conditions.The approach used in our work is somewhat similar to that in [32] where the authors consider building time integration schemes which can be composed as convex combinations of forward Euler and Taylor series time steps,where they aim to find methods which are optimized for the largest SSP coefficients.However,there are several differences between our approach and the one of [32],which results in the fact that in this paper we are able to find more methods,of higher order,and with better SSP coefficients.In addition,in the present work we find and prove an order barrier for SSP-TS methods.

The first difference between our approach and the approach in [32] is that we allow computations ofof the intermediate values,rather than only(un) .Another way of saying this is that we consider SSSP-TS methods that are not of type M3,while the methods considered in [32] are all of type M3.In some cases,when we restrict our search to M3 methods andK=1 ,we find methods with the same SSP coefficient as in [32].For example,HBT34 matches our SSP-TS M3(3,4,1) method with an SSP coefficient of CTS =1 ,HBT44 matches our SSP-TS M3(4,4,1) method withHBT54 matches our SSP-TS M3(5,4,1) method with CTS =2.441 ,and HBT55 matches our SSP-TS M3(5,5,1) method with an SSP coefficient of CTS =1.062 .While methods of type M3 have their advantages,they are sometimes sub-optimal in terms of efficiency,as we point out in the tables.

The second difference between the SSP-TS methods in this paper and the methods in [32] is that in [32] only one method of orderp>4 is reported,while we have many fifth- and sixth-order methods of various types and stages,optimized for a variety ofKvalues.

The most fundamental difference between our approach and the approach in [32] is that our methods are optimized for the relationship between the forward Euler restriction and the Taylor series restriction while the time-step restriction in the methods of [32] is defined as the most restrictive of the forward Euler and Taylor series time-step conditions.Respecting the minimum of the two cases will still satisfy the nonlinear stability property,but this approach does not allow for a balance between the restrictions considered,which can lead to severely more restrictive conditions.In our approach we use the relationship between the two time-step restrictions to select optimal methods.For this reason,the methods we find have larger allowable time steps in many cases.To understand this a little better consider the case where the forward Euler condition is ΔtFE≤Δxand the Taylor series condition isIn the approach used in [32],the base time-step restriction is thenThe HBT23 method in [32] is a third-order scheme with two stages which has an SSP coefficient of CTS =1 ,so the allowable time step with this scheme will be the sameOn the other hand,using our optimal SSP-TS M2(2,3,0.5) scheme,which has an SSP coeficient=0.75 ,the allowable time step isa 50% increase.This is not only true whenK <1 : consider the case whereand ΔtTS≤Δx.Once again the HBT23 method in [32] will have a time-step restriction ofwhile our M2(2,3,2) method has an SSP coefficient CTS =1.88 ,so that the overall time-step restriction would bewhich is 88% larger.Even when the two base conditions are the same (i.e.,K=1 ) and we have ΔtFE≤Δxand ΔtTS≤Δx,the HBT23 method in [32] gives an allowable time step of CTS =1 while our SSP-TS M2(2,3,1) has an SSP coefficient CTS =1.5 ,so that our method allows a time step that is 50% larger.3These efficiency measures do not account for the fact that the methods in [32] are of type SSP-TS M3 and so require fewer funding evaluations.Correcting for this,our methods are still 10%–40% more efficient.These simple cases demonstrate that our methods,which are optimized for the value ofK,will usually allow a larger SSP coefficient that the methods obtained in [32].

4 Numerical Results

4.1 Overview of Numerical Tests

We wish to test our methods on what are now considered standard benchmark tests in the SSP community.In this subsection we preview our results,which we then present in more detail throughout the remainder of the section.

First,in the tests in Sect.4.3 we focus on how the strong stability properties of these methods are observed in practice,by considering the total variation of the numerical solution.We focus on two scalar PDEs: the linear advection equation and Burgers’ equation,using simple first-order spatial discretizations which are known to satisfy a total variation-diminishing property over time for the forward Euler and Taylor series building blocks.We want to ensure that our numerical approximation to these solutions observe similar properties as long as the predicted SSP time-step restriction,Δt≤CTSΔtFE,is respected.These scalar one-dimensional partial differential equations are chosen for their simplicity so we may understand the behavior of the numerical solution,but the discontinuous initial conditions may lead to instabilities if standard time discretization techniques are employed.Our tests show that the methods we design here preserve these properties as expected by the theory.

In Example 2,we extend the results from Example 1 to the case where we use the higher order weighted essentially non-oscillatory (WENO) method,which is not probably TVD but gives results that have very small increases in total variation.We demonstrate that our methods out-perform other methods,such as the SSP-SD MDRK methods in [5],and that non-SSP methods that are standard in the literature do not preserve the TVD property for any time step.

In many of these examples we are concerned with the total variation-diminishing property.To measure the sharpness of the SSP condition we compute the maximal observed rise in total variation over each step,defined by as well as the maximal observed rise in total variation over each stage,defined by

wherey(s+1)corresponds toun+1.The quantity ofinterest is the time step Δtobs,or the SSP coefficientat which this rise becomes significant,as defined by a maximal increase of 10−10.

It is important to notice that the SSP-TS methods we designed depend on the value ofKin (14).However,in practice we often do not know the exact value ofK.In Example 3 we investigate what happens when we use spatial discretizations with a given value ofKwith time discretization methods designed for an incorrect value ofK.We conclude that although in some cases a smaller step size is required,for methods of type M3 there is generally no adverse result from selecting the wrong value ofK.

In Example 4 we investigate the increased flexibility in the choice of spatial discretization that results from relying on the (12) and (14) base conditions.The only constraint in the choice of differentiation operatorsDxand(described at the end of Sect.1.2) is that the resulting building blocks must satisfy the monotonicity conditions (12) and (14) in the desired convex functional ‖⋅‖ .As noted above,this constraint is less restrictive than requiring that (12) and (13) are satisfied: any spatial discretizations for which (12) and (13) are satisfied will also satisfy (14).However,there are some spatial discretizations that satisfy (12) and (14) that do not satisfy (13).In Example 4 we find that choosing spatial discretizations that satisfy (12) and (14) but not (13) allows for larger time steps before the rise in total variation.And finally,in Example 5,we demonstrate the positivity-preserving behavior of our methods when applied to a nonlinear system of equations.

4.2 On the Numerical Implementation of the Second Derivative

In the following numerical test cases the spatial discretization is performed as follows: at each iteration we take the known valueunand compute the fluxf(un)=−unin the linear case andfor Burgers’ equation.Now to compute the spatial derivativef(un)xwe use an operatorDxand compute

In the numerical examples below the differential operatorDxwill represent,depending on the problem,a first-order upwind finite difference scheme and the fifth-order finite difference WENO method [20].In our scalar test cases(u) does not change sign,so we avoid flux splitting.

Now we have the approximation toUtat timetn,and wish to compute the approximation toUtt.For the linear advection problem,this is very straightforward asUtt=Uxx.To compute this,we takeuxas computed before,and differentiate it again.For Burgers’ equation,we haveWe take the approximation toUtthat we obtained above,and we multiply it byun,then differentiate in space once again.In pseudocode,the calculation takes the form

Using these,we can now construct our two building blocks

In choosing the spatial discretizarionsDxandit is important that these building blocks satisfy (12) and (14) in the desired convex functional ‖⋅‖.

4.3 Example 1: TVD First‑Order Finite Difference Approximations

In this section we use first-order spatial discretizations,that are probably total variation diminishing (TVD),coupled with a variety of time-stepping methods.We look at the maximal rise in total variation.

Example 1a: Linear advectionAs a first test case,we consider a linear advection problem

on a domainx∈[−1,1] ,with step-function initial conditions

and periodic boundary conditions.This simple example is chosen as our experience has shown [13] that this problem often demonstrates the sharpness of the SSP time step.

For the spatial discretization we use a first-order forward difference for the first and second derivative:

These spatial discretizations satisfy

So that ΔtFE=Δxand in this case we haveK=1 in (14).Note that the second-derivative discretization used above does not satisfy the second-derivative condition (13),so that most of the methods we devised in [5] do not guarantee strong stability preservation for this problem.

For all of our simulations for this example,we use a fixed grid ofM=601 points,for a grid sizeand a time stepwhere we varyuntil beyond the point where the TVD property is violated.We step each method forward byN=50 timesteps and compare the performance of the various time-stepping methods constructed earlier in this work,forK=1 .We define the observed SSP coefficientas the multiple of ΔtFEfor which the maximal rise in total variation exceeds 10−10.

We verify that the observed values of ΔtFEandKmatch the predicted values,and test this problem to see how well the observed SSP coefcfiientmatches the predicted SSP coefcfiientfor the fourth-,fifth-,and sixth-order methods.The results are listed in the upper half of Table 6.

Example 1b: Burgers’ equationWe repeat the example above with all the same parameters but for the problem

Table6 Example 1: and for SSP-TS M2 and M3 methods

Table6 Example 1: and for SSP-TS M2 and M3 methods

Method Cpred TS Cobs TS Cpred eff Cobs eff Linear advection FE 1.000 0 1.000 0 1.00 1.00 TS 1.000 0 1.000 0 0.50 0.50 M2(3,4,1) 1.878 8 1.878 8 0.31 0.31 M3(3,4,1) 1.000 0 1.000 0 0.25 0.25 M2(4,4,1) 2.666 8 2.666 8 0.33 0.33 M3(4,4,1) 1.818 1 1.818 1 0.36 0.36 M2(5,4,1) 3.538 1 3.629 1 0.35 0.36 M3(5,4,1) 2.440 6 2.440 6 0.40 0.40 M2(4,5,1) 2.186 4 2.223 9 0.27 0.27 M2(5,5,1) 2.928 0 3.168 1 0.29 0.31 M3(5,5,1) 1.062 5 1.571 0 0.17 0.26 M2(6,5,1) 3.874 9 3.874 9 0.32 0.32 M3(6,5,1) 1.820 7 1.956 2 0.26 0.27 M2(5,6,1) 0.350 0 1.939 8 0.03 0.19 M2(6,6,1) 1.522 5 2.354 8 0.12 0.19 M2(7,6,1) 2.115 0 2.369 5 0.15 0.19 M3(7,6,1) 0.894 6 1.320 7 0.11 0.16 M3(8,6,1) 1.736 9 1.986 1 0.19 0.22 Burgers’ FE 1.000 0 1.000 0 1.00 1.00 TS 1.000 0 1.000 0 0.50 0.50 M2(3,4,1) 1.878 8 1.878 8 0.31 0.31 M3(3,4,1) 1.000 0 1.000 0 0.25 0.25 M2(4,4,1) 2.666 8 2.666 8 0.33 0.33 M3(4,4,1) 1.818 1 1.818 1 0.36 0.36 M2(5,4,1) 3.538 1 3.610 2 0.35 0.36 M3(5,4,1) 2.440 6 2.440 6 0.40 0.40 M2(4,5,1) 2.186 4 2.213 0 0.27 0.27 M2(5,5,1) 2.928 0 3.100 9 0.29 0.31 M3(5,5,1) 1.062 5 1.543 6 0.17 0.25 M2(6,5,1) 3.874 9 3.874 9 0.32 0.32 M3(6,5,1) 1.820 7 2.000 3 0.26 0.28 M2(5,6,1) 0.350 0 1.923 9 0.03 0.19 M2(6,6,1) 1.522 5 2.287 5 0.12 0.19 M2(7,6,1) 2.115 0 2.318 9 0.15 0.16 M3(7,6,1) 0.894 6 1.289 3 0.11 0.16 M3(8,6,1) 1.736 9 1.973 4 0.19 0.21

onx∈(−1,1) .Here we use the spatial derivatives

and

Using Harten’s lemma we can easily show that these definitions ofFandcause the Taylor series condition to be satisfied for Δt≤Δx.The results are quite similar to those of the linear advection equation in Example 1a,as can be seen in the lower half of Table 6.

The results from these two studies show that the SSP-TS methods provide a reliable guarantee of the allowable time step for which the method preserves the strong stability condition in the desired norm.For methods of orderp=4 ,we observe that the SSP coefficient is sharp: the predicted and observed values of the SSP coefficient are identical for all the fourth-order methods tested.For methods of higher order (p=5,6 ) the observed SSP coefficient is often significantly higher than the minimal value guaranteed by the theory.

4.4 Example 2: Weighted Essentially Non‑oscillatory (WENO) Approximations

In this section we re-consider the nonlinear Burgers’ equation (35)

onx∈(−1,1) .We use the step function initial conditions (34),and periodic boundaries.We useM=201 points in the spatial domain,so thatand we step forward forN=50 time steps and measure the maximal rise in total variation for each case.

For the spatial discretization,we use the fifth-order finite difference WENO method [20] in space,as this is a high-order method that can handle shocks.We describe this method in Appendix 3.Recall that the motivation for the development of SSP multistage multi-derivative time-stepping is for use in conjunction with high-order methods for problems with shocks.Ideally,the specially designed spatial discretizations satisfy (12) and (14).Although the weighted essentially non-oscillatory (WENO) methods do not have a theoretical guarantee of this type,in practice we observe that these methods do control the rise in total variation,as long as the step-size is below a certain threshold.

Fig.4 Example 2: Comparison of the maximal rise in total variation (on the y-axis) as a function of(on the x-axis) for a selection of time-stepping methods for evolving Burgers’ equation with WENO spatial discretizations.(a) Fifth-order methods.(b) Sixth-order methods

Below,we refer to the WENO method on a flux withas WENO+defined in (41) and to the corresponding method on a flux with(u)≤0 as WENO−defined in (42).Because(u) is strictly non-negative in this example,we do not need to use flux splitting,and useD=WENO+.For the second derivative we have the freedom to use=WENO+or=WENO−.In this example,we use=Dx=WENO+.In Example 4 below we show that this is more efficient.

In Fig.4(a),we compare the performance of our SSP-TS M3(7,5,1) and SSP-TS M2(4,5,1) methods,which both have eight function evaluations per time step,and our SSP-TS M3(5,5,1),which has six function evaluations per time step,to the SSP-SD MDRK(3,5,2) in [5] and non-SSP RK(6,5) Dormand-Prince method [8],which also have six function evaluations per time step.We note that we use the SSP-SD MDRK(3,5,2) (designed forK=2 ) because this method performs best compared to other explicit two-derivative multistage methods designed for different values ofK.Clearly,the non-SSP method is not safe to use on this example.The M3 methods are most efficient,allowing the largest time step per function evaluation before the total variation begins to rise.

This conclusion is also the case for the sixth-order methods.In Fig.4(b),we compare our SSP-TS M3(9,6,1) and M2(5,6,1) methods,which both have ten function evaluations per time step,and our M3(7,6,1),which has eight function evaluations per time step,to the SSP-SD MDRK(4,6,1) and non-SSP RK(8,6) method given in Verner’s paper table [52],which also have eight function evaluations per time step.Clearly,the non-SSP method is not safe to use on this example.The M3 methods are most efficient,allowing the largest time step per function evaluation before the total variation begins to rise.

This example demonstrates the need for SSP methods: classical non-SSP methods do not control the rise in total variation.We also observe that the methods of type M3 are efficient,and may be the preferred choice of methods for use in practice.

4.5 Example 3: Testing Methods Designed with Various Values of K

In general,the value ofKis not exactly known for a given problem,so we cannot choose a method that is optimized for the correctK.We wish to investigate how methods with different values ofKperform for a given problem.In this example,we re-consider the linear advection Eq.(33)

with step function initial conditions (34),and periodic boundary conditions onx∈(−1,1) .We use the fifth-order WENO method withM=201 points in the spatial domain,so thatand we step forward forN=50 time steps and measure the maximal rise in total variation for each case.Using this example,we investigate how time-stepping methods optimized for differentKvalues perform on the linear advection with finite difference spatial approximation test case above,where it is known thatK=1 .We use a variety of fifth- and sixth-order methods,designed for 0.1 ≤K≤2 and give the value offor which the maximal rise in total variation becomes large,when applied to the linear advection problem.

Fig.5 Example 3: The observed value of such that the method is TVD (y-axis) when methods designed for different K values (on the x-axis) are applied to the problem with K =1 .For each method,the observed value(solid line) is higher than the predicted value (dashed line)

In Fig.5(a) we give the observed value (solid lines) offor a number of SSP-TS methods,M2(4,5,K),M2(5,5,K),M2(6,5,K),M3(5,5,K),and M3(6,5,K),and the corresponding predicted value (dotted lines) that a method designed forK=1 should give.In Fig.5(b) we repeat this study with sixth-order methods M2(5,6,K),M2(6,6,K),M3(7,6,K),and M3(8,6,K).We observe that while choosing the correctKvalue can be beneficial,and is certainly important theoretically,in practice using methods designed for differentKvalues often makes little difference,particularly when the method is optimized for a value close to the correctKvalue.For the sixth-order methods in particular,the observed values of the SSP coefficient are all larger than the predicted SSP coefficient.

4.6 Example 4: The Benefit of Different Base Conditions

In [5] we use the choice ofDx=WENO+defined in (41),followed by=WENO−defined in (42),by analogy to the first-order finite difference for the linear advection caseUt=Ux,where we use a differentiation operatorD+xfollowed by the downwind differentiation operatorD−xto produce a centered difference for the second derivative.In fact,this approach makes sense for these cases because it respects the properties of the flux for the second derivative and consequently satisfies the second-derivative condition (13).However,if we simply wish the Taylor series formulation to satisfy a TVD-like condition,we are free to use the same operator (WENO+or WENO−,as appropriate) twice,and indeed this gives a larger allowable Δt.

In Fig.6 we show how using the repeated upwind discretizationD=WENO−and=WENO−(solid lines) which satisfy the Taylor Series Condition (14) but not the second-derivative condition (13) to approximate the higher order derivative allows for a larger time step than the spatial discretizations (dashed lines) used in 9.We see that for the fifth-order methods the rise in total variation always occurs for largerfor the solid lines (=Dx=WENO−) than for the dashed lines (Dx=WENO−and=WENO+),even for the method designed in [5] to be SSP for the second case but not the first case.For the sixth-order methods the results are almost the same,though the SSP-SD MDRK(4,6,1) method that is SSP for base conditions of the type in [5] performs identically in both cases.These results demonstrate that requiring that the spatial discretizations only satisfy (12) and (14) (but not necessarily (13)) results in methods with larger allowable time steps.

Fig.6 Example 4: The maximal rise in total variation (on the y-axis) for values of (on the x-axis).Simulations using the repeated upwind discretization Dx =WENO− and =WENO− (solid lines) are more efficient than those using Dx =WENO− and =WENO+ (dashed lines).This demonstrates the enhanced allowable time step afforded by the SSP-TS methods

4.7 Example 5: Nonlinear Shallow Water Equations

As a final test case we consider the shallow water equations,where we are concerned with the preservation of positivity in the numerical solution.The shallow water equations [2] are a nonlinear system of hyperbolic conservation laws defined by

whereh(x,t) denotes the water height at locationxand timet,v(x,t) the water velocity,gis the gravitational constant,andU=(h,hv)Tis the vector of unknown conserved variables.In our simulations,we setg=1 .To discretize this problem,we use the standard Lax-Friedrichs splitting

and define the (conservative) approximation to the first derivative as

We discretize the spatial gridx∈(0,1) withM=201 points.To approximate the second derivative,we start with element-wise first derivativeand then approximate the second derivative (consistent with (11)) as

Table7 The predicted and observed values offor which positivity of the height of the water is preserved in the shallow water equations in Example 5

Table7 The predicted and observed values offor which positivity of the height of the water is preserved in the shallow water equations in Example 5

?

In problems such as the shallow water equations,the non-negativity of the numerical solution is important as a height ofh<0 is not physically meaningful,and the system loses hyperbolicity when the height becomes negative.For a positivity-preserving test case,we consider a Riemann problem with zero initial velocity,but with a wet and a dry state [2,53]:

In our numerical simulations,we focus on the impact of the numerical scheme on the positivity of the solver for the the water heighth(x,t).This quantity is ofinterest from a numerical perspective because if the heighth(x,t)<0 for anyxort,the code will crash due to square-root of height.

First,we investigate the behavior of the base methods in terms of the positivity-preserving time step.In other words,we want to get a numerical value for ΔtFE,K.To do so,we numerically study the positivity behavior of the forward Euler and Taylor series approach.To do this,we evolve the solution forward for more time steps with different values ofto identify the predicted positivity-preserving value.Using the approach,we see that as we increase the number of steps the predicted value of the positivity preserving value,for both forward Euler and Taylor series.We are not able to numerically identifyresulting from the second-derivative condition,which cannot be evolved forward as it does not approximate the solution to the ODE at all.

In Table 7 we compare the positivity-preserving time step of a variety of numerical time integrators.We consider the fifth-order SSP-TS methods M2(4,5,1),M3(5,5,1),and M3(6,5,1),and compare their performance to the SSP-SD MDRK(3,5,2) method in [5],and the non-SSP Dormand-Prince method.We also consider the sixth-order SSP-TS methods M2(5,6,1),M3(7,6,1),and M3(9,6,1),as well as the SSP-SD MDRK(4,6,1) from [5] and the non-SSPRK(8,6) method.Positivity of the water height is measured at each stage for a total ofN=60 time steps.We report the largest allowable value of(is the maximal wavespeed for the domain) for which the solution remains positive.For each method,the predicted valuesare obtained by multiplying the SSP coefficientof that method byFor the SSP-SD MDRK methods we do not make a prediction as we are not able to identifyresulting from the second-derivative condition.

In Table 7 we show that all of our SSP-TS methods preserve the positivity of the solution for values larger than those predicted by the theoryand that even for the SSP MSRK methods there is a large region of valuesfor which the solution remains positive.However,the non-SSP methods permit no positive time step that retains positivity of the solution,highlighting the importance of SSP methods.

5 Conclusions

In [5] we introduced a formulation and base conditions to extend the SSP framework to multistage multi-derivative time-stepping methods,and the resulting SSP-SD methods.While the choice of base conditions we used in [5] gives us more flexibility in finding SSP time-stepping schemes,it limits the flexibility in the choice of the spatial discretization.In the current paper we introduce an alternative SSP formulation based on the conditions (12) and (14) and investigate the resulting explicit two-derivative multistage SSP-TS time integrators.These base conditions are relevant because some commonly used spatial discretizations may not satisfy the second-derivative condition (13) which we required in [5],but do satisfy the Taylor series condition (14).This approach decreases the flexibility in our choice of time discretization because some time discretizations that can be decomposed into convex combinations of (12) and (13) cannot be decomposed into convex combinations of (12) and (14).However,it increases the flexibility in our choice of spatial discretizations,as we may now consider spatial methods that satisfy (12) and (14) but not (13).In the numerical tests we showed that this increased flexibility allowed for more efficient simulations in several cases.

In this paper,we proved that explicit SSP-TS methods have a maximum obtainable order ofp=6 .Next we formulated the proper optimization procedure to generate SSP-TS methods.Within this new class we were able to organize our schemes into three sub categories that reflect the different simplifications used in the optimization.We obtained methods up to and including orderp=6 thus breaking the SSP order barrier for explicit SSP Runge-Kutta methods.Our numerical tests show that the SSP-TS explicit two-derivative methods perform as expected,preserving the strong stability properties satisfied by the base conditions (12) and (14) under the predicted time-step conditions.Our simulations demonstrate the sharpness of the SSP-TS condition in some cases,and the need for SSP-TS time-stepping methods.Furthermore the numerical results indicate that the added freedom in the choice of spatial discretization results in larger allowable time steps.The coefficients of the SSP-TS methods described in this work can be downloaded from [14].

AcknowledgementsThe work of D.C.Seal was supported in part by the Naval Academy Research Council.The work of S.Gottlieb and Z.J.Grant was supported by the AFOSR Grant #FA9550-15-1-0235.A part of this research is sponsored by the Office of Advanced Scientific Computing Research; US Department of Energy,and was performed at the Oak Ridge National Laboratory,which is managed by UT-Battelle,LLC under Contract no.De-AC05-00OR22725.This manuscript has been authored by UT-Battelle,LLC,under contract DE-AC05-00OR22725 with the US Department of Energy.The United States Government retains and the publisher,by accepting the article for publication,acknowledges that the United States Government retains a non-exclusive,paid-up,irrevocable,world-wide license to publish or reproduce the published form of this manuscript,or allow others to do so,for United States Government purposes.

Appendix 1 Order Conditions

Any method of the form (9) must satisfy the order conditions for allp≤Pto be of orderP.

Appendix 2 Coefficients of Selected Methods

All the time-stepping methods in this work can be downloaded as Matlab files from [14].In this appendix we present selected methods.

SSP-TS M2(4,5,1) This method has=2.186 48,

SSP-TS M3(8,6,1) This method has=1.736 9,

Appendix 3: Fifth‑Order WENO Method of Jiang and Shu

To solve the PDE

we approximate the spatial derivative to obtainF(u)≈−f(u)x,and then use a time-stepping method to solve the resulting system of ODEs.In this section we describe the fifthorder WENO scheme presented in [20].

First,we split the flux into the positive and negative parts

This can be accomplished in a variety of ways,e.g.,the Lax-Friedrichs flux splitting

wherem=maxIn this way we ensure that

To calculate the numerical fluxeswe begin by calculating the smoothness measurements to determine if a shock lies within the stencil.For our fifth-order scheme,these are

and

Next,we use the smoothness measurements to calculate the stencil weights

and

Finally,the numerical fluxes are

and

Finally,we compute

and put it all together