Strong Stability Preserving Integrating Factor Two-step Runge--Kutta Methods

TL;DR

This paper introduces explicit SSP two-step Runge-Kutta integrating factor methods that achieve higher order and larger SSP coefficients, overcoming the order barrier of traditional SSP Runge-Kutta methods.

Contribution

The authors develop and analyze explicit SSP two-step Runge-Kutta integrating factor methods that surpass the fourth order barrier and have larger SSP coefficients.

Findings

Methods up to eighth order are constructed.

Two-step methods can break the order barrier of p=4.

Two-step methods often have larger SSP coefficients than one-step counterparts.

Abstract

Problems that feature significantly different time scales, where the stiff time-step restriction comes from a linear component, implicit-explicit (IMEX) methods alleviate this restriction if the concern is linear stability. However, where the SSP property is needed, IMEX SSP Runge-Kutta (SSP-IMEX) methods have very restrictive time-steps. An alternative to SSP-IMEX schemes is to adopt an integrating factor approach to handle the linear component exactly and step the transformed problem forward using some time-evolution method. The strong stability properties of integrating factor Runge--Kutta methods were previously established, where it was shown that it is possible to define explicit integrating factor Runge-Kutta methods that preserve strong stability properties satisfied by each of the two components when coupled with forward Euler time-stepping. It was proved that the solution will be SSP if the transformed problem is stepped forward with an explicit SSP Runge-Kutta method that has non-decreasing abscissas. However, explicit SSP Runge-Kutta methods have an order barrier of p=4, and sometimes higher order is desired. In this work we consider explicit SSP two-step Runge--Kutta integrating factor methods to raise the order. We show that strong stability is ensured if the two-step Runge-Kutta method used to evolve the transformed problem is SSP and has non-decreasing abscissas. We find such methods up to eighth order and present their SSP coefficients. Adding a step allows us to break the fourth order barrier on explicit SSP Runge-Kutta methods; furthermore, our explicit SSP two-step Runge--Kutta methods with non-decreasing abscissas typically have larger SSP coefficients than the corresponding one-step methods.

Tables

Table 1. Table 1: SSP coefficients of the optimized eSSP-TSRK(s,p) methods [ 22 ] .

	2	3	4	5	6	7
1	-	-	-	-	-	-
2	1.4142	0.7320	-	-	-	-
3	2.4495	1.6506	0.8588	-	-	-
4	3.4641	2.3027	1.5926	0.8542	-	-
5	4.4721	2.9879	2.3605	1.6481	-	-
6	5.4772	3.7768	3.0559	2.3093	0.5957	-
7	6.4807	4.4836	3.7405	2.9278	1.2719	-
8	7.4833	5.2227	4.4921	3.5794	1.9384	0.5666
9	8.4853	6.0498	5.2705	3.9415	2.5826	1.1199
10	9.4868	6.8274	6.1039	4.2544	3.1992	1.7857

Table 2. Table 2: Effective SSP coefficients of the optimized eSSP-TSRK(s,p) methods [ 22 ] .

	2	3	4	5	6	7
1	-	-	-	-	-	-
2	0.7071	0.3660	-	-	-	-
3	0.8165	0.5502	0.2863	-	-	-
4	0.8660	0.5757	0.3982	0.2135	-	-
5	0.8944	0.5976	0.4721	0.3296	-	-
6	0.9129	0.6295	0.5093	0.3849	0.0993	-
7	0.9258	0.6405	0.5343	0.4183	0.1817	-
8	0.9354	0.6528	0.5615	0.4474	0.2423	0.0708
9	0.9428	0.6722	0.5856	0.4379	0.2869	0.1244
10	0.9487	0.6827	0.6104	0.4254	0.3199	0.1786

Table 3. Table 3: SSP coefficients of the optimized eSSP-TSRK + (s,p) methods.

	2	3	4	5	6	7
1	-	-	-	-	-	-
2	1.4142	0.7320	-	-	-	-
3	2.4495	1.6506	0.8588	-	-	-
4	3.4641	2.3027	1.5926	0.8542	-	-
5	4.4721	2.9807	2.3523	1.6481	-	-
6	5.4772	3.7672	3.0140	2.3093	0.5958	-
7	6.4807	4.4533	3.6751	2.9173	1.2671	-
8	7.4833	5.2134	4.4178	3.5477	1.8728	0.5666
9	8.4853	6.0012	5.2120	3.9426	2.4784	1.0715
10	9.4868	6.7916	6.0626	4.2362	3.1646	1.6892

Table 4. Table 4: Effective SSP coefficients of the optimized eSSP-TSRK + (s,p) methods.

	2	3	4	5	6	7
1	-	-	-	-	-	-
2	0.7071	0.3660	-	-	-	-
3	0.8165	0.5502	0.2863	-	-	-
4	0.8660	0.5757	0.3982	0.2135	-	-
5	0.8944	0.5961	0.4705	0.3296	-	-
6	0.9129	0.6279	0.5023	0.3849	0.0993	-
7	0.9258	0.6362	0.5250	0.4168	0.1810	-
8	0.9354	0.6517	0.5522	0.4435	0.2341	0.0708
9	0.9428	0.6668	0.5791	0.4381	0.2754	0.1191
10	0.9487	0.6792	0.6063	0.4236	0.3165	0.1689

Table 5. Table 5: The observed SSP coefficient of a variety of SSPIF-TSRK ( s , p ) 𝑠 𝑝 (s,p) methods compared to theirpredicted SSP coefficient, shown for various wavespeeds a 𝑎 a . The value of a 𝑎 a does not negatively impact the observed SSP coefficient.

Method	$𝒞$	${\lambda }_{obs}^{TVD}$
		for $a=0$	$a=1.0$	$a=5$
SSPIF-TSRK(3,4)	0.8588	1.0454	1.2550	1.2621
SSPIF-TSRK(5,4)	2.3523	2.3523	2.3523	2.4123
SSPIF-TSRK(9,4)	5.2120	5.2120	5.2120	6.4010
SSPIF-TSRK(4,5)	0.8542	1.1852	1.3388	1.3389
SSPIF-TSRK(6,5)	2.3093	2.3093	2.3093	2.3094
SSPIF-TSRK(9,5)	3.9426	3.9426	3.9426	4.1173
SSPIF-TSRK(6,6)	0.5958	1.7771	1.7891	1.7893
SSPIF-TSRK(7,6)	1.2671	2.0239	2.0239	2.0261
SSPIF-TSRK(9,6)	2.4784	2.8038	2.8038	2.8204
SSPIF-TSRK(8,7)	0.5666	1.6624	2.6737	2.7788
SSPIF-TSRK(9,7)	1.0715	2.1626	2.4053	2.4053
SSPIF-TSRK(11,8)	0.2743	2.3871	3.1137	3.1271

Table 6. Table 6: The observed SSP coefficients for integrating factor methods, an SSP-IMEX method, and an explicit Runge–Kutta method for Example 2b with various wavespeeds a 𝑎 a . The value of a 𝑎 a does not negatively impact the observed SSP coefficient for the integrating factor methods, but it does for the IMEX and the explicit Runge–Kutta methods.

Method	${\lambda }_{obs}^{TVD}$
	$a=0$	$a=1.0$	$a=2$	$a=10$
SSPIF-TSRK(4,3)	2.303	2.303	2.303	2.775
SSPIF-TSRK(4,4)	1.593	1.593	1.593	1.639
SSPIFRK(4,3)	1.818	1.818	1.818	1.818
SSP-IMEX(4,3,K)	2.000	1.476	1.192	0.310
eSSPRK(4,3)	2.000	1.000	0.666	0.181

Equations

$$U_t+f(U{)}_x=0,$$

$$u_t=F(u),$$

$$\Vert u^{n+1}\Vert =\Vert u^n+\Delta tF(u^n)\Vert \le \Vert u^n\Vert ,$$

$$0\le \Delta t\le \Delta t_{\text{FE}}.$$

$$u^{(0)}$$

$$u^{(i)}$$

$$u^{n+1}$$

$$\Vert u^{(i)}\Vert =\Vert \sum _{j=0}^{i-1}\left({\alpha }_{i,j}{u}^{(j)}+\Delta t{\beta }_{i,j}F(u^{(j)})\right)\Vert =\Vert \sum _{j=0}^{i-1}{\alpha }_{i,j}\left(u^{(j)}+\Delta t\frac{{\beta }_{i,j}}{{\alpha }_{i,j}}F(u^{(j})\right)\Vert .$$

$$\Vert u^{(i)}\Vert \le \sum _{j=0}^{i-1}{\alpha }_{i,j}\Vert u^{(j)}+\Delta t\frac{{\beta }_{i,j}}{{\alpha }_{i,j}}F(u^{(j})\Vert \le \Vert u^n\Vert ,$$

$$\Delta t\le \underset{i,j}{\min }\frac{{\alpha }_{i,j}}{{\beta }_{i,j}}\Delta t_{\text{FE}}.$$

$$u_t=Lu+N(u)$$

$$\Vert u^n+\Delta tN(u^n)\Vert \le \Vert u^n\Vert \text{mbox}for\Delta t\le \Delta t_{\text{FE}}$$

$$\Vert u^n+\Delta tL{u}^n\Vert \le \Vert u^n\Vert \text{mbox}for\Delta t\le {\widetilde{\Delta t}}_{\text{FE}}.$$

$$e^{-Lt}{u}_t-e^{-Lt}Lu=e^{-Lt}N(u)$$

$${\left(e^{-Lt}u\right)}_t=e^{-Lt}N(u).$$

$$w_t=e^{-Lt}N(e^{Lt}w)=G(w),$$

$$u_t+a{u}_x+{\left(\frac{1}{2}{u}^2\right)}_x$$

$$u^{(1)}$$

$$u^{(2)}$$

$$u^{n+1}$$

$$u^{(1)}$$

$$u^{(2)}$$

$$u^{n+1}$$

$$u^{(i)}$$

$$\Vert e^{\tau L}{u}^n\Vert \le \Vert u^n\Vert \forall \tau \ge 0.$$

$$\Vert e^{\tau L}(u^n+\Delta tN(u^n))\Vert \le \Vert u^n\Vert \forall \Delta t\le \Delta t_{\text{FE}},\text{mbox}providedthat\tau \ge 0.$$

$$u^{(0)}$$

$$u^{(i)}$$

$$u^{n+1}$$

$$\Vert u^{n+1}\Vert \le \Vert u^n\Vert \forall \Delta t\le \mathcal{C}\Delta t_{\text{FE}}.$$

$$\mathcal{C}=\underset{i,j}{\min }\frac{{\alpha }_{i,j}}{{\beta }_{i,j}}.$$

$$\Vert e^{-\tau \tilde{L}}{u}^n\Vert \le \Vert u^n\Vert \forall \tau \ge 0,$$

$$y_1^n$$

$$y_i^n$$

$$u^{n+1}$$

$$\mathbf{w}$$

Full text

This paper is dedicated to the memory of Saul Abarbanel. His wisdom, humor, and kindness were a gift to all who knew him.

1 Introduction

The behavior of the numerical solution of a hyperbolic partial differential equation (PDE) of the form

[TABLE]

depends on properties of the spatial discretization and of the time discretization. When the solution is smooth, stability is guaranteed by analyzing the $L_{2}$ stability properties of the discretization applied to the linear problem. However, when dealing with a non-smooth solution, the numerical solution may contain non-physical oscillations that prevent the approximation from converging uniformly and thus $L_{2}$ linear stability is not sufficient to ensure convergence [26]. To ensure that the numerical solution does not form stability-destroying oscillations, we require that the numerical method satisfies nonlinear non-inner-product stability properties such as stability in the maximum norm or in the total variation (TV) semi-norm.

For nonlinear hyperbolic problems with discontinuous solutions we must, therefore, analyze the nonlinear non-inner-product stability properties of a highly nonlinear complex spatial discretization combined with a high order time discretization. Instead of this difficult task a method-of-lines formulation is generally followed: we develop a spatial discretization that satisfies nonlinear non-inner-product stability properties when coupled with the forward Euler time stepping method. Next, we use a high order strong stability preserving (SSP) time discretization [29, 30, 28, 34, 14, 19, 18, 22, 9, 23] which preserves the properties of the spatial discretization coupled with forward Euler. Explicit strong stability preserving (SSP) Runge–Kutta methods were developed in [29, 30] to preserve the properties of total variation diminishing (TVD) spatial discretizations for hyperbolic conservation laws (1) with discontinuous solutions. They have since been developed extensively and have been widely used to preserve different numerical stability properties needed in a variety of application areas. Furthermore, many classes of SSP time-stepping methods have been studied, including SSP explicit and implicit linear multi-step methods [9], Runge–Kutta methods [29, 30, 18, 19], multi-stage multi-step methods [2, 22, 23], and multi-stage multi-derivative methods [3, 12].

1.1 SSP methods

The key to developing SSP time-stepping methods is ensuring that the methods can be re-written as convex combinations of forward Euler steps. To illustrate this concept, we begin with the PDE (1) above, and use a spatial discretizations of $f(U)_{x}$ that ensures that when we evolve the resulting semi-discretized system of ordinary differential equations (ODEs)

[TABLE]

using the forward Euler method, the strong stability property in the convex functional $\|\cdot\|$ is satisfied

[TABLE]

under a step size restriction

[TABLE]

Next, we use a higher order time integrator that can be written as a convex combination of forward Euler steps, so that we ensure that any convex functional strong stability property that is satisfied by the forward Euler method will still be satisfied by the higher order time discretization, perhaps under a different time-step. For example, an $s$-stage explicit Runge–Kutta method (denoted eSSPRK$(s,p)$ where $p$ is the order) can be written in Shu-Osher form [9]:

[TABLE]

(Where $\sum_{j=0}^{i-1}\alpha_{i,j}=1$ is required for consistency). If $\alpha_{i,j}$ and $\beta_{i,j}$ are non-negative, and any $\alpha_{i,j}$ is zero only if its corresponding $\beta_{i,j}$ is zero, then we can rearrange each stage into a convex combination of forward Euler steps:

[TABLE]

Clearly, then, (3) tells us that

[TABLE]

while (4) imposes the time-step restriction

[TABLE]

(We use the convention that the ratio in considered infinite in the case where a $\beta_{i,j}$ is equal to zero). This convex combination ensures that the internal stages also satisfy the strong stability property under the same time-step restriction; this can be important when pressure, density, or water height, are simulated, because in these cases a negative value at an intermediate or final stage may prevent the simulation from continuing [13]. Clearly, the convex combination condition is a sufficient condition for strong stability preservation for explicit Runge–Kutta methods; in [9, 24, 33] it was shown to be a necessary condition as well.

This shows that when a higher order time discretization method is written as a convex combination of forward Euler steps, it will preserve the forward Euler condition (3), under a modified time-step restriction $\Delta t\leq\mathcal{C}\Delta t_{\textup{FE}}$. When $\mathcal{C}>0$, the method is called strong stability preserving (SSP) with SSP coefficient $\mathcal{C}$ [29]. While in the original papers [29, 30], the term $\|\cdot\|$ was the total variation semi-norm, the strong stability preservation property holds for any convex functional $\|\cdot\|$ as long as the spatial discretization satisfies (3) in that convex functional for some $\Delta t_{\textup{FE}}$ as in (4).

Not all methods can be decomposed into convex combinations of forward Euler steps with $\mathcal{C}>0$. Explicit SSP Runge–Kutta methods cannot exist for order $p>4$ [24, 28]. Furthermore, the SSP coefficient $\mathcal{C}$ is restricted too, as all explicit $s$-stage Runge–Kutta methods have a SSP bound $\mathcal{C}\leq s$ [9]. Unlike for smooth problems, where implicit methods (or implicit treatment of stiff terms) can eliminate the time-step restriction needed for linear or nonlinear inner-product norm stability, a strong stability preserving general linear method (GLM) of order greater than one has a finite SSP time-step [32]. In fact, it has been shown [25, 6, 19, 4] that the SSP time-step restrictions for implicit and implicit-explicit (IMEX) SSP Runge–Kutta methods have a restrictive observed bound of $\mathcal{C}\leq 2s$ . The limitation on the SSP coefficient of implicit and IMEX SSP Runge–Kutta methods led to the study of SSP integrating factor Runge–Kutta methods in [15].

1.2 Overview of current paper

The explicit SSP integrating factor Runge–Kutta methods and the associated theory developed in our prior work [15], reviewed in Section 2, allow us to alleviate the time-step restriction while preserving the strong stability property. However, because explicit SSP Runge–Kutta methods have an order barrier of $p\leq 4$, we cannot hope to get higher order explicit SSP integrating factor Runge–Kutta methods. We could try to work with an integrating factor multi-step approach, which is SSP as long as the multi-step method is SSP (as we show in Subsection 4.1) but these methods have small SSP coefficients and require many steps for high order. Instead, we follow the process in [22, 23] and examine the SSP properties of integrating factor methods based on explicit SSP two-step Runge–Kutta (eSSP-TSRK) methods such as those in [22]. When these are used as a basis for integrating factor TSRK methods, the result is SSP provided that the abscissas are non-decreasing. In Section 3 we discuss explicit SSP two-step Runge–Kutta methods, review the optimization problem, and review the resulting eSSP-TSRK methods presented in [22]. In Section 4 we develop the SSP theory for explicit SSP integrating factor two-step Runge–Kutta (SSPIF-TSRK) methods, which is very similar to that of the explicit SSP integrating factor Runge–Kutta (SSPIFRK) methods in [15], and show that as long as these methods have non-decreasing abscissas the result is SSP. We show that to find appropriate explicit SSP two-step Runge–Kutta for pairing with integrating factor methods, the optimization problem needs to be augmented by a simple condition on the abscissas. In Section 5 we present the optimized methods found using this approach, and discuss their features. In Section 6 we test these methods for convergence and for their SSP performance on numerical test cases.

1.3 Efficient computation of the matrix exponential

It is important to note that in the current work as well as our prior work [15, 16], the cost of computation of the matrix exponential will be a critical factor in the ability to efficiently implement the proposed methods. If the computation of the matrix exponential is only needed once per simulation (i.e. if $L$ is a constant coefficient operator), the cost may be reasonable, but in some cases the exponential must be computed at every step, and low-storage, matrix-free approaches are needed. New approaches to efficiently compute the matrix exponential have been recently considered in [1, 31, 27, 7], and such approaches, as well as others, will be critical for bringing the integrating factor methods proposed in [15, 16] and the current work into practical use.

2 Review of SSP integrating factor Runge–Kutta methods

In [15, 16] we considered problems of the form

[TABLE]

that result from a semi-discretization of a hyperbolic PDE. We focus on the case where the problem has a linear constant coefficient component $Lu$ and a nonlinear component $N(u)$, such that

[TABLE]

and

[TABLE]

The notation $\|\cdot\|$ here refers to some convex functional (not usually a norm) needed for nonlinear non-inner-product stability. Of particular interest is the case where $L$ is a linear operator that significantly restricts the allowable time-step due to strong stability concerns, i.e. where $\tilde{\Delta t_{\textup{FE}}}<<{\Delta t_{\textup{FE}}}$. For such cases, where nonlinear non-inner-product stability properties are of concern, an implicit or implicit-explicit (IMEX) SSP scheme doesn’t significantly alleviate the allowable time-step [4] and such methods will result in severe constraints on the allowable time-step [25, 19, 6, 9, 4].

This motivated our investigation of integrating factor methods, where the linear component $Lu$ is handled exactly and then the time-step restriction is dependent only upon the step size restriction $\Delta t\leq{\Delta t_{\textup{FE}}}$ coming from the the nonlinear component $N(u)$. The methodology we considered in [15] in order to alleviate the restriction on the allowable time-step begins with an integrating factor approach:

[TABLE]

Defining $w=e^{-Lt}u$ gives us the modified ODE system

[TABLE]

which can then be evolved forward in time. In [15] we used an explicit SSP Runge–Kutta method of the form (1.1). We observed that this approach is not sufficient to ensure strong stability of the system. However, as we showed, if the transformed problem is stepped forward using a SSP Runge–Kutta method that is carefully chosen to have non-decreasing abscissas the resulting method will preserve the desired strong stability property.

2.1 Motivating example

Before we repeat the formal results from [15], which explain how non-decreasing abscissas in the Runge–Kutta method are needed to preserve strong stability when using an integrating factor approach, we present the following example, presented in [15], which illustrates this fact:

Consider the partial differential equation

[TABLE]

where $x\in[0,1]$ and the boundary conditions are periodic. In this example we consider $a=10$. A first-order upwind difference with $400$ points in space is used to semi-discretize the linear term $Lu\approx-10u_{x}$. The nonlinear terms $N(u)\approx-\left(\frac{1}{2}u^{2}\right)_{x}$ are approximated using a fifth order WENO finite difference [17].

The transformed problem (10) is stepped forward in time by the explicit SSP Runge–Kutta method with non-decreasing abscissas, denoted eSSPRK${}^{+}(3,3)$ in [15], which has $\mathcal{C}=\frac{3}{4}$. The resulting SSP integrating factor Runge–Kutta method is:

[TABLE]

For comparison, the transformed problem (10) is also stepped forward in time by the explicit eSSPRK(3,3) Runge–Kutta method in [29] which has $\mathcal{C}=1$ (this method is widely known as the Shu-Osher method). The resulting (non-SSP) integrating factor Runge–Kutta method becomes:

[TABLE]

Exponentials with negative exponents appear in this formulation due to the fact that the SSP Shu-Osher method it is based on has decreasing abscissas. Using different values of $\Delta t$ we evolved the solution 25 time steps using the integrating factor Runge–Kutta methods (2.1) and (2.1), and calculated the maximal rise in total variation over each stage for all time steps.

In Figure 1 we show the $\log_{10}$ of the maximal rise in total variation vs. the value of $\lambda=\frac{\Delta t}{\Delta x}$ of the evolution using (2.1) (in blue) and using (2.1) (in red). The results from (2.1) maintain a small maximal rise in total variation up to $\lambda\approx 0.8$, while in the results from (2.1) we observe that the maximal rise in total variation is large even for very small values of $\lambda$. This behavior is due to the fact that basing an integrating factor Runge–Kutta method on an explicit SSP Runge–Kutta method is not enough to ensure the preservation of a strong stability property – the method must also satisfy the non-decreasing abscissa condition to ensure that the strong stability property is preserved.

2.2 SSP analysis of integrating factor Runge–Kutta method

To understand the results in [15] recall that each stage of the explicit Runge–Kutta method (1.1) applied to the transformed problem (10) becomes

[TABLE]

In the following, we consider each piece of this formula. Each of the following results is taken directly from [15] and the proofs appear in that work. We first consider the behavior of the exponential operator $e^{\tau L}$. The following theorem ensures that this term is strongly stable as long as $\tau\geq 0$:

Theorem 2.1**.**

If a linear operator $L$ satisfies (9) for some value of $\tilde{\Delta t}_{\textup{FE}}>0$, then

[TABLE]

Note that in the following results we no longer require the condition (9) but rather the weaker condition (14). Pairing this exponential with a forward Euler type step still preserves strong stability:

Corollary 1**.**

Given a linear operator $L$ that satisfies (14) and a (possibly nonlinear) operator $N(u)$ that satisfies (8) for some value of ${\Delta t}_{FE}>0$, we have

[TABLE]

Finally, we put the pieces together:

Theorem 2.2**.**

Given a linear operator $L$ that satisfies (14) and a (possibly nonlinear) operator $N(u)$ that satisfies (8) for some value of ${\Delta t}_{FE}>0$, and a Runge–Kutta integrating factor method of the form

[TABLE]

where $0=c_{1}\leq c_{2}\leq...\leq c_{s}$, then $u^{n+1}$ obtained from (2.2) satisfies

[TABLE]

where

[TABLE]

In [15] we used these results to motivate the search for explicit SSP Runge–Kutta methods with non-decreasing abscissas, that would pair well with an integrating factor approach to produce SSPIFRK methods. We found explicit SSP Runge–Kutta methods of up to $s=10$ stages and order $p=4$ and demonstrated the performance of the corresponding explicit SSP integrating factor Runge–Kutta methods on widely used test cases.

In [16] we show that it is not absolutely necessary to have non-decreasing abscissas: even if there are some decreasing abscissas one can preserve the strong stability property by replacing the operator $L$ by an operator $\tilde{L}$ that satisfies

[TABLE]

whenever the difference of abscissas $c_{i}-c_{j}$ is negative. For hyperbolic partial differential equations, $\tilde{L}$ is simply the spatial discretization that is stable for the backwards in time version of the equation. This approach is employed in the classical SSP literature, called downwinding. In that case, negative coefficients $\beta_{i,j}$ preserve the SSP property as long as the operator is replaced by a downwinded operator [10, 11, 21].

3 A review of explicit SSP two-step Runge–Kutta methods

The order barrier of $p\leq 4$ on explicit SSP Runge–Kutta methods and the need for higher order explicit methods in simulations led to the investigation of the SSP properties of explicit two-step Runge–Kutta methods (TSRK) [22], and later of explicit multi-step Runge–Kutta methods (MSRK) [23]. In [22] it was shown that explicit SSP TSRK methods do not use values of stages from the previous steps, so that they can be written as:

[TABLE]

An extension of this to multiple steps was studied in [23].

We can write the method in the matrix form

[TABLE]

where $\mathbf{S}$ and $\mathbf{T}$ are the matrices defined by:

[TABLE]

and $\mathbf{w}$, $\mathbf{x}$, and $\mathbf{f}$ are defined by:

[TABLE]

The $\mathbf{D}$ is a $2\times s$ matrix that contains the elements $\mathbf{D}_{i,1}=d_{i}$ and $\mathbf{D}_{i,2}=1-d_{i}$, for $i=1,\dots,s$. The vector $\boldsymbol{\theta}$ contains the elements $\boldsymbol{\theta}_{1}=\theta$ and $\boldsymbol{\theta}_{2}=1-\theta$. Note that $\mathbf{S}\mathbf{e}=\mathbf{e}$, where $\mathbf{e}$ is a vector ones of the proper length.

3.1 Formulating the optimization problem

In [22] a SSP optimization problem was defined for TSRK methods. The derivation of the optimization problem begins with the matrix form (19) of the TSRK:

[TABLE]

We then add $r\mathbf{T}\mathbf{w}$ to each side of the method

[TABLE]

and rearrange to obtain

[TABLE]

Where $\mathbf{P}=r\left(\mathbf{I}+r\mathbf{T}\right)^{-1}\mathbf{T}$ and $\mathbf{R}=\left(\mathbf{I}+r\mathbf{T}\right)^{-1}\mathbf{S}=\left(\mathbf{I}-\mathbf{P}\right)\mathbf{S}$. Note that the row sums of $[\mathbf{R}~{}\mathbf{P}]$ are each equal to one

[TABLE]

It is easy to see that if $\mathbf{R}$ and $\mathbf{P}$ have no negative values, then the method is a convex combination of forward Euler steps, and so is SSP with SSP coefficient defined by

[TABLE]

This leads directly to the optimization problem: We wish to find coefficients $\mathbf{T}$ and $\mathbf{S}$ for explicit SSP two-step Runge–Kutta methods that maximize the value of $r$ while satisfying

[TABLE]

This requires that $\left(\mathbf{I}+r\mathbf{T}\right)^{-1}$ exists. The inequalities (22a) and (22b) are understood componentwise. The $\tau_{\rho}$ in (22c) are the order conditions listed in Appendix A.

3.2 Optimized SSP two-step Runge–Kutta methods

In [22] a SSP optimization problem was defined for the TSRK methods, and optimized explicit TSRK methods of up to eighth order were found numerically. The SSP coefficients $\mathcal{C}$ of these methods up to seventh order are in Table 1. Of course, a method with more stages will have a higher SSP coefficient, but also require more computations. To account for this we define the effective SSP coefficient $\mathcal{C}_{\textup{eff}}=\frac{1}{s}\mathcal{C}$, which is the SSP coefficient normalized by the number of stages, in Table 2. An eleven stage eighth order method has an SSP cefficient $\mathcal{C}=0.341$ and an effective SSP Coefficient $\mathcal{C}_{\textup{eff}}=0.031$.

These methods break the fourth-order barrier of (one-step) SSP Runge–Kutta methods and have significantly larger SSP coefficients than the corresponding order multi-step methods. Numerical tests in [22] showed that these high order SSP two-step Runge-Kutta methods are beneficial for preserving the order and strong stability properties when used with high order spatial discretizations for the time integration of a variety of hyperbolic PDEs. It was proved in [22] that explicit SSP TSRK methods have order at most eight. In [23] it was shown that adding more steps overcomes this order barrier as well.

4 Explicit SSP two-step Runge–Kutta schemes for use with integrating factor methods

In this section we consider, as we did in [15, 16], a problem of the form (7) with a linear constant coefficient component $Lu$ that satisfies the strong stability condition (9) and a nonlinear component $N(u)$ that satisfies the strong stability condition (8). When $\tilde{\Delta t}_{\textup{FE}}<<\Delta t_{\textup{FE}}$ using an explicit SSP Runge–Kutta method will result in severe constraints on the allowable time-step that is not alleviated by using an implicit or an implicit-explicit (IMEX) SSP Runge–Kutta method [25, 19, 6, 9, 4]. This is in contrast to the case with linear $L_{2}$ stability, where using an implicit or an IMEX method may completely alleviate the time-step restriction coming from the stiff component.

To alleviate the restriction on the allowable time-step we can solve the linear part exactly by an integrating factor approach as in (10). This new ODE can be evolved forward in time using standard methods. As described in Section 2, in [15] we considered stepping the transformed problem (10) forward using an explicit SSP Runge–Kutta method of the form (1.1) and found that the result was SSP provided that the method had non-decreasing abscissas. In the following, we consider extending the SSP integrating factor approach by using multi-step and two-step Runge–Kutta integrating factor methods to obtain order $p\geq 4$

4.1 SSP integrating factor linear multi-step methods

It can be easily shown that if we use an explicit SSP multi-step method

[TABLE]

which satisfies

[TABLE]

with SSP coefficient $\mathcal{C}=\max_{l}\frac{\alpha_{l}}{\beta_{l}}$, the result is SSP as well, with the same SSP coefficient:

Theorem 4.1**.**

Given a linear operator $L$ that satisfies (14) and a (possibly nonlinear) operator $N(u)$ that satisfies (8) for some value of ${\Delta t}_{FE}>0$, and a $k$-step explicit SSP integrating factor multi-step method of the form

[TABLE]

then the numerical solution satisfies

[TABLE]

*where $\mathcal{C}=\max_{l}\frac{\alpha_{l}}{\beta_{l}}$. *

Notice that for multi-step methods the notion of SSP includes all the previous steps considered.

Proof 4.2**.**

We write the method in the form

[TABLE]

the last inequality holds provided that $\Delta t\frac{\beta_{l}}{\alpha_{l}}\leq\Delta t_{\textup{FE}}$. Using the fact that by consistency we have $\sum_{l=1}^{k}\alpha_{l}=1$, we can conclude that

[TABLE]

However, SSP multi-step methods require many steps for high order and have small SSP coefficients as shown in [9]. Now that we see that building explicit SSP multi-step integrating factor methods is trivial, we are ready to approach the main question of this paper: is it possible to extend our results in [15] to two-step Runge–Kutta methods?

4.2 SSP integrating factor two-step Runge–Kutta methods

Any explicit SSP TSRK method of the form (18) that is SSP with SSP coefficient $\mathcal{C}=r$ can be re-written in the matrix-vector form (19), which can also be represented as

[TABLE]

where $0\leq v_{i}\leq 1$, $\hat{\alpha}_{i}\geq 0$, and $\alpha_{i,j}\geq 0$. Here, the consistency condition $\mathbf{R}\mathbf{e}+\mathbf{P}\mathbf{e}=\mathbf{e}$ can be written as $v_{i,1}+v_{i,2}+\hat{\alpha}_{i}+\sum_{j=1}^{i-1}\alpha_{i,j}=1$ for each row $i$.

We note that the time-levels for each stage, also called the abscissas,

[TABLE]

can be computed (in the matrix form) by $\mathbf{D}$, $\mathbf{A}$, and $\hat{\mathbf{A}}$ where $\mathbf{c}=\hat{\mathbf{A}}\mathbf{e}+\mathbf{A}\mathbf{e}-\mathbf{D}_{i,1}$. The abscissa corresponding to $u^{n}$ is [math], the ones corresponding to all the internal stages $j=2,...,s$ are $c_{j}$ and the one corresponding to the final stage $u^{n+1}$ is $1$.

Applying (26) to the transformed equation becomes

[TABLE]

The following theorem tells us what conditions an explicit SSP TSRK method (26) has to satisfy to preserve the strong stability of a numerical solution when coupled with an integrating factor approach as in (27).

Theorem 4.3**.**

Given a linear operator $L$ that satisfies (14) and a (possibly nonlinear) operator $N(u)$ that satisfies (8) for some value of ${\Delta t}_{FE}>0$, and an explicit SSP two-step Runge–Kutta integrating factor method of the form (27) based on the SSP method (26) with non-decreasing abscissas $0=c_{1}\leq c_{2}\leq c_{3}.......\leq c_{s}\leq c_{s+1}=1$, then the numerical solution satisfies

[TABLE]

*where the SSP coefficient $\mathcal{C}=r$. *

Proof 4.4**.**

Consider each stage of (27)

[TABLE]

Now noting that $c_{i}-c_{j}\geq 0$ and applying the result in Corollary 1, we obtain

[TABLE]

*for any $\Delta t\leq\mathcal{C}{\Delta t}_{FE}$. *

This theorem shows us that to find optimal SSP two-step Runge–Kutta methods that are suitable for use with the integrating factor approach, we must augment the optimization problem with additional conditions reflecting the need for non-decreasing abscissas.

To see how this theorem matters in practice, consider, once again, the motivating example in Subsection 2.1 where we solve Equation (11) using an integrating factor approach. In this case, we use $a=2$ so $Lu\approx-2u_{x}$. Otherwise, everything is the same as in Subsection 2.1. We evolve the transformed problem forward 25 time-steps for various time-steps using the eSSP-TSRK$(10,5)$ method in [22] and a new eSSP-TSRK${}^{+}(10,5)$ method that has non-decreasing abscissas. The $\log_{10}$ of the maximal rise in total variation at each stage vs. the value $\lambda=\frac{\Delta t}{\Delta x}$ is shown in Figure 2. We observe that the new method preserves the TVD behavior of the spatial discretization up to a large $\lambda$, while the eSSP-TSRK$(10,5)$ method in [22] does not preserve the TVD property of the spatial discretization.

5 Optimal and optimized methods with non-decreasing abscissas

We augmented the optimization problem from [22] described in Subsection 3.1 by imposing the non-decreasing abscissa constraint

[TABLE]

in addition to (22). We then implemented in Matlab (as in [18, 20, 19]), and used to find optimized explicit eSSP-TSRK+ methods of $s\leq 10$ stages and order $p\leq 7$, as well as an eighth order method with $s=11$. The new methods have non-decreasing abscissas, and are denoted eSSP-TSRK${}^{+}(s,p)$. As we saw above, these methods preserve the SSP properties of a transformed problem (10) and so can be used as a basis for explicit SSP integrating factor two-step Runge–Kutta (SSPIF-TSRK) methods.

The SSP coefficients and effective SSP coefficients of the optimized methods with non-decreasing abscissas (eSSP-TSRK+) are listed in Tables 3 and 4. We boldface the coefficients of the eSSP-TSRK+ methods that are the same as the optimized eSSP-TSRK methods (previously found in [22, 23]) shown in Tables 1 and 2. In Figure 3 we show the effective SSP coefficients of the eSSP-TSRK+ methods and compare them to those of the corresponding eSSP-TSRK methods. For orders $p\leq 4$ we also show the effective SSP coefficients of the eSSPRK+, for comparison. These tables and graphs include all the second order $p=2$ methods, which already have non-decreasing coefficients.

The third order eSSP-TSRK${}^{+}(s,3)$ methods with $s=2,3,4$ also have the same effective SSP coefficients as the eSSP-TSRK$(s,3)$ methods in [22]. Although the eSSP-TSRK$(3,3)$ method in [22] has decreasing coefficients, we were able to find an eSSP-TSRK${}^{+}(3,3)$ non-decreasing coefficient method with the same SSP coefficient. The eSSP-TSRK+(s,3) methods with $s>4$ have smaller SSP coefficients than the corresponding eSSP-TSRK methods, but the loss in size of effective SSP coefficient is under $1$%. These methods have a significant advantage over the eSSPRK+ methods, but this advantage disappears as we allow a larger number of stages $s$, as can be seen in Figure 3.

The fourth order eSSP-TSRK${}^{+}(s,4)$ methods with $s=3,4$ also have the same SSP coefficients as the corresponding eSSP-TSRK in [22]. These two methods are significant because allowing an additional step makes possible fourth order SSP methods with $s=3,4$, which cannot be achieved for eSSPRK methods. Once we look at more stages, $s>4$, we see that the eSSP-TSRK${}^{+}(s,4)$ methods suffer some loss over the eSSP-TSRK$(s,4)$ in the magnitude of the effective SSP coefficient, but it is still under $2$%, so the differences are not obvious in Figure 3. Once again, these eSSP-TSRK${}^{+}(s,4)$ methods offer a significantly larger effective SSP coefficient than the corresponding eSSPRK+ methods.

The fifth order eSSP-TSRK${}^{+}(s,5)$ methods with $s=4,5,6,9$ also have the same SSP coefficients as the corresponding eSSP-TSRK methods in [22]. We note that our optimization code found a slightly better eSSP-TSRK${}^{+}(9,5)$ than the eSSP-TSRK$(9,5)$ previously found in the literature. For $s=7,8,10$ the eSSP-TSRK${}^{+}(s,5)$ cannot match the SSP coefficients of the eSSP-TSRK$(s,5)$ methods in [22]. The loss that we see in the effective SSP coefficients is, once again, under $1$%, so the difference cannot be seen in Figure 3. It is important to note that there are no explicit SSP Runge–Kutta methods of fifth order or above, so the eSSP-TSRK methods are particularly useful if one needs to go to higher order.

The sixth order eSSP-TSRK${}^{+}(s,6)$ methods we found all have smaller SSP coefficients than the eSSP-TSRK$(s,6)$ in [22]. For $s\leq 7$ the loss is under $0.5$% which cannot be seen in Figure 3. However, once we get to $s=8,9,10$ we see more significant losses in the size of the SSP coefficients, however these are still below $5$%. It is interesting to notice that the explicit eSSP-TSRK${}^{+}(10,6)$ method has an effective SSP coefficient $\mathcal{C}_{\textup{eff}}\approx 0.31$ that is about half that of the eSSPRK$(10,4)$ (which is $\mathcal{C}_{\textup{eff}}=0.6$) and about $60$% of that of the eSSPRK${}^{+}(10,4)$ (which is $\mathcal{C}_{\textup{eff}}\approx 0.5299$).

The eSSP-TSRK${}^{+}(8,7)$ method we found has the same SSP coefficient as the eSSP-TSRK$(8,7)$ in [22], while the eSSP-TSRK${}^{+}(s,7)$ methods with $s=9,10$ have a loss in the effective SSP coefficient of $4.4$% and $5.55$% compared with the corresponding eSSP-TSRK methods in [22]. Although it is possible to find eighth order methods, these require more than ten stages. Using our optimization code, we found an eleven stage eighth order method with non-decreasing abscissas eSSP-TSRK+(11,8) with an effective SSP coefficient $\mathcal{C}_{\textup{eff}}\approx 0.0249$. The SSP coefficient of the eSSP-TSRK${}^{+}(11,8)$ methods is significantly smaller (by close to 20%) than the eleven stage eighth order method eSSP-TSRK$(11,8)$ in [22] which has $\mathcal{C}_{\textup{eff}}\approx 0.031$.

In general, we see that for lower number of stages, the SSP coefficients of the methods with non-decreasing abscissas are no smaller than the methods with decreasing abscissas. Perhaps this reflects the fact that the solution space for methods with $s$ close to $p$ is already limited, and so the additional constraint of non-decreasing abscissas poses no major restriction. However, it is interesting to notice that despite the constrained solution space when $s$ is close to $p$, in every case where an eSSP-TSRK$(s,p)$ method was found it was also possible to find an eSSP-TSRK${}^{+}(s,p)$ method.

6 Numerical Results

In this section we demonstrate the performance of the SSP integrating factor two-step Runge–Kutta (SSPIF-TSRK) methods based on the eSSP-TSRK+ methods reported in Section 5. The SSP coefficients of these methods are presented in Table 3; the coefficients of these methods can be downloaded from [8].

Our tests focus on convergence and the SSP properties of the methods. In Subsection 6.1 we verify convergence at the design order for a selection of these methods on the van der Pol problem. In Subsections 6.2 and 6.3 we study the behavior of these methods in terms of their allowable TVD time-step on linear and nonlinear problems with spatial discretizations that are TVD. The simple TVD test in these sections have been used extensively because they demonstrate the behavior of the SSP time-step, and provide evidence that there are some cases in which the SSP property is necessary to preserve the TVD property.

6.1 Example 1: Convergence study

Consider the van der Pol problem, a nonlinear system of ODEs of the form

[TABLE]

The problem can be written as $\mathbf{u}_{t}=L\mathbf{u}+N(\mathbf{u})$ where $\mathbf{u}=(u_{1};u_{2})$ and

[TABLE]

We initialize the problem with $\mathbf{u}_{0}=(2;0)$. Using a variety of SSPIF-TSRK methods, we run the problem to final time $T_{final}=2$, with $\Delta t=0.01,0.02,0.04,0.05,0.08,0.10$. The initial values and the exact solution (for error calculation) was calculated by Matlab’s ODE45 routine with tolerances set to AbsTol=$2\times 10^{-14}$ and RelTol=$2\times 10^{-14}$. We tested the SSPIF-TSRK methods based on the eSSP-TSRK+ methods whose SSP coefficients are presented in Table 3 above. We calculated the orders by finding the slopes of the $\log_{10}(errors)$ using Matlab’s polyfit function. The results, shown in Figure 4 for several selected methods, validate that the methods exhibit the expected order of convergence.

6.2 Example 2: Sharpness of SSP time-step for a linear problem

While the importance of SSP methods is not limited to cases where we need to preserve the TVD property, TVD studies tend to exhibit the sharpness of the SSP time-step and so are traditionally used to test SSP time-stepping methods. We begin with a case where the TVD property of the spatial discretization can be proven, and investigate the TVD behavior of the numerical solutions resulting from evolving the semi-discretized problem using our SSP integrating factor two-step Runge–Kutta methods.

We use our SSPIF-TSRK methods to evolve in time the linear advection equation on the domain $x\in[0,1]$ with periodic boundary conditions and a step function initial condition

[TABLE]

Using a simple first-order forward difference to discretize the spatial derivatives with a spatial grid of $1000$ points. We split the problem as with $Lu\approx-au_{x}$ and $N(u)\approx-u_{x}$. The spatial discretization is TVD when coupled with forward Euler in the sense that

[TABLE]

and

[TABLE]

We evolve the numerical solution forward ten time-steps for different values of $\lambda=\frac{\Delta t}{\Delta x}$, and measure the total variation at each stage. To ensure internal stage monotonicity we compare the total variation at any stage to the total variation at the previous stage. We are interested in the observed TVD time-step $\Delta t^{TVD}_{obs}$, which is the size of time-step $\Delta t$ at which the maximal rise in total variation between any two stages is greater than $10^{-12}$. We compare the observed TVD time-step with the theoretical TVD time-step condition. The SSP coefficient that corresponds to the value of the observed TVD time-step is called the observed SSP coefficient $\mathcal{C}_{obs}$, defined by

[TABLE]

Since the forward Euler TVD time-step for this problem is $\Delta t_{\textup{FE}}=\Delta x$, we have

[TABLE]

6.2.1 Example 2a: Comparison of integrating methods for wavespeeds $\boldsymbol{a=5}$ and $\boldsymbol{a=2}$

In this example we consider a variety of SSPIF-TSRK methods for the time-evolution of (30) with wavespeeds $a=5$ and $a=2$.

In Figure 5 (left) we show the observed maximal rise in total variation when this equation is evolved forward by several SSPIF-TSRK methods using different values of the CFL number $\lambda=\frac{\Delta t}{\Delta x}$. We consider a selection of fourth order methods. First, we explore the performance of the SSPIF-TSRK$(s,p)$ methods with $(s,p)=(3,4),(4,4)$ (the corresponding SSPIFRK methods do not exist). If we consider the effective observed time-step

[TABLE]

(i.e. the observed time-step normalized by the number of stages $s$) we see that the SSPIFRK$(9,4)$ has $\mathcal{C}_{obs}^{eff}=0.47$ while the SSPIF-TSRK$(3,4)$ has a smaller $\mathcal{C}_{obs}^{eff}=0.42$ and SSPIF-TSRK$(4,4)$ has an even smaller $\mathcal{C}_{obs}^{eff}=0.4$. However, we observe that the allowable TVD time-step of the SSPIFRK methods with $(s,p)=(5,4),(9,4)$ is smaller than that of the corresponding SSPIF-TSRK methods. The SSPIF-TSRK methods allow us to obtain higher order for fewer stages than the SSPIFRK methods, and larger observed time-steps for the same $(s,p)$.

Figure 5 (right) shows us the performance of the higher order SSPIF-TSRK methods, with orders $p=5,6,7,8$. In this case we use $a=2$ in (30). We selected the best performing methods in terms of observed SSP time-step, the SSPIF-TSRK$(s,p)$ methods for $(s,p)=(10,5),(10,6),(8,7),(11,8)$. We note that the $(s,p)=(11,8)$ method has the same SSP time-step performance as the $(s,p)=(10,6)$, showing the eighth order method to be a viable and efficient method.

6.2.2 Example 2b: Considering different wavespeeds

In [15] we showed that unlike fully explicit SSP Runge–Kutta methods or implicit-explicit SSP Runge–Kutta (SSP-IMEX) methods, the SSP integrating factor Runge–Kutta schemes have an allowable SSP time-step that is not impacted by the restriction on the time-step resulting from the linear component $L$. This is also true for the SSP integrating factor two-step Runge–Kutta methods. In this section, we investigate how different wavespeeds in $L$ impact the allowable SSP time-step of different methods.

In Table 5 we show the observed SSP coefficient $\lambda^{TVD}_{obs}$ for values of $a=0,1,5$ in (30) for a variety of SSPIF-TSRK methods. We notice that for the methods with $(s,p)=(5,4),(9,4),(6,5),(9,5)$ the observed SSP coefficient for $a=0,1$ was sharp, i.e. exactly the SSP coefficient $\mathcal{C}$ predicted by the theory. As $a$ gets larger ($a=5$) the observed SSP coefficient becomes larger than predicted by the theory, which only provides a guaranteed lower bound. This occurs because the damping effect of the exponential increases as $a$ gets larger and so rises in total variation are damped out and fall below the threshold. For the other methods, the observed SSP coefficient for $a=0$ is larger than predicted by the theory. For $(s,p)=(7.6),(9,6)$ the observed SSP coefficient for $a=0,1$ was the same, while for $a=5$ the observed SSP coefficient is larger. For the other methods, we see an increase in the observed SSP coefficient as $a$ increases.

These results are in sharp distinction to explicit SSP Runge–Kutta methods and SSP-IMEX methods. In Table 6 we compare the observed SSP coefficient of some four stage explicit SSP integrating factor methods to those of the explicit SSP Runge–Kutta method eSSPRK(4,3) and the SSP-IMEX(4,3,K) methods.

When the eSSPRK(4,3) method is applied to (30) with wavespeed $a+1$, where $a=0,1,2,10$, Table 6 shows that the observed SSP coefficient exactly matches the predicted

[TABLE]

This means that as the wavespeed $a$ increases, the observed SSP coefficient goes down as expected.

The SSP-IMEX(4,3,K) we use are from [4], and have SSP explicit and implicit parts that were optimized for the SSP step size for each value of $K=\frac{1}{a}$. The results presented in Table 6 show that, as expected from SSP theory, the observed SSP coefficient decays linearly as the wavespeed $a$ rises. This matches with the results in [4] that, in contrast to the fact that IMEX methods may be A-stable, they cannot be unconditionally SSP past first order. This limitation highlights the need for SSP integrating factor methods.

Finally, we see that the observed SSP coefficients for the explicit SSP integrating factor methods do not decay as the wavespeed $a$ goes up. For the SSPIFRK(4,3) method the observed SSP coefficient is exactly the same for $a=0,1,2,10$, while the SSPIF-TSRK(4,3) and SSPIF-TSRK(4,4) methods have the same observed SSP coefficient for $a=0,1,2$ and a larger observed SSP coefficient for $a=10$.

6.3 Example 3: Sharpness of SSP time-step for a nonlinear problem

Once again we consider the motivating problem (11) in Subsection 2.1:

[TABLE]

on the domain $[0,1]$ with periodic boundary conditions.

We used a first-order upwind difference to semi-discretize, with $400$ spatial points, the linear term $Lu\approx-au_{x}$, and a fifth order WENO finite difference for the nonlinear terms $N(u)\approx-\left(\frac{1}{2}u^{2}\right)_{x}$. We note that the WENO scheme is not TVD. However, its total variation is generally well-controlled in typical simulations.

We evolve the semi-discrerized problem forward 25 time steps using $\Delta t=\lambda\Delta x$, and measure the total variation at each stage. We then calculate the maximal rise in total variation over each stage. Once again, the quantity of interest is the observed TVD time-step, but in the absence of a forward Euler condition (as WENO is not TVD) we cannot determine the observed SSP coefficient, only the observed TVD time-step.

We first show the results for a selection of third order methods on problem (11) with the value $a=5$. In Figure 6 (left) we plot $\log_{10}$ of the maximal rise in total variation versus the ratio $\lambda=\frac{\Delta t}{\Delta x}$. Our SSPIF-TSRK(3,3) method out-performs the SSPIFRK(3,3) method in [15] by more than 10%. In contrast, the fully explicit three stage third order Shu-Osher method applied begins to feature a large rise in total variation for a much smaller value of $\lambda=.15\approx\frac{1}{1+a}$, as expected. Evolving the transformed problem with the Shu-Osher method (2.1) (denoted IF with eSSPRK(3,3)) and with the three-stage third order Cox and Matthews exponential time differencing Runge–Kutta method (ETDRK3) [5] results in a maximal rise in total variation that increases rapidly with $\lambda$.

In Figure 6 (right) we show a similar study using fourth order methods. Our SSPIF-TSRK(5,4) method allows a TVD time-step of $\lambda_{obs}\approx 1.3$, while the corresponding integrating factor Runge–Kutta method SSPIFRK(5,4) method maintains a small maximal rise in total variation until $\lambda_{obs}\approx 1$. The fourth order ETDRK4 method of [5] does not have good TVD performance: the maximal rise in total variation when using the ETDRK4 rises rapidly from the smallest value of $\lambda$. The SSP-IMEX(5,4,0.1) method features an observed value of $\lambda_{obs}\approx 0.5$, pointing once again to the fact that IMEX methods are not as suitable when the SSP property is desired.

Finally, we compare TVD time-step of the SSPIF-TSRK$(s,p)$ methods with $(s,p)=(4,3),(5,4),(6,5),(7,6),(8,7)$. In Figure 7 we show the $\log_{10}$ of the maximal rise in total variation at each stage compared to the previous stage, plotted against the ratio $\lambda=\frac{\Delta t}{\Delta x}$. These integrating factor methods all perform well for this problem. Notably, the seventh order method has the largest observed TVD time-step and the sixth order method has the smallest observed TVD time-step for both $a=1$ and $a=5$. The fourth order method behaves consistently for both $a=1$ and $a=5$. These results are all consistent with those of the linear test case. Once again, we observe that as the exponent increases we generally have a larger allowable TVD time-step, as we noted above.

7 Conclusions

In this work we extended the SSP theory for integrating factor Runge–Kutta methods established in [15] to integrating factor two step Runge–Kutta methods. SSP integrating factor methods are needed for problems that require a nonlinear non-inner product stability property, and where the linear component is severely restricting the time-step. Unlike the case where linear $L_{2}$ stability is of interest, implicit or IMEX methods do not alleviate this restriction. However, integrating factor Runge–Kutta methods based on applying a SSP Runge–Kutta method with nondecreasing abscissas to the transformed equation completely alleviate the time-step restriction. SSP Runge–Kutta methods only exist up to fourth order. To extend these results to higher order methods, we considered herein SSP two-step Runge–Kutta methods applied in conjunction with an integrating factor approach. We showed that, just as for the Runge–Kutta case, if the two-step Runge–Kutta methods are SSP and have non-decreasing abscissas then they can be used to time-step the transformed problem. We enhanced the optimization problem in [22] with the condition (29) and found SSP two-step Runge–Kutta methods with non-decreasing abscissas that have optimized SSP coefficients (presented in Table 3). We tested the resulting SSPIF-TSRK methods on a variety of problems and showed that the new high order methods perform as expected in terms of convergence and the strong stability property.

Appendix A Order Conditions

The order conditions, (22c) for a $k$-step SSP Runge-Kutta were derived by Ketcheson and used in [23]. We can write the generalized form for a multi-step Runge–Kutta method as

[TABLE]

and rewrite it into a matrix form by defining

[TABLE]

The method then becomes

[TABLE]

and by leting $\boldsymbol{l}$ be the vector $\boldsymbol{l}=\left(k-1,k-2,...,1,0\right)^{T}$ to compute the abscissas $\mathbf{c}=\tilde{\mathbf{A}}\mathbf{e}-\tilde{\mathbf{D}}\boldsymbol{l}$ to get the the following expressions

[TABLE]

Below are the order conditions for methods from first to eighth order. As the order increases the methods must satisfy all previous order conditions as well as the additional ones for each order.

First:

[TABLE]

Second:

[TABLE]

Third:

[TABLE]

Fourth:

[TABLE]

Fifth:

[TABLE]

Sixth:

[TABLE]

Seventh:

[TABLE]

Eighth:

[TABLE]

Where $\mathbf{C}$ is a diagonal matrix of the abscissas $\mathbf{C}=diag(\mathbf{c})$ and the exponentiation is considered to be element wise i.e. $\mathbf{c}^{3}=\left(\mathbf{c}\cdot\mathbf{c}\cdot\mathbf{c}\right)$.

Acknowledgment. This publication is based on work supported by AFOSR grant FA9550-18-1-0383. 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.