# Resolving Entropy Growth from Iterative Methods

Viktor Linders${}^{\mbox{\tiny\rm 1}}$, Hendrik Ranocha${}^{\mbox{\tiny\rm
2}}$, Philipp Birken${}^{\mbox{\tiny\rm 1}}$

###### 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.

${}^{\mbox{\tiny\rm 1}}$Centre for mathematical sciences, Lund University,
Lund, Sweden
${}^{\mbox{\tiny\rm 2}}$Applied Mathematics, University of Hamburg, Hamburg,
Germany
email<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
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:

$u_{t}={\cal L}(u),\quad x\in(x_{\mathrm{min}},x_{\mathrm{max}}],\quad
t\in[t_{0},t_{e}].$ (1)

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,

$u_{t}+6uu_{x}=0,\quad u(x,t=0)=u_{0}.$ (2)

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

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

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

$\mathbf{\boldsymbol{\mathsf{u}}}_{t}+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}}})=\mathbf{\boldsymbol{\mathsf{0}}}.$
(3)

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

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\frac{1}{2}\|\mathbf{\boldsymbol{\mathsf{u}}}\|^{2}$
$\displaystyle=-2\Delta
x(\mathbf{\boldsymbol{\mathsf{u}}}^{\top}\mathbf{\boldsymbol{\mathsf{D}}}\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{u}}})\mathbf{\boldsymbol{\mathsf{u}}}+\mathbf{\boldsymbol{\mathsf{u}}}^{\top}\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{u}}})\mathbf{\boldsymbol{\mathsf{D}}}\mathbf{\boldsymbol{\mathsf{u}}})$
$\displaystyle=-2\Delta
x(\mathbf{\boldsymbol{\mathsf{u}}}^{\top}\mathbf{\boldsymbol{\mathsf{D}}}\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{u}}})\mathbf{\boldsymbol{\mathsf{u}}}+\mathbf{\boldsymbol{\mathsf{u}}}^{\top}\mathbf{\boldsymbol{\mathsf{D}}}^{\top}\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{u}}})\mathbf{\boldsymbol{\mathsf{u}}})$
$\displaystyle=-2\Delta
x(\mathbf{\boldsymbol{\mathsf{u}}}^{\top}(\mathbf{\boldsymbol{\mathsf{D}}}+\mathbf{\boldsymbol{\mathsf{D}}}^{\top})\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{u}}})\mathbf{\boldsymbol{\mathsf{u}}}=\mathbf{\boldsymbol{\mathsf{0}}}.$

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:

$\displaystyle\mathbf{\boldsymbol{\mathsf{U}}}$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\frac{\Delta
t_{n}}{2}\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{U}}})$
$\displaystyle\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\Delta
t_{n}\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{U}}})\equiv
2\mathbf{\boldsymbol{\mathsf{U}}}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}.$

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,

$\displaystyle\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}})$
$\displaystyle:=\mathbf{\boldsymbol{\mathsf{U}}}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\Delta
t_{n}(\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}}})=\mathbf{\boldsymbol{\mathsf{0}}},$
(4) $\displaystyle\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}$
$\displaystyle\phantom{:}=2\mathbf{\boldsymbol{\mathsf{U}}}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}.$

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

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

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

$\displaystyle\mathbf{\boldsymbol{\mathsf{U}}}^{\top}\mathbf{\boldsymbol{\mathsf{U}}}-\mathbf{\boldsymbol{\mathsf{U}}}^{\top}\mathbf{\boldsymbol{\mathsf{u}}}^{n}$
$\displaystyle=-\Delta
t_{n}\mathbf{\boldsymbol{\mathsf{U}}}^{\top}\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{U}}})$
(6) $\displaystyle=-\Delta
t_{n}\mathbf{\boldsymbol{\mathsf{U}}}^{\top}(\mathbf{\boldsymbol{\mathsf{D}}}\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}})+\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}})\mathbf{\boldsymbol{\mathsf{D}}})\mathbf{\boldsymbol{\mathsf{U}}}$
$\displaystyle=-\Delta
t_{n}\mathbf{\boldsymbol{\mathsf{U}}}^{\top}((\mathbf{\boldsymbol{\mathsf{D}}}+\mathbf{\boldsymbol{\mathsf{D}}}^{\top})\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}}))\mathbf{\boldsymbol{\mathsf{U}}}=0,$

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

$\displaystyle\mathbf{\boldsymbol{\mathsf{F}}}^{\prime}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\Delta\mathbf{\boldsymbol{\mathsf{U}}}+\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{0}}},$ (7)
$\displaystyle\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}+\Delta\mathbf{\boldsymbol{\mathsf{U}}},\quad
k=0,1,\dots$

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

$\mathbf{\boldsymbol{\mathsf{F}}}^{\prime}({\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}})=\mathbf{\boldsymbol{\mathsf{I}}}+\Delta
t_{n}\left(\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\mathbf{\boldsymbol{\mathsf{D}}}+\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{D}}}\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})+2\mathbf{\boldsymbol{\mathsf{D}}}\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\right).$
(8)

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)}$:

$\displaystyle\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}$
$\displaystyle+\Delta
t_{n}\left(\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\mathbf{\boldsymbol{\mathsf{D}}}\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}+\mathbf{\boldsymbol{\mathsf{D}}}\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}\right)$
(9) $\displaystyle+\Delta
t_{n}\underbrace{\left[\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{D}}}\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})+\mathbf{\boldsymbol{\mathsf{D}}}\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\right]}_{\mathbf{\boldsymbol{\mathsf{M}}}}\Delta\mathbf{\boldsymbol{\mathsf{U}}}=\mathbf{\boldsymbol{\mathsf{0}}}.$

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

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

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,

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

and hence

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

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

$\displaystyle\mathbf{\boldsymbol{\mathsf{M}}}^{\top}\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}$
$\displaystyle=\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{D}}}\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}+\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\mathbf{\boldsymbol{\mathsf{D}}}^{\top}\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}$
$\displaystyle=\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{D}}}\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\mathbf{\boldsymbol{\mathsf{1}}}-\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{D}}}\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\mathbf{\boldsymbol{\mathsf{1}}}=\mathbf{\boldsymbol{\mathsf{0}}}.$

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

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

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

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

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

$u_{0}=\frac{c}{2}\operatorname{sech}^{2}\left(\frac{\sqrt{c}}{2}x\right),$
(11)

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.

(a) Entropy dissipative 3-stage Lobatto IIIC.

(b) Entropy conservative implicit midpoint rule.

Figure 1: Entropy
$\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n})=\|\mathbf{\boldsymbol{\mathsf{u}}}^{n}\|/2$
and residual
$\|\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\|$
after $n=1$ time steps with $k$ Newton iterations applied to Burgers’ equation
with $\Delta x=0.1$ and $\Delta t=0.1$.

##  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

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

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.

Figure 2: Entropy $\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n})$ and residual
$\|\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\|$
after $n=1$ time steps with $k\in\\{0,\dots,14\\}$ iterations of method of
Newton-type for the implicit midpoint rule applied to Burgers’ equation with
$\Delta x=0.1$ and $\Delta t=0.5$. The convergence is linear.

###  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

$\displaystyle\mathbf{\boldsymbol{\mathsf{F}}}^{\prime}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})(\tilde{\mathbf{\boldsymbol{\mathsf{U}}}}-\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})+\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{0}}},$ (12)
$\displaystyle\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}$
$\displaystyle=\alpha_{k}\tilde{\mathbf{\boldsymbol{\mathsf{U}}}}+(1-\alpha_{k})\mathbf{\boldsymbol{\mathsf{U}}}^{(k)},$

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

$\displaystyle\mathbf{\boldsymbol{\mathsf{F}}}^{\prime}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\Delta\mathbf{\boldsymbol{\mathsf{U}}}+\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})$
$\displaystyle=(1-\alpha_{k})\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}),$
$\displaystyle\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}+\Delta\mathbf{\boldsymbol{\mathsf{U}}}.$

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

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

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

$\alpha_{k}=-\frac{\langle\Delta\mathbf{\boldsymbol{\mathsf{U}}},\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}\rangle+\langle\tilde{\mathbf{\boldsymbol{\mathsf{U}}}},\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}\rangle}{\|\Delta\mathbf{\boldsymbol{\mathsf{U}}}\|^{2}}.$

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.

Figure 3: Entropy $\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n})$ and residual
$\|\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})\|$
after $n=1$ time steps with $k\in\\{0,\dots,6\\}$ inexact Newton iterations
for the implicit midpoint rule applied to Burgers’ equation with $\Delta
x=0.1$ and $\Delta t=0.5$. The convergence is linear.

##  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

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

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

$\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}t\\\
\mathbf{\boldsymbol{\mathsf{u}}}(t)\\\
\eta\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}(t)\bigr{)}\end{pmatrix}=\begin{pmatrix}1\\\
f\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}(t)\bigr{)}\\\
\eta^{\prime}\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}(t)\bigr{)}f\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}(t)\bigr{)}\end{pmatrix}.$

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

$\begin{pmatrix}t^{n+1}_{\gamma}\\\
\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}_{\gamma}\\\
\eta\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}_{\gamma}\bigr{)}\end{pmatrix}=\begin{pmatrix}t^{n}\\\
\mathbf{\boldsymbol{\mathsf{u}}}^{n}\\\
\eta\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}^{n}\bigr{)}\end{pmatrix}+\gamma\Delta
t\begin{pmatrix}1\\\
\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}\\\
\eta^{\mathrm{new}}-\eta\bigl{(}\mathbf{\boldsymbol{\mathsf{u}}}^{n}\bigr{)}\end{pmatrix}$
(13)

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

$\displaystyle\mathbf{\boldsymbol{\mathsf{y}}}^{i}$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\Delta
t\sum_{j=1}^{s}a_{ij}\,f\bigl{(}\mathbf{\boldsymbol{\mathsf{y}}}^{j}\bigr{)},\qquad
i\in\\{1,\dots,s\\},$ $\displaystyle\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\Delta
t\sum_{i=1}^{s}b_{i}\,f\bigl{(}\mathbf{\boldsymbol{\mathsf{y}}}^{i}\bigr{)},$

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

$\eta^{\mathrm{new}}=\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n})+\Delta
t\sum_{i=1}^{s}b_{i}(\eta^{\prime}f)(\mathbf{\boldsymbol{\mathsf{y}}}^{i}).$

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

$\eta^{\prime}(\mathbf{\boldsymbol{\mathsf{u}}}^{n+1})\frac{\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}}{\|\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}\|}=c\Delta
t+\mathcal{O}(\Delta t^{2}),$

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

$\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}_{\gamma}=\mathbf{\boldsymbol{\mathsf{u}}}(t^{n+1}_{\gamma})+\mathcal{O}(\Delta
t^{p+1}).$

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

$u_{t}+6uu_{x}+u_{xxx}=0,\quad u(x,t=0)=u_{0}.$ (14)

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

$u(x,t)=\frac{c}{2}\operatorname{sech}^{2}\left(\frac{\sqrt{c}}{2}(x-ct)\right).$
(15)

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

$\gamma=\frac{\|\mathbf{\boldsymbol{\mathsf{u}}}^{n}\|^{2}-\langle\mathbf{\boldsymbol{\mathsf{U}}},\mathbf{\boldsymbol{\mathsf{u}}}^{n}\rangle}{\|\mathbf{\boldsymbol{\mathsf{U}}}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}\|^{2}}.$

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.

Figure 4: Error of numerical solutions and entropy error of the KdV equation
(14) with periodic boundary conditions using the implicit midpoint rule with
time step size $\Delta t=0.05$. The implicit equations are solved with Newton-
GMRES and different relative tolerances.

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: Error of numerical solutions and entropy error of the KdV equation
(14) with periodic boundary conditions using the fourth-order, three-stage
Lobatto IIIC method with time step size $\Delta t=0.1$. The implicit equations
are solved with Newton-GMRES and different absolute and relative tolerances.

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

$u_{t}+u_{x}+uu_{x}-u_{txx}=0.$ (16)

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

$\displaystyle J^{\text{BBM}}_{1}(u)$ $\displaystyle=\int u,$ (17)
$\displaystyle J^{\text{BBM}}_{2}(u)$
$\displaystyle=\frac{1}{2}\int(u^{2}+(u_{x})^{2})=\frac{1}{2}\int
u(\operatorname{I}-\partial_{x}^{2})u,$ $\displaystyle J^{\text{BBM}}_{3}(u)$
$\displaystyle=\int(u+1)^{3},$

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

$\mathbf{\boldsymbol{\mathsf{u}}}_{t}+(\mathbf{\boldsymbol{\mathsf{I}}}-\mathbf{\boldsymbol{\mathsf{D}}}_{2})^{-1}\left(\frac{1}{3}\mathbf{\boldsymbol{\mathsf{D}}}_{1}\mathbf{\boldsymbol{\mathsf{u}}}^{2}+\frac{1}{3}\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{u}}})\mathbf{\boldsymbol{\mathsf{D}}}_{1}\mathbf{\boldsymbol{\mathsf{u}}}+\mathbf{\boldsymbol{\mathsf{D}}}_{1}\mathbf{\boldsymbol{\mathsf{u}}}\right)=\mathbf{\boldsymbol{\mathsf{0}}}$
(18)

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

$J^{\text{BBM}}_{1}(\mathbf{\boldsymbol{\mathsf{u}}})=\mathbf{\boldsymbol{\mathsf{1}}}^{\top}\mathbf{\boldsymbol{\mathsf{M}}}\mathbf{\boldsymbol{\mathsf{u}}},\qquad
J^{\text{BBM}}_{2}(\mathbf{\boldsymbol{\mathsf{u}}})=\mathbf{\boldsymbol{\mathsf{u}}}^{\top}(\mathbf{\boldsymbol{\mathsf{I}}}-\mathbf{\boldsymbol{\mathsf{D}}}_{2})\mathbf{\boldsymbol{\mathsf{u}}}.$

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

$\mathbf{\boldsymbol{\mathsf{u}}}_{t}+(\mathbf{\boldsymbol{\mathsf{I}}}-\mathbf{\boldsymbol{\mathsf{D}}}_{2})^{-1}\mathbf{\boldsymbol{\mathsf{D}}}_{1}\left(\frac{1}{2}\mathbf{\boldsymbol{\mathsf{u}}}^{2}+\mathbf{\boldsymbol{\mathsf{u}}}\right)=\mathbf{\boldsymbol{\mathsf{0}}}$
(19)

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

$J^{\text{BBM}}_{3}(\mathbf{\boldsymbol{\mathsf{u}}})=\mathbf{\boldsymbol{\mathsf{1}}}^{\top}\mathbf{\boldsymbol{\mathsf{M}}}(\mathbf{\boldsymbol{\mathsf{1}}}+\mathbf{\boldsymbol{\mathsf{u}}})^{3}$

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

$\mathbf{\boldsymbol{\mathsf{1}}}^{\top}\mathbf{\boldsymbol{\mathsf{M}}}\mathbf{\boldsymbol{\mathsf{u}}}_{t}=-\mathbf{\boldsymbol{\mathsf{1}}}^{\top}\mathbf{\boldsymbol{\mathsf{M}}}(\mathbf{\boldsymbol{\mathsf{I}}}-\mathbf{\boldsymbol{\mathsf{D}}}_{2})^{-1}\mathbf{\boldsymbol{\mathsf{D}}}_{1}\left(\frac{1}{2}\mathbf{\boldsymbol{\mathsf{u}}}^{2}+\mathbf{\boldsymbol{\mathsf{u}}}\right)=-\mathbf{\boldsymbol{\mathsf{1}}}^{\top}\mathbf{\boldsymbol{\mathsf{D}}}_{1}\left(\frac{1}{2}\mathbf{\boldsymbol{\mathsf{u}}}^{2}+\mathbf{\boldsymbol{\mathsf{u}}}\right)=0.$

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

$\mathcal{H}(\mathbf{\boldsymbol{\mathsf{u}}})=\mathbf{\boldsymbol{\mathsf{1}}}^{\top}\mathbf{\boldsymbol{\mathsf{M}}}\left(\frac{1}{6}\mathbf{\boldsymbol{\mathsf{u}}}^{3}+\frac{1}{2}\mathbf{\boldsymbol{\mathsf{u}}}^{2}\right).$

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

$\mathcal{H}_{t}(\mathbf{\boldsymbol{\mathsf{u}}})=\mathcal{H}^{\prime}(\mathbf{\boldsymbol{\mathsf{u}}})\mathbf{\boldsymbol{\mathsf{u}}}_{t}=-\left(\frac{1}{6}\mathbf{\boldsymbol{\mathsf{u}}}^{3}+\frac{1}{2}\mathbf{\boldsymbol{\mathsf{u}}}^{2}\right)^{\top}\mathbf{\boldsymbol{\mathsf{M}}}(\mathbf{\boldsymbol{\mathsf{I}}}-\mathbf{\boldsymbol{\mathsf{D}}}_{2})^{-1}\mathbf{\boldsymbol{\mathsf{D}}}_{1}\left(\frac{1}{6}\mathbf{\boldsymbol{\mathsf{u}}}^{3}+\frac{1}{2}\mathbf{\boldsymbol{\mathsf{u}}}^{2}\right)=0.$

∎

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

$\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}=\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\Delta
t\int_{0}^{1}f\bigl{(}s\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}+(1-s)\mathbf{\boldsymbol{\mathsf{u}}}^{n}\bigr{)}\mathrm{d}s.$

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

$\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}=\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\frac{\Delta
t}{6}\left(f(\mathbf{\boldsymbol{\mathsf{u}}}^{n})+4f\biggl{(}\frac{\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}+\mathbf{\boldsymbol{\mathsf{u}}}^{n}}{2}\biggr{)}+f(\mathbf{\boldsymbol{\mathsf{u}}}^{n+1})\right)$
(20)

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

$u(t,x)=A\operatorname{sech}\bigl{(}K(x-ct)\bigr{)}^{2},\quad A=3(c-1),\quad
K=\frac{1}{2}\sqrt{1-1/c},$ (21)

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.

(a) Discrete $L^{2}$ and entropy error for the quadratic invariant
$J_{2}^{BBM}(\mathbf{\boldsymbol{\mathsf{u}}})$. The implicit midpoint rule is
applied to the semidiscretization (18).

(b) Discrete $L^{2}$ and entropy error for the cubic invariant
$J_{3}^{BBM}(\mathbf{\boldsymbol{\mathsf{u}}})$. The AVF method is applied to
the semidiscretization (19).

Figure 6: Discrete $L^{2}$ error of numerical methods for the BBM equation
(16) with periodic boundary conditions. The implicit equations are solved with
Newton-GMRES and different relative tolerances.

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

$\displaystyle\mathbf{\boldsymbol{\mathsf{y}}}^{i}$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\Delta
t_{n}\sum_{j=1}^{s}a_{ij}\,\mathbf{\boldsymbol{\mathsf{f}}}\bigl{(}\mathbf{\boldsymbol{\mathsf{y}}}^{j}\bigr{)},\qquad
i\in\\{1,\dots,s\\},$ (22)
$\displaystyle\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\Delta
t_{n}\sum_{i=1}^{s}b_{i}\,\mathbf{\boldsymbol{\mathsf{f}}}\bigl{(}\mathbf{\boldsymbol{\mathsf{y}}}^{i}\bigr{)}.$

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

$\begin{array}[]{c|ccc}0&1/6&-1/3&\phantom{-}1/6\\\
1/2&1/6&\phantom{-}5/12&-1/12\\\ 1&1/6&\phantom{-}2/3&\phantom{-}1/6\\\
\hline\cr&1/6&\phantom{-}2/3&\phantom{-}1/6\end{array}.$

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

$\displaystyle\mathbf{\boldsymbol{\mathsf{F}}}(\mathbf{\boldsymbol{\mathsf{U}}}):$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{U}}}-\mathbf{\boldsymbol{1}}\otimes\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\Delta
t_{n}(\mathbf{\boldsymbol{A}}\otimes\mathbf{\boldsymbol{\mathsf{I}}})\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{U}}})=\mathbf{\boldsymbol{\mathsf{0}}},$
(23) $\displaystyle\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}$
$\displaystyle=\mathbf{\boldsymbol{\mathsf{y}}}^{s}.$

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

$\mathbf{\boldsymbol{B}}\mathbf{\boldsymbol{A}}^{-1}+\mathbf{\boldsymbol{A}}^{-\top}\mathbf{\boldsymbol{B}}=\mathrm{diag}(1,0,\dots,0,1)=\mathbf{\boldsymbol{e}}_{1}\mathbf{\boldsymbol{e}}_{1}^{\top}+\mathbf{\boldsymbol{e}}_{s}\mathbf{\boldsymbol{e}}_{s}^{\top}.$
(24)

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

$\Delta
x\mathbf{\boldsymbol{\mathsf{U}}}^{\top}(\mathbf{\boldsymbol{B}}\mathbf{\boldsymbol{A}}^{-1}\otimes\mathbf{\boldsymbol{\mathsf{I}}})\mathbf{\boldsymbol{\mathsf{U}}}=\Delta
x\mathbf{\boldsymbol{\mathsf{U}}}^{\top}(\mathbf{\boldsymbol{B}}\mathbf{\boldsymbol{A}}^{-1}\mathbf{\boldsymbol{1}}\otimes\mathbf{\boldsymbol{\mathsf{u}}}^{n}).$
(25)

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

$\Delta
x\mathbf{\boldsymbol{\mathsf{U}}}^{\top}\left(\frac{\mathbf{\boldsymbol{B}}\mathbf{\boldsymbol{A}}^{-1}+\mathbf{\boldsymbol{A}}^{-\top}\mathbf{\boldsymbol{B}}}{2}\otimes\mathbf{\boldsymbol{\mathsf{I}}}\right)\mathbf{\boldsymbol{\mathsf{U}}}=\Delta
x\mathbf{\boldsymbol{\mathsf{U}}}^{\top}(\mathbf{\boldsymbol{e}}_{1}\otimes\mathbf{\boldsymbol{\mathsf{u}}}^{n}).$

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

$\frac{1}{2}\|\mathbf{\boldsymbol{\mathsf{y}}}^{1}\|^{2}+\frac{1}{2}\|\mathbf{\boldsymbol{\mathsf{y}}}^{s}\|^{2}=\Delta
x(\mathbf{\boldsymbol{\mathsf{y}}}^{1})^{\top}\mathbf{\boldsymbol{\mathsf{u}}}^{n}.$

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

$\displaystyle\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n+1})$
$\displaystyle=\frac{1}{2}\|\mathbf{\boldsymbol{\mathsf{u}}}^{n}\|^{2}-\frac{1}{2}\|\mathbf{\boldsymbol{\mathsf{y}}}^{1}\|^{2}+\Delta
x(\mathbf{\boldsymbol{\mathsf{y}}}^{1})^{\top}\mathbf{\boldsymbol{\mathsf{u}}}^{n}-\frac{1}{2}\|\mathbf{\boldsymbol{\mathsf{u}}}^{n}\|^{2}$
$\displaystyle=\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n})-\eta(\mathbf{\boldsymbol{\mathsf{y}}}^{1}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}).$

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

$\mathbf{\boldsymbol{\mathsf{F}}}^{\prime}({\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}})=\mathbf{\boldsymbol{\mathsf{I}}}+\Delta
t_{n}(\mathbf{\boldsymbol{A}}\otimes\mathbf{\boldsymbol{\mathsf{I}}})\mathbf{\boldsymbol{\mathsf{f}}}^{\prime}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)}),$
(26)

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

$\mathbf{\boldsymbol{\mathsf{U}}}^{(k+1)}-\mathbf{\boldsymbol{1}}\otimes\mathbf{\boldsymbol{\mathsf{u}}}^{n}+\Delta
t_{n}(\mathbf{\boldsymbol{A}}\otimes\mathbf{\boldsymbol{\mathsf{I}}})\mathbf{\boldsymbol{\mathsf{f}}}(\mathbf{\boldsymbol{\mathsf{U}}}^{(k)})+2\Delta
t_{n}(\mathbf{\boldsymbol{A}}\otimes\mathbf{\boldsymbol{\mathsf{I}}})\tilde{\mathbf{\boldsymbol{\mathsf{M}}}}\Delta\mathbf{\boldsymbol{\mathsf{U}}}=\mathbf{\boldsymbol{\mathsf{0}}}.$
(27)

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

$\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n+1})=\eta(\mathbf{\boldsymbol{\mathsf{u}}}^{n})-\eta(\mathbf{\boldsymbol{\mathsf{y}}}^{1}-\mathbf{\boldsymbol{\mathsf{u}}}^{n})-2\Delta
t_{n}\Delta
x\Delta\mathbf{\boldsymbol{\mathsf{U}}}^{\top}\mathbf{\boldsymbol{\mathsf{M}}}\Delta\mathbf{\boldsymbol{\mathsf{U}}},$

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.

## References

  * [1] A. Araújo and A. Durán, Error propagation in the numerical integration of solitary waves. The regularized long wave equation, Appl. Numer. Math., 36 (2001), pp. 197–217.
  * [2] M. J. Bencomo and J. Chan, Discrete adjoint computations for relaxation Runge-Kutta methods, J. Sci. Comput., 94 (2023), p. 59.
  * [3] T. B. Benjamin, J. L. Bona, and J. J. Mahony, Model equations for long waves in nonlinear dispersive systems, Philos. Trans. Royal Soc. A PHILOS T R SOC A, 272 (1972), pp. 47–78.
  * [4] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Rev., 59 (2017), pp. 65–98.
  * [5] P. Birken, Numerical Methods for Unsteady Compressible Flow Problems, Chapman and Hall/CRC, New York, 2021.
  * [6] P. Birken and V. Linders, Conservation properties of iterative methods for implicit discretizations of conservation laws, J. Sci. Comput., 92 (2022), pp. 1–32.
  * [7] P. D. Boom and D. W. Zingg, High-order implicit time-marching methods based on generalized summation-by-parts operators, SIAM J. Sci. Comput., 37 (2015), pp. A2682–A2709.
  * [8] J. C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons Ltd, Chichester, 2016.
  * [9] M. Calvo, D. Hernández-Abreu, J. I. Montijano, and L. Rández, On the preservation of invariants by explicit Runge-Kutta methods, SIAM J. Sci. Comput., 28 (2006), pp. 868–885.
  * [10] E. Celledoni, R. I. McLachlan, D. I. McLaren, B. Owren, G. R. W. Quispel, and W. Wright, Energy-preserving Runge-Kutta methods, ESAIM-Math. Model. Num., 43 (2009), pp. 645–649.
  * [11] J. Chan and C. G. Taylor, Efficient computation of Jacobian matrices for entropy stable summation-by-parts schemes, J. Comput. Phys., 448 (2022), p. 110701.
  * [12] J. De Frutos and J. M. Sanz-Serna, Accuracy and conservation properties in numerical integration: the case of the Korteweg-de Vries equation, Numer. Math., 75 (1997), pp. 421–445.
  * [13] K. Dekker and J. G. Verwer, Stability of Runge-Kutta methods for stiff nonlinear differential equations, vol. 2 of CWI Monographs, North-Holland, Amsterdam, 1984.
  * [14] S. C. Eisenstat and H. F. Walker, Choosing the forcing terms in an inexact Newton method, SIAM J. Sci. Comput., 17 (1996), pp. 16–32.
  * [15] D. C. D. R. Fernández, P. D. Boom, and D. W. Zingg, A generalized framework for nodal first derivative summation-by-parts operators, J. Comput. Phys., 266 (2014), pp. 214–239.
  * [16] T. C. Fisher and M. H. Carpenter, High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains, J. Comp. Phys., 252 (2013), pp. 518–557.
  * [17] L. Friedrich, G. Schnücke, A. R. Winters, D. C. Del Rey Fernández, G. J. Gassner, and M. H. Carpenter, Entropy Stable Space–Time Discontinuous Galerkin Schemes with Summation-by-Parts Property for Hyperbolic Conservation Laws, J. Sci. Comput., 80 (2019), pp. 175–222.
  * [18] G. J. Gassner and A. R. Winters, A novel robust strategy for discontinuous Galerkin methods in computational fluid mechanics: why? when? what? where?, Front. Phys., 8 (2021), p. 500690.
  * [19] E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, vol. 31 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin Heidelberg, 2006.
  * [20] J. Jackaman and S. MacLachlan, Preconditioned Krylov solvers for structure-preserving discretisations, 2022, arXiv:2212.05127.
  * [21] S. Kang and E. M. Constantinescu, Entropy-preserving and entropy-stable relaxation IMEX and multirate time-stepping methods, J. Sci. Comput., 93 (2022), p. 23.
  * [22] C. T. Kelley, Iterative methods for linear and nonlinear equations, SIAM, Philadelphia, 1995.
  * [23] D. I. Ketcheson, Relaxation Runge-Kutta methods: Conservation and stability for inner-product norms, SIAM J. Numer. Anal., 57 (2019), pp. 2850–2870.
  * [24] D. A. Kopriva, Implementing spectral methods for partial differential equations: Algorithms for scientists and engineers, Springer Science & Business Media, Dordrecht, 2009.
  * [25] D. Kröner and M. Ohlberger, A posteriori error estimates for upwind finite volume schemes for nonlinear conservation laws in multi dimensions, Math. Comp., 69 (1999), pp. 25–39.
  * [26] D. Li, X. Li, and Z. Zhang, Implicit-explicit relaxation Runge-Kutta methods: construction, analysis and applications to PDEs, Math. Comput., (2022).
  * [27] , Linearly implicit and high-order energy-preserving relaxation schemes for highly oscillatory Hamiltonian systems, J. Comput. Phys., (2023), p. 111925.
  * [28] V. Linders, On an eigenvalue property of summation-by-parts operators, J. Sci. Comput., 93 (2022), p. 82.
  * [29] V. Linders and P. Birken, Locally conservative and flux consistent iterative methods, 2022, arXiv:2206.10943.
  * [30] V. Linders, T. Lundquist, and J. Nordström, On the order of accuracy of finite difference operators on diagonal norm based summation-by-parts form, SIAM J. Numer. Anal., 56 (2018), pp. 1048–1063.
  * [31] V. Linders, J. Nordström, and S. H. Frankel, Properties of Runge-Kutta-summation-by-parts methods, J. Comput. Phys., 419 (2020), p. 109684.
  * [32] V. Linders, H. Ranocha, and P. Birken, Reproducibility repository for ”Resolving entropy growth from iterative methods”.  https://gitlab.maths.lu.se/Linders/iteration-entropy, 2023.
  * [33] T. Lundquist and J. Nordström, The SBP-SAT technique for initial value problems, J. Comput. Phys., 270 (2014), pp. 86–104.
  * [34] R. I. McLachlan, G. Quispel, and N. Robidoux, Geometric integration using discrete gradients, Philos. Trans. Royal Soc. A PHILOS T R SOC A, 357 (1999), pp. 1021–1045.
  * [35] A. Montoison, D. Orban, and contributors, Krylov.jl: A Julia basket of hand-picked Krylov methods.  https://github.com/JuliaSmoothOptimizers/Krylov.jl, June 2020.
  * [36] J. Nordström and T. Lundquist, Summation-by-parts in time, J. Comput. Phys., 251 (2013), pp. 487–499.
  * [37] P. J. Olver, Euler operators and conservation laws of the BBM equation, in Mathematical Proceedings of the Cambridge Philosophical Society, vol. 85, Cambridge University Press, 1979, pp. 143–160.
  * [38] G. Quispel and D. I. McLaren, A new class of energy-preserving numerical integration methods, J. Phys. A Math. Theor., 41 (2008), p. 045206.
  * [39] H. Ranocha, Some notes on summation by parts time integration methods, Results Appl. Math., 1 (2019), p. 100004.
  * [40] , SummationByPartsOperators.jl: A Julia library of provably stable semidiscretization techniques with mimetic properties, J. Open Source Softw., 6 (2021), p. 3454.
  * [41] H. Ranocha, L. Dalcin, and M. Parsani, Fully-discrete explicit locally entropy-stable schemes for the compressible Euler and Navier-Stokes equations, Comput. Math. Appl., 80 (2020), pp. 1343–1359.
  * [42] H. Ranocha and D. I. Ketcheson, Relaxation Runge-Kutta methods for Hamiltonian problems, J. Sci. Comput., 84 (2020).
  * [43] H. Ranocha, L. Lóczi, and D. I. Ketcheson, General relaxation methods for initial-value problems with application to multistep schemes, Numer. Math., 146 (2020), pp. 875–906.
  * [44] H. Ranocha, D. Mitsotakis, and D. I. Ketcheson, A broad class of conservative numerical methods for dispersive wave equations, Commun. Comput. Phys., 29 (2021), pp. 979–1029.
  * [45] H. Ranocha, M. Sayyari, L. Dalcin, M. Parsani, and D. I. Ketcheson, Relaxation Runge-Kutta methods: Fully-discrete explicit entropy-stable schemes for the compressible Euler and Navier-Stokes equations, SIAM J. Sci. Comput., 42 (2020), pp. A612–A638.
  * [46] J. Revels, M. Lubin, and T. Papamarkou, Forward-mode automatic differentiation in Julia, 07 2016, arXiv:1607.07892.
  * [47] R. D. Richtmyer and K. W. Morton, Difference Methods for Boundary-Value Problems, John Wiley & Sons, New York, London, Sydney, 1967.
  * [48] J. M. Sanz-Serna, An explicit finite-difference scheme with exact conservation properties, J. Comput. Phys., 47 (1982), pp. 199–210.
  * [49] J. M. Sanz-Serna and V. Manoranjan, A method for the integration in time of certain partial differential equations, J. Comput. Phys., 52 (1983), pp. 273–289.
  * [50] L. M. Versbach, V. Linders, R. Klöfkorn, and P. Birken, Theoretical and practical aspects of space-time DG-SEM implementations, 2022, arxiv:2201.05800.