# Well-balanced adaptive compact approximate Taylor methods
for systems of balance laws

H. Carrillo${}^{\hskip 1.42262pt\textrm{a}}$, E. Macca${}^{\hskip
1.42262pt\textrm{b}}$, Carlos Parés${}^{\hskip 1.42262pt\textrm{c}}$, G.
Russo${}^{\hskip 1.42262pt\textrm{d}}$
<EMAIL_ADDRESS>, TecNm region
carbonifera,Coahuila México.
<EMAIL_ADDRESS>, Department of Mathematics and
Computer Scince, University of Catania, Catania Italy.
<EMAIL_ADDRESS>, Department of Applied Mathematics,
University of Malaga, Malaga, Spain.
<EMAIL_ADDRESS>, Department of Mathematics and Computer
Scince, University of Catania, Catania Italy.

###### Abstract

Compact Approximate Taylor (CAT) methods for systems of conservation laws were
introduced by Carrillo and Parés in 2019. These methods, based on a strategy
that allows one to extend high-order Lax-Wendroff methods to nonlinear systems
without using the Cauchy-Kovalevskaya procedure, have arbitrary even order of
accuracy $2p$ and use $(2p+1)$-point stencils, where $p$ is an arbitrary
positive integer. More recently in 2021 Carrillo, Macca, Parés, Russo and
Zorío introduced a strategy to get rid of the spurious oscillations close to
discontinuities produced by CAT methods. This strategy led to the so-called
Adaptive CAT (ACAT) methods, in which the order of accuracy – and thus the
width of the stencils – is adapted to the local smoothness of the solution.
The goal of this paper is to extend CAT and ACAT methods to systems of balance
laws. To do this, the source term is written as the derivative of its
indefinite integral that is formally treated as a flux function. The well-
balanced property of the methods is discussed and a variant that allows in
principle to preserve any stationary solution is presented. The resulting
methods are then applied to a number of systems going from a linear scalar
conservation law to the 2D Euler equations with gravity, passing by the
Burgers equations with source term and the 1D shallow water equations: the
order and well-balanced properties are checked in several numerical tests.

††journal: Journal of LaTeX Templates

Keywords: High order fully-discrete schemes; High order reconstruction for
systems of balance laws; Finite difference schemes; Well-balanced methods.

##  1 Introduction.

This paper deals with the design of high-order well-balanced methods for
hyperbolic quasi-linear systems of balance laws

$U_{t}+F(U)_{x}=S(U)H_{x},$ (1.1)

with initial condition $U(x,0)=U_{0}(x)$, where
$U:\mathbb{R}\times[0,+\infty)\rightarrow\mathbb{R}^{d}$ is the unknown vector
field; $F:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is the flux function;
$S:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is the source term; and
$H:\mathbb{R}\rightarrow\mathbb{R}$ is a known function. PDE systems of this
form appear in many fluid models in different contexts: shallow water models,
multiphase flow models, gas dynamic, elastic wave equations, etc.

More precisely, we focus on the extension of high-order Lax-Wendroff methods
to systems (1.1). For linear systems of conservation laws, these methods are
based on Taylor expansions in time in which the time derivatives are
transformed into spatial derivatives using the governing equations [22, 32,
15]. The discretization of the spatial derivatives by means of centered high-
order differentiation formulas leads then to numerical methods with arbitrary
order of accuracy.

The main difficulty to extend Lax-Wendroff methods to nonlinear problems comes
from the transformation of time derivatives into spatial derivatives: the use
of the equations through the Cauchy-Kovalesky (CK) procedure may be
impractical from the computational point of view because it often requires
extended symbolic calculus, ended up with inefficient codes. In the context of
ADER methods introduced by Toro and collaborators (see [23, 31, 29]), this
difficulty has been circumvented by replacing the CK procedure by local space-
time problems that are solved with a Galerkin method: see [12], [11].

For systems of conservation laws

$U_{t}+F(U)_{x}=0,$ (1.2)

the CK procedure was avoided in [36] by computing time derivatives in a
recursive way using high-order centered differentiation formulas combined with
Taylor expansions in time. Nevertheless, if $(2p+1)$-point differentiation
formulas are used to compute spatial and temporal derivatives, the resulting
method use $(4p+1)$-point stencils while Lax-Wendroff methods for linear
systems use $(2p+1)$-point ones. In [4] a variant of these methods that use
$(2p+1)$-point stencils, the so-called Compact Approximated Taylor methods
(CAT), was introduced. CAT methods were shown to reduce to the standard high-
order Lax-Wendroff methods when applied to linear problems. The technique used
to reduce the length of the stencils increases the computational cost of a
time step compared to the methods introduced in [36]: the Taylor expansions
are computed locally, so that the total number of expansions needed to update
the numerical solution is multiplied by $(2p+1)$. Nevertheless, CAT methods
have better stability properties allowing larger time steps, thus compensating
the extra cost per time iteration: see [4].

As it happens with high-order Lax-Wendroff methods for linear systems, their
extensions to nonlinear problems produce spurious oscillations close to
discontinuities and a strategy is needed to get rid of them. In [36] they were
combined with WENO reconstructions to compute the first time derivatives. CAT
methods were also combined with WENO in [4] and [5] to avoid oscillations near
discontinuities. Nevertheless this combination was not optimal: while the best
CAT methods are those of even order, WENO methods have odd accuracy order.
Moreover, the restriction on the time step imposed by WENO methods may spoil
the advantages of the better stability property of CAT methods. To avoid this,
in [3] a new version of CAT methods, the so-called Adaptive CAT (ACAT)
methods, was introduced in which the oscillations near discontinuities are
cured by adapting the order of accuracy – and thus the width of the stencils –
to the smoothness of the solution. To do this, a class of smoothness
indicators was introduced.

The main goal of this article is to extend ACAT methods to balance laws (1.1)
and to its 2D counterpart. To do this, we follow the strategy in [14] (see
also [10]) that consists in writing the source term as the derivative of the
indefinite integral

$\int_{-\infty}^{x}S(U(\sigma,t))H_{x}(\sigma)\,d\sigma,$

that is then formally treated as a new flux function. Please note that,
although this technique is applied here to derive ACAT methods for (1.1), it
can be applied in principle to other extensions of high-order Lax-Wendroff
methods to nonlinear problems using any technique to remove the spurious
oscillations.

Systems of balance laws (1.1) have non-trivial stationary solutions that
satisfy the ODE system

$F(U)_{x}=S(U)H_{x}.$ (1.3)

The objective of well balanced schemes is to preserve exactly or with enhanced
accuracy some of these steady state solutions. In the context of shallow water
equations, Bermúdez and Vázquez-Cendón introduced in [1] the condition called
C-property: a scheme is said to satisfy this condition if it preserves the
water at rest solutions. Since then, many different numerical methods that
satisfy this property have been introduced in the literature: see [2], [33]
and their references. In the framework of finite difference methods, high-
order schemes that satisfy the C-property were introduced in [6] and [35]:
while the former was based on the formal writing of the system in conservative
form based on the above mentioned technique, the latter relied on the
expression of the source term as a function of variables that are constants
for the stationary solutions to be preserved: see [34]. A similar technique
based on the reconstruction of equilibrium variables has been proposed in
[28]. In [24] a general technique to derive high-order well-balanced finite-
difference methods for systems of balance-laws was introduced. The strategy,
inspired on the general technique for finite volume methods discussed in [8],
was as follows: let $U_{i}$ be the numerical approximation of the solution
$U(x_{i},t)$ at the node $x_{i}$ at time $t$ and let $U^{*}_{i}$ be the
stationary solution satisfying the Cauchy problem:

$\begin{cases}&\displaystyle F(U^{*}_{i})_{x}=S(U^{*}_{i})H_{x},\\\\[5.0pt]
&\displaystyle U^{*}_{i}(x_{i})=U_{i}.\end{cases}$ (1.4)

Then, if $U^{*}_{i}$ can be found, one has trivially

$S(U_{i})H_{x}(x_{i})=S(U^{*}_{i}(x_{i}))H_{x}(x_{i})=F(U^{*}_{i}(x_{i}))_{x}.$
(1.5)

Therefore, locally the system of balance laws can be written in conservation
form as follows

$U_{t}+(F(U)-F(U^{*}_{i}(x)))_{x}=0.$

The numerical method is obtained then by discretizating this conservative form
by means of high-order WENO reconstruction of this extended ’flux’ function.
We will follow this strategy here to derive well-balanced ACAT methods.

This paper is organized as follows: in Section 2, CAT methods for systems of
conservation laws (1.2) are recalled. In Section 3, these methods are extended
to systems of balance laws (1.1): after obtaining a first high-order version
of the methods for systems of balance laws, the well-balanced property is
discussed and a second version is introduced that preserves in principle any
stationary solution. In Section 4 the technique introduced in [3] based on the
order-adaption of the methods is recalled and it is applied to the 2 versions
of CAT methods derived in Section 3, which leads to ACAT and well-balanced
ACAT methods for systems of balance laws. Section 5 is devoted to the
extension of the methods to 2D problems. In Section 6 the methods are applied
to a number of problems: a 1D scalar balance law, Burgers equation with source
term, the 1D shallow water model, and 2D Euler equations with gravitational
potential. Finally, in the last section, we draw some conclusions.

##  2 Compact Approximate Taylor Methods for systems of conservation laws

First of all, let us recap the expression of CAT methods for a 1D system of
conservation laws (1.2) with initial condition $U(x,0)=U_{0}(x)$. In the
CAT$2P$ method, a Taylor expansion in time is used to update the numerical
solution:

$U_{i}^{n+1}=U_{i}^{n}+\sum_{k=1}^{2P}\frac{(\Delta t)^{k}}{k!}U^{(k)}_{i},$
(2.1)

where $\\{x_{i}\\}$ are the nodes of a uniform mesh of step $\Delta x;$
$U_{i}^{n}$ is an approximation of the value of the exact solution
$U(x_{i},t_{n})$ at time $t_{n}=n\Delta t$ at $x_{i}$; and $U_{i}^{(k)}$ is an
approximation of $\partial^{k}_{t}U(x_{i},t_{n})$. Following the strategy
introduced in [36] to derive Approximate Taylor methods, the Cauchy-
Kovalevskaya procedure is avoided using the equality

$\partial_{t}^{k}U=-\partial_{x}\partial_{t}^{k-1}F(U),$ (2.2)

satisfied by smooth enough solutions.

In the Approximate Taylor methods introduced in [36], global approximations
that have to be computed only once at every node were used. Although the use
of local approximations may increase the number of calculations, it allows us
to prevent the increase of the stencil observed in the methods introduced in
this reference: indeed, CAT$2P$ methods are written in conservative form

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Bigl{(}F_{i-\frac{1}{2}}^{P}-F_{i+\frac{1}{2}}^{P}\Bigr{)},$ (2.3)

and only the values

$U_{i-P+1}^{n},\dots,U_{i+P}^{n}$

at the stencil $\mathcal{S}^{P}_{i+1/2}$

$\mathcal{S}^{P}_{i+1/2}=\\{x_{i-P+1},\dots,x_{i+P}\\},$

are used to compute the numerical flux $F_{i+\frac{1}{2}}^{P}$, so that the
stencil used to update the solution at the $i$th node is
$\\{x_{i-P},\dots,x_{i+P}\\}$, while in [36] the stencil for the same order of
accuracy required $4P$ points. More precisely, the time derivatives of the
solution will be approximated by applying a formula of numerical
differentiation for first order spatial derivatives to some approximations

${F}^{(k-1)}_{i,j}\approx\partial_{t}^{k-1}F(U)(x_{i+j},t_{n}),\quad
j=-P+1,\dots,P$ (2.4)

that will be computed using recursively Taylor expansions in time. In this
notation, $P$ is an arbitrary positive integer and $j$ is a local coordinate
in the stencil $\mathcal{S}^{P}_{i+1/2}\equiv\\{-P+1,\dots,P\\}$ i.e.
$F^{(k)}_{i,j}$ is the approximation of $\partial_{t}^{k}F(U)$ at time $t_{n}$
at the node of local coordinate $j$ of the stencil $\mathcal{S}^{P}_{i+1/2}$
which is $x_{i+j}$. These approximations are local in the following sense: let
us suppose that $i_{1}+j_{1}=i_{2}+j_{2}=l$, i.e. $x_{l}$ belongs to
$\mathcal{S}^{P}_{i_{1}+1/2}$ and $\mathcal{S}^{P}_{i_{2}+1/2}$ with local
coordinates $j_{1}$ and $j_{2}$ respectively. Then $F^{(k)}_{i_{1},j_{1}}$ and
$F^{(k)}_{i_{2},j_{2}}$ are, in general, two different approximations of
$\partial_{t}^{k}F(U)(x_{l},t_{n})$. Moreover, in [4] it has been shown that
CAT$2P$ methods reduce to the $2P$-order Lax-Wendroff method for linear
systems. This implies the linear stability for these methods under the usual
CFL requirement (which is not the case for the methods introduced in [36]).

Since numerical differentiation plays a fundamental role in the algorithm,
before giving the expression of the numerical flux, let us introduce the
notation to describe the formulas that will be used. Given two positive
integers $p$, $k$, an index $i$, and a real number $q$, we will represent by

$f^{(k)}(x_{i}+q\Delta x)\approx A^{k,q}_{p}(f_{i},\Delta x)=\frac{1}{\Delta
x^{k}}\sum_{j=-p+1}^{p}\gamma^{k,q}_{p,j}f(x_{i+j}),\quad\quad
k=0,\ldots,2p-1,$ (2.5)

the interpolatory formula that approximates the $k$-th derivative of a
function $f$ at the point $x_{i}+q\Delta x$ using its values at the $2p$
points $x_{i-p+1},\dots,x_{i+p}$. For $k=0$,

$f(x_{i}+q\Delta x,t_{n})=A^{0,q}_{p}(f_{i},\Delta
x)=\sum_{j=-p+1}^{p}\gamma^{0,q}_{p,j}f(x_{i+j}),$

represents the value at $x_{i}+q\Delta x$ of the Lagrange polynomial that
interpolates the values of $f$ at the points $x_{i-p+1},\dots,x_{i+p}$.

######  Remark 2.1.

The coefficients $\gamma^{k,q}_{p,j}$ of the differentiation formulas can be
recursively computed using the algorithm introduced in [13]. See also [4].

The following notation

$\displaystyle\partial_{x}^{k}f(x_{i}+q\Delta x,t_{n})\approx
A^{k,q}_{p}(f_{i,{*}}^{n},\Delta x)$ $\displaystyle=\frac{1}{\Delta
x^{k}}\sum_{j=-p+1}^{p}\gamma^{k,q}_{p,j}f_{i+j},$ (2.6)
$\displaystyle\partial_{t}^{k}f(x_{i},t_{n})\approx
A^{k,0}_{p}(f^{n,{*}}_{i},\Delta t)$ $\displaystyle=\frac{1}{\Delta
t^{k}}\sum_{j=-p+1}^{p}\gamma^{k,q}_{p,j}f_{i+j},$ (2.7)

will be used to indicate that the formula is applied to some approximations
$f_{i}$ of $f$ and not to its exact point values $f(x_{i})$. In cases where
there are two or more indices, the symbol ${*}$ will be used to indicate to
which index (space or time) the differentiation is applied. From now on, since
all the formulas are computed at time $t=t_{n}$ we avoid the extra index $n$
on equations (2.6) and (2.7). For instance, the following approximations will
be used in the algorithm adopted to compute the numerical fluxes:

$\displaystyle\partial^{k}_{t}U(x_{i+j},t_{n})\approx-
A_{P}^{1,j}\Bigl{(}F_{i,{*}}^{(k-1)},\Delta x\Bigr{)}=\frac{1}{\Delta
x}\sum_{l=-P+1}^{P}\gamma^{1,j}_{P,l}F_{i,l}^{(k)},$ (2.8)
$\displaystyle\partial^{k}_{t}F(U)(x_{i+j},t_{n})\approx
A_{P}^{k,0}\Bigl{(}F_{i,j}^{k,{*}},\Delta t\Bigr{)}=\frac{1}{\Delta
t^{k}}\sum_{r=-P+1}^{P}\gamma^{k,0}_{P,r}F_{i,j}^{k,n+r},$ (2.9)
$\displaystyle\partial^{k}_{t}F(U)\left(x_{i}+\frac{\Delta
x}{2},t_{n}\right)\approx A^{0,1/2}_{P}\left({F}_{i,{*}}^{(k)},\Delta
x\right)=\sum_{j=-P+1}^{P}\gamma^{0,1/2}_{P,j}{F}_{i,j}^{(k)}.$ (2.10)

In (2.8), that is the discrete version of (2.2), numerical differentiation in
space is used to approximate the time derivative of the solution at $x_{i+j}$
from the local approximations $F_{i,l}^{(k)}$, $l=-P+1,\dots,P$. In (2.9),
numerical differentiation in time is used to approximate the $k$-th time
derivative of $F(U)$ at $x_{i+j}$ from some approximations $F_{i,j}^{k,n+r}$
of $F(U)(x_{i+j},t_{n+r})$, $r=-P+1,\dots,P$. Finally, in (2.10), Lagrange
interpolation is used to approximate the value of the $k$-th time derivative
of $F(U)$ at $x_{i}+\Delta x/2$ at time $t_{n}$ from $F_{i,l}^{(k)}$,
$l=-P+1,\dots,P$.

Using this notation, the expression of the numerical flux is as follows:

$F_{i+\frac{1}{2}}^{P}=\sum_{k=1}^{2P}\frac{(\Delta
t)^{k-1}}{k!}A^{0,1/2}_{P}({F}_{i,{*}}^{(k-1)},\Delta
x)=\sum_{k=1}^{2P}\frac{(\Delta
t)^{k-1}}{k!}\sum_{j=-P+1}^{P}\gamma_{P,j}^{0,\frac{1}{2}}F_{i,j}^{(k-1)},$
(2.11)

where the time derivatives of the flux are computed by the following iterative
algorithm (see [4] for details):

  * 1.

Define $F_{i,j}^{(0)}:=F(U_{i+j}^{n})$ for all $j=-P+1,\ldots,P.$

  * 2.

For $k=1,\ldots,2P-1:$

    * (a)

Compute for all $j=-P+1,\ldots,P$

$U_{i,j}^{(k)}=-A_{P}^{1,j}\Bigl{(}F_{i,{*}}^{(k-1)},\Delta x\Bigr{)}.$

    * (b)

Define for all $j,r=-P+1,\ldots,P$

$F_{i,j}^{k,n+r}:=F\Bigl{(}U_{i,j}^{k,n+r}\Bigr{)},$

where $U_{i,j}^{k,n+r}$ is the approximation of $U(x_{i+j},t_{n+r})$ given by
the approximate Taylor expansion in time:

$U_{i,j}^{k,n+r}=U_{i+j}^{n}+\sum_{m=1}^{k}\frac{(\Delta
t)^{m}}{m!}U_{i,j}^{(m)}.$

    * (c)

Compute for all $j=-P+1,\ldots,P$

$F_{i,j}^{(k)}=A_{P}^{k,0}\Bigl{(}F_{i,j}^{k,{*}},\Delta t\Bigr{)}.$

######  Remark 2.2.

Observe that the computation of the numerical flux $F_{i+\frac{1}{2}}^{P}$
requires the approximation of $U$ at the nodes of a space-time grid of
$2P\times 2P$ points: $U_{i,j}^{k,n+r}$, $j,r=-P+1,\ldots,P$ (see Figure 2.1).
The approximations of the solution $U$ at times $(n-P+1)\Delta t$, …,
$(n-1)\Delta t$ are different from the ones already computed in the previous
steps: $U_{i+j}^{n-p}$, …, $U_{i+j}^{n-1}$. In other words, the discretization
in time is not based on a multistep method but in a one-step one: in fact it
can be interpreted as a RK method whose stages are $\tilde{U}^{n+r}_{i,j}$,
$r=-p,\dots,p$: see [3].

Figure 2.1: Local space-time grid where approximations of $U$ are computed to
calculate $F^{P}_{i+1/2}$ with $P=2$. For simplicity a pair $j,r$ represents
the point $(x_{i+j},t_{n+r})$. Taylor expansions in time are used to obtain
these approximations following the blue lines. These Taylor expansions are
centered in the points lying on the black line.

In [4] it has been shown that CAT$2P$ has order of accuracy $2P$ in time and
it is linearly stable under the usual CFL condition.

For $P=1$, the expression of CAT2 numerical flux reduces to:

$F_{i+\frac{1}{2}}^{1}=\frac{1}{4}\Bigl{(}F_{i,0}^{1,n+1}+F_{i,1}^{1,n+1}+F(U_{i}^{n})+F(U_{i+1}^{n})\Bigr{)}.$
(2.12)

where, for $j=0,1,$

$F_{i,j}^{1,n+1}=F\Bigl{(}U_{i,j}^{1,n+1}\Bigr{)}=F\Bigl{(}U_{i+j}^{n}-\frac{\Delta
t}{\Delta x}(F(U_{i+1}^{n})-F(U_{i}^{n}))\Bigr{)}$

Consequently, the corresponding numerical scheme (2.3) writes as follows:

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{4\Delta
x}\Bigl{(}F_{i-1,0}^{1,n+1}+F_{i-1,1}^{1,n+1}-F_{i,0}^{1,n+1}-F_{i,1}^{1,n+1}+F(U_{i-1}^{n})-F(U_{i+1}^{n})\Bigr{)}.$
(2.13)

######  Remark 2.3.

The second order CAT2 numerical flux (2.12) could be seen as a new one-step
Jacobian-free extension of Lax-Wendroff scheme that is even different from the
two-step MacCormack and Richtmyer method.

##  3 Compact Approximate methods for systems of balance laws

###  3.1 CAT2P methods

The goal of this section is to extend CAT$2P$ methods to systems of balance
laws (1.1) with initial condition $U(x,0)=U_{0}(x)$. Many authors have treated
the source term as divergence of flux introducing an elliptic equation, see
[26, 27]. To avoid the introduction of a new equation we follow the strategy
in [14] (see also [10]) in which the sistem (1.1) is first written in
conservative form through the definition of a ‘combined flux’ formed by the
sum of flux function $F$ and the indefinite integral of the source term: more
precisely, let us introduce the function $\mathcal{F}$ given by

$\mathcal{F}(U)(x,t)=F(U(x,t))-\int_{-\infty}^{x}S(U(\sigma,t))H_{x}(\sigma)\,d\sigma,$
(3.1)

assuming that the integral is finite. Then, the equality

$\mathcal{F}(U)_{x}=F(U)_{x}-S(U)H_{x},$

allows one to write the system of balance laws (1.1) in the form

$U_{t}+\mathcal{F}(U)_{x}=0.$ (3.2)

Then, the formal expression of CAT$2P$ is given by:

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Big{(}\mathfrak{F}_{i-\frac{1}{2}}^{P}-\mathfrak{F}_{i+\frac{1}{2}}^{P}\Big{)},$
(3.3)

where

$\mathfrak{F}_{i+\frac{1}{2}}^{P}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A_{P}^{0,1/2}\Big{(}\mathcal{F}_{i,{*}}^{(k-1)}\Big{)}.$ (3.4)

Here, $\mathcal{F}^{(k)}_{i,j}$ are local approximation of
$\partial_{t}^{(k)}\mathcal{F}(U)(x_{i+j},t_{n})$ that are computed by
adapting the algorithm described in Section 2. Formally, the algorithm is as
follows:

  * 1.

Define

$\displaystyle F_{i,j}^{(0)}$ $\displaystyle:=$ $\displaystyle
F(U_{i+j}^{n}),\quad j=-P+1,\ldots,P;$ $\displaystyle I_{i,j}^{(0)}$
$\displaystyle:=$
$\displaystyle\int_{-\infty}^{x_{i+j}}S(U(x,t_{n}))H_{x}(x)\,dx,\quad
j=-P+1,\dots,P.$

  * 2.

For $k=1,\ldots,2P-1:$

    * (a)

Compute for all $j=-P+1,\ldots,P$

$U_{i,j}^{(k)}=-A_{P}^{1,j}\Bigl{(}F_{i,{*}}^{(k-1)},\Delta
x\Bigr{)}+A_{P}^{1,j}\Bigl{(}I_{i,{*}}^{(k-1)},\Delta x\Bigr{)}.$

    * (b)

Define for all $j,r=-P+1,\ldots,P$

$\displaystyle I_{i,j}^{n+r}$ $\displaystyle:=$
$\displaystyle\int_{-\infty}^{x_{i+j}}S(U(x,t_{n+r}))H_{x}(x)\,dx,$
$\displaystyle F_{i,j}^{k,n+r}$ $\displaystyle:=$ $\displaystyle
F\Bigl{(}U_{i,j}^{k,n+r}\Bigr{)},$

where $U_{i,j}^{k,n+r}$ is the approximation of $U(x_{i+j},t_{n+r})$ given by
the Taylor expansion in time:

$U_{i,j}^{k,n+r}=U_{i+j}^{n}+\sum_{m=1}^{k}\frac{(\Delta
t)^{m}}{m!}U_{i,j}^{(m)}.$

    * (c)

Compute for all $j=-P+1,\ldots,P$

$F_{i,j}^{(k)}=A_{P}^{k,0}\Bigl{(}F_{i,j}^{k,{*}},\Delta t\Bigr{)},\quad
I_{i,j}^{(k)}=A_{P}^{k,0}\Bigl{(}I_{i,j}^{{*}},\Delta t\Bigr{)}.$

The ’numerical fluxes’ are then defined by:

$\mathfrak{F}_{i+\frac{1}{2}}^{P}=F^{P}_{i+1/2}-I^{P}_{i+1/2}$ (3.5)

where

$\displaystyle F^{P}_{i+1/2}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,1/2}_{P}(F_{i,{*}}^{(k-1)},\Delta x),$ (3.6) $\displaystyle
I^{P}_{i+1/2}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,1/2}_{P}(I_{i,{*}}^{(k-1)},\Delta x).$ (3.7)

This algorithm is formal, since it requires the computation of integrals that
depend on the exact solution in intervals of the form $(-\infty,x_{i+j}]$. In
order to be computationally implementable, let us first rewrite it using only
integrals in bounded intervals. To do that, the key point is the following
chain of equalities:

$\displaystyle A_{P}^{1,j}\Bigl{(}I_{i,{*}}^{(k-1)},\Delta x\Bigr{)}$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta
x}\sum_{s=-P+1}^{P}\gamma_{P,s}^{1,j}I_{i,s}^{(k-1)}$ $\displaystyle=$
$\displaystyle\frac{1}{\Delta
x}\sum_{s=-P+1}^{P}\gamma_{P,s}^{1,j}A_{P}^{k-1,0}\Bigl{(}I_{i,s}^{{*}},\Delta
t\Bigr{)}$ $\displaystyle=$ $\displaystyle\frac{1}{\Delta x\Delta
t^{k-1}}\sum_{s=-P+1}^{P}\gamma_{P,s}^{1,j}\sum_{r=-P+1}^{P}\gamma_{P,r}^{k-1,0}I_{i,s}^{k,n+r}$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta x\Delta
t^{k-1}}\sum_{r=-P+1}^{P}\gamma_{P,r}^{k-1,0}\sum_{s=-P+1}^{P}\gamma_{P,s}^{1,j}\int_{-\infty}^{x_{i+s}}S(U(x,t_{n+r}))H_{x}(x)\,dx$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta x\Delta
t^{k-1}}\sum_{r=-P+1}^{P}\gamma_{P,r}^{k-1,0}\sum_{s=-P+1}^{P}\gamma_{P,s}^{1,j}\Bigl{(}\int_{-\infty}^{x_{i+s}}S(U(x,t_{n+r}))H_{x}(x)\,dx$
$\displaystyle\qquad-\int_{-\infty}^{x_{i-P+1}}S(U(x,t_{n+r}))H_{x}(x)\,dx\Bigr{)}$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta x\Delta
t^{k-1}}\sum_{r=-P+1}^{P}\gamma_{P,r}^{k-1,0}\sum_{s=-P+1}^{P}\gamma_{P,s}^{1,j}\int_{x_{i-P+1}}^{x_{i+s}}S(U(x,t_{n+r}))H_{x}(x)\,dx$

where the identity

$\sum_{s=-P+1}^{P}\gamma^{1,j}_{P,s}=0,\quad j=-P+1,\ldots,P,$

has been used: remember that an interpolatory formula of numerical
differentiation that uses $2P$ points is exact at least for polynomials of
degree $2P-1$ and thus it is exact for constant polynomials. Therefore, if the
formula is applied to the constant polynomial $p\equiv 1$, we have

$0=p^{\prime}(x_{i+j})=A^{1,j}_{P}(p_{i,{*}},\Delta
x)=\sum_{s=-P+1}^{P}\gamma^{1,j}_{P,s}p(x_{i+s})=\sum_{s=-P+1}^{P}\gamma^{1,j}_{P,s}.$

By introducing the notation

$\displaystyle I_{i,j,l}^{m}$ $\displaystyle:=$
$\displaystyle\int_{x_{i+j}}^{x_{i+l}}S(U(x,t_{m}))H_{x}(x)\,dx,$
$\displaystyle I_{i,j,l}^{(k-1)}$ $\displaystyle:=$ $\displaystyle
A_{P}^{k-1,0}\Bigl{(}I_{i,j,l}^{{*}},\Delta t\Bigr{)},$

we obtain

$\displaystyle A_{P}^{1,j}\Bigl{(}I_{i,{*}}^{(k-1)},\Delta x\Bigr{)}$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta x\Delta
t^{k-1}}\sum_{r=-P+1}^{P}\gamma_{P,r}^{k-1,0}\sum_{s=-P+1}^{P}\gamma_{P,s}^{1,j}I_{i,-P+1,s}^{n+r}$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta
x}\sum_{r=-P+1}^{P}\gamma_{P,r}^{k-1,0}I_{i,-P+1,s}^{(k-1)}$ $\displaystyle=$
$\displaystyle A_{P}^{1,j}\Bigl{(}I_{i,-P+1,{*}}^{(k-1)},\Delta x\Bigr{)},$

where only integrals in intervals of the form $[x_{i-P+1},x_{i+s}]$ appear.
Observe that $I_{i-P+1,i-P+1}^{n+r}=0$ for all $k$.

Concerning the expression of the numerical method, observe that:

$\displaystyle I^{P}_{i+1/2}-I^{P}_{i-1/2}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}\left(A^{0,1/2}_{P}(I_{i,{*}}^{(k-1)},\Delta
x)-A^{0,1/2}_{P}(I_{i-1,{*}}^{(k-1)},\Delta x)\right)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}\sum_{j=-P+1}^{P}\gamma_{P,j}^{0,1/2}\left(I_{i,j}^{(k-1)}-I_{i-1,j}^{(k-1)}\right)$
$\displaystyle=$ $\displaystyle\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}\sum_{j=-P+1}^{P}\gamma_{P,j}^{0,1/2}A_{P}^{k-1,0}\left(I_{i,j}^{{*}}-I_{i-1,j}^{{*}},\Delta
t\right).$

Since

$I_{i,j}^{n+r}-I_{i-1,j}^{n+r}=I_{i,j-1,j}^{n+r}=\int_{x_{i+j-1}}^{x_{i+j}}S(U(x,t_{n+r}))H_{x}(x)\,dx,$
(3.8)

if we define

$\mathcal{I}_{i,j}^{(k-1)}=A_{P}^{k-1,0}\Bigl{(}{I}_{i,j-1,j}^{{*}},\Delta
t\Bigr{)},$ (3.9)

we have

$I^{P}_{i+1/2}-I^{P}_{i-1/2}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,1/2}_{P}(\mathcal{I}_{i,{*}}^{(k-1)},\Delta x),$

so that (3.3) can be written in equivalent form

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Big{(}F_{i-\frac{1}{2}}^{P}-{F}_{i+\frac{1}{2}}^{P}+S_{i}^{P}\Big{)},$
(3.10)

where

$S_{i}^{P}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,1/2}_{P}(\mathcal{I}_{i,{*}}^{(k-1)},\Delta x).$ (3.11)

Observe that only integrals (3.8) in intervals of length $\Delta x$ appear in
the expression of the numerical source term.

Finally, in order to have an implementable algorithm, all the integrals
appearing in it are approximated using quadrature formulas combined with the
approximations $U_{i,j}^{k,n+r}$ of the exact solution that are available at
every stage. To do this, given $i$ and $j=-P+2,\dots,P$, we consider at
$[x_{i+j-1},x_{i+j}]$ the interpolatory quadrature formula

$\int_{x_{i+j-1}}^{x_{i+j}}f(x)\,dx\approx\Delta
x\sum_{s=-P+1}^{P}a^{i,j}_{P,s}f(x_{i+s})$

whose nodes are $x_{i+s}$, $s=-P+1,\dots,P$. This formula will be used to
approximate the integrals appearing at the $k$-th stage of the algorithm as
follows: given two indices $j_{1}<j_{2}$

$\displaystyle
I_{i,j_{1},j_{2}}^{m}\approx\widetilde{I}^{k,m}_{i,j_{1},j_{2}}:=\Delta
x\sum_{s=j_{1}+1}^{j_{2}}\sum_{l=-P+1}^{P}a^{i,s}_{P,l}S(U^{k,m}_{i,l})H_{x}(x_{i+l}).$

Taking into account these approximations of the integral terms, the algorithm
is finally as follows, see Figure 2.1 for the case $P=2:$

  * 1.

Compute

$\displaystyle F_{i,j}^{(0)}=F(U_{i+j}^{n}),\quad j=-P+1,\ldots,P;$
$\displaystyle\widetilde{I}_{i,j-1,j}^{(0)}=\Delta
x\sum_{l=-P+1}^{P}a^{i,j}_{P,l}S(U^{n}_{i+l})H_{x}(x_{i+l}),\quad
j=-P+2,\dots,P;$ $\displaystyle\widetilde{I}_{i,-P+1,-P+1}^{(0)}=0;$
$\displaystyle\widetilde{I}_{i,-P+1,j}^{(0)}=\sum_{s=-P+2}^{j}\widetilde{I}_{i,s-1,s}^{(0)},\quad
j=-P+2,\dots,P.$

  * 2.

For $k=1,\ldots,2P-1:$

    * (a)

Compute for all $j=-P+1,\ldots,P$

$U_{i,j}^{(k)}=-\mathcal{A}_{P}^{1,j}\Bigl{(}F_{i,{*}}^{(k-1)},\Delta
x\Bigr{)}+\mathcal{A}_{P}^{1,j}\Bigl{(}I_{i,-P+1,{*}}^{(k-1)},\Delta
x\Bigr{)}.$

    * (b)

Compute for all $j,r=-P+1,\dots,P$

$U_{i,j}^{k,n+r}=U_{i+j}^{n}+\sum_{m=1}^{k}\frac{(r\Delta
t)^{m}}{m!}U_{i,j}^{(m)}.$

    * (c)

Compute for all $j,r=-P+1,\ldots,P$

$F_{i,j}^{k,n+r}=F\Bigl{(}U_{i,j}^{k,n+r}\Bigr{)},$

    * (d)

Compute for all $r=-P+1,\dots,P$, $j=-P+2,\dots,P$

$\widetilde{I}^{k,n+r}_{i,j-1,j}=\Delta
x\sum_{l=-P+1}^{P}a^{i,j}_{P,l}S(U^{k,n+r}_{i,l})H_{x}(x_{i+l}).$

    * (e)

Compute for all $j=-P+2,\dots,P$

$\widetilde{I}^{(k)}_{i,j-1,j}=A_{P}^{k,0}\Bigl{(}\widetilde{I}_{i,j-1,j}^{k,{*}},\Delta
t\Bigr{)}.$

    * (f)

Compute

$\displaystyle F_{i,j}^{(k)}=A_{P}^{k,0}\Bigl{(}F_{i,j}^{k,{*}},\Delta
t\Bigr{)},\quad j=-P+1,\ldots,P;$
$\displaystyle\widetilde{I}_{i,-P+1,-P+1}^{(k)}=0;$
$\displaystyle\widetilde{I}_{i,-P+1,j}^{(k)}=\sum_{s=-P+2}^{j}\widetilde{I}^{(k)}_{i,s-1,s}\quad
j=-P+2,\ldots,P.$

Once the algorithm has been executed, the integrals already computed can be
used to approximate the source term as follows:

  * 1.

For $k=1,\dots,2P$ define

$\widetilde{\mathcal{I}}_{i,j}^{(k-1)}=\begin{cases}\widetilde{I}^{(k-1)}_{i-1,j,j+1}&\text{if
$j=-P+1,\ldots,0$;}\\\ &\\\ \widetilde{I}^{(k-1)}_{i,j-1,j}&\text{if
$j=1,\dots,P$. }\end{cases}$

  * 2.

Compute

$\widetilde{S}_{i}^{P}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,1/2}_{P}(\widetilde{\mathcal{I}}_{i,{*}}^{(k-1)},\Delta x).$

Observe that the first $P$ integral terms $\mathcal{I}_{i,j}^{(k-1)}$
appearing in the expression of the numerical source term (3.11) are
approximated with the values $I_{i-1,j,j+1}^{(k-1)}$, used to compute the flux
at the intercell $i-1/2$, and the $P$ last ones by $I_{i,j-1,j}^{(k-1)}$, used
to compute the flux at the intercell $i+1/2$.

The final expression of the numerical method is then

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Big{(}F_{i-\frac{1}{2}}^{P}-{F}_{i+\frac{1}{2}}^{P}+\widetilde{S}_{i}^{P}\Big{)},$
(3.12)

where ${F}_{i+\frac{1}{2}}^{P}$ is given by (3.6).

###  3.2 CAT2 for system of balance laws

Let us illustrate the above numerical method in the easiest case $P=1$. In
this case, the quadrature formula used to compute integrals in intervals of
length $\Delta x$ is the trapezoidal rule:

$\int_{x_{i}}^{x_{i+1}}f(x)\,dx\approx\frac{\Delta
x}{2}\Bigl{(}f(x_{i})+f(x_{i+1})\Bigr{)}.$

The numerical method is then as follows: for every $i$

  * 1.

Compute

$U_{i,j}^{(1)}=-\frac{1}{\Delta
x}\left(F(U_{i+1}^{n})-F(U_{i}^{n})\right)+\frac{1}{2}\left(S(U_{i}^{n})H_{x}(x_{i})+S(U_{i+1}^{n})H_{x}(x_{i+1})\right),\quad
j=0,1.$

Observe that, there is no dependence on $j$ in the right term because the
derivative are computed in the same way.

  * 2.

Compute

$U_{i,j}^{1,n+1}=U_{i+j}^{n}+\Delta t\;U_{i,j}^{(1)},\quad j=0,1.$

Then, define

${F}_{i+\frac{1}{2}}^{1}:=\frac{1}{4}\left(F(U_{i}^{n})+F(U_{i+1}^{n})+F(U_{i,0}^{1,n+1})+F(U_{i,1}^{1,n+1})\right)$
(3.13) $\displaystyle\widetilde{S}_{i}^{1}$ $\displaystyle:=$
$\displaystyle\frac{\Delta
x}{8}\Bigl{(}(S(U_{i-1}^{n})+S(U_{i-1,0}^{1,n+1}))H_{x}(x_{i-1})+(S(U_{i}^{n})+S(U_{i-1,1}^{1,n+1}))H_{x}(x_{i})$
$\displaystyle+(S(U_{i}^{n})+S(U_{i,0}^{1,n+1}))H_{x}(x_{i})+(S(U_{i+1}^{n})+S(U_{i,1}^{1,n+1}))H_{x}(x_{i+1})\Bigr{)}.$

The numerical method is then (3.12) with $P=1.$

######  Remark 3.1.

The CAT procedure applied to systems of balance law introduces spurious
oscillations near discontinuity points. In Section 5 we introduce the adaptive
version of those methods, that prevent formation of spurious oscillations.

##  4 Well-balanced CAT2P for systems of balance laws

###  4.1 WBCAT2P methods

The goal of this section is to derive a well-balanced version of the CAT$2P$
methods introduced in the previous section. The idea is as follows: let us
suppose that the initial condition is given by

$U(x,0)=U^{*}(x),$

were $U^{*}$ is a stationary solution of (1.1). Let us introduce then the
function $\tilde{\mathcal{F}}$ given by

$\displaystyle\tilde{\mathcal{F}}(U)(x,t)=$
$\displaystyle\mathcal{F}(U)(x,t))-\mathcal{F}(U^{*})(x)=$ $\displaystyle=$
$\displaystyle
F(U(x,t))-F(U^{*}(x))-\int_{-\infty}^{x}\Bigl{(}S(U(\sigma,t))-S(U^{*}(\sigma))\Bigr{)}H_{x}(\sigma)d\sigma.$
(4.1)

Hence, observing that

$\tilde{\mathcal{F}}(U)_{x}=F(U)_{x}-F(U^{*})_{x}-(S(U)-S(U^{*}))H_{x}=F(U)_{x}-S(U)H_{x},$

the system of balance laws (1.1) can be formally written in the form

$U_{t}+\tilde{\mathcal{F}}(U)_{x}=0.$ (4.2)

Obviously $\tilde{\mathcal{F}}(U^{*})=0$, therefore a numerical method based
on the discretization of this conservative form is expected to preserve
$U^{*}$ exactly.

In practice, this strategy is applied as follows: once the approximation
$U_{i}^{n}$ has been obtained, we consider the stationary solution $U^{*}_{i}$
that satisfies

$U^{*}_{i}(x_{i})=U_{i}^{n},$

i.e. $U^{*}_{i}$ solves the Cauchy problem

$\left\\{\begin{array}[]{l}F(U)_{x}=S(U)H_{x}\\\
U(x_{i})=U_{i}^{n}.\end{array}\right.$ (4.3)

Let us assume for simplicity that this Cauchy problem has a unique solution
that is explicitly known. Then, the system of balance laws is rewritten in the
form

$U_{t}+\tilde{\mathcal{F}}_{i}(U)_{x}=0.$ (4.4)

where

$\displaystyle\tilde{\mathcal{F}}_{i}(U)(x,t)=$
$\displaystyle\mathcal{F}(U)(x,t))-\mathcal{F}(U^{*}_{i})(x)=$
$\displaystyle=$ $\displaystyle
F(U(x,t))-F(U^{*}_{i}(x))-\int_{-\infty}^{x}\Bigl{(}S(U(\sigma,t))-S(U^{*}_{i}(\sigma))\Bigr{)}H_{x}(\sigma)\,d\sigma$
(4.5)

and the CAT$2P$ method is then applied:

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Big{(}\mathfrak{\tilde{F}}_{i,i-\frac{1}{2}}^{P}-\mathfrak{\tilde{F}}_{i,i+\frac{1}{2}}^{P}\Big{)},$
(4.6)

where,

$\displaystyle\mathfrak{\tilde{F}}_{i,i+\frac{1}{2}}^{P}$
$\displaystyle=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A_{P}^{0,1/2}\Big{(}\tilde{\mathcal{F}}_{i;i,{*}}^{(k-1)}\Big{)},$
(4.7) $\displaystyle\mathfrak{\tilde{F}}_{i,i-\frac{1}{2}}^{P}$
$\displaystyle=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A_{P}^{0,1/2}\Big{(}\tilde{\mathcal{F}}_{i;i-1,{*}}^{(k-1)}\Big{)}.$
(4.8)

Here $\tilde{\mathcal{F}}^{(k)}_{i;l,j}$ is an approximation of
$\partial_{t}^{(k)}\tilde{\mathcal{F}}_{i}(U)(x_{l+j},t_{n}).$ Observe that,
in this case, two numerical fluxes have to be computed at every inter-cell
$x_{i+1/2}$, $\mathfrak{\tilde{F}}_{i,i+\frac{1}{2}}^{P}$ and
$\mathfrak{\tilde{F}}_{i+1,i+\frac{1}{2}}^{P}$, whose computation are based
respectively on the stationary solutions $U^{*}_{i}$ (that satisfies
$U^{*}_{i}(x_{i})=U_{i}^{n}$) and $U^{*}_{i+1}$ (that satisfies
$U^{*}_{i+1}(x_{i+1})=U_{i+1}^{n}$).

Observe that, if the initial condition $U_{0}$ is a stationary solution, then
at time $t=0$, $U^{*}_{i}=U_{0}$ for all $i$, so that

$\tilde{\mathcal{F}}_{i}(U_{0})=\mathcal{F}(U_{0})-\mathcal{F}(U^{*}_{i})=\mathcal{F}(U_{0})-\mathcal{F}(U_{0})=0,\quad\forall
i,$

and the numerical method is expected to preserve the initial condition.

The algorithm is then as follows: for every $i$

  * 1.

Compute the solution $U^{*}_{i}(x)$ of the Cauchy problem (4.3).

  * 2.

Compute

$\displaystyle F_{i;i,j}^{(0)}=F(U_{i+j}^{n})-F(U^{*}_{i}(x_{i+j})),\quad
j=-P+1,\ldots,P;$ $\displaystyle
F_{i;i-1,j}^{(0)}=F(U_{i-1+j}^{n})-F(U^{*}_{i}(x_{i-1+j})),\quad
j=-P+1,\ldots,P;$ $\displaystyle\widetilde{I}_{i;i,j-1,j}^{(0)}=\Delta
x\sum_{l=-P+1}^{P}a^{i,j}_{P,l}\left(S(U^{n}_{i+l})-S(U^{*}_{i}(x_{i+l}))\right)H_{x}(x_{i+l}),\quad
j=-P+2,\dots,P;$ $\displaystyle\widetilde{I}_{i;i,-P+1,-P+1}^{(0)}=0;$
$\displaystyle\widetilde{I}_{i;i,-P+1,j}^{(0)}=\sum_{s=-P+2}^{j}\widetilde{I}_{i;i,s-1,s}^{(0)},\quad
j=-P+3,\dots,P;$ $\displaystyle\widetilde{I}_{i;i-1,j-1,j}^{(0)}=\Delta
x\sum_{l=-P+1}^{P}a^{i-1,j}_{P,l}\left(S(U^{n}_{i-1+l})-S(U^{*}_{i}(x_{i-1+l}))\right)H_{x}(x_{i-1+l}),\quad
j=-P+3,\dots,P;$ $\displaystyle\widetilde{I}_{i;i-1,-P+1,-P+1}^{(0)}=0;$
$\displaystyle\widetilde{I}_{i;i-1,-P+1,j}^{(0)}=\sum_{s=-P+2}^{j}\widetilde{I}_{i;i-1,s-1,s}^{(0)},\quad
j=-P+2,\dots,P.$

  * 3.

For $k=1,\ldots,2P-1:$

    * (a)

Compute for all $j=-P+1,\ldots,P$

$\displaystyle U_{i;i,j}^{(k)}=-A_{P}^{1,j}\Bigl{(}F_{i;i,{*}}^{(k-1)},\Delta
x\Bigr{)}+A_{P}^{1,j}\Bigl{(}I_{i;i,-P+1,{*}}^{(k-1)},\Delta x\Bigr{)};$
$\displaystyle
U_{i;i-1,j}^{(k)}=-A_{P}^{1,j}\Bigl{(}F_{i;i-1,{*}}^{(k-1)},\Delta
x\Bigr{)}+A_{P}^{1,j}\Bigl{(}I_{i;i-1,-P+1,{*}}^{(k-1)},\Delta x\Bigr{)}.$

    * (b)

Compute for all $j,r=-P+1,\dots,P$

$\displaystyle U_{i;i,j}^{k,n+r}=U_{i+j}^{n}+\sum_{m=1}^{k}\frac{(\Delta
t)^{m}}{m!}U_{i;i,j}^{(m)},$ $\displaystyle
U_{i;i-1,j}^{k,n+r}=U_{i+j-1}^{n}+\sum_{m=1}^{k}\frac{(\Delta
t)^{m}}{m!}U_{i;i-1,j}^{(m)}.$

    * (c)

Compute for all $j,r=-P+1,\ldots,P$

$F_{i;i,j}^{k,n+r}=F\Bigl{(}U_{i;i,j}^{k,n+r}\Bigr{)},\quad
F_{i;i-1,j}^{k,n+r}=F\Bigl{(}U_{i;i-1,j}^{k,n+r}\Bigr{)}.$

    * (d)

Compute for all $r=-P+1,\dots,P$, $j=-P+2,\dots,P$

$\displaystyle\widetilde{I}^{k,n+r}_{i;i,j-1,j}=\Delta
x\sum_{l=-P+1}^{P}a^{i,j}_{P,l}S(U^{k,n+r}_{i;i,l})H_{x}(x_{i+l}),$
$\displaystyle\widetilde{I}^{k,n+r}_{i;i-1,j-1,j}=\Delta
x\sum_{l=-P+1}^{P}a^{i-1,j}_{P,l}S(U^{k,n+r}_{i;i-1,l})H_{x}(x_{i-1+l}).$

    * (e)

Compute for all $j=-P+2,\dots,P$

$\widetilde{I}^{(k)}_{i;i,j-1,j}=A_{P}^{k,0}\Bigl{(}\widetilde{I}_{i;i,j-1,j}^{k,{*}},\Delta
t\Bigr{)},\quad\widetilde{I}^{(k)}_{i;i-1,j-1,j}=A_{P}^{k,0}\Bigl{(}\widetilde{I}_{i;i-1,j-1,j}^{k,{*}},\Delta
t\Bigr{)}.$

    * (f)

Compute

$\displaystyle F_{i;i,j}^{(k)}=A_{P}^{k,0}\Bigl{(}F_{i;i,j}^{k,{*}},\Delta
t\Bigr{)},\quad j=-P+1,\ldots,P;$
$\displaystyle\widetilde{I}_{i;i,-P+1,-P+1}^{(k)}=0;$
$\displaystyle\widetilde{I}_{i;i,-P+1,j}^{(k)}=\sum_{s=-P+2}^{j}\widetilde{I}^{(k)}_{i;i,s-1,s}\quad
j=-P+2,\ldots,P;$ $\displaystyle
F_{i;i-1,j}^{(k)}=A_{P}^{k,0}\Bigl{(}F_{i;i-1,j}^{k,{*}},\Delta
t\Bigr{)},\quad j=-P+1,\ldots,P;$
$\displaystyle\widetilde{I}_{i;i-1,-P+1,-P+1}^{(k)}=0;$
$\displaystyle\widetilde{I}_{i;i-1,-P+1,j}^{(k)}=\sum_{s=-P+2}^{j}\widetilde{I}^{(k)}_{i;i-1,s-1,s}\quad
j=-P+2,\ldots,P.$

Once the algorithm has been executed, the integrals already computed can be
used to approximate the source term as follows:

  * 1.

For $k=1,\dots,2P$ define

$\widetilde{\mathcal{I}}_{i,j}^{(k-1)}=\begin{cases}\widetilde{I}^{(k-1)}_{i;i-1,j,j+1}&\text{if
$j=-P+1,\ldots,0$;}\\\ &\\\ \widetilde{I}^{(k-1)}_{i;i,j-1,j}&\text{if
$j=1,\dots,P$. }\end{cases}$

  * 2.

Compute

$\widetilde{S}_{i}^{P}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,1/2}_{P}(\widetilde{\mathcal{I}}_{i,{*}}^{(k-1)},\Delta x).$
(4.9)

The final expression of the numerical method is then

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Big{(}F_{i;i-\frac{1}{2}}^{P}-{F}_{i;i+\frac{1}{2}}^{P}+\widetilde{S}_{i}^{P}\Big{)},$
(4.10)

where ${F}_{i,i\pm\frac{1}{2}}^{P}$ are given by

$\displaystyle F^{P}_{i;i+1/2}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,1/2}_{P}(F_{i;i,{*}}^{(k-1)},\Delta x),$ (4.11)
$\displaystyle F^{P}_{i;i-1/2}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,1/2}_{P}(F_{i;i-1,{*}}^{(k-1)},\Delta x).$ (4.12)

######  Remark 4.1.

Apparently scheme 4.10 is not exactly conservative because the fluxes at the
two sides of a cell edge may be different, since they depend on the local
reconstruction. However, the schemes are exactly well balanced (see Theorem
4.1 below), and as a consequence, in case the source is identically zero, the
numerical fluxes at the two sides of a cell edges are equal, and the scheme is
conservative.

######  Remark 4.2.

Observe that this algorithm can be used to update the solution at the point
$x_{i}$ at time $t_{n}$ only if the Cauchy problem (4.3) has a solution that
is defined in the cells of the stencils $S^{P}_{i\pm 1/2}$ whose analytic
expression is known. Therefore:

  * 1.

If (4.3) has no solution, the CAT$2P$ method will be used instead. Please note
that this choice does not spoil the well-balanced character of the numerical
method: in this case, the cell values in the stencil cannot be the point
values of a stationary solution (otherwise there would be at least one
solution of (4.3)) and thus there is no local equilibrium to preserve.

  * 2.

If (4.3) has more than one solution, a criterion is needed to select one of
them: this is the case for the shallow water system that will be discussed in
Section 4.4.

  * 3.

If (4.3) has a solution defined in the stencils but it is not possible to find
its expression by analytic procedures, it is possible to apply an ODE solver
to approximate it, like it has been done in [9] for finite-volume methods. In
all the problems considered in Section 7 the analytic expression of the
stationary solutions is available either in explicit or implicit form.

###  4.2 Well-balanced property

Numerical method (4.10) is fully well-balanced in the following sense:

######  Theorem 4.1.

Let $U^{*}$ be a continuous stationary solution of (1.1). Then, if numerical
method (4.10) is applied to the initial condition

$U^{0}_{i}=U^{*}(x_{i}),\quad\forall i,$

then

$U_{i}^{n}=U_{i}^{0},\quad\forall i,n.$

Proof: The proof is based on induction on the order of the proximate time
derivative. Observe first that $U^{*}$ solves any Cauchy problem (4.3) for
$n=0$. Therefore, at the first step the solution of (4.3) is given by

$U^{*}_{i}=U^{*},\quad\forall i.$

Therefore, for every $i$:

$\displaystyle F_{i;i,j}^{(0)}=F_{i;i-1,j}^{(0)}=0,\quad j=-P+1,\ldots,P;$
$\displaystyle\widetilde{I}_{i;i,j-1,j}^{(0)}=\widetilde{I}_{i;i-1,j-1,j}^{(0)}=0,\quad
j=-P+2,\ldots,P;$
$\displaystyle\widetilde{I}_{i;i,-P+1,j}^{(0)}=\widetilde{I}_{i;i-1,-P+1,j}^{(0)}=0,\quad
j=-P+1,\dots,P;$

and thus

$U_{i;i,j}^{(1)}=U_{i;i-1,j}^{(1)}=0,\quad j=-P+1,\ldots,P.$

As a consequence:

$\displaystyle U_{i;i,j}^{1,r}=U_{i+j}^{0},\quad
F_{i;i,j}^{1,r}=F\Bigl{(}U_{i+j}^{0}\Bigr{)},\quad j,r=-P+1,\dots,P;$
$\displaystyle U_{i;i-1,j}^{1,r}=U_{i-1+j}^{0},\quad
F_{i;i-1,j}^{1,r}=F\Bigl{(}U_{i-1+j}^{0}\Bigr{)},\quad j,r=-P+1,\dots,P$
$\displaystyle\widetilde{I}^{1,r}_{i;i,j-1,j}=\Delta
x\sum_{l=-P+1}^{P}a^{i,j}_{P,l}S(U^{0}_{i+l})H_{x}(x_{i+l}),$
$\displaystyle\widetilde{I}^{1,r}_{i;i-1,j-1,j}=\Delta
x\sum_{l=-P+1}^{P}a^{i-1,j}_{P,l}S(U^{0}_{i-1+l})H_{x}(x_{i-1+l}).$

Notice that the values of all these quantities do not depend on $r$.
Therefore, when numerical differentiation in time is applied we obtain:

$\displaystyle F_{i;i,j}^{(1)}=F_{i;i-1,j}^{(1)}=0,\quad j=-P+1,\ldots,P;$
$\displaystyle\widetilde{I}_{i;i,j-1,j}^{(1)}=\widetilde{I}_{i;i-1,j-1,j}^{(1)}=0,\quad
j=-P+2,\ldots,P;$
$\displaystyle\widetilde{I}_{i;i,-P+1,j}^{(1)}=\widetilde{I}_{i;i-1,-P+1,j}^{(1)}=0,\quad
j=-P+1,\ldots,P.$

Therefore

$U_{i;i,j}^{(2)}=U_{i;i-1,j}^{(2)}=0,\quad j=-P+1,\ldots,P.$

Repeating the reasoning we obtain

$\displaystyle F_{i;i,j}^{(k)}=F_{i;i-1,j}^{(k)}=0,\quad j=-P+1,\ldots,P,\quad
k=0,\dots,2P;$
$\displaystyle\widetilde{I}_{i;i,j-1,j}^{(k)}=\widetilde{I}_{i;i-1,j-1,j}^{(k)}=0,\quad
j=-P+2,\ldots,P,\quad k=0,\dots,2P;$
$\displaystyle\widetilde{I}_{i;i,-P+1,j}^{(1)}=\widetilde{I}_{i;i-1,-P+1,j}^{(1)}=0,\quad
j=-P+1,\ldots,P,\quad k=0,\dots,2P.$

Therefore

$F^{P}_{i;i+1/2}=F^{P}_{i,i-1/2}=S^{P}_{i+1/2}=0,$

and we obtain

$U^{1}_{i}=U^{0}_{i},\quad\forall i,$

as we wanted to prove.

###  4.3 WBCAT2 for system of balance laws

Let us illustrate again this numerical method in the case $P=1$:

For every $i$

  * 1.

Compute the solution $U^{*}_{i}$ of the Cauchy problem (4.3).

  * 2.

Compute

$\displaystyle U_{i;i,j}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{1}{\Delta
x}\Bigl{(}F(U_{i+1}^{n})-F(U_{i}^{n})-F(U_{i}^{*}(x_{i+1}))+F(U_{i}^{*}(x_{i}))\Bigr{)}$
$\displaystyle+\frac{1}{2}\Bigl{(}(S(U_{i}^{n})-S(U_{i}^{*}(x_{i})))H_{x}(x_{i})+(S(U_{i+1}^{n})-S(U_{i}^{*}(x_{i+1})))H_{x}(x_{i+1})\Bigr{)},\quad
j=0,1;$ $\displaystyle U_{i;i-1,j}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{1}{\Delta
x}\Bigl{(}F(U_{i}^{n})-F(U_{i-1}^{n})-F(U_{i}^{*}(x_{i}))+F(U_{i}^{*}(x_{i-1}))\Bigr{)}$
$\displaystyle+\frac{1}{2}\Bigl{(}(S(U_{i-1}^{n})-S(U_{i}^{*}(x_{i-1})))H_{x}(x_{i-1})+(S(U_{i}^{n})-S(U_{i}^{*}(x_{i})))H_{x}(x_{i})\Bigr{)},\quad
j=0,1.$

  * 3.

Compute

$\displaystyle U_{i;i,j}^{1,n+1}=U_{i+j}^{n}+\Delta t\;U_{i;i,j}^{(1)},\quad
j=0,1,$ $\displaystyle U_{i;i-1,j}^{1,n+1}=U_{i+j}^{n}+\Delta
t\;U_{i-1;i,j}^{(1)},\quad j=0,1.$

  * 4.

Define

$\displaystyle{F}_{i,i+\frac{1}{2}}^{1}:=\frac{1}{4}\left(F(U_{i}^{n})+F(U_{i+1}^{n})+F(U_{i;i,0}^{1,n+1})+F(U_{i;i,1}^{1,n+1})-2F(U^{*}_{i}(x_{i}))+\right.$
(4.13) $\displaystyle\qquad\qquad\quad-\left.2F(U^{*}_{i}(x_{i+1}))\right),$
$\displaystyle{F}_{i,i-\frac{1}{2}}^{1}:=\frac{1}{4}\left(F(U_{i-1}^{n})+F(U_{i}^{n})+F(U_{i;i-1,0}^{1,n+1})+F(U_{i;i-1,1}^{1,n+1})-2F(U^{*}_{i}(x_{i-1}))\right.$
(4.14) $\displaystyle\qquad\qquad\quad-\left.-2F(U^{*}_{i}(x_{i}))\right),$
$\displaystyle\widetilde{S}_{i}^{1}:=\frac{\Delta
x}{8}\Bigl{(}(S(U_{i-1}^{n})+S(U_{i-1,0}^{1,n+1})-2S(U^{*}_{i}(x_{i-1})))H_{x}(x_{i-1})$
(4.15)
$\displaystyle\qquad\qquad+(S(U_{i}^{n})+S(U_{i;i-1,1}^{1,n+1})-2S(U^{*}_{i}(x_{i})))H_{x}(x_{i})$
$\displaystyle\qquad\qquad+(S(U_{i}^{n})+S(U_{i;i,0}^{1,n+1})-2S(U^{*}_{i}(x_{i})))H_{x}(x_{i})$
$\displaystyle\qquad\qquad+(S(U_{i+1}^{n})+S(U_{i;i,1}^{1,n+1})-2S(U^{*}_{i}(x_{i+1})))H_{x}(x_{i+1})\Bigr{)}.$

The numerical method is then (4.10) with $P=1$.

###  4.4 CAT2P methods that preserve a family of stationary solutions

The strategy described in Subsection 4.1 can be easily adapted to obtain
schemes that only preserve a prescribed set of stationary solutions: this
would be the case if, for instance, the set to be preserved is a $k$-parameter
family of stationary solutions,

$U^{*}(x;C_{1},\dots,C_{k}),$

with $k<d$ where $d$ is the number of variables. If it is the case, instead of
looking for a solution of (4.3), one looks for a solution of the following
nonlinear system:

Find $C^{i}_{1},\dots,C^{i}_{k}$ such that:

$u_{j_{l}}^{*}(x_{i};C^{i}_{1},\dots,C^{i}_{k})=u_{i,j_{l}},\quad
l=1,\dots,k,$ (4.16)

where $u^{*}_{j}$, $u_{i,j}$ denote respectively the $j$-th component of
$U^{*}$ and $U_{i}$ and $\\{j_{1},\dots,j_{k}\\}$ is a set of $k$ indices that
is predetermined in order to have the same number of unknowns and equations in
(4.16). These indices $j_{1},\dots,j_{k}$ are chosen so that systems of
equations (4.16) have a unique solution, if it is possible. Once the problem
has been solved, the numerical fluxes and source terms are computed as in
Section 4.1 with the choice

$U^{*}_{i}(x)=U^{*}(x,C^{i}_{1},\dots,C^{i}_{k}).$

##  5 Adaptive CAT2P for systems of balance law

Despite the fact that Compact Approximate Taylor (CAT) schemes are linearly
stable in the $L^{2}$ sense under the usual CFL condition, spurious
oscillations may appear next to a discontinuity of the solution, as it happens
for the Lax-Wendroff method: see [4]. In order to get rid of them, we consider
here the shock-capturing technique introduced in [3] based on a family of
high-order smoothness indicators. The idea is as follows: once the
approximations at time $t^{n}$ have been computed, the candidate stencils to
compute $\mathfrak{F}_{i+\frac{1}{2}}$ are

$\mathcal{S}^{p}_{i+1/2}=\\{x_{i-p+1},\ldots,x_{i+p}\\},\quad p=1,\ldots,P.$

The selected stencil is the one with maximal length among those in which the
solution at time $t^{n}$ is smooth, according to some smoothness indicators
$\psi^{p}_{i+1/2}$ for $p=1,\dots,P$. If a discontinuity is detected in the
stencil $\mathcal{S}^{1}_{i+1/2}$ a robust first-order numerical method is
used. The ingredients of this strategy are described below:

###  5.1 First-order numerical method

As first-order robust scheme to be combined with CAT2P methods for balance
laws, we select the Lax-Friedrichs method applied to (3.2) which leads to the
formal expression:

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Bigl{(}F^{LF}_{i-1/2}-F^{LF}_{i+1/2}+\widetilde{S}^{LF}_{i}\Bigr{)}$ (5.1)

where

$\displaystyle F^{LF}_{i+1/2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(F(U^{n}_{i})+F(U^{n}_{i+1})\right)-\frac{\Delta
x}{2\Delta t}\left(U^{n}_{i+1}-U^{n}_{i}\right),$ (5.2)
$\displaystyle\widetilde{S}^{LF}_{i}$ $\displaystyle=$ $\displaystyle\Delta
xS(U_{i}^{n})H_{x}(x_{i}),$ (5.3)

where the mid-point rule has been used to approximate the integral
corresponding to the source term.

In the case of the WBCAT2$P$ methods, the Lax-Friedrichs method is formally
applied to (4.4) and the mid-point formula is again used to compute the
integral terms, which leads again to a method of the form (5.1) where
$F^{LF}_{i+1/2}$ is again the standard Lax-Friedrichs numerical flux (5.2) but
now

$\widetilde{S}^{LF}_{i}={F}^{LF,*}_{i;i+1/2}-{F}^{LF,*}_{i;i-1/2},$ (5.4)

where

${F}^{LF,*}_{i;i+1/2}=\frac{1}{2}\left(F(U_{i}^{*}(x_{i}))+F(U_{i}^{*}(x_{i+1}))\right)-\frac{\Delta
x}{2\Delta t}\left(U^{*}_{i}(x_{i+1})-U^{*}_{i}(x_{i})\right)$ (5.5)

and $U^{*}_{i}$ is the stationary solution that satisfies

$U^{*}_{i}(x_{i})=U_{i}^{n}.$

(5.4) is a consistent approximation of the integral of the source term, since

${F}^{LF,*}_{i;i+1/2}-{F}^{LF,*}_{i;i-1/2}\approx\Delta
xF(U^{*})_{x}(x_{i})=\Delta xS(U^{*}_{i}(x_{i}))H_{x}(x_{i})=\Delta
xS(U^{n}_{i})H_{x}(x_{i}).$

On the other hand, the well-balanced property of the method can be trivially
checked.

###  5.2 ACAT2 and WBACAT2 methods

The expression of the ACAT2 numerical method is based on a flux limiter (see
[21, 22, 32]). Its expression is as follows:

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Bigl{(}F^{*}_{i-1/2}-F^{*}_{i+1/2}+\widetilde{S}^{*}_{i}\Bigr{)}$ (5.6)

where

$\displaystyle{F}^{*}_{i\pm 1/2}$ $\displaystyle=$
$\displaystyle\varphi^{1}_{i}\,{F}^{1}_{i\pm
1/2}+(1-\varphi^{1}_{i})\,{F}^{LF}_{i\pm 1/2},$ (5.7)
$\displaystyle\widetilde{S}_{i}^{*}$ $\displaystyle=$
$\displaystyle\varphi^{1}_{i}\widetilde{S}^{1}_{i}+(1-\varphi^{1}_{i})\widetilde{S}^{LF}_{i},$
(5.8)

where ${F}^{1}_{i\pm 1/2}$ and $\widetilde{S}^{1}$ are given by (3.13) and
(3.2) respectively; $F^{LF}_{i\pm 1/2}$ and $\widetilde{S}^{LF}_{i}$ are given
by (5.2) and (5.3) respectively;

$\varphi^{1}_{i}=\min(\varphi^{1}_{i-1/2},\varphi^{1}_{i+1/2}),$

where $\varphi^{1}_{i+1/2}$ is computed by a flux limiter, and has the
property

$\varphi^{1}_{i+1/2}\approx\begin{cases}1&\text{if
$\\{U_{i-1}^{n},\dots,U_{i+2}^{n}\\}$ are 'smooth';}\\\
0&\text{otherwise.}\end{cases}$ (5.9)

For scalar problems, standard flux limiter functions $\varphi^{1}(r)$, such as
minmod, superbee, van Leer [30, 19], may be used:

$\varphi_{i+1/2}^{1}=\varphi^{1}(r_{i+1/2}),$ (5.10)

where

$r_{i+1/2}=\frac{\Delta_{\rm upw}}{\Delta_{\rm
loc}}=\left\\{\begin{array}[]{cl}\displaystyle
r^{-}_{i+1/2}:=\frac{u^{n}_{i}-u^{n}_{i-1}}{u^{n}_{i+1}-u^{n}_{i}}&\mbox{if
}a_{i+1/2}>0,\\\ \displaystyle
r^{+}_{i+1/2}:=\frac{u^{n}_{i+2}-u^{n}_{i+1}}{u^{n}_{i+1}-u^{n}_{i}}&\mbox{if
}a_{i+1/2}\leq 0;\end{array}\right.$ (5.11)

and $a_{i+1/2}$ is an estimate of the wave speed such as for instance Roe’s
intermediate speed:

$a_{i+1/2}=\begin{cases}\displaystyle\frac{f(u^{n}_{i+1})-f(u^{n}_{i})}{u^{n}_{i+1}-u^{n}_{i}}&\text{
if $|u^{n}_{i}-u^{n}_{i+1}|>{\rm tol};$}\\\
f^{\prime}(u^{n}_{i})&\text{otherwise.}\end{cases}$

where tol is a given small tolerance. An alternative that avoids the
computation of an intermediate speed was introduced in [32] and it consists on
defining

$\varphi^{1}_{i+1/2}=\min(\varphi^{1}(r_{i+1/2}^{+}),\varphi^{1}(r^{-}_{i+1/2})).$
(5.12)

This strategy can be easily extended to systems by computing the flux limiter
component by component.

The well-balanced version of the ACAT2 method has the form

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Bigl{(}F^{*}_{i;i-1/2}-F^{*}_{i;i+1/2}+\widetilde{S}^{*}_{i}\Bigr{)}$
(5.13)

where

$\displaystyle{F}^{*}_{i;i\pm 1/2}$ $\displaystyle=$
$\displaystyle\varphi^{1}_{i}\,{F}^{1}_{i;i\pm
1/2}+(1-\varphi^{1}_{i})\,{F}^{LF}_{i;i\pm 1/2},$ (5.14)
$\displaystyle\widetilde{S}_{i}^{*}$ $\displaystyle=$
$\displaystyle\varphi^{1}_{i}\widetilde{S}^{1}_{i}+(1-\varphi^{1}_{i})\widetilde{S}^{LF,*}_{i}=\varphi^{1}_{i}\widetilde{S}^{1}_{i}-(1-\varphi^{1}_{i})\left(F^{LF,*}_{i;i+\frac{1}{2}}-F^{LF,*}_{i;i-\frac{1}{2}}\right),$
(5.15)

where ${F}^{1}_{i;i\pm 1/2}$ are given by (4)-(4); $F^{LF,*}_{i;i\pm 1/2}$ is
given by (5.5); and $\widetilde{S}_{i}^{1}$ is given by (4.15).

###  5.3 Smoothness Indicators

The smoothness indicators introduced in [3] will be used here: given a set of
the point values $f_{j}$ of a function $f$ in the nodes of the stencil
$S^{p}_{i+1/2}$, $p\geq 2$ we define $\psi_{i+1/2}^{p}$ as follows: first
define the lateral weights $I_{p,L}$ and $I_{p,R}$ as

$I_{p,L}:=\sum_{j=-p+1}^{-1}(f_{i+1+j}-f_{i+j})^{2}+\varepsilon,\quad
I_{p,R}:=\sum_{j=1}^{p-1}(f_{i+1+j}-f_{i+j})^{2}+\varepsilon,$ (5.16)

where $\varepsilon$ is a small quantity that is added to prevent the lateral
weights to vanish when the function is constant. Next, compute the half
harmonic mean

$I_{p}=\frac{I_{p,L}I_{p,R}}{I_{p,L}+I_{p,R}}.$ (5.17)

Figure 5.1: Left and right part to compute the high order smoothness
indicators where $P=2.$

Finally, we define the high order smoothness indicator of the stencil $S_{p}$
by

$\psi_{i+1/2}^{p}:=\left(\frac{I_{p}}{I_{p}+\tau_{p}}\right),$ (5.18)

where

$\displaystyle\tau_{p}:=$
$\displaystyle\left(\Delta^{2p-1}_{i-p+1}f\right)^{2}.$ (5.19)

The script $\Delta^{2p-1}_{i-p+1}f$ represents the undivided difference of
$\\{f_{i-p+1},\dots,f_{i+p}\\}$:

$\displaystyle\Delta^{2p-1}_{i-p+1}f=$
$\displaystyle(2p-1)!\sum^{p}_{j=-p+1}\,\gamma^{2p-1,1/2}_{p,j}\,f^{n}_{i+j}.$
(5.20)

These indicators are such that

$\psi^{p}_{i+1/2}\approx\left\\{\begin{array}[]{cl}1&\mbox{ if $\\{f_{j}\\}$
are 'smooth' in $\mathcal{S}^{p}_{i}$;}\\\
0&\mbox{otherwise;}\end{array}\right.$ (5.21)

see [3] for a precise statement of this property and its proof.

###  5.4 ACAT2P and WBACAT2P methods

Using all these ingredients, the final expression of the ACAT2$P$ method is as
follows:

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Bigl{(}{F}^{\mathcal{A}_{i}}_{i;i-\frac{1}{2}}-{F}^{\mathcal{A}_{i}}_{i;i+\frac{1}{2}}+\widetilde{S}_{i}^{\mathcal{A}_{i}}\Bigr{)},$
(5.22)

where

${F}^{\mathcal{A}_{i}}_{i;i\pm 1/2}=\begin{cases}{F}^{*}_{i\pm 1/2}&\text{if
$\mathcal{A}_{i}=\emptyset$;}\\\ {F}^{p_{s}}_{i\pm 1/2}&\text{otherwise, where
$p_{s}=\max(\mathcal{A}_{i})$;}\end{cases}$ (5.23)

and

$\widetilde{S}^{\mathcal{A}_{i}}_{i}=\begin{cases}\widetilde{S}^{*}_{i}&\text{if
$\mathcal{A}_{i}=\emptyset$;}\\\ \widetilde{S}^{p_{s}}_{i}&\text{otherwise,
where $p_{s}=\max(\mathcal{A}_{i}).$}\end{cases}$ (5.24)

Here, $\mathcal{A}_{i}$ is the set of indices given by

$\mathcal{A}_{i}=\\{p\in\\{2,\dots,P\\}\ s.t.\ \psi^{p}_{i-1/2}\approx 1\text{
and }\psi^{p}_{i+1/2}\approx 1\\};$ (5.25)

$F^{*}_{i+1/2}$ and $\widetilde{S}_{i}^{*}$ are the ACAT2 numerical flux and
source terms given by (5.7), (5.8); $F^{p_{s}}_{i+1/2}$ and
$\widetilde{S}_{i}^{p_{s}}$ are the ACAT2$p_{s}$ numerical fluxes and source
terms defined in (3.5), (3.11). Throughout the paper we assume $\psi\approx 1$
if $\psi\geq 0.9$.

######  Remark 5.1.

Notice that in Equation 5.25 index $p$ starts from 2, since it is not possible
to determine the smoothness of the data in the two-point stencil
$S^{1}_{i+1/2}$: see [3].

Analogously, the final expression of the WBACAT2$P$ method is as follows:

$U_{i}^{n+1}=U_{i}^{n}+\frac{\Delta t}{\Delta
x}\Bigl{(}{F}^{\mathcal{A}_{i}}_{i;i-\frac{1}{2}}-{F}^{\mathcal{A}_{i}}_{i;i+\frac{1}{2}}+\widetilde{S}_{i}^{\mathcal{A}_{i}}\Bigr{)},$
(5.26)

where

${F}^{\mathcal{A}_{i}}_{i;i+1/2}=\begin{cases}{F}^{*}_{i;i+1/2}&\text{if
$\mathcal{A}_{i}=\emptyset$;}\\\ {F}^{p_{s}}_{i;i+1/2}&\text{otherwise, where
$p_{s}=\max(\mathcal{A}_{i});$}\end{cases}$ (5.27)

and

$\widetilde{S}^{\mathcal{A}_{i}}_{i}=\begin{cases}\widetilde{S}^{*}_{i}&\text{if
$\mathcal{A}_{i}=\emptyset$;}\\\ \widetilde{S}^{p_{s}}_{i}&\text{otherwise,
where $p_{s}=\max(\mathcal{A}_{i});$}\end{cases}$ (5.28)

where ${F}^{*}_{i;i\pm 1/2}$ and $\widetilde{S}^{*}_{i}$ are the WBACAT2
numerical fluxes and source terms given by (5.14), (5.15);
${F}^{p_{s}}_{i;i+1/2}$ and $\widetilde{S}^{p_{s}}_{i}$ are WBCAT$2p_{s}$
numerical fluxes and source term given by (4.12), (4.11), (4.9) with
$P=p_{s}$.

##  6 Two-dimensional problems

In this section we will focus on the extension of ACAT methods to non-linear
two-dimensional systems of hyperbolic balance laws

${U}_{t}+{F}({U})_{x}+{G}({U})_{y}=S_{1}(U)H_{x}+S_{2}(U)H_{y}.$ (6.1)

As we did for the 1D systems of balance laws (1.1), let us introduce the
functions $\mathcal{F}$ and $\mathcal{G}$ as:

$\displaystyle\mathcal{F}(U)((x,y),t)$
$\displaystyle=F(U(x,y,t))-\int_{-\infty}^{x}S_{1}(U)H_{\sigma}d\sigma;$ (6.2)
$\displaystyle\mathcal{G}(U)((x,y),t)$
$\displaystyle=G(U(x,y,t))-\int_{-\infty}^{y}S_{2}(U)H_{\tau}d\tau,$ (6.3)

assuming that the integrals are finite. Then, the identities

$\mathcal{F}(U)_{x}=F(U(x,y,t))_{x}-S_{1}(U)H_{x}$

and

$\mathcal{G}(U)_{y}=G(U(x,y,t))_{y}-S_{2}(U)H_{y}$

allows one to write the 2D system of balance laws (6.1) in the equivalent
conservative form

$U_{t}+\mathcal{F}(U)_{x}+\mathcal{G}(U)_{y}=0.$ (6.4)

The idea is now extend the ACAT$2P$ schemes to 2D systems written in the form
(6.4).

The following multi-index notation will be used:

$\mathbf{i}=(i_{1},i_{2})\in\mathbb{Z}\times\mathbb{Z},$

and

$\textbf{0}=(0,0)\quad\mathbf{1}=(1,1)\quad\mathbf{1/2}=(1/2,1/2),\quad\mathbf{e}_{1}=(1,0),\quad\mathbf{e}_{2}=(0,1).$

We consider Cartesian meshes with nodes

$\mathbf{x}_{\mathbf{i}}=(i_{1}\Delta x,i_{2}\Delta y).$

Using this notation, CAT$2P$ methods can be extended as follows:

$U_{\mathbf{i}}^{n+1}=U_{\mathbf{i}}^{n}+\frac{\Delta t}{\Delta
x}\left[{F}_{\mathbf{i}-\frac{1}{2}\mathbf{e}_{1}}^{P}-{F}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{1}}^{P}+\widetilde{S}^{P}_{1,\mathbf{i}}\right]+\frac{\Delta
t}{\Delta
y}\left[{G}_{\mathbf{i}-\frac{1}{2}\mathbf{e}_{2}}^{P}-{G}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{2}}^{P}+\widetilde{S}^{P}_{2,\mathbf{i}}\right]$
(6.5)

where the numerical fluxes ${F}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{1}}^{P}$,
${G}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{2}}^{P}$ and the source terms
$\widetilde{S}^{P}_{j,\mathbf{i}}$, $j=1,2$ are computed using the values of
the numerical solution $U_{\mathbf{i}}^{n}$ in the $P^{2}-$point stencil
centered at $\mathbf{x}_{\mathbf{i}+\mathbf{1/2}}=((i_{1}+\frac{1}{2})\Delta
x,(i_{2}+\frac{1}{2})\Delta y)$

$\mathcal{S}^{P}_{\mathbf{i}+\mathbf{1/2}}=\\{\mathbf{x}_{\mathbf{i}+\mathbf{j}},\quad\mathbf{j}\in\mathfrak{I}_{P}\\},$

where

$\mathfrak{I}_{P}=\\{\mathbf{j}=(j_{1},j_{2})\in\mathbb{Z}\times\mathbb{Z},\quad-P+1\leq
j_{k}\leq P,\quad k=1,2\\}.$

See Figure 6.1 for an example.

Figure 6.1: Stencil $S_{2}$ centered in $\mathbf{x}_{\mathbf{1/2}}=(0.5\Delta
x,0.5\Delta y$). The rectangular black grid represents the space dimension
while the blue line means the time direction of Taylor approximation.

For the sake of simplicity let us introduce index sets which will be used in
the approximation of integrals of the source term

$\mathfrak{M}_{P}=\\{\mathbf{j}\in\mathfrak{I}_{P}\quad{\rm
such\;that}\quad-P+2\leq j_{1}\leq P\\}$

and

$\mathfrak{N}_{P}=\\{\mathbf{j}\in\mathfrak{I}_{P}\quad{\rm
such\;that}\quad-P+2\leq j_{2}\leq P\\}.$

The following algorithm will be used to compute the numerical fluxes
${F}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{1}}^{P}$,
${G}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{2}}^{P}$ and the source terms
$\widetilde{S}^{P}_{j,\mathbf{i}}$, $j=1,2$ of 2D CAT$2P$ method:

  * 1.

Define

$\displaystyle
F_{\mathbf{i},\mathbf{j}}^{(0)}:=F(U_{\mathbf{i}+\mathbf{j}}^{n}),\quad\mathbf{j}\in\mathfrak{I}_{P};$
$\displaystyle
G_{\mathbf{i},\mathbf{j}}^{(0)}:=G(U_{\mathbf{i}+\mathbf{j}}^{n}),\quad\mathbf{j}\in\mathfrak{I}_{P};$
$\displaystyle
I_{\mathbf{i},\mathbf{j}-\mathbf{e}_{1},\mathbf{j}}^{(0)}:=\Delta
x\sum_{q=-P+1}^{P}a_{P,q}^{i_{1},j_{1}}S_{1}(U^{n}_{\mathbf{i}+q\mathbf{e}_{1}})H_{x}(x_{\mathbf{i}+q\mathbf{e}_{1}}),\quad\mathbf{j}\in\mathfrak{M}_{P};$
$\displaystyle I_{\mathbf{i},(-P+1,j_{2}),(-P+1,j_{2})}^{(0)}=0,\quad
j_{2}=-P+1,\dots,P;$ $\displaystyle
I_{\mathbf{i},(-P+1,j_{2}),\mathbf{j}}^{(0)}=\sum_{s=-P+2}^{j_{1}}I_{\mathbf{i},(s-1,j_{2}),(s,j_{2})}^{(0)}\quad\mathbf{j}\in\mathfrak{M}_{P};$
$\displaystyle
J_{\mathbf{i},\mathbf{j}-\mathbf{e}_{2},\mathbf{j}}^{(0)}:=\Delta
y\sum_{q=-P+1}^{P}a_{P,q}^{i_{2},j_{2}}S_{2}(U^{n}_{\mathbf{i}+q\mathbf{e}_{2}})H_{y}(x_{\mathbf{i}+q\mathbf{e}_{2}}),\quad\mathbf{j}\in\mathfrak{N}_{P};$
$\displaystyle J_{\mathbf{i},(j_{1},-P+1),(j_{1},-P+1)}^{(0)}=0,\quad
j_{1}=-P+1,\dots,P;$ $\displaystyle
J_{\mathbf{i},(j_{1},-P+1),\mathbf{j}}^{(0)}=\sum_{s=-P+2}^{j_{2}}J_{\mathbf{i},(j_{1},s-1),(j_{1},s)}^{(0)}\quad\mathbf{j}\in\mathfrak{N}_{P};$

  * 2.

For $k=1,\ldots,2P-1:$

    * (a)

Compute for all $\mathbf{j}\in\mathfrak{I}_{P}$

$\displaystyle U_{\mathbf{i},\mathbf{j}}^{(k)}$ $\displaystyle=$
$\displaystyle-
A_{P}^{1,j_{1}}\Bigl{(}F_{\mathbf{i},({*},j_{2})}^{(k-1)},\Delta
x\Bigr{)}+A_{P}^{1,j_{1}}\Bigl{(}I_{\mathbf{i},(-P+1,j_{2}),({*},j_{2})}^{(k-1)},\Delta
x\Bigr{)}$ $\displaystyle-
A_{P}^{1,j_{2}}\Bigl{(}G_{\mathbf{i},(j_{1},{*})}^{(k-1)},\Delta
y\Bigr{)}+A_{P}^{1,j_{2}}\Bigl{(}J_{\mathbf{i},(j_{1},-P+1),(j_{1},{*})}^{(k-1)},\Delta
y\Bigr{)}.$

    * (b)

Compute for all $\mathbf{j}\in\mathfrak{I}_{P}$ and for all $r=-P+1,\ldots,P$

$U_{\mathbf{i},\mathbf{j}}^{k,n+r}=U_{\mathbf{i}+\mathbf{j}}^{n}+\sum_{m=1}^{k}\frac{(\Delta
t)^{m}}{m!}U_{\mathbf{i},\mathbf{j}}^{(m)}.$

    * (c)

Compute for all $\mathbf{j}\in\mathfrak{I}_{P}$ and for all $r=-P+1,\ldots,P$

$F_{\mathbf{i},\mathbf{j}}^{k,n+r}=F(U_{\mathbf{i},\mathbf{j}}^{k,n+r})\quad{\rm{and}}\quad
G_{\mathbf{i},\mathbf{j}}^{k,n+r}=G(U_{\mathbf{i},\mathbf{j}}^{k,n+r}).$

    * (d)

Compute for all $\mathbf{j}\in\mathfrak{M}_{P}$ and for all $r=-P+1,\ldots,P$

$\displaystyle I_{\mathbf{i},\mathbf{j}-\mathbf{e}_{1},\mathbf{j}}^{k,n+r}$
$\displaystyle=$ $\displaystyle\Delta
x\sum_{q=-P+1}^{P}a_{P,q}^{i_{1},j_{1}}S_{1}(U^{n}_{\mathbf{i}+q\mathbf{e}_{1}})H_{x}(x_{\mathbf{i}+q\mathbf{e}_{1}});$
$\displaystyle I_{\mathbf{i},\mathbf{j}-\mathbf{e}_{1},\mathbf{j}}^{(k)}$
$\displaystyle=$ $\displaystyle
A^{k,0}_{P}\Bigl{(}I_{\mathbf{i},\mathbf{j}-\mathbf{e}_{1},\mathbf{j}}^{k,{*}},\Delta
t\Bigr{)}.$

    * (e)

Compute for all $\mathbf{j}\in\mathfrak{N}_{P}$ and for all $r=-P+1,\ldots,P$

$\displaystyle J_{\mathbf{i},\mathbf{j}-\mathbf{e}_{2},\mathbf{j}}^{k,n+r}$
$\displaystyle=$ $\displaystyle\Delta
y\sum_{q=-P+1}^{P}a_{P,q}^{i_{2},j_{2}}S_{2}(U^{n}_{\mathbf{i}+q\mathbf{e}_{2}})H_{y}(x_{\mathbf{i}+q\mathbf{e}_{2}});$
$\displaystyle J_{\mathbf{i},\mathbf{j}-\mathbf{e}_{2},\mathbf{j}}^{(k)}$
$\displaystyle=$ $\displaystyle
A^{k,0}_{P}\Bigl{(}J_{\mathbf{i},\mathbf{j}-\mathbf{e}_{2},\mathbf{j}}^{k,{*}},\Delta
t\Bigr{)}.$

    * (f)

Compute

$\displaystyle
F_{\mathbf{i},\mathbf{j}}^{(k)}=A^{k,0}_{P}\Bigl{(}F_{\mathbf{i},\mathbf{j}}^{k,{*}},\Delta
t\Bigr{)},\quad\mathbf{j}\in\mathfrak{I}_{P};$ $\displaystyle
G_{\mathbf{i},\mathbf{j}}^{(k)}=A^{k,0}_{P}\Bigl{(}G_{\mathbf{i},\mathbf{j}}^{k,{*}},\Delta
t\Bigr{)},\quad\mathbf{j}\in\mathfrak{I}_{P};$ $\displaystyle
I_{\mathbf{i},(-P+1,j_{2}),(-P+1,j_{2})}^{(k)}=0,\quad j_{2}=-P+1,\dots,P;$
$\displaystyle
I_{\mathbf{i},(-P+1,j_{2}),\mathbf{j}}^{(k)}=\sum_{s=-P+2}^{j_{1}}I_{\mathbf{i},(s-1,j_{2}),(s,j_{2})}^{(k)}\quad\mathbf{j}\in\mathfrak{M}_{P};$
$\displaystyle J_{\mathbf{i},(j_{1},-P+1),(j_{1},-P+1)}^{(k)}=0,\quad
j_{1}=-P+1,\dots,P;$ $\displaystyle
J_{\mathbf{i},(j_{1},-P+1),\mathbf{j}}^{(k)}=\sum_{s=-P+2}^{j_{2}}I_{\mathbf{i},(j_{1},s-1),(j_{1},s)}^{(k)}\quad\mathbf{j}\in\mathfrak{N}_{P};$

  * 3.

Compute

$\displaystyle
F^{P}_{\mathbf{i}+\mathbf{1/2}\mathbf{e}_{1}}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,\frac{1}{2}}_{P}\Bigl{(}F_{\mathbf{i},({*},0)}^{(k-1)},\Delta
x\Bigr{)};$ (6.6) $\displaystyle
G^{P}_{\mathbf{i}+\mathbf{1/2}\mathbf{e}_{2}}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,\frac{1}{2}}_{P}\Bigl{(}G_{\mathbf{i},(0,{*})}^{(k-1)},\Delta
y\Bigr{)}$ (6.7)

Once the algorithm has finished, the integrals will have already been computed
and can be used to approximate the source term as follows:

  * 4.

For $k=1,\ldots,2P$ define

$\displaystyle\mathcal{I}_{\mathbf{i},j_{1}}^{(k-1)}$ $\displaystyle=$
$\displaystyle\begin{cases}I_{\mathbf{i}-\mathbf{e}_{1},j_{1}\mathbf{e}_{1},(j_{1}+1)\mathbf{e}_{1}}^{(k-1)}\quad{\rm{if}}\quad
j_{1}=-P+1,\ldots,0;\\\
I_{\mathbf{i},(j_{1}-1)\mathbf{e}_{1},j_{1}\mathbf{e}_{1}}^{(k-1)}\quad{\rm{for}}\quad
j_{1}=1,\ldots,P.\end{cases}$
$\displaystyle\mathcal{J}_{\mathbf{i},j_{2}}^{(k-1)}$ $\displaystyle=$
$\displaystyle\begin{cases}J_{\mathbf{i}-\mathbf{e}_{2},j_{2}\mathbf{e}_{2},(j_{2}+1)\mathbf{e}_{2}}^{(k-1)}\quad{\rm{for}}\quad
j_{2}=-P+1,\ldots,0;\\\
J_{\mathbf{i},(j_{2}-1)\mathbf{e}_{2},j_{2}\mathbf{e}_{2}}^{(k-1)}\quad{\rm{if}}\quad
j_{2}=1,\ldots,P.\end{cases}$

  * 5.

Compute

$\displaystyle\widetilde{S}^{P}_{1,\mathbf{i}}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,\frac{1}{2}}_{P}\Bigl{(}\mathcal{I}_{\mathbf{i},{*}}^{(k-1)},\Delta
x\Bigr{)};$ (6.8)
$\displaystyle\widetilde{S}^{P}_{2,\mathbf{i}}=\sum_{k=1}^{2P}\frac{\Delta
t^{k-1}}{k!}A^{0,\frac{1}{2}}_{P}\Bigl{(}\mathcal{J}_{\mathbf{i},{*}}^{(k-1)},\Delta
y\Bigr{)}.$ (6.9)

In the case of the well-balanced methods WBCAT$2P$ there is an important
difference: if the algorithm described in Subsection 4.1 (adopting the 2D
above notation) had to be used, the first step to update the numerical
solution at the point $\mathbf{x}_{\mathbf{i}}$ at a time $t_{n}$ would be to
find a solution to the problem

$\left\\{\begin{array}[]{l}F(U)_{x}+G(U)_{y}=S_{1}(U)H_{x}+S_{2}(U)H_{y}\\\
U(\mathbf{x}_{\mathbf{i}})=U_{\mathbf{i}}^{n}.\end{array}\right.$ (6.10)

This problem is obviously much more difficult to solve, either exactly or
numerically, than (4.3) since it is about a nonlinear PDE system instead of an
ODE system. Moreover in this case there may exist infinitely many stationary
solutions satisfying the condition at only one point
$\mathbf{x}_{\mathbf{i}}$: some extra conditions have to be imposed to
determine one of them.

Nevertheless, if the stationary solutions to be preserved constitute a
$k$-parameter family

$U^{*}(x,y;C_{1},\dots,C_{k}),$

with $k<d$, then the numerical strategy described in Subsection 4.4 can be
followed: this strategy will be used in Subsection 7.4 to preserve a family of
stationary solutions of the 2D Euler system with gravity.

The extension of the adaptive CAT2$P$ and WBCAT2$P$ to (6.1) is similar.
Following what has been done for ACAT2$P$ and WBACAT2$P$ in 1D we define
$\mathcal{A},$ the index set used to select the square stencils according with
the smoothness of numerical data, is defined as:

$\displaystyle\mathcal{A}_{\mathbf{i}}$ $\displaystyle=$
$\displaystyle\\{p\in\\{2,\dots,P\\}\text{ s.t.
}\psi^{p}_{\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{1}}\approx
1\quad\rm{and}\quad\psi^{p}_{\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{2}}\approx
1\\}$ (6.11)

where $\psi^{p}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{1}}$,
$\psi^{p}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{2}}$ are the smoothness
indicators introduced in Subsection 5.3 computed direction by direction. Then
the 2D adaptive CAT2P are so defined:

${F}_{\mathbf{i};\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{1}}^{\mathcal{A}_{\mathbf{i}}}=\begin{cases}{F}^{*}_{\mathbf{i};\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{1}}&\text{if
$\mathcal{A}_{\mathbf{i}}=\emptyset$;}\\\
{F}^{p_{s}}_{\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{1}}&\text{otherwise, where
$p_{s}=\max(\mathcal{A}_{\mathbf{i}});$ }\end{cases}$ (6.12)
${G}^{\mathcal{A}_{\mathbf{i}}}_{\mathbf{i};\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{2}}=\begin{cases}{G}^{*}_{\mathbf{i};\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{2}}&\text{if
$\mathcal{A}_{\mathbf{i}}=\emptyset$;}\\\
{G}^{p_{s}}_{\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{2}}&\text{otherwise, where
$p_{s}=\max(\mathcal{A}_{\mathbf{i}});$}\end{cases}$ (6.13)

and

$\widetilde{S}^{\mathcal{A}_{\mathbf{i}}}_{j\mathbf{i}}=\begin{cases}\widetilde{S}^{*}_{j,\mathbf{i}}&\text{if
$\mathcal{A}_{\mathbf{i}}=\emptyset$;}\\\
\widetilde{S}^{p_{s}}_{j,\mathbf{i}}&\text{otherwise, where
$p_{s}=\max(\mathcal{A}_{\mathbf{i}}).$}\end{cases}$ (6.14)

$F^{*}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{1}},$
$G^{*}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{2}},$ and
$\widetilde{S}_{j,\mathbf{i}}^{*}$ are the ACAT2 numerical flux and source
terms given by 1D (5.7), (5.8);
$F^{p_{s}}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{1}},$
$G^{p_{s}}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{2}}$ and
$\widetilde{S}_{j,\mathbf{i}}^{p_{s}}$ are the ACAT2$p_{s}$ numerical fluxes
and source terms defined in 1D (3.5), (3.11).

Using this notation, ACAT$2P$ methods can be extended as follows:

$U_{\mathbf{i}}^{n+1}=U_{\mathbf{i}}^{n}+\frac{\Delta t}{\Delta
x}\left[{F}_{\mathbf{i}-\frac{1}{2}\mathbf{e}_{1}}^{\mathcal{A}_{\mathbf{i}}}-{F}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{1}}^{\mathcal{A}_{\mathbf{i}}}+\widetilde{S}^{\mathcal{A}_{\mathbf{i}}}_{1,\mathbf{i}}\right]+\frac{\Delta
t}{\Delta
y}\left[{G}_{\mathbf{i}-\frac{1}{2}\mathbf{e}_{2}}^{\mathcal{A}_{\mathbf{i}}}-{G}_{\mathbf{i}+\frac{1}{2}\mathbf{e}_{2}}^{\mathcal{A}_{\mathbf{i}}}+\widetilde{S}^{\mathcal{A}_{\mathbf{i}}}_{2,\mathbf{i}}\right]$

In a similar way, the 2D adaptive well-balanced WBACAT2P writes as follows:

$U_{\mathbf{i}}^{n+1}=U_{\mathbf{i}}^{n}+\frac{\Delta t}{\Delta
x}\left[{F}_{\mathbf{i};\mathbf{i}-\frac{1}{2}\mathbf{e}_{1}}^{\mathcal{A}_{\mathbf{i}}}-{F}_{\mathbf{i};\mathbf{i}+\frac{1}{2}\mathbf{e}_{1}}^{\mathcal{A}_{\mathbf{i}}}+\widetilde{S}^{\mathcal{A}_{\mathbf{i}}}_{1,\mathbf{i}}\right]+\frac{\Delta
t}{\Delta
y}\left[{G}_{\mathbf{i};\mathbf{i}-\frac{1}{2}\mathbf{e}_{2}}^{\mathcal{A}_{\mathbf{i}}}-{G}_{\mathbf{i};\mathbf{i}+\frac{1}{2}\mathbf{e}_{2}}^{\mathcal{A}_{\mathbf{i}}}+\widetilde{S}^{\mathcal{A}_{\mathbf{i}}}_{2,\mathbf{i}}\right],$

where

${F}_{\mathbf{i};\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{1}}^{\mathcal{A}_{\mathbf{i}}}=\begin{cases}{F}^{*}_{\mathbf{i};\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{1}}&\text{if
$\mathcal{A}_{\mathbf{i}}=\emptyset$;}\\\
{F}^{p_{s}}_{\mathbf{i};\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{1}}&\text{otherwise,
where $p_{s}=\max(\mathcal{A}_{\mathbf{i}});$}\end{cases}$ (6.15)
${G}^{A}_{\mathbf{i};\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{2}}=\begin{cases}{G}^{*}_{\mathbf{i};\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{2}}&\text{if
$\mathcal{A}_{\mathbf{i}}=\emptyset$;}\\\
{G}^{p_{s}}_{\mathbf{i};\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{2}}&\text{otherwise,
where $p_{s}=\max(\mathcal{A}_{\mathbf{i}});$}\end{cases}$ (6.16)

and

$\widetilde{S}^{\mathcal{A}_{\mathbf{i}}}_{j\mathbf{i}}=\begin{cases}\widetilde{S}^{*}_{j,\mathbf{i}}&\text{if
$\mathcal{A}_{\mathbf{i}}=\emptyset$;}\\\
\widetilde{S}^{p_{s}}_{j,\mathbf{i}}&\text{otherwise, where
$p_{s}=\max(\mathcal{A}_{\mathbf{i}}).$}\end{cases}$ (6.17)

Here, $F^{*}_{\mathbf{i};\mathbf{i}+\frac{1}{2}\mathbf{e}_{1}},$
$G^{*}_{\mathbf{i};\mathbf{i}+\frac{1}{2}\mathbf{e}_{2}},$ and
$\widetilde{S}_{j,\mathbf{i}}^{*}$ are the ACAT2 numerical flux and source
terms given by 1D (5.14),(5.15);
$F^{p_{s}}_{\mathbf{i};\mathbf{i}+\frac{1}{2}\mathbf{e}_{1}},$
$G^{p_{s}}_{\mathbf{i};\mathbf{i}+\frac{1}{2}\mathbf{e}_{2}}$ and
$\widetilde{S}_{j,\mathbf{i}}^{p_{s}}$ are the ACAT2$p_{s}$ numerical fluxes
and source terms defined in 1D (4.12), (4.11), (4.9).

##  7 Numerical results

In this section we apply ACAT$2P$ and WBACAT$2P$, $P=1,2$ methods to several
1D and 2D problems: the 1D linear transport equation and Burgers equation with
source term, the 1D shallow water system, and 2D Euler equations with a
gravitational potential. The Minmod flux limiter [25] is used in ACAT2 and the
smoothness indicators (5.18) are used for ACAT4: no loss of precision for
first order critical points has been observed in any of the test problems
considered here due to the use of $\psi^{2}_{i+\frac{1}{2}}$. Fornberg’s
algorithm [13] is used to compute the coefficients of the numerical
differentiation formulas.

###  7.1 Linear Equation

We consider the linear scalar balance law

$u_{t}+u_{x}=u,$ (7.1)

that has the form (1.1) with $H(x)=x$. The analytic solution of the initial
value problem with condition

$u(x,0)=u_{0}(x)$

is given by:

$u(x,t)=u_{0}(x-t)e^{t}.$ (7.2)

The stationary solutions solve the ODE

$u_{x}=u.$

Hence, the set of stationary solutions is

$u^{*}(x)=Ce^{x},\quad C\in\mathbb{R}.$

####  7.1.1 Order test

Following [24], we consider (7.1) with initial condition

$u_{0}(x)=\begin{cases}0\quad\quad\;\;\;\mathrm{if}\quad x<0;\\\
p(x)\quad\;\;\mathrm{if}\quad 0\leq x\leq 1;\\\
1\quad\quad\;\;\;\mathrm{if}\quad x>1;\end{cases}$ (7.3)

Figure 7.1: Test 7.1.1. Initial condition and exact solution (left);
difference between the numerical and the exact solutions computed with CAT2,
CAT4, WBCAT2 and WBCAT4 at $t=1$ using a mesh of 41 points (right) on the
interval $[-0.2,2]$ and CFL$=0.9.$

where $p(x)$ is the polynomial that satisfies $p(0)=0$, $p(1)=1$,
$p^{k}(0)=p^{k}(1)=0$, $k=1,\ldots,5$:

$p(x)=x^{6}\Bigl{(}\sum_{k=0}^{5}(-1)^{k}\binom{5+k}{k}(x-1)^{k}\Bigr{)}$

(see Figure 7.1). The methods ACAT2, ACAT4, WBACAT2, WBACAT4 have been applied
to (7.1) with initial condition (7.3) in the spatial interval $[-0.2,2],$ with
CFL$=0.9.$ Dirichlet boundary conditions are considered to the left and free
boundary conditions to the right based on the use of ghost cells. Figure 7.1
shows the numerical solutions obtained at time $t=1$ on the interval
$[-0.2,2]$ using $40$ mesh points.

Tables 1-2 show the $L^{1}$-errors and the empirical order of convergence
corresponding to the standard and Adaptive CAT2P and WBCAT2P with $p=1,2.$ As
it can be seen, all the schemes keep the expected order and the errors
corresponding to methods of the same order are almost identical. In the first
case, the smoothness indicators of ACAT2P and WBACAT2P have been fixed to 1,
hence, the Adaptive CAT2P coincides exactly with the standard CAT2P scheme;
regarding the second case, no restrictions are imposed on the smoothness
indicators, requiring an increase in the number of points to capture the
theoretical order. No further restrictions are required for the time step.

| CAT2 | WBCAT2 | CAT4 | WBCAT4
---|---|---|---|---
Points | Order | Error | Order | Error | Order | Error | Order | Error
6 | - | 1.79E-1 | - | 1.80E-1 | - | 6.03E-2 | - | 6.03E-2
11 | 1.57 | 6.00E-2 | 1.56 | 6.05E-2 | 2.86 | 8.31E-3 | 2.86 | 8.31E-3
21 | 1.91 | 1.60E-2 | 1.92 | 1.60E-2 | 3.58 | 6.94E-4 | 3.58 | 6.94E-4
41 | 2.02 | 3.93E-3 | 2.02 | 3.94E-3 | 3.88 | 4.69E-5 | 3.89 | 4.69E-5
81 | 2.03 | 9.63E-4 | 2.02 | 9.66E-4 | 4.01 | 2.90E-6 | 4.00 | 2.90E-7
Table 1: Test 7.1.1: Errors in $L^{1}$ norm and convergence rates for CAT2, CAT4, WBCAT2 and WBCAT4 at time $t=1.$ | ACAT2 | WBACAT2 | ACAT4 | WBACAT4
---|---|---|---|---
Points | Order | Error | Order | Error | Order | Error | Order | Error
16 | - | 2.00E-1 | - | 1.99E-1 | - | 1.57E-2 | - | 1.61E-2
31 | 1.51 | 7.01E-2 | 1.50 | 7.04E-2 | 3.01 | 1.95E-3 | 2.99 | 2.02E-3
61 | 1.85 | 1,94E-2 | 1.86 | 1.94E-2 | 3.62 | 1.58E-4 | 3.60 | 1.67E-4
121 | 2.19 | 4.25E-3 | 2.19 | 4.25E-3 | 8.10 | 5.71E-7 | 8.20 | 5.71E-7
241 | 2.00 | 1.05E-3 | 2.00 | 1.06E-3 | 3.99 | 3.57E-8 | 4.00 | 3.56E-8
481 | 2.00 | 2.64E-4 | 2.00 | 2.64E-4 | 4.00 | 2.22E-9 | 4.00 | 2.22E-9
Table 2: Test 7.1.1: Errors in $L^{1}$ norm and convergence rates for ACAT2,
ACAT4, WBACAT2 and WBACAT4 at time $t=1.$

###  7.2 Burgers Equation

In this section we consider the scalar Burgers equation with source term:

$u_{t}+\Bigl{(}\frac{u^{2}}{2}\Bigr{)}_{x}=u^{2}H_{x}(x).$ (7.4)

The stationary solutions solve now the ODE

$u_{x}=uH_{x}(x)$

whose general solution is

$u^{*}(x)=Ce^{H(x)},\quad C\in\mathbb{R}.$

####  7.2.1 Preservation of a stationary solution

As first test let us check numerically the well-balanced property. For this
reason, we consider an oscillatory

$H(x)=x+0.1\sin(10x)$

and consequently the stationary solutions become $u^{*}(x)=Ce^{H(x)}.$ We
solve (7.4) with initial condition

$u_{0}(x)=e^{H(x)}$

in the interval $[-1,1]$ using $100$ mesh points and CFL$=0.9.$ As boundary
conditions, the stationary solution is imposed at ghost cells.

Table 3 shows the $L^{1}$-errors and the empirical order of convergence
corresponding to ACAT$2P$, and the error of WBACAT$2P$ methods, with $P=1,2$
at time $t=8.$ While the non well-balanced schemes give accurate solutions
according to their order, the well-balanced methods capture the stationary
solution with machine precision.

| ACAT2 | WBACAT2 | ACAT4 | WBACAT4
---|---|---|---|---
Points | Order | Error | Error | Order | Error | Error
100 | - | 1.82E-3 | 2.66E-17 | - | 1.93E-5 | 3.99E-17
200 | 2.05 | 4.32E-4 | 3.11E-17 | 4.09 | 1.13E-6 | 6.21E-17
400 | 2.03 | 1.07E-4 | 2.44E-17 | 4.05 | 6.84E-8 | 2.22E-17
800 | 2.01 | 2.63E-5 | 2.78E-17 | 4.03 | 4.19E-9 | 3.55E-17
1600 | 2.00 | 6.55E-6 | 1,99E-17 | 4.01 | 2.59E-10 | 5.21E-17
Table 3: Test 7.2.1. Test of well-balance property. Errors in $L^{1}$ norm and
convergence rates at time $t=8$ for ACAT2-4. The errors for WBACAT2-4 are due
to round-off.

Checked numerically the well-balanced property for the Burgers equation with
source we will now focus on some experiment to test the accuracy of the
numerical solutions where initial conditions are: a perturbation of the
stationary solution; and a smooth perturbation of Test 7.1.1.

####  7.2.2 Perturbation of a stationary solution

Let us consider (7.4) with oscillatory $H$ given by

$H(x)=x+0.1\sin(10x)$

and initial condition

$u_{0}(x)=e^{H(x)}+0.02e^{-200(x+0.7)^{2}},$ (7.5)

that is a small smooth perturbation of the stationary solution
$u^{*}(x)=e^{H(x)}$: see Figure 7.2.

Figure 7.2: Test 7.2.2. Initial condition and stationary solutions (left).
Difference between reference and stationary solution at initial and final time
(right). The perturbation of the initial condition (left) is amplified by 10
times in order to see clearly the perturbation. The reference solution is
computed with WBACAT4 using 1000 mesh points and CFL$=0.9$ at time $t=1.2.$

We solve the problem in the interval $[-1,1]$ using 200 mesh points and
CFL$=0.9.$ As boundary conditions the stationary solution is imposed at left
ghost point and free boundary at right.

Figure 7.3: Test 7.2.2. Difference between numerical solutions computed with
ACAT$2P$ (top) and WBCAT$2P,$ (bottom) $P=1,2,$ and reference solution at
$t=1.2$ using a mesh of 200 points and CFL$=0.9.$ For the reference solution
the WBACAT4 is adopted with a mesh of 1000 points.

Figure 7.3 shows that, all the schemes are able to evolve the perturbation in
according with the order. Nevertheless, the well-balanced methods, WBACAT2 and
WBACAT4, are able to capture more precisely the evolution of the perturbation
with a smaller error than the relative non well-balanced schemes.
For the errors in $L^{1}-$norm and convergence rates, we adopt as initial
condition

$u_{0}(x)=e^{H(x)}+0.2e^{-200(x+0.7)^{2}},$

in other word, a bigger perturbation is considered reducing the final time to
$t=0.2$.
Table 4 shows that the non well-balanced methods introduce a bigger error in
comparison with the well-balanced approach but an increasing of points is
necessary to achieve the theoretical order. This phenomenon is partly
attributable to reconstruction partly to smoothness indicators, because they
fail to detect the theoretical regularity.

| ACAT2 | WBACAT2 | ACAT4 | WBACAT4
---|---|---|---|---
Points | Order | Error | Order | Error | Order | Error | Order | Error
80 | - | 6.44E-2 | - | 1.07E-2 | - | 1.53E-3 | - | 1.17E-4
160 | 1.15 | 2.88E-2 | 1.74 | 3.21E-3 | 2.11 | 3.10E-4 | 2.03 | 2.86E-5
320 | 1.45 | 1.05E-2 | 1.89 | 8.61E-4 | 4.28 | 1.59E-5 | 4.45 | 1.31E-6
640 | 1.62 | 3.42E-3 | 1.96 | 2.21E-4 | 4.14 | 8.99E-7 | 3.85 | 8.99E-8
1280 | 1.86 | 9.41E-4 | 1.99 | 5.58E-5 | 4.01 | 5.57E-8 | 3.96 | 5.76E-9
Table 4: Test 7.2.2: Errors in $L^{1}$ norm and convergence rates for ACAT2,
ACAT4, WBACAT2 and WBACAT4 at time $t=0.2$ and CFL$=0.9.$

This experiment is very interesting because, on the one hand it shows that the
well-balanced methods work well near the stationary solution, on the other
hand, it happens that a small perturbation of the initial state, although
rather smooth, may result in a loss of accuracy in the adaptive order
reconstruction.

####  7.2.3 Preservation of a stationary solution with oscillatory $H$

Following [24], we consider (7.4) with

$H(x)=x+\frac{1}{10}\sin(100x),$ (7.6)

and we take as initial condition the stationary solution

$u^{*}(x)=e^{H(x)}$

(see Figure 7.4).

Figure 7.4: Test 7.2.3. Initial condition (top) and $H$ (down).

We solve the problem in the interval $[-1,1]$ using 100 mesh points and
$CFL=0.9.$ With this choice of mesh points, the period of the oscillations of
$H$ is close to $\Delta x.$ As boundary conditions the stationary solution is
imposed again at ghost points.

Figure 7.5: Test 7.2.3. Exact and numerical stationary solutions computed with
ACAT$2P$ and WBCAT$2P$ at $t=1$ using a mesh of 100 points and CFL$=0.9$: non
well-ballanced (top) and well-balanced (bottom). Figure 7.6: Test 7.2.3. Top:
differences between the exact and the numerical stationary solutions computed
with ACAT$2P$, $P=1,2$, at time $t=1$ using 100 mesh points and $CFL=0.9$.
Bottom: differences between the exact and the numerical solutions computed
with WBACAT$2P$, $P=1,2$, at time $t=1$ using 100 mesh points and $CFL=0.9$.

Figure 7.5 shows that, while WBACAT2 and WBACAT4 capture the stationary
solution with machine precision, this is not the case for ACAT2 and ACAT4.
Figure 7.6 displays the differences between the numerical solution and the
stationary solution obtained at time $t=1$ (top). In this case, the results
provided by WBACAT2 and WBACAT4 (bottom) are very similar and they are able to
capture the machine precision; while, the results provided by ACAT2 and ACAT4
(top) show that the non well-balanced methods are not able to detect the
stationary solution with high precision even if a 4 order method is applied.

As a final check we consider the behaviour of the methods in the case of an
initial condition of class $\mathcal{C}^{5}$ which is far from the stationary
solution.

####  7.2.4 Order Test

Let us consider (7.4) with $H(x)=x$ and initial condition (7.3) (see Figure
7.7)

$u_{0}(x)=\begin{cases}0\quad\quad\;\;\;\mathrm{if}\quad x<0;\\\
p(x)\quad\;\;\mathrm{if}\quad 0\leq x\leq 1;\\\
1\quad\quad\;\;\;\mathrm{if}\quad x>1;\end{cases}$ (7.7)

where

$p(x)=x^{6}\Bigl{(}\sum_{k=0}^{5}(-1)^{k}\binom{5+k}{k}(x-1)^{k}\Bigr{)}.$

We solve the problem in the interval $[-0.2,2]$ using 80 mesh point and
CFL$=0.9.$ As boundary condition free boundary is imposed at ghost points.

Figure 7.7: Test 7.2.4. Initial condition and Reference solution obtained with WBACAT4 using a 2560 mesh points and CFL$=0.9$ at time $t=0.5.$ Figure 7.8: Test 7.2.4. Differences between numerical solutions computed with ACAT$2P$ and WBCAT$2P$, $P=1,2$, and the reference solution at $t=0.5$ using a mesh of 80 points and CFL$=0.9.$ For the reference solution a mesh of 2560 has been adopted. | ACAT2 | WBACAT2 | ACAT4 | WBACAT4
---|---|---|---|---
Points | Order | Error | Order | Error | Order | Error | Order | Error
80 | - | 1.74E-3 | - | 1.91E-3 | - | 2.41E-4 | - | 2.42E-4
160 | 1.93 | 4.11E-4 | 1.92 | 5.05E-4 | 2.40 | 4.56E-5 | 2.41 | 4.57E-5
320 | 1.93 | 1.08E-4 | 1.94 | 1.32E-4 | 3.17 | 5.05E-6 | 3.18 | 5.06E-6
640 | 1.98 | 2.74E-5 | 1.98 | 3.34E-5 | 3.68 | 3.95E-7 | 3.68 | 3.94E-7
1280 | 1.99 | 6.91E-6 | 1.99 | 8.41E-6 | 3.91 | 2.60E-8 | 3.92 | 2.61E-8
2560 | 2.00 | 1.73E-6 | 2.00 | 2-11E-6 | 3.98 | 1.65E-9 | 3.98 | 1.65E-9
Table 5: Test 7.2.4: Errors in $L^{1}$ norm and convergence rates for ACAT2,
ACAT4, WBACAT2 and WBACAT4 at time $t=0.5$ and CFL$=0.9.$

Figure 7.8 and Table 5 show us how all the methods manage to produce solutions
in agreement with each other obtaining the expected order.

This experiment shows that, when the initial condition is far from the
stationary solution, well-balanced and non well-balanced methods produce
essentially the same results, with the expected order of accuracy.

###  7.3 Shallow water model

In this section we will focus on the one-dimensional hyperbolic shallow water
model

$\begin{cases}h_{t}+q_{x}=0\\\ \displaystyle
q_{t}+\left(\frac{q^{2}}{h}+\frac{g}{2}h^{2}\right)_{x}=ghH_{x},\end{cases}$
(7.8)

that can be written in the form (1.1) with

$U=\begin{bmatrix}h\\\ q\end{bmatrix},\quad F(U)=\begin{bmatrix}q\\\
\displaystyle\frac{q^{2}}{h}+\frac{g}{2}h^{2}\end{bmatrix},\quad
S(U)=\begin{bmatrix}0\\\ gh\end{bmatrix}.$

The variable $x$ refers to the axis of the channel and $t$ is time; $q(x,t)$
and $h(x,t)$ represent the discharge and the water thickness; $g,$ the
acceleration due to gravity; $H(x),$ the depth measured from a fixed level of
reference; furthermore, the following relation is verified
$q(x,t)=h(x,t)u(x,t),$ with $u$ the depth average horizontal velocity. The
eigenvalues of the Jacobian matrix $J(U)$ of the flux $F(U)$ are

$\lambda_{1}=u-\sqrt{gh}\quad{\rm{and}}\quad\lambda_{2}=u+\sqrt{gh}.$

The local Froude number is defined by

$Fr=\frac{|u|}{\sqrt{gh}}.$

The flow is said to be supercritical if $Fr>1$ for all $x\in[a,b],$ critical
if $Fr=1$ for all $x\in[a,b]$ and subcritical if $Fr<1$ for all $x\in[a,b].$

The stationary solution of the shallow water system (7.8) are implicitly given
by

$q=Q\quad{\rm{and}}\quad\frac{1}{2}\frac{q^{2}}{h^{2}}+gh-gH=C,$ (7.9)

where $Q$ and $C$ are arbitrary constants [7]. In order to implement the well-
balanced methods, given $U^{n}_{i}=[h^{n}_{i},q^{n}_{i}]^{T}$ one has to find
the stationary solution $U^{*}_{i}=[q^{*}_{i},h^{*}_{i}(x)]^{T}$ that solves
(4.3): it is implicitly given by

$q^{*}_{i}(x)=q^{n}_{i},\quad\frac{1}{2}\frac{(q^{n}_{i})^{2}}{{h^{*}_{i}(x)}^{2}}+gh^{*}_{i}(x)-gH=C_{i},\quad\forall
x,$

with

$C_{i}=\frac{1}{2}\frac{(q^{n}_{i})^{2}}{(h_{i}^{n})^{2}}+gh^{n}_{i}-gH(x_{i}).$

Therefore, at a point $x_{j}$ of the stencil, one has
$q_{i}^{*}(x_{j})=q^{n}_{i}$ and $h_{i}^{*}(x_{j})$ has to be a positive root
of the polynomial:

$P_{i,j}(h)=h^{3}-\left(\frac{C_{i}}{g}+H(x_{j})\right)h^{2}+\frac{1}{2g}(q_{i}^{n})^{2}.$

This polynomial can have two, one, or zero positive roots. In the first case,
one of the roots corresponds to a supercritical state and the other one to a
subcritical state: a criterion is necessary to select one root or the other.
We follow here a criterion similar to the one chosen in [7] in the context of
finite volume methods: the solution whose regime (sub or supercritical) is the
same as the one of $U^{n}_{i}$ is selected. A careful implementation is needed
to capture transcritical stationary solutions: see for instance the discussion
in [7] or [24].

####  7.3.1 Preservation of a subcritical stationary solution

Let us consider the shallow water model in the space interval $[-3,3]$ with
bottom depth given by

$H(x)=\begin{cases}-0.25(1+\cos(5\pi x))\quad\rm{if}\;-0.2\leq x\leq 0.2;\\\
0\quad\rm{otherwise};\end{cases}$ (7.10)

and initial condition given by the subcritical stationary solution $U^{*}$
that satisfies

$q^{*}=2.5,\quad h^{*}(-3)=2$

(see Figure 7.9). The numerical methods are applied to this problem using 200
mesh points and CFL$=0.8$. At the boundaries, the stationary solution is
imposed at ghost points.

Figure 7.9: Test 7.3.1. Discrete initial condition with 100 mesh points. Free
surface and bathymetry. Figure 7.10: Test 7.3.1. Difference between the
numerical solution for $h$ (top) and $q$ (bottom) obtained with ACAT methods
and the exact stationary one, at time $t=4$ using 200 mesh points and
CFL$=0.8.$

Figures 7.10 shows the differences between the exact and the non well-balanced
numerical solutions obtained at time $t=4.$ As it can be seen, the well-
balanced methods capture the stationary solution to machine accuracy. This
behaviour is confirmed by Table 6 that shows the $L^{1}-$errors corresponding
to WBACAT$2P$, $P=1,2$, using $50,$ $100,$ $200$ and $400$ mesh points at time
$t=4.$

| WBACAT2 | WBACAT4
---|---|---
Points | h | q | u | h | q | u
50 | 2.93E-16 | 1.07E-16 | 2.66E-16 | 2.39E-16 | 5.32E-17 | 1.87E-16
100 | 3.46E-16 | 7.99E-17 | 1.86E-16 | 2.13E-16 | 0 | 1.20E-16
200 | 3.40E-16 | 0 | 2.46E-16 | 3.99E-17 | 0 | 1.99E-17
400 | 1.77E-16 | 0 | 1.20E-16 | 0 | 5.99E-17 | 2.98E-17
Table 6: Test 7.3.1. Errors in $L^{1}$ norm for WBACAT$2P$, $P=1,2$, at time
$t=4.$

The introduction of spurious oscillations with the not well-balanced schemes
involves a order reduction since the high order smoothness indicators are not
able to detect a priori the real smoothness of the solution. This behaviour is
highlighted in the next experiments.

####  7.3.2 Perturbation of a subcritical stationary solution

The setting of this test is similar to the previous one but now the initial
condition is a smooth perturbation of the subcritical stationary solution
$U^{*}$ (see Figure 7.11) considered there:

$U_{0}=\begin{bmatrix}h^{*}+0.006e^{(-20(x+1)^{2})}\\\ q^{*}\end{bmatrix}.$

Figure 7.11: Test 7.3.2. Initial condition and reference solution obtained
with WBACAT4 computed at time $t=0.4$ using $2000$ mesh points and CFL$=0.8:$
$h$ (left); $q$ (right). In the plot of $q$ there appear the left and right
traveling waves, as well as a small left moving reflected wave. Figure 7.12:
Test 7.3.2. Difference between Reference and numerical solutions obtained with
WBACAT$2P$, $P=1,2$, computed at time $t=0.4$ using $200$ mesh points and
CFL$=0.8:$ $h$ (top); $q$ (bottom). The Reference solution is computed with
WBACAT4 adopting a $2000$ mesh points. Figure 7.13: Test 7.3.2. Difference
between Reference and numerical solutions obtained with ACAT$2P$, $P=1,2$,
computed at time $t=0.4$ using $200$ mesh points and CFL$=0.8:$ $h$ (top); $q$
(bottom). The Reference solution is computed with WBACAT4 adopting a $2000$
mesh points

The numerical solutions are computed on the interval $[-3,3]$ using $200$ mesh
points at time $t=0.4$ with CFL $=0.8.$ As boundary conditions the subcritical
stationary solution is imposed at ghost points.

Figures 7.12 and 7.13 show the errors obtained by the differences between
reference solution and the numerical solutions computed with well-balanced and
not well-balanced methods at time $t=0.4.$ The Reference solution considered
is WBACAT4 adopting a $2000$ mesh points. As expected, WBACAT$2P$, $P=1,2$,
capture better the waves generated by the initial perturbation than ACAT$2P$,
$P=1,2$.

| ACAT2 | WBACAT2 | ACAT4 | WBACAT4
---|---|---|---|---
Points | Order | Error | Order | Error | Order | Error | Order | Error
50 | - | 2.18E-3 | - | 1.64E-5 | - | 2.07E-4 | - | 1.81E-5
100 | 0.98 | 1.10E-3 | 1.95 | 4.22E-6 | 1.11 | 9.61E-5 | 2.25 | 3.79E-6
200 | 1.23 | 4.96E-4 | 1.94 | 1.09E-6 | 1.59 | 3.19E-5 | 4.72 | 1.44E-7
400 | 1.48 | 1.68E-4 | 1.97 | 2.77E-7 | 1.84 | 8.87E-6 | 4.33 | 8.22E-9
800 | 1.53 | 5.82E-5 | 1.97 | 7.07E-8 | 1.93 | 2.31E-6 | 4.07 | 6.54E-10
Table 7: Test 7.3.2: Errors in $L^{1}$ norm and convergence rates related to $h$ for ACAT2, ACAT4, WBACAT2 and WBACAT4 at time $t=0.15$ and CFL$=0.8.$ | ACAT2 | WBACAT2 | ACAT4 | WBACAT4
---|---|---|---|---
Points | Order | Error | Order | Error | Order | Error | Order | Error
50 | - | 8.21E-3 | - | 5.54E-5 | - | 7.85E-4 | - | 6.15E-5
100 | 1.00 | 4.11E-3 | 1.63 | 1.79E-5 | 1.22 | 3.35E-4 | 2.86 | 1.69E-5
200 | 1.25 | 1.73E-3 | 1.94 | 4.68E-6 | 1.67 | 1.05E-4 | 4.78 | 6.17E-7
400 | 1.51 | 6.06E-4 | 1.98 | 1.19E-6 | 1.87 | 2.86E-5 | 4.31 | 4.65E-8
800 | 1.54 | 2.08E-4 | 1.98 | 3.01E-7 | 1.98 | 7.31E-6 | 4.05 | 1-41E-9
Table 8: Test 7.3.2: Errors in $L^{1}$ norm and convergence rates related to
$q$ for ACAT2, ACAT4, WBACAT2 and WBACAT4 at time $t=0.15$ and CFL$=0.8.$

In addition, Tables 7-8 show how well-balanced methods manage to reach the
expected order, behavior not respected by non-well-balanced methods. In this
case, as seen above, this phenomenon is partly attributable to not well-
balanced reconstruction partly to smoothness indicators. In fact, the not
well-balanced method at first step introduces a spurious error which implies a
loss of numerical smoothness resulting in degradation of the order.

####  7.3.3 Smooth initial condition with flat bottom

We now check that in the case of flat bottom and smooth solution, well-
balanced and non well-balanced schemes give the same result, all with the
expected order of accuracy. In order to obtain these results we consider the
Shallow water equation (7.8) with flat bottom and smooth initial condition
(7.3) (see Figure 7.14):

$U_{0}(x)=\begin{bmatrix}h_{0}(x)\\\ q_{0}(x)\end{bmatrix},$ (7.11)

where

$h_{0}(x)=q_{0}(x)=\begin{cases}0\quad\quad\;\;\;\mathrm{if}\quad x<0;\\\
p(x)\quad\;\;\mathrm{if}\quad 0\leq x\leq 1;\\\
1\quad\quad\;\;\;\mathrm{if}\quad x>1;\end{cases}$

and

$p(x)=x^{6}\Bigl{(}\sum_{k=0}^{5}(-1)^{k}\binom{5+k}{k}(x-1)^{k}\Bigr{)}.$

Figure 7.14: Test 7.3.3. Initial condition and reference solution obtained by WBACAT4 at time $t=0.2$ using $3200$ mesh points and CFL$=0.8;$ $h$ (top); $q$ (bottom). Figure 7.15: Test 7.3.3. Differences between numerical solutions obtained with ACAT$2P$ and WBACAT$2P$, $P=1,2$, computed at time $t=0.2$ using $100$ mesh points and CFL$=0.8$ and the reference solution. $h$ (top); $q$ (bottom). For the reference solution a 3200 mesh points has been adopted.  | ACAT2 | WBACAT2 | ACAT4 | WBACAT4
---|---|---|---|---
Points | Order | Error | Order | Error | Order | Error | Order | Error
200 | - | 7.27E-5 | - | 7.25E-5 | - | 2.25E-6 | - | 2.25E-6
400 | 2.01 | 1.80E-5 | 2.01 | 1.80E-5 | 3.98 | 1.42E-7 | 3.98 | 1.42E-7
800 | 2.00 | 4.50E-6 | 2.00 | 4.50E-6 | 3.99 | 8.91E-9 | 3.99 | 8,91E-9
1600 | 2.00 | 1.13E-6 | 2.00 | 1.13E-6 | 4.00 | 5.56E-10 | 4.00 | 5.57E-10
Table 9: Test 7.3.3: Errors in $L^{1}$ norm and convergence rates related to $h$ for ACAT2, ACAT4, WBACAT2 and WBACAT4 at time $t=0.2$ and CFL$=0.9.$ | ACAT2 | WBACAT2 | ACAT4 | WBACAT4
---|---|---|---|---
Points | Order | Error | Order | Error | Order | Error | Order | Error
200 | - | 8.75E-4 | - | 8.74E-4 | - | 2.67E-5 | - | 2.67E-5
400 | 2.01 | 2.17E-4 | 2.01 | 2.17E-4 | 3.98 | 1.69E-6 | 3.98 | 1.69E-6
800 | 2.00 | 5.42E-5 | 2.00 | 5.42E-5 | 3.99 | 1.06E-7 | 3.99 | 1.06E-7
1600 | 2.00 | 1.35E-5 | 2.00 | 1.35E-5 | 4.00 | 6.61E-9 | 4.00 | 6.61E-9
Table 10: Test 7.3.3: Errors in $L^{1}$ norm and convergence rates related to
$q$ for ACAT2, ACAT4, WBACAT2 and WBACAT4 at time $t=0.2$ and CFL$=0.9.$

The numerical solutions are computed on the interval $[-2,4]$ using $100$ mesh
points and CFL $=0.9$ at time $t=0.1,$ while for the reference solution a
$3200$ mesh points is adopted. As boundary conditions free boundary is imposed
at ghost points.

Figure 7.15 and Table 9-10 show that all methods have a similar behavior by
reproducing similar results. In particular, all schemes are accurate as
expected.

With this experiment we have proven that all methods, both well-balanced and
not, have a similar behavior when they are very far from the stationary
condition; while, well-balanced reconstructions reproduce better results, both
in accuracy and numerical convergence, when we are close to the stationary
solution.

###  7.4 2D Euler system with gravity

As last test, let us consider the 2D system of compressible Euler equations
with a gravity

$\begin{cases}\rho_{t}+(\rho u)_{x}+(\rho v)_{y}=0,\\\ (\rho u)_{t}+(\rho
u^{2}+p)_{x}+(\rho uv)_{y}=-\rho H_{x},\\\ (\rho v)_{t}+(\rho uv)_{x}+(\rho
v^{2}+p)_{y}=-\rho H_{y},\\\ E_{t}+(u(E+p))_{x}+(v(E+p))_{y}=-\rho uH_{x}-\rho
vH_{y}.\end{cases}$ (7.12)

Here, $\rho$ is the density; $u$, the velocity in $x-$direction; $v$, the
velocity in $y-$direction; $p$, the pressure; $E$, the energy per unit volume
excluding the gravitational energy; and $H(x,y)$, the gravitational potential
[20]. The pressure is supposed to satisfy the equation of state

$p=(\gamma-1)\Bigl{(}E-\frac{1}{2}\rho(u^{2}+v^{2})\Bigr{)},$

where $\gamma$ is the ratio between specific heats at constant pressure and
volume, which is taken to be constant. System (7.12) can be written in the
form (6.1) with

$U=\begin{bmatrix}\rho\\\ \rho u\\\ \rho v\\\ E\end{bmatrix},\quad
F(U)=\begin{bmatrix}\rho u\\\ \rho u^{2}+p\\\ \rho uv\\\
u(E+p)\end{bmatrix},\quad G(U)=\begin{bmatrix}\rho v\\\ \rho uv\\\ \rho
v^{2}+p\\\ v(E+p)\end{bmatrix},\quad S(U)=\begin{bmatrix}0\\\ -\rho H_{x}\\\
-\rho H_{y}\\\ -\rho uH_{x}-\rho vH_{y}\end{bmatrix}.$

Hydrostatic stationary solutions satisfy

$u=0,\quad v=0,\quad\nabla p=-\rho\nabla H.$

A two-parameter family of isothermal stationary solution is given by

$\rho^{*}(\mathbf{x})=C_{1}e^{-H(\mathbf{x})}\geq 0;\quad
p^{*}(\mathbf{x})=C_{2}\rho^{*}(\mathbf{x})\geq 0;\quad u^{*}=v^{*}=0;\quad
E^{*}(x)=\frac{p^{*}(x)}{\gamma-1}.$ (7.13)

Given

$U_{\mathbf{i}}^{n}=[\rho_{\mathbf{i}}^{n},\rho_{\mathbf{i}}^{n}u_{\mathbf{i}}^{n},\rho_{\mathbf{i}}^{n}v_{\mathbf{i}}^{n},E_{\mathbf{i}}^{n}]^{T},$

the stationary solution $U^{*}_{\mathbf{i}}$ selected by applying the
technique described in Subsection 4.4 is then

$\rho_{\mathbf{i}}^{*}(\mathbf{x})=\rho_{\mathbf{i}}^{n}e^{-(H(\mathbf{x})-H(\mathbf{x}_{\mathbf{i}}))};\;p^{*}_{\mathbf{i}}(\mathbf{x})=\rho_{\mathbf{i}}^{n}e^{-(H(\mathbf{x})-H(\mathbf{x}_{\mathbf{i}}))};\;u_{\mathbf{i}}^{*}=v_{\mathbf{i}}^{*}=0;\;E_{\mathbf{i}}^{*}(x)=\frac{p^{*}_{\mathbf{i}}(x)}{\gamma-1}.$
(7.14)

###  7.5 Preservation of a continuous stationary solution

Following [18]-[20]-[16]-[17] we consider Euler equations in the 2D domain
$[0,1]\times[0,1]$ with two different gravitational potentials

$H_{1}(x,y)=x+y,\quad
H_{2}(x,y)=\dfrac{1}{\sqrt{(x-\frac{1}{3})^{2}+(y+\frac{1}{2})^{2}}}$

and initial condition

$\displaystyle\rho(\mathbf{x},0)=e^{-H(\mathbf{x})};\quad
p(\mathbf{x},0)=e^{-H(\mathbf{x})};\quad u(\mathbf{x},0)=v(\mathbf{x},0)=0.$
(7.15)

| 2D density
---|---
| 2D ACAT2 | 2D ACAT4 | 2D WBACAT2 | 2D WBACAT4
Points | Error | Order | Error | Order | Error | Error
20$\times$20 | 4.87E-6 | - | 7.85E-9 | - | 2.72E-17 | 2.96E-18
40$\times$40 | 1.91E-6 | 1.35 | 1.01E-9 | 2.95 | 2.19E-17 | 2.39E-18
80$\times$80 | 5.62E-7 | 1.76 | 8.54E-11 | 3.56 | 1.82E-17 | 2.39E-18
160$\times$160 | 1.43E-7 | 1.98 | 5.61E-12 | 3.93 | 2.37E-18 | 2.64E-18
Table 11: Test 7.5: 2D Euler equation with gravity and gravitational potential $H_{1}$. Errors in $L^{1}$ norm for density at time $t=0.3.$ | 2D density
---|---
| 2D ACAT2 | 2D ACAT4 | 2D WBACAT2 | 2D WBACAT4
Points | Error | Order | Error | Order | Error | Error
20$\times$20 | 3.85E-5 | - | 3.87E-5 | - | 2.50E-17 | 3.77E-17
40$\times$40 | 1.58E-5 | 1.28 | 5.16E-6 | 2.91 | 3.23E-17 | 3.59E-17
80$\times$80 | 4.78E-6 | 1.72 | 4.45E-7 | 3.53 | 3.33E-17 | 3.33E-17
160$\times$160 | 1.23E-6 | 1.96 | 2.89E-8 | 3.95 | 3.15E-17 | 3.23E-17
Table 12: Test 7.5: 2D Euler equation with gravity and gravitational potential
$H_{2}$. Errors in $L^{1}$ norm for density at time $t=0.3.$

We solve numerically the equations using a $(21\times 21)-$point mesh and
CFL$=0.9$. As boundary condition the exact solution is imposed to all sides
through the ghost points. Tables 11 and 12 exhibit the errors in $L^{1}-$norm
for ACAT$2P$, WBACAT$2P$, $P=1,2$ corresponding to $H=H_{1}$ and $H=H_{2}$
respectively. As it can be seen, the differences between the solutions given
by well-balanced and no well-balanced methods are bigger for $H=H_{2}$: please
note that, in the case of the linear potential $H_{1},$ the stationary
solution is essentially 1D while this is not true in the case $H=H_{2}.$

###  7.6 Perturbation of the stationary solution

We consider now Euler equations in the 2D domain $[0,1]\times[0,1]$ with the
gravitational potential $H_{2}$ and an initial condition that represents a
perturbation of the hydrostatic stationary considered in the previous test
case:

$\displaystyle\rho(\mathbf{x},0)$
$\displaystyle=e^{-H(\mathbf{x})}+0.008e^{-200(x-0.5)^{2}-200(y-0.5)^{2}};\quad
p(\mathbf{x},0)=e^{-H(\mathbf{x})}+0.008e^{-200(x-0.5)^{2}-200(y-0.5)^{2}};$
$\displaystyle u(\mathbf{x},0)$ $\displaystyle=v(\mathbf{x},0)=0.$ (7.16)

| 2D density
---|---
| 2D ACAT2 | 2D ACAT4 | 2D WBACAT2 | 2D WBACAT4
Points | Error | Order | Error | Order | Error | Order | Error | Order
20$\times$20 | 4.49E-5 | - | 5.93E-6 | - | 8.27E-6 | - | 4.28E-7 | -
40$\times$40 | 2.47E-5 | 0.86 | 1.71E-6 | 1.79 | 4.41E-6 | 0.91 | 4.92E-8 | 2.71
80$\times$80 | 1.21E-5 | 1.03 | 4.37E-7 | 1.97 | 2.08E-6 | 1.08 | 7.05E-9 | 2.80
160$\times$160 | 5.45E-6 | 1.15 | 9.85E-8 | 2.15 | 8.13E-7 | 1.36 | 9.47E-10 | 2.90
320$\times$320 | 2.43E-6 | 1.17 | 2.13E-8 | 2.21 | 2.45E-7 | 1.73 | 1.22E-10 | 2.96
Table 13: Test 7.6: 2D Euler equation with gravity and gravitational potential
$H_{2}$. Errors in $L^{1}$ norm for density at time $t=0.2.$

Table 13 shows error in $L^{1}-$norm and convergence rates for the numerical
solutions obtained with ACAT$2P$ and WBACAT$2P$ and the reference solution at
time $t=0.2,$ with $P=1,2.$ As happened for Shallow water, in case that a
small perturbation of the stationary solution is considered as initial
condition, the well-balanced schemes manage to capture the solution with a
better accuracy than standard methods. This phenomena is shown on Tables 13.

###  7.7 Acoustic propagation

As last experiment we consider the Euler equations in the 2D domain
$[0,2]\times[0,2]$ with the gravitational potential $H_{3},$

$H_{3}(x,y)=\dfrac{1}{\sqrt{(x-0.4)^{2}+(y+0.1)^{2}}},$

and an initial condition that represents a very small perturbation of the
hydrostatic stationary considered in the previous test case:

$\displaystyle\rho(\mathbf{x},0)=e^{-H(\mathbf{x})}+0.000001e^{-200(x-1)^{2}-200(y-1)^{2}};\quad
p(\mathbf{x},0)=\rho(\mathbf{x},0);\quad u(\mathbf{x},0)=v(\mathbf{x},0)=0.$
(7.17)

Figure 7.16: Test 7.7: Euler equations with gravitational potential $H_{3}$.
Initial perturbation using a $101\times 101$ mesh points. Figure 7.17: Test
7.7: Euler equations with gravitational potential $H_{3}$. Difference between
the numerical solutions and the stationary solution computed at time $t=0.75$
with WBACAT2 using $101\times 101$ mesh points and CFL$=0.8.$ Figure 7.18:
Test 7.7: Euler equations with gravitational potential $H_{3}$. Difference
between the density solutions and the stationary solution computed at time
$t=0.75$ with ACAT2 using $101\times 101$ mesh points and CFL$=0.8.$

Figures 7.17 shows the difference between the numerical solutions and the
stationary solution computed at time $t=0.75$ with WBACAT2 a using $101\times
101$ mesh points and CFL$=0.8.$ As expected, the singularity of the
gravitational potential modifies the thickness of the corona relative to the
signal propagation, thinning it close to the singularity. $H_{3}$ has a
singularity on position $(0.4,-0.1).$
As we can see in Figure 7.18, the non well-balanced method ACAT2 is not able
to capture the evolution of the wave generated by the initial perturbation,
since the numerical errors are much bigger than the wave amplitude: one would
need a space step at least two orders of magnitude lower in order to have a
truncation error of the same order of the signal, making computation with non
well-balanced method absolutely impractical.

##  8 Conclusion

This paper deals with the construction, analysis, implementation and testing
of well-balanced
ACAT2$P$ methods to the treatment of hyperbolic systems of balance laws. The
starting point is to re-write the systems as conservation laws, by subtracting
to the flux a primitive of the the source term. Meanwhile, the well-balanced
property has been obtained rewriting the systems as conservation laws, by
subtracting to the flux the flux of the stationary solution, and adding to the
primitive of the source the source computed at stationary solution. The
methods are developed for systems in one and two space dimensions, and can be
extended to 3D. In principle the procedure allows the construction of well-
balanced schemes of arbitrary order, although the computational complexity
quickly increases with the order of accuracy. We prove that the constructed
schemes are exactly well-balanced.

The use of suitable limiters allow an effective treatment of discontinuous
solutions. Several test cases have been performed for the scalar equation and
for systems in one and two space dimensions. In all cases we observe that
stationary solutions are preserved within machine precision, allowing very
accurate results when the solution is a small deviation from equilibrium.

The main advantage of the method consists in its generality: it allows the
automatic construction of very high order well-balanced schemes.

There are still a few things that require improvement and generalization.
First, a careful complexity analysis of the methods is needed, together with
an improvement of the computational efficiency, possibly by parallel computing
techniques, which should be possible thanks to the local computations required
by these methods.

Second, we shall explore new limiting strategies, which will allow more
accurate results when the solution is smooth.

An open problem is how to couple ACAT methodology with IMEX method for the
treatment of problems with stiff source.

All these issues are subject of current investigation.

## Acknowledgements

This research has received funding from the European Union’s Horizon 2020
research and innovation program, under the Marie Sklodowska-Curie grant
agreement No 642768. E. Macca was partially supported by GNCS Research Project
"Approssimazione numerica di problemi di natura iperbolica ed applicazioni".
E. Macca and G.Russo would like to thank the Italian Ministry of Instruction,
University and Research (MIUR) to support this research with funds coming from
PRIN Project 2017 (No. 2017KKJP4X entitled “Innovative numerical methods for
evolutionary partial differential equations and applications”). The research
of C. Parés was partially supported by the Spanish Government(SG), the
European Regional Development Fund(ERDF), the Regional Government of
Andalusia(RGA), and the University of Málaga(UMA) through the projects of
reference RTI2018-096064-B-C21 (SG-ERDF), UMA18-Federja-161 (RGA-ERDF-UMA),
and P18-RT-3163 (RGA-ERDF). E. Macca and G. Russo are members of the INdAM
Research group GNCS.

## References

  * [1] Bermúdez, A., and Vázquez, M. E. Upwind methods for hyperbolic conservation laws with source terms.  Computers & Fluids 23(8) (1994), 1049–1071.
  * [2] Bouchut, F. Non-linear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources.  Frontiers in Mathematics Birkhauser (2004).
  * [3] Carrillo, H., Macca, E., Parés, C., Russo, G., and Zorío, D. An order-adaptive compact approximate taylor method for systems of conservation law.  Journal of Computational Physics 438 (2021), 31.
  * [4] Carrillo, H., and Parés, C. Compact Approximate Taylor methods for systems of conservation laws.  Journal of Scientific Computing 80 (2019), 1832–1866.
  * [5] Carrillo, H., Parés, C., and Zorío, D. Lax-Wendroff Approximate Taylor methods with fast and optimized weighted essentially non-oscillatory reconstructions.  J. Sci. Comput. 86 (2021), 41.
  * [6] Caselles, V., Donat, R., and Haro, G. Flux-gradient and source-term balancing for certain high resolution shock-capturing schemes.  Computers & Fluids 38 (2009), 16–36.
  * [7] Castro, M., López-García, J., and Parés, C. High order exactly well-balanced numerical methods for shallow water systems.  Journal of Computational Physics 246 (2013), 242–264.
  * [8] Castro, M., and Parés, C. Well-balanced high-order finite volume methods for systems of balance laws.  J. Sci. Comput. 82 (2020), 48.
  * [9] Castro, M. J., Gómez-Bueno, I., and Parés, C. High-order well-balanced methods for systems of balance laws: a control-based approach.  Appl. Math. Comput. 394 (2021), 31.
  * [10] Donat, R., and Martínez-Gavera, A. Hybrid second order schemes for scalar balance laws.  Journal of Scientific Computing 48 (2011), 52–69.
  * [11] Dumbser, M., Balsara, D., Toro, E., and Munz, C. A unified framework for the construction of one-step finite-volume and discontinuous Galerkin schemes.  Journal of Computational Physics 227 (2008), 8209–8253.
  * [12] Dumbser, M., Enaux, C., and Toro, E. F. Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws.  J. Comput. Phys. 227 (2008), 3971–4001.
  * [13] Fornberg, B. Generation of finite difference formulas on arbitrarily space grids.  Mathematics of Computation 51 (1988), 699–706.
  * [14] Gascón, L., and Corderán, J. M. Construction of second-order TVD schemes for nonhomogeneous hyperbolic conservation laws.  Journal of Computational Physics 172 (2001), 261–297.
  * [15] Gideon, R. A. Some methods of approximating distribution functions.  Thesis (Ph.D.)–The University of Wisconsin - Madison (1971), 323.
  * [16] Grosheintz-Laval, L., and Kappeli, R. High-order well-balanced finite volume schemes for the Euler equations with gravitation.  Journal of Computational Physics 378 (2019), 324–343.
  * [17] Grosheintz-Laval, L., and Kappeli, R. Well-balanced finite volume schemes for nearly steady adiabatic flows.  Journal of Computational Physics 423 (2020), 28.
  * [18] Kappeli, R., and S.Mishra.  Well-balanced schemes for the Euler equations with gravitation.  Journal of Computational Physics 259 (2014), 199–219.
  * [19] Kemm, F. A comparative study of tvd-limiters - well-known limiters and an introduction of new ones.  Internationale Journal of Numerical Methods in Fluids 67 (2010), 404–440.
  * [20] Klingenberg, C., Puppo, G., and Semplice, M. Arbitrary order finite volume well-balanced schemes for the Euler equations with gravity.  SIAM J. Sci. Comput. 41(2) (2019), 695–721.
  * [21] LeVeque, R. Finite volume methods for hyperbolic problems, 1 ed.  Cambridge Texts in Applied Mathematics. Cambridge University Press., 2002\.
  * [22] LeVeque, R. Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems (Classics in Applied Mathematics), 1 ed.  Society for Industrial and Applied Mathematics, Philadelpia, PA. USA., 2007.
  * [23] Millington, R. C., Titarev, V. A., and Toro, E. F. Ader: arbitrary-order non-oscillatory advection schemes. hyperbolic problems: theory, numerics, applications.  Internat. Ser. Numer. Math., 140,141, Birkhauser, Basel I (2001), 723–732.
  * [24] Parés, C., and Parés-Pulido, C. Well-balanced high-order finite difference methods for systems of balance laws.  J. Comput. Phys. 45 (2021), 35.
  * [25] Roe, P. Characteristic-based schemes for Euler equations.  Annu. Rev. Fluid Mech. 18 (1986), 337–365.
  * [26] Russo, G. Deterministic diffusion of particles.  Communications on Pure and Applied Mathematics 43, 6 (1990), 697–733.
  * [27] Russo, G. A particle method for collisional kinetic equations. i. basic theory and one-dimensional results.  Journal of Computational Physics 87, 2 (1990), 270–300.
  * [28] Russo, G., and Khe, A. High order well-balanced schemes based on numerical reconstruction of the equilibrium variables.  In Proc. of the 15th Internatinal Conference on Waves and Stability on Continuous Media (2010), World Scientific, pp. 230–241.
  * [29] Schwartzkopff, T., Munz, C., and Toro, E. F. Ader: a high-order approach for linear hyperbolic systems in 2d. proceedings of the fifth international conference on spectral and high order methods (icosahom-01) (uppsala).  J. Sci. Comput. 17 (2002), 231–240.
  * [30] Sweby, P. K. High resolution schemes using flux limiters for hyperbolic conservation laws.  SIAM J. Numer. Anal. 21 (5) (1984), 995–1011.
  * [31] Titarev, V. A., and Toro, E. T. Ader: arbitrary high order godunov approach. proceedings of the fifth international conference on spectral and high order methods (icosahom-01) (uppsala).  J. Sci. Comput. 17 (2002), 609–618.
  * [32] Toro, E. Riemann Solvers and Numerical Methods for Fluid Dynamics, third ed.  Springer, 2009.
  * [33] Xing, Y. Numerical methods for the nonlinear Shallow water equations.  Handbook of Numerical Methods, Springer 18 (2017), 361–384.
  * [34] Xing, Y., and Shu, C. W. High-order well-balanced finite difference weno schemes for a class of hyperbolic systems with source terms.  Journal of Scientific Computing 27 (2006), 477–494.
  * [35] Xing, Y., and Shu, C. W. High-order well-balanced finite difference weno schemes with the exact conservation property for the shallow water equations.  Journal of Computational Physics 208 (2006), 206–227.
  * [36] Zorío, D., Baeza, A., and Mulet, P. An approximate Lax-Wendroff-type procedure for high order accurate scheme for hyperbolic conservation laws.  J. Sci. Comput. 71, 1 (2017), 246–273.