Based on the fact that integration by parts plays a major role in the development of energy and entropy estimates for initial boundary value problems, one may conjecture that the summation by parts (SBP) property [svard2014review, fernandez2014review] is a key factor in provably stable schemes. Although it is complicated to formulate such a conjecture mathematically, there are several attempts to unify stable methods in the framework of summation by parts schemes, starting from the origin of SBP operators in finite difference methods [kreiss1974finite, strand1994summation] and ranging from finite volume [nordstrom2001finite, nordstrom2003finite] and discontinuous Galerkin methods [gassner2013skew] to flux reconstruction schemes [ranocha2016summation].
Turning to SBP methods in time [nordstrom2013summation, lundquist2014sbp, boom2015high], a class of linearly and nonlinearly stable SBP schemes has been constructed and studied in this context, see also [ruggiu2018pseudo, ruggiu2019eigenvalue, linders2019convergence]. If the underlying quadrature is chosen as Radau or Lobatto quadrature, these Runge–Kutta schemes are exactly the classical Radau IA, Radau IIA, and Lobatto IIIC methods [ranocha2019some]. Having the conjecture “stability results require an SBP structure” in mind, this article provides additional insights to this topic by constructing new classes of SBP schemes, which reduce to the classical Lobatto IIIA and Lobatto IIIB methods if that quadrature rule is used. Consequently, all stable classical Runge–Kutta methods based on Radau and Lobatto quadrature rules can be formulated in the framework of SBP operators. Notably, instead of using simultaneous approximation terms (SATs) [carpenter1994time, carpenter1999stable] to impose initial conditions weakly, these new schemes use a strong imposition of initial conditions in combination with a projection method [olsson1995summationI, olsson1995summationII, mattsson2018improved].
By mimicking integration by parts at a discrete level, the stability of SBP methods can be obtained by mimicking the continuous analysis. All known SBP time integration methods are implicit and the stability does not depend on the size of the time step. In contrast, the stability analysis of explicit time integration methods can also use techniques similar to summation by parts, but the analysis is in general more complicated and restricted to sufficiently small time steps [ranocha2018L2stability, sun2017stability, sun2019strong]. Since there are strict limitations for stability of explicit methods, especially for nonlinear problems [ranocha2019strong, ranocha2019energy], an alternative to stable fully implicit methods is to modify less expensive (explicit or not fully implicit) time integration schemes to get the desired stability results [ketcheson2019relaxation, ranocha2019relaxation, ranocha2020general, ranocha2020relaxationHamiltonian, offner2018artificial, sun2019enforcing].
This article is structured as follows. At first, the existing class of SBP time integration methods is introduced in Section 2, including a description of the related stability properties. Thereafter, the novel SBP time integration methods are proposed in Section 3. Their stability properties are studied and the relation to Runge–Kutta methods is described. In particular, the Lobatto IIIA and Lobatto IIIB methods are shown to be recovered using this framework. Afterwards, results of numerical experiments demonstrating the established stability properties are reported in Section 4. Finally, the findings of this article are summed up and discussed in Section 5.
2 Known Results for SBP Schemes
Consider an ordinary differential equation (ODE)
with solution in a Hilbert space. Summation by parts schemes approximate the solution on a finite grid pointwise as and . The SBP operators can be defined as follows, cf. [svard2014review, fernandez2014review, fernandez2014generalized].
A first derivative SBP operator of order on consists of
a discrete operator approximating the derivative with order of accuracy ,
a symmetric and positive definite discrete quadrature matrix approximating the scalar product ,
There are analogous definitions of SBP operators for second or higher order derivatives. In this article, only first derivative SBP operators are considered.
The quadrature matrix is sometimes called norm matrix (since it induces a norm via a scalar product) or mass matrix (in a finite element context).
Because of the SBP property (2), SBP operators mimic integration by parts discretely via
However, this mimetic property does not suffice for the derivations to follow. Nullspace consistency will be used as an additional required mimetic property. This novel property was introduced in [svard2019convergence] and has been a key factor in [linders2019convergence, ranocha2020discrete].
A first derivative SBP operator is nullspace consistent, if the nullspace (kernel) of satisfies .
Here, denotes the discrete grid function with value unity at every node.
Every first derivative operator (which is at least first order accurate) maps constants to zero, i.e. . Hence, the kernel of always satisfies . Here and in the following, denotes the subspace relation of vector spaces. If is not nullspace consistent, there are more discrete grid functions besides constants which are mapped to zero (which makes it inconsistent with ). Then, and undesired behavior can occur, cf. [svard2019convergence, linders2019convergence, ranocha2019some].
An SBP time discretization of (1) using a simultaneous approximation term (SAT) to impose the initial condition weakly is [nordstrom2013summation, lundquist2014sbp, boom2015high]
The numerical solution at is given by , where solves (4).
The interval can be partitioned into multiple subintervals/blocks such that multiple steps of this procedure can be used sequentially [lundquist2014sbp].
In order to guarantee that (4) can be solved for a linear scalar problem, must be invertible, where is a real parameter usually chosen as . The following result has been obtained in [linders2019convergence, Proposition 4].
If is a first derivative SBP operator, is invertible if and only if is nullspace consistent.
In [ruggiu2018pseudo, ruggiu2019eigenvalue], it was explicitly shown how to prove that is invertible in the pseudospectral/polynomial and finite difference case.
As many other one-step time integration schemes, SBP-SAT schemes (4) can be characterized as Runge–Kutta methods, given by their Butcher coefficients [hairer2008solving, butcher2016numerical]
where and . For (1), a step from to is given by
Here, are the stage values of the Runge-Kutta method. The following characterization of (4) as Runge–Kutta method was given in [boom2015high].
Consider a first derivative SBP operator . If is invertible, (4) is equivalent to an implicit Runge–Kutta method with the Butcher coefficients
The factor is needed since the Butcher coefficients of a Runge–Kutta method are normalized to the interval .
Next, we recall some classical stability properties of Runge–Kutta methods for linear problems, cf. [hairer2010solving, Section IV.3]. The absolute value of solutions of the scalar linear ordinary differential equation (ODE)
cannot increase if . The numerical solution after one time step of a Runge–Kutta method with Butcher coefficients is , where
is the stability function of the Runge–Kutta method. The stability property of the ODE is mimicked discretely as if .
A Runge–Kutta method with stability function for all with is stable. The method is stable, if it is stable and .
Hence, stable methods are stable for every time step and stable methods damp out stiff components as .
The following stability properties have been obtained in [lundquist2014sbp, boom2015high].
Consider a first derivative SBP operator . If is invertible, then the SBP-SAT scheme (4) is both and stable.
The SBP-SAT scheme (4) is both and stable if is a nullspace consistent SBP operator.
3 The New Schemes
The idea behind the novel SBP time integration scheme introduced in the following is to mimic the reformulation of the ODE (1) as an integral equation
Taking the time derivative on both sides yields . The initial condition is satisfied because . Hence, the solution of (1) can be written implicitly as the solution of the integral equation (10). Note that the integral operator is the inverse of the derivative operator with a vanishing initial condition at . Hence, a discrete inverse (an integral operator) of the discrete derivative operator with a vanishing initial condition will be our target.
In the space of discrete grid functions, the scalar product induced by is used throughout this article. The adjoint operators with respect to this scalar product will be denoted by , i.e. . The adjoint of a discrete grid function is denoted by .
By definition, the adjoint operator of satisfies
for all grid functions . The adjoint is a discrete representation of the inverse Riesz map applied to a grid function [roman2008advanced, Theorem 9.18] and satisfies
The following lemma and definition were introduced in [ranocha2020discrete].
For a nullspace consistent first derivative SBP operator , .
A fixed but arbitrarily chosen basis vector of for a nullspace consistent SBP operator is denoted as .
The name is intended to remind the reader of (grid) oscillations, since the kernel of is orthogonal to the image of [roman2008advanced, Theorem 10.3] which contains all sufficiently resolved functions. Several examples are given in [ranocha2020discrete]. To prove , choose any and and compute
Consider the SBP operator of order defined by the Lobatto Legendre nodes and in . Then,
and , where . Here, represents the highest resolvable grid oscillation on and is orthogonal to , since .
The following technique has been used in [ranocha2020discrete] to analyze properties of SBP operators in space. Here, it will be used to create new SBP schemes in time. Consider a nullspace consistent first derivative SBP operator on the interval using grid points and the corresponding subspaces
Here and in the following, denotes the value of the discrete function at the initial time . For example, is the first coefficient of if and .
is the vector space of all grid functions which vanish at the left boundary point, i.e. . is the vector space of all grid functions which can be represented as derivatives of other grid functions, i.e. is the image of .
From this point in the paper, denotes a nullspace consistent first derivative SBP operator.
The mapping is bijective, i.e. one-to-one and onto, and hence invertible.
Given , there is a such that . Hence,
To prove that is injective (i.e. one-to-one), consider an arbitrary and assume there are such that . Then, . Because of nullspace consistency, for a scalar . Since ,
Thus, . ∎
The inverse operator of is denoted as .
The inverse operator is a discrete integral operator, with zero initial data.
Continuing Example 3.4, the vector spaces and are
This can be seen as follows. For , using implies that the first component of is zero and that the second one can be chosen arbitrarily. For , note that both rows of are identical. Hence, every must have the same first and second component.
At the level of , the inverse of can be represented as
Indeed, if , then
Similarly, if , then
Hence, and , where is the identity on .
Now, we have introduced the inverse of , which is a discrete integral operator . However, the integral operator is only defined for elements of the space . Hence, one has to make sure that a generic right hand side vector is in the range of the derivative operator in order to apply the inverse . To guarantee this, components in the direction of grid oscillations must be removed. For this, the discrete projection/filter operator
will be used.
The projection/filter operator defined in (24) is an orthogonal projection onto the range of , i.e. onto . It is symmetric and positive semidefinite with respect to the scalar product induced by .
Clearly, is the usual orthogonal projection onto [roman2008advanced, Theorems 9.14 and 9.15]. Hence, is the orthogonal projection onto the orthogonal complement . In particular, for a (real or complex valued) discrete grid function ,
because of the Cauchy–Schwarz inequality [roman2008advanced, Theorem 9.3]. ∎
Now, all ingredients to mimic the integral equation (10) have been provided. Applying at first the discrete projection operator and second the discrete integral operator to a generic right hand side results in , which is a discrete analog of the integral . Additionally, the initial condition has to be imposed, which is done by adding the constant initial value as . Putting it all together, a new class of SBP schemes mimicking the integral equation (10) discretely is proposed as
Hence, the projection/filter operator (24) is
Thus, is a smoothing filter operator that removes the highest grid oscillations and maps a grid function into the image of the derivative operator . Hence, the inverse , the discrete integral operator, can be applied after , resulting in
Finally, for an arbitrary ,
3.1 Summarizing the Development
As stated earlier, the SBP time integration scheme (26) mimics the integral reformulation (10) of the ODE (1). Instead of using as discrete analog of the integral operator directly, the projection/filter operator defined in (24) must be applied first in order to guarantee that a generic vector is in the image of . Finally, the initial condition is imposed strongly.
The second summand vanishes because returns a vanishing value at , i.e. , since maps onto .
Moreover, taking the discrete derivative on both sides of (26) results in
since . The orthogonal projection operator on the right hand side ensures that this equation can be solved for for any given right hand side . Such a projection onto is necessary in the discrete case because of the finite dimensions.
In the context of SBP operators, two essentially different interpretations of integrals arise. Firstly, the integral gives the scalar product, approximated by the mass matrix which maps discrete functions to scalar values. Secondly, the integral is the inverse of the derivative with vanishing values at . This operator is discretized as on its domain of definition and maps a discrete grid function in to a discrete grid function in .
3.2 Linear Stability
In this section, linear stability properties of the new scheme (26) are established.
For nullspace consistent SBP operators, the scheme (26) is stable.
For the scalar linear ODE (8) with , the energy method will be applied to the scheme (26). Using to denote the complex conjugate, the difference of the energy at the final and initial time (see (31)) is
The second factor is non-negative because is positive semidefinite with respect to the scalar product induced by , cf. Lemma 3.9. Therefore, , implying that the scheme is stable. ∎
3.3 Characterization as a Runge–Kutta Method
Unsurprisingly, the new method (26) can be characterized as a Runge–Kutta method.
For nullspace consistent SBP operators that are at least first order accurate, the method (26) is a Runge–Kutta method with Butcher coefficients
First, note that one step (6) from zero to of a Runge–Kutta method with coefficients can be written as
where the right hand side vector is given by . Comparing this expression with
where the right hand side is given by , the form of and is immediately clear. The new value of the new SBP method (26) is
Using the SBP property (2) and ,
Inserting and from (32) results in
Because of , we have . Hence,
Comparing this expression with (36) yields the final assertion for . ∎
By definition, yields a vector with vanishing initial condition at the left endpoint, i.e. . Because of , we have , which is the first row of , where a notation as in Julia [bezanson2017julia] has been used. ∎
3.4 Operator Construction
To implement the SBP scheme (26), the product has to be computed, which is (except for a scaling by ) the matrix of the corresponding Runge–Kutta method, cf. Theorem 3.13. Since the projection operator maps vectors into , the columns of are in the image of the nullspace consistent SBP derivative operator . Hence, the matrix equation can be solved for , which is a matrix of the same size as , e.g. via a QR factorization, yielding the least norm solution. Then, we have to ensure that the columns of are in , since maps into . This can be achieved by subtracting from each column of , , where a notation as in Julia [bezanson2017julia] has been used. After this correction, we have .
3.5 Lobatto IIIA Schemes
General characterizations of the SBP-SAT scheme (4) on Radau and Lobatto nodes as classical collocation Runge–Kutta methods (Radau IA and IIA or Lobatto IIIC, respectively) have been obtained in [ranocha2019some]. A similar characterization will be obtained in this section.
If the SBP operator is given by the nodal polynomial collocation scheme on Lobatto Legendre nodes, the SBP method (26) is the classical Lobatto IIIA method.
The Lobatto IIIA methods are given by the nodes and weights of the Lobatto Legendre quadrature, just as the SBP method (26). Hence, it remains to prove that the classical condition is satisfied [butcher2016numerical, Section 344], where
In other words, all polynomials of degree must be integrated exactly by with vanishing initial value at . By construction of , see (35), this is satisfied for all polynomials of degree , since the grid oscillations are given by , where is the Legendre polynomial of degree , cf. [ranocha2020discrete, Example 3.6].
Finally, it suffices to check whether is integrated exactly by . The left hand side of (41) yields
since . On the right hand side, transforming the time domain to the standard interval , the analytical integrand of is
since the Legendre polynomials satisfy Legendre’s differential equation
Hence, vanishes exactly at the Legendre nodes for polynomials of degree , which are and the roots of . Thus, the analytical integral of vanishes at all grid nodes. ∎
Finally, the remaining Butcher coefficients are given by
which are exactly the coefficients of the Lobatto IIIA method with stages.
Since the Lobatto IIIA methods are neither nor stable, the new SBP method (26) is in general not or stable, too.
Because of Lemma 3.14, the classical Gauss, Radau IA, Lobatto IIIB, and Lobatto IIIC methods cannot be expressed in the form (26). The classical Radau I, Radau II, and Radau IIA methods are also not included in the class (26). For example, for two nodes, these methods have the matrices
while the methods (26) on the left and right Radau nodes yield the matrices
This also shows that the condition , which is common in the Runge–Kutta literature, is in general not satisfied by the schemes (26).
We develop a related SBP time integration scheme that includes the classical Lobatto IIIB collocation method in the appendix. Additionally, we mention why it seems to be difficult to describe Gauss collocation methods in a general SBP setting, cf. Appendix A.
4 Numerical Experiments
Numerical experiments corresponding to the ones in [nordstrom2013summation] will be conducted. The novel SBP methods (26) have been implemented in Julia v1.3 [bezanson2017julia] and Matplotlib [hunter2007matplotlib] has been used to generate the plots. The source code for all numerical examples is available online [ranocha2020classRepro].
4.1 Non-stiff Problem
The non-stiff test problem
with analytical solution is solved in the time interval using the SBP method (26) with the diagonal norm operators of [mattsson2004summation]. The errors of the numerical solutions at the final time are shown in Figure 1. As can be seen, they converge with an order of accuracy equal to the interior approximation order of the diagonal norm operators. For the operator with interior order eight, the error reaches machine precision for nodes and does not decrease further.
4.2 Stiff Problem
The stiff test problem
with analytical solution and parameter is solved in the time interval . The importance of such test problems for stiff equations has been established in [prothero1974stability]. Using the diagonal norm operators of [mattsson2004summation] for the method (26) yields the convergence behavior shown in Figure 2. As for the SBP-SAT schemes of [nordstrom2013summation], the order of convergence is reduced to the approximation order at the boundaries.
5 Summary and Discussion
A novel class of stable summation by parts time integration methods has been proposed. Instead of using simultaneous approximation terms to impose the initial condition weakly, the initial condition is imposed strongly.
Compared to the SAT approach, some linear and nonlinear stability properties such as and stability are lost in general. On the other hand, well-known stable methods such as the Lobatto IIIA schemes are included in this new SBP framework. Additionally, a related SBP time integration method has been proposed which includes the classical Lobatto IIIB schemes.
This article provides new insights into the relations of numerical methods and contributes to the discussion of whether SBP properties are necessarily involved in numerical schemes for differential equations which are provably stable.
Appendix A Lobatto IIIB Schemes
Another scheme similar to (26) can be constructed by considering as bijective operator acting on functions that vanish at the right endpoint. The corresponding matrix of the Runge–Kutta method becomes
The Lobatto IIIB methods are given by the nodes and weights of the Lobatto Legendre quadrature, just as the SBP method. Hence, it remains to prove that the classical condition is satisfied [butcher2016numerical, Section 344], where
In matrix vector notation, this can be written as
where the exponentiation is performed pointwise.
In other words, all polynomials of degree must be integrated exactly by with vanishing final value at . By construction of , this is satisfied for all polynomials of degree , and the proof can be continued as the one of Theorem 3.15. ∎
Up to now, it is unclear whether the classical collocation Runge–Kutta schemes on Gauss nodes can be constructed as special members of a family of schemes that can be formulated for (more) general SBP operators. The problem seems to be that SBP schemes rely on differentiation, while the conditions and describing the Runge–Kutta methods rely on integration. Thus, special compatibility conditions as in the proof of Theorem 3.15 are necessary. To the authors’ knowledge, no insights in this direction have been achieved for Gauss methods.
Research reported in this publication was supported by the King Abdullah University of Science and Technology (KAUST).