Some Notes on Summation by Parts Time Integration Methods

Some properties of numerical time integration methods using summation by parts operators and simultaneous approximation terms are studied. These schemes can be interpreted as implicit Runge-Kutta methods with desirable stability properties such as $A$-, $B$-, $L$-, and algebraic stability. Here, insights into the necessity of certain assumptions, relations to known Runge-Kutta methods, and stability properties are provided by new proofs and counterexamples. In particular, it is proved that a) a technical assumption is necessary since it is not fulfilled by every SBP scheme, b) not every Runge-Kutta scheme having the stability properties of SBP schemes is given in this way, c) the classical collocation methods on Radau and Lobatto nodes are SBP schemes, and d) nearly no SBP scheme is strong stability preserving.


Known Results on SBP SAT Schemes
In order to solve an ordinary differential equation (ODE) ∀t ∈ (0, T) : u (t) f (t, u(t)), a grid 0 ≤ τ 1 < · · · < τ s ≤ T is introduced and the numerical solution is approximated pointwise as u i u(τ i ) and f i f (τ i , u i ). Summation by parts (SBP) operators can be defined as follows, cf. [2,3,10]. • a symmetric and positive definite discrete mass/norm matrix M approximating the L 2 scalar product u T Mv ≈ ∫ T 0 u(t)v(t) dt, • and interpolation vectors t L , t R approximating the interpolation to the boundary as t T L u ≈ u(0), t T R u ≈ u(T) with order of accuracy at least p, such that SBP operators mimic integration by parts discretely via the summation by parts property (2). An SBP time discretisation using a simultaneous approximation term (SAT) of (1) with parameter σ ∈ R is [1,7,8] Du f + σM −1 t L u 0 − t T L u .
Most stability results have been achieved for the choice σ 1, i.e.
Hence, this discretisation will be considered in the following. The numerical solution at t T is given by t T R u, where u solves (4). The interval [0, T] can also be partitioned into multiple subintervals/blocks such that multiple steps of this procedure are used sequentially.
In order to guarantee that (4) can be solved for a dissipative linear scalar problem, the following assumption is introduced [8]. Assumption 1.2. For σ > 1 2 , all eigenvalues of D + σM −1 t L t T L have strictly positive real part. The following characterisation of (4) as Runge-Kutta method has been developed in [1]. Theorem 1.3. If assumption 1.2 is satisfied, (4) is equivalent to an implicit Runge-Kutta method with the following Butcher coefficients, where 1 denotes also the vector (1, . . . , 1) T ∈ R s .
The factor 1 T is needed since Runge-Kutta coefficients are normalised to the interval [0, 1]. In order to make this article sufficiently self-contained, some classical stability properties of Runge-Kutta methods will be recalled briefly, cf.
Re λ ≤ 0. The numerical solution after one time step of a Runge-Kutta method with Butcher is the stability function of the Runge-Kutta method. The stability property is mimicked discretely as |u + | ≤ |u 0 | if R(λ ∆t) ≤ 1. Definition 1.4. A Runge-Kutta method with stability function R(z) is A-stable, if R(z) ≤ 1 for all z ∈ C with Re(z) ≤ 0. The method is L-stable, if it is A-stable and lim z→∞ R(z) 0.
Hence, A-stable methods are stable for every time step ∆t > 0 and L-stable methods damp out stiff components corresponding to λ −x with large x ∈ R sufficiently fast.
Another classical stability property is connected with possibly nonlinear problems (1) in Hilbert spaces satisfying a one-sided Lipschitz condition ∀t, u, v : where ν ∈ R is the one-sided Lipschitz constant of f . This condition gives some bounds on the growth rate of the difference between two solutions. In particular, the distance between two solutions cannot increase if ν ≤ 0.
The following stability properties have been obtained in [1,7]. Theorem 1.6. Suppose that assumption 1.2 holds. Then, the SBP SAT scheme (4) is A-and L-stable. If the mass matrix M is diagonal, the scheme is also B-stable.

Assumptions and Algebraic Stability
In this section, the new results of this short note concerning the necessity of assumption 1.2 and the necessity of an SBP SAT form for stability properties guaranteed by Theorem 1.6 are presented.

Assumption on Eigenvalues of
Assumption 1.2 has been proved for classical second order SBP operators in [8] and for SBP operators on Gauss, Radau, and Lobatto quadrature nodes in [9]. It has been examined numerically for other classical finite difference SBP operators in [8]. Since assumption 1.2 holds for all known SBP SAT schemes investigated in [1,[7][8][9], it is interesting to know whether it follows from properties of SBP operators. Theorem 2.1. There are SBP operators that do not satisfy assumption 1.2.
Proof. Consider the operators on the uniform grid with four nodes 0, 1 3 (2) is satisfied, t L and t R are exact, and D is a first order accurate SBP derivative operator. However,

Algebraic Stability
Many stability properties such as Aand B-stability are satisfied if the following algebraic criterion is fulfilled by the coefficients of a Runge-Kutta method [5,Theorem 12.4].

Definition 2.2. A Runge-Kutta method with Butcher coefficients
It has been noted in [1] that an SBP SAT scheme (4) with diagonal M is algebraically stable, since the nodes τ i are pairwise distinct, i.e. the corresponding Runge-Kutta method is nonconfluent. In that case, Band algebraic stability are equivalent [5,Corollary 12.14]. This can also be proved directly, cf. [1,Theorem 5.8].
It is interesting to know whether all Runge-Kutta methods with stability properties guaranteed by Theorem 1.6 can be constructed as SBP SAT schemes. Since those schemes are L-stable, the classical Gauss collocation schemes (which are not L-stable) cannot be constructed in this way, cf. [1]. However, there is ii) The Runge-Kutta method is given via Theorem Then, the algebraic stability matrix diag(b)A+A T diag(b)−bb T has the eigenvalues 5 8 , 3 8 , and zero (twofold). Hence, the Runge-Kutta method is algebraically stable (because b i > 0 is satisfied additionally) and therefore also Aand B-stable. Its stability functions fulfils lim z→∞ R(z) 0. Thus, the scheme is also L-stable.
It suffices to consider T 1. If the scheme is given by an SBP SAT method (4) Similarly, consistency of D and t L implies t R defined by (11) is first order accurate, i.e. t T R 1 1 and t T R c 1. The same accuracy of t L requires t T L 1 1, t T L c 0.
Because of (12), D can be written as Since M ∈ R 4×4 should be symmetric, it is determined by ten real parameters, e.g.

Classical Collocation Methods
In [1], it has been shown that the SBP SAT scheme with Lobatto quadrature on four nodes corresponds to the classical Lobatto IIIC method with s 4. It has been mentioned that this is similar for the Radau IA and Radau IIA schemes. However, to the authors knowledge, no general proof of this result has been given up to now. To prove it, the classical conditions will be used.
where the exponentiation c q is performed pointwise. Inserting A from (5) yields This is equivalent to where v is any polynomial of degree ≤ s − 1, evaluated at the nodes c i . Since the left endpoint 0 The Radau quadrature is exact for polynomials of degree ≤ 2s − 2. Hence, for every q ∈ {1, . . . , s}, and (using integration by parts) proving D(s).
The weights and nodes of the right Radau quadrature (right endpoint 1 included) are the weights b i and nodes c i of the Radau IIA method. The matrix A of the Radau IIA method is determined uniquely by the condition C(s), i.e. C(η) with η s in (15) [5, section IV.5]. Hence, it suffices to prove that the SBP SAT method satisfies C(s), which can be written using M diag(b) as where the exponentiation c q is again performed pointwise. Inserting A from (5), this is equivalent to where v is any polynomial of degree ≤ s − 1, evaluated at the nodes c i . Using the SBP property (2), this can be rewritten as Since the right endpoint 1 is included, Using the exactness of the Radau quadrature for polynomials of degree ≤ 2s − 2, for every q ∈ {1, . . . , s}, and (using integration by parts) proving C(s). Finally, the weights and nodes of the Lobatto quadrature (left and right endpoints 0, 1 included) are the weights b i and nodes c i of the Lobatto IIIC method. The matrix A of the Lobatto IIIC method is determined uniquely by the condition C(s − 1) and a i,1 b 1 , i ∈ {1, . . . , s} [5,section IV.5]. Since the order of accuracy of the SBP operator is s − 1, C(s − 1) is satisfied [1,Lemma 5.3]. This can also be proved using similar manipulations as above. Hence, it remains to show a i,1 b 1 , i ∈ {1, . . . , s}. Since D is exact for constants, t L (1, 0, . . . , 0) T , and M diag(b 1 , . . . , b s ), Therefore, (a i,1 ) s i 1 Another desirable stability property of time integration methods is that they are strong stability preserving (SSP), i.e. that they preserve convex stability properties of the explicit Euler method [4]. Definition 4.1. A numerical time integration method is called strongly stable for a given convex functional η if η(u + ) ≤ η(u 0 ), possibly using some time step restriction of the form 0 < ∆t ≤ ∆t max . A numerical time integration method is called strong stability preserving with SSP coefficient c > 0, if η(u + ) ≤ η(u 0 ) for all time steps 0 < ∆t ≤ c ∆t E whenever the explicit Euler method is strongly stable for the convex functional η and time steps 0 < ∆t ≤ ∆t E .
Typical convex functionals η considered for SSP methods are the norm in a Hilbert space for dissipative operators or the total variation seminorm for semidiscretisations of scalar conservation laws.
If a i, j were non-negative, the left hand side would be non-positive for i 1 (since c j ≥ c 1 ) and thus zero. Hence, the first row of A would be zero, which is impossible, because A is invertible. If the SBP operator is at least first order accurate, the corresponding Runge-Kutta method satisfies C(1) and D(1) [1, Lemma 5.3 and Lemma 5.4], i.e.

Remark 4.3.
Classical finite difference SBP operators and those based on Radau or Lobatto quadrature include at least one endpoint and can thus not result in SSP schemes. The SBP SAT scheme (4) on two Gauss nodes does not contain an endpoint and has a first order accurate derivative operator. Nevertheless, the scheme is not SSP, since the corresponding matrix A has a negative entry. Example 4.4. There is a first order accurate SBP operator with diagonal norm matrix not including any boundary node such that the resulting Runge-Kutta method given by Theorem 1.3 is The operators D, t L , t R are exact for polynomials of degree one, assumption 1.2 has been verified numerically for σ ∈ (1/2, 2), A and b have only non-negative entries, and the scheme is strong stability preserving with SSP coefficient ≈ 1.35, computed using NodePy [6].