DIRK Schemes with High Weak Stage Order

Runge-Kutta time-stepping methods in general suffer from order reduction: the observed order of convergence may be less than the formal order when applied to certain stiff problems. Order reduction can be avoided by using methods with high stage order. However, diagonally-implicit Runge-Kutta (DIRK) schemes are limited to low stage order. In this paper we explore a weak stage order criterion, which for initial boundary value problems also serves to avoid order reduction, and which is compatible with a DIRK structure. We provide specific DIRK schemes of weak stage order up to 3, and demonstrate their performance in various examples.


Introduction
Runge-Kutta (RK) methods achieve high-order accuracy in time by means of combining approximations to the solution at multiple stages. An s-stage RK scheme can be represented via the Butcher tableau c A b T = c 1 a 11 · · · a 1s . . . . . . . . . c s a s1 · · · a ss b 1 · · · b s . Throughout the whole paper we assume that c = A e, where e is the vector of all ones. The scheme's stability function [12] R(ζ) = 1 + ζ b T (I − ζA) −1 e measures the growth u n+1 /u n per step ∆t, when applying the scheme to the linear model equation u (t) = λu, with ζ = λ∆t.
A particular interest lies in the accuracy of the RK scheme for stiff problems, i.e., problems in which a larger time step is chosen than the fastest time scale of the problem's dynamics. A standard stiff model problem [8] is the scalar linear ordinary differential equation (ODE) with i.c. u(0) = φ(0) and Re λ ≤ 0. The true solution y(t) = φ(t) evolves on an O(1) time scale. Hence, λ-values with large negative real part result in stiffness. Considering a family of test problems (parametrized by λ), one can now establish the scheme's convergence via two different limits: (a) the non-stiff limit ∆t → 0 and ζ → 0; and (b) the stiff limit ∆t → 0 and ζ → −∞. A characteristic property of most RK schemes is that, while the non-stiff limit recovers the scheme's order (as given by the order conditions [2,5]), the error decays at a reduced order in the stiff limit. This phenomenon is called "order reduction" (OR) [11,10,7,3,1] and it manifests in various ways for more complex problems, including numerical boundary layers [6]. The OR phenomenon can be seen by studying the RK scheme applied to (1). The approximation error at time t n+1 reads [12, Chapter IV.15] where R(ζ) is the growth factor, and are the truncation errors incurred at the intermediate stages and at the end of the step, respectively. Here, φ (j) denotes the j-th derivative of the solution, and the vectors . . we call the stage order residuals or stage order vectors. The condition τ (η) = 0 for 0 ≤ η ≤ j appears often in the literature and is also referred to as the simplifying assumption C(η) [12]. In (2), the step error δ n+1 is of the formal order (in ∆t) of the scheme (due to the order conditions). Moreover, the growth factor carries over (more or less, see [4]) the accuracy from one to the next step. Hence, the critical expression for OR is the term involving the stage error δ n+1 s . Specifically, the asymptotic behavior of the expression matters. In the non-stiff limit (ζ 1), a Neumann expansion yields ζ(I −ζA) −1 = ζI +ζ 2 A+ ζ 3 A 2 + . . . , leading to expressions b T A τ (j) with > 0. And in fact the order conditions guarantee that b T A τ (j) = 0 for 0 ≤ + j ≤ p − 1 to ensure the formal order of the scheme. Conversely, in the stiff limit we can treat ζ −1 as the small parameter and expand ζ( . , leading to expressions b T A τ (j) with < 0. The order conditions do not imply that these quantities vanish, and in general one may observe a reduced rate of convergence.
A key question is therefore whether additional conditions can be imposed on the RK scheme that recover the scheme's order in the stiff regime. A well-known answer to the question is: Definition 1. Letp denote the order of the quadrature rule of an RK scheme. Letq denote the largest integer such that τ (j) = 0 for 1 ≤ j ≤q. The stage order of a RK scheme is q = min(p,q).
Having stage order q implies that the error decays at an order of (at least) q in the stiff regime (see also [12]). This work focuses particularly on diagonally-implicit Runge-Kutta (DIRK) schemes, for which A is lower diagonal. A known drawback of DIRK schemes is that they cannot have high stage order: Theorem 1. The stage order of an irreducible DIRK scheme is at most 2. The stage order of a DIRK scheme with non-singular A is at most 1.
Proof. Since c = A e, we have τ Thus if A is non-singular, one has τ (2) = 0, so q ≤ 1. Consider now the case that a 11 = c 1 = 0, and suppose that the method has stage order 3. The conditions τ Hence, while DIRK schemes possess an implementation-friendly structure (each stage is a backward-Euler-type solve), their potential to avoid OR by means of high stage order is limited. We therefore move to a weaker condition that can avoid OR in some situations for higher order in the context of DIRK schemes.

Weak Stage Order
To avoid order reduction, the expressions g (j) in (3) need to vanish in the stiff limit. In line with [9], we define the following criteria: Definition 2 (weak stage order). A RK scheme has weak stage order (WSO)q if there is an A-invariant subspace that is orthogonal to b and that contains the stage order vectors τ (j) for 1 ≤ j ≤q.
Theorem 2. (WSO is the most general condition that ensures g (j) = 0 for all ζ > 0) Let coefficients A, b be given. Then g (j) = 0 for all ζ > 0 and 1 ≤ j ≤q if and only if the corresponding RK scheme has weak stage orderq.
Proof. Let C(G) denote the column space of

From the Cayley-Hamilton theorem it follows that WSOq is equivalent to
Differentiating both sides of this equation -times, with respect to ζ, and taking the limit as ζ → 0 + , yields the conditions in equation (4).
Definition 3 (weak stage order eigenvector criterion). A RK scheme satisfies the WSO eigenvector criterion of orderq e if for each 1 ≤ j ≤q e , there exists µ j such that A τ (j) = µ j τ (j) , and moreover, b T τ (j) = 0.
The WSO eigenvector criterion of orderq e implies WSO (of at least)q e . For a given scheme, let p denote the classical order, q the stage order, andq the weak stage order. Then we havẽ q ≥ q and p ≥ q. Note however that a method with WSOq ≥ 1 need not even be consistent; order conditions must be imposed separately.
The WSO eigenvector criterion may serve to avoid OR because it implies that i.e., it allows one to "push" the stage order residuals past the matrix (1 − ζA) −1 , and then use b T τ (j) = 0. Note that the condition b T τ (j) = 0 that is required in Def. 3 is actually automatically satisfied (due to the order conditions) if p >q e (or p ≥q e for stiffly accurate schemes).
It must be stressed that the concept of WSO (both criteria) is based on the linear test equation (1), hence it is not clear to what extent WSO will remedy OR for nonlinear problems or problems with time-dependent coefficients. In Section 4 we numerically investigate some nonlinear test problems.
Finally, we present a limitation theorem on the WSO eigenvector criterion.
Proof. Because the τ (j) only depend on A, the eigenvector relation in Def. 3 depends only on A, not on b. With A lower triangular, the first k components of τ (j) depend only on the upper k rows of A; and the same is true for the eigenvector relation as well. Hence, for a scheme to have an A that allows for the WSO eigenvector criterion of orderq e , all upper sub-matrices of A must admit the same, too. We can therefore study A row by row. The first component of τ (j) equals (1 − 1 j )a j 11 , which is nonzero for j > 1. Hence, the first row of the equation A τ (j) = µ j τ (j) is equivalent to µ j = a 11 . With that, we can move to the second row of the equation, which reads To determine the set of solutions (a 11 , a 21 , a 22 ) of (5), we first observe that (5) is homogeneous, i.e., if (a 11 , a 21 , a 22 ) solves (5), then (µa 11 , µa 21 , µa 22 ) solves (5) as well for any µ ∈ R. It therefore suffices to consider the solutions of (5) in the 2D-plane ( a 11 a 21 , a 22 a 21 ). Figure 1 shows the resulting solution curves for j ∈ {2, 3, 4}.
One class of solutions lies on the straight line of slope 1 passing through (1, 0). Those schemes are equal-time methods, i.e., RK schemes that have c = ν e, where ν ∈ R is a constant. In fact, equal-time schemes satisfy the eigenvector relation for all j. However, they are not particularly useful RK methods, because-among other limitations-they are restricted to second order. This follows because the order 1 and 2 conditions require b T e = 1 and b T c = 1 2 . Thus ν = 1 2 , and b T c 2 = ν 2 = 1 4 , which contradicts the order 3 condition b T c 2 = 1 3 . Note that the equal-time scenario also covers the points at infinity in Fig. 1, i.e., the schemes with a 21 = 0.
Non-equal-time schemes that satisfy (5) for j = 2 and j = 3 are the following two points in the ( a 11 a 21 , a 22 a 21 ) plane: None of these two points satisfies (5) for j = 4 (green curve in Fig. 1). Thereforeq e ≤ 3.
Among the two sets of solutions found in the proof, P 1 implies that a 11 , a 21 , and a 22 all have the same sign, which is a desirable property. In contrast, P 2 implies that a 21 < 0. Both WSO 3 schemes presented below correspond to the P 1 solution.

DIRK Schemes With High Weak Stage Order
Imposing the classical order conditions [2,5], together with the WSO eigenvector relation (Def. 3), we determine RK schemes by searching the parameter space of DIRK schemes (with all diagonal entries non-zero). A stiffly accurate structure ( b T equals the last row of A) is imposed, as is A-stability (verified by evaluating the stability function R(ζ) along the imaginary axis). Together this implies that the resulting scheme is L-stable; i.e., it ensures that unresolved stiff modes decay [5]. The number of stages is chosen so that the constraints admit solutions. The optimization itself is carried out using MATLAB's optimization toolbox, using multiple local optimization algorithms included in the function fmincon. An effort was made to minimize the L 2 norm of the local truncation error coefficients. However, in multiple cases the solver exhibited bad convergence properties; so while the schemes below yield reasonable truncation errors, it should not be expected that they are optimal. We find an order 3 scheme with WSO 2 (see also [9]

Numerical Results
In this section we verify the order of accuracy of the schemes above and demonstrate that WSO remedies order reduction for linear problems. We confirm that WSO p is required for ODEs, and WSO p − 1 is required for PDE IBVPs. In addition, we study the effect of WSO for two nonlinear problems.

4.1.
Linear ODE test problem. We consider the linear ODE test problem (1) with the true solution φ(t) = sin(t + π 4 ), the stiffness parameter λ = −10 4 , and the initial condition u(0) = sin( π 4 ). The problem is solved using three 3rd order DIRK schemes (with WSO 1, 2, and 3) and two 4th order DIRK schemes (with WSO 1 and 3) 1 up to the final time T = 10. The convergence results are shown in Fig. 2. In the stiff regime where |ζ| = |λ|∆t 1, first order convergence is observed for the WSO 1 schemes as expected, the WSO 2 scheme improves the convergence rate to 2, and the WSO 3 schemes exhibit 3rd order convergence. In addition to yielding better convergence orders in the stiff regime, the schemes with higher WSO also turn out to yield substantially smaller error constants in the non-stiff regime (∆t 1/|λ|). For comparison, we also display a DIRK scheme with explicit first stage (EDIRK), that is, a 11 = 0, of stage order 2 (see Thm. 1). The left panel of Fig. 2 shows that the WSO 2 scheme exhibits the same convergence behavior as the stage order 2 EDIRK scheme and performs equally well in terms of accuracy.

4.2.
Linear PDE test problem: Schrödinger equation. As a linear PDE test problem, we study the dispersive Schrödinger equation. The method of manufactured solutions is used, i.e., the forcing, the boundary conditions (b.c.) and initial conditions (i.c.) are selected to generate a desired true solution. The spatial approximation is carried out using 4th order centered differences on a fixed spatial grid of 10000 cells. This renders spatial approximation errors negligible and thus isolates the temporal errors due to DIRK schemes. The errors are measured in the maximum norm in space.  We consider with the true solution u(x, t) = e i(kx−ωt) , ω = 2π and k = 5. Figure 3 shows the convergence orders of u, u x and u xx for 3rd order DIRK schemes with WSO 1 (left), WSO 3 (middle) and a 4th order DIRK scheme with WSO 3 (right). For IBVPs, spatial boundary layers are produced by RK methods, thus limiting the convergence order in u toq + 1, with an additional half an order loss per derivative whenq < p [9]. As a result, the 4th order WSO 3 scheme recovers 4th order convergence in u and improves the convergence in u x and u xx . Whenq = p, the full convergence order in u, u x and u xx is achieved, as seen in the middle panel in Fig. 3 Here ν = 0.1 and u(x, t) = cos(2+10t) sin(0.2+20x). The nonlinear implicit equations arising at each time step are solved using a standard Newton iteration. The choice of Neumann b.c. distinguishes this example from the one given in [9]. With Neumann b.c., the convergence order in u is limited toq + 1.5 (half an order better than with Dirichlet b.c.). Figure 4 shows that order reduction arises with the stage order 1 scheme, and that the WSO2 scheme recovers 3rd order convergence for u and u x , and the 3rd order WSO 3 scheme yields 3rd order convergence for u, u x and u xx .
with i.c. (x(0), y(0)) = (2, 0), stiffness parameter µ = 500, and final time T = 10. The nonlinear system at each time step is solved via MATLAB's built-in nonlinear system solver. The "exact" solution is computed using explicit RK4 with a time step ∆t = 10 −6 . In this case, the presented DIRK schemes with high WSO do not improve the convergence rates in the stiff regime and they perform worse than the WSO 1 scheme in terms of accuracy (see Fig. 5). On the other hand, an EDIRK with stage order 2 improves the rate of convergence in the stiff regime (see right panel in Fig. 5). However, it does so, interestingly, by yielding larger errors for large time steps.

Conclusions and Outlook
This study demonstrates that it is possible to overcome order reduction (OR) for certain classes of problems in the context of DIRK schemes, even though these are limited to low stage order. A specific weak stage order (WSO) "eigenvector" criterion has been presented, analyzed, and applied to determine DIRK schemes with WSO up to 3. The numerical results confirm that the schemes avoid OR for linear problems and for some nonlinear problems  Figure 5. Error convergence for Van der Pol's equation. Left: 3rd order DIRK schemes with WSO 1 (blue circles), WSO 2 (red triangles) and WSO 3 (black squares). Right: 4th order DIRK schemes with WSO 1 (blue circles) and WSO 3 (red triangles), and a 3rd order EDIRK scheme with stage order 2 (black squares).
in which the mechanism for order reduction is linear (i.e., boundary conditions). The key limitation found herein is that the eigenvector criterion cannot go beyond WSO 3 for DIRK schemes. Hence, a key question of future research is how high WSO is admitted by the general criterion in Def. 2. Another important future research task is to devise further DIRK schemes that are truly optimized in terms of truncation error coefficients or other criteria.