Relaxation Runge-Kutta Methods: Conservation and stability for Inner-Product Norms

TL;DR

This paper introduces a modification to Runge-Kutta methods that ensures conservation or stability under any inner-product norm, maintaining original accuracy and stability, with demonstrated effectiveness in numerical examples and hyperbolic PDEs.

Contribution

It presents a simple modification to Runge-Kutta methods that guarantees conservation or stability for any inner-product norm while preserving their original properties.

Findings

Modified methods ensure stability and conservation under any inner-product norm.

The methods retain the accuracy and stability of the original Runge-Kutta methods.

Numerical examples demonstrate effectiveness, including in hyperbolic PDE applications.

Abstract

We further develop a simple modification of Runge--Kutta methods that guarantees conservation or stability with respect to any inner-product norm. The modified methods can be explicit and retain the accuracy and stability properties of the unmodified Runge--Kutta method. We study the properties of the modified methods and show their effectiveness through numerical examples, including application to entropy-stability for first-order hyperbolic PDEs.

Figures

Tables

Table 1. Table 1 : Values of γ ∗ subscript 𝛾 \gamma_{*} from Theorem 9 for some well-known SSP methods. For 0 ≤ γ ≤ γ ∗ 0 𝛾 subscript 𝛾 0\leq\gamma\leq\gamma_{*} , the relaxation method has the same SSP coefficient as the original method.

Method	${\gamma }_{\ast }$
SSPRK(2,2)	2
SSPRK(s,2)	s/(s-1)
SSPRK(3,3)	3/2
SSPRK(4,3)	1
SSPRK(5,3)	1
SSPRK(9,3)	1
SSPRK(5,4)	1.312
SSPRK(10,4)	25/24

Equations

$$u^{\prime }(t)$$

$$u(t_0)$$

$$\frac{d}{dt}\Vert u(t){\Vert }^2=2⟨u,f(t,u)⟩\le 0.$$

$$\Vert u^{n+1}\Vert \le \Vert u^n\Vert .$$

$$⟨u,f(t,u)⟩=0.$$

$$\Vert u^{n+1}\Vert =\Vert u^n\Vert .$$

$$y_i$$

$$u(t_n+\Delta t)\approx u^{n+1}$$

$$f_i=f(t_n+c_i\Delta t,y_i)$$

$$\Vert u^{n+1}{\Vert }^2-\Vert u^n{\Vert }^2$$

$$=2\Delta t\sum _{j=1}^s{b}_j⟨u^n,f_j⟩+{\Delta t}^2\sum _{i,j=1}^s{b}_i{b}_j⟨f_i,f_j⟩$$

$$=2\Delta t\sum _{j=1}^s{b}_j⟨y_j,f_j⟩-2\Delta t\sum _{j=1}^s{b}_j⟨u^n-y_j,f_j⟩+{\Delta t}^2\sum _{i,j=1}^s{b}_i{b}_j⟨f_i,f_j⟩,$$

$$\Vert u^{n+1}{\Vert }^2-\Vert u^n{\Vert }^2$$

$$-{\Delta t}^2\sum _{i,j=1}^s{m}_{ij}⟨f_i,f_j⟩,$$

$$M=BA+A^{T}B-b{b}^T.$$

$$u(t_n+{\gamma }_n\Delta t)\approx u_{\gamma }^{n+1}$$

$$\begin{array}{rl}\Vert u_{\gamma }^{n+1}{\Vert }^2-\Vert u^n{\Vert }^2=2{\gamma }_n\Delta t\sum _{j=1}^s{b}_j⟨y_j,f_j⟩& -2{\gamma }_n{\Delta t}^2\sum _{i,j=1}^s{b}_i{a}_{ij}⟨f_j,f_i⟩\\ & +{\gamma }_n^2{\Delta t}^2\sum _{i,j=1}^s{b}_i{b}_j⟨f_i,f_j⟩.\end{array}$$

$${\gamma }_n=\frac{2{\sum }_{i,j=1}^s{b}_i{a}_{ij}⟨f_i,f_j⟩}{{\sum }_{i,j=1}^s{b}_i{b}_j⟨f_i,f_j⟩},$$

$$\Vert u_{\gamma }^{n+1}{\Vert }^2-\Vert u^n{\Vert }^2$$

$${\gamma }_n$$

$$\Delta E(\gamma ):=\Vert u_{\gamma }^{n+1}{\Vert }^2-\Vert u^{n+1}{\Vert }^2=0,$$

$$\Delta E(1)=\Vert u(t_n+\Delta t){\Vert }^2-\Vert u(t_n){\Vert }^2=2\Delta t\sum _j{b}_j⟨y_j,f_j⟩+\mathcal{O}({\Delta t}^{p+1}).$$

$${\gamma }_n=1+\mathcal{O}({\Delta t}^{p-1}).$$

$$b^T{c}^{k-1}-1/k$$

$$b^{T}v$$

$${\gamma }_n{b}^T{c}^{k-1}-1/k$$

$${\gamma }_n{b}^{T}v$$

$${\gamma }_n{b}^T{c}^{k-1}-{\gamma }_n^k/k$$

$${\gamma }_n{b}^{T}v$$

$${\gamma }_n(b^T{c}^{k-1}-{\gamma }_n^{k-1}/k)=\frac{{\gamma }_n}{k}(1-(1+{\Delta t}^{(k-1)(p-1)}))=\mathcal{O}({\Delta t}^{(k-1)(p-1)}).$$

$${\gamma }_n=1+\mathcal{O}({\Delta t}^{p-1}).$$

$$\Vert u_{\gamma }^1-u({\gamma }_1\Delta t)\Vert =\mathcal{O}({\Delta t}^{p+1}).$$

$$\rho (t^{\prime })=\rho (t)+1$$

$$\gamma ([t_1,t_2])=\gamma (t_1)\gamma (t_2)(\rho (t_1)+\rho (t_2)+1).$$

Full text

1 Motivation and background

Consider the initial value problem

[TABLE]

where111We consider real spaces for simplicity, but the methods developed here are also applicable in complex spaces. $u:{\mathbb{R}}\to{\mathbb{R}}^{m}$ and $f:{\mathbb{R}}\times{\mathbb{R}}^{m}\to{\mathbb{R}}^{m}$. In this work we focus on problems that are dissipative with respect to some inner-product norm:

[TABLE]

Here and throughout, $\langle\cdot,\cdot\rangle$ denotes an inner product and $\|\cdot\|$ the corresponding norm; we will sometimes refer to $\|u\|^{2}$ as energy. For dissipative problems, it is desirable that the numerical solution mimic (2):

[TABLE]

Herein, a method is called monotonicity preserving if it guarantees (3) for all problems satisfying (2).

We say the problem (1) is conservative if

[TABLE]

For such problems it is desirable to discretely conserve energy:

[TABLE]

A method is also called conservative if it guarantees (5) for all problems satisfying (4).

In many applications, numerical conservation or monotonicity preservation are of great importance. Not only do these properties guarantee that the solution remains bounded, but violation of these properties can lead to solutions that are unphysical and qualitatively wrong. Nevertheless, most numerical methods do not enforce these properties exactly, but only up to truncation errors. This is particularly true for explicit Runge–Kutta methods. Even if one considers only linear autonomous systems of ODEs, it has been shown that many explicit Runge–Kutta methods — including the classical 4th-order method — are not even conditionally monotonicity preserving; see [26, 24, 21, 25]. For general nonlinear autonomous systems, no methods of order greater than one are known to be monotonicity preserving [19]. For conservative problems, no explicit Runge–Kutta method of any order can enforce discrete conservation, even for linear autonomous problems. Formulas for the production of a general convex entropy function (which includes the special case of an inner-product norm) have been derived in [18].

Monotonicity preservation for Runge–Kutta methods with respect to an inner-product norm was studied by Higueras [14], with results closely connected to an earlier study of contractivity preservation by Dahlquist & Jeltsch[6]. However, to obtain results for explicit RK methods in that framework, it is necessary to require strict inequality in (2), making the results inapplicable for conservative problems or for typical high-order semi-discretizations of hyperbolic PDEs.

Unconditionally conservative and monotonicity-preserving methods exist, for instance, among the classes of implicit Runge–Kutta and (for Hamiltonian systems) partitioned Runge–Kutta methods; see [12] and references therein.

While classical explicit Runge–Kutta or linear multistep methods cannot preserve general quadratic invariants, some explicit energy-conservative methods have been developed by going outside these traditional classes (as is done also in the present work). Indeed, it is possible to modify any Runge–Kutta method to preserve energy or other first integrals by a technique known as projection; see e.g. [10, 4]. In this approach, at the end of each time step the solution is projected onto a desired set in order to ensure some property like conservation or monotonicity. The approach used in the current work can be viewed as a projection method where the projection is performed along a direction corresponding to the next time step update, but with an additional modification that the size of the time step is also modified. This approach was originally proposed by Dekker & Verwer [7, pp. 265-266], who noted that the classical four-stage RK method could be modified slightly to conserve energy while maintaining its accuracy. The idea was extended in [8] to a restricted class of fourth order methods. This was developed further in [5] by giving a general proof that applying the technique (without the step size adjustment) to a RK method of order $p$ results in a method of order at least $p-1$. Subsequent development in [5, 4, 17] focuses on embedded projection methods; i.e. methods that project in a direction determined by an embedded Runge–Kutta pair. The approach described in the present work could be viewed as a variant in which the “embedded” method is simply the identity map, but with an additional twist that requires reinterpretation of the new step solution as an approximation at a slightly different time.

Like embedded projection methods, the methods proposed here are:

• •

explicit

• •

conditionally conservative and monotonicity-preserving for general nonlinear ODEs

• •

arbitrarily high order accurate

• •

linearly covariant

Furthermore, they do not require partitioning or temporal staggering, and they inherit other useful properties (such as strong stability preservation) of a selected Runge–Kutta method. Preservation of more general (non-inner-product) functionals is also possible with modification similar to that described herein; see [22].

The main contributions of the present work are: first, to further develop these methods that, while not entirely new, seem to have been overlooked; second, to put them on a rigorous footing in terms of accuracy and stability properties; and third, to explore their properties through analysis and numerical experiments.

1.1 Energy evolution by Runge–Kutta methods

A Runge–Kutta method applied to (1) takes the form

[TABLE]

We make the usual assumption that $c_{j}=\sum_{i}a_{ij}$. For convenience we introduce the shorthand

[TABLE]

for the $i$th stage derivative. The change in energy from one step to the next is

[TABLE]

which can be rewritten using (6a) as

[TABLE]

The first sum on the right side of (7) is zero for conservative systems, and it is negative for dissipative systems if $b_{j}\geq 0$ for all $j$. However, the remaining two terms may lead to violation of the conservation or monotonicity property. Those two terms can be written together as the bilinear form

[TABLE]

where, letting $B$ denote the diagonal matrix with entries $b_{j}$,

[TABLE]

This is the traditional analysis used in studying symplecticity, algebraic stability, and other properties of RK methods (see e.g. [3, 7, 12]). If the matrix $M$ is positive semidefinite and the weights are nonnegative, then the method is said to be algebraically stable. Clearly, such methods are unconditionally monotonicity-preserving. If $M=0$, the method is said to be symplectic; clearly such methods are unconditionally conservative. Certain well-known implicit methods have these properties. However, explicit methods cannot be algebraically stable or symplectic.

2 Relaxation Runge–Kutta methods

The relaxation version of the method (6) is obtained by replacing (6b) with the update formula

[TABLE]

The only difference between (9) and (6b) is the factor $\gamma_{n}$ that multiplies the step size. We can think of the original Runge–Kutta method (6) as determining only the direction in which the solution will be updated, while the choice of $\gamma_{n}$ determines how far to step in that direction. From this point of view $\gamma_{n}$ is similar to the relaxation parameter used in some iterative algebraic solvers, and for this reason we refer to these methods as relaxation Runge–Kutta (RRK) methods.

With this change, (7) becomes

[TABLE]

We can eliminate the last two terms by setting

[TABLE]

so that

[TABLE]

In case the the denominator of (11) vanishes, we have $u^{n+1}=u^{n}$, so we can achieve conservation or monotonicity by taking simply $\gamma_{n}=1$. We thus define in place of (11)

[TABLE]

Because we will interpret $u_{\gamma}^{n+1}$ as an approximation to the solution at time $u(t_{n}+\gamma_{n}{\Delta t})$, it is important that $\gamma_{n}>0$. It is straightforward to show the following:

Lemma 1**.**

Let $\sum b_{i}a_{ij}>0$, let $f$ be sufficiently smooth, and let $\gamma_{n}$ be defined by (12). Then $\gamma_{n}>0$ for sufficiently small ${\Delta t}>0$.

Note that the condition in Lemma 1 holds for all methods of order two or higher.

Remark 1**.**

A more detailed analysis indicates that small enough here means that ${\Delta t}$ should be no more than about $1/L$, where $L$ is the Lipschitz constant of $f$. This is also the order of the absolutely stable step size when using any explicit method, so one may expect that using an absolutely stable step size will also yield $\gamma_{n}>0$.

For conservative systems this gives exact energy conservation; for dissipative systems it preserves dissipativity as long as all weights are non-negative.

Theorem 2**.**

Let $(A,b)$ be the coefficients of a Runge–Kutta method of order at least two. The corresponding relaxation Runge–Kutta method defined by (6a) and (9) with $\gamma_{n}$ defined by (12) is conservative. If the weights are nonnegative, then the relaxation method is also monotonicity preserving as long as ${\Delta t}$ is chosen so that $\gamma_{n}\geq 0$.

Remark 2**.**

The relaxation RK method still conserves linear invariants, which is important for instance in the semi-discretization of hyperbolic conservation laws.

Remark 3**.**

A version of Theorem 2 applicable only to conservative systems appears in [5, Thm. 2.1], along with a formula for $\gamma_{n}$ that is more computationally efficient (but correct only for conservative systems).

Remark 4**.**

The denominator of the expression for $\gamma_{n}$ in (11) is simply the norm of the step update, and can be computed with a single inner product. It has been pointed out by Hendrik Ranocha that by solving (7) for the term that is the numerator of the expression for $\gamma_{n}$, it can be computed using only $s$ inner products [20]. Thus $\gamma_{n}$ can be computed with just $s+1$ inner products.

The update formula (9) is equivalent to replacing the coefficients $b_{j}$ in the Runge–Kutta method with $\gamma_{n}b_{j}$. It can also be viewed roughly as taking a Runge–Kutta step of size $\gamma_{n}{\Delta t}$ in place of ${\Delta t}$, but notice that the original step size ${\Delta t}$ is still used in the computation of the stages $y_{i}$. Both viewpoints (rescaling $b$ and rescaling ${\Delta t}$) will be useful in our analysis.

We will see below in Lemma 4 that, for reasonable values of ${\Delta t}$, $\gamma_{n}$ is close to unity. Although the proof of this fact is technical, the result itself is not surprising, at least for conservative systems. For such systems, the value of $\gamma_{n}$ given by (12) is a solution of

[TABLE]

which is a quadratic function of $\gamma$ with one root at zero. Given that $\Delta E(1)={\mathcal{O}}({\Delta t}^{p+1})$ and that according to (10) $\Delta E\propto{\Delta t}^{2}$ for conservative problems, it is natural to expect that $\Delta E(\gamma)$ has a zero within ${\mathcal{O}}({\Delta t}^{p-1})$ of unity. The method described here can be viewed as using a line search to determine the step size that solves $\Delta E(\gamma)=0$ and thus exactly conserves energy.

For dissipative systems also, we will see that the last two terms in (7) are not important (in the sense of accuracy) to the numerical approximation of the energy evolution; i.e. the term $2{\Delta t}\sum_{j}b_{j}\langle y_{j},f_{j}\rangle$ approximates the energy evolution over one step to order $p$:

[TABLE]

2.1 Accuracy

At first glance, the RRK method (given by (6a) with (9)) seems to be not even consistent, since $\sum_{j}\gamma_{n}b_{j}=\gamma_{n}\neq 1$ in general. For a classical RK method, the condition $\sum_{j}b_{j}=1$, along with higher order conditions, is necessary for local consistency of a given order. But for an RRK method, because the coefficients depend on ${\Delta t}$, we can still obtain high order accuracy if the order conditions are nearly satisfied, which is true if $\gamma_{n}$ is sufficiently close to unity.

Theorem 3**.**

Let $(a_{ij},b_{j})$ be the coefficients of a Runge–Kutta method of order $p$. Consider the RRK method defined by (6a), (9) and suppose that

[TABLE]

Then:

(IDT method) If the solution $u^{n+1}_{\gamma}$ is interpreted as an approximation to $u(t_{n}+{\Delta t})$, the method has order $p-1$. 2. 2.

(RRK method) If the solution $u^{n+1}_{\gamma}$ is interpreted as an approximation to $u(t_{n}+\gamma_{n}{\Delta t})$, the method has order $p$.

Proof.

First we note that by taking appropriate linear combinations of the usual order conditions, the conditions for a RK method to have order $p$ can be written as (see [1])

[TABLE]

where $V_{p}$ is a set of vectors depending only on $a_{ij}$ whose specific elements are not important here. For the RRK method, we replace $b$ with $\gamma_{n}b$, obtaining the conditions

[TABLE]

Given a method that satisfies (14), clearly (15b) is satisfied as well, while the left hand side of (15a) is ${\mathcal{O}}({\Delta t}^{p-1})$. In the error expansion, this value gets multiplied by ${\mathcal{O}}({\Delta t}^{k})$, so the leading error term is ${\mathcal{O}}({\Delta t}^{p})$ and the method has order $p-1$.

To prove the second part of Theorem 3, we view the solution given by (9) as an interpolant for the RK solution, evaluated at $\gamma_{n}{\Delta t}$. The order conditions for this interpolated solution are

[TABLE]

We see that the conditions (16b) and the first bushy-tree condition ((16a) with $k=1$) are still exactly fulfilled. The remaining bushy-tree conditions (16a) are not exactly satisfied; for $k\geq 2$ we have, using (14a) and (13),

[TABLE]

In the error expansion, each of these residuals is multiplied by ${\mathcal{O}}({\Delta t}^{k})$, so that the overall error incurred is ${\mathcal{O}}({\Delta t}^{(k-1)p+1})$. The exponent in this expression is at least $p+1$ for $k\geq 2$. ∎

It turns out that $\gamma_{n}$ satisfies the condition required in Theorem 3.

Lemma 4**.**

Let $a_{ij},b_{j}$ denote the coefficients of a Runge–Kutta method (6) of order $p$, let $f$ be a sufficiently smooth function, and let $\gamma_{n}$ be defined by (12). Then

[TABLE]

We illustrate Lemma 4 in Figure 1, which shows the convergence of $\gamma_{n}$ to 1 as the step size is reduced for the problem (28). Due to the symmetry of the problem, even faster convergence is observed for some methods.

Remark 5**.**

Lemma 4 was proved in [8, Proposition 4] for the special case of a certain family of four-stage, fourth-order methods.

Before proving Lemma 4 we state the main consistency result, which follows immediately from Theorem 3 and Lemma 4.

Corollary 5**.**

Let $(A,b)$ be the coefficients of a Runge–Kutta method of order $p$, and consider the RRK method defined by (6a) and (9) with $\gamma_{n}$ defined by (12).

• •

(IDT method) If the solution $u^{n+1}_{\gamma}$ is interpreted as an approximation to $u(t_{n}+{\Delta t})$, the method has order $p-1$.

• •

(RRK method) If the solution $u^{n+1}_{\gamma}$ is interpreted as an approximation to $u(t_{n}+\gamma_{n}{\Delta t})$, the method has order $p$; i.e. the solution after one step satisfies

[TABLE]

Remark 6**.**

The first part of Corollary 5 was proved by different means in [5, Thm. 2.1 (ii)].

In order to prove Lemma 4 we use the theory of B-series and follow the notation of [13]. In the remainder of this section, $t$ denotes a tree rather than time and $\gamma(t)$ denotes the density of a tree rather than the relaxation step length (see [13, Dfn. 2.10]). Let $\rho(t)$ denote the order of tree $t$, and let $t^{\prime}$ denote the tree obtained by attaching a new root node to the root of $t$. Let $[t_{1},t_{2}]$ denote the tree obtained by attaching a new root node to the root nodes of $(t_{1},t_{2})$ as in [13, Dfn. 2.12]. Then

[TABLE]

and

[TABLE]

Finally, let $\Phi_{j}(t)$ be defined as in [13, Dfn. 2.9]. We recall [13, Thm. 2.13]:

Theorem 6**.**

A Runge–Kutta method is of order $p$ iff

[TABLE]

for all trees of order $\leq p$.

Proof of Lemma 4.

We can write $\gamma=1-\delta/\eta$ where $\eta={\mathcal{O}}(1)$ and

[TABLE]

Thus it suffices to show that

[TABLE]

From [13, Thm. II.2.11] we have the Taylor series for the $i$th stage derivative:

[TABLE]

where $F(t)$ is the elementary differential corresponding to $t$. Thus $\langle f_{i},f_{j}\rangle$ can be expressed as a linear combination of inner products of elementary differentials:

[TABLE]

where $\rho(t)$ is the order (number of nodes) of tree $t$, and $t_{1},t_{2}$ range over the set of all labelled rooted trees. Here

[TABLE]

Because of the symmetry of the inner product, it is sufficient to show that

[TABLE]

for all pairs of trees $(t_{1},t_{2})$ satisfying

[TABLE]

We have

[TABLE]

Here we have applied (17) in various places, using the fact that our method is of order $p$ and that (due to (19)) the trees in question all have order less than or equal to $p$. ∎

Remark 7**.**

It is possible to prove a generalization of Lemma 4 without the use of B-series; see [22]. The proof above is included here to show a direct approach that is very different from that used in [22].

2.2 Comparison with projection methods

A projection Runge–Kutta method [12, Section IV.4] consists of the traditional Runge–Kutta formula (6) followed by a projection step:

[TABLE]

where $\Phi$ is the projection direction and $\lambda$ is chosen so that $u_{\lambda}^{n+1}$ lies on a desired manifold. A relaxation Runge–Kutta method can be viewed as a projection method along the direction $\Phi=u^{n+1}-u^{n}$ with step length $\lambda=1-\gamma_{n}$. However, the existing formalism for such methods does not include the possibility of interpreting the new solution as an approximation at a time different from $t_{n}+{\Delta t}$. In [5], this projection perspective was applied and it was shown that the resulting method is indeed of order $p-1$ when the original RK method is order $p$. We follow the terminology of [5] and refer to this interpretation as the incremental direction technique, or IDT.

To properly write a relaxation method as a projection method we can write (1) as the equivalent autonomous system of $m+1$ ODEs $v^{\prime}=g(v)$ in the standard way with

[TABLE]

where $g(v)=f(v_{m+1},v_{1:m})$. The projection of the solution of this problem in the direction $\Phi=v^{n+1}-v^{n}$ requires also projecting the updated value of $t$ to obtain the value $t_{n}+\gamma_{n}{\Delta t}$.

3 Stability properties of explicit relaxation RK methods

All of the results in the previous section apply to general Runge–Kutta methods. In the rest of this work, we focus on explicit methods.

An important question not explicitly answered in the foregoing analysis is how large the step size ${\Delta t}$ can be taken in practice. Theorem 2 guarantees unconditional stability for RRK methods applied to any conservative problem, and guarantees stability for dissipative problems as long as the step size is small enough. Due to the overall explicit nature of the methods, we should not expect in either case to obtain accurate results using step sizes much larger than what the original explicit RK method allows; in this respect RRK methods behave similarly to rational Runge–Kutta methods, to which their form is also very similar [27, 11]. In this section we investigate linear stability of RRK methods and also study how the strong stability preserving (SSP) property is affected by the use of RRK methods.

The behavior of an RRK method may be quite challenging to analyze, since $\gamma$ depends in a nonlinear way on the method coefficients and the numerical solution itself. On the other hand, over a single step we can think of $\gamma$ as a fixed value close to unity, and study the properties of the Runge–Kutta method

[TABLE]

3.1 The stability function

The stability function of a Runge–Kutta method is

[TABLE]

where $e$ is vector with all entries equal to unity. Thus the stability function corresponding to one step of the relaxation method (22) is

[TABLE]

Letting $\alpha_{k}$ denote the coefficients of $R(z)$ for an explicit method:

[TABLE]

we have

[TABLE]

Let $S(A,b)\subset{\mathbb{C}}$ denote the region of absolute stability of RK method $(A,b)$. It turns out that $S(A,\gamma b)$ grows as $\gamma$ decreases.

Theorem 7**.**

Let $\gamma_{1},\gamma_{2}$ be given such that $0\leq\gamma_{1}\leq\gamma_{2}$. Then $S(A,\gamma_{2}b)\subseteq S(A,\gamma_{1}b)$.

Proof.

It is sufficient to show that if $|R_{\gamma_{2}}(z_{*})|\leq 1$ for some $z_{*}\in{\mathbb{C}}$ then $|R_{\gamma_{1}}(z_{*})|\leq 1$. Let $w=\sum_{k=1}^{s}\alpha_{k}z_{*}^{k}$; then clearly $|R_{\gamma_{2}}(z_{*})|=|1+\gamma_{2}w|\leq 1$ implies $|R_{\gamma_{1}}(z_{*})|=|1+\gamma_{1}w|\leq 1$. ∎

Theorem 7 is illustrated in Figure 2, which shows the stability region for relaxation RK methods with 2–4 stages and order equal to the stage number. Stability region boundaries are shown for different values of $\gamma$ in the range $[0.7,1.3]$. This is an exaggerated range compared to values of $\gamma$ that are used in practice; for practical values of $\gamma$, the change in the stability region is visually too small to notice. It is particularly interesting to note that small reductions in $\gamma$ lead to significantly enhanced stability along the imaginary axis. Note also that as $\gamma\to 0$, the RRK method tends to the identity map and $S(A,\gamma b)\to{\mathbb{C}}$.

Consideration of absolute stability along the imaginary axis leads to the so-called E-polynomial, which for explicit RK methods is

[TABLE]

Clearly the method is stable for $z=iy$ such that $E(y)\geq 0$. Direct calculation shows that for any RRK method of order at least $p$ (where $p\geq 2$), we have

[TABLE]

The leading terms up to $y^{p}$ (which vanish for a standard RK method) are positive for $0\leq\gamma<1$, so we see that for every method of order two or higher $S(A,\gamma b)$ includes a segment of the imaginary axis containing $z=0$ when $\gamma\in[0,1)$. This can be observed for instance in Figure 2(a), where the RK method is unstable over the whole imaginary axis, but for $\gamma<1$ the stability regions include part of the imaginary axis.

3.2 Strong stability preservation

Because the RRK method is only a small perturbation of the original RK method, desirable properties of the original method may remain in effect when using the RRK version. We illustrate this idea by studying strong stability preserving RRK methods.

In the following, ${\mathcal{C}}(A,b)$ denotes the SSP coefficient (or radius of absolute monotonicity) of the Runge–Kutta method with coefficients $(A,b)$. Recall that ${\mathcal{C}}(A,b)$ is equal to the largest value $r\geq 0$ such that the method $(A,b)$ is absolutely monotonic at $-r$ [9].

For $0\leq\gamma\leq 1$, $u_{\gamma}^{n+1}$ is a convex combination of $u^{n}$ and $u^{n+1}$, which implies that ${\mathcal{C}}(A,\gamma b)\geq{\mathcal{C}}(A,b)$ for $0\leq\gamma\leq 1$. For many methods, ${\mathcal{C}}(A,\gamma b)$ also does not decrease when $\gamma$ is taken a little larger than 1.

Lemma 8**.**

Let the RK method with coefficients $(A,b)$ be absolutely monotonic at $z=-r$. Then the method $(A,\gamma b)$ with $\gamma\geq 0$ is also absolutely monotonic at $z=-r$ iff $R_{\gamma}(-r)\geq 0$.

Proof.

The conditions for absolute monotonicity of method $(A,b)$ at $z$ can be written (see, e.g. [9, p. 211]

[TABLE]

Conditions (24a) do not depend on $b$, while the first condition of (24b) will hold for any positive multiple of $b$ if it holds for $b$. ∎

Theorem 9**.**

Given a RK method with coefficients $(A,b)$, we have ${\mathcal{C}}(A,\gamma b)\geq{\mathcal{C}}(A,b)$ for $0\leq\gamma\leq\gamma_{*}$, where

[TABLE]

Proof.

Combining Lemma 8 with a standard result on absolute monotonicity (see e.g. [16, Lemma 3.1]) we have that $(A,\gamma b)$ is absolutely monotonic on $[-{\mathcal{C}}(A,b),0]$ as long as $R_{\gamma}(-{\mathcal{C}}(A,b))\geq 0$. Since $R_{\gamma}(z)-1=\gamma(R(z)-1)$, we obtain the condition stated in the theorem. ∎

Values of $\gamma_{*}$ are given in Table 1 for some well-known SSP methods. For all of these methods, direct computation shows that ${\mathcal{C}}(A,\gamma b)={\mathcal{C}}(A,b)$ for $0\leq\gamma\leq\gamma_{*}$, while taking $\gamma>\gamma_{*}$ leads to a decrease in the SSP coefficient.

4 Numerical examples

For conservative systems, stability is guaranteed under any step size. However, one expects the accuracy to deteriorate significantly if $\gamma$ is not close to unity. In the dissipative case, it is possible that $\gamma$ becomes negative and stability is lost. A careful examination of the analysis in the previous sections suggests that step sizes on the order of what linearized stability analysis predicts to be stable should be acceptable.

The following Runge–Kutta methods will be used in the numerical experiments.

• •

SSPRK(2,2): Two-stage, second-order SSP method of [23].

• •

SSPRK(3,3): Three-stage, third-order SSP method of [23].

• •

SSPRK(10,4): Ten-stage, fourth-order SSP method of [15].

• •

RK(4,4): Classical four-stage, fourth-order method.

• •

BSRK(8,5): Eight-stage, fifth-order method of [2]. A fixed step size ${\Delta t}$ is used and the embedded method is not used.

All of these methods have non-negative weights. We refer to the relaxation version of a method by replacing “RK” with “RRK”; e.g. RRK(4,4).

4.1 Linear, skew-Hermitian system

With this example we investigate the behavior of RRK methods for linear skew-Hermitian problems. For concreteness, we consider the advection equation

[TABLE]

with periodic boundary conditions over a spatial interval of length $2\pi$ discretized in space by a Fourier spectral collocation method with $m$ points. This results in a linear, constant-coefficient system of ODEs:

[TABLE]

where $u_{j}(t)\approx U(x_{j},t)$ at evenly spaced points $x_{j}$ and $D$ is the $m\times m$ skew-Hermitian Fourier spectral differentiation matrix. Since $D$ is normal, the behavior of any RK method on this problem can be characterized simply in terms of the eigenvalues of $D$ and the corresponding eigenvectors, which are discrete Fourier modes:

[TABLE]

Let us express the initial data in terms of these modes:

[TABLE]

where $\hat{u}$ is the discrete Fourier transform (DFT) of $u$. Then the exact solution of the semi-discrete system is given by

[TABLE]

Thus the energy associated with each mode is constant in time. Applying a Runge–Kutta method we instead obtain the solution

[TABLE]

where $R(z)$ is again the stability function of the method. The energy associated with mode $\xi$ is modified by the factor $|R(i\xi{\Delta t})|$ at each step. The maximum stable step size is the value that guarantees $i\xi{\Delta t}\in S(A,b)$ for all $\xi$. This is just

[TABLE]

where $I(A,b)$ is the length of the method’s imaginary axis stability interval. In the following experiments, we use a step size ${\Delta t}=\mu{\Delta t}_{\max}$. Using a given (standard) RK method and $0\leq\mu\leq 1$, we have absolute stability and the energy of each mode decays. This is illustrated in Figure 3(a), where we solve with $m=128$ and plot the relative change in amplitude

[TABLE]

for each mode, for a range of time step sizes $\mu{\Delta t}_{\max}$. Due to symmetry we plot only the positive wavenumbers. This figure does not depend on the initial data but only on $|R(iy)|^{N}$, where $N$ is the total number of time steps taken. For larger values of $\mu$, the high wavenumber modes are strongly damped.

Now let us consider what happens when applying a relaxation Runge–Kutta method with $\gamma_{n}$ chosen according to (12) so that energy is conserved. At each step, the energy in mode $\xi$ is modified by the factor $|R_{\gamma}(i\xi{\Delta t})|$. If the initial data is chosen to consist of a single wavenumber $k_{\xi}$, then $\gamma_{n}$ will take a value such that $|R_{\gamma}(i\xi{\Delta t})|=1$, and the same value of $\gamma_{n}$ will be used at every step. For more general initial data, $\gamma_{n}$ is chosen precisely so that the change in energy when summed over all modes is zero. This value depends on the data, so $\gamma_{n}$ will be different at each step. Furthermore, at each step some modes will be diminished while others will grow. This is illustrated in Figure 3(b), which is analogous to Figure 3(a) but for the energy-conserving RRK(4,4) method, with initial data taken as white noise; i.e. $\hat{u}^{0}_{j}=e^{i\theta_{j}}$ where the phases $\theta_{j}$ are random. We see that high-wavenumber modes are again damped, especially for step sizes close to the stability limit. Meanwhile, some lower-wavenumber modes are amplified in order to preserve the total energy. If we instead take initial data that is reasonably well-resolved, such as

[TABLE]

the resulting amplification curves are nearly indistinguishable from those of the standard RK4 method (Figure 3(a)), as shown in Figure 3(c). This is because most of the energy is in the low-wavenumber modes, which are propagated fairly accurately by the standard RK4 method, so little compensation is need in order to conserve energy. In the latter figure we also plot the amplification for $\mu=1.02$, which is beyond the absolute stability limit. We see that in this case the RRK method greatly amplifies the highest-wavenumber mode.

In Figure 4, we compare the behavior of RK(4,4) and relaxation RK(4,4) (RRK(4,4)) for different values of $\mu$. For $\mu<1$, the methods give similar solutions. Unlike RK(4,4), RRK(4,4) remains stable for $\mu>1$. However, taking $\mu=1.016$ leads to highly inaccurate and oscillatory approximations.

Remark 8**.**

Even for the linearly unstable value $\mu=1.016$, we have found that $\gamma$ remains within less than $10^{-2}$ of unity and the solution never blows up. For much larger values ($\mu\gtrsim 1.25$) the value of $\gamma$ tends to zero after a few steps and the calculation is never completed.

4.2 A non-normal linear autonomous problem

Here we consider a linear autonomous problem $u^{\prime}(t)=Au(t)$ with non-normal right-hand side introduced by Sun & Shu [24]:

[TABLE]

This problem is dissipative, but — as shown in [24] — the classical 4th-order Runge–Kutta method RK(4,4) is not monotone for this problem, no matter how small one takes ${\Delta t}$. To provide a concrete example, we compute $R({\Delta t}A)$, where $R(z)$ is the stability polynomial of RK(4,4) and we choose step size ${\Delta t}=0.5$. The first singular value of the matrix $R({\Delta t}A)$ is approximately $1.001$. Taking a single step with RK(4,4) and initial condition equal to the first right singular vector thus leads to an increase in the energy, as shown in Figure 5. The increase is even larger for the step size 0.7. In contrast, RRK(4,4) preserves monotonicity with either step size.

4.3 Nonlinear oscillator

Here we consider the problem

[TABLE]

with analytical solution

[TABLE]

Although energy is conserved in the exact solution, the widely used SSPRK(3,3) method of Shu & Osher [23] produces a solution whose energy is monotonically increasing for every positive step size [19]. Similar behavior is observed for several other explicit RK methods we have tested, as shown in Figure 6(a). By applying the modification described in the present work, we obtain instead Figure 6(b), showing that energy is conserved up to roundoff error for all methods.

Figure 7 shows the convergence for each of the standard RK methods (solid lines) and its energy-conserving modification (dashed lines). We see that the relaxation method is in each case on par with or more accurate than the original method.

4.4 Burgers’ equation

We solve the inviscid Burgers’ equation

[TABLE]

on the interval $-1\leq x\leq 1$ with periodic boundary conditions and initial data

[TABLE]

with the flux-differencing discretization

[TABLE]

where $F_{i\pm 1/2}$ is the numerical flux, defined below. The spatial domain is discretized with 50 equally-spaced points. In the convergence tests below, this spatial discretization is held fixed while the time step is varied, in order to investigate only the temporal convergence.

4.4.1 Energy-conservative semi-discretization

We take the second-order accurate symmetric flux [26]

[TABLE]

which yields a conservative semi-discrete system. Because of the lack of numerical viscosity, the semi-discrete solution is not the vanishing-viscosity solution, but instead develops a dispersive shock.

Energy evolution up to $t=0.2$ is shown in Figures 8. As expected, standard RK methods all exhibit significant growth or dissipation of energy, while RRK methods preserve energy up to roundoff error. Convergence results at $t=0.03$ (just before shock formation) are shown in Figure 9, where we compare the IDT approach with the full relaxation approach. All RRK methods achieve the order of accuracy of the corresponding RK method, whereas for IDT methods the convergence rate is reduced by one, as predicted by Theorem 3.

4.4.2 Energy-dissipative semi-discretization

We obtain a dissipative system by adding a centered difference to the flux:

[TABLE]

The amount of dissipation is controlled by $\epsilon>0$. The scheme is still consistent with (29) since the amount of dissipation is proportional to $\Delta x$. We take $\epsilon=1/100$. With this dissipative flux, the solution develops a viscous shock.

Results are shown in Figures 10 and 11. What is most interesting is that applying the relaxation approach dramatically improves the numerical approximation of the global dissipation. Convergence results are similar to those obtained with the conservative flux.

5 Conclusions

The relaxation approach we have proposed seems to be a simple and effective way to make any Runge–Kutta method preserve conservation or dissipativity with respect to an inner-product norm. This can be extended to more general convex functionals; see [22]. While we have focused here on explicit methods, the technique applies to implicit methods as well. Like Runge–Kutta methods, relaxation Runge–Kutta methods automatically preserve linear first integrals. If the original method is equipped with an embedded error estimator or dense output formula, these can be used without modification (other than the rescaling of the step size when determining the dense output times). As we have shown, many SSP RK methods retain the same SSP coefficient when used as RRK methods. The linear stability properties of a method are slightly modified depending on the choice of $\gamma$, but in general the allowable step size for an RRK method is essentially the same as that allowed for the original RK method.

Together these properties make RRK methods an attractive choice for symmetric hyperbolic systems; in combination with entropy-stable spatial discretizations they give a fully-discrete, explicit scheme that is provably entropy stable (for quadratic entropies).

In several experiments we have compared results obtained by viewing $u_{\gamma}^{n+1}$ as an approximation to $u(t_{n}+{\Delta t})$ (the so-called incremental direction technique, or IDT) versus those obtained by viewing it as an approximation to $u(t_{n}+\gamma_{n}{\Delta t})$ (the relaxation approach). The main purpose of this comparison is to illustrate our theoretical convergence estimates; it is clear that in practice one should always use the latter interpretation. Comparison with embedded projection methods is planned as future work.

For an RK method with $s$ stages, determination of $\gamma_{n}$ via the formula (12) requires the evaluation of $s+1$ inner products (see Remark 4). For typical high-order discretizations of nonlinear PDEs this cost is negligible compared to the $s$ evaluations of $f$ required for each step. For simpler applications (such as linear wave equations) the cost may be an important factor. A careful comparison of RRK schemes relative to other time discretizations for linear wave equations is the subject of ongoing work.

6 Acknowledgments

The author is grateful to Hendrik Ranocha for helpful comments on drafts of this work. This work was supported by funding from King Abdullah University of Science & Technology.