Resolving Entropy Growth from Iterative Methods

TL;DR

This paper investigates how iterative solvers affect entropy properties in discretized nonlinear conservation laws, revealing that relaxation techniques can preserve entropy dissipation and improve solution accuracy.

Contribution

It demonstrates that Newton's method can reverse entropy dissipation in implicit schemes and proposes relaxation as an effective remedy to maintain entropy properties.

Findings

Newton's method can cause anti-dissipative behavior in entropy schemes

Relaxation techniques effectively preserve entropy dissipation in iterative solvers

Entropy conservative schemes yield more accurate solutions even with larger tolerances

Abstract

We consider entropy conservative and dissipative discretizations of nonlinear conservation laws with implicit time discretizations and investigate the influence of iterative methods used to solve the arising nonlinear equations. We show that Newton's method can turn an entropy dissipative scheme into an anti-dissipative one, even when the iteration error is smaller than the time integration error. We explore several remedies, of which the most performant is a relaxation technique, originally designed to fix entropy errors in time integration methods. Thus, relaxation works well in consort with iterative solvers, provided that the iteration errors are on the order of the time integration method. To corroborate our findings, we consider Burgers' equation and nonlinear dispersive wave equations. We find that entropy conservation results in more accurate numerical solutions than non-conservative schemes, even when the tolerance is an order of magnitude larger.

Equations

$$u_t=\mathcal{L}(u),x\in (x_{\min },x_{\max }],t\in [t_0,t_e].$$

$$u_t+6u{u}_x=0,u(x,t=0)=u_0.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\frac{1}{2}\Vert u{\Vert }^2=0,$$

$${\mathsf{u}}_t+2(\mathsf{D}\mathrm{diag}(\mathsf{u})\mathsf{u}+\mathrm{diag}(\mathsf{u})\mathsf{D}\mathsf{u})=0.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\frac{1}{2}\Vert \mathsf{u}{\Vert }^2$$

$$=-2\Delta x({\mathsf{u}}^{\top }\mathsf{D}\mathrm{diag}(\mathsf{u})\mathsf{u}+{\mathsf{u}}^{\top }{\mathsf{D}}^{\top }\mathrm{diag}(\mathsf{u})\mathsf{u})$$

$$=-2\Delta x({\mathsf{u}}^{\top }(\mathsf{D}+{\mathsf{D}}^{\top })\mathrm{diag}(\mathsf{u})\mathsf{u}=0.$$

$$\mathsf{U}$$

$${\mathsf{u}}^{n+1}$$

$$\mathsf{F}(\mathsf{U})$$

$${\mathsf{u}}^{n+1}$$

$$\eta ({\mathsf{u}}^{n+1})=\frac{\Delta x}{2}(2\mathsf{U}-{\mathsf{u}}^n{)}^{\top }(2\mathsf{U}-{\mathsf{u}}^n)=\eta ({\mathsf{u}}^n)+2\Delta x({\mathsf{U}}^{\top }\mathsf{U}-{\mathsf{U}}^{\top }{\mathsf{u}}^n).$$

$${\mathsf{U}}^{\top }\mathsf{U}-{\mathsf{U}}^{\top }{\mathsf{u}}^n$$

$$=-\Delta t_n{\mathsf{U}}^{\top }(\mathsf{D}\mathrm{diag}(\mathsf{U})+\mathrm{diag}(\mathsf{U})\mathsf{D})\mathsf{U}$$

$$=-\Delta t_n{\mathsf{U}}^{\top }((\mathsf{D}+{\mathsf{D}}^{\top })\mathrm{diag}(\mathsf{U}))\mathsf{U}=0,$$

$${\mathsf{F}}^{\prime }({\mathsf{U}}^{(k)})\Delta \mathsf{U}+\mathsf{F}({\mathsf{U}}^{(k)})$$

$${\mathsf{U}}^{(k+1)}$$

$${\mathsf{F}}^{\prime }({\mathsf{U}}^{(k)})=\mathsf{I}+\Delta t_n\left(\mathrm{diag}({\mathsf{U}}^{(k)})\mathsf{D}+\mathrm{diag}(\mathsf{D}{\mathsf{U}}^{(k)})+2\mathsf{D}\mathrm{diag}({\mathsf{U}}^{(k)})\right).$$

$${\mathsf{U}}^{(k+1)}-{\mathsf{u}}^n$$

$$+\Delta t_n\underset{\mathsf{M}}{\underbrace{\left[\mathrm{diag}(\mathsf{D}{\mathsf{U}}^{(k)})+\mathsf{D}\mathrm{diag}({\mathsf{U}}^{(k)})\right]}}\Delta \mathsf{U}=0.$$

$$\eta ({\mathsf{u}}^{n+1})=\eta ({\mathsf{u}}^n)+2\Delta x\left(({\mathsf{U}}^{(k+1)}{)}^{\top }{\mathsf{U}}^{(k+1)}-({\mathsf{U}}^{(k+1)}{)}^{\top }{\mathsf{u}}^n\right).$$

$$({\mathsf{U}}^{(k+1)}{)}^{\top }{\mathsf{U}}^{(k+1)}-({\mathsf{U}}^{(k+1)}{)}^{\top }{\mathsf{u}}^n=-\Delta t_n({\mathsf{U}}^{(k+1)}{)}^{\top }\mathsf{M}\Delta \mathsf{U},$$

$$\eta ({\mathsf{u}}^{n+1})=\eta ({\mathsf{u}}^n)-2\Delta x\Delta t_n({\mathsf{U}}^{(k+1)}{)}^{\top }\mathsf{M}\Delta \mathsf{U}.$$

$${\mathsf{M}}^{\top }{\mathsf{U}}^{(k)}$$

$$=\mathrm{diag}(\mathsf{D}{\mathsf{U}}^{(k)})\mathrm{diag}({\mathsf{U}}^{(k)})1-\mathrm{diag}({\mathsf{U}}^{(k)})\mathrm{diag}(\mathsf{D}{\mathsf{U}}^{(k)})1=0.$$

$$\eta ({\mathsf{u}}^{n+1})=\eta ({\mathsf{u}}^n)-2\Delta x\Delta t_n\Delta {\mathsf{U}}^{\top }\mathsf{M}\Delta \mathsf{U}.$$

$$\Delta {\mathsf{U}}^{\top }\mathsf{M}\Delta \mathsf{U}=\Delta {\mathsf{U}}^{\top }\left[\mathrm{diag}(\mathsf{D}{\mathsf{U}}^{(k)})+\mathsf{D}\mathrm{diag}({\mathsf{U}}^{(k)})\right]\Delta \mathsf{U}.$$

$$u_0=\frac{c}{2}{\mathrm{sech}}^2\left(\frac{\sqrt{c}}{2}x\right),$$

$${\tilde{\mathsf{F}}}^{\prime }({\mathsf{U}}^{(k)}):=\mathsf{I}+\Delta t_n\left(\mathrm{diag}({\mathsf{U}}^{(k)})\mathsf{D}+\mathsf{D}\mathrm{diag}({\mathsf{U}}^{(k)})\right)$$

$${\mathsf{F}}^{\prime }({\mathsf{U}}^{(k)})(\tilde{\mathsf{U}}-{\mathsf{U}}^{(k)})+\mathsf{F}({\mathsf{U}}^{(k)})$$

$${\mathsf{U}}^{(k+1)}$$

$${\mathsf{F}}^{\prime }({\mathsf{U}}^{(k)})\Delta \mathsf{U}+\mathsf{F}({\mathsf{U}}^{(k)})$$

$${\mathsf{U}}^{(k+1)}$$

$$({\mathsf{U}}^{(k+1)}{)}^{\top }{\mathsf{U}}^{(k+1)}-({\mathsf{U}}^{(k+1)}{)}^{\top }{\mathsf{u}}^n=0,$$

$${\alpha }_k=-\frac{⟨\Delta \mathsf{U},{\mathsf{U}}^{(k)}⟩+⟨\tilde{\mathsf{U}},{\mathsf{U}}^{(k)}-{\mathsf{u}}^n⟩}{\Vert \Delta \mathsf{U}{\Vert }^2}.$$

$${\mathsf{u}}_{\gamma }^{n+1}={\mathsf{u}}^n+\gamma ({\mathsf{u}}^{n+1}-{\mathsf{u}}^n),$$

Full text

${}^{\mbox{\tiny\rm 1}}$*Centre for mathematical sciences, Lund University, Lund, Sweden

${}^{\mbox{\tiny\rm 2}}$Applied Mathematics, University of Hamburg, Hamburg, Germany

email:   viktor.linders@math.lu.se

mail@ranocha.de

    philipp.birken@na.lu.se *

Keywords: iterative methods, entropy conservation, implicit methods, dispersive wave equations

1 Introduction

For many partial differential equations (PDEs) in computational fluid dynamics (CFD), the notion of mathematical entropy plays an important role. It can be that entropy is preserved or that the solution satisfies an entropy inequality. Alternatively, for nonlinear systems such as the Euler equations, an entropy inequality can be used to select a unique weak solution, at least in 1D. Enforcing such an entropy inequality on the discrete level has, for certain equations, led to schemes that are provably convergent to weak entropy solutions [25].

In recent years it has been observed that the robustness of high-order discontinuous Galerkin (DG) discretizations for compressible turbulent flows, greatly benefits from discrete entropy stability; see [18] and the references therein for an overview. Following the groundbreaking paper [16], a space-time version of the DG spectral element method (DG-SEM) with a specific choice of fluxes, was proven to be entropy-stable for systems of hyperbolic conservation laws [17], giving the first fully discrete high order scheme with this property. This method is inherently implicit.

Proofs of entropy stability for schemes involving implicit time integration assume that the arising systems of nonlinear equations are solved exactly. In practice, iterative solvers are used, which are terminated when a tolerance is reached. These in turn, have rarely been constructed with the intention of preserving physical invariants. Indeed, the impact of iterative solvers on the conservation of linear invariants was analyzed in [6, 29], and it was found that many iterative methods violate global or local conservation. A recent study of Jackaman and MacLachlan [20] focuses on Krylov subspace methods for linear problems with quadratic invariants, and is related in spirit to this paper. Given the importance of entropy estimates in discrete schemes, there is a need to analyze the role played by iterative solvers in this context.

We consider PDEs posed in a single spatial dimension with periodic boundary conditions:

[TABLE]

Here, ${\cal L}$ is a nonlinear differential operator. The cases considered are such that there is a convex nonlinear functional $\eta(u)$, here called the entropy, which is conserved. By splitting the spatial terms in certain ways, entropy conservative semi-discretizations are obtained with the aid of skew-symmetric difference operators. Next, we apply well-known implicit time integration methods that are able to either conserve or dissipate entropy, resulting in nonlinear algebraic systems with entropy bounded solutions.

Using convergent iterative solvers within implicit time integration, the iteration error (and thus the entropy error) can be made arbitrarily small by iterating long enough. However, such an approach is ill-advised since the number of iterations dictates the efficiency of the scheme. On the other hand, it is desirable to iterate long enough to prevent the iteration error from affecting the time integration error. The iteration error should thus be kept on the order of the local time integration error, but smaller in magnitude. The interested reader will find an overview of the use of iterative solvers in computational fluid dynamics in [5].

In this setting, we analyze the entropy behavior of Newton’s method, which forms the basis of many iterative solvers. As it turns out, Newton’s method can generate undesired growth (or decay) of the entropy error. We consider Burgers’ equation in Section 2 and demonstrate that Newton’s method can turn an entropy dissipative scheme into an anti-dissipative method when the tolerance is comparable to the error in the time step. A detailed analysis reveals that the source of the entropy error is the Jacobian of the spatial discretization. In Section 3, strategies are evaluated that recover the entropy conservation after each Newton iteration. It is seen that these strategies come with a large cost to the efficiency of the solver, since they significantly reduce the convergence rate.

To get the best of both worlds — entropy conservation and fast convergence — we introduce relaxation methods in Section 4. Under mild assumptions, relaxation can be used to recover any convex entropy. We apply relaxation to nonlinear dispersive wave equations in Sections 5 and 6. Experiments reveal that entropy conservation leads to numerical solutions of considerably higher quality than non-conservative schemes, even when larger tolerances on the iterations are used. We thus conclude that relaxation works well in combination with iterative methods. Finally, we summarize the developments and discuss our conclusions in Section 7.

1.1 Software

The numerical experiments discussed in this article are implemented in Julia [4]. We use the Julia packages SummationByPartsOperators.jl [40], ForwardDiff.jl [46], and Krylov.jl [35]. All source code required to reproduce the numerical experiments is available in our repository [32].

2 Burgers’ equation

As a starting point, we consider Burgers’ equation as a classical model for nonlinear conservation laws. It is well known how to design both entropy conservative and entropy dissipative discretizations for this problem. The purpose of this section is to demonstrate that Newton’s method can destroy the entropic behavior of the discretization and cause undesired entropy growth (or decay). This example serves as an illustration of a general truth: Iterative methods may destroy the design principles upon which a discretization has been built.

2.1 The continuous case

Consider the inviscid Burgers’ equation defined on a 1D periodic domain,

[TABLE]

The factor 6 has been added for consistency with the Korteweg-de Vries equation discussed in Section 5. Multiplying (2) by the solution $u$ and integrating in space results in the identity

[TABLE]

where periodicity has been used to eliminate the boundary terms. Here, $\|\cdot\|$ denotes the $L^{2}$ norm. Evidently, the quantity $\eta(u)=\frac{1}{2}\|u\|^{2}$ is conserved. Throughout, we refer to $\eta$ as an entropy for Burgers’ equation (2).

2.2 The semi-discrete case

It is well known how to design discretizations of (2) that mimic the entropy conservation of the continuous problem. Let $\mathbf{\boldsymbol{\mathsf{D}}}$ be a skew-symmetric matrix that approximates the spatial derivative operator. In this section, we use classical fourth-order accurate central finite differences with periodic boundaries. Consider the semi-discretization

[TABLE]

Here and elsewhere a sans serif font is used to denote quantities evaluated on a discrete, uniform computational grid with grid spacing $\Delta x$. The formulation (3) constitutes a so-called skew-symmetric split form [47, eq. (6.40)]. It arises from the identity 6uu_{x}=2\bigl{(}(u^{2})_{x}+uu_{x}\bigr{)}.

Multiplying (3) by $\Delta x\mathbf{\boldsymbol{\mathsf{u}}}^{\top}$ and using the skew-symmetry of $\mathbf{\boldsymbol{\mathsf{D}}}$ yields

[TABLE]

The final equality follows from the skew-symmetry of $\mathbf{\boldsymbol{\mathsf{D}}}$. Here we have defined the discrete norm $\|\mathbf{\boldsymbol{\mathsf{u}}}\|^{2}\equiv\Delta x\mathbf{\boldsymbol{\mathsf{u}}}^{\top}\mathbf{\boldsymbol{\mathsf{u}}}$ with a slight notational abuse. Evidently the semi-discrete entropy $\eta(\mathbf{\boldsymbol{\mathsf{u}}})=\frac{1}{2}\|\mathbf{\boldsymbol{\mathsf{u}}}\|^{2}$ is conserved.

2.3 The fully discrete case

Since the entropy $\eta$ is a quadratic functional of the solution, a fully discrete scheme that conserves entropy is obtained by discretizing (3) in time using the implicit midpoint rule. Further, any $B$-stable Runge-Kutta method will be entropy dissipative [8, Section 357]. Here, we consider both the (conservative) implicit midpoint rule and the (dissipative) fourth-order, three-stage Lobatto IIIC method. The entropy analyses in this section and the next are for simplicity reserved for the midpoint rule. The analogous results for Lobatto IIIC are found in Appendix A. They heavily utilize the fact that Lobatto IIIC is equivalent to a Summation-By-Parts (SBP) method in time [39]; see [50, 36, 33, 7, 31] for further developments of this topic.

For a generic system of ordinary differential equations $\mathbf{\boldsymbol{\mathsf{u}}}_{t}=\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{u}}})$, the midpoint rule can be expressed as a Runge-Kutta method as follows:

[TABLE]

Here, $\mathbf{\boldsymbol{\mathsf{U}}}$ denotes an intermediate stage used to define the numerical solution $\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}\approx\mathbf{\boldsymbol{\mathsf{u}}}(t_{n+1})$ in terms of the solution at the previous time step; $\mathbf{\boldsymbol{\mathsf{u}}}^{n}\approx\mathbf{\boldsymbol{\mathsf{u}}}(t_{n})$. It is understood that $\Delta t_{n}=t_{n+1}-t_{n}$.

Applied to the semi-discretization (3), the fully discrete scheme becomes, after slight rearrangement,

[TABLE]

From the second line in (4), the entropy $\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n+1})$ is given by

[TABLE]

Left multiplication of the first line in (4) by $\mathbf{\boldsymbol{\mathsf{U}}}^{\top}$ reveals that

[TABLE]

where the skew-symmetry of $\mathbf{\boldsymbol{\mathsf{D}}}$ has been used in the same way as in the semi-discrete analysis. Consequently, $\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n+1})=\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n})$, hence entropy is conserved.

2.4 Newton’s method

The analysis above assumes that the stage equations in the first line of (4) are solved exactly. In practice this is infeasible. Instead, the stage vector $\mathbf{\boldsymbol{\mathsf{U}}}$ is approximated using iterative methods. Here we consider Newton’s method as an illustrative example. If Newton iterates are computed until the residual is sufficiently small, then entropy conservation is effectively retained. However, for efficiency reason it is desirable to terminate the iterates when the residual is smaller than some tolerance, chosen to reflect the overall expected accuracy of the scheme. We will show that Newton’s method is not entropy conservative from one iteration to the next. Consequently, schemes utilizing Newton’s method (and other iterative methods) tend to lose the design principle upon which the discrete scheme is built, namely entropy conservation.

The iterates produced by Newton’s method applied to the stage equation in (4) are obtained as

[TABLE]

Here, $\mathbf{\boldsymbol{\mathsf{F}}}^{\prime}$ denotes the Jacobian of $\mathbf{\boldsymbol{\mathsf{F}}}$ and can be explicitly evaluated; see [11]:

[TABLE]

Here, $\mathbf{\boldsymbol{\mathsf{I}}}$ is the identity matrix. Inserting this expression into (7), explicitly writing out $\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})$ from (4), and utilizing the fact that $\Delta\mathbf{\boldsymbol{\mathsf{U}}}=\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}-\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}$ leads to the following equation for $\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}$:

[TABLE]

This equation replaces the stage equations in (4) when $k+1$ Newton iterations are performed.

Suppose that Newton’s method is terminated after $k+1$ iterations and that the solution is updated as $\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}=2\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}$. As in (5), the entropy is given by

[TABLE]

Upon computing the parenthesized term, note that the first line of (9) is precisely $\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\Delta t_{n}\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})$. By the same derivation as in (6) it follows that $(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})^{\top}\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})=0$. Consequently,

[TABLE]

and hence

[TABLE]

A slightly more elegant expression is obtained by replacing $(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})^{\top}$ with $\Delta\mathbf{\boldsymbol{\mathsf{U}}}^{\top}$. This can be done since the vector $\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}$ lies in the kernel of $\mathbf{\boldsymbol{\mathsf{M}}}^{\top}$. To see this, note that

[TABLE]

Here, the skew-symmetry of $\mathbf{\boldsymbol{\mathsf{D}}}$ has been used together with the fact that any vector $\mathbf{\boldsymbol{\mathsf{v}}}$ satisfies $\mathbf{\boldsymbol{\mathsf{v}}}=\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{v}}})\mathbf{\boldsymbol{\mathsf{1}}}$, where $\mathbf{\boldsymbol{\mathsf{1}}}$ is the vector of ones. The final equality follows by the commutativity of diagonal matrices.

We summarize these observations in the following:

Proposition 1**.**

Consider the discretization (4) of Burgers’ equation (2) and apply $(k+1)$ iterations with Newton’s method to the stage equations. The entropy of the resulting numerical solution satisfies

[TABLE]

Thus, the entropy error $\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n+1})-\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n})$ depends on the indefinite quadratic form

[TABLE]

Consequently, Newton’s method may cause both entropy growth and decay.

With Lobatto IIIC in place of the midpoint rule, the entropy error has an identical form, although $\Delta\mathbf{\boldsymbol{\mathsf{U}}}$ and $\mathbf{\boldsymbol{\mathsf{M}}}$ incorporate all intermediate stages in this case. Assuming that the iterates converge, $\Delta\mathbf{\boldsymbol{\mathsf{U}}}$ will decrease until the entropy error is vanishingly small. However, in practical applications, we terminate the iterations at some tolerance matching the accuracy of the scheme. In such circumstances, the entropy error must be corrected by other means.

As an illustrative example we compute a single time step for Burgers’ equation (2). The spatial domain is set to $(-10,10]$ and the initial condition is given by

[TABLE]

where $c=2$. In the spatial discretization we set $\Delta x=0.1$. Both Lobatto IIIC and the midpoint rule are used in time with $\Delta t=0.5$. The entropy $\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n})$ (with $n=1$) and residual $\|\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\|$ are shown in Fig. 1 when $k$ Newton iterations have been used to approximate the solution to the stage equation in (4). Throughout, $\mathbf{\boldsymbol{\mathsf{U}}}^{(0)}=\mathbf{\boldsymbol{\mathsf{u}}}^{n}$ is used as initial guess. This choice is known to conserve linear invariants [6, 29].

The entropy visibly grows when too few Newton iterations are used. As expected, it drops back to its correct value with further iterations, testifying to the indefiniteness of the quadratic form in (10). The residual converges quadratically as expected for Newton’s method.

These examples show that if Newton’s method is applied with a tolerance around $10^{-3}$, then entropy growth is to be expected in these simulations. Note that entropy growth may occur even for the provably entropy dissipative discretization.

3 Strategies for recovering entropy conservation

There are several ways in which we can recover entropy conservation within Newton’s method. The most obvious is to simply keep iterating until the entropy error is negligible. However, this process might be costly since we may be forced to iterate to residuals that are smaller than what is motivated by the accuracy of the discretization.

The entropy error in (10) opens up several avenues for modifying Newton’s method to recover entropy conservation in each iteration. In the following subsections we describe two strategies that are independent of the choice of tolerance. However, as we will see, they both converge slower than Newton’s method. The two methods are, respectively

• •

Method of Newton-type: By modifying the Jacobian, the entropy error is removed in its entirety.

• •

Inexact Newton: Entropy conservation is recovered using a line search.

Throughout this section, we confine our attention to the entropy conservative discretization that utilizes the implicit midpoint rule.

3.1 Method of Newton-type

The entropy error (10) caused by Newton’s method originates from the Jacobian of the spatial discretization. A simple remedy is thus to modify the Jacobian by removing the problematic terms, thereby obtaining a method of Newton-type. In other words, replacing the exact Jacobian $\mathbf{\boldsymbol{\mathsf{F}}}^{\prime}$ in (8) with the approximation

[TABLE]

should result in a vanishing entropy error. However, this will simultaneously reduce the convergence speed from quadratic to linear [22, Chapter 5].

Repeating the experiment in section 2.4 reveals that this is indeed the case. Fig. 2 shows the entropy and residual. Clearly the entropy is conserved as desired. However, the convergence is no longer quadratic but linear. After 14 iterations, the residual is comparable to four iterations with Newton’s method. It is thus questionable if anything has been gained.

3.2 Inexact Newton

As a second alternative, entropy conservation can be recovered through a line search. We replace Newton’s method as stated in (7) by

[TABLE]

where $\alpha_{k}\in[0,1]$ is a sequence of scalar parameters. This formulation is equivalent to the inexact Newton’s method

[TABLE]

If the sequence $\{\alpha_{k}\}$ is such that $(1-\alpha_{k})=\mathcal{O}(\|\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\|)$, then quadratic convergence is retained under standard assumptions. Under the much less severe condition $0\leq 1-\alpha_{k}<1$, covergence is generally linear [14].

The essential property used in the proof that the fully discrete scheme (4) is entropy conservative is that the stage vector $\mathbf{\boldsymbol{\mathsf{U}}}$ satisfies $\mathbf{\boldsymbol{\mathsf{U}}}^{\top}\mathbf{\boldsymbol{\mathsf{U}}}-\mathbf{\boldsymbol{\mathsf{U}}}^{\top}\mathbf{\boldsymbol{\mathsf{u}}}^{n}=0$. Proposition 1 shows that the same relation does not hold for the Newton iterations. To ensure entropy conservation for the inexact Newton method, $\alpha_{k}$ must therefore be chosen such that

[TABLE]

where $\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}$ is given as in (12). This is a quadratic equation in $\alpha_{k}$. Assuming that $\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}$ has correct entropy, the two roots are $\alpha_{k}=0$ and

[TABLE]

Here, $\langle\cdot,\cdot\rangle$ denotes the Euclidean inner product scaled by $\Delta x$. The trivial root $\alpha_{k}=0$ results in $\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}=\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}$ and is of no use, hence only the second root needs to be considered.

Repeating the previous experiment with inexact Newton leads to the results displayed in Fig. 3. Entropy is conserved as expected. The convergence is faster than it was with the method of Newton-type, however it is still linear. This happens since $\alpha_{k}\in[0,1]$ for each $k$ but does not approach unity, as is necessary for quadratic convergence. Instead, its value appears to plateau around $0.95$ in this particular experiment.

4 Theory of relaxation methods

The idea of a line search discussed in the previous section can be applied, not necessarily after each Newton iteration, but alternatively just once after a full time step computed using multiple Newton iterations. The resulting schemes are known as relaxation methods [23, 45, 43]. As we will see, relaxation solves the problem with entropy growth without impacting the convergence of Newton’s method, or that of any other iterative method.

While a general theory of relaxation methods that includes multistep methods is available [43], we restrict the following discussion to one-step methods for simplicity. Thus, we consider a system of ODEs \mathbf{\boldsymbol{\mathsf{u}}}^{\prime}(t)=f\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}(t)\bigr{)} and a one-step method. We again assume that there is a (nonlinear and sufficiently smooth) functional $\eta$ of interest, which we call an entropy.

In this setting, the basic idea of relaxation methods is to perform a line search along the secant line connecting the new approxiation $\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}$ produced by a given time integration method, and the previous value $\mathbf{\boldsymbol{\mathsf{u}}}^{n}$. This results in the relaxed value

[TABLE]

where $\gamma\approx 1$ is the relaxation parameter chosen to enforce the desired conservation or dissipation of $\eta$. This idea goes back to Sanz-Serna [48, 49] and Dekker & Verwer [13, pp. 265–266], who considered entropies $\eta$ given by inner product norms. However, the first approaches resulted in an order reduction [9], which has been fixed in [23] for inner product norms and extended to general entropies in [45, 43]. Some further developments can be found in [42, 41, 2, 21, 26, 27].

To motivate the approach, we augment the ODE by the equations $t^{\prime}=1$ and the entropy evolution given by the chain rule, leading to

[TABLE]

Given a suitable estimate of the entropy \eta^{\mathrm{new}}=\eta\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}(t^{n+1})\bigr{)}+\mathcal{O}(\Delta t^{p+1}), where $p$ is the order of the time integration method, relaxation methods enforce

[TABLE]

by inserting the second equation into the third one, resulting in a scalar equation for the relaxation parameter $\gamma$. Next, the numerical solution is updated according to the second equation and the current simulation time is set to $t^{n+1}_{\gamma}$. The last step is required to avoid order reduction. Finally, the time integration is continued using $(t^{n+1}_{\gamma},\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}_{\gamma})$ instead of $(t^{n+1},\mathbf{\boldsymbol{\mathsf{u}}}^{n+1})$.

If entropy conservation is to be enforced, the canonical choice of the entropy estimate is \eta^{\mathrm{new}}=\eta\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}^{n}\bigr{)}. In case of entropy dissipation, $\eta^{\mathrm{new}}$ can be estimated [23, 45, 43]: For an RK method

[TABLE]

with non-negative weights $b_{i}\geq 0$, choose

[TABLE]

This leads to an entropy dissipative scheme.

The general theory of relaxation methods yields the following result [45, 43]:

Theorem 1**.**

Consider the system of ODEs \mathbf{\boldsymbol{\mathsf{u}}}^{\prime}(t)=f\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}(t)\bigr{)}. Assume $\mathbf{\boldsymbol{\mathsf{u}}}^{n}=\mathbf{\boldsymbol{\mathsf{u}}}(t^{n})$ and $\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}=\mathbf{\boldsymbol{\mathsf{u}}}(t^{n+1})+\mathcal{O}(\Delta t^{p+1})$ with $p\geq 2$ and $t^{n+1}=t^{n}+\Delta t$. If \eta^{\mathrm{new}}=\eta\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}(t^{n+1})\bigr{)}+\mathcal{O}(\Delta t^{p+1}), $\Delta t>0$ is small enough, and the non-degeneracy condition

[TABLE]

is satisfied with $c\neq 0$, then there is a unique $\gamma=1+\mathcal{O}(\Delta t^{p-1})$ that satisfies the relaxation condition (13) and the resulting relaxation method is of order $p$ such that

[TABLE]

The crucial observation for our application in this article is that there is no further assumption on the preliminary time step update $\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}$ other than the accuracy constraint $\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}=\mathbf{\boldsymbol{\mathsf{u}}}(t^{n+1})+\mathcal{O}(\Delta t^{p+1})$ in Theorem 1. In particular, the results can be applied to implicit Runge-Kutta methods where the stage equations are solved inexactly with Newton’s method (or another iterative method) up to some tolerance as long as the tolerance of the nonlinear solver does not affect the accuracy of the time integration method. A perturbation analysis shows that an approximation $\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}+\mathbf{\boldsymbol{\mathsf{\varepsilon}}}$ results in a perturbed relaxation parameter that can be estimated. In general, we need that the size of the perturbation $\mathbf{\boldsymbol{\mathsf{\varepsilon}}}$ is at most of the same order as the local error of the baseline time integration method, i.e. $\mathcal{O}(\Delta t^{p+1})$. This is precisely the situation desired for iterative methods in practice.

5 Korteweg-de Vries equation

We have motivated our study by looking at the influence of inexact solutions obtained by Newton’s method applied to provably entropy conservative and entropy dissipative time discretization methods for Burgers’ equation. Next, we will look at a more complicated equation where i) implicit methods are required due to stiffness constraints and ii) entropy conservative methods have a significant advantage. Indeed, the classical Korteweg-de Vries (KdV) equation

[TABLE]

is stiff due to the linear dispersive term $u_{xxx}$. We consider solutions of the form

[TABLE]

In the context of the KdV equation, (15) describes a soliton propagating with speed $c$. The behavior of numerical time integration methods applied to soliton solutions has been analyzed in [12]: If the time integration method conserves the linear invariant $\int u$ and the quadratic invariant $\int u^{2}$, the error of the numerical solution has a leading-order term that grows linearly with time. Otherwise, the error grows quadratically in time at leading order.

De Frutos and Sanz-Serna [12] verified this numerically using the implicit midpoint rule and an entropy non-conservative third order SDIRK method. They used essentially exact solutions of the stage equations. Under such conditions, relaxation has been demonstrated to yield the same reduced error growth when applied to the SDIRK method [42].

Here, we extend the investigations and focus on the influence of non-negligible tolerances of the nonlinear solver. In other words, we approximate the stage solutions to the stage equations with a tolerance that matches that of the time discretization.

5.1 Numerical methods

We consider the soliton solution (15) of the KdV equation (14) in the periodic domain $x\in(-10,10]$. The spatial discretization is the same as for Burgers’ equation, i.e., five-point centered (i.e. skew-symmetric) differences with a spatial increment $\Delta x=0.1$. In time, the implicit midpoint rule is used with time step $\Delta t=0.05$. The final time is set to $t=1000$.

The solution to the nonlinear system arising within each time step is approximated with Newton-GMRES. Here, the absolute tolerance is set to zero and the relative tolerance is chosen as $tol\in\{10^{-3},10^{-4},10^{-5}\}$.

Applying relaxation to preserve a quadratic invariant yields a linear equation for the relaxation parameter that we solve analytically. For the implicit midpoint rule, it is given by

[TABLE]

The tolerances for GMRES are chosen using the procedure suggested by Eisenstat and Walker [14], which is designed to recover quadratic convergence without over-iterating the linear solver in the early Newton iterations. This procedure requires the choice of two scalar parameters, $\eta_{\mathrm{max}}$ and $\gamma$ (not to be confused with the relaxation parameter). Here, we follow the implementation described in [22, Chapter 6.3] with the parameter choices $(\gamma,\eta_{\mathrm{max}})=(0.9,0.9)$.

5.2 Numerical results

Fig. 4 shows the discrete $L^{2}$ error and the error in the quadratic invariant $\eta(\mathbf{\boldsymbol{\mathsf{u}}})=\frac{1}{2}\|\mathbf{\boldsymbol{\mathsf{u}}}\|^{2}$ for four cases: The unrelaxed solver with relative tolerance chosen from the set $\{10^{-3},10^{-4},10^{-5}\}$ and the relaxed solver with relative tolerance $10^{-3}$. The absolute tolerance is set to zero.

The discrete $L^{2}$ error after the first time step is roughly $10^{-3}$. Newton’s method with a relative tolerance of $10^{-3}$ results in a superlinear error growth in time. The error eventually plateaus due to the fact that the numerical solution drifts completely out of phase with the true solution. It subsequently drifts in and out of phase, which gives rise to the oscillating nature of the error. Subsequent growth is also a consequence of the soliton losing its initial shape.

Shrinking the tolerance to $10^{-4}$ reduces the error growth rate significantly. However, the growth is still superlinear and reaches the same plateau seen with the larger tolerance. Yet another tolerance reduction to $10^{-5}$ leads to considerably improved results. The error grows linearly throughout the simulation time. However, the quadratic invariant is still not conserved.

By instead applying relaxation after every time step and using the original tolerance $10^{-3}$, the error growth rate is reduced even further. Additionally, the entropy is conserved up to machine precision. Thus, entropy conservation yields a better numerical solution even with a two orders of magnitude larger tolerance compared to the non-conservative scheme.

Figure 5 shows the results using the same space discretization but the fourth-order, three-stage Lobatto IIIC method in time with $\Delta t=0.1$. Since this method is $B$-stable, it dissipates the quadratic entropy when the stage equations are solved exactly. Here, we again apply Newton-GMRES with different tolerances. In this case, we set the absolute and relative tolerances to the same value, again chosen from the set $\{10^{-3},10^{-4},10^{-5}\}$. The $L^{2}$ error after one time step is approximately $10^{-3}$.

The resulting scheme is anti-dissipative with tolerances $10^{-3}$ and $10^{-4}$. The error displays similar characteristics to that of the midpoint rule with an initially superlinear growth.

With the tolerance $10^{-5}$ the entropy error is negative. The $L^{2}$ error shows a hump, suggesting that the numerical solution first drifts out of phase in one direction, then turns around and drifts the other way. Eventually the error starts to grow superlinearly as seen for the other tolerances.

As for the midpoint rule, relaxation leads to entropy conservation, a slower error growth and an overall better numerical solution. Again, this holds even when the tolerance is two orders of magnitude larger than the non-conservative scheme.

6 Benjamin-Bona-Mahony equation

An alternative to the stiff dispersive term of the KdV equation has been proposed by Benjamin, Bona, and Mahony (BBM) in [3], leading to the BBM equation

[TABLE]

Due to the mixed-derivative term $-u_{txx}$, the system is not stiff, but a linear elliptic equation must be solved to evaluate the time derivative $u_{t}$. Assuming periodic boundary conditions, the functionals

[TABLE]

are invariants of solutions to the BBM equation [37]. The error growth under time discretization has been analyzed in [1]: For methods conserving the linear invariant and one of the nonlinear invariants, the error of solitary waves grows linearly in time while general methods have a quadratically growing error (both at leading order). We use this example to investigate the behavior of the methods also for non-quadratic entropies.

6.1 Numerical methods

We introduce numerical methods conserving important invariants of the BBM equation (16) with periodic boundary conditions. To broaden the scope of the following derivations, we assume that a discrete inner product is given by a diagonal, symmetric, and positive-definite matrix $\mathbf{\boldsymbol{\mathsf{M}}}$. For the classical finite-difference methods described above, $\mathbf{\boldsymbol{\mathsf{M}}}=\Delta x\mathbf{\boldsymbol{\mathsf{I}}}$. Next, we need derivative operators $\mathbf{\boldsymbol{\mathsf{D}}}_{1,2}$ approximating the first and second derivative operators. We require compatibility with integration by parts, i.e., we need that $\mathbf{\boldsymbol{\mathsf{D}}}_{1}$ is skew-symmetric with respect to $\mathbf{\boldsymbol{\mathsf{M}}}$ and that $\mathbf{\boldsymbol{\mathsf{D}}}_{2}$ is symmetric and negative semidefinite with respect to $\mathbf{\boldsymbol{\mathsf{M}}}$. This is satisfied by classical central finite differences and Fourier collocation methods but also by appropriate continuous and discontinuous Galerkin methods; see e.g. [44].

For brevity, let $\mathbf{\boldsymbol{\mathsf{u}}}^{2}=\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{u}}})\mathbf{\boldsymbol{\mathsf{u}}}$ denote the elementwise square of $\mathbf{\boldsymbol{\mathsf{u}}}$. The convention can similarly be extended to other elementwise powers as necessary.

Using such periodic derivative operators, the semidiscretization

[TABLE]

conserves both the linear and the quadratic invariant [44], which are discretely represented as

[TABLE]

All symplectic Runge-Kutta methods, such as the implicit midpoint rule, conserve these linear and quadratic invariants [19].

Next, we construct a method conserving the cubic invariant.

Theorem 2**.**

The semidiscretization

[TABLE]

conserves the linear invariant $J^{\text{BBM}}_{1}(\mathbf{\boldsymbol{\mathsf{u}}})$ and the cubic invariant

[TABLE]

for commuting periodic derivative operators $\mathbf{\boldsymbol{\mathsf{D}}}_{1}$, $\mathbf{\boldsymbol{\mathsf{D}}}_{2}$.

Proof.

Conservation of the linear invariant follows from Lemma 2.1 and Lemma 2.2 of [44] since

[TABLE]

Given conservation of the linear invariant, conservation of the cubic invariant is equivalent to conservation of the Hamiltonian

[TABLE]

This Hamiltonian is conserved, since $\mathbf{\boldsymbol{\mathsf{M}}}(\mathbf{\boldsymbol{\mathsf{I}}}-\mathbf{\boldsymbol{\mathsf{D}}}_{2})^{-1}\mathbf{\boldsymbol{\mathsf{D}}}_{1}$ is skew-symmetric [44, Lemma 2.3] and

[TABLE]

∎

The average vector field (AVF) method [34] for an ODE \mathbf{\boldsymbol{\mathsf{u}}}^{\prime}(t)=f\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}(t)\bigr{)} is given by

[TABLE]

It conserves the Hamiltonian $\mathcal{H}$ of Hamiltonian systems \mathbf{\boldsymbol{\mathsf{u}}}^{\prime}(t)=S\,\mathcal{H}^{\prime}\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}(t)\bigr{)} with a constant (in time) skew-symmetric operator $S$ [38]. For quadratic Hamiltonians, it is equivalent to the implicit midpoint rule; for (up to) quartic Hamiltonians, it is equivalent to

[TABLE]

see [10]. Thus, we obtain

Proposition 2**.**

The time integration method (20) applied to the semidiscretization (19) of the BBM equation with periodic boundary conditions conserves the linear and cubic invariants under the conditions given in Theorem 2.

6.2 Numerical results

We consider the traveling wave solution

[TABLE]

with speed $c=1.2$ in the periodic domain $(-90,90]$. We apply the semidiscretizations described above with Fourier collocation methods [24, Chapter 4] using $2^{6}$ nodes. The time integration methods use fixed time steps $\Delta t=0.25$. The nonlinear systems are solved with the same Newton-GMRES method as for the KdV equation. The relative tolerances are, as before, chosen as $tol\in\{10^{-3},10^{-4}\}$.

Applying relaxation to preserve a quadratic invariant yields a linear equation for the relaxation parameter that we solve analytically. For a cubic invariant, the relaxation parameter is determined by a quadratic equation that we solve analytically; we always choose the solution closer to unity.

The results of these numerical experiments are shown in Fig. 6. The discrete $L^{2}$ error after the first time step is of the order $10^{-3}$. Using a relative tolerance of $10^{-3}$ for Newton-GMRES results in a discrete $L^{2}$ error that grows superlinearly in time. Reducing the relative tolerance of Newton’s method to $10^{-4}$ reduces the error growth rate to linear, resulting in significantly better results for long-time simulations. Applying relaxation after each time step with a relative tolerance of $10^{-3}$ for Newton’s method also yields a linear error growth rate and even slightly smaller discrete $L^{2}$ errors for long-time simulations.

7 Conclusion

We have analyzed entropy properties of iterative solvers in the context of time integration methods for nonlinear conservation laws. Continuing recent work on linear invariants in [6, 29], we have focused on the conservation of nonlinear functionals. In particular, we have considered combinations of space and implicit time discretization that result in entropy conservative and entropy dissipative schemes when the arising equation systems were solved exactly. In practice, the iterative solver is terminated once a tolerance is reached. This tolerance is chosen such that the iteration error is smaller than the time integration error. We have demonstrated that, in this situation, Newton’s method can result in a qualitatively wrong behavior of the entropy, both for entropy conservative and dissipative schemes.

Based on an analysis for Burgers’ equation, we have explored several possible entropy fixes for Newton’s method. Of these, an idea stemming from the recently developed relaxation methods, is most performant. These methods are designed as small modifications of time integrations schemes that are able to preserve the correct evolution of nonlinear functionals. Here we have shown that as long as the iteration error is small enough in the sense described above, they can also be used within implicit time integrators with inexact solves. We have demonstrated that Newton’s method with inexact linear solves and reasonable tolerances combines well with the relaxation approach, in particular for nonlinear dispersive wave equations. The numerical results show that, for the problems considered here, entropy conservation leads to smaller errors than non-conservative methods, even when the tolerance of the iterative method is an order of magnitude larger.

Acknowledgments

We thank Gregor Gassner for stimulating discussions about this research topic and comments on an early draft of the manuscript.

Statements and Declarations

Funding

Viktor Linders was partially funded by The Royal Physiographic Society in Lund.

Hendrik Ranocha was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, project number 513301895) and the Daimler und Benz Stiftung (Daimler and Benz foundation, project number 32-10/22).

Code availability

We have set up a reproducibility repository [32] for this article, containing all Julia source code required to fully reproduce the numerical experiments discussed in this article.

Appendix A Entropy analysis for Lobatto IIIC

The purpose of this appendix is to provide details of the entropic behavior of the Lobatto IIIC method used in the numerical experiments in Section 2. The analysis will be kept general enough to encompass methods of arbitrary order of accuracy. It utilizes the Summation-By-Parts (SBP) property [36, 7] satisfied by the RK method. The same argument can be made for any RK method associated with the SBP framework such as Radau IA and IIA [39]. See [31, 50] for details about the properties and implementation of such methods.

A.1 The fully discrete scheme

For a system of ordinary differential equations $\mathbf{\boldsymbol{\mathsf{u}}}_{t}=\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{u}}})$, the Runge-Kutta stage equations and update are given by

[TABLE]

In the experiment in Section 2, the 3-stage Lobatto IIIC method is used, which is associated with the Butcher tableau

[TABLE]

It will be convenient for our purposes to express (22) in a vector format as

[TABLE]

Here, $\mathbf{\boldsymbol{A}}=\{a_{ij}\}_{i,j=1}^{s}$ is the Butcher coefficient matrix, $\mathbf{\boldsymbol{1}}$ is the $s$-element vector of all ones, $\mathbf{\boldsymbol{\mathsf{U}}}$ is the (column) vector containing the stacked stage vectors $\mathbf{\boldsymbol{\mathsf{y}}}^{1},\dots,\mathbf{\boldsymbol{\mathsf{y}}}^{s}$ and $\otimes$ denotes the Kronecker product.

The update $\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}$ can be computed as in (23) due to the fact that the final row of $\mathbf{\boldsymbol{A}}$ is identical to the vector $\mathbf{\boldsymbol{b}}^{\top}$. This property can be generalized by considering methods for which there is a vector $\mathbf{\boldsymbol{\mathsf{v}}}\in\mathbb{R}^{s}$ such that $\mathbf{\boldsymbol{A}}^{\top}\mathbf{\boldsymbol{v}}=\mathbf{\boldsymbol{b}}$, in which case $\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}=(\mathbf{\boldsymbol{v}}^{\top}\otimes\mathbf{\boldsymbol{\mathsf{I}}})\mathbf{\boldsymbol{\mathsf{U}}}$. For the 3-stage Lobatto IIIC method, we have $\mathbf{\boldsymbol{v}}=(0,0,1)^{\top}$.

We now revisit the spatial discretization of Burgers’ equation in (3) by setting $\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{u}}})=-2(\mathbf{\boldsymbol{\mathsf{D}}}\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{u}}})\mathbf{\boldsymbol{\mathsf{u}}}+\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{u}}})\mathbf{\boldsymbol{\mathsf{D}}}\mathbf{\boldsymbol{\mathsf{u}}})$. The entropy behavior of the fully discrete scheme can be analyzed with the aid of the SBP property, which can be expressed as

[TABLE]

Here, $\mathbf{\boldsymbol{B}}=\mathrm{diag}(\mathbf{\boldsymbol{b}})$ and $\mathbf{\boldsymbol{e}}_{j}$ denotes the $j$th column of the $s\times s$ identity matrix. In this context, $\mathbf{\boldsymbol{A}}^{-1}$ can be viewed as a difference operator adjoined with an initial condition, and $\mathbf{\boldsymbol{B}}$ as a quadrature rule [30]. The SBP property (24) is thus a discrete version of integration by parts.

We begin by left-multiplying the stage equations in (23) by $\Delta x\mathbf{\boldsymbol{\mathsf{U}}}^{\top}(\mathbf{\boldsymbol{B}}\mathbf{\boldsymbol{A}}^{-1}\otimes\mathbf{\boldsymbol{\mathsf{I}}})$. This is a well-defined operation since $\mathbf{\boldsymbol{A}}$ is invertible for any $s$ [28]. By a derivation identical to (6), it holds that $(\mathbf{\boldsymbol{\mathsf{y}}}^{i})^{\top}\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{y}}}^{i})=\mathbf{\boldsymbol{\mathsf{0}}}$ for each $i=1,\dots,s$. Since $\mathbf{\boldsymbol{B}}$ is diagonal it follows that $\Delta x\mathbf{\boldsymbol{\mathsf{U}}}^{\top}(\mathbf{\boldsymbol{B}}\otimes\mathbf{\boldsymbol{\mathsf{I}}})\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{U}}})=\mathbf{\boldsymbol{\mathsf{0}}}$, and consequently

[TABLE]

The left-hand side of (25) is a quadratic form and is therefore equal to its symmetric part. The right-hand side is simplified by the relation $\mathbf{\boldsymbol{B}}\mathbf{\boldsymbol{A}}^{-1}\mathbf{\boldsymbol{1}}=\mathbf{\boldsymbol{e}}_{1}$; see [31, Lemma 3]. The identity (25) therefore reduces to

[TABLE]

Simplification using the SBP property (24) allows us to express this in terms of the first and last stages as

[TABLE]

By adding and subtracting $\frac{1}{2}\|\mathbf{\boldsymbol{\mathsf{u}}}^{n}\|^{2}$ from the right-hand side, the entropy, given by $\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n+1})=\frac{1}{2}\|\mathbf{\boldsymbol{\mathsf{y}}}^{s}\|^{2}$, can be expressed as

[TABLE]

Lobatto IIIC therefore dissipates entropy. This result is independent of the number of stages and holds more generally for RK methods with the SBP property. A slight generalization (to the so called gSBP property [15]) is necessary for the analysis to hold for some RK methods such as the Radau family.

A.2 Newton’s method

The Jacobian $\mathbf{\boldsymbol{\mathsf{F}}}^{\prime}$ is explicitly given by

[TABLE]

where $\mathbf{\boldsymbol{\mathsf{f}}}^{\prime}$ is the Jacobian of the spatial discretization. The form of $\mathbf{\boldsymbol{\mathsf{f}}}^{\prime}$ is identical to that seen for the midpoint rule, except that it is now repeated for each stage $\mathbf{\boldsymbol{\mathsf{y}}}^{i}$ in a block-diagonal fashion. This leads to an equation for $\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}$ of the form

[TABLE]

The matrix $\tilde{\mathbf{\boldsymbol{\mathsf{M}}}}$ is block-diagonal with each block identical to the case for the midpoint rule, but evaluated at the individual stages. As before, (27) represents a linearization of the fully discrete scheme around the iterate $\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}$, perturbed by a term arising from the spatial Jacobian.

By an analysis completely analogous to that of in the previous subsection, the entropy relation evaluates to

[TABLE]

where $\mathbf{\boldsymbol{\mathsf{M}}}=(\mathbf{\boldsymbol{B}}\otimes\mathbf{\boldsymbol{\mathsf{D}}})\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})+\mathrm{diag}((\mathbf{\boldsymbol{B}}\otimes\mathbf{\boldsymbol{\mathsf{D}}})\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})$. The entropy error induced by Newton’s method is thus of the same form as for the midpoint rule (for which $\mathbf{\boldsymbol{B}}=1$), but includes all $s$ stages of the Runge-Kutta scheme.