Implicit and Implicit-Explicit Strong Stability Preserving Runge-Kutta Methods with High Linear Order

TL;DR

This paper develops high linear order implicit and implicit-explicit SSP Runge-Kutta methods that optimize stability properties for solving hyperbolic PDEs, extending the order beyond previous limitations for implicit methods.

Contribution

It introduces new high linear order implicit and IMEX SSP Runge-Kutta methods with optimized stability regions, surpassing previous order restrictions for implicit SSP methods.

Findings

Implicit SSP Runge-Kutta methods of any linear order were constructed.

Optimized IMEX SSP methods with high linear and nonlinear orders were developed.

Numerical tests confirmed the high order convergence and stability properties.

Abstract

When evolving in time the solution of a hyperbolic partial differential equation, it is often desirable to use high order strong stability preserving (SSP) time discretizations. These time discretizations preserve the monotonicity properties satisfied by the spatial discretization when coupled with the first order forward Euler, under a certain time-step restriction. While the allowable time-step depends on both the spatial and temporal discretizations, the contribution of the temporal discretization can be isolated by taking the ratio of the allowable time-step of the high order method to the forward Euler time-step. This ratio is called the strong stability coefficient. The search for high order strong stability time-stepping methods with high order and large allowable time-step had been an active area of research. It is known that implicit SSP Runge-Kutta methods exist only up to sixth order. However, if we restrict ourselves to solving only linear autonomous problems, the order conditions simplify and we can find implicit SSP Runge-Kutta methods of any linear order. In the current work we aim to find very high linear order implicit SSP Runge-Kutta methods that are optimal in terms of allowable time-step. Next, we formulate an optimization problem for implicit-explicit (IMEX) SSP Runge-Kutta methods and find implicit methods with large linear stability regions that pair with known explicit SSP Runge-Kutta methods of orders plin=3,4,6 as well as optimized IMEX SSP Runge-Kutta pairs that have high linear order and nonlinear orders p=2,3,4. These methods are then tested on sample problems to verify order of convergence and to demonstrate the sharpness of the SSP coefficient and the typical behavior of these methods on test problems.

Figures

Tables

Table 1. Table 1: SSP coefficients of optimal explicit and implicit Runge–Kutta methods and the ratio of the implicit coefficient to the explicit coefficient. A dash indicates that SSP methods of this type cannot exist, a blank space indicates none were found.

	Explicit Methods			Implicit Methods					Ratio of Im/Ex
$s$ $\backslash$ $p$	2	3	4	2	3	4	5	6	2	3	4
1	-	-	-	2	-	-	-	-
2	1	-	-	4	2.73	-	-	-	4
3	2	1	-	6	4.83	2.05	-	-	3	4.83
4	3	2	-	8	6.87	4.42	1.14		2.67	3.44
5	4	2.65	1.51	10	8.90	6.04	3.19		2.50	3.36	4.00
6	5	3.54	2.28	12	10.92	7.80	4.97	0.18	2.40	3.08	3.42
7	6	4.27	3.29	14	12.93	9.19	6.21	0.26	2.33	3.03	2.79
8	7	5.12	4.15	16	14.94	10.67	7.56	2.25	2.29	2.92	2.57
9	8	6.03	4.86	18	16.94	12.04	8.90	5.80	2.25	2.81	2.48
10	9	6.80	6.00	20	18.95	13.64	10.13	8.10	2.22	2.79	2.27
11	10	7.59	6.50	22	20.95	15.18	11.33	8.85	2.20	2.76	2.34

Table 2. Table 2: SSP coefficients of LNL explicit Runge–Kutta methods from [ 10 ] . A dash indicates that SSP methods of this type cannot exist, and the boldface indicates that the SSP coefficients are equal to the corresponding ones of lower order.

	Explicit $p=2,p=3$ methods
$s\backslash p_{lin}$	5	6	7	8	9
5	1.00	–	–	–	–
6	2.00	1.00	–	–	–
7	2.65	2	1.00	–	—
8	3.37	2.65	2.00	1.00	–
9	4.1	3.37	2.65	2.00	1.00

Table 3. Table 3: The value of 𝒞 𝒞 {\cal{C}} for LNL iRK with linear order p l ​ i ​ n subscript 𝑝 𝑙 𝑖 𝑛 p_{lin} , and nonlinear order p = 2 𝑝 2 p=2 on left and p = 4 𝑝 4 p=4 on right. Note along the diagonal these the SSP coefficients are close to six times larger than the corresponding explicit LNL methods.

	Implicit $p=2,p=3$ methods						Implicit $p=4$ methods
$s\backslash p_{lin}$	4	5	6	7	8	9	4	5	6	7	8	9
3	3.24	–	–	–	–	–	2.05	-	-	-	-	-
4	4.56	3.64	–	–	–	–	4.42	3.54	–	–	–	–
5	6.15	5.01	3.97	–	–	–	6.04	4.63	3.81	–	–	–
6	7.85	6.66	5.17	4.27	–	–	7.80	6.49	5.14	4.14	–	–
7	9.60	8.42	6.54	5.51	4.54	–	9.19	7.86	6.42	5.42	4.46	–
8	11.23	10.21	7.92	6.91	5.82	4.79	10.67	9.25	7.86	6.82	5.80	4.74
9	12.81	11.82	9.48	8.33	7.03	6.10	12.04	11.15	9.46	8.32	6.98	6.08

Table 4. Table 4: The difference between the theoretical and observed SSP coefficient at which the maximal rise in total variation is about the threshold of 10 − 10 superscript 10 10 10^{-10} for the p = 2 𝑝 2 p=2 methods.

$s\backslash p_{lin}$	2	3	4	5	6	7	8	9
2	2.5$\times {10}^{-12}$	1.5$\times {10}^{-12}$
3	4.0$\times {10}^{-12}$	2.5$\times {10}^{-12}$	3.8$\times {10}^{-2}$
4	8.1$\times {10}^{-12}$	4.5$\times {10}^{-12}$	2.2$\times {10}^{-11}$	6.4$\times {10}^{-2}$
5	1.6$\times {10}^{-11}$	9.1$\times {10}^{-12}$	1.2$\times {10}^{-11}$	7.6$\times {10}^{-12}$	1.7 $\times {10}^{-1}$
6	3.2$\times {10}^{-11}$	1.8$\times {10}^{-11}$	1.1$\times {10}^{-11}$	6.6$\times {10}^{-12}$	2.6$\times {10}^{-11}$	3.0$\times {10}^{-1}$
7	6.4$\times {10}^{-11}$	3.7$\times {10}^{-11}$	1.4$\times {10}^{-11}$	7.6$\times {10}^{-12}$	2.6$\times {10}^{-10}$	7.1$\times {10}^{-12}$	3.6$\times {10}^{-1}$
8	1.2$\times {10}^{-10}$	7.4$\times {10}^{-11}$	3.7$\times {10}^{-10}$	1.0$\times {10}^{-11}$	5.5$\times {10}^{-11}$	5.5$\times {10}^{-11}$	1.4$\times {10}^{-2}$	4.7$\times {10}^{-1}$
9	2.5$\times {10}^{-10}$	1.5$\times {10}^{-10}$	1.4$\times {10}^{-10}$	1.7$\times {10}^{-10}$	3.4$\times {10}^{-11}$	2.1$\times {10}^{-11}$	2.2$\times {10}^{-10}$	3.4$\times {10}^{-2}$

Table 5. Table 5: the SSP coefficient of optimal SSP IMEX methods for nonlinear order p = p e = p i = 2 , 3 , 4 formulae-sequence 𝑝 superscript 𝑝 𝑒 superscript 𝑝 𝑖 2 3 4 p=p^{e}=p^{i}=2,3,4 with K = ∞ 𝐾 K=\infty . These match with the SSP coefficients of the explicit methods previously reported in [ 11 , 10 ] , which verifies that the coupling conditions do not adversely affect the SSP coefficient in the case K = ∞ 𝐾 K=\infty . Note that the SSP coefficients in bold correspond to methods that have negative values in 𝐒 ~ ~ 𝐒 \tilde{\mathbf{S}} , while the other methods have only nonnegative coefficient matrices.

nonlinear order	stages	$p_{lin}=3$	$p_{lin}=4$	$p_{lin}=5$	$p_{lin}=6$	$p_{lin}=7$
$p=2,3$	$s=p_{lin}$	1.0000	1.0000	1.0000	1.0000	1.0000
$p=2,3$	$s=p_{lin}+1$	2.0000	2.0000	2.0000	2.0000	2.000
$p=2,3$	$s=p_{lin}+2$	2.6505	2.6505	2.6505	2.6505	2.6505
$p=4$	$s=p_{lin}$	—	—	0.7603	0.8677	1.0000
$p=4$	$s=p_{lin}+1$	—	1.5082	1.8091	1.8269	1.9293
$p=4$	$s=p_{lin}+2$	—	2.2945	2.5753	2.5629	2.6192

Table 6. Table 6: SSP coefficients of some optimized IMEX methods with s = p l ​ i ​ n 𝑠 subscript 𝑝 𝑙 𝑖 𝑛 s=p_{lin} and s = p l ​ i ​ n + 1 𝑠 subscript 𝑝 𝑙 𝑖 𝑛 1 s=p_{lin}+1 stages optimized for different values of K 𝐾 K .

$s=p_{lin}$
$\frac{1}{K}$	p	$p_{lin}$
		3	4	5	6	7
10	2	2.030$\times {10}^{-1}$	1.729$\times {10}^{-1}$	1.520$\times {10}^{-1}$	1.323$\times {10}^{-1}$	1.109$\times {10}^{-1}$
	3	1.492$\times {10}^{-1}$	1.727$\times {10}^{-1}$	1.520 $\times {10}^{-1}$	1.323$\times {10}^{-1}$	0.932$\times {10}^{-1}$
	4	–	–	–	1.132$\times {10}^{-1}$	0.849$\times {10}^{-1}$
100	2	2.23$\times {10}^{-2}$	2.06$\times {10}^{-2}$	1.58$\times {10}^{-2}$	1.51$\times {10}^{-2}$	0.99$\times {10}^{-2}$
	3	1.63$\times {10}^{-2}$	2.05$\times {10}^{-2}$	1.58$\times {10}^{-2}$	1.51$\times {10}^{-2}$	0.99$\times {10}^{-2}$
	4	–	–	–	1.18$\times {10}^{-2}$	0.89$\times {10}^{-2}$

Table 7. Table 7: Comparison of the theoretical and observed allowable time-step before an SSP violation of 10 − 12 superscript 10 12 10^{-12} occurs for the linear problem in Example 4.1 with wavespeed ω = 100 𝜔 100 \omega=100 and for method optimized for the value of K 𝐾 K above.

Type	p	$p_{lin}$	s	$\frac{1}{K}$	Predicted $\Delta t$	Observed $\Delta t$	ratio
Explicit	2	5	6	–	1.98$\times {10}^{-2}$	1.98$\times {10}^{-2}$	–
IMEX				1	1.25$\times {10}^{-2}$	3.07$\times {10}^{-2}$	1.55
IMEX				10	2.89$\times {10}^{-2}$	3.33$\times {10}^{-2}$	1.68
IMEX				100	3.26$\times {10}^{-2}$	3.37$\times {10}^{-2}$	1.70
Implicit				–	6.59$\times {10}^{-2}$	6.59$\times {10}^{-2}$	3.32
Explicit	2	6	6	–	9.90$\times {10}^{-3}$	9.90$\times {10}^{-3}$	–
IMEX				1	6.00$\times {10}^{-3}$	1.47$\times {10}^{-2}$	1.48
IMEX				10	1.32$\times {10}^{-2}$	1.94$\times {10}^{-2}$	1.96
IMEX				100	1.51$\times {10}^{-2}$	1.97$\times {10}^{-2}$	1.99
Implicit				–	5.12$\times {10}^{-2}$	5.12$\times {10}^{-2}$	5.17
Explicit	3	5	5	–	9.90$\times {10}^{-3}$	9.90$\times {10}^{-3}$	–
IMEX				1	6.34$\times {10}^{-3}$	2.00$\times {10}^{-2}$	2.02
IMEX				10	1.52$\times {10}^{-2}$	2.24$\times {10}^{-2}$	2.26
IMEX				100	1.58$\times {10}^{-2}$	2.38$\times {10}^{-2}$	2.40
Implicit				–	4.96$\times {10}^{-2}$	4.96$\times {10}^{-2}$	5.01
Explicit	3	5	6	–	1.98$\times {10}^{-2}$	1.98$\times {10}^{-2}$	–
IMEX				1	1.23$\times {10}^{-2}$	3.10$\times {10}^{-2}$	1.56
IMEX				10	2.77$\times {10}^{-2}$	3.30$\times {10}^{-2}$	1.66
IMEX				100	3.20$\times {10}^{-2}$	3.50$\times {10}^{-2}$	1.76
Implicit				–	6.59$\times {10}^{-2}$	6.59$\times {10}^{-2}$	3.32
Explicit	4	6	6	–	8.59$\times {10}^{-3}$	9.90$\times {10}^{-3}$	–
IMEX				1	5.12$\times {10}^{-3}$	1.91$\times {10}^{-2}$	1.93
IMEX				10	1.13$\times {10}^{-2}$	2.14$\times {10}^{-2}$	2.16
IMEX				100	1.18$\times {10}^{-2}$	2.09$\times {10}^{-2}$	2.11
Implicit				–	5.09$\times {10}^{-2}$	5.09$\times {10}^{-2}$	5.92

Table 8. Table 8: Comparison of the theoretical and observed SSP coefficients that preserve the nonlinear stability properties in Example 4.2 with the methods optimized for the value 1 K 1 𝐾 \frac{1}{K} to match the wavespeed ω 𝜔 \omega .

p	$p_{lin}$	s	$\omega =\frac{1}{K}$	Predicted $\Delta t$	Observed $\Delta t$
2	4	5	10	3.099$\times {10}^{-1}$	3.099$\times {10}^{-1}$
			100	3.65$\times {10}^{-2}$	3.65$\times {10}^{-2}$
3	4	5	10	3.084$\times {10}^{-1}$	3.084$\times {10}^{-1}$
			100	3.56$\times {10}^{-2}$	3.56$\times {10}^{-2}$
3	5	5	10	1.520$\times {10}^{-1}$	2.135$\times {10}^{-1}$
			100	1.58$\times {10}^{-2}$	2.39$\times {10}^{-2}$
2	4	6	10	4.088$\times {10}^{-1}$	4.088$\times {10}^{-1}$
			100	4.67$\times {10}^{-2}$	4.67$\times {10}^{-2}$
3	4	6	10	4.171$\times {10}^{-1}$	4.171$\times {10}^{-1}$
			100	4.59$\times {10}^{-2}$	4.59$\times {10}^{-2}$
3	6	7	10	2.29$\times {10}^{-1}$	3.00$\times {10}^{-1}$
			100	2.98$\times {10}^{-2}$	3.50$\times {10}^{-2}$
4	6	6	10	1.13$\times {10}^{-1}$	2.01$\times {10}^{-1}$
			100	1.18$\times {10}^{-2}$	2.10$\times {10}^{-2}$

Equations

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

$$\Vert u^n+\Delta tF(u^n)\Vert \le \Vert u^n\Vert ,0\le \Delta t\le \Delta t_{FE},$$

$$u_t=F(u),$$

$$y^{(i)}$$

$$u^{n+1}$$

$$y^{(i)}=v_i{u}^n+\sum _{j=1}^s{\alpha }_{i,j}\left(y^{(j)}+\Delta t\frac{{\beta }_{i,j}}{{\alpha }_{i,j}}F(y^{(j)})\right),$$

$$u^{n+1}=v_{s+1}{u}^n+\sum _{j=1}^s{\alpha }_{s+1,j}\left(y^{(j)}+\Delta t\frac{{\beta }_{s+1,j}}{{\alpha }_{s+1,j}}F(y^{(j)})\right).$$

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

$$\Vert u^{n+1}\Vert \le v_{s+1}\Vert u^n\Vert +\sum _{j=1}^s{\alpha }_{s+1,j}\Vert y^{(j)}+\Delta t\frac{{\beta }_{s+1,j}}{{\alpha }_{s+1,j}}F(y^{(j})\Vert .$$

$$\Vert y^{(i)}\Vert \le \Vert u^n\Vert \text{mbox}and\Vert u^{n+1}\Vert \le \Vert u^n\Vert \text{mbox}forany\Delta t\le \mathcal{C}\Delta t_{\text{FE}}$$

$$u_t=F(u)+G(u)$$

$$u^{(i)}$$

$$u^{n+1}$$

$$\displaystyle\left(\begin{array}[]{ll}\mathbf{A}&0\\ {\bf b}&0\\ \end{array}\right)\left(\mathbf{I}+r\left(\begin{array}[]{ll}\mathbf{A}&0\\ {\bf b}&0\\ \end{array}\right)\right)^{-1}\geq 0$$

$$\displaystyle\left\|r\left(\begin{array}[]{ll}\mathbf{A}&0\\ {\bf b}&0\\ \end{array}\right)\left(\mathbf{I}+r\left(\begin{array}[]{ll}\mathbf{A}&0\\ {\bf b}&0\\ \end{array}\right)\right)^{-1}\right\|_{\infty}\leq 1$$

$${\tau }_k(\mathbf{A},b)=0\text{mbox}fork=1,...,P,$$

$${\tau }_q^{lin}(\mathbf{A},\mathbf{b})={\mathbf{b}}^T{\mathbf{A}}^{q-2}\mathbf{c}={\mathbf{b}}^T{\mathbf{A}}^{q-1}\mathbf{e}=\frac{1}{q!}\forall q=1,...,p_{lin}.$$

$${\tau }_k(\mathbf{A},\mathbf{b})\text{mbox}for1\le k\le p$$

$${\tau }_q^{lin}(\mathbf{A},\mathbf{b})\text{mbox}forp\le q\le p_{lin}$$

$$y^{(i)}$$

$$u^{n+1}$$

$$u_1^{\prime }=u_2$$

$$u_2^{\prime }=\frac{1}{ϵ}(-u_1+(1-u_1^2)u_2)$$

$$u_t+f(u{)}_x$$

$$U_t+U_x=0,$$

$$U_{0}(x)=\left\{\begin{array}[]{rc}1&-0.1\leq x\leq 0.1\\ 0&otherwise\\ \end{array}\right.$$

$$\underset{0\le n\le N-1}{\max }\left(\Vert u^{n+1}{\Vert }_{TV}-\Vert u^n{\Vert }_{TV}\right).$$

$$U^{\prime }(t)=f(U(t))+g(U(t)),U(t_0)=U_{0}t\ge t_0$$

$$u^{\prime }(t)=F(u)+G(u).$$

$$u^{(i)}$$

$$u^{n+1}$$

$$\Vert u^n+\Delta tF(u^n)\Vert \le \Vert u^n\Vert ,0\le \Delta t\le \Delta t_{FE},$$

$$\Vert u^n+\Delta tG(u^n)\Vert \le \Vert u^n\Vert ,0\le \Delta t\le K\Delta t_{FE}.$$

$$Y=\mathbf{e}{u}^n+\Delta t\mathbf{S}F(Y)+\Delta t\tilde{\mathbf{S}}G(Y).$$

$$\left(I+r\mathbf{S}+\tilde{r}\tilde{\mathbf{S}}\right)Y$$

Full text

1 Strong Stability Preserving Runge–Kutta methods

The exact solution of a hyperbolic conservation law of the form

[TABLE]

frequently develops sharp gradients or discontinuities, which may cause significant difficulties in numerical simulations. For this reason, the development of high order spatial discretizations that can handle such discontinuities or sharp gradients has been an area of active research for the past few decades. Such methods have special nonlinear non-inner-product stability properties that mimic some significant physical properties of the exact solution, such as total variation stability, a maximum principle, or positivity. These properties ensure that when (1) is discretized in space, the spatial discretization $F(u)$ satisfies

[TABLE]

where $u^{n}$ is a discrete approximation to $U$ at time $t^{n}$ and $\|\cdot\|$ is the desired norm, semi-norm, or convex functional. In other words, when the semi-discretized system

[TABLE]

is evolved forward in time using a first order explicit Euler method, the numerical solution satisfies the desired strong stability property, as long as the time-step is suitably limited. In practice we want to use a higher order time integrator instead of the first order forward Euler method (2), but we still need the strong stability property $\|u^{n+1}\|\leq\|u^{n}\|$ to be satisfied, perhaps under a modified time-step restriction.

Some higher order Runge–Kutta methods can be decomposed into convex combinations of forward Euler steps [31], so that any convex functional property satisfied by (2) is automatically preserved, usually under a different time-step restriction. For example, if the $s$-stage implicit Runge–Kutta method is written in the form [32, 11],

[TABLE]

it is obvious that each stage can be written as a linear combination of the solution at the previous step and forward Euler steps of the stages

[TABLE]

If the coefficients $v_{i}$, $\alpha_{i,j}$ and $\beta_{i,j}$ are all non-negative and $\beta_{i,j}$ is zero whenever the corresponding $\alpha_{i,j}$ is zero, the stages decompose into a convex combination of forward Euler steps of the form (2), with each time step replaced by $\frac{\beta_{i,j}}{\alpha_{i,j}}\Delta t$. Each stage is then bounded by

[TABLE]

Recall that $v_{i}+\sum_{j=1}^{s}\alpha_{i,j}=1$ for consistency, and that each $\|y^{(j)}+\Delta t\frac{\beta_{i,j}}{\alpha_{i,j}}F(y^{(j)})\|\leq\|y^{(j)}\|$ for $\frac{\beta_{i,j}}{\alpha_{i,j}}\Delta t\leq\Delta t_{\textup{FE}}$ from (2), so we have

[TABLE]

where ${\cal{C}}=\min_{i,j}\frac{\alpha_{i,j}}{\beta_{i,j}}$. (If any of the $\beta_{i,j}$’s are equal to zero, we consider the corresponding ratios to be infinite.) Any method that can be written in this form with ${\cal{C}}>0$ is called a strong stability preserving (SSP) method.

Strong stability preserving second and third order explicit Runge–Kutta methods [32] and fourth order methods [33, 21] were developed using this convex combination approach. These methods ensure that any strong stability property satisfied by the spatial discretization when using the forward Euler condition (2) is preserved by the higher order strong stability preserving Runge–Kutta method. Furthermore, the convex combination decomposition above ensures that the strong stability property is also satisfied by the intermediate stages in a Runge–Kutta method. This may be desirable in many applications, notably in simulations that require positivity. Clearly, the condition above is sufficient for preservation of strong stability; in [8, 9, 15, 16] it was shown that this condition is necessary, as well.

While the time-step depends on both the spatial and temporal discretizations, we isolate the contribution of the temporal discretization to the time-step restriction by considering the ratio of the allowable time-step of the high order method to the forward Euler time-step. This ratio is called the strong stability preserving coefficient. Using this approach, we view the time-step restriction (7) as a combination of two factors: the forward Euler time-step $\Delta t_{\textup{FE}}$ that comes from the spatial discretization, and the SSP coefficient ${\cal{C}}$ that is a property only of the time-discretization. For this reason, the research on SSP methods focuses on optimizing the allowable time-step $\Delta t\leq{\cal{C}}\Delta t_{\textup{FE}}$ by maximizing the SSP coefficient ${\cal{C}}$ of the method. Among methods of similar types, a more relevant quantity is the effective SSP coefficient ${\cal{C}}_{\textup{eff}}=\frac{{\cal{C}}}{s}$, which takes into account the computational cost of the method at each iteration, defined by the number of stages $s$ (typically also the number of function evaluations per time-step). However, this measure is not applicable when comparing explicit and implicit methods. Table 1 gives SSP coefficients for explicit SSP Runge–Kutta methods of orders $p\leq 4$ and implicit methods of orders $p\leq 6$ and different stages, and the ratio of the SSP methods of corresponding implicit and explicit methods. Explicit Runge–Kutta methods (and in fact all explicit general linear methods) have a bound on the SSP coefficient ${\cal{C}}\leq s$ [11], while all optimal implicit Runge–Kutta methods have been observed to have ${\cal{C}}\ \leq 2s$ [24], which is only twice that of the explicit method. However, an implicit methods of a given number of stages and order typically has an SSP coefficient that is greater than twice that of the corresponding explicit Runge–Kutta method, as shown in Table 1.

In addition to bounds on the effective SSP coefficients, SSP Runge–Kutta methods also suffer from barriers on their order: explicit Runge–Kutta methods with positive SSP coefficient cannot be more than fourth-order accurate [26, 30] and implicit Runge–Kutta methods with positive SSP coefficient cannot be more than sixth-order accurate [24, 11]. These restrictive order barriers of SSP Runge–Kutta methods is a result of the nonlinearity of the ODEs. However, if we are only interested in the order of accuracy on linear autonomous ODE systems, Runge–Kutta methods need only satisfy a smaller set of order conditions. If we only want the method to have high linear order ($p_{lin}$), then the order barrier is broken and Runge–Kutta methods with positive SSP coefficients exist for arbitrarily high linear orders. These methods are of interest because their SSP coefficients serve as upper bounds for nonlinear methods. However, they are also useful in their own right for linear problems where the strong stability preserving property is required, such as Maxwell’s equations and the equations of linear elasticity.

In [10] the observation that the order barrier is only applicable to the nonlinear order led to the study of explicit SSP Runge–Kutta methods with high linear order and optimal nonlinear order $p=4$, which demonstrated that where high linear order is desired, going to higher nonlinear order costs little or nothing in terms of the SSP coefficient. We refer to methods that may have a higher linear order than nonlinear order ($p_{lin}\geq p$) as linear/nonlinear (LNL) methods. Table 2 shows the optimized SSP coefficients of these explicit SSP LNL Runge–Kutta methods. These methods may be of particular value for problems that require high linear order throughout but also have some nonlinear features in some regions which benefit from order greater than two.

The explicit LNL methods found in [10] suggest that implicit methods in this class could be beneficial as well. Implicit SSP Runge–Kutta methods with very high linear order have not been widely studied, and these methods could be very useful for linear problems. In addition, requiring that these methods have a nonlinear order $p>2$ would allow these methods to be more useful if the problem is generally linear but has regions that feature nonlinearities. In this work, we expand upon the explicit methods in [10] and seek optimal implicit SSP methods with very high linear order and nonlinear orders $p\leq 6$. We first consider implicit SSP Runge–Kutta methods that have high order for linear problems and nonlinear order $p=2$. We then consider optimal SSP methods with nonlinear orders $p=3,4,5,6$ and compare the SSP coefficients of these optimal methods. These methods are presented in Section 2. These implicit SSP LNL Runge–Kutta methods have SSP coefficients that are up to six times the size of the corresponding explicit SSP LNL Runge–Kutta methods. In Section 2.4 we verify the convergence of these methods as well as the sharpness of the SSP coefficient on sample problems. Next, in Section 3 we consider implicit-explicit Runge–Kutta methods, for use with problems of the form

[TABLE]

where we wish to treat $F$ explicitly and $G$ implicitly. We formulate an optimization problem that depends on the ratio of the forward Euler condition of the two components $F$ and $G$. Using the optimization routine we developed, we found SSP IMEX Runge–Kutta methods where the explicit component has nonlinear order $p^{e}=3$ and $p^{e}=4$ and higher linear order $p_{lin}>4$, and the implicit component has nonlinear orders $p^{i}=2,3,4$ with higher linear order $p_{lin}$ that matches the linear order of the explicit part. In Section 3.2 we present the optimal methods of this type and in Section 3.3 we verify the convergence of these methods and show the behavior of these methods on sample problems that require the SSP property.

2 Implicit Runge–Kutta with high linear order

2.1 The SSP optimization problem

Our goal is to find implicit SSP Runge–Kutta methods of a given order and number of stages with the largest possible SSP coefficient ${\cal{C}}$. While the most convenient formulation for directly observing the SSP coefficient of a Runge–Kutta method is the Shu-Osher form (4), for formulating the optimization problem it is more convenient to use the Butcher form

[TABLE]

(where the coefficients $a_{ij}$ are place into the matrix $\mathbf{A}$ and $b_{j}$ into the row vector ${\bf b}$) [21]. The reason this approach is preferable is that the Shu-Osher form of a Runge–Kutta method is not unique, while the Butcher form is, so that rather than perform a search for an optimal convex combination of the Shu-Osher form, we seek the unique optimal Butcher coefficients (4).

Following the approach developed by Ketcheson [21],and used in [21, 24, 22, 11, 23, 10] to find optimal SSP methods, we aim to maximize the value of $r$ under the following constraints:

[TABLE]

where $\tau_{k}$ are the order conditions, described below. This optimization gives the Butcher coefficients $\mathbf{A}$ and ${\bf b}$ and an optimal value of the SSP coefficient ${\cal{C}}=r$.

The conversion from the optimal Butcher form to the canonical Shu-Osher form is given in [11], and given the matrix $\mathbf{A}$, the vector ${\bf b}$ and the SSP coefficient $r$, the Shu-Osher coefficients $\alpha$ and $\beta$ can be easily accomplished in MATLAB by

s=size(A,1);
K=[A;b’];
G=eye(s)+r*A;
beta=K/G;
alpha=r*beta;

where the coefficients $v_{i}$ are then computed by the consistency condition $v_{i}+\sum_{j=1}^{s}\alpha_{i,j}=1$.

Some of the major constraints in the optimization problem come from the order conditions $\tau_{k}(\mathbf{A},{\bf b})$ above. For to demonstrate the correct order of accuracy for nonlinear problems, it must satisfy the order conditions:

$\tau_{1}(\mathbf{A},{\bf b})$:   ${\bf b}^{T}{\bf e}=1$

$\tau_{2}(\mathbf{A},{\bf b})$:   ${\bf b}^{T}{\bf c}=\frac{1}{2}$

$\tau_{3}(\mathbf{A},{\bf b})$:   ${\bf b}^{T}\mathbf{A}{\bf c}=\frac{1}{6}$   and   ${\bf b}^{T}{\bf c}^{2}=\frac{1}{3}$

$\tau_{4}(\mathbf{A},{\bf b})$:   ${\bf b}^{T}\mathbf{A}^{2}{\bf c}=\frac{1}{24}$   and   ${\bf b}^{T}\mathbf{A}{\bf c}^{2}=\frac{1}{12}$   and   ${\bf b}^{T}{\bf c}\mathbf{A}{\bf c}=\frac{1}{8}$   and   ${\bf b}^{T}{\bf c}^{3}=\frac{1}{4}$ ,

where ${\bf c}=\mathbf{A}{\bf e}$ and ${\bf e}$ is a vector of ones. Thus, for first order there is only one condition ($P=1$ in Equation (9k) above), for second order two ($P=2$), for third order four ($P=4$), and for fourth order eight conditions are needed ($P=8$). For fifth order, we require a total of $P=17$ conditions, and for sixth order $P=37$ conditions must be satisfied.

It is well-known that there are no implicit SSP Runge–Kutta methods greater than sixth order. In this work, we will consider only implicit methods up to order $p=6$, but we will allow the methods to satisfy higher order $p_{lin}\geq p$ when applied to a linear problem. In this case, the order conditions simplify, and can be expressed as

[TABLE]

In the following sections we give the SSP conditions of the optimized methods resulting from the optimization problem given in (9) with order conditions

[TABLE]

and

[TABLE]

which we implemented in a MATLAB code [5]. The coefficients of these optimized methods can be downloaded from [5].

2.2 Optimal SSP LNL implicit Runge–Kutta methods

We started our search by considering fully implicit Runge–Kutta methods. However, as in [10] these methods had SSP coefficients that were identical to those of the optimized diagonally implicit Runge–Kutta (DIRK) methods, so we proceeded to consider only DIRK methods.

We first found SSP DIRK methods of up to $s=9$ stages with with linear order $2\leq p_{lin}\leq s+1$ and nonlinear order $p=2$. These methods tend to have nice low-storage forms when written in their Shu-Osher arrays. Our first observation is that the methods with nonlinear order $p=2$ and linear order $p_{lin}=3$ the same SSP coefficient as the third order ($p=3$) optimal implicit SSP Runge–Kutta methods listed [24, 11]. Clearly, the additional nonlinear order condition imposed did not constrain the methods any further.

A similar behavior is observed for SSP DIRK methods with $p_{lin}\geq 4$ and nonlinear orders $p=2$ and $p=3$, which has SSP coefficients that are the same up to two decimal places, except for the case of $s=3$ and $p_{lin}=4$ for which there was a difference of $\approx 1.43\times 10^{-2}$, and the $s=4$ and $p_{lin}=5$ where there was a difference of $\approx 2\times 10^{-2}$. The SSP coefficients are given in Table 3 (on left). Finding optimized DIRK SSP Runge–Kutta methods with $p=4$ and $\leq 4p_{lin}\geq 9$, we see that the SSP coefficients in this case are smaller than the methods with nonlinear orders $p=2,3$, but as the number of stages increases this difference diminishes, as can be seen on Table 3 (on right). This is also clear in Figure 2, which plots the SSP coefficients of methods of nonlinear order $p=2$ (and therefore of $p=3$ as well) and $p=4$ with different linear orders $3\leq p_{in}\leq 9$. In this figure we observe that if we compare methods with a given $p_{lin}$, the SSP coefficient of the methods with nonlinear order $p=2$ in does not have significantly higher SSP coefficient than the methods with nonlinear order $p=4$. This would seem to suggest that the linear order is the main constraint on the SSP coefficient, as was the case for the explicit methods in [10]. However, Figure 2, which looks at methods with higher nonlinear orders $p=5$ and $p=6$ tells a different story. In this figure we compare three methods compare three methods with linear order $p_{lin}=5$ and nonlinear orders $p=2,4,5$ (in blue solid, dot-dashed with circle marker, and dashed with + markers) and similarly, three methods with linear order $p_{lin}=6$ and nonlinear orders $p=2,4,6$ (in red solid, dot-dashed with circle marker, and dashed with + markers) We say that in this case, as we go to higher nonlinear orders the SSP coefficient is significantly reduced. For example, a six stage ($s=6$) sixth nonlinear order method has an SSP coefficient ${\cal{C}}=0.18$ whereas the six stage LNL method with nonlinear order $p=4$ and linear order $p_{lin}=6$ has an SSP coefficient ${\cal{C}}=5.14$.

The SSP coefficients of the implicit SSP LNL Runge–Kutta methods in 3 should be compared to those of the explicit SSP LNL Runge–Kutta methods 2. we see that the implicit methods have SSP coefficients that are significantly larger than the corresponding explicit methods: up to six times larger. This suggests that the cost of the implicit solver may be offset by the larger allowable time-step in some cases. While the methods in this section are not necessarily suitable for every application, they are valuable in that they provide an upper bound on the possible SSP coefficient for a given number of stages and linear and nonlinear orders, and demonstrate the effect of increasing the nonlinear order.

2.2.1 Coefficients of selected optimal methods

The coefficients of all the methods listed in the section above can be downloaded as .mat files from [5]. In this section we list, for the user’s convenience, the non-zero coefficients of three methods of particular interest. We use the canonical Shu–Osher form here, where all the forward-Euler steps are of the same size $\frac{\Delta t}{r}$ and $r={\cal{C}}$:

[TABLE]

SSP iRK LNL $p=4$, $s=p_{lin}=6$: One of the most efficient method produced was the $s=p_{lin}=6$ LNL method of nonlinear order $p=4$ that has SSP coefficient ${\cal{C}}=5.138$. The corresponding implicit Runge–Kutta method in [11] which has $s=p=p_{lin}=6$ has a much smaller SSP coefficient of ${\cal{C}}=0.18$. The non-zero coefficients are:

[TABLE]

SSP iRK LNL $p=4$, $s=8$, $p_{lin}=9$: A very high order LNL method with $s=8$, $p_{lin}=9$, and $p=4$ has SSP coefficient ${\cal{C}}=4.735$. This method has the following non-zero coefficients given in the canonical Shu-Osher form (11):

[TABLE]

**SSP iRK LNL $p=2$, $s=10$, $p_{lin}=11$: ** A very high order LNL method with $s=10$, $p_{lin}=11$, and $p=2$ has SSP coefficient ${\cal{C}}=5.2306$. This method has a very low storage form, where the non-zero coefficients are given in the canonical Shu-Osher form (11):

[TABLE]

2.3 Optimizing for Additional Properties

In the section above we optimized the implicit Runge–Kutta methods for the largest possible SSP coefficients. However, in many cases the motivation for using an implicit method is not only the SSP coefficient; this is especially true in cases where the additional time-step allowed for the implicit SSP methods is not enough to offset the cost of the implicit solver needed. If we wish to optimize for other properties, such as linear stability regions, alongside the SSP property, we can add a condition to the inequality constraints in the optimization routine. The methods we found above did not have large linear stability regions, and so we wish to explore whether one can find methods that have a large SSP coefficient as well as large linear stability region. In this section we present our optimized methods using this approach. The resulting methods are suboptimal in terms of SSP coefficients but benefit, as desired, from larger regions of absolute stability. In Figures (4) and (4) we show the linear stability regions of methods found from this co-optimization approach compared to those found in the subsection above.

Figure (4) shows several five stage ($s=5$) methods with $p=3$ and $p_{lin}=4$ resulting from optimization runs that required increasing linear stability regions that included more of the real axis. We show four methods: in blue is the linear stability of the optimal SSP method in the subsection above. This method has SSP coefficient ${\cal{C}}=6.1472$, and crosses the negative real axis at $x=-26.3$. In red is the linear stability region of a method that has SSP coefficient ${\cal{C}}=5.9537$ but allows slightly more of the negative real axis, crossing it at $x=-29.766$. Next, in green is the linear stability of a method with a slightly smaller SSP coefficient of ${\cal{C}}=5.6249$ but which allows much more of the negative real axis, crossing it at $x=-56.247$. Finally, in black is the linear stability region of a method with ${\cal{C}}=5.4459$ but which allows significantly more of the negative real axis, crossing it at $x=-81.68$. Thus we see that we can co-optimize for additional linear stability properties while balancing the need for a large SSP coefficient.

Another frequently desirable feature of time-stepping methods is that they include the imaginary axis or points near the imaginary axis. This is particularly desirable when solving hyperbolic PDEs. In our case, the optimal SSP five stage ($s=5$) method with $p=3$ and $p_{lin}=4$ does not include the imaginary axis at all ($y=0$). Even when we consider the close neighborhood of the imaginary axis (points which are closer than $10^{-5}$ to the imaginary axis) it only allows these values up to a value of $y=.9099$. Figure (4) shows the linear stability of this region in blue. If we wish to increase these values, we can obtain a method whose linear stability region is shown in red. This method clearly captures much less of the real axis, but it includes the imaginary axis up to the value of $y=5.28$. The zoomed image of the region of linear stability in Figure (4) shows this difference clearly. The SSP coefficient of the second method is significantly less, at ${\cal{C}}=4.1322$, but this may be a worthwhile trade-off where needed.

The methods in this subsection are meant to represent a small slice of the range of possible properties for which one may co-optimize. As in [27], it is possible to optimize SSP methods for particular linear stability regions, and it is also possible to consider other desirable properties with respect to which we can to co-optimize. Furthermore, other optimization routes are possible, such as starting with a known linear stability polynomial or a particular form of a method that is known to have desirable properties.

2.4 Numerical Results

The methods above were found by the optimization code [5], where the order conditions are imposed as equality constraints. It is a good idea to check that the methods indeed perform as designed on a series of linear and nonlinear test cases. These convergence studies are reported in Section 2.4.1. Next, we wish to test how the methods perform in terms of the predicted SSP time-step. Clearly, the SSP coefficient gives a guarantee that the desired property is preserved by the time-stepping method at the corresponding time-step. In Section 2.4.2 we study the difference between the guaranteed time-step and the actual time-step at which the desired nonlinear stability property (in this case the total variation diminishing property) still holds. A close agreement between these two time-steps demonstrates the relevance of the SSP property for typical cases.

2.4.1 Verification of the linear and nonlinear orders of convergence

Example 1.1: Study of the convergence rate for a nonlinear ODE To verify the nonlinear order of the implicit Runge–Kutta methods with nonlinear order $p$ and linear order $p_{lin}$ we us a nonlinear system of ODEs

[TABLE]

known as the van der Pol problem. We use $\epsilon=10$ and initial conditions $u_{0}=(0.5;0)$.

This problem was tested with methods with number of stages $s$, linear order $p_{lin}$ and nonlinear order $p$. Each methods is designated by $(s,p_{lin},p)$. Of each nonlinear order we test two methods: one with low stages and $p_{lin}=p$ and one with high number of stages and linear order. The methods we test are $(s,p_{lin},p)=(2,2,2),(7,7,2),(3,3,3),(7,7,3),(4,4,4),(7,7,4)$. We use $\Delta t=\frac{1}{250},\frac{1}{350},\frac{1}{450},\frac{1}{550},\frac{1}{650}$ and step forward to time $T_{final}=1.0$. The errors are calculated by using a very accurate solution calculated by MATLAB’s ODE45 routine with tolerances set to AbsTol=RelTol= $10^{-14}$. In Figure 6 we show that the $log_{10}$ of the errors in the first component vs. the $log_{10}$ of the time-step $\Delta t$. The orders (slope of the line) are taken by taking a linear fit using MATLAB’s polyfit. We observe the rate of convergence expected for the nonlinear order of the method.

Example 1.2: Study of the convergence rate for a linear PDE To verify the linear order of convergence of these methods we solve a linear advection problem $U_{t}+U_{x}=0$ with periodic boundary conditions and initial conditions $U_{0}(x)=\sin(x)$ on the spatial domain $x\in[0,2\pi]$. We discretize the spatial grid with $N=11$ equidistant points and use the Fourier pseudospectral differentiation matrix $\mathbf{D}$ [14] to compute $F\approx-U_{x}$. In this case, the solution is a sine wave, so that the pseudospectral method is exact, and the spatial discretization contributes no errors. For this reason, we use grid refinement in time only, and use a range of time steps, $\Delta t=\lambda\Delta x$ where $\lambda=\frac{1}{10},\frac{2}{10},\frac{3}{10},\frac{4}{10},\frac{5}{10},\frac{6}{10},\frac{7}{10},\frac{8}{10},\frac{9}{10}$ to compute the solution to final time $T_{f}=5.0$. The errors are measured in the $\ell_{2}$ norm.

Two methods from each order were tested for convergence on this problem. In Figure 7 we show the convergence plots for the $s=p_{lin}$ methods the $s=p_{lin}-1$ methods, both with nonlinear order $p=4$. The slopes were measured in the region before round-off error dominates. We observe that the design-order $p_{lin}$ of each method for linear problems is apparent.

Example 1.3: Study of the convergence rate for a nonlinear PDE To verify the nonlinear order on PDE, we solve the Buckley-Leverett equation which is commonly used to model two-phase flow through porous media:

[TABLE]

on $x\in[0,2\pi]$, with periodic boundary conditions. We take $a=\frac{1}{3}$ and initial condition $u(x,0)=\sin(x)$. For the spatial discretization, we use a Fourier pseudospectral method as above. $\Delta t=\lambda\Delta x$ with $\lambda=\frac{1}{32},\frac{1}{16},\frac{1}{8},\frac{1}{4},\frac{1}{2}$, and $\Delta x=\frac{\pi}{4}$, evolved to final time $T=2.0$. The errors are calculated by using a very accurate solution calculated by MATLAB’s ODE45 routine with tolerances set to AbsTol=RelTol= $10^{-14}$. In Figure 6 we show that the $log_{10}$ of the errors in the first component vs. the $log_{10}$ of the time-step $\Delta t$. The orders (slope of the line) are taken by taking a linear fit using MATLAB’s polyfit. We observe the rate of convergence expected for the nonlinear order of the method.

2.4.2 Verification of the SSP property

Example 2: The methods above are selected to optimize the SSP coefficient ${\cal{C}}$. This value guarantees that the desired strong stability property that is satisfied by the forward Euler method will be preserved under the condition $\Delta t\leq{\cal{C}}\Delta t_{\textup{FE}}$. It is interesting to see whether for typical problems this value is predictive of the actual time-step at which the desired strong stability properties are violated. As an example, we consider the linear advection equation

[TABLE]

on the domain $x\in[-1,1]$ with periodic boundary conditions and the step function initial condition

[TABLE]

We approximate the spatial derivative with a first order finite difference method, which is total variation diminishing (TVD) for all values $\Delta t\leq\Delta x$. This simple example is chosen as our experience has shown [11] that this problem often demonstrates the sharpness of the SSP time-step. For all of our simulations, we use a fixed grid of size $\Delta x=\frac{1}{300}$, and a time-step $\Delta t=\lambda\Delta x$, where we choose values of $\lambda$ starting from the minimum of ${\cal{C}}/10$ or $0.1$. We then increase the value of $\lambda$ until the TVD property is violated. To measure the effectiveness of these methods, we consider the maximum observed rise in total variation defined by

[TABLE]

over the first $N=20$ steps. We are interested in the time-step in which this rise becomes evident, so we consider a violation if the TV rises $10^{-10}$ or more at any time-step. The value of $\lambda$ for which this violation in TV occurs is refined to the level of $10^{-12}$. In Figure 8 we plot the maximal rise in total variation over different values of $\lambda=\frac{\Delta t}{\Delta x}$ for a selection of methods. It is interesting to observe that the change in behavior of the methods is quite sharp – once a certain threshold is reached, the SSP coefficient goes from being very small ($<10^{-14}$) to being much larger ($>10^{-4}$) with small changes in the time-step.

We are also interested in difference between the predicted and observed value of ${\cal{C}}$ (or equivalently, $\Delta t$) at which this violation happens. Table 4 provides these values for the case $p=2$. We observe that for all the methods with $p=2$ and $p_{lin}<s$, the observed time-step for which the TVD property is preserved matches the predicted time-step up to $\approx 10^{-10}$. For these cases, the SSP coefficient is an excellent predictor of the actual value for which the TVD property breaks down.

For the cases where $p_{lin}=s\leq 7$, and $p=2$, the observed time step is also within $\approx 10^{-10}$ of the predicted value. This is true also for the methods with $p=3$. However, when $p_{lin}=s=8,9$, and $p=2,3,4$, the observed time step is much larger than the predicted time-step ($\approx 10^{-1}$ or $10^{-2}$). Also, for $p=4$, the $s=p_{lin}=5$, we also have a larger difference between the observed and predicted time-step. Finally, in the cases where $s=p_{lin}-1$, the observed time step is usually much larger than the predicted time-step ($\approx 10^{-1}$ or $10^{-2}$), and the SSP coefficient is not sharp at all.

3 Optimal SSP IMEX Runge–Kutta methods with $p_{lin}\geq p$

For cases in which the problem has a linear component which restricts the time-step significantly, and a nonlinear component which does not, it may be advantageous to treat the two components differently by using an implicit-explicit method. Consider an initial value problem of the form

[TABLE]

which is semi-discretized to give an ODE of the form

[TABLE]

To step the two components forward in time, one explicitly and one implicitly, we use an Implicit-Explicit (IMEX) Runge–Kutta method, which is a particular kind of additive Runge–Kutta method. An $s$-stage Additive Runge–Kutta (ARK) method is defined by two $s\times s$ real matrices $\mathbf{A},\tilde{\mathbf{A}}$ and two real vectors ${\bf b},{\tilde{{\bf b}}}$ such that:

[TABLE]

We select $\mathbf{A}$ to be a strictly lower triangular matrix, so that this additive Runge–Kutta method contains an explicit method (represented by ($\mathbf{A},{\bf b}$)) used for the non-stiff part $F$. As above, we limit ourselves to diagonally implicit methods and require $\tilde{\mathbf{A}}$ to be a lower triangular matrix, where ($\tilde{\mathbf{A}},{\tilde{{\bf b}}}$) represents an implicit Runge–Kutta (DIRK) method used for the stiff part $G$, as in [18, 2, 20].

IMEX methods were first introduced by Crouziex in 1980 [7] for evolving parabolic equations. In 1997, Ascher, Ruuth, and Wetton [1] introduced IMEX multi-step methods for time dependent PDEs, notably convection-diffusion equations, and in the same year Ascher, Ruuth and Spiteri [2] presented the IMEX Runge–Kutta schemes for such problems. Although the authors were not focused on designing SSP pairs, some of methods are in fact known SSP methods. For example, the method in [2, Section 2.4] is the midpoint explicit and implicit SSP methods. In this work, all the implicit methods are SDIRK methods, and most have nice properties such as L-stability, with diagonal values $\gamma$ that are those mentioned in the books of Hairer, Norsett, and Wanner [12, 13].

Implicit methods are often particularly desirable when applied to a linear component. In this case, the order conditions simplify. In [4], Calvo, Frutos, and Novo developed IMEX pairs for the case where the implicit component $G$ is linear. This was the first work where the linear and nonlinear orders were separated. The methods they produced had nonlinear implicit order $p_{im}=2$, nonlinear explicit order and linear order $p_{ex}=p_{lin}=3,4$.

Kennedy and Carpenter [20] derived IMEX Runge–Kutta methods based on singly diagonally implicit Runge–Kutta (SDIRK) methods. This work introduced sophisticated IMEX methods with good accuracy and stability properties, as well as high quality embedded methods for error control and other features that make these methods usable in complicated applications.

The first IMEX methods that had SSP properties were considered by Pareschi and Russo in [29]. In this work, explicit SSP Runge–Kutta methods L-stable implicit components. These schemes were designed to be asymptotically preserving. The authors listed a $p_{ex}=p_{im}=p=3$ order conditions, and also presented a table depicting the number of coupling conditions under each underlying assumptions. In this work, the authors designed their methods while enforcing the condition that the abscissas of the explicit and implicit methods were equal ${\bf c}=\tilde{{\bf c}}$. They observed that under the assumption ${\bf c}={\tilde{{\bf c}}}$ and ${\bf b}={\tilde{{\bf b}}}$ the coupling conditions were redundant for the $p=3$. It is worth noting the SDIRK implicit pairs from [29] have the same $\gamma$ value listed in the books of Hairer, Norsett, and Wanner [12, 13]. Further work by Boscarino, Pareschi, and Russo [3] observed that the order reduction phenomenon disappears when ${\bf b}={\tilde{{\bf b}}}$. In this work they presented a globally stiffly accurate, third order method with non-negative explicit coefficients (named BPR(3,5,3).

Higueras [17] was first to consider SSP IMEX methods where both explicit and implicit pairs of the methods were SSP. In this work, she presented conditions for an IMEX Runge–Kutta methods to be SSP, and listed some previously known methods, such as methods from [2], that satisfy the SSP conditions. The methods that appear in this work are of order $p\leq 3$, and there is no distinction made between the explicit, implicit and linear order of the methods. In [18] Higueras presented order barriers and other characteristics of strong stability preserving additive Runge–Kutta methods. This work contains several very important properties of SSP IMEX pairs and is necessary to the understanding of the structure of these methods. In later work, [28] several known methods were analyzed in terms of their linear stability region and SSP performance on astrophysics-type problems. The implicit second order methods considered are L-stable (and thus A-stable as well), but the third order method was not A-stable. In their work in [19], the authors mainly focused on second order, three-stage explicit methods that implicit SSP pairs. In this work, the main focus was on methods with suboptimal SSP coefficients that offered significant additional benefits (including large linear stability regions).

Our approach to deriving optimal SSP methods is slightly different from [17, 19]. First, we use a slightly different formulation of the SSP condition, which facilitates the construction of an optimization problem. Next, we introduce a coefficient which relates the strong stability condition of $F$ to that of $G$, and make an observation about time step of optimal methods using this coefficient, which also simplifies the optimization routine. But the main distinction between our work and prior work is the fact that our investigations focus on higher order methods, and especially methods with higher linear order.

To analyze an IMEX method in the SSP context, we assume as in [17, 19] that when each component ($F$ or $G$) of the ODE is stepped forward individually with a forward Euler method, the new solution will satisfy a strong stability property in some desired convex functional $\|\cdot\|$, but under very different time-step restrictions:

[TABLE]

and

[TABLE]

If $K$ is small, the $G$ component will make the overall allowable time-step of the method very small, so it may be worthwhile to treat this component implicitly. We know [24, 11] that the allowable time-step for an $s$-stage SSP implicit Runge–Kutta method of order $p>1$ can only be expected to be, at most, $2s$ times the allowable explicit forward Euler time-step $\Delta t_{\textup{FE}}$, which does not typically offset the additional cost of solving the implicit system. However, in the previous discussions we show that the allowable time-step of an implicit method is in fact more than twice, and in the new methods presented in Section 2 up to six times, as large as that of an explicit method with the same order and number of stages and order. Furthermore, in the case where $G$ is a linear operator, the cost of the implicit solver is much smaller and under such circumstances it may be worthwhile to use an implicit SSP Runge–Kutta method for this part of the problem. However, if $F$ is nonlinear, or if for any other reason the cost of implicitly solving the system for $F$ is large, we wish to use an explicit time-stepping method for this part. The IMEX approach may also be desirable if the value of $K$ is not small (perhaps even infinitely large) but other factors, such as linear stability requirements, greatly limit the step size if $G$ is treated explicitly.

The main contribution of our work is to present SSP IMEX methods with $p_{lin}\geq 4$. These methods have not appeared in the literature, and while it is not clear yet how useful they will be in actual applications, we believe that understanding the SSP bounds on such methods and its dependence on the linear and nonlinear orders is of value.

3.1 Formulating the Optimization Problem

To formulate the optimization method for this problem, we stack the matrix $\mathbf{A}$ and the vector ${\bf b}$, padded with zeroes, into a square matrix $\mathbf{S}$ (as we did in Equation (9)). Similarly, we convert $\tilde{\mathbf{A}}$ and ${\tilde{{\bf b}}}$ into a matrix $\tilde{\mathbf{S}}$. We now rewrite the method (16) as

[TABLE]

Now, we add the terms

[TABLE]

where

[TABLE]

From this formulation we see that if $R{\bf e}$, $P$, and $Q$ are all positive component-wise, then the resulting method is simply a convex combinations of forward Euler steps $Y+\frac{\Delta t}{r}F(Y)$ and $Y+\frac{\Delta t}{\tilde{r}}G(Y)$, and therefore will preserve the strong stability property in the desired convex functional $\|\cdot\|$, under the time-step restriction $\Delta t\leq\min\left(r\Delta t_{\textup{FE}},\tilde{r}K\Delta t_{\textup{FE}}\right)$. We observe that the optimal methods will have $r=\tilde{r}K$, so the optimization problem becomes:

Maximize $r$ such that

[TABLE]

where the first three inequalities are understood component-wise, and $\tau_{k}(\mathbf{A},{\bf b},\tilde{\mathbf{A}},{\tilde{{\bf b}}})$ are the order conditions, described in the sections below. This optimization gives the Butcher coefficients $\mathbf{A}$, ${\bf b}$, $\tilde{\mathbf{A}}$, and ${\tilde{{\bf b}}}$, and an optimal value of the SSP coefficient ${\cal{C}}=r$. Of course, for each value of $K$, a different method will be optimal.

Note that although this process defines a sufficient condition for the resulting method to be SSP with SSP coefficient ${\cal{C}}=r$, there is no reason to expect this condition to be necessary. In particular, this formulation ignores the possible interactions between $F$ and $G$ that would result in more relaxed SSP conditions.

Once again, the order conditions (20d) are a key piece of the optimization, and serve as equality constraints. These conditions will be described in the next two subsections.

3.1.1 Linear Order Conditions

First order $p_{lin}=1$: The order conditions for an additive method (16) to be first order are

[TABLE]

This has two conditions: one for the explicit part and one for the implicit part. To generalize these order conditions, we define $\Phi_{1}=\{\phi_{1,1},\;\phi_{1,2}\}=\{{\bf b}^{T},\;{\tilde{{\bf b}}}^{T}\},$ and the order conditions become ${\bf b}^{T}{\bf e}=1,\;\;\;{\tilde{{\bf b}}}^{T}{\bf e}=1.$

Second order $p_{lin}=2$: Define $\Phi_{2}=\Phi_{1}\otimes\{\mathbf{A},\;\tilde{\mathbf{A}}\}=\{\phi_{1,1}\mathbf{A},\;\phi_{1,1}\tilde{\mathbf{A}},\;\phi_{1,2}\mathbf{A},\;\phi_{1,2}\tilde{\mathbf{A}}\},$ and the order conditions are each of these right-multiplied by ${\bf e}$ and set equal to $\frac{1}{2}$. Recall that ${\bf c}=\mathbf{A}{\bf e}$ and ${\tilde{{\bf c}}}=\tilde{\mathbf{A}}{\bf e}$ The order conditions then become

[TABLE]

Third order $p_{lin}=3$: use $\Phi_{3}=\Phi_{2}\otimes\{\mathbf{A},\;\tilde{\mathbf{A}}\}$ and right-multiply by ${\bf e}$ to obtain:

[TABLE]

[TABLE]

To obtain higher linear order, we define

[TABLE]

and the resulting new order conditions for this order take the form

[TABLE]

To go from linear order $p_{lin}=q-1$ to linear order $p_{lin}=q$ we require an additional $2^{q}$ order conditions. We observe that unlike in the implicit case, in the IMEX case the number of order conditions grows very rapidly even when both $F$ and $G$ are linear. This is due to the many coupling order conditions that are required for the two components to interact properly.

3.1.2 Nonlinear Order Conditions

The order conditions for first and second order are the same for both linear and nonlinear problems. In this section, we look at the order conditions for nonlinear orders $p^{e}=3,4$ for the explicit part and $p^{i}=3,4$ for the implicit part.

Nonlinear order $p=3$: For third order, we have the linear order conditions for the explicit and implicit parts, and the nonlinear coupling conditions given in the section above. In addition, we have the following nonlinear condition on the explicit part:

[TABLE]

When we wish to have $p^{e}=3$ but use a linear $G$ (i.e $p^{i}=2$) we need to include all the linear order conditions for $p_{lin}=3$, as well as the nonlinear equations for the explicit part (22a), and we must include the coupled order conditions in (22c). We can neglect the nonlinear implicit conditions (22b) and (22d).

Nonlinear order $p=4$: We include all the linear and nonlinear order conditions above, for the explicit and implicit parts as well as the coupled conditions. In addition, we have the nonlinear order condition Equation (23a) for the explicit method:

[TABLE]

And the nonlinear coupled order conditions

[TABLE]

[TABLE]

In the case where $G$ is linear and we wish to have $p^{i}=2$, we must satisfy all the linear order conditions, the nonlinear explicit conditions (23a), and the coupling conditions (24), but we can neglect the nonlinear implicit conditions (23b) and the nonlinear coupling conditions (25)

3.2 Optimal SSP IMEX Methods

We formulate the optimization problem above in a MATLAB routine following [25] and use it to generate optimized methods for different choices of $K$. We focus primarily on cases in which the number of stages $s$ is the same as the linear order $p_{lin}$ or a little larger. The nonlinear orders of interest are $p=2,3,4$. We note that in finding optimal methods, we frequently observed convergence to higher order than we required or expected. For example, in the case where we only require $p^{i}=2$ and $p^{e}>2$ the optimal SSP methods generally have ${\bf b}={\tilde{{\bf b}}}$ and ${\bf c}={\tilde{{\bf c}}}$, as we would expect from [18] and due to this the nonlinear implicit order is higher than required, $p^{i}=3$. However, we also find that in some cases we required only $p^{i}=2$ and $p^{e}=4$, but the methods converged to $p^{i}=p^{e}=4$. For this reason, in the following sections we describe mostly optimal methods that have $p=p^{i}=p^{e}$ (though some exceptions exist as is noted below). In some cases, especially where the number of stages $s$ is large and where there are more order constraints, the optimization routine had difficulties converging and we are not comfortable stating that the methods found are optimal. For this reason, we refer to the methods as optimized rather than optimal.

IMEX pairs for $K=\infty$ First, we look at the case where there is no constraint resulting from the $G$ component, which is equivalent to setting $K=\infty$. We reformulate the optimization problem to account for this by shutting off the $\tilde{\mathbf{S}}$ terms in (20a) – (20c). The case of $K=\infty$ is equivalent to the case where the implicit method is non-SSP while the explicit method is SSP. In many such cases, researchers have focused on finding methods where the explicit component is SSP and the implicit component is A-stable, L-stable, or stiffly stable. These methods were mostly of order $p\leq 3$ with few $p=4$ methods. It is known from the order barrier on explicit SSP Runge–Kutta methods that methods of order $p>4$ cannot exist. In our case, we focus primarily on methods that have $p_{lin}\geq 4$.

We first study methods with $p_{lin}\geq p$, and we do not optimize for linear stability: our rationale for studying these methods is to demonstrate that where there is no SSP constraint on $G$ we are able to find implicit pairs for optimal explicit SSP methods for $F$. Using our optimization routine we found methods with $s$ stages, linear order $p_{lin}$ and nonlinear orders $p=p^{i}=p^{e}=2,3$ with $K=\infty$ with $s=p_{lin},p_{lin}+1,p_{lin}+2$. Initially, we imposed a non-negativity conditions on all the coeffiicents, and were able to find many good methods. Later, we relaxed this condition and allowed negative coefficients in $\tilde{\mathbf{S}}$, which allowed us to find additional methods. Table 5 contains the SSP coefficients of these methods, which clearly match those of the explicit LNL SSP Runge–Kutta methods found in [10] and repeated in Table 2. The SSP coefficients of methods that have negative values in $\tilde{\mathbf{S}}$ are listed in bold. These methods do not have large linear stability regions, as they were not optimized for this purpose. However, these results show that the addition of an implicit component that does not have its own SSP constraint does not have an adverse effect on the size of the SSP IMEX method despite the addition of many new order conditions (coupling conditions) which must be satisfied.

Next, we look for methods that pair with popular SSP explicit Runge–Kutta methods and have large linear stability regions. We found a family of methods that pair with the three-stage third order Shu-Osher method:

[TABLE]

These implicit pair methods are defined by ${\bf b}={\tilde{{\bf b}}}$, and

[TABLE]

with

[TABLE]

The methods are third order and for $\beta>\frac{1}{2}$ are all A-stable. In particular, one member of this family is the nice looking method

[TABLE]

Another appealing method in this family occurs for the value $\beta=\frac{\sqrt{3}}{6}+\frac{1}{2}$. In this case the two non-zero values on the diagonal are equal, and so we have a type of SDIRK method (though with the first diagonal equal to zero). Note that this value of the diagonal may look familiar: in fact, it is the same as the diagonal value of the $s=2$, $p=3$ SDIRK method value in [12, Table 7.2], but that method cannot be paired with the Shu-Osher method.

This family of methods can be compared with two of Pareschi and Russo’s methods in [29] that also pair with the explicit SSPRK(3,3), Shu Osher method. The first is a three stage method that is singly implicit and L-stable and has both nonzero elements on the diagonal of $\tilde{A}$ but is only of overall order $p=2$. The second method is a four stage method that is singly implicit, has four nonzero elements on the diagonal of $\tilde{A}$ but is of overall order $p=3$. Depending on the cost of inverting the operator and the need for L-stability, the family of methods we present, which require only two inverse solves (and possibly with the same diagonal value) and are A-stable but not L-stable, may be preferable.

A useful and efficient method is Ketcheson’s explicit SSP Runge–Kutta method with ten stages and nonlinear order $p^{e}=4$. We found an SDIRK method which pairs well with Ketcheson’s SSPRK(10,4) method and has a very large linear stability region. The pair has linear order $p_{lin}=4$, and although it was produced by the optimizer under the assumption that $G$ is linear, we obtained better-than-expected nonlinear order $p^{i}=3$ due to the equality ${\bf c}={\tilde{{\bf c}}}$. The coefficients of this method are given in the appendix and the linear stability region is plotted in Figure 10. Note that the stability region this method crosses the real axis at $r_{1}\approx-102,775$ and the imaginary axis at $r_{2}\approx\pm 138,891$. The optimization code was quickly able to match any value of $r=\min\{|r_{1}|,|r_{2}|\}$ that we requested. This leads us to suspect that there may well be an A-stable method that pairs with Ketcheson’s explicit SSPRK(10,4) method; however, the number of coefficients here is prohibitive and unlike the case of the Shu-Osher SSPRK(3,3) method solving for this analytically seems difficult.

It is important to note that the explicit methods in the pair above all have an optimal SSP coefficient among methods of its number of stages and order, and the coupling with an implicit method under the assumption $K=\infty$ does not adversely affect the SSP coefficient, so these methods are the best possible in their class in terms of SSP coefficient.

Finally, we found an IMEX pair with $s=10$ stages, linear order $p_{lin}=6$, and nonlinear order $p=3$. The explicit part has SSP coefficient ${\cal{C}}=3.3733$. This method has an SDIRK implicit part, with a good linear stability region as shown in in Figure 10. The stability region crosses the real axis at $r_{1}\approx-1,350$ and the imaginary axis at $r_{2}\approx\pm 600$.

IMEX pairs for small $K$ An interesting case that has not been considered extensively in the literature (except indirectly in [17]), is the situation where the implicit component $G$ introduces a very tight time-step restriction for the strong stability property to be satisfied. This is the case where the value of the parameter $K$ is very small. We expect that in these cases the SSP coefficient of the IMEX pair will be correspondingly limited, and indeed the coefficients are small. Table 6 gives the SSP coefficients of some methods we found using our optimization routine. We see in these tables that for IMEX SSP Runge–Kutta methods the SSP coefficient is reduced as both the linear and nonlinear orders rise; however, in general the SSP coefficients for $p=2$ and for $p=3$ are very similar, and for these cases the increasing linear order causes the decrease in the SSP coefficient. In addition, we see that there may be some modest benefit to tailoring the method to the actual value of $K$ in the problem; whether this is true in practice is investigated in Section 3.3.2. The coefficients $\mathbf{A},\tilde{\mathbf{A}},{\bf b}$ and ${\tilde{{\bf b}}}$ of these methods can be downloaded at [6].

These methods are not at all optimized in terms of the linear stability region, as the main constraint of interest for these methods comes from the SSP condition, and indeed their linear stability region is similar to explicit methods. We believe these methods may be useful where the cost of the implicit solution is negligible while the SSP condition is needed for small values of $K$. However, even in cases where these methods are not directly useful, they give us an upper bound on the possible SSP coefficients of SSP IMEX methods for typical values of $K$, and so serve as a guide to what we can expect when we co-optimize these methods with other desirable properties. The numerical tests, particularly those in Section 3.3.2, show the performance of some of these methods for the types of problem that require the SSP property for both the explicit $F$ and implicit $G$.

3.3 Numerical Experiments

3.3.1 Convergence studies

To test the accuracy of the SSP IMEX Runge–Kutta methods, we consider a linear or nonlinear explicit part. As our main motivation for these methods are where the implicit component is linear, we test the methods on problems where $G$ is linear. In these tests we confirm that our methods give us the expected linear and nonlinear orders for the explicit and implicit parts.

Example 3.1: convergence study with explicit linear advection and implicit linear diffusion. We approximate the solution to the equation

[TABLE]

where $\epsilon=0.01$ with periodic boundary conditions and sine wave initial conditions. We discretize the spatial grid with $N=8$ equidistant points and use the Fourier pseudospectral differentiation matrix for both the first and second derivative [14]. The advection part is dealt with explicitly and the diffusion implicitly (and exactly using MATLAB’s backslash operator). Temporal grid refinement with pseudospectral approximation of the spatial derivative. We use a range of time steps, $\Delta t=\lambda\Delta x$, where $\Delta x=\frac{\pi}{4}$ and we pick $\lambda=\frac{1}{16},\frac{1}{8},\frac{1}{4},\frac{1}{2}$ to compute the solution to final time $T_{f}=5.0$. The methods we test are our IMEX SSP methods with $s$ stages, nonlinear order $p=p^{i}=p^{e}$ and linear order $p_{lin}$ $(s,p,p_{lin})=(3,3,3),(4,3,4),(5,4,4),(6,2,5),(6,3,6),(7,2,7)$ designed for $K=\infty$ and for $K=10^{-2}$. The $\ell_{2}$ errors are measured compared to the exact solution and shown in Figure 11. We note that the methods designed for $K=\infty$ and those designed for $K=10^{-2}$ perform essentially the same on convergence studies. The orders (slope $m$ of the line) are taken by taking a linear fit using MATLAB’s polyfit and given in the caption. We see that, as expected, the linear design order is apparent in this linear problem. It is interesting to note that the $(s,p,p_{lin})=(3,4,4)$ and the $(s,p,p_{lin})=(4,4,5)$ have the same slope, because their linear order is the same, but the latter method has a smaller error constant so the errors are smaller.

Example 3.2: convergence study with explicit Burgers’ and implicit linear advection. We approximate the solution to the equation

[TABLE]

with periodic boundary conditions and sine wave initial conditions. We discretize the spatial grid with $N=24$ equidistant points and use the Fourier pseudospectral differentiation matrix [14]. The Burgers’ flux $\left(\frac{1}{2}U^{2}\right)_{x}$ is dealt with explicitly and the linear advection implicitly (and exactly using MATLAB’s backslash operator). We perform grid refinement with pseudospectral approximation of the spatial derivative, we use a range of time steps, $\Delta t=\lambda\Delta x$ where we pick $\lambda=\frac{1}{128},\frac{1}{64},\frac{1}{32},\frac{1}{16},\frac{1}{8}$ to compute the solution to final time $T_{f}=0.8$ before the shock forms. The methods we test are our IMEX SSP methods with $(s,p,p_{lin})=(3,3,3),(4,3,4),(5,4,4),(6,2,5),(6,3,6),(7,2,7)$ designed for $K=\infty$ and for $K=10^{-2}$. The $\ell_{2}$ errors are measured compared to the the results from MATLAB’s ode45 and graphed in Figure 12. The orders (slope $m$ of the line) are taken by taking a linear fit using MATLAB’s polyfit and given in the caption. We see that the nonlinear order $p$ typically dominates: this occurs since the linear order $p_{lin}$ is no smaller than the nonlinear orders. In the case $(7,2,7)$ for $K=10^{-2}$ we get a better order than expected by one, due to a small error constant ($\approx 10^{-8}$) for this method.

3.3.2 Numerical verification of the SSP properties of these methods

In this section we focus on the behavior of our time-stepping methods on problems which require the SSP property. In particular, we use out IMEX SSP methods for problems where the desired property to be preserved is the total variation diminishing (TVD) property and where the spatial discretization is simple enough that the theoretical bound on the TVD condition is easy to see.

**Example 4.1: Explicit linear advection with implicit linear advection. ** Consider the linear equation

[TABLE]

with periodic boundary conditions and step function initial conditions

[TABLE]

We use a first order upwind differencing the linear advection terms, with $N=301$ points in space. The time-step is set to $\Delta t=\lambda\Delta x$, where different values of $\lambda$ are tested to observe the value at which the discretization is no longer TVD (i.e. the total variation rises by more than $10^{-10}$ at any given time-step). We first compared the Predicted $\Delta t$ and Observed $\Delta t$ for which the TVD property is violated. The results are very similar to those in Table 8 and so we do not report them here.

We consider several IMEX methods that were optimized for $K=1,\frac{1}{10},\frac{1}{100}$ and see how their observed time-step for the TVD to be satisfied compares to the time-step predicted by the theory, even when the methods were designed for different values of $K$. We also compare the performance of these methods to the situation where both components are advected by the optimal explicit LNL SSP Runge–Kutta method of the corresponding number of stages, nonlinear order, and linear order. For completeness, we also include a comparison to treating both components implicitly using an optimal implicit LNL SSP Runge–Kutta method.

Table 7 contains the observed and predicted time-step for which the TVD property is violated. Clearly, the observed $\Delta t$ is never smaller than the predicted $\Delta t$, as the SSP property provides a guarantee of this property. We note that while the predicted and observed TVD time-step are always exact when both components are treated implicitly, and frequently so when they are both treated explicitly, this is not as often the case where we use an IMEX method (though in many cases, particularly those noted in Table 8, we do have sharp agreement). It is notable that using an IMEX method that was optimized for a value of $K$ close to the actual value of $\frac{1}{\omega}$ we generally see larger allowable time-step, though this difference may not be large enough to generate optimal methods for many different values. The main question we wish to answer is whether these methods are useful as replacements to the fully explicit method. In these cases we see that using an IMEX scheme allows us to take a time-step that is between $1.5$ and $2$ times the fully explicit time-step (see the ratio column). There are very limited circumstances in which this larger time-step is truly enough to offset the cost of the implicit solver; however, this is what one would have expected from the fact that the SSP coefficient of the explicit methods are bounded by the number of stages while the SSP coefficient of the implicit methods are bounded by twice the number of stages. The more interesting observation in this table may be the ratio of the fully implicit allowable time-step to the fully explicit allowable time-step: these can be quite large, possibly large enough to offset the implicit solver cost in some cases.

**Example 4.2: Explicit Burgers’ with implicit linear advection. ** We consider the equation

[TABLE]

with periodic boundary conditions and step function initial conditions

[TABLE]

We use a first order upwind differencing for both the Burgers’ term and the linear advection term, with $N=301$ points in space. The methods used are those optimized for the value of $K$ that corresponds to the wavespeed $\omega=\frac{1}{K}$, and these IMEX methods have $s$ stages, nonlinear order $p^{e}=p^{i}=p$ and nonlinear order $p_{lin}$. The time-step is set to $\Delta t=\lambda\Delta x$, where different values of $\lambda$ are tested to observe the value at which the discretization is no longer TVD. This value of $\Delta t$ is then called the “Observed $\Delta t$” while the value we would expect from the SSP coefficient is called the “Predicted $\Delta t$”.

Table 8 contains the observed and predicted time-step for which the TVD property is violated. Clearly, the observed $\Delta t$ is never smaller than the predicted $\Delta t$, as the SSP property provides a guarantee of this property. As seen in the table, in some cases, the methods feature close agreement between the predicted and observed SSP coefficient, while in others the actual time-step allowed for this problem is greater than predicted. We note that these methods performed very similarly on the linear problem in Example 4.2 above, so the results were not repeated in a separate table. These results, and in particular the sharpness of the predicted allowable time-step in several cases, demonstrate that the theoretical SSP property is a good predictor of the actual behavior of the method.

4 Conclusions

In this work, we investigated implicit and IMEX SSP methods with very high linear order. We first considered implicit SSP Runge–Kutta methods that have high order for linear problems and order $p=2$ for nonlinear problems. We show that as is the case with explicit methods we are able to find implicit methods with $p_{lin}>6$ with no linear order barrier. We were further able to show that (as in [10]) the optimal methods had no better SSP coefficient that optimal methods that are diagonally implicit with a low storage formulation. Thus we can say that optimal implicit SSP Runge–Kutta methods for any linear order $p_{lin}$ and nonlinear order $p=2$ are diagonally implicit and low storage.

We continued our study of implicit linear/nonlinear (LNL) methods for larger nonlinear orders $p>2$ following the approach used in [10]. We designed SSP implicit Runge–Kutta methods where the linear order exceeded the nonlinear order (called LNL methods), and observed that the SSP coefficient for methods with linear order $p_{lin}\leq 9$ was similar whether we chose nonlinear order $p=2$ and $p=3$, and that if we went up to $p=4$, the SSP coefficient was not typically significantly reduced as long as we had sufficiently many stages. This implies that just as the case for the explicit LNL methods, as we increase the number of stages and maintain a moderate nonlinear order $p=3,4$, the linear order is the main constraint on the SSP coefficient. However, if we increase the nonlinear order to $p=5,6$, this did have a strong adverse effect on the size of the SSP coefficient. We also observed that the SSP coefficients of the implicit LNL SSP Runge–Kutta methods were up to six times as large as those of the corresponding explicit LNL SSP Runge–Kutta of the same number of stages and linear and nonlinear orders. We report the coefficients of methods with high linear order and moderate number of stages that have reasonable SSP coefficients, for example the LNL implicit RK $(s,p,p_{lin})=(6,4,6)$ method with SSP coefficient ${\cal{C}}=0.856$. These methods may be useful when a high linear order method is desired, while still reasonable for use with nonlinear problems. We then verified the convergence of these methods as well as the sharpness of the SSP coefficient on sample problems.

The second part of this paper focused on implicit-explicit methods, where as above the methods have higher linear order than nonlinear order. These methods are of value when we desire both the explicit and implicit parts to have the SSP property, and are particularly beneficial when the linear part of a problem is treated implicitly, and the nonlinear part explicitly. We found diagonally implicit non-SSP Runge–Kutta methods with large regions of linear stability that pair with well-known explicit SSP methods. Next, we found pairs of IMEX SSP methods, with high linear order, that were optimized for the relationship between the forward Euler conditions of each component. We verified the order of these methods on a variety of problems, and the sharpness of the SSP coefficient on a typical test cases, and conclude that these IMEX LNL methods demonstrate high order convergence and perform as desired in terms of SSP.

This work shows that it is possible to produce implicit and IMEX SSP methods of very high order if we consider only the linear order. For implicit methods, we found methods of up to order $p_{lin}=9$ and for the IMEX methods up to order $p_{lin}=7$, which show that the order barriers of the implicit and IMEX SSP Runge–Kutta methods apply only to the nonlinear order. The approach described in this work can be used to produce optimal SSP implicit and IMEX methods with additional desired properties, as well. The SSP coefficients presented in this work give a baseline with which to compare any methods that are SSP and also have other properties, and give a value to strive toward. As expected at the outset, the allowable time-step for these methods is very restricted, especially for the IMEX schemes and therefore the cost of the implicit solvers is still the major issue that arises in the solution of these methods. However, in particular cases (such as where the implicit term is a constant coefficient linear term) where the cost of the implicit solver is controllable, these methods may be of more than theoretical value.

Acknowledgements This work was supported by AFOSR grant FA9550-15-1-0235, and was partially supported by DOE NNSA ASC Algorithms effort, the DOE Office of Science AMR program at Sandia National Laboratory under contract DE-AC04-94AL85000

Appendix A Coefficients of some optimized methods

A.1 Methods for $K=\infty$

A.1.1 An IMEX pair based on Ketcheson’s SSPRK10,4 method

Ketcheson’s SSPRK10,4 method is given by a vector ${\bf b}$ where $({\bf b})_{i}=\frac{1}{10}$ for $1\leq i\leq 10$ and the coefficients in $\mathbf{A}$ are

[TABLE]

The SDIRK method described above that pairs with Ketcheson’s 10,4 method has coefficients ${\tilde{{\bf b}}}={\bf b}$ and ${\tilde{{\bf c}}}={\bf c}$ and the non-zero coefficients in $\tilde{\mathbf{A}}$ are:

[TABLE]

A.2 Methods for $K=10$

A.2.1 An IMEX pair with $p^{e}=p^{i}=3$, $s=p_{lin}=5$

The five stage IMEX pair with

The nonzero coefficients in $\mathbf{A}$ are given by:

[TABLE]

[TABLE]

[TABLE]

A.2.2 An IMEX pair with $p^{e}=p^{i}=4$, $p_{lin}=6$

The nonzero coefficients in $\mathbf{A}$ are given by:

[TABLE]

The nonzero coefficients in $\tilde{\mathbf{A}}$ are given by:

[TABLE]

[TABLE]

A.3 Methods for $K=100$

A.3.1 An IMEX pair with $p^{e}=p^{i}=3$, $s=p_{lin}=5$

The nonzero coefficients in $\mathbf{A}$ are given by:

[TABLE]

[TABLE]

[TABLE]

A.3.2 An IMEX pair with $s=10$, $p^{e}=4$, $p^{i}=3$, $p_{lin}=6$

The six stage IMEX pair with fourth order for nonlinear problems and sixth order on linear problems is given by the nonzero coefficients in $\mathbf{A}$:

[TABLE]

The nonzero coefficients in $\tilde{\mathbf{A}}$ are given by:

[TABLE]

[TABLE]