# Locally conservative and flux consistent iterative methods

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

###### Abstract

Conservation and consistency are fundamental properties of discretizations of
systems of hyperbolic conservation laws. Here, these concepts are extended to
the realm of iterative methods by formally defining _locally conservative_ and
_flux consistent_ iterations. These concepts are of both theoretical and
practical importance: Based on recent work by the authors, it is shown that
pseudo-time iterations using explicit Runge-Kutta methods are locally
conservative but not necessarily flux consistent. An extension of the Lax-
Wendroff theorem is presented, revealing convergence towards weak solutions of
a temporally retarded system of conservation laws. Each equation is modified
in the same way, namely by a particular scalar factor multiplying the spatial
flux terms. A technique for enforcing flux consistency, and thereby recovering
convergence, is presented. Further, local conservation is established for all
Krylov subspace methods, with and without restarts, and for Newton’s method
under certain assumptions on the discretization. Thus it is shown that Newton-
Krylov methods are locally conservative, although not necessarily flux
consistent. Numerical experiments with the 2D compressible Euler equations
corroborate the theoretical results. Further numerical investigations of the
impact of flux consistency on Newton-Krylov methods indicate that its effect
is case dependent, and diminishes as the number of iterations grow.

${}^{\mbox{\tiny\rm 1}}$Centre for mathematical sciences, Lund University,
Lund, Sweden.
email<EMAIL_ADDRESS>
<EMAIL_ADDRESS>

Keywords: Iterative methods, Conservation laws, Lax-Wendroff theorem, Pseudo-
time iterations, Newton-Krylov methods

##  1 Introduction

Conservation laws arise ubiquitously in the modelling of physical phenomena
and their discretizations remain the subject of intense research. Fundamental
properties of successful schemes include conservation, consistency and
convergence. These concepts are well defined for both space and time
discretizations; explicit and implicit.

Implicit discretizations typically result in a large, sparse systems of
nonlinear equations to be solved in each time step. The solution is usually
approximated through the application of iterative methods; see e.g. [2,
Chapters 5 & 6]. Yet, discussions about conservation, consistency and
convergence (in the sense of grid refinement) for schemes involving iterative
methods are rare. In [9, 1], studys were conducted of particular implicit
finite volume schemes applied to the steady Euler equations, solved using a
variety of modified Newton-type methods. The results indicate that the choice
of iterative method has a significant impact on the convergence of the scheme.
Based on these results, a study of similar schemes applied to the Reynolds-
Averaged Navier-Stokes (RANS) equations were carried out in [10], where it was
found that the less performant methods violate mass conservation.

In [4], the authors considered general finite volume discretizations of 1D
scalar conservation laws, discretized in time with the implicit Euler method.
Global (i.e. mass) conservation was proven for many methods, including pseudo-
time iterations, Krylov subspace methods, Newton’s method and certain
multigrid techniques. On the other hand, the Jacobi and Gauss-Seidel
iterations were shown to violate mass conservation in general, corroborating
the observations in [10]. A stronger notion is that of local conservation,
defined formally below, which loosely means that mass is not only conserved
but also not teleported. This notion is important for physical correcteness,
and allows to prove extensions of the Lax-Wendroff theorem [12], thus giving a
much stronger mathematical backing of such nonlinear schemes. In [4], a start
was made in this vein for pseudo-time iterations. It was shown that in case of
convergence, the resulting scheme converges to a solution of a conservation
law, where the flux is multiplied by a scheme dependent factor, unless
particular care is taken. We say that the iterative method lacks flux
consistence, which manifests as a temporal retardation.

Throughout this article, we work with systems of conservation laws. After
introducing relevant notation and terminology, we formally define _locally
conservative_ and _flux consistent_ iterative methods in 2. We extend the
results on pseudotime iterations from [4] in 3, while considering a large
class of implicit Runge-Kutta (RK) methods in place of Euler’s method. As it
turns out, even for systems flux inconsistency manifests through a scheme
dependent scalar factor. Thus, a method to enforce flux consistency, first
introduced in [4], applies here too.

The second focus of this work is local conservation of the important class of
Newton-Krylov methods. In 4, we first prove that Newton’s method is both
locally conservative and flux consistent under certain assumptions on the
spatial discretization and when solving all linear systems exactly. We can
thus establish that if there are problems with conservation within an implicit
solver using Newton’s method, they stem from the iterative solver for the
linear systems. Secondly, by relating Krylov subspace methods to pseudo-time
iterations, local conservation is shown also for these in 5. This subsequently
leads to a proof of local conservation for Newton-Krylov methods.

Numerical examples corroborate the theoretical findings in 6. We further
explore the impact of flux consistency on Newton-Krylov methods by applying
the aforementioned technique for ensuring flux consistency of pseudo-time
iterations. The results indicate that role of flux consistency is case
dependent, and that its effect diminishes as the number of iterations grow.

##  2 Preliminaries

This section introduces relevant notation and the theoretical background upon
which the remaining paper rests.

###  2.1 Notation

Scalar quantities are denoted by letters in normal font. Vectors and matrices
are bold, with vectors being lower case and matrices upper case. Vectors and
matrices of several different dimensions are treated in the manuscript:

  * •

We consider systems of $m$ conservation laws. Vectors in $\mathbb{R}^{m}$ are
written with an underline e.g. $\underline{\mathbf{\boldsymbol{u}}}$.

  * •

We consider $s$-stage explicit Runge-Kutta methods. Vectors in
$\mathbb{R}^{s}$ are written with a right-pointing arrow, e.g.
$\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}}$, and
similarly for matrices in $\mathbb{R}^{s\times s}$.

  * •

We consider $\tilde{s}$-stage implicit Runge-Kutta methods. Vectors in
$\mathbb{R}^{\tilde{s}}$ are written with a left-pointing arrow e.g.
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{b}}}$, and similarly
for matrices in $\mathbb{R}^{\tilde{s}\times\tilde{s}}$.

  * •

Vectors whose dimension is the product of the dimensions listed above are
written with a combination of attributes, e.g.
$\accentset{\text{\tiny$\leftrightarrow$}}{\mathbf{\boldsymbol{x}}}\in\mathbb{R}^{\tilde{s}s}$,
$\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{y}}}}}}\in\mathbb{R}^{sm}$,
$\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftrightarrow$}}{\mathbf{\boldsymbol{z}}}}}}\in\mathbb{R}^{\tilde{s}sm}$.

  * •

Vectors representing quantities on spatial grids are expressed in a bold sans
serif font, e.g.

$\mathbf{\boldsymbol{\mathsf{u}}}^{\top}=(\dots,u_{i-1},u_{i},u_{i+1},\dots).$

This notation is combined with the accents above if the evaluated quantity in
question is a vector. Thus, a vector of $m$-element vectors is represented as

$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{\top}=(\dots,\underline{\mathbf{\boldsymbol{u}}}_{i-1}^{\top},\underline{\mathbf{\boldsymbol{u}}}_{i}^{\top},\underline{\mathbf{\boldsymbol{u}}}_{i+1}^{\top},\dots).$

Matrices operating on these vectors are denoted similarly with capital letters
and are constructed block-diagonally:

$\underline{\mathbf{\boldsymbol{\mathsf{A}}}}=\text{blkdiag}(\dots,\underline{\mathbf{\boldsymbol{A}}}_{i-1},\underline{\mathbf{\boldsymbol{A}}}_{i},\underline{\mathbf{\boldsymbol{A}}}_{i+1},\dots).$

  * •

A flux function that takes $(p+q+1)$ arguments, e.g.
$\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{u}}}_{i-p},\dots,\underline{\mathbf{\boldsymbol{u}}}_{i+q})$,
is sometimes denoted with the abbreviated argument
$\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}})$.

###  2.2 Conservation laws and the Lax-Wendroff theorem

Consider the system of 1D conservation laws

$\underline{\mathbf{\boldsymbol{u}}}_{t}+\underline{\mathbf{\boldsymbol{f}}}_{x}=\underline{\mathbf{\boldsymbol{0}}},\quad\underline{\mathbf{\boldsymbol{u}}}(x,0)=\underline{\mathbf{\boldsymbol{u}}}_{0}(x),\quad
x\in\Omega,\quad t>0.$ (1)

Here,
$\underline{\mathbf{\boldsymbol{u}}}(x,t),\,\underline{\mathbf{\boldsymbol{f}}},\,\underline{\mathbf{\boldsymbol{u}}}_{0}\in\mathbb{R}^{m}$.
Throughout, it is assumed that (1) is posed either as a Cauchy problem or on a
periodic domain. Under these circumstances, the quantity
$\int_{\Omega}\underline{\mathbf{\boldsymbol{u}}}\text{d}x$ is conserved. A
space-time discretization of (1) that discretely mimics this property is said
to be _globally conservative_.

We consider discretizations of (1) that may be expressed in the form

$\frac{\underline{\mathbf{\boldsymbol{u}}}_{i}^{n+1}-\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}}{\Delta
t}+\frac{1}{\Delta
x}\left(\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}-\underline{\mathbf{\boldsymbol{f}}}_{i-\frac{1}{2}}\right)=\underline{\mathbf{\boldsymbol{0}}},\quad
i=\dots,-1,0,1,\dots$ (2)

Here,
$\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}\approx\underline{\mathbf{\boldsymbol{u}}}(x_{i},t_{n})\equiv\underline{\mathbf{\boldsymbol{u}}}(i\Delta
x,n\Delta t)$. In this paper we consider implicit discretizations and
therefore restrict our attention to numerical fluxes of the type

$\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}\equiv\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{n+1})=\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{u}}}_{i-p}^{n+1},\dots,\underline{\mathbf{\boldsymbol{u}}}_{i+q}^{n+1}),$

where $p$ and $q$ are nonnegative integers with $p+q>0$. Throughout, it is
assumed that these numerical fluxes are consistent:

######  Definition 1.

The numerical flux
$\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}})=\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{u}}}_{i-p},\dots,\underline{\mathbf{\boldsymbol{u}}}_{i+q})$
is said to be consistent with
$\underline{\mathbf{\boldsymbol{f}}}(\underline{\mathbf{\boldsymbol{u}}})$ if
it it Lipschitz continuous in each argument and if
$\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{u}}},\dots,\underline{\mathbf{\boldsymbol{u}}})=\underline{\mathbf{\boldsymbol{f}}}(\underline{\mathbf{\boldsymbol{u}}})$.

The concept of local conservation will be central in the remainder. It applies
to explicit and implicit discretizations alike:

######  Definition 2.

A discretization of (1) that can be expressed in the form (2) is said to be
locally conservative.

Local conservation is a useful property for both physical and mathematical
reasons (1). It enforces that solution components leaving one computational
cell necessarily enter the neighbouring one. Conservation of ”total mass”,
$\sum_{i}\Delta x\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}$ is thereby
ensured, in analogy with the continuous problem (i.e. global conservation).
Further, it is an essential ingredient in the ubiquitous Lax-Wendroff theorem.

The Lax-Wendroff theorem applies to the Cauchy problem for (1) and considers
locally conservative discretizations with consistent numerical flux. If the
numerical solution of such a scheme converges to a function
$\underline{\mathbf{\boldsymbol{u}}}$ in the limit of vanishing $\Delta x$ and
$\Delta t$, the theorem provides sufficient conditions for
$\underline{\mathbf{\boldsymbol{u}}}$ to be a weak solution of the
conservation law (1) [13, Chapter 12]. More precisely, consider a sequence of
grids $(\Delta x_{\ell},\Delta t_{\ell})$ such that $\Delta x_{\ell},\Delta
t_{\ell}\rightarrow 0$ as $\ell\rightarrow\infty$. Let
$\underline{\mathbf{\boldsymbol{\mathcal{U}}}}_{\ell}(x,t)$ denote the
piecewise constant function that takes the solution value
$\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}$ in
$(x_{i},x_{i+1}]\times(t_{n-1},t_{n}]$ on the $\ell$th grid. We make the
following assumptions:

######  Assumption 1.

  1. 1.

There is a function $\underline{\mathbf{\boldsymbol{u}}}(x,t)$ such that over
every bounded set $\Omega=[a,b]\times[0,T]$ in $x$-$t$ space,

$\|\underline{\mathbf{\boldsymbol{\mathcal{U}}}}_{\ell}(x,t)-\underline{\mathbf{\boldsymbol{u}}}(x,t)\|_{1,\Omega}\rightarrow
0\quad\text{as}\quad\ell\rightarrow\infty.$

  2. 2.

For each $T\geq 0$ there is a constant $R>0$ such that the total variation

$TV(\underline{\mathbf{\boldsymbol{\mathcal{U}}}}_{\ell}(\cdot,t))<R\quad\text{for
all}\quad 0\leq t\leq T,\quad\ell=1,2,\dots$

The Lax-Wendroff theorem can then be stated as follows:

######  Theorem 1.

Consider a sequence of grids $(\Delta x_{\ell},\Delta t_{\ell})$ such that
$\Delta x_{\ell},\Delta t_{\ell}\rightarrow 0$ as $\ell\rightarrow\infty$.
Suppose that the numerical flux
$\underline{\mathbf{\boldsymbol{f}}}_{i\pm\frac{1}{2}}$ in (2) is consistent
with $\underline{\mathbf{\boldsymbol{f}}}$ and that 1 is satisfied. Then,
$\underline{\mathbf{\boldsymbol{u}}}(x,t)$ is a weak solution of (1).

1 is not strictly speaking necessary for the Lax-Wendroff theorem. Both
conditions can be relaxed somewhat. Indeed, the original proof due to Lax and
Wendroff instead assumes that
$\underline{\mathbf{\boldsymbol{\mathcal{{U}}}}}_{\ell}$ converges boundedly
almost everywhere to $\underline{\mathbf{\boldsymbol{u}}}$ [12].

The discretization (2) is implicit, hence the solution generally must be
approximated using iterative methods. Since we are interested in the
convergence properties of schemes involving iterative methods, the concepts of
local conservation and consistency must be extended to this setting. To this
end, we propose the following definition:

######  Definition 3.

Suppose that the solution to (2) is approximated by a sequence of iterates
$\underline{\mathbf{\boldsymbol{u}}}_{i}^{(k)},\,k=0,\dots,N$ and set
$\underline{\mathbf{\boldsymbol{u}}}^{n+1}=\underline{\mathbf{\boldsymbol{u}}}_{i}^{(N)}$.
If there is a numerical flux function
$\underline{\mathbf{\boldsymbol{h}}}_{i+\frac{1}{2}}^{(N)}$ such that the
approximate numerical solution $\underline{\mathbf{\boldsymbol{u}}}^{n+1}$
satisfies

$\frac{\underline{\mathbf{\boldsymbol{u}}}_{i}^{n+1}-\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}}{\Delta
t}+\frac{1}{\Delta
x}\left(\underline{\mathbf{\boldsymbol{h}}}_{i+\frac{1}{2}}^{(N)}-\underline{\mathbf{\boldsymbol{h}}}_{i-\frac{1}{2}}^{(N)}\right)=\underline{\mathbf{\boldsymbol{0}}},\qquad
i=\dots,-1,0,1,\dots$ (3)

then the iterative method is said to be locally conservative. If
$\underline{\mathbf{\boldsymbol{h}}}_{i+\frac{1}{2}}^{(N)}$ further satisfies
1, then the iterative method is said to be flux consistent.

##  3 Conservation and consistency of pseudo-time iterations

In order to approximate the solution of the nonlinear system (2) using pseudo-
time iterations, we introduce a pseudo-time derivative,

$\frac{\partial\underline{\mathbf{\boldsymbol{u}}}_{i}}{\partial\tau}+\underline{\mathbf{\boldsymbol{g}}}_{i}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}})=0,\qquad\underline{\mathbf{\boldsymbol{u}}}_{i}(0)=\underline{\mathbf{\boldsymbol{u}}}_{0_{i}},\qquad
i=\dots,-1,0,1,\dots$

where the nonlinear function $\underline{\mathbf{\boldsymbol{g}}}_{i}$ is
given by

$\underline{\mathbf{\boldsymbol{g}}}_{i}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}})=\frac{\underline{\mathbf{\boldsymbol{u}}}_{i}-\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}}{\Delta
t}+\frac{1}{\Delta
x}\left(\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}})-\underline{\mathbf{\boldsymbol{f}}}_{i-\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}})\right).$
(4)

Several different methods are available for iterating in pseudo-time [19, 3].
Herein, we use an $s$-stage explicit Runge-Kutta (ERK) method. Let
$(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}},\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}},\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{c}}})$
denote the coefficient matrix and vectors of the ERK method. We denote the
$k$th pseudo-time iterate by $\underline{\mathbf{\boldsymbol{u}}}_{i}^{(k)}$.
The subsequent iterate $\underline{\mathbf{\boldsymbol{u}}}_{i}^{(k+1)}$ is
computed from $\underline{\mathbf{\boldsymbol{u}}}_{i}^{(k)}$ as

$\underline{\mathbf{\boldsymbol{u}}}_{i}^{(k+1)}=\underline{\mathbf{\boldsymbol{u}}}_{i}^{(k)}-\Delta\tau_{k}\sum_{j=1}^{s}b_{j}\underline{\mathbf{\boldsymbol{g}}}_{i}\left(\underline{\mathbf{\boldsymbol{\mathsf{U}}}}_{j}^{(k)}\right),\qquad
i=\dots,-1,0,1,\dots$ (5)

where the stage vectors
$\underline{\mathbf{\boldsymbol{\mathsf{U}}}}_{j}^{(k)}$, $j=1,\dots,s$ have
elements

$\underline{\mathbf{\boldsymbol{U}}}_{j_{\iota}}^{(k)}=\underline{\mathbf{\boldsymbol{u}}}_{\iota}^{(k)}-\Delta\tau_{k}\sum_{l=1}^{j-1}a_{j,l}\underline{\mathbf{\boldsymbol{g}}}_{\iota}\left(\underline{\mathbf{\boldsymbol{\mathsf{U}}}}_{l}^{(k)}\right),\quad\iota=i-p,\dots,i+q.$
(6)

As previously, $p$ and $q$ determine the bandwidth of the finite volume
stencil.

###  3.1 Systems of conservation laws

In [4] the scalar version (i.e. $m=1$) of (1) was considered. There, it was
shown that the scheme (5)–(6) preserves the local conservation of the space-
time discretization. Further, an extension of the Lax-Wendroff theorem was
provided that incorporates a fixed number $N$ of pseudo-time iterations.

Here, we present two theorems that generalize the results in [4] to systems of
conservation laws (i.e. $m\geq 1$). We define a step in physical time by
setting
$\underline{\mathbf{\boldsymbol{u}}}_{i}^{n+1}=\underline{\mathbf{\boldsymbol{u}}}_{i}^{(N)}$.
Throughout,
$\underline{\mathbf{\boldsymbol{u}}}_{i}^{(0)}=\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}$
is chosen as initial guess.

Recall that the _stability function_ $\phi(z)$ of an RK method
$(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}},\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}},\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{c}}})$
is given by

$\phi(z)=1+z\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}}^{\top}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}-z\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}})^{-1}\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}},$
(7)

where $\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}$ is the
$s\times s$ identity matrix and
$\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\in\mathbb{R}^{s}$
is the vector of all ones; see e.g. [21, Chapter IV.3]. The _stability region_
of the RK method is defined as the subset of the complex plane for which
$|\phi(z)|<1$.

The proofs of the following theorems are very similar to those presented for
the scalar case in [4]. The details are therefore omitted. The following is a
generalization to systems of conservation laws, which also allows each pseudo-
time step to be taken with different ERK methods.

######  Theorem 2.

Choose the initial guess
$\underline{\mathbf{\boldsymbol{u}}}_{i}^{(0)}=\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}$.
Apply $N$ pseudo-time iterations to (2), where the $k$th iteration is
performed with an ERK method
$(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}_{k},\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}}_{k},\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{c}}}_{k})$
with stability function $\phi_{k}(z)$ and pseudo-time step $\Delta\tau_{k}$.
Let $\mu_{k}=\Delta\tau_{k}/\Delta t$ for $k=0,\dots,N-1$. The pseudo-time
iterations are locally conservative with numerical flux

$\underline{\mathbf{\boldsymbol{h}}}_{i+\frac{1}{2}}^{(N)}=\sum_{k=0}^{N-1}\left(\mu_{k}\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}}_{k}^{\top}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}+\mu_{k}\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}_{k})^{-1}\otimes\underline{\mathbf{\boldsymbol{I}}}\right)\left(\prod_{l=k+1}^{N-1}\phi_{l}(-\mu_{l})\right)\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{f}}}}}}_{i+\frac{1}{2}}^{(k)}.$
(8)

Here, $\underline{\mathbf{\boldsymbol{I}}}$ is the $m\times m$ identity matrix
and

$\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{f}}}}}}_{i+\frac{1}{2}}^{(k)}=\left(\underline{\mathbf{\boldsymbol{f}}}^{\top}_{i+\frac{1}{2}}\left(\underline{\mathbf{\boldsymbol{\mathsf{U}}}}_{1}^{(k)}\right),\dots,\underline{\mathbf{\boldsymbol{f}}}^{\top}_{i+\frac{1}{2}}\left(\underline{\mathbf{\boldsymbol{\mathsf{U}}}}_{s}^{(k)}\right)\right)^{\top}\in\mathbb{R}^{sm}.$

The numerical flux is consistent with
$c\underline{\mathbf{\boldsymbol{f}}}(u)$, where

$c\equiv c(\mu_{0},\dots,\mu_{N-1})=1-\prod_{l=0}^{N-1}\phi_{l}(-\mu_{l}).$
(9)

Thus, pseudo-time iterations are flux consistent if and only if $c=1$.

######  Remark 1.

The product in (8) is empty when $k=N-1$. To handle this case we use the
convention

$\prod_{l=N}^{N-1}\phi_{l}(-\mu_{l})=1.$

###### Proof.

The proof is step by step the same as those of [4, Lemma 3 & Theorem 2]. The
only changes necessary are to replace
$\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}$ by
$\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{I}}}$,
and to swap each multiplication of the form
$\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\phi$, where
$\phi$ is a scalar, to a corresponding Kronecker product
$\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\phi}}}$.
∎

2 reveals that pseudo-time iterations are locally conservative but not
necessarily flux consistent, as characterized by the scalar $c$ in (9). To
remove the inconsistency, it suffices to make a single iteration with any
explicit RK method for which $\phi_{l}(-\mu_{l})=0$. In [4] it was suggested
to iterate once with the explicit Euler method, choosing
$\Delta\tau_{0}=\Delta t$ so that $\mu_{0}=1$. The stability function is
$\phi(-\mu_{0})=1-\mu_{0}$, hence $\phi(-1)=0$. This rids the iterative method
of its flux inconsistency and the remaining iterations can be made with any
other ERK method as preferred. Numerical experiments in [4] showed that
enforcing flux consistency can have a profound impact on the convergence of
the pseudo-time iterations.

The coefficient $c$ causes the numerical fluxes in (8) to be consistent with
the modified system of conservation laws

$\underline{\mathbf{\boldsymbol{u}}}_{t}+c(\mu_{0},\dots,\mu_{N-1})\underline{\mathbf{\boldsymbol{f}}}_{x}=0.$
(10)

Note that $c$ affects all components of the system identically. An
interpretation of its presence is that the pseudo-time iterations alter the
rate of flow of time. Defining the modified time $t_{c}=ct$ it follows that
$\underline{\mathbf{\boldsymbol{u}}}_{t}=c\underline{\mathbf{\boldsymbol{u}}}_{t_{c}}$
so that
$\underline{\mathbf{\boldsymbol{u}}}_{t_{c}}+\underline{\mathbf{\boldsymbol{f}}}_{x}=\underline{\mathbf{\boldsymbol{0}}}$.
Hence, $t_{c}$ rather than $t$ is the governing time variable in (10).

Convergent pseudo-time iterations require that the pseudo-time steps are
restricted to the stability domain of the explicit RK method, at least for
most of the iterations. Thus, the product in (9) is generally taken over
factors bounded by unity, and consequently $c\leq 1$. Thus, $t_{c}$ represents
a time retardation, i.e. time flows slower in (10) than in the original
conservation law (1). This is in line with the experimental observations made
in [4].

###  3.2 Higher order implicit Runge-Kutta methods

Let us return to the system of conservation laws (1) and introduce a spatial
semi-discretization

$\underline{\mathbf{\boldsymbol{u}}}_{t}+\frac{1}{\Delta
x}\left(\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}})-\underline{\mathbf{\boldsymbol{f}}}_{i-\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}})\right)=\underline{\mathbf{\boldsymbol{0}}}.$
(11)

If the implicit Euler method is used to discretize (11) in time, then (2) is
recovered. However, suppose that we instead wish to discretize in time using
an $\tilde{s}$-stage implicit Runge-Kutta (IRK) method with Butcher matrix and
vectors
$(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{A}}},\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{b}}},\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{c}}})$.
The resulting scheme can be expressed as

$\displaystyle\underline{\mathbf{\boldsymbol{U}}}_{i}^{(j)}$
$\displaystyle=\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}-\frac{\Delta
t}{\Delta
x}\sum_{l=1}^{\tilde{s}}\accentset{\text{\tiny$\leftarrow$}}{a}_{jl}\left(\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{U}}}}^{(j)})-\underline{\mathbf{\boldsymbol{f}}}_{i-\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{U}}}}^{(j)})\right),\quad
j=1,\dots,\tilde{s},$ (12)
$\displaystyle\underline{\mathbf{\boldsymbol{u}}}_{i}^{n+1}$
$\displaystyle=\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}-\frac{\Delta
t}{\Delta
x}\sum_{j=1}^{\tilde{s}}\accentset{\text{\tiny$\leftarrow$}}{b}_{j}\left(\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{U}}}}^{(j)})-\underline{\mathbf{\boldsymbol{f}}}_{i-\frac{1}{2}}(\underline{\mathbf{\boldsymbol{\mathsf{U}}}}^{(j)})\right).$

Define the quantities

$\displaystyle\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{U}}}}}}_{i}$
$\displaystyle=\left(\underline{\mathbf{\boldsymbol{U}}}_{i}^{(1),\top},\dots,\underline{\mathbf{\boldsymbol{U}}}_{i}^{(\tilde{s}),\top}\right)^{\top}\in\mathbb{R}^{\tilde{s}m},$
$\displaystyle\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{f}}}}}}_{i+\frac{1}{2}}(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})$
$\displaystyle=\left(\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}^{\top}(\underline{\mathbf{\boldsymbol{\mathsf{U}}}}^{(1)}),\dots,\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}^{\top}(\underline{\mathbf{\boldsymbol{\mathsf{U}}}}^{(\tilde{s})})\right)^{\top}\in\mathbb{R}^{\tilde{s}m}.$

Then, (12) can be reformulated as

$\displaystyle\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{U}}}}}}_{i}$
$\displaystyle=\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}-\frac{\Delta
t}{\Delta
x}\left(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{I}}}\right)\left(\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{f}}}}}}_{i+\frac{1}{2}}(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})-\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{f}}}}}}_{i-\frac{1}{2}}(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})\right),$
(13) $\displaystyle\underline{\mathbf{\boldsymbol{u}}}_{i}^{n+1}$
$\displaystyle=\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}-\frac{\Delta
t}{\Delta
x}\left(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{b}}}^{\top}\otimes\underline{\mathbf{\boldsymbol{I}}}\right)\left(\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{f}}}}}}_{i+\frac{1}{2}}(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})-\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{f}}}}}}_{i-\frac{1}{2}}(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})\right),$

where $\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{1}}}$ is the
vector of ones in $\mathbb{R}^{\tilde{s}}$.

We seek an approximation of the solution
$\underline{\mathbf{\boldsymbol{u}}}_{i}^{n+1}$. To this end we must find an
approximate solution of the nonlinear equation system in the second line of
(13). Rewriting this equation as

$\frac{\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{U}}}}}}_{i}-\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}}{\Delta
t}+\frac{1}{\Delta
x}\left(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{I}}}\right)\left(\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{f}}}}}}_{i+\frac{1}{2}}(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})-\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{f}}}}}}_{i-\frac{1}{2}}(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})\right)=\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{0}}}}}},$
(14)

we see that it is of precisely the same form as the discretization (2),
although with
$\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{U}}}}}}_{i}$
taking the place of $\underline{\mathbf{\boldsymbol{u}}}_{i}^{n+1}$,
$\accentset{\text{\tiny$\leftarrow$}}{1}\otimes\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}$
replacing $\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}$ and
$(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{I}}})\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{f}}}}}}_{i+\frac{1}{2}}$
in place of $\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}$. Thus, if
pseudo-time iterations are used to approximate a solution, then 2 applies and
we can immediately conclude that the resulting scheme can be written in the
conservative form

$\frac{\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{U}}}}}}_{i}-\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}}{\Delta
t}+\frac{1}{\Delta
x}\left(\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{h}}}}}}_{i+\frac{1}{2}}^{(N)}-\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{h}}}}}}_{i-\frac{1}{2}}^{(N)}\right)=\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{0}}}}}},$
(15)

where the numerical flux is given by

$\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{h}}}}}}_{i+\frac{1}{2}}^{(N)}=\sum_{k=0}^{N-1}\left(\mu_{k}\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}}^{\top}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}+\mu_{k}\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}})\otimes\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{I}}}\right)\left(\prod_{l=k+1}^{N-1}\phi_{l}(\mu_{l})\right)\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftrightarrow$}}{f}}}}_{i+\frac{1}{2}}^{(k)},$

and

$\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftrightarrow$}}{f}}}}_{i+\frac{1}{2}}^{(k)}=\left(\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{f}}}}^{\top}\left(\underline{\mathbf{\boldsymbol{\mathsf{\accentset{\text{\tiny$\leftarrow$}}{U}}}}}_{1}^{(k)}\right),\dots,\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{f}}}}^{\top}\left(\underline{\mathbf{\boldsymbol{\mathsf{\accentset{\text{\tiny$\leftarrow$}}{U}}}}}_{s}^{(k)}\right)\right)^{\top}.$

At this point we note that if the scheme is such that a vector
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}$ exists,
satisfying

$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}^{\top}\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{A}}}=\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{b}}}^{\top},\qquad\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}^{\top}\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{1}}}=1,$
(16)

then the second line in (13) can be evaluated by left-multiplying the first
line by
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}^{\top}\otimes\underline{\mathbf{\boldsymbol{I}}}$.
In other words, it follows that

$\underline{\mathbf{\boldsymbol{u}}}_{i}^{n+1}=\left(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}^{\top}\otimes\underline{\mathbf{\boldsymbol{I}}}\right)\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{U}}}}}}_{i}.$
(17)

If we adopt this principle and use it to compute
$\underline{\mathbf{\boldsymbol{u}}}_{i}^{n+1}$ based on the approximation of
$\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{U}}}}}}_{i}$
obtained by applying pseudo-time iterations to (14), then the resulting scheme
is conservative:

######  Theorem 3.

Apply $N$ pseudo-time iterations to the stage equations (14) with the same
assumptions as in 2. Suppose that a vector
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}$ exists that
satisfies conditions (16) and compute
$\underline{\mathbf{\boldsymbol{u}}}_{i}^{n+1}$ using (17). Then the pseudo-
time iterations are locally conservative with the numerical flux

$\underline{\mathbf{\boldsymbol{h}}}_{i+\frac{1}{2}}^{(N)}=\sum_{k=0}^{N-1}\left(\mu_{k}\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}+\mu_{k}\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}})\otimes\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{b}}}^{\top}\otimes\underline{\mathbf{\boldsymbol{I}}}\right)\left(\prod_{l=k+1}^{N-1}\phi(\mu_{l})\right)\underline{\mathbf{\boldsymbol{\accentset{\text{\tiny$\leftrightarrow$}}{f}}}}_{i+\frac{1}{2}}^{(k)}.$
(18)

The numerical flux is consistent with
$c(\mu_{0},\dots,\mu_{N-1})\underline{\mathbf{\boldsymbol{f}}}$. Thus, the
pseudo-time iterations are flux consistent if and only if $c=1$.

###### Proof.

Local conservation with the flux (18) follows from left-multiplying (14) by
$\left(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}^{\top}\otimes\underline{\mathbf{\boldsymbol{I}}}\right)$
and using (16) to conclude that
$\left(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}^{\top}\otimes\underline{\mathbf{\boldsymbol{I}}}\right)(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{u}}}_{i}^{n})=\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}$.
Consistency with $c\underline{\mathbf{\boldsymbol{f}}}$ follows from the fact
that
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{b}}}^{\top}\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{1}}}=1$
for every consistent RK method. ∎

It should be noted that we cannot find a vector
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}$ that satisfies
conditions (16) for all IRK methods. For example, the implicit midpoint rule
is given by
$(\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{A}}},\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{b}}},\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{c}}})=(1,1/2,1/2)$.
Thus, the first condition in (16) gives
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}=1/2$ whereas
the second one gives
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}=1$, both of
which cannot be satisfied. However, we remark that the important class of IRK
methods associated with Summation-By-Parts (SBP) methods all have such a
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}$ by
construction [5, 14]. This is a broad class of methods with the ability to
preserve $L^{2}$-type estimates of the solution to systems of differential
equations [15] and encompass the ubiquitous Radau IA and IIA and Lobatto IIIC
methods [16] as special cases. For further details about the theoretical and
practical aspects of SBP methods for time marching, see [20] and the
references therein.

A generalization of the Lax-Wendroff theorem can be obtained if we restrict
ourselves to considering a fixed number of iterations on a sequence of ever
finer grids:

######  Theorem 4.

Consider a sequence of grids $(\Delta x_{\ell},\Delta t_{\ell})$ such that
$\Delta x_{\ell},\Delta t_{\ell}\rightarrow 0$ as $\ell\rightarrow\infty$. Fix
$N$ independently of $\ell$, set
$\underline{\mathbf{\boldsymbol{u}}}_{i}^{(0)}=\underline{\mathbf{\boldsymbol{u}}}_{i}^{n}$.
Apply $N$ pseudo-time iterations to the conservative discretization (2), or to
the stage equations (6) followed by the RK step (5). Let
$\Delta\tau_{k,\ell}/\Delta t_{\ell}=\mu_{k,\ell}=\mu_{k}$ be constants
independent of $\ell$ for each $k=0,\dots,N-1$. Suppose that the numerical
flux $\underline{\mathbf{\boldsymbol{f}}}_{i\pm\frac{1}{2}}$ in (2) is
consistent with $\underline{\mathbf{\boldsymbol{f}}}$ and that 1 is satisfied.
Then, $\underline{\mathbf{\boldsymbol{u}}}(x,t)$ is a weak solution of the
conservation law (10).

###### Proof.

The proof is identical to that of [4, Theorem 3], with the same changes as
those in the proof of 2. ∎

##  4 Newton’s method

We now return to the nonlinear system (2), or equivalently to the implicit
Runge-Kutta stage equations in (13) and consider Newton’s method, which
replaces the nonlinear system with a sequence of linear ones. The solution to
each linear system can then be approximated using pseudo-time iterations, or
more commonly,using a Krylov subspace method.

In this section we limit our attention to the case when the numerical flux is
bivariate, i.e. when
$\underline{\mathbf{\boldsymbol{f}}}_{i+\frac{1}{2}}=\hat{\underline{\mathbf{\boldsymbol{f}}}}(\underline{\mathbf{\boldsymbol{u}}}_{i},\underline{\mathbf{\boldsymbol{u}}}_{i+1})$
for some function
$\hat{\underline{\mathbf{\boldsymbol{f}}}}(\underline{\mathbf{\boldsymbol{\theta}}},\underline{\mathbf{\boldsymbol{\phi}}})$.
It will be convenient to split the flux into a symmetric (or _convective_) and
an anti-symmetric (or _dissipative_) component;

$\hat{\underline{\mathbf{\boldsymbol{f}}}}=\hat{\underline{\mathbf{\boldsymbol{f}}}}^{(+)}+\hat{\underline{\mathbf{\boldsymbol{f}}}}^{(-)},\qquad\hat{\underline{\mathbf{\boldsymbol{f}}}}^{(+)}(\underline{\mathbf{\boldsymbol{\theta}}},\underline{\mathbf{\boldsymbol{\phi}}})=\hat{\underline{\mathbf{\boldsymbol{f}}}}^{(+)}(\underline{\mathbf{\boldsymbol{\phi}}},\underline{\mathbf{\boldsymbol{\theta}}}),\qquad\hat{\underline{\mathbf{\boldsymbol{f}}}}^{(-)}(\underline{\mathbf{\boldsymbol{\theta}}},\underline{\mathbf{\boldsymbol{\phi}}})=-\hat{\underline{\mathbf{\boldsymbol{f}}}}^{(-)}(\underline{\mathbf{\boldsymbol{\phi}}},\underline{\mathbf{\boldsymbol{\theta}}}).$

Any bivariate function can be expressed in this way by setting

$\hat{\underline{\mathbf{\boldsymbol{f}}}}^{(+)}(\underline{\mathbf{\boldsymbol{\theta}}},\underline{\mathbf{\boldsymbol{\phi}}})=\frac{1}{2}(\hat{\underline{\mathbf{\boldsymbol{f}}}}(\underline{\mathbf{\boldsymbol{\theta}}},\underline{\mathbf{\boldsymbol{\phi}}})+\hat{\underline{\mathbf{\boldsymbol{f}}}}(\underline{\mathbf{\boldsymbol{\phi}}},\underline{\mathbf{\boldsymbol{\theta}}})),\qquad\hat{\underline{\mathbf{\boldsymbol{f}}}}^{(-)}(\underline{\mathbf{\boldsymbol{\theta}}},\underline{\mathbf{\boldsymbol{\phi}}})=\frac{1}{2}(\hat{\underline{\mathbf{\boldsymbol{f}}}}(\underline{\mathbf{\boldsymbol{\theta}}},\underline{\mathbf{\boldsymbol{\phi}}})-\hat{\underline{\mathbf{\boldsymbol{f}}}}(\underline{\mathbf{\boldsymbol{\phi}}},\underline{\mathbf{\boldsymbol{\theta}}})).$

By anti-symmetry, the dissipative component satisfies
$\hat{\underline{\mathbf{\boldsymbol{f}}}}^{(-)}(\underline{\mathbf{\boldsymbol{u}}},\underline{\mathbf{\boldsymbol{u}}})=\underline{\mathbf{\boldsymbol{0}}}$.
Consistency of the numerical flux is therefore equivalent to consistency of
the convective component.

In the following subsections, we demonstrate that Newton’s method is locally
conservative and flux consistent when bivariate fluxes are used. To simplify
the presentation, we begin by proving these results for scalar conservation
laws. The extension to systems is straightforward but notationally
complicated. We therefore postpone this to a separate subsection.

###  4.1 Scalar conservation laws

Consider the scalar conservation law

$u_{t}+f_{x}=0,$

posed on a periodic spatial domain and adjoined with appropriate initial data.
Let $\hat{f}^{(\pm)}(\theta,\phi)$ be a consistent bivariate numerical flux.
Without loss of generality we may assume that the flux is either symmetric or
anti-symmetric. The analysis will proceed in the same way in both cases and we
may thereafter form linear combinations of such fluxes as we like. We
discretize with a finite volume method,

$\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta
t}+\frac{\hat{f}^{(\pm)}(u_{i}^{n+1},u_{i+1}^{n+1})-\hat{f}^{(\pm)}(u_{i-1}^{n+1},u_{i}^{n+1})}{\Delta
x}=0.$ (19)

The discretization (19) may equivalently be expressed in vector form as

$\frac{\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}}{\Delta
t}+\frac{1}{\Delta
x}(\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}\circ\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}(\mathbf{\boldsymbol{\mathsf{u}}}^{n+1}))\mathbf{\boldsymbol{\mathsf{1}}}=\mathbf{\boldsymbol{\mathsf{0}}}.$
(20)

The dimensions of the vectors and matrices match the number of cells in the
computational grid. Here, $\circ$ denotes the Hadamard product. The elements
of the matrix
$\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}(\mathbf{\boldsymbol{\mathsf{u}}})$
are given by
$(\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)})_{ij}=\hat{f}^{(\pm)}(u_{i},u_{j})$
and $\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}$ is given by

$\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}=\begin{bmatrix}0&1&&\dots&\mp 1\\\
\mp 1&0&1&&\\\ &\mp 1&0&1&\\\ \vdots&&&\ddots&1\\\ 1&&0&\mp 1&0\end{bmatrix}.$

We define the function
$\mathbf{\boldsymbol{\mathsf{g}}}(\mathbf{\boldsymbol{\mathsf{v}}})$ as

$\mathbf{\boldsymbol{\mathsf{g}}}(\mathbf{\boldsymbol{\mathsf{v}}}):=\frac{\mathbf{\boldsymbol{\mathsf{v}}}-\mathbf{\boldsymbol{\mathsf{u}}}^{n}}{\Delta
t}+\frac{1}{\Delta
x}(\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}\circ\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}(\mathbf{\boldsymbol{\mathsf{v}}}))\mathbf{\boldsymbol{\mathsf{1}}}.$
(21)

Newton’s method applied to the nonlinear system (20) is then given by

$\mathbf{\boldsymbol{\mathsf{g}}}^{\prime}(\mathbf{\boldsymbol{\mathsf{v}}}^{(k)})\Delta\mathbf{\boldsymbol{\mathsf{v}}}+\mathbf{\boldsymbol{\mathsf{g}}}(\mathbf{\boldsymbol{\mathsf{v}}}^{(k)})=\mathbf{\boldsymbol{\mathsf{0}}},\qquad\mathbf{\boldsymbol{\mathsf{v}}}^{(k+1)}=\mathbf{\boldsymbol{\mathsf{v}}}^{(k)}+\Delta\mathbf{\boldsymbol{\mathsf{v}}},$
(22)

where $\mathbf{\boldsymbol{\mathsf{g}}}^{\prime}$ is the Jacobian of
$\mathbf{\boldsymbol{\mathsf{g}}}$. An explicit expression for
$\mathbf{\boldsymbol{\mathsf{g}}}^{\prime}$ is given in [6, Theorems 2.1 &
5.1]. Let $\hat{f}^{(\pm)}_{\phi}=\partial\hat{f}^{(\pm)}/\partial\phi$ and
introduce the matrix
$\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{\phi}(\mathbf{\boldsymbol{\mathsf{v}}})$
with elements
$(\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{\phi})_{ij}=\hat{f}^{(\pm)}_{\phi}(v_{i},v_{j})$.
Defining the matrix

$\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}):=\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}\circ\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{\phi}(\mathbf{\boldsymbol{\mathsf{v}}})-\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{1}}}^{\top}(\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}\circ\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{\phi}(\mathbf{\boldsymbol{\mathsf{v}}}))),$

the Jacobian is given by
$\mathbf{\boldsymbol{\mathsf{g}}}^{\prime}(\mathbf{\boldsymbol{\mathsf{v}}}^{(k)})=\Delta
t^{-1}\mathbf{\boldsymbol{\mathsf{I}}}+\Delta
x^{-1}\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}(\mathbf{\boldsymbol{\mathsf{v}}}^{(k)})$.
This is a tridiagonal matrix, which we may explicitly write as

$\displaystyle\mathbf{\boldsymbol{\mathsf{g}}}^{\prime}(\mathbf{\boldsymbol{\mathsf{v}}}^{(k)})=$
$\displaystyle\text{tri}\left(\mp\frac{\hat{f}^{(\pm)}_{\phi}(v_{i}^{(k)},v_{i-1}^{(k)})}{\Delta
x},\frac{1}{\Delta
t}-\frac{\hat{f}^{(\pm)}_{\phi}(v_{i-1}^{(k)},v_{i}^{(k)})\mp\hat{f}^{(\pm)}_{\phi}(v_{i+1}^{(k)},v_{i}^{(k)})}{\Delta
x},\frac{\hat{f}^{(\pm)}_{\phi}(v_{i}^{(k)},v_{i+1}^{(k)})}{\Delta x}\right)$

Inserting this expression for
$\mathbf{\boldsymbol{\mathsf{g}}}^{\prime}(\mathbf{\boldsymbol{\mathsf{v}}}^{(k)})$
into (22) and collecting terms leads to a system of equations of the form

$\frac{v_{i}^{(k+1)}-u_{i}^{n}}{\Delta t}+\frac{1}{\Delta
x}\left(h_{i+\frac{1}{2}}^{(k+1)}-h_{i-\frac{1}{2}}^{(k+1)}\right)=0,\qquad
i=\dots,-1,0,1,\dots,$ (23)

where the numerical flux function $h_{i+\frac{1}{2}}^{(k+1)}$ is given by

$h_{i+\frac{1}{2}}^{(k+1)}=\hat{f}_{\phi}^{(\pm)}(v_{i}^{(k)},v_{i+1}^{(k)})\Delta
v_{i+1}\pm\hat{f}_{\phi}^{(\pm)}(v_{i+1}^{(k)},v_{i}^{(k)})\Delta
v_{i}+\hat{f}^{(\pm)}(v_{i+1}^{(k)},v_{i}^{(k)}).$ (24)

######  Theorem 5.

Choose $v_{i}^{(0)}=u_{i}^{n}$ and suppose that
$\mathbf{\boldsymbol{\mathsf{g}}}^{\prime}(u\mathbf{\boldsymbol{\mathsf{1}}})$
is nonsingular for any non-zero scalar $u$. Then Newton’s method, applied to
the discretization (20), is locally conservative and flux consistent.

###### Proof.

Setting $u_{i}^{n+1}=v_{i}^{(k+1)}$ in (23) shows that Newton’s method applied
to (20) is locally conservative. To establish flux consistency, let
$\mathbf{\boldsymbol{\mathsf{v}}}^{(k)}$ and
$\mathbf{\boldsymbol{\mathsf{u}}}^{n}$ both be given as
$u\mathbf{\boldsymbol{\mathsf{1}}}$ for some nonzero scalar $u$ and some $k$
(e.g. $k=0$). Observe that, by the consistency of $\hat{f}^{(\pm)}$, we have
$(\mathbf{\boldsymbol{\mathsf{F}}}^{(+)}(u\mathbf{\boldsymbol{\mathsf{1}}}))_{ij}=f(u)$
and
$(\mathbf{\boldsymbol{\mathsf{F}}}^{(-)}(u\mathbf{\boldsymbol{\mathsf{1}}}))_{ij}=0$.
Consequently,
$\mathbf{\boldsymbol{\mathsf{Q}}}^{(+)}\circ\mathbf{\boldsymbol{\mathsf{F}}}^{(-)}(u\mathbf{\boldsymbol{\mathsf{1}}})=\mathbf{\boldsymbol{\mathsf{0}}}$
and
$\mathbf{\boldsymbol{\mathsf{Q}}}^{(-)}\circ\mathbf{\boldsymbol{\mathsf{F}}}^{(+)}(u\mathbf{\boldsymbol{\mathsf{1}}})=f(u)\mathbf{\boldsymbol{\mathsf{Q}}}^{(-)}$.
Since
$\mathbf{\boldsymbol{\mathsf{Q}}}^{(-)}\mathbf{\boldsymbol{\mathsf{1}}}=\mathbf{\boldsymbol{\mathsf{0}}}$
it therefore follows from (21) that
$\mathbf{\boldsymbol{\mathsf{g}}}(u\mathbf{\boldsymbol{\mathsf{1}}})=\mathbf{\boldsymbol{\mathsf{0}}}$.
Since
$\mathbf{\boldsymbol{\mathsf{g}}}^{\prime}(u\mathbf{\boldsymbol{\mathsf{1}}})$
is nonsingular by assumption, (22) implies that
$\Delta\mathbf{\boldsymbol{\mathsf{v}}}=\mathbf{\boldsymbol{\mathsf{0}}}$.
Thus, from (24) we have

$h_{i+\frac{1}{2}}^{(k+1)}=\hat{f}^{(\pm)}_{\phi}(u,u)\cdot
0\pm\hat{f}^{(\pm)}_{\phi}(u,u)\cdot
0+\hat{f}^{(\pm)}(u,u)=\hat{f}^{(\pm)}(u,u),$

the latter of which is consistent by assumption. Note that
$h_{i+\frac{1}{2}}^{(k+1)}$ exclusively depends on the current iterate
$v_{\,}^{(k+1)}$ via $\Delta v$ and is therefore a linear function. The
numerical flux is thus Lipschitz continuous, hence consistent.

Finally, since any bivariate function can be written as a linear combination
of symmetric and anti-symmetric components, the argument above extends by
linearity to arbitrary bivariate numerical fluxes. ∎

###  4.2 Systems of conservation laws

We now return to the system of conservation laws (1) and the discretization
(2), which may correspond either to the implicit Euler method or to the stage
equations of a higher order RK method. The analysis proceeds very similarly to
the scalar case, although we are now dealing with a flux function
$\underline{\mathbf{\boldsymbol{f}}}:\mathbb{R}^{m}\mapsto\mathbb{R}^{m}$.
Suppose that the corresponding numerical flux is of the form

$\hat{\underline{\mathbf{\boldsymbol{f}}}}^{(\pm)}(\underline{\mathbf{\boldsymbol{\theta}}},\underline{\mathbf{\boldsymbol{\phi}}})=(\hat{f}^{(\pm)}_{1}(\underline{\mathbf{\boldsymbol{\theta}}},\underline{\mathbf{\boldsymbol{\phi}}}),\dots,\hat{f}^{(\pm)}_{m}(\underline{\mathbf{\boldsymbol{\theta}}},\underline{\mathbf{\boldsymbol{\phi}}}))^{\top},$

where each $\hat{f}^{(\pm)}_{\ell}$, $\ell=1,\dots,m$, is bivariate,
consistent and either symmetric or anti-symmetric.

This time, Newton’s method is given by

$\underline{\mathbf{\boldsymbol{\mathsf{g}}}}^{\prime}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})\Delta\underline{\mathbf{\boldsymbol{\mathsf{v}}}}+\underline{\mathbf{\boldsymbol{\mathsf{g}}}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})=\underline{\mathbf{\boldsymbol{\mathsf{0}}}},\qquad\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k+1)}=\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)}+\Delta\underline{\mathbf{\boldsymbol{\mathsf{v}}}},$
(25)

where we define
$\underline{\mathbf{\boldsymbol{\mathsf{g}}}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}})$
as

$\displaystyle\underline{\mathbf{\boldsymbol{\mathsf{g}}}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}):$
$\displaystyle=\frac{\underline{\mathbf{\boldsymbol{\mathsf{v}}}}-\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{n}}{\Delta
t}+\frac{1}{\Delta
x}((\underline{\mathbf{\boldsymbol{I}}}\otimes\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)})\circ\underline{\mathbf{\boldsymbol{\mathsf{F}}}}^{(\pm)}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}))\underline{\mathbf{\boldsymbol{\mathsf{1}}}}$
(26)
$\displaystyle=\frac{\underline{\mathbf{\boldsymbol{\mathsf{v}}}}-\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{n}}{\Delta
t}+\frac{1}{\Delta
x}\begin{bmatrix}(\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}\circ\mathbf{\boldsymbol{\mathsf{F}}}_{1})\mathbf{\boldsymbol{\mathsf{1}}}\\\
\vdots\\\
(\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}\circ\mathbf{\boldsymbol{\mathsf{F}}}_{m})\mathbf{\boldsymbol{\mathsf{1}}}\end{bmatrix}.$

Here we have introduced the matrix
$\underline{\mathbf{\boldsymbol{\mathsf{F}}}}^{(\pm)}=\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{1}\oplus\dots\oplus\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{m}$
with
$(\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{\ell})_{ij}=\hat{f}^{(\pm)}_{\ell}(\underline{\mathbf{\boldsymbol{v}}}_{i},\underline{\mathbf{\boldsymbol{v}}}_{j})$.
We denote the partial derivatives of the numerical fluxes as
$\hat{f}^{(\pm)}_{\ell,\phi_{\kappa}}=\partial\hat{f}^{(\pm)}_{\ell}/\partial\phi_{\kappa}$
and define the matrices
$\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{\ell,\phi_{\kappa}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}})$
with elements
$(\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{\ell,\phi_{\kappa}})_{ij}=\hat{f}^{(\pm)}_{\ell,\phi_{\kappa}}(\underline{\mathbf{\boldsymbol{v}}}_{i},\underline{\mathbf{\boldsymbol{v}}}_{j})$.
Finally we obtain

$\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{\ell,\phi_{\kappa}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}):=\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}\circ\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{\ell,\phi_{\kappa}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}})-\mathrm{diag}(\mathbf{\boldsymbol{\mathsf{1}}}^{\top}\left(\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}\circ\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{\ell,\phi_{\kappa}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}})\right)).$

In [6], the Jacobian
$\underline{\mathbf{\boldsymbol{\mathsf{g}}}}^{\prime}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})$
is shown to be given by

$\underline{\mathbf{\boldsymbol{\mathsf{g}}}}^{\prime}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})=\frac{\underline{\mathbf{\boldsymbol{\mathsf{I}}}}}{\Delta
t}+\frac{1}{\Delta
x}\begin{bmatrix}\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{1,\phi_{1}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})&\dots&\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{1,\phi_{m}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})\\\
\vdots&\ddots&\vdots\\\
\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{m,\phi_{1}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})&\dots&\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{m,\phi_{m}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})\end{bmatrix}.$

This is a block matrix formed by $m^{2}$ submatrices, each tridiagonal as in
the scalar case (although the off-diagonal blocks do not depend on $\Delta
t$). If we write
$\Delta\underline{\mathbf{\boldsymbol{\mathsf{v}}}}=(\Delta\mathbf{\boldsymbol{\mathsf{v}}}_{1}^{\top},\dots,\Delta\mathbf{\boldsymbol{\mathsf{v}}}_{m}^{\top})^{\top}$
and collect terms, then the Newton iteration (25) takes the form

$\displaystyle\frac{\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k+1)}-\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{n}}{\Delta
t}$ $\displaystyle+\frac{1}{\Delta
x}\begin{bmatrix}\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{1,\phi_{1}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})&\dots&\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{1,\phi_{m}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})\\\
\vdots&\ddots&\vdots\\\
\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{m,\phi_{1}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})&\dots&\mathbf{\boldsymbol{\mathsf{\partial}}}\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{m,\phi_{m}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)})\end{bmatrix}\begin{bmatrix}\Delta\mathbf{\boldsymbol{\mathsf{v}}}_{1}\\\
\vdots\\\ \Delta\mathbf{\boldsymbol{\mathsf{v}}}_{m}\end{bmatrix}$
$\displaystyle+\frac{1}{\Delta
x}\begin{bmatrix}(\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}\circ\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{1}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)}))\mathbf{\boldsymbol{\mathsf{1}}}\\\
\vdots\\\
(\mathbf{\boldsymbol{\mathsf{Q}}}^{(\mp)}\circ\mathbf{\boldsymbol{\mathsf{F}}}^{(\pm)}_{m}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(k)}))\mathbf{\boldsymbol{\mathsf{1}}}\end{bmatrix}=\underline{\mathbf{\boldsymbol{\mathsf{0}}}}.$

Looking at the $i$th row of the $\ell$th block in this linear system we find
that the scheme once again may be expressed on the locally conservative form

$\frac{v_{\ell_{i}}^{(k+1)}-u_{\ell_{i}}^{n}}{\Delta t}+\frac{1}{\Delta
x}\left(h_{\ell_{i}+\frac{1}{2}}^{(k+1)}-h_{\ell_{i}-\frac{1}{2}}^{(k+1)}\right)=0,\qquad
i=\dots,-1,0,1,\dots,$ (27)

where $\ell=1,\dots,m$. The numerical flux is given by

$h_{\ell_{i}+\frac{1}{2}}^{(k+1)}=\sum_{\kappa=1}^{m}\left(\hat{f}^{(\pm)}_{\ell,\phi_{\kappa}}(\underline{\mathbf{\boldsymbol{v}}}_{\kappa_{i}}^{(k)},\underline{\mathbf{\boldsymbol{v}}}_{\kappa_{i+1}}^{(k)})\Delta
v_{\kappa_{i+1}}\pm\hat{f}^{(\pm)}_{\ell,\phi_{\kappa}}(\underline{\mathbf{\boldsymbol{v}}}_{\kappa_{i+1}}^{(k)},\underline{\mathbf{\boldsymbol{v}}}_{\kappa_{i}}^{(k)})\Delta
v_{\kappa_{i}}\right)+\hat{f}^{(\pm)}_{\ell}(\underline{\mathbf{\boldsymbol{v}}}_{\ell_{i+1}}^{(k)},\underline{\mathbf{\boldsymbol{v}}}_{\ell_{i}}^{(k)}).$

Further, if
$\underline{\mathbf{\boldsymbol{\mathsf{g}}}}^{\prime}(\underline{\mathbf{\boldsymbol{u}}}\otimes\mathbf{\boldsymbol{\mathsf{1}}})$
is nonsingular for any vector $\underline{\mathbf{\boldsymbol{u}}}$ with non-
zero elements, then by the same argument as in the scalar case, the numerical
flux is consistent. We may thus conclude:

######  Theorem 6.

Choose $v_{\ell_{i}}^{(0)}=u_{\ell_{i}}^{n}$ for each $i$ and $\ell=1,\dots,m$
and suppose that
$\underline{\mathbf{\boldsymbol{\mathsf{g}}}}^{\prime}(\underline{\mathbf{\boldsymbol{u}}}\otimes\mathbf{\boldsymbol{\mathsf{1}}})$
is nonsingular for any vector $\underline{\mathbf{\boldsymbol{u}}}$ with non-
zero elements. Then Newton’s method, applied to the discretization (26), is
locally conservative and flux consistent.

Motivated by the conservation and flux consistency of Newton’s method we
conjecture that under conditions similar to 1, the numerical approximation
converges to a weak solution of the original conservation law (1). In 6 we
present numerical results supporting this conjecture.

The solution $\Delta\underline{\mathbf{\boldsymbol{v}}}$ to the linear systems
arising within Newton’s method will in practice be approximated using a second
iterative method. If we for this purpose consider pseudo-time iterations, then
by the conservation and flux consistency of Newton’s method, 2 applies. The
only change necessary is to choose the initial guess
$\Delta\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}$ as the null-vector.
We summarize these observations in the following:

######  Theorem 7.

Apply $K$ iterations of Newton’s method to the implicit scheme (2) with
initial guess
$\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(0)}=\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{n}$.
Approximate the solution to the linear system in the $j$th Newton iteration
using $N_{j}$ pseudo-time iterations with ERK methods and initial guess
$\Delta\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}=\underline{\mathbf{\boldsymbol{\mathsf{0}}}}$.
The Newton-pseudo-time iterations are locally conservative with numerical
fluxes consistent with $c\underline{\mathbf{\boldsymbol{f}}}$, where

$c=1-\prod_{l=0}^{\sum_{j}^{K}N_{j}-1}\phi_{l}(-\mu_{l}).$ (28)

###### Proof.

Since Newton’s method is locally conservative and flux consistent, the result
is a direct consequence of 2. ∎

##  5 Krylov subspace methods

Krylov subspace methods are a more common choice than pseudo-time iterations
for approximating the solutions of the linear systems arising within Newton’s
method. The goal of this section is to establish local conservation of Krylov
subspace methods, and thereby also of Newton-Krylov methods, by relating them
to pseudo-time iterations.

Let us consider a system of linear conservation laws and a correponding
linearly implicit discretization in the form (2). Due to the linearity of the
numerical flux, there is a matrix
$\underline{\mathbf{\boldsymbol{\mathsf{D}}}}$ such that the scheme can be
expressed as

$\underbrace{\left(\frac{\underline{\mathbf{\boldsymbol{\mathsf{I}}}}}{\Delta
t}+\frac{\underline{\mathbf{\boldsymbol{\mathsf{D}}}}}{\Delta
x}\right)}_{\underline{\mathbf{\boldsymbol{\mathsf{M}}}}}\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{n+1}=\underbrace{\frac{\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{n}}{\Delta
t}}_{\underline{\mathbf{\boldsymbol{\mathsf{d}}}}}.$

We assume that $\underline{\mathbf{\boldsymbol{\mathsf{M}}}}$ is an invertible
matrix.

By definition, the $(k+1)$st iteration of a Krylov subspace method finds an
approximate solution $\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(k+1)}$ to
$\underline{\mathbf{\boldsymbol{\mathsf{M}}}}\underline{\mathbf{\boldsymbol{\mathsf{u}}}}=\underline{\mathbf{\boldsymbol{\mathsf{d}}}}$
such that
$\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(k+1)}-\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(0)}\in\mathcal{K}_{k+1}(\underline{\mathbf{\boldsymbol{\mathsf{M}}}},\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)})$.
The Krylov subspace
$\mathcal{K}_{k+1}(\underline{\mathbf{\boldsymbol{\mathsf{M}}}},\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)})$
is defined as

$\mathcal{K}_{k+1}(\underline{\mathbf{\boldsymbol{\mathsf{M}}}},\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)}):=\mathrm{span}(\\{\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)},\underline{\mathbf{\boldsymbol{\mathsf{M}}}}\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)},\dots,\underline{\mathbf{\boldsymbol{\mathsf{M}}}}^{k}\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)}\\}),$
(29)

and
$\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)}=\underline{\mathbf{\boldsymbol{\mathsf{d}}}}-\underline{\mathbf{\boldsymbol{\mathsf{M}}}}\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(0)}$
is the initial residual. In other words, there are coefficients
$\alpha_{0},\dots,\alpha_{k}$ such that

$\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(k+1)}=\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(0)}+\left(\sum_{\ell=0}^{k}\alpha_{\ell}\underline{\mathbf{\boldsymbol{\mathsf{M}}}}^{\ell}\right)\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)}.$
(30)

We will now show that all Krylov subspace methods are locally conservative. To
this end, we first establish a connection between Krylov subspace methods and
pseudo-time iterations with explicit RK methods:

######  Lemma 8.

Let $\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(k+1)}$ denote the
$(k+1)$st iterate of a Krylov subspace method. There exists an ERK method with
$s=k+1$ stages such that the first iterate
$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(1)}$ of a pseudo-time
iteration using this method satisfies
$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(1)}=\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(k+1)}$.

###### Proof.

We prove the lemma in two steps: The first step is to show that

$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(1)}-\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}\in\mathcal{K}_{k+1}(\underline{\mathbf{\boldsymbol{\mathsf{M}}}},\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)})$

when
$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}=\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(0)}$.
In other words, a single pseudo-time iteration with a $(k+1)$-stage explicit
RK methods yields an approximation from the appropriate Krylov subspace. The
second step is to choose
$(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}},\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}},\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{c}}})$
for the Runge-Kutta method such that the coefficients
$\alpha_{0},\dots,\alpha_{k}$ in (30) are recovered.

For the first step, recall that pseudo-time iterations approximate solutions
to the steady state problem

$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}_{\tau}=\underline{\mathbf{\boldsymbol{\mathsf{d}}}}-\underline{\mathbf{\boldsymbol{\mathsf{M}}}}\underline{\mathbf{\boldsymbol{\mathsf{u}}}}=:\underline{\mathbf{\boldsymbol{\mathsf{r}}}}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}}),\qquad\underline{\mathbf{\boldsymbol{\mathsf{u}}}}(0)=\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}.$

A single iteration with an explicit RK method is obtained as

$\displaystyle\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}}$
$\displaystyle=\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}+\Delta\tau(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{I}}}})\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{R}}}}}}}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}}),$
(31) $\displaystyle\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(1)}$
$\displaystyle=\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}+\Delta\tau(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}}^{\top}\otimes\underline{\mathbf{\boldsymbol{\mathsf{I}}}})\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{R}}}}}}}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}}),$

where

$\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{R}}}}}}}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})=((\underline{\mathbf{\boldsymbol{\mathsf{d}}}}-\underline{\mathbf{\boldsymbol{\mathsf{M}}}}\underline{\mathbf{\boldsymbol{\mathsf{U}}}}_{1})^{\top},\dots,(\underline{\mathbf{\boldsymbol{\mathsf{d}}}}-\underline{\mathbf{\boldsymbol{\mathsf{M}}}}\underline{\mathbf{\boldsymbol{\mathsf{U}}}}_{s})^{\top})^{\top}=(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{d}}}})-(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{M}}}})\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}}.$

Inserting
$\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}}$
from (31) into the latter expression leads to

$\displaystyle\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{R}}}}}}}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})$
$\displaystyle=(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{d}}}})-(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{M}}}})\left[\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}+\Delta\tau(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{I}}}})\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{R}}}}}}}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})\right]$
$\displaystyle=(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{d}}}})-(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{M}}}}\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)})-\Delta\tau(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{M}}}})\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{R}}}}}}}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})$
$\displaystyle=(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)})-\Delta\tau(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{M}}}})\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{R}}}}}}}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}}).$

Solving for
$\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{R}}}}}}}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})$
results in

$\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{R}}}}}}}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{\underline{\mathbf{\boldsymbol{\mathsf{U}}}}}}})=\left[\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{I}}}}+\Delta\tau(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{M}}}})\right]^{-1}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)}).$

Inserting this expression into
$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(1)}$ from (31) yields

$\displaystyle\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(1)}$
$\displaystyle=\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}+\Delta\tau(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{I}}}})\left[\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{I}}}}+\Delta\tau(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{M}}}})\right]^{-1}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)})$
$\displaystyle=\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}+\underline{\mathbf{\boldsymbol{\mathsf{M}}}}^{-1}\left(\Delta\tau(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{b}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{M}}}})\left[\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{I}}}}+\Delta\tau(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{M}}}})\right]^{-1}(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{1}}}\otimes\underline{\mathbf{\boldsymbol{\mathsf{I}}}})\right)\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)}.$

We identify the paranthesised portion of this expression as
$\underline{\mathbf{\boldsymbol{\mathsf{I}}}}-\phi(-\Delta\tau\underline{\mathbf{\boldsymbol{\mathsf{M}}}})$,
where $\phi(z)$ is the stability function of the RK method. Thus,

$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(1)}=\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}+\underline{\mathbf{\boldsymbol{\mathsf{M}}}}^{-1}(\underline{\mathbf{\boldsymbol{\mathsf{I}}}}-\phi(-\Delta\tau\underline{\mathbf{\boldsymbol{\mathsf{M}}}}))\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)}.$
(32)

Now, $\phi(z)$ is a polynomial of degree $s=k+1$ and from (7) we see that its
constant term is $1$. Hence,
$\underline{\mathbf{\boldsymbol{\mathsf{M}}}}^{-1}(\underline{\mathbf{\boldsymbol{\mathsf{I}}}}-\phi(-\underline{\mathbf{\boldsymbol{\mathsf{M}}}}))$
is a polynomial in $\underline{\mathbf{\boldsymbol{\mathsf{M}}}}$ of degree
$k$. Consequently,
$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(1)}-\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}\in\mathcal{K}_{k+1}(\underline{\mathbf{\boldsymbol{\mathsf{M}}}},\underline{\mathbf{\boldsymbol{\mathsf{r}}}}^{(0)})$,
which completes the first step of the proof.

For the second step we must find an explicit RK method whose stability
polynomial satisfies
$z^{-1}(1-\phi(-z))=\sum_{\ell=0}^{k}\alpha_{\ell}z^{\ell}$. We do not require
this RK method to be accurate or even consistent; it suffices that it exists.

Consider the $(k+1)$-stage explicit RK method defined by the Butcher tableau

$\begin{array}[]{c|ccccc}0&0\\\ c_{1}&-a&0\\\ c_{2}&&-a&0\\\
\vdots&&&\ddots&\ddots\\\ c_{k}&&&&-a&0\\\
\hline\cr&b_{0}&b_{1}&\dots&b_{k-1}&b_{k}\end{array}$

By direct computation we can verify that
$(\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{I}}}-z\accentset{\text{\tiny$\rightarrow$}}{\mathbf{\boldsymbol{A}}})^{-1}$
is a lower triangular Toeplitz matrix with elements $(-az)^{\ell}$ on the
$\ell$th subdiagonal ($\ell=0$ being the main diagonal and $\ell=k$ being the
bottom left element). From (7), the stability function of this method is
therefore

$\phi(z)=1+z\sum_{i=0}^{k}\sum_{j=i}^{k}b_{j}(-az)^{i}.$

Consequently we are looking for coefficients $b_{j}$ that satisfy the relation

$\frac{1-\phi(-z)}{z}=\sum_{i=0}^{k}\sum_{j=i}^{k}b_{j}(az)^{i}=\sum_{\ell=0}^{k}\alpha_{\ell}z^{\ell}.$

Matching exponents of $z$, this amounts to solving the triangular system

$\begin{bmatrix}1&1&\dots&1&1\\\ 0&1&\dots&1&1\\\ &&\ddots&&\\\
0&0&\dots&0&1\end{bmatrix}\begin{pmatrix}b_{0}\\\ b_{1}\\\ \vdots\\\
b_{k}\end{pmatrix}=\begin{pmatrix}a^{-0}\alpha_{0}\\\ a^{-1}\alpha_{1}\\\
\vdots\\\ a^{-k}\alpha_{k}\end{pmatrix}.$

This system is uniquely solvable for any coefficients
$\alpha_{\ell},\,\ell=0,\dots,k$ and $a\neq 0$. Hence, an explicit RK method
can always be found such that a single pseudo-time iteration reproduces the
$(k+1)$th Krylov vector when the initial guess is given by
$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}=\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(0)}$.
∎

With 8 established, the following is a direct consequence of 2:

######  Theorem 9.

With the initial guess
$\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(0)}=\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{n}$,
Krylov subspace methods applied to the locally conservative linear
discretization (2) are locally conservative.

###### Proof.

For any iteration of any Krylov subspace method, 8 shows that there is an
equivalent explicit RK method that yields the same numerical solution in one
pseudo-time iteration. By 2, all explicit RK methods are locally conservative,
hence so are Krylov subspace methods. ∎

We do not know what the coefficients $\alpha_{0},\dots,\alpha_{k}$ are and
consequently we cannot identify a corresponding RK method. Further, each
$\alpha_{j}$ will in general change with every iteration. We can therefore not
know a priori if a given Krylov subspace method is flux consistent. We can
also not apply the extension of the Lax-Wendroff theorem in 4 since
$\alpha_{j}$ in general will not remain constant upon grid refinement.

On the other hand, 9 suggests that we can apply one pseudo-time iteration with
the explicit Euler method prior to applying the Krylov subspace method in
order to enforce consistency, without violating local conservation. Further
yet, since 2 is indifferent to the RK method used in previous iterations, we
can conclude that 9 also applies to _restarted_ Krylov methods, with each
restart corresponding to a pseudo-time iteration with its own RK method. In
fact, local conservation will be retained even if we swap Krylov method mid-
solve, or if we mix Krylov subspace methods and pseudo-time iterations.

Finally we note that 9 together with 6 implies that Newton-Krylov methods are
locally conservative:

######  Theorem 10.

Apply Newton’s method to the implicit scheme (2) with initial guess
$\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(0)}=\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{n}$.
Approximate the solution to the each linear system within the Newton
iterations using any Krylov subspace method with initial guess
$\Delta\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(0)}=\underline{\mathbf{\boldsymbol{\mathsf{0}}}}$.
Then the Newton-Krylov iterations are locally conservative.

###### Proof.

Since Newton’s method is locally conservative and flux consistent by 6, the
result is a direct consequence of 9. ∎

##  6 Numerical experiments

As target problem for the experiments in this section we use the 2D
compressible Euler equations,

$\begin{bmatrix}\rho\\\ \rho u\\\ \rho v\\\ \rho
E\end{bmatrix}_{t}+\begin{bmatrix}\rho u\\\ \rho u^{2}+p\\\ \rho uv\\\ (\rho
E+p)u\end{bmatrix}_{x}+\begin{bmatrix}\rho v\\\ \rho uv\\\ \rho v^{2}+p\\\
(\rho E+p)v\end{bmatrix}_{y}=0,$ (33)

posed on the domain $(x,y)\in(-5,15]\times(-5,5]$. Here, $\rho,u,v,E$ and $p$
respectively denote density, horizontal and vertical velocity components,
total energy per unit mass and pressure. The pressure is related to the other
variables through the equation of state

$p=(\gamma-1)\rho e,$

where $\gamma=1.4$ and $e=E-(u^{2}+v^{2})/2$ is the internal energy density.
The domain is taken to be periodic in both spatial coordinates. The setting is
the isentropic vortex problem [17] with initial conditions

$\displaystyle\rho_{0}$
$\displaystyle=\left(\frac{1-\epsilon^{2}(\gamma-1)M_{\infty}^{2}}{8\pi^{2}}\exp{(r)}\right)^{\frac{1}{\gamma-1}},$
$\displaystyle u_{0}$ $\displaystyle=1-\frac{\epsilon y}{2\pi}\exp{(r/2)},$
$\displaystyle v_{0}$ $\displaystyle=\frac{\epsilon x}{2\pi}\exp{(r/2)},$
$\displaystyle p_{0}$ $\displaystyle=\frac{\rho_{0}^{\gamma}}{\gamma
M_{\infty}^{2}},$

where $r=1-x^{2}-y^{2}$. Here, $\epsilon=5$ is the circulation and
$M_{\infty}=0.5$ is the Mach number. As the solution evolves in time, the
initial vortex propagates in the horizontal direction with unit speed.

This problem does not adhere to the assumptions made throughout the paper
since it is posed in 2D. Nevertheless, we will see shortly that the
theoretical results presented so far appear to hold also in this setting.

###  6.1 Convergence tests

#### Pseudo-time iterations: SSPRK3

We begin by corroborating the extension of the Lax-Wendroff theorem presented
in 4. With implicit Euler used in time, experimental verification of this
theorem was provided in [4]. However, 3 implies that the result holds also
when other implicit RK methods are used in time, so long as they fulfill
condition (16).

In the following experiments we discretize the Euler equations in space using
finite volumes with central fluxes and the three-stage Lobatto IIIC method in
time. This fully implicit RK method satisfies (16) with the vector
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}=(0,0,1)^{\top}$.
The simulations are run to time $t=0.1$ with space and time grids satisfying
by $\Delta y=\Delta x/2$ and $\Delta t=\Delta x/4$. In each time step, the
discrete solution is approximated using three pseudo-time iterations using the
explicit strong stability preserving method SSPRK3 from [18]. The pseudo-time
steps are chosen to be $\Delta\tau_{i}=\Delta t/\sqrt{i},\,i=1,2,3$. The
stability function of SSPRK3 is given by

$\phi(z)=1+z+z^{2}/2+z^{3}/6.$

4 therefore predicts that the solution of the scheme, if convergent,
approaches the solution to a modified conservation law with modification
constant

$c(\mu_{0},\mu_{1},\mu_{2})=1-\phi(\mu_{0})\phi(\mu_{1})\phi(\mu_{2})\approx
0.9101.$

For this particular problem, the presence of $c$ causes a corresponding
reduction of the vortex propagation speed. 1a shows the $L^{2}$-error in the
density component of the solution, as measured with respect to the exact
solution of the original system of conservation laws (33) and that of the
modified version. Clearly the numerical solution approaches the solution of
the modified equations, not the original ones.

Shown also is the error, measured with respect to the original conservation
law, of the same scheme but where flux consistency has been enforced using a
single pseudo-time iteration with explicit Euler and $\mu=1$. As expected,
this ensures convergence towards the correct solution.

#### Newton-SSPRK3

We now repeat the experiment but add two Newton iterations per time step. The
solutions to the linear systems are then approximated using the same three
pseudo-time iterations as before. In addition to the discretization described
previously, we also perform the experiment with a different one: In space, the
entropy conservative and kinetic energy preserving finite volume scheme of
Chandrashekar [7] is used and in time, the two-stage Radau IIA method, which
satisfies (16) with the vector
$\accentset{\text{\tiny$\leftarrow$}}{\mathbf{\boldsymbol{v}}}=(0,1)^{\top}$.

7 suggests that the resulting schemes are consistent with a modified system of
conservation laws with modification constant

$c=1-(\phi(\mu_{0})\phi(\mu_{1})\phi(\mu_{2}))^{2}\approx 0.9919.$

Rather than computing the Jacobian explicitly, it is approximated using finite
differences such that for any two vectors
$\underline{\mathbf{\boldsymbol{\mathsf{u}}}}$ and
$\underline{\mathbf{\boldsymbol{\mathsf{v}}}}$,

$\underline{\mathbf{\boldsymbol{\mathsf{g}}}}^{\prime}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}})\underline{\mathbf{\boldsymbol{\mathsf{v}}}}\approx\frac{\underline{\mathbf{\boldsymbol{\mathsf{g}}}}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}}+\epsilon\underline{\mathbf{\boldsymbol{\mathsf{v}}}})-\underline{\mathbf{\boldsymbol{\mathsf{g}}}}(\underline{\mathbf{\boldsymbol{\mathsf{u}}}})}{\epsilon},\quad\epsilon=\frac{10^{-7}}{\|\underline{\mathbf{\boldsymbol{\mathsf{v}}}}\|}.$
(34)

Although an extension of the Lax-Wendroff theorem incorporating Newton’s
method is not yet available, 1b indicates convergence towards the solution of
the modified equations, not of the original ones. Once again, this is remedied
by enforcing flux consistency with a single explicit Euler iteration. Here,
the errors corresponding to the two different space-time discretizations are
similar enough that the lines are on top of each other.

(a) SSPRK3

(b) Newton-SSPRK3

Figure 1: Density error upon grid refinement for the compressible Euler
equations. Convergence is seen towards the modified conservation laws
(Modified), not the original ones (Original), unless flux consistency is
enforced (Consistent). (a) Pseudo-time iterations using SSPRK3. (b) Newton’s
method with SSPRK3 as subsolver for the linear systems.

### Newton-Krylov: Fixed iterations

Newton-Krylov methods are a more common choice of solver than Newton with
pseudo-time iterations. 10 establishes that they are all locally conservative.
Conceivably, a Lax-Wendroff type result might be available for these methods,
although we have not presented one here. It is of interest to explore
experimentally if convergence is observed, and if so, towards what solution.

We repeat the previous experiment, this time with the three-stage Lobatto IIIC
method in time and Chandrashekar’s finite volume scheme in space, but replace
SSPRK3 with GMRES. Here we consider four cases: Using one Newton and one
Krylov iteration per time step (N1K1); one Newton and two Krylov iterations
(N1K2); one Newton iteration with GMRES run to a tolerance of $10^{-14}$ (N1),
i.e. effectively with ’exact’ linear solves; Newton-GMRES run until the Newton
residual is smaller than $10^{-15}$ (Exact), i.e. effectively with ’exact’
nonlinear solves. 2a shows the error with respect to the original conservation
law of the four schemes. The errors for N1K2, N1 and the exact solver are on
top of each other, suggesting that the discretization dominates the error. The
numerical solutions appear to converge to the correct solution. However, N1K1
displays a different behaviour, suggesting that the iteration error dominates.
The error curve appears to flatten as the grid is refined, although it cannot
be deduced whether it will continue towards zero or if it reaches a plateau.
Thus, it remains unclear if GMRES is flux consistent.

### Newton-Krylov: Tolerance

In practice, the number of Newton iterations will not be preset but rather
governed by a relative and/or an absolute tolerance. Here, we set both of
these to a value $tol$ and once again explore the convergence behaviour. To
set the tolerances of the GMRES iterations we follow the procedure described
by Eisenstat and Walker [8] with parameters $\gamma=\eta_{max}=0.9$; see [11,
Chapter 6] for details. The density errors using
$tol\in\\{10^{-3},10^{-4},10^{-5}\\}$ are shown in 2b.

Two distinct phenomena can be observed: Firstly, when $tol=10^{-5}$, the
convergence behaviour changes from one similar to the exact solver in 2a to
one resembling N1K1. Thus, the error is seen to change from being
discretization dominated to being iteration dominated as the grid is refined.
With $tol=10^{-4}$ the iteration error appears to dominate throughout.

Secondly, when $tol=10^{-3}$ the error eventually stops converging. This
behavior is explained by the observation that with a fixed tolerance, the
initial guess will be a sufficiently accurate approximation of the solution if
$\Delta t$ is small enough. In that case, the Newton and GMRES iterations are
terminated without updating the solution.

(a) Fixed iterations

(b) Tolerance

Figure 2: Density error upon grid refinement for the compressible Euler
equations. (a) The number of iterations is fixed: One Newton and one GMRES
iteration per time step (N1K1); One Newton and two GMRES iterations (N1K2);
one Newton with a nearly exact linear solver (N1); nearly exact Newton-GMRES
(Exact). (b) Tolerance governed Newton-GMRES with the Eisenstat-Walker
procedure.

###  6.2 Acceleration experiments

In [4], numerical experiments were performed indicating that considerable
efficiency gains can be made with pseudo-time iterations by enforcing flux
consistency. It is not clear whether Newton-Krylov methods are flux
consistent. It is therefore worthwhile exploring if enforced flux consistency
leads to efficiency gains also in this case. We use the tolerance governed
Newton-GMRES solver with the Eisenstat-Walker procedure and compare the
standard solver with one where flux consistency is enforced with explicit
Euler before each call to GMRES.

In this experiment we use Chandrashekar’s finite volume scheme with the
implicit Euler method in time to compute a single time step with $\Delta
t=0.1$. The running cost of the two solvers is measured in terms of the number
of evaluations of the full space discretization. Each Newton iteration
requires a single such evaluation; see (22). Each GMRES iteration also needs a
single function evaluation during the computation of the approximate Jacobian
matrix-vector product (34).

The $k$th pseudo-time step with explicit Euler applied within the $j$th Newton
iteration takes the form

$\frac{\Delta\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(k+1)}-\Delta\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(k)}}{\Delta\tau_{k}}+\underline{\mathbf{\boldsymbol{\mathsf{g}}}}^{\prime}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(j)})\Delta\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(k)}+\underline{\mathbf{\boldsymbol{\mathsf{g}}}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(j)})=\underline{\mathbf{\boldsymbol{\mathsf{0}}}}.$

Consider the case $k=0$ and recall from 7 that local conservation follows if
the initial guess is
$\Delta\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(0)}=\underline{\mathbf{\boldsymbol{\mathsf{0}}}}$.
Thus, enforcing flux consistency by one pseudo-time iteration with explicit
Euler simply amounts to setting
$\Delta\underline{\mathbf{\boldsymbol{\mathsf{u}}}}^{(1)}=-\Delta
t\underline{\mathbf{\boldsymbol{\mathsf{g}}}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(j)})$.
This single function evaluation is already computed within Newton’s method,
hence flux consistency comes at no additional cost. In fact, the only change
necessary to the solver is to alter the initial guess for GMRES from
$\Delta\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(0)}=\underline{\mathbf{\boldsymbol{\mathsf{0}}}}$
to $\Delta\underline{\mathbf{\boldsymbol{\mathsf{w}}}}^{(0)}=-\Delta
t\underline{\mathbf{\boldsymbol{\mathsf{g}}}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(j)})$.

Three cases with different CFL numbers are considered. 3a shows the number of
function evaluations required to reach a particular residual
$\|\underline{\mathbf{\boldsymbol{\mathsf{g}}}}(\underline{\mathbf{\boldsymbol{\mathsf{v}}}}^{(j)})\|$,
where the $L^{2}$-norm is used. 3b shows the function evaluations distributed
accross the Newton iterations. In all cases, the flux consistent initial guess
(dotted lines) reduces the number of necessary iterations for residuals
greater than roughly $10^{-3}$, compared to the regular solver (solid lines).
However, for smaller residuals the situation varies with the CFL number. At
the smallest CFL, flux consistency remains beneficial even for finer
tolerances. However, for the largest CFL the opposite trend is seen.

A possible explanation for these observations is that the explicit Euler
method introduces significant errors to the numerical solution when the CFL
number is large. Its stability region is small, hence problems with large CFL
numbers are unsurprising. It is possible to find other explicit Runge-Kutta
methods that enforce flux consistency, in principle with much larger stability
regions. However, such methods will necessarily have more stages and thus
impart additional costs on the solver.

Finally we note that at the smallest tolerances, no discernable difference is
seen between the two solvers. Presumably this happens because flux consistency
is achieved, either exactly or very nearly, by the standard Newton-GMRES
solver when the tolerance is small and the number of iterations is large. In
conclusion, whether enforced flux consistency is beneficial for Newton-GMRES
is case dependent.

(a)

(b)

Figure 3: Efficiency study for Newton-GMRES with standard (solid lines) and
flux consistent (dotted lines) initial guesses. (a) Function evaluations
needed to reach a given residual. (b) Function evaluations per Newton
iteration.

##  7 Conclusions

In this paper, the concepts of locally conservative and flux consistent
iterative methods have been introduced and shown to be of both theoretical and
practical interest. Based on earlier work in [4], it was shown that pseudo-
time iterations using explicit Runge-Kutta methods are locally conservative
but not necessarily flux consistent, when applied to conservative
discretizations of finite volume-type with a broad class of implicit Runge-
Kutta methods. For 1D problems, an extension of the Lax-Wendroff theorem
reveals convergence towards weak solutions of a temporally retarded system of
conservation laws. Each equation is modified in the same way, namely by a
particular scalar factor multiplying the spatial flux terms. Flux consistency,
and thereby convergence, is recovered through a technique based on using the
explicit Euler method.

Local conservation has further been established for all Krylov subspace
methods, with and without restarts, as well as for Newton’s method under the
assumption of bivariate fluxes. Thus it follows that Newton-Krylov methods are
locally conservative, although not necessarily flux consistent. Numerical
experiments with the 2D compressible Euler equations suggest that the role
enforced flux consistency is case dependent. Its effect diminishes as the
number of GMRES iterations grow, presumably because flux consistency is
achieved automatically.

## References

  * [1] T. Barth, Analysis of implicit local linearization techniques for upwind and TVD algorithms, in 25th AIAA Aerospace Sciences Meeting, 1987, p. 595.
  * [2] P. Birken, Numerical methods for unsteady compressible flow problems, Chapman and Hall/CRC, 2021.
  * [3] P. Birken, J. Bull, and A. Jameson, Preconditioned smoothers for the full approximation scheme for the RANS equations, Journal of Scientific Computing, 78 (2019), pp. 995–1022.
  * [4] P. Birken and V. Linders, Conservative iterative methods for implicit discretizations of conservation laws, arXiv preprint arXiv:2106.10088, (2021).
  * [5] P. D. Boom and D. W. Zingg, High-order implicit time-marching methods based on generalized Summation-By-Parts operators, SIAM Journal on Scientific Computing, 37 (2015), pp. A2682–A2709.
  * [6] 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.
  * [7] P. Chandrashekar, Kinetic energy preserving and entropy stable finite volume schemes for compressible Euler and Navier-Stokes equations, Communications in Computational Physics, 14 (2013), pp. 1252–1286.
  * [8] S. C. Eisenstat and H. F. Walker, Choosing the forcing terms in an inexact Newton method, SIAM Journal on Scientific Computing, 17 (1996), pp. 16–32.
  * [9] D. Jespersen and T. Pulliam, Flux vector splitting and approximate Newton methods, in 6th Computational Fluid Dynamics Conference Danvers, 1983, p. 1899.
  * [10] C. Junqueira-Junior, L. C. Scalabrin, E. Basso, and J. L. F. Azevedo, Study of conservation on implicit techniques for unstructured finite volume Navier–Stokes solvers, Journal of Aerospace Technology and Management, 6 (2014), pp. 267–280.
  * [11] C. T. Kelley, Iterative methods for linear and nonlinear equations, SIAM, 1995.
  * [12] P. Lax and B. Wendroff, Systems of conservation laws, tech. rep., LOS ALAMOS NATIONAL LAB NM, 1959.
  * [13] R. J. LeVeque, Numerical methods for conservation laws, vol. 3, Springer, 1992.
  * [14] V. Linders, J. Nordström, and S. H. Frankel, Properties of Runge-Kutta-Summation-By-Parts methods, Journal of Computational Physics, 419 (2020), p. 109684.
  * [15] J. Nordström and T. Lundquist, Summation-By-Parts in time, Journal of Computational Physics, 251 (2013), pp. 487–499.
  * [16] H. Ranocha, Some notes on Summation By Parts time integration methods, Results in Applied Mathematics, 1 (2019), p. 100004.
  * [17] C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, in Advanced numerical approximation of nonlinear hyperbolic equations, Springer, 1998, pp. 325–432.
  * [18] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, Journal of computational physics, 77 (1988), pp. 439–471.
  * [19] R. C. Swanson, E. Turkel, and C.-C. Rossow, Convergence acceleration of Runge–Kutta schemes for solving the Navier–Stokes equations, Journal of Computational Physics, 224 (2007), pp. 365–388.
  * [20] L. M. Versbach, V. Linders, R. Klöfkorn, and P. Birken, Theoretical and practical aspects of space-time DG-SEM implementations, arXiv preprint arXiv:2201.05800, (2022).
  * [21] G. Wanner and E. Hairer, Solving ordinary differential equations II, Springer Berlin Heidelberg, 1996.