∎

11institutetext: Jeremy Chouchoulis , Jochen Schuetz , Jonas Zeifang
11email<EMAIL_ADDRESS>
22institutetext: Hasselt University, Faculty of Sciences & Data Science
Institute, Agoralaan Gebouw D, BE-3590 Diepenbeek, Belgium

# Jacobian-free explicit multiderivative Runge-Kutta methods for hyperbolic
conservation laws

Jeremy Chouchoulis  Jochen Schütz  Jonas Zeifang  0000-0003-3451-2568
0000-0002-6355-9130 0000-0002-7388-923X

(Received: date / Accepted: date)

###### Abstract

Based on the recent development of Jacobian-free Lax-Wendroff (LW) approaches
for solving hyperbolic conservation laws [Zorio, Baeza and Mulet, Journal of
Scientific Computing 71:246-273, 2017], [Carrillo and Parés, Journal of
Scientific Computing 80:1832-1866, 2019], a novel collection of explicit
Jacobian-free multistage multiderivative solvers for hyperbolic conservation
laws is presented in this work. In contrast to Taylor time-integration
methods, multiderivative Runge-Kutta (MDRK) techniques achieve higher-order of
consistency not only through the excessive addition of higher temporal
derivatives, but also through the addition of Runge-Kutta-type stages. This
adds more flexibility to the time integration in such a way that more stable
and more efficient schemes could be identified. The novel method permits the
practical application of MDRK schemes. In their original form, they are
difficult to utilize as higher-order flux derivatives have to be computed
analytically. Here we overcome this by adopting a Jacobian-free approximation
of those derivatives. In this paper, we analyze the novel method with respect
to order of consistency and stability. We show that the linear CFL number
varies significantly with the number of derivatives used. Results are verified
numerically on several representative testcases.

###### Keywords:

Hyperbolic conservation laws Multiderivative Runge-Kutta Lax-Wendroff Finite
differences

###### MSC:

65M06 65M08 65M12 35L65

##  1 Introduction

In this work, we present a novel discretization method for the numerical
approximation of one-dimensional hyperbolic conservation laws on domain
$\Omega\subset\mathbb{R}$,

$\displaystyle w_{t}+{f}({w})_{x}$
$\displaystyle=0\,,\qquad\qquad\text{on}\quad(x,t)\in\Omega\times(0,T_{\text{end}}]\,,$
(1) $\displaystyle{w}(x,0)$ $\displaystyle={w}_{0}(x)\,.$

Our primary interest is on temporal integration. In recent years, there has
been quite some progress on the further development of the _multiderivative_
paradigm for temporal integration, see, e.g., TC10 ; Seal2015b ;
2021_Gottlieb_EtAl ; SSJ2017 ; Seal13  and the references therein. Assume that
one is given a scalar ODE, e.g.,

$\displaystyle{y}^{\prime}(t)$ $\displaystyle=\Phi({y})$ (2)

for some flux function $\Phi$. Multiderivative schemes make use of not only
$\Phi$, but also of the quantities
$y^{\prime\prime}(t)\equiv\Phi^{\prime}(y)\Phi$,
$y^{\prime\prime\prime}(t)\equiv\ldots$ and so on. Using this approach, one
can derive stable, high-order and storage-efficient schemes very easily
SealSchuetzZeifang21 . This can be extended to partial differential equations
(PDEs) with a time-component, such as Eq. (1), depending on the method either
directly through the method of lines-discretization SSJ2017  or through a Lax-
Wendroff procedure, see, e.g., CarrilloPares2019 ; CarrilloParesZorio2021 ;
LaxWend1960 ; LiDu2016 ; Qiu08 ; ZorioEtAl . The Lax-Wendroff method expresses
temporal derivatives of the unknown function $w$ in terms of the fluxes
through the Cauchy-Kowalevskaya procedure. As an example, we consider – for
simplicity given that $f$ is scalar – the second time-derivative of $w$. Due
to Eq. (1), there holds

$\displaystyle w_{tt}=-(f(w)_{x})_{t}=-(f(w)_{t})_{x},$ (3)

and

$\displaystyle f(w)_{t}=f^{\prime}(w)w_{t}=-f^{\prime}(w)f(w)_{x};$

hence

$\displaystyle
w_{tt}=\left(f^{\prime}(w)f(w)_{x}\right)_{x}=2f^{\prime}(w)f^{\prime\prime}(w)w_{x}^{2}+f^{\prime}(w)^{2}w_{xx}.$
(4)

Already at this stage, one can see that this approach is very tedious as it
necessitates highly complex symbolic calculations.

Still, the potential LW-methods bear is very well recognized among
researchers. Over the last two decades, plenty of authors have put effort into
developing high-order variants of the LW-method for nonlinear systems.
Particularly the ADER (Arbitrary order using DERivatives) methods, see, e.g.,
PnPm0 ; 2008_Dumbser_Enaux_Toro ; axioms7030063 ; Schwartzkopff2002 ;
Titarev2002 ; 2005_Titarev_Toro  and the references therein, gained a lot of
interest. Also, higher-order extensions of the LW-method using WENO and
discontinuous Galerkin (DG) reconstructions were investigated GuoQiuQiu15 ;
MultiDerHDG2015 ; 2011_Lu_Qiu ; Qiu08 ; Qiu2005 ; QiuShu03 .

Our essential intent of this paper is to make explicit multistage
multiderivative solvers more accessible as a means to solve PDEs. Although
such solvers have been theoretically studied since the early 1940’s (see
Seal13  for an extensive review), the schemes have not been put much to
practice, which is most likely due to the necessary cumbersome calculation of
flux derivatives. In CarrilloPares2019 , Carrillo and Parés have, based on the
earlier work ZorioEtAl , developed the _compact approximate Taylor_ (CAT)
method to circumvent having to symbolically compute flux derivatives. By means
of an automatic procedure, the higher-order temporal derivatives of $w$, such
as in Eq. (4), are approximated. Their work is based on Taylor methods, i.e.,
time integration is given by

${w}(x,t^{n+1})={w}(x,t^{n})+\sum\limits_{k=1}^{\mathtt{r}}\frac{{\Delta
t}^{k}}{k!}\partial^{k}_{t}{w}(x,t^{n})+\mathcal{O}({\Delta
t}^{\mathtt{r+1}})\,.$

In this work, we extend their approach to more general multiderivative
integration methods, more precisely, to multiderivative Runge-Kutta (MDRK)
methods.

The paper is structured in the following manner: In Sect. 2 multiderivative
Runge-Kutta (MDRK) time integrators for ODEs are introduced, given that they
form the central mechanism of this work. Thereafter, in Sect. 3 we shortly
revisit the Jacobian-free approach of the CAT method and introduce the
explicit Jacobian-free MDRK solver for hyperbolic conservation laws, termed
MDRKCAT. After describing the numerical scheme, in Sect. 4 we prove
consistency, and in Sect. 5 analyze linear stability. Via several numerical
cases we verify and expand on the theoretical results in Sect. 6. At last, we
draw our conclusions and discuss future perspectives in Sect. 7.

##  2 Explicit multiderivative Runge-Kutta solvers

We start by considering the system of ODEs defined by Eq. (2) in which $\Phi$
is a function of the solution variable ${y}\in\mathbb{R}^{\mathtt{m}}$. In
order to apply a time-marching scheme, we discretize the temporal domain with
a fixed timestep $\Delta t$ by iterating $N$ amount of steps such that $\Delta
t=T_{\text{end}}/N$. Consequently, we define the time levels by

$t^{n}:=n\Delta t\,\qquad 0\leq n\leq N.$

###### Remark 1

Note that, although the fully space-time-discrete algorithm (Alg. 3.2) _seems_
to have a multistep flavour, this is ultimately not the case. It is therefore
of no necessity to consider a uniform timestep, which is also demonstrated
numerically in Sect. 6.

The central class of time integrators in this work are _explicit_
multiderivative Runge-Kutta (MDRK) methods. These form a natural
generalization of classical explicit Runge-Kutta methods by adding extra
temporal derivatives of $\Phi(w)$. The additional time derivatives can be
recursively calculated via the chain rule, there holds

$\frac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}{\Phi}\left({y}\right)=\frac{\mathrm{d}^{k-1}}{\mathrm{d}t^{k-1}}\left(\Phi^{\prime}(y)\dot{y}\right)=\frac{\mathrm{d}^{k-1}}{\mathrm{d}t^{k-1}}\left(\Phi^{\prime}(y)\Phi(y)\right).$

For a more detailed description, we refer to Seal13 . To present our ideas,
let us formally define the MDRK scheme as follows:

######  Definition 1 ((Seal13, , Def. 2))

An explicit $q$-th order accurate $\mathtt{r}$-derivative Runge-Kutta scheme
using $\mathtt{s}$ stages ($\mathtt{r}\text{DRK}q\text{-}\mathtt{s}$) is any
method which can be formalized as

$\hskip
64.00003pt{y}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}:={y}^{n}+\sum\limits_{k=1}^{\mathtt{r}}{\Delta
t}^{k}\sum\limits_{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}=1}^{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1}a_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}^{(k)}\frac{\mathrm{d}^{k-1}}{\mathrm{d}t^{k-1}}{\Phi}\left({y}^{n,{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}\right)\quad{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1,\dots,\mathtt{s},$
(5) where
${y}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>is a stage approximation at time
$t^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}:=t^{n}+c_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>t$. The update is given by
${y}^{n+1}:={y}^{n}+\sum\limits_{k=1}^{\mathtt{r}}{\Delta
t}^{k}\sum\limits_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1}^{\mathtt{s}}b^{(k)}_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}\frac{\mathrm{d}^{k-1}}{\mathrm{d}t^{k-1}}{\Phi}({y}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}})\,.$

The given coefficients
$a_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}^{(k)}$
and
$b^{(k)}_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>determine the scheme; they are typically summarized in an extended Butcher
tableau.

###### Remark 2

Note that standard Taylor methods can be cast in the framework of Def. 1
through setting $\mathtt{s}=1$, $a_{11}^{(k)}=0$ and $b_{1}^{(k)}=1/k!$
$(k=1,\dots,\mathtt{r})$. The multiderivative Runge-Kutta schemes used in this
work can be found through their extended Butcher tableaux in Appendix A.

###### Remark 3

The stability regions of the used Runge-Kutta methods are visualized in Fig.
1, see TC10 ; HaiWan ; TurTur2017  for more details. Note that except for
$3\text{DRK}5\text{-}2$ and $4\text{DRK}6\text{-}2$, all schemes contain parts
of the imaginary axis.

Scheme | $R(z)$
---|---
$2\text{DRK}3\text{-}2$ | $\sum_{k=0}^{3}\frac{z^{k}}{k!}$
$2\text{DRK}4\text{-}2$ | $\sum_{k=0}^{4}\frac{z^{k}}{k!}$
$2\text{DRK}5\text{-}3$ | $\sum_{k=0}^{5}\frac{z^{k}}{k!}+\frac{z^{6}}{600}$
$3\text{DRK}5\text{-}2$ | $\sum_{k=0}^{5}\frac{z^{k}}{k!}+\frac{z^{6}}{900}$
$3\text{DRK}7\text{-}3$ | $\sum_{k=0}^{7}\frac{z^{k}}{k!}+c_{8}z^{8}+c_{9}z^{9}$
$4\text{DRK}6\text{-}2$ | $\sum_{k=0}^{6}\frac{z^{k}}{k!}+\frac{z^{7}}{6480}+\frac{z^{8}}{77760}$

Figure 1: Regions of absolute stability
$\mathcal{R}:=\left\\{z\in\mathbb{C}~{}|~{}|R(z)|\leq 1\right\\}$ and
stability functions $R(z)$ for the two-derivative schemes
$2\text{DRK}3\text{-}2$, $2\text{DRK}4\text{-}2$, $2\text{DRK}5\text{-}3$ TC10
, the three-derivative schemes $3\text{DRK}5\text{-}2$ and
$3\text{DRK}7\text{-}3$ TurTur2017  and the four-derivative scheme
$4\text{DRK}6\text{-}2$ (see Appendix A). Except for $3\text{DRK}5\text{-}2$
and $4\text{DRK}6\text{-}2$, all schemes contain parts of the imaginary axis.
Note that in the stability function of $3\text{DRK}7\text{-}3$, we have
defined $c_{8}:=\frac{1}{23520}-\frac{\sqrt{2}}{70560}$ and
$c_{9}:=\frac{11}{1481760}-\frac{\sqrt{2}}{246960}$.

##  3 Multiderivative Runge-Kutta solvers for hyperbolic conservation laws

Discretizing the spatial part of the hyperbolic conservation law (1)
necessitates a discretization of the domain $\Omega$. Hence, consider

$\displaystyle\\{x_{1},\ldots,x_{M}\\}$

to be a uniform partition of $\Omega$ into $M$ cells of size $\Delta x$. A
natural extension of Def. 1 applied to Eq. (1) can then be expressed as

$\displaystyle
w_{i}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>$\displaystyle:=w_{i}^{n}-\sum\limits_{k=1}^{\mathtt{r}}{\Delta
t}^{k}\sum\limits_{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}=1}^{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1}a_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}^{(k)}D_{x}D^{k-1}_{t}f(w_{i}^{n,{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}})\,,$
(6a) $\displaystyle w_{i}^{n+1}$
$\displaystyle:=w_{i}^{n}-\sum\limits_{k=1}^{\mathtt{r}}{\Delta
t}^{k}\sum\limits_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1}^{\mathtt{s}}b^{(k)}_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}D_{x}D^{k-1}_{t}f(w_{i}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}),$
(6b)

for
${\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1,\dots,\mathtt{s}$;
with $D_{x}$ and $D_{t}$ being suitable approximations to $\partial_{x}$ and
$\partial_{t}$ to be explained in the sequel. Contrary to the complete Cauchy-
Kovalevskaya procedure as outlined for the second derivative in (4), only
_one_ time derivative of the solution is transformed into a spatial derivative
(see (3)).

The core focus of this paper is to avoid the explicit use of Jacobians of the
flux function $f$. Jacobians of $f$ arise due to the usage of higher temporal
derivatives, see, e.g., Eq. (4). In this, we follow the compact approximate
Taylor (CAT) approach outlined in CarrilloPares2019 . Since the CAT method
heavily relies on discrete differentiation, first a small part is devoted to
introducing the fundamental notation. Thereafter the method is described and
applied to Eqs. (6).

###  3.1 Discrete differentiation

In this short section, we fix the notation on using finite differencing
1964_Abramowitz_Stegun ; QuarteroniNumA ; 2001_Ralston_Rabinowitz .
Considering central differences, the $(2p+1)$-point Lagrangian polynomials are
given by

$\hskip 40.00006ptL_{p,j}(\omega):=\prod\limits_{\begin{subarray}{c}r=-p\\\
r\neq j\end{subarray}}^{p}\dfrac{\omega-r}{j-r}\qquad j=-p,\dots,p\,.$ (7)

It is well-known that these polynomials can be used to interpolate
$\varphi:\mathbb{R}\to\mathbb{R}$ in the points $x_{i-p},\dots,x_{i+p}$
through

$\mathcal{P}_{i}\varphi(x):=\sum\limits_{j=-p}^{p}L_{p,j}\left(\frac{x-x_{i}}{\Delta
x}\right)\varphi(x_{i+j}).$ (8)

Similarly, from the $2p$-point Lagrangian polynomials

$\hskip
40.00006pt\ell_{p,j}(\omega):=\prod\limits_{\begin{subarray}{c}r=-p+1\\\ r\neq
j\end{subarray}}^{p}\dfrac{\omega-r}{j-r}\qquad j=-p+1,\dots,p\,,$ (9)

(note that the index of the product begins at $r=-p+1$) we obtain the unique
polynomial of degree $2p-1$ interpolating $\varphi$ in the points
$x_{i-p+1},\dots,x_{i+p}$ through

$\mathcal{Q}_{i}\varphi(x):=\sum\limits_{j=-p+1}^{p}\ell_{p,j}\left(\frac{x-x_{i}}{\Delta
x}\right)\varphi(x_{i+j}).$ (10)

In the sequel, we use a similar notation as in CarrilloPares2019 :

######  Definition 2 (CarrilloPares2019 )

For the $k$-th derivative $0\leq k\leq 2p$, we define the following
quantities:

$\displaystyle\delta^{k}_{p,j}$ $\displaystyle:=L_{p,j}^{(k)}(0),\qquad j$
$\displaystyle=\ $ $\displaystyle-p,\ldots,p,$
$\displaystyle\gamma^{k,m}_{p,j}$ $\displaystyle:=\ell_{p,j}^{(k)}(m),\qquad
j$ $\displaystyle\ =\ $ $\displaystyle-p+1,\ldots,p,\quad m=-p+1,\ldots,p.$

Note that the $k$-th derivatives $L_{p,j}^{(k)}(0)$ and $\ell_{p,j}^{(k)}(m)$
are derived analytically from Eq. (7) and Eq. (9), respectively.

Approximate derivatives can thus be derived from

$\displaystyle(\mathcal{P}_{i}\varphi)^{(k)}(x_{i})$
$\displaystyle=\frac{1}{{\Delta
x}^{k}}\sum\limits_{j=-p}^{p}\delta^{k}_{p,j}\varphi(x_{i+j})\,.$ (11a)
$\displaystyle(\mathcal{Q}_{i}\varphi)^{(k)}(x_{i+m})$
$\displaystyle=\frac{1}{{\Delta
x}^{k}}\sum\limits_{j=-p+1}^{p}\gamma^{k,m}_{p,j}\varphi(x_{i+j})\,,$ (11b)

with $m=-p+1,\dots,p$. Since we are working in a discrete context, we define
the linear operator counterparts of Eq. (11a) and Eq. (11b) as

$\displaystyle P^{(k)}\colon\mathbb{R}^{2p+1}$
$\displaystyle\rightarrow\mathbb{R},$
$\displaystyle\mathbf{v}\mapsto\frac{1}{{\Delta
x}^{k}}\sum\limits_{j=-p}^{p}\delta^{k}_{p,j}v_{j},$ $\displaystyle
Q_{m}^{(k)}\colon\mathbb{R}^{2p}$ $\displaystyle\rightarrow\mathbb{R},$
$\displaystyle\mathbf{w}\mapsto\frac{1}{{\Delta
x}^{k}}\sum\limits_{j=-p+1}^{p}\gamma^{k,m}_{p,j}w_{j}\,.$

###### Remark 4

The spatial index $i$ is neglected for the linear operators as its direct
dependency on the node $x_{i}$ is lost, cf. Eq. (8) and Eq. (10). Notice also
that $Q_{m}^{(k)}$ takes vectors in $\mathbb{R}^{2p}$ as input, whereas
$P^{(k)}$ takes vectors in $\mathbb{R}^{2p+1}$ as input.

A non-centered $2p$-point finite difference method to approximate
$\partial^{k}_{t}w(x_{i},t^{n+m})$ for $m=-p+1,\dots,p$ can therefore be
written as

$Q_{m}^{(k)}\mathbf{w}^{\langle n\rangle}_{i}=\frac{1}{{\Delta
t}^{k}}\sum\limits_{r=-p+1}^{p}\gamma^{k,m}_{p,j}w_{i}^{n+r}\,,$

with vector notation

$\mathbf{w}^{\langle n\rangle}_{i}:=\begin{pmatrix}w_{i}^{n-p+1}\\\ \vdots\\\
w_{i}^{n+p}\end{pmatrix}\,.$ (12)

The angled brackets represent the local stencil function

$\langle\cdot\rangle\colon\mathbb{Z}\to\mathbb{Z}^{2p}\colon
n\mapsto\begin{pmatrix}n-p+1,&\dots,&n+p\end{pmatrix}^{T}$ (13)

throughout this paper, and will be considered for both the spatial index $i$
as the temporal index $n$. Note that the position of the angled bracket (top
or bottom) determines whether derivation is w.r.t. time (top) or space
(bottom).

To put the scheme into conservation form, in ZorioEtAl  auxiliary centered
coefficients have been introduced. Here, the operators $P^{(k)}$ for $k\geq 1$
are written as differences of new ‘half-way point’ interpolation operators.

######  Definition 3 (ZorioEtAl )

Define $\lambda^{k-1}_{p,j}$ via the relations

$\displaystyle\delta^{k}_{p,p}$ $\displaystyle=:\lambda^{k-1}_{p,p}\,,$ (14a)
$\displaystyle\delta^{k}_{p,j}$
$\displaystyle=:\lambda^{k-1}_{p,j}-\lambda^{k-1}_{p,j+1}\,,\quad
j=-p+1,\dots,p-1\,,$ (14b) $\displaystyle\delta^{k}_{p,-p}$
$\displaystyle=:-\lambda^{k-1}_{p,-p+1}\,.$ (14c)

###### Remark 5

The relations given in Eq. (14) make up an overdetermined system, yet provide
a unique solution $\lambda^{k-1}_{p,j}$ obtained from Eq. (14a) and Eq. (14b),
see (ZorioEtAl, , Theorem 2).

Notice the shift between $k$ and $k-1$. This is justified because we enforce a
first order derivative relation for the approximation
$(\mathcal{P}_{i}\varphi)^{(k)}(x_{i})$ by splitting the operator as

$(\mathcal{P}_{i}\varphi)^{(k)}(x_{i})=\frac{\left(\Uplambda^{(k-1)}\varphi\right)(x_{i+1/2})-\left(\Uplambda^{(k-1)}\varphi\right)(x_{i-1/2})}{\Delta
x}\,,$ (15)

in which $\Uplambda^{(k-1)}$ is an operator mapping to $\mathbb{P}_{2p-1}$ so
that

$\left(\Uplambda^{(k-1)}\varphi\right)(x_{i+1/2}):=\frac{1}{{\Delta
x}^{k-1}}\sum\limits_{j=-p+1}^{p}\lambda^{k-1}_{p,j}\varphi(x_{i+j})\,.$

The linear operator alternative is defined by

$\displaystyle\Lambda^{(k-1)}\colon\mathbb{R}^{2p}$
$\displaystyle\longrightarrow\mathbb{R},\qquad\mathbf{v}\mapsto\frac{1}{{\Delta
x}^{k-1}}\sum\limits_{j=-p+1}^{p}\lambda^{k-1}_{p,j}v_{j}.$ (16)

An overview of all the defined interpolation operators is given in Tbl. 1.

Table 1: A summary of the defined interpolation operators scaled and shifted to fit the uniform mesh locally at $x_{i}$.  Functional | Linear |
---|---|---
$(\mathcal{P}_{i}\varphi)^{(k)}$ | $P^{(k)}\mathbf{v}$ |  $k$-th derivative of the Lagrangian interpolation
|  | polynomial in the nodes $x_{i-p},\dots,x_{i+p}$
$(\mathcal{Q}_{i}\varphi)^{(k)}$ | $Q_{m}^{(k)}\mathbf{v}$ |  $k$-th derivative of the Lagrangian interpolation
|  | polynomial in the nodes $x_{i-p+1},\dots,x_{i+p}$
$\left(\Uplambda^{(k-1)}\varphi\right)(x_{i+1/2})$ | $\Lambda^{(k-1)}\mathbf{v}$ |  $(k-1)$-th derivative of a half-way interpolation
|  | at $x_{i+1/2}$ using the nodes $x_{i-p+1},\dots,x_{i+p}$
|  | with the difference coefficients defined by Eq. (14)

###  3.2 A Jacobian-free MDRK scheme

With all the building blocks at our disposal, we can now describe how the
final class of methods, that we call MDRKCAT, is assembled. Starting from Eq.
(6), we define the conservative updates of the solution via

$\displaystyle
w^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>$\displaystyle:=w^{n}_{i}-\frac{\Delta t}{\Delta
x}\big{(}F_{i+1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}-F_{i-1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}\big{)}\qquad{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1,\dots,\mathtt{s}\,,$
(17a) $\displaystyle w^{n+1}_{i}$ $\displaystyle:=w^{n}_{i}-\frac{\Delta
t}{\Delta x}\big{(}F_{i+1/2}^{n}-F_{i-1/2}^{n}\big{)}\,,$ (17b)

in which the numerical fluxes are given by,

$\displaystyle
F_{i+1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>$\displaystyle=\sum\limits_{k=1}^{\mathtt{r}}{\Delta
t}^{k-1}\sum\limits_{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}=1}^{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1}a_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}^{(k)}\Lambda^{(0)}(\widetilde{\mathbf{f}}^{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}})^{(k-1)}_{i,\langle
0\rangle}\qquad{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1,\dots,\mathtt{s}\,,$
(18a) $\displaystyle F_{i+1/2}^{n}$
$\displaystyle=\sum\limits_{k=1}^{\mathtt{r}}{\Delta
t}^{k-1}\sum\limits_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1}^{\mathtt{s}}b^{(k)}_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}\Lambda^{(0)}(\widetilde{\mathbf{f}}^{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l})^{(k-1)}_{i,\langle
0\rangle}\,.$ (18b)

For the calculation of $\widetilde{\mathbf{f}}^{(k-1)}_{i,\langle 0\rangle}$
the compact approximate Taylor (CAT) procedure CarrilloPares2019  is used and
the flux derivatives can be calculated according to Eq. (16) by

$\Lambda^{(0)}\widetilde{\mathbf{f}}^{(k-1)}_{i,\langle
0\rangle}:=\sum\limits_{j=-p+1}^{p}\lambda^{0}_{p,j}\widetilde{f}^{(k-1)}_{i,j}\,.$

$\widetilde{f}^{(k-1)}_{i,j}\approx\partial^{k-1}_{t}f(w)(x_{i+j},t^{n})$
indicates the local approximations for the time-derivatives of the flux and
are given by

$\displaystyle\hskip 85.35826pt\widetilde{f}^{(k-1)}_{i,j}$
$\displaystyle:=Q_{0}^{(k-1)}\left(\bm{\mathfrak{f}_{T}}\right)^{k-1,\langle
n\rangle}_{i,j}\qquad j=-p+1,\dots,p\,.$

They rely on the approximate flux values
$\left(\mathfrak{f}_{T}\right)^{k-1,n+r}_{i,j}\approx f(w(x_{i+j},t^{n+r}))$.
In other words, we take the $(k-1)$-st discrete temporal derivative in
$x_{i+j}$ using approximate fluxes

$\displaystyle\hskip 22.76228pt\left(\mathfrak{f}_{T}\right)^{k-1,n+r}_{i,j}$
$\displaystyle:=f\left(w^{n}_{i+j}+\sum\limits_{m=1}^{k-1}\frac{(r\Delta
t)^{m}}{m!}\widetilde{w}^{(m)}_{i,j}\right),$

for $j,r=-p+1,\dots,p$. The only thing that is left to define are the
quantities $\widetilde{w}^{(m)}_{i,j}\approx\frac{\partial^{m}}{\partial
t^{m}}w(x_{i+j},t^{n})$. Their approximation makes heavy use of the Cauchy-
Kovalevskaya identity $\partial^{m}_{t}w=-\partial_{x}\partial^{m-1}_{t}f(w)$,
they are hence approximated by

$\displaystyle\hskip 85.35826pt\widetilde{w}^{(m)}_{i,j}$
$\displaystyle:=-Q_{j}^{(1)}\widetilde{\mathbf{f}}^{(m-1)}_{i,\langle
0\rangle}\qquad j=-p+1,\dots,p.$

Via the described steps, the vectors
$\widetilde{\mathbf{f}}^{(k-1)}_{i,\langle 0\rangle}$ are recursively
obtained, see also CarrilloPares2019 .

###### Definition 4

For a more precise terminology, we define the specific
$\mathtt{r}$-derivative, $q$-th order, $\mathtt{s}$-stage MDRKCAT method as
$\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$.

A summary of the $\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$ procedure to
obtain the stage values is given in Alg. 3.2. Note that the flux at the left
half-way point is obtained either from a shift of the index, i.e.
$F^{n}_{i-1/2}=F^{n}_{i-1+1/2}$ or is given by the boundary condition.

Stage solution
(${\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=2,\dots,\mathtt{s}$):

for  _ $j=-p+1$ to $p$ _ do

$(\widetilde{f}^{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1})^{(0)}_{i,j}=f(w^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1}_{i+j})$

end for

$F_{i+1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}=\sum\limits_{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}=1}^{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1}a_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}^{(1)}\Lambda^{(0)}(\widetilde{\mathbf{f}}^{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}})^{(0)}_{i,\langle
0\rangle}$

for  _ $k=2$ to $\mathtt{r}$ _ do

Get
$(\widetilde{f}^{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1})^{(k-1)}_{i,j}$
via CAT procedure.

$F_{i+1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>t}^{k-1}\sum\limits_{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}=1}^{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1}a_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}^{(k)}\Lambda^{(0)}(\widetilde{\mathbf{f}}^{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}})^{(k-1)}_{i,\langle
0\rangle}$

end for

$w^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>t}{\Delta
x}\big{(}F_{i+1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}-F_{i-1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}\big{)}$

CAT procedure CarrilloPares2019  ($k=2,\dots,\mathtt{r}$):

for  _ $j=-p+1$ to $p$ _ do

$\widetilde{w}^{(k-1)}_{i,j}=-Q_{j}^{(1)}\widetilde{\mathbf{f}}^{(k-2)}_{i,\langle
0\rangle}$

$\hskip 28.45274pt=-\frac{1}{{\Delta
x}}\sum\limits_{r=-p+1}^{p}\gamma^{1,j}_{p,r}\widetilde{f}^{(k-2)}_{i,r}$

for  _ $r=-p+1$ to $p$ _ do

$\left(\mathfrak{f}_{T}\right)^{k-1,n+r}_{i,j}=f\left(w^{n}_{i+j}+\sum\limits_{m=1}^{k-1}\frac{(r\Delta
t)^{m}}{m!}\widetilde{w}^{(m)}_{i,j}\right)$

end for

$\widetilde{f}^{(k-1)}_{i,j}=Q_{0}^{(k-1)}\left(\bm{\mathfrak{f}_{T}}\right)^{k-1,\langle
n\rangle}_{i,j}$

$\hskip 27.0301pt=\frac{1}{{\Delta
t}^{k-1}}\sum\limits_{r=-p+1}^{p}\gamma^{k-1,0}_{p,r}\left(\mathfrak{f}_{T}\right)^{k-1,n+r}_{i,j}$

end for

Algorithm 1 Stages of $\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$, an
$\mathtt{r}$-derivative, $q$-th order, $\mathtt{s}$-stage MDRKCAT method

##  4 Consistency analysis

In this section, we show that the $\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$
methods are consistent. The order of consistency is, as to be expected, the
minimum of the underlying Runge-Kutta order ($q$) and the order of the
interpolation ($2p$). Let us make the following two important assumptions:

###### Assumption 1

We assume both $f$ and $w$ to be smooth functions in $C^{\infty}$.
Furthermore, we assume that $\Delta t$ and $\Delta x$ are asymptotically
comparable in size, i.e.,

$\displaystyle\mathcal{O}(\Delta t)=\mathcal{O}(\Delta x).$

Throughout this section, we use the following notation to reduce the number of
function arguments:

${\partial_{t}^{k}f(w)}_{i+j}\equiv\partial_{t}^{k}f(w(x_{i+j},t^{n})).$

Whenever possible, similar notation is used for other functions. The time
index $n$ is only mentioned when necessary. We immediately state the main
result and thereafter, in a successive form, deduce the necessary lemmas upon
which its proof relies.

###### Theorem 4.1

The consistency order of an explicit
$\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$ method is given by $\min(2p,q)$.
Here, $q$ is the consistency order of the underlying MDRK method, while the
stencil to update $w(x_{i},t^{n})$ is given by $\\{i-p,i-p+1,\ldots,i+p\\}$.

###### Proof

The proof relies on Lemmas that will be proven in the sequel. In La. 1, it is
shown that the numerical flux difference gives the correct flux up to an order
of $2p$. We can hence substitute the exact solution $w(x,t)$ into Eq. (17b),
which immediately gives the requested result due to the fact that the Runge-
Kutta update is an integration scheme of order $q+1$:

$\displaystyle\hphantom{=}w(x_{i},t^{n+1})-w(x_{i},t^{n})+\frac{\Delta
t}{\Delta x}\big{(}F_{i+1/2}^{n}-F_{i-1/2}^{n}\big{)}$
$\displaystyle=w(x_{i},t^{n+1})-w(x_{i},t^{n})+\sum\limits_{k=1}^{\mathtt{r}}{\Delta
t}^{k}\sum\limits_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1}^{\mathtt{s}}b^{(k)}_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}\partial_{x}\partial^{k-1}_{t}f(w)_{i}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}+\mathcal{O}({\Delta
x}^{2p+1})$ $\displaystyle=\mathcal{O}({\Delta t}^{q+1})+\mathcal{O}({\Delta
x}^{2p+1})\,.$

$~{}\hfill\square$

###### Remark 6

Since the convergence order is $\min(2p,q)$, the optimal choice w.r.t.
computational efficiency is to set $p=\lceil q/2\rceil$. Hence,
“$\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$” does not contain the variable
$p$.

###### Lemma 1

The update numerical flux (18b) satisfies

$\frac{F_{i+1/2}^{n}-F_{i-1/2}^{n}}{\Delta
x}=\sum\limits_{k=1}^{\mathtt{r}}{\Delta
t}^{k-1}\sum\limits_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1}^{\mathtt{s}}b^{(k)}_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}\partial_{x}\partial^{k-1}_{t}f(w)_{i}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}+\mathcal{O}({\Delta
x}^{2p})\,.$

An analogous result holds for the stage flux (18a).

###### Proof

From La. 3 and La. 4, we obtain that for $k>1$ and
${\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1,\dots,\mathtt{s}$
there holds

$\displaystyle(\widetilde{f}^{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l})^{(k-1)}_{i,j}={\partial_{t}^{k-1}f(w)}_{i+j}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>$\displaystyle+\eta_{p}^{(k-1)}\big{(}R^{(k-1)}_{f,t}\big{)}_{i+j}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>t}^{2p-k+1}$
$\displaystyle+\xi_{p,j}^{(k-1)}\big{(}R^{(k-1)}_{f,x}\big{)}_{i}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>x}^{2p-k+1}+\mathcal{O}({\Delta x}^{2p-k+2})\,,$

with $\eta_{p}^{(k-1)}$ and $\xi_{p,j}^{(k-1)}$ real-valued coefficients; and
$R^{(k-1)}_{f,x}$, $R^{(k-1)}_{f,t}$ smooth functions of space and time. The
above formula is put into use by substituting it into the numerical flux
(18b):

$\displaystyle\frac{\Lambda^{(0)}(\widetilde{\mathbf{f}}^{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l})^{(k-1)}_{i,\langle
0\rangle}-\Lambda^{(0)}(\widetilde{\mathbf{f}}^{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l})^{(k-1)}_{i-1,\langle
0\rangle}}{\Delta x}$
$\displaystyle\quad=\frac{\left(\Uplambda^{(0)}\partial_{t}^{k-1}f(w)\right)_{i+1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}-\left(\Uplambda^{(0)}\partial_{t}^{k-1}f(w)\right)_{i-1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>x}$ $\displaystyle\quad\quad+{\Delta
t}^{2p-k+1}\cdot\eta_{p}^{(k-1)}\frac{\big{(}\Uplambda^{(0)}R^{(k-1)}_{f,t}\big{)}_{i+1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}-\big{(}\Uplambda^{(0)}R^{(k-1)}_{f,t}\big{)}_{i-1/2}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>x}$ $\displaystyle\quad\quad+{\Delta
x}^{2p-k+1}\cdot\frac{\big{(}R^{(k-1)}_{f,x}\big{)}_{i}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}-\big{(}R^{(k-1)}_{f,x}\big{)}_{i-1}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>x}\underbrace{\sum\limits_{j=-p+1}^{p}\lambda^{0}_{p,j}\xi_{p,j}^{(k-1)}}_{\mathcal{O}(1)}$
$\displaystyle\quad\quad+\frac{1}{\Delta x}\mathcal{O}({\Delta x}^{2p-k+2})$
$\displaystyle\quad\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny\eqref{eq:SplittingIntPolyCent}}}}{=}}\left(\mathcal{P}_{i}\partial_{t}^{k-1}f(w)\right)^{(1)}(x_{i},t^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}})$
$\displaystyle\quad\quad+{\Delta
t}^{2p-k+1}\cdot\eta_{p}^{(k-1)}\big{(}\mathcal{P}_{i}R^{(k-1)}_{f,t}\big{)}^{(1)}(x_{i},t^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}})+\mathcal{O}({\Delta
x}^{2p-k+1})$
$\displaystyle\quad={\partial_{x}\partial_{t}^{k-1}f(w)}_{i}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}+\mathcal{O}({\Delta
x}^{2p-k+1})\,.$

Please note that the term
$\left(\big{(}R^{(k-1)}_{f,x}\big{)}_{i}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}-\big{(}R^{(k-1)}_{f,x}\big{)}_{i-1}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}\right)/\Delta
x$ can be interpreted as a finite difference approximation of the derivative
of the smooth function $R^{k-1}_{f,x}$ and therefore remains bounded, i.e., is
$\mathcal{O}(1)$.
$~{}\hfill\square$

Before we regard the consistency analysis of the CAT steps in La. 3 and La. 4,
the follwing identity on the difference coefficients is described.

###### Lemma 2

Consider a local stencil index $j=-p+1,\dots,p$. Then there holds for
$k=1,\dots,2p-1$:

$\hskip
40.00006pt\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{k,j}(r-j)^{s}=k!\delta_{s,k}\,,\quad
s=0,\dots,2p-1\,.$

The symbol $\delta_{s,k}$ here represents the Kronecker-delta function.

###### Proof

We consider a mesh centered around $x_{i}=0$ with spatial size $\Delta x=1$.
The operator $\mathcal{Q}_{i}$ exactly interpolates the polynomial function
$\varphi(x):=(x-j)^{s}$ for $s=0,\dots,2p-1$ such that Eq. (10) becomes

$(x-j)^{s}=\mathcal{Q}_{i}\varphi(x)=\sum\limits_{r=-p+1}^{p}\ell_{p,r}(x)\varphi(r)=\sum\limits_{r=-p+1}^{p}\ell_{p,r}(x)(r-j)^{s}\,.$

Deriving the above relation $k$ times in $x$ and thereafter evaluating in
$x=j$ gives the result. $\newline ~{}\hfill\square$

Now we can provide a consistency proof for the CAT procedure in Alg. 3.2. To
this purpose we establish the following notation for the exact Taylor
approximation in time of order $k\in\mathbb{N}$,

$T^{k}_{i,j}(\tau):=\sum_{m=0}^{k}\frac{\tau^{m}}{m!}\partial_{t}^{m}w(x_{i+j},t^{n})\,.$
(19)

It is the $k$-th order approximation of $w(x_{i+j},t^{n}+\tau)$. The proof
itself is built in a similar fashion as (CarrilloPares2019, , Theorem 2) and
(ZorioEtAl, , Proposition 1).

###### Lemma 3

For $k=1$ the steps of the CAT algorithm (Alg. 3.2, right side) satisfy

$\displaystyle\widetilde{w}^{(1)}_{i,j}$
$\displaystyle={\partial_{t}w}_{i+j}+\xi_{p,j}^{(1)}\big{(}R^{(1)}_{w}\big{)}_{i}\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p})\,,$
$\displaystyle\widetilde{f}^{(1)}_{i,j}$
$\displaystyle={\partial_{t}f(w)}_{i+j}+\eta_{p}^{(1)}\big{(}R^{(1)}_{f,t}\big{)}_{i+j}\cdot{\Delta
t}^{2p-1}+\xi_{p,j}^{(1)}\big{(}R^{(1)}_{f,x}\big{)}_{i}\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p})\,,$

with real-valued coefficients $\xi_{p,j}^{(1)}$, $\eta_{p}^{(1)}$; and smooth
functions $\big{(}R^{(1)}_{w}\big{)}_{i}$,
$\big{(}R^{(1)}_{f,t}\big{)}_{i+j}$, $\big{(}R^{(1)}_{f,x}\big{)}_{i}$.
(Please note again that $(\cdot)_{i}$ stands for function-evaluation at
$x=x_{i}$.)

###### Proof

A straightforward computation on $\widetilde{w}^{(1)}_{i,j}$, using the Taylor
expansion of $f(w)$, reveals that

$\displaystyle\widetilde{w}^{(1)}_{i,j}$
$\displaystyle=-Q_{j}^{(1)}\widetilde{\mathbf{f}}^{(0)}_{i,\langle
0\rangle}=-\frac{1}{{\Delta
x}}\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,j}\widetilde{f}^{(0)}_{i,r}$
$\displaystyle=-\frac{1}{{\Delta
x}}\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,j}f(w(x_{i+r},t^{n}))$
$\displaystyle=-\frac{1}{{\Delta
x}}\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,j}\left(\sum\limits_{s=0}^{2p}\frac{\left((r-j)\Delta
x\right)^{s}}{s!}{\partial_{x}^{s}f(w)}_{i+j}+\mathcal{O}({\Delta
x}^{2p+1})\right)$ $\displaystyle=-\frac{1}{{\Delta
x}}\sum\limits_{s=0}^{2p-1}\frac{{\Delta
x}^{s}}{s!}{\partial_{x}^{s}f(w)}_{i+j}\left(\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,j}(r-j)^{s}\right)$
$\displaystyle\qquad-\left[\frac{1}{(2p)!}\left(\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,j}(r-j)^{2p}\right){\partial_{x}^{2p}f(w)}_{i}+\mathcal{O}({\Delta
x})\right]{\Delta x}^{2p-1}$ $\displaystyle\qquad+\mathcal{O}({\Delta
x}^{2p})$ $\displaystyle\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny La.
\ref{lem:LagrSumLemma}}}}{=}}-{\partial_{x}f(w)}_{i+j}+\xi_{p,j}^{(1)}\big{(}R^{(1)}_{w}\big{)}_{i}\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p})$
$\displaystyle={\partial_{t}w}_{i+j}+\xi_{p,j}^{(1)}\big{(}R^{(1)}_{w}\big{)}_{i}\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p})\,,$

in which we define $\xi_{p,j}^{(1)}$ and $\big{(}R^{(1)}_{w}\big{)}_{i}$ by

$\displaystyle\xi_{p,j}^{(1)}:=\frac{-1}{(2p)!}\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,j}(r-j)^{2p}\,,\quad\big{(}R^{(1)}_{w}\big{)}_{i}:={\partial_{x}^{2p}f(w)}_{i}.$

The succeeding step of the method evaluates the flux $f$ in an approximate
Taylor series. We find,

$\displaystyle\left(\mathfrak{f}_{T}\right)^{1,n+r}_{i,j}$
$\displaystyle=f\left(w(x_{i+j},t^{n})+(r\Delta
t)\widetilde{w}^{(1)}_{i,j}\right)$
$\displaystyle\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny\eqref{eq:exactTaylorTime}}}}{=}}f\left(T^{1}_{i,j}(r{\Delta
t})+(r\Delta t)\xi_{p,j}^{(1)}\big{(}R^{(1)}_{w}\big{)}_{i}\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p+1})\right)$ $\displaystyle=(f\circ
T^{1}_{i,j})(r{\Delta t})$ $\displaystyle\qquad+(f^{\prime}\circ
T^{1}_{i,j})(r{\Delta t})\cdot\left[(r\Delta
t)\xi_{p,j}^{(1)}\big{(}R^{(1)}_{w}\big{)}_{i}\right]\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p+1})$ $\displaystyle=(f\circ
T^{1}_{i,j})(r{\Delta t})$ $\displaystyle\qquad+(f^{\prime}\circ
T^{1}_{i,j})(0)\cdot\left[(r\Delta
t)\xi_{p,j}^{(1)}\big{(}R^{(1)}_{w}\big{)}_{i}\right]\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p+1})$ $\displaystyle=(f\circ
T^{1}_{i,j})(r{\Delta t})+{f^{\prime}(w)}_{i}\left[(r\Delta
t)\xi_{p,j}^{(1)}\big{(}R^{(1)}_{w}\big{)}_{i}\right]\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p+1})\,.$

Consequently, we can find for the temporal interpolation of the fluxes:

$\displaystyle\widetilde{f}^{(1)}_{i,j}$
$\displaystyle=Q_{0}^{(1)}\left(\bm{\mathfrak{f}_{T}}\right)^{1,\langle
n\rangle}_{i,j}=\frac{1}{{\Delta
t}}\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,0}\left(\mathfrak{f}_{T}\right)^{1,n+r}_{i,j}$
$\displaystyle=\frac{1}{{\Delta
t}}\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,0}\left[\sum\limits_{s=0}^{2p}\left(\frac{(r{\Delta
t})^{s}}{s!}\frac{\mathrm{d}^{s}(f\circ
T^{1}_{i,j})}{{\mathrm{d}\tau}^{s}}(0)\right)+\mathcal{O}({\Delta
x}^{2p+1})\right]$
$\displaystyle\qquad+{f^{\prime}(w)}_{i}\underbrace{\left(\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,0}r\right)}_{=1,\text{~{}La.~{}}\ref{lem:LagrSumLemma}}\xi_{p,j}^{(1)}\big{(}R^{(1)}_{w}\big{)}_{i}\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p})$ $\displaystyle=\frac{1}{{\Delta
t}}\sum\limits_{s=0}^{2p-1}\frac{{\Delta
t}^{s}}{s!}\underbrace{\left(\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,0}r^{s}\right)}_{=\delta_{1,s},\text{~{}La.~{}}\ref{lem:LagrSumLemma}}\frac{\mathrm{d}^{s}(f\circ
T^{1}_{i,j})}{{\mathrm{d}\tau}^{s}}(0)$
$\displaystyle\qquad+\underbrace{\frac{1}{(2p)!}\left(\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,0}r^{2p}\right)}_{=:\eta_{p}^{(1)}}\underbrace{\frac{\mathrm{d}^{2p}(f\circ
T^{1}_{i,j})}{{\mathrm{d}\tau}^{2p}}(0)}_{=:\big{(}R^{(1)}_{f,t}\big{)}_{i+j}}\cdot{\Delta
t}^{2p-1}$
$\displaystyle\qquad+\xi_{p,j}^{(1)}\underbrace{{f^{\prime}(w)}_{i}\big{(}R^{(1)}_{w}\big{)}_{i}}_{=:\big{(}R^{(1)}_{f,x}\big{)}_{i}}\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p})$ $\displaystyle=\frac{\mathrm{d}(f\circ
T^{1}_{i,j})}{{\mathrm{d}\tau}}(0)+\eta_{p}^{(1)}\big{(}R^{(1)}_{f,t}\big{)}_{i+j}\cdot{\Delta
t}^{2p-1}+\xi_{p,j}^{(1)}\big{(}R^{(1)}_{f,x}\big{)}_{i}\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p})$
$\displaystyle={\partial_{t}f(w)}_{i+j}+\eta_{p}^{(1)}\big{(}R^{(1)}_{f,t}\big{)}_{i+j}\cdot{\Delta
t}^{2p-1}+\xi_{p,j}^{(1)}\big{(}R^{(1)}_{f,x}\big{)}_{i}\cdot{\Delta
x}^{2p-1}+\mathcal{O}({\Delta x}^{2p}).$

$~{}\hfill\square$

Via induction one can generalize this result.

###### Lemma 4

For $k=2,\dots 2p-1$ the steps of the CAT algorithm (Alg. 3.2, right side)
satisfy

$\displaystyle\widetilde{w}^{(k)}_{i,j}$
$\displaystyle={\partial_{t}^{k}w}_{i+j}+\xi_{p,j}^{(k)}\big{(}R^{(k)}_{w}\big{)}_{i}\cdot{\Delta
x}^{2p-k}+\mathcal{O}({\Delta x}^{2p-k+1})\,,$
$\displaystyle\widetilde{f}^{(k)}_{i,j}$
$\displaystyle={\partial_{t}^{k}f(w)}_{i+j}+\eta_{p}^{(k)}\big{(}R^{(k)}_{f,t}\big{)}_{i+j}\cdot{\Delta
t}^{2p-k}$ $\displaystyle\hskip
58.04343pt+\xi_{p,j}^{(k)}\big{(}R^{(k)}_{f,x}\big{)}_{i}\cdot{\Delta
x}^{2p-k}+\mathcal{O}({\Delta x}^{2p-k+1})\,,$

with

$\displaystyle\eta_{p}^{(k)}$
$\displaystyle:=\frac{1}{(2p)!}\left(\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{k,0}r^{2p}\right),\qquad\xi_{p,j}^{(k)}:=-\sum\limits_{r=-p+1}^{p}\gamma_{p,r}^{1,j}\xi_{p,r}^{(k-1)}$
$\displaystyle\big{(}R^{(k)}_{w}\big{)}_{i}$
$\displaystyle:=\big{(}{f^{\prime}(w)}_{i}\big{)}^{k-1}{\partial_{x}^{2p}f(w)}_{i},\qquad\big{(}R^{(k)}_{f,t}\big{)}_{i+j}:=\frac{\mathrm{d}^{2p}(f\circ
T^{k}_{i,j})}{{\mathrm{d}\tau}^{2p}}(0),$
$\displaystyle\big{(}R^{(k)}_{f,x}\big{)}_{i}$
$\displaystyle:={f^{\prime}(w)}_{i}\big{(}R^{(k)}_{w}\big{)}_{i}\,.$

##  5 von Neumann stability of $\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$
methods

The original CAT procedure was developed with the intention to create a scheme
which linearly reduces back to high-order Lax-Wendroff methods. As a
consequence, the scheme is CFL-1 stable for linear equations
(CarrilloPares2019, , Theorem 1). In this section, we discuss the stability
properties of the $\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$ scheme
presented in this work.

######  Theorem 5.1 (MDRK-LW scheme)

Explicit $\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$ methods for the linear
advection flux $f(w)=\alpha w$ reduce to the numerical scheme

$\displaystyle
w^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>$\displaystyle=w^{n}_{j}+\sum\limits_{k=1}^{\mathtt{r}}(-1)^{k}\alpha^{k}{\Delta
t}^{k}\left(\sum\limits_{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}=1}^{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1}a_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}^{(k)}P^{(k)}\mathbf{w}^{n,{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color<EMAIL_ADDRESS>j\rangle}\right)\qquad{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1,\dots,\mathtt{s}\,,$
(20a) $\displaystyle w^{n+1}_{j}$
$\displaystyle=w^{n}_{j}+\sum\limits_{k=1}^{\mathtt{r}}(-1)^{k}\alpha^{k}{\Delta
t}^{k}\left(\sum\limits_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1}^{\mathtt{s}}b^{(k)}_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}P^{(k)}\mathbf{w}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>j\rangle}\right)\,,$ (20b)

with
$P^{(k)}\mathbf{w}^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>j\rangle}$ the centered difference approximation of the $k$-th spatial
derivative.

The proof is similar to (CarrilloPares2019, , Theorem 1) and is hence left
out. In this form it is possible to perform a von Neumann stability analysis
(see for example QuarteroniNumA ; 2004_Strikwerda ). That is, we fill in the
Fourier mode $w^{n}_{j}=\mathcal{W}(t^{n})e^{\mathrm{i}\kappa x_{j}}$ with
wave number $\kappa\in\mathbb{Z}$ and search for the the amplification factors
via the relations

$\mathcal{W}(t^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}):=g^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}(\kappa)\mathcal{W}(t^{n})\quad\text{and}\quad\mathcal{W}(t^{n+1}):=g(\kappa)\mathcal{W}(t^{n})\,.$

Doing so gives an additional recurrence relation.

###### Proposition 1

The amplification factors obtained from a von Neumann analysis on the MDRK-LW
scheme Eq. (20) is defined by the recurrence relations

$\displaystyle
g^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}(\kappa)$
$\displaystyle=1+\sum\limits_{k=1}^{\mathtt{r}}(-1)^{k}\sigma^{k}~{}\mathtt{P}^{(k)}(\kappa)\left(\sum\limits_{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}=1}^{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1}a_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}^{(k)}g^{n,{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}(\kappa)\right)\qquad{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1,\dots,\mathtt{s}\,,$
(21a) $\displaystyle g(\kappa)$
$\displaystyle=1+\sum\limits_{k=1}^{\mathtt{r}}(-1)^{k}\sigma^{k}~{}\mathtt{P}^{(k)}(\kappa)\left(\sum\limits_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1}^{\mathtt{s}}b^{(k)}_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}g^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}(\kappa)\right)\,,$
(21b)

with wave number $\kappa\in\mathbb{Z}$, $\sigma:=\frac{\alpha{\Delta
t}}{{\Delta x}}$ the corresponding CFL number and

$\mathtt{P}^{(k)}(\kappa):=\sum_{r=-p}^{p}\delta^{k}_{p,r}e^{\mathrm{i}r\kappa{\Delta
x}}\,.$ (22)

The term $\mathtt{P}^{(k)}(\kappa)$ can be interpreted as the $k$-th
derivative of the centered $(2p+1)$-point Lagrangian interpolation Eq. (11a)
using Fourier basis
$\left\\{e^{\mathrm{i}l\chi}~{}|~{}l\in\mathbb{Z}\right\\}$ with the grid
frequency $\chi=\kappa\Delta x$.

In Tbl. 2 we display the amplification factors $g(\kappa)$ for the considered
MDRK schemes summarized earlier in Fig. 1. Of main interest w.r.t. to these
amplification factors is to obtain a critical CFL value
$\sigma^{*}\in\mathbb{R}^{+}$ such that

$|g(\kappa)|\leq 1\quad\Leftrightarrow\quad\frac{\alpha{\Delta t}}{{\Delta
x}}\leq\sigma^{*}\,.$

We have used the Symbolic Math Toolbox in MATLAB 2020_matlabsymbolic  along
with a bisection method on the CFL variable $\sigma$ in Eq. (21b) to
numerically obtain a $\sigma^{*}$:

  * •

In our approach we take $\Delta x=1$; only the frequency of $g(\kappa)$ is
influenced by $\Delta x$, there is no change in absolute value. Assume we
would compare mesh sizes $\Delta x$ and $\Delta\tilde{x}$, then

$\mathtt{P}^{(k)}(\kappa)=\sum\limits_{r=-p}^{p}\delta^{k}_{p,r}e^{\mathrm{i}r\kappa{\Delta
x}}=\sum\limits_{r=-p}^{p}\delta^{k}_{p,r}e^{\mathrm{i}r\left(\kappa\frac{\Delta
x}{\Delta\tilde{x}}\right)\Delta\chi}=\mathtt{P}^{(k)}(\widetilde{\kappa})\,,$

where $\widetilde{\kappa}:=\kappa\frac{\Delta x}{\Delta\tilde{x}}$. Thus
behavior of $g(\kappa)$ is the same up to a recalibration of the frequency
space.

  * •

Via the toolbox, the value $\max|g(\kappa)|$ is calculated on a uniform
1000-cell mesh of $\kappa\in[-\pi,\pi]$. This domain suffices for this
purpose, since for $\Delta x=1$ all terms $\mathtt{P}^{(k)}(\kappa)$ in Eq.
(22) are $2\pi$-periodic.

Table 2: Amplification factors $g(\kappa)$ of the considered MDRK schemes (see
Appendix A). The CFL criteria (up to four decimals) obtained from
$|g(\kappa)|\leq 1$ are also shown.

Scheme | $g(\kappa)$ | CFL
---|---|---
$2\text{DRKCAT}3\text{-}2$ | $1-\sigma~{}\mathtt{P}^{(1)}+\frac{1}{2}\sigma^{2}\left(\frac{2}{3}\mathtt{P}^{(1)}\mathtt{P}^{(1)}+\frac{1}{3}\mathtt{P}^{(2)}\right)-\frac{1}{6}\sigma^{3}~{}\mathtt{P}^{(1)}\mathtt{P}^{(2)}$ | 1.2954
$2\text{DRKCAT}4\text{-}2$ | $1-\sigma~{}\mathtt{P}^{(1)}+\frac{1}{2}\sigma^{2}~{}\mathtt{P}^{(2)}-\frac{1}{6}\sigma^{3}~{}\mathtt{P}^{(1)}\mathtt{P}^{(2)}+\frac{1}{24}\sigma^{4}~{}\mathtt{P}^{(2)}\mathtt{P}^{(2)}$ | 1.4718
$2\text{DRKCAT}5\text{-}3$ | $1-\sigma~{}\mathtt{P}^{(1)}+\frac{1}{2}\sigma^{2}~{}\mathtt{P}^{(2)}-\frac{1}{6}\sigma^{3}~{}\mathtt{P}^{(1)}\mathtt{P}^{(2)}+\frac{1}{24}\sigma^{4}~{}\mathtt{P}^{(2)}\mathtt{P}^{(2)}$ | 1.0619
| $-\frac{1}{120}\sigma^{5}~{}\mathtt{P}^{(1)}\mathtt{P}^{(2)}\mathtt{P}^{(2)}+\frac{1}{600}\sigma^{6}~{}\mathtt{P}^{(2)}\mathtt{P}^{(2)}\mathtt{P}^{(2)}$ |
$3\text{DRKCAT}5\text{-}2$ | $1-\sigma~{}\mathtt{P}^{(1)}+\frac{1}{2}\sigma^{2}~{}\mathtt{P}^{(2)}-\frac{1}{6}\sigma^{3}~{}\mathtt{P}^{(3)}+\frac{1}{24}\sigma^{4}~{}\mathtt{P}^{(1)}\mathtt{P}^{(3)}$ | 0.4275
| $-\frac{1}{120}\sigma^{5}~{}\mathtt{P}^{(2)}\mathtt{P}^{(3)}+\frac{1}{900}\sigma^{6}~{}\mathtt{P}^{(3)}\mathtt{P}^{(3)}$ |
$3\text{DRKCAT}7\text{-}3$ | $1-\sigma~{}\mathtt{P}^{(1)}+\frac{1}{2}\sigma^{2}~{}\mathtt{P}^{(2)}-\frac{1}{6}\sigma^{3}~{}\mathtt{P}^{(3)}+\frac{1}{24}\sigma^{4}~{}\mathtt{P}^{(1)}\mathtt{P}^{(3)}$ | 0.2300
| $-\frac{1}{120}\sigma^{5}~{}\mathtt{P}^{(2)}\mathtt{P}^{(3)}+\frac{1}{720}\sigma^{6}~{}\mathtt{P}^{(3)}\mathtt{P}^{(3)}-\frac{1}{5040}\sigma^{7}~{}\mathtt{P}^{(1)}\mathtt{P}^{(2)}\mathtt{P}^{(3)}$ |
| $+\left(\frac{1}{23520}-\frac{\sqrt{2}}{70560}\right)\sigma^{8}~{}\mathtt{P}^{(2)}\mathtt{P}^{(3)}\mathtt{P}^{(3)}-\left(\frac{11}{1481760}-\frac{\sqrt{2}}{246960}\right)\sigma^{9}~{}\mathtt{P}^{(3)}\mathtt{P}^{(3)}\mathtt{P}^{(3)}$ |
$4\text{DRKCAT}6\text{-}2$ | $1-\sigma~{}\mathtt{P}^{(1)}+\frac{1}{2}\sigma^{2}~{}\mathtt{P}^{(2)}-\frac{1}{6}\sigma^{3}~{}\mathtt{P}^{(3)}+\frac{1}{24}\sigma^{4}~{}\mathtt{P}^{(4)}$ | 0.8563
| $-\frac{1}{120}\sigma^{5}~{}\mathtt{P}^{(1)}\mathtt{P}^{(4)}+\frac{1}{720}\sigma^{6}~{}\mathtt{P}^{(2)}\mathtt{P}^{(4)}-\frac{1}{6480}\sigma^{7}~{}\mathtt{P}^{(3)}\mathtt{P}^{(4)}+\frac{1}{77760}\sigma^{8}~{}\mathtt{P}^{(4)}\mathtt{P}^{(4)}$ |

The critical CFL values $\sigma^{*}$ up to four decimals are shown in Tbl. 2.
Notice that the two-derivative schemes improve the linear stability, whereas
the other schemes reduce the stability compared to the original CAT method. To
put this observation into perspective, let us point out that the CAT algorithm
uses $2p$ derivatives, and thus is based on high-order Lax-Wendroff methods
with $\mathtt{r}$ even. If, only for the sake of discussion, we take _uneven_
$\mathtt{r}$, we get another picture. For $\mathtt{r}=1$, we obtain the
forward-time central-space scheme, which is infamous for being unconditionally
unstable LEV  and for odd $\mathtt{r}>1$, we obtain CFL numbers smaller than
one. Hence, we can make some important observations:

  * •

Most importantly we can conclude that a method of lines (MOL) viewpoint is
inadequate. Solely regarding the stability regions $\mathcal{R}$ in Fig. 1
would give the idea that the $3\text{DRKCAT}7\text{-}3$ scheme provides the
best stability. This is clearly not the case since the even-derivative schemes
are shown to be better in terms of stability.

  * •

Choosing even order derivatives gives better results than choosing an odd
number of derivatives, i.e. the two- and four-derivative schemes show better
stability than the three-derivative schemes. This is in very good agreement
with the observations on the original CAT method.

  * •

In contrast to Taylor-methods, where only the number of derivatives can be
prescribed, MDRK methods have two free parameters, namely the number of
derivatives _and_ the number of stages. We observe that stages and derivatives
highly influence the stability properties of the MDRKCAT method. This allows
the identification of well-suited MDRK schemes and gives more flexibility
compared to original Taylor-methods.

##  6 Numerical results

In this section, we show numerical results validating our analytical findings.
By means of several continuous test cases ranging from scalar PDEs to the
system of Euler equations, we show that the expected orders of convergence are
obtained. For brevity, we do not include flux limiting techniques to this
work; hence, we avoid setups where shock formation occurs. Note that flux
limiting can be incorporated in a straightforward way as in CarrilloPares2019
.

The measure for the accuracy in this section is the scaled $l_{1}$-error at
time $t^{N}\equiv T_{\text{end}}$

$\|\mathbf{w}(T_{\text{end}})-\mathbf{w}^{N}\|_{1}:=\Delta
x\sum\limits_{i=1}^{M}|w(x_{i},T_{\text{end}})-w_{i}^{N}|\,,$

$\mathbf{w}:t\mapsto\mathbf{w}(t)$ being a function of time returning a vector
of exact solution values (or a reference solution) in the nodes
$x_{1},\dots,x_{M}$, and $\mathbf{w}^{N}$ being the vector of approximations
$w^{N}_{i}$ at $t^{N}$. For systems, the sum of the $l_{1}$-errors
corresponding to the separate solution variables is considered. All displayed
convergence plots begin at $M=8$ cells and double the amount of cells with
each iteration. In order to enforce the adopted CFL value $\sigma$, the local
eigenvalues $\lambda^{n}_{\text{eig},i}$ of the Jacobian w.r.t. $w_{i}^{n}$
are computed. By means of the relation

$\Delta t^{n}:=\frac{\sigma\Delta x}{\max_{i}|\lambda^{n}_{\text{eig},i}|}\,$

the timestep is then computed. As the maximum eigenvalue in the computational
domain varies over time, a non-constant timestep is prescribed. This
highlights the ability of the novel method to use varying timestep sizes, see
Rem. 1.

###  6.1 Burgers equation

First, we consider Burgers equation

$\partial_{t}w+\partial_{x}\left(\frac{w^{2}}{2}\right)=0,$

with the cosine-wave initial condition with periodic boundary conditions

$w(x,0)=\frac{1}{4}\cos(\pi x)\qquad\text{on}\quad x\in[0,2]\,,$ (23)

to certify the accuracy $\min(2p,q)$ obtained in Thm. 4.1. Since the cosine-
wave Eq. (23) has both positive and negative values, the characteristic lines
must cross and a shock is formed at some point. The breaking time of the wave
LEV ; 2011_Whitham  is at $t^{*}=\frac{4}{\pi}$. Hence we set the final time
to $T_{\text{end}}=0.8$, well before shock formation.

Figure 2: Convergence order of the CAT and explicit MDRKCAT methods applied to
Burgers equation on the cosine wave $w_{0}(x)=\frac{1}{4}\cos(\pi x)$ up to
$T_{\text{end}}=0.8$ with CFL $\sigma=0.5$. The $2\text{DRKCAT}5\text{-}3$,
$3\text{DRKCAT}5\text{-}2$ and $4\text{DRKCAT}6\text{-}2$ schemes behave in a
very alike manner. For large $\Delta t$ and large $p$, we have observed
stability problems for the CAT method. Those divergent results of the CAT
method have been omitted in the convergence plot.

Using the characteristic lines solution, $l_{1}$-errors have been calculated
for the CAT-$p$ methods ($p=1,\dots,6$) and the MDRKCAT methods (Appendix A).
The results with CFL $\sigma=0.5$ are visualized in Fig. 2. All expected
convergence orders are obtained; order $2p$ for the CAT-$p$ methods and at
least order $q$ of the MDRK schemes. The schemes $2\text{DRKCAT}5\text{-}3$,
$3\text{DRKCAT}5\text{-}2$ and $3\text{DRKCAT}7\text{-}3$ behave better than
expected. This is caused by the fact that the spatial order of accuracy is
higher than the temporal one; and spatial errors dominate the overall behavior
at least for ‘large’ $\Delta t$. We have observed similar behavior also for
other schemes where $2p>q$. Given that CAT methods have been designed as a
natural generalization of Lax-Wendroff methods with an even-order accuracy,
odd-order MDRK schemes take advantage here. We expect this behavior to become
more apparent when computing with finer machine precision. Convergence plots
such as in Fig. 2 will then manifest as a stretched out S-curve.

Further, we notice that for higher values of $p$ and larger values $\Delta t$,
the CAT methods tend to be less stable. Unstable results have been left out in
Fig. 2 for $p=4,5$ and $6$. Even though CAT methods have been shown to be
linearly stable under a CFL-$1$ condition CarrilloPares2019 , rapid divergence
is observed well before shocks are formed for larger values of $p$ at regions
where the derivative of the solution $w^{\prime}(x)$ is large in absolute
value. The MDRKCAT methods used in this work do not suffer this fate.

Next, we perform a similar study having the initial solution,

$w(x,0)=\frac{1}{4}e^{\cos(\pi x)+\sin(\pi x)}\qquad\text{on}\quad
x\in[0,2]\,,$

with periodic boundary conditions. Having a steeper peak than the cosine (23),
the breaking time will be earlier. We find $t^{*}=\frac{4}{\pi e}$, and choose
$T_{\text{end}}=0.3$ accordingly.

Figure 3: Burgers equation applied to $w_{0}(x)=\frac{1}{4}e^{\cos(\pi
x)+\sin(\pi x)}$ up to $T_{\text{end}}=0.3$ with CFL $\sigma=0.5$. Left:
$3\text{DRKCAT}7\text{-}3$ solution using $M=8192$ cells, less nodal points
are shown for better visual distinctness. Right: convergence order of the
explicit MDRKCAT methods.

In Fig. 3 the final solution at $T_{\text{end}}=0.3$ is visualized and
accuracy is studied in the convergence plots, solely focused on the MDRKCAT
method. The behavior is very similar to the previous case, except that the
$3\text{DRKCAT}7\text{-}3$ scheme is driven back faster to order $7$.

###  6.2 Buckley-Leverett equation

Next, we consider the Buckley-Leverett flux LEV ,

$f(w)=\frac{4w^{2}}{4w^{2}+(1-w)^{2}}\,.$

This flux is non-convex and introduces more nonlinearities compared to the
Burgers flux. We consider the initial condition

$w(x,0)=1-\frac{3}{4}\cos^{2}\left(\frac{\pi}{2}x\right)\qquad\text{on}\quad
x\in[-1,1]\,.$

The typical Buckley-Leverett profile consist of a shock wave followed directly
by a rarefaction wave. We set $T_{\text{end}}=0.1$ to remain continuous and be
able to calculate the exact solution via its characteristics. The solution and
convergence plots are visualized in Fig. 4 with CFL $\sigma=0.5$. We notice
that the numerical solution tends toward the expected Buckley-Leverett
profile. All schemes converge with the expected accuracy in a similar way as
for Burgers equation.

Figure 4: Buckley-Leverett equation applied to
$w_{0}(x)=1-\frac{3}{4}\cos^{2}\left(\frac{\pi}{2}x\right)$ up to
$T_{\text{end}}=0.1$ with CFL $\sigma=0.5$. Left: $3\text{DRKCAT}7\text{-}3$
solution using $M=8192$ cells, less nodal points are shown for better visual
distinctness. Right: convergence order of the explicit MDRKCAT methods.

###  6.3 One-dimensional Euler equations

Finally, we consider the Euler equations of gas dynamics

$\partial_{t}{w}+\partial_{x}{f}({w})=0\,,$

in which

${w}=\begin{pmatrix}\rho\\\ \rho u\\\
E\end{pmatrix}\,,\quad{f}({w})=\begin{pmatrix}\rho u\\\ \rho u^{2}+p\\\
u(E+p)\end{pmatrix}\,,$

with $\rho$ the density, $u$ the velocity, $E$ the energy of the system and
$p$ the pressure. The system is closed via the equation of state for an ideal
gas

$p=(\gamma-1)\left(E-\frac{1}{2}\rho u^{2}\right)\,,$

with $\gamma$ being the ratio of specific heats, assumed to be $1.4$ LEV ;
2011_Whitham .

First we initialize the primitive variables $(\rho,u,p)$ such that the Euler
equations describe the linear advection of a density profile. To this end, we
take

$\rho(x,0)=1+0.3\sin(\pi x)\qquad\text{on}\quad x\in[0,4]\,,$

and set both $u_{0}(x,0)$ and $p_{0}(x,0)$ to be one. Periodic boundary
conditions are used and $T_{\text{end}}$ is set to $0.8$. In Fig. 5 the
corresponding convergence plots are displayed with CFL $\sigma=0.5$.
Immediately starting from the coarsest meshes the expected convergence orders
are obtained.

Figure 5: Convergence order of the explicit MDRKCAT methods applied to
$\rho_{0}(x)=1+0.3\sin(\pi x)$, $u_{0}(x)=p_{0}(x)=1$ on $x\in[0,4]$ up to
$T_{\text{end}}=0.8$ with $\sigma=0.5$.

Secondly, we consider the initial condition

${w}(x,0)=\frac{1}{4}\begin{pmatrix}3\\\ 1\\\ 3\end{pmatrix}+\frac{\sin(\pi
x)}{2}\begin{pmatrix}1\\\ 1\\\ 1\end{pmatrix}\qquad\text{on}\quad
x\in[0,2]\,,$

with periodic boundary conditions and $T_{\text{end}}=0.2$ for ${w}$ to remain
continuous. In Fig. 6 the solution is visualized for the
$3\text{DRKCAT}7\text{-}3$ scheme with $p=4$, CFL $\sigma=0.5$ on $M=4096$
cells.

In order to inspect accuracy, a reference solution has been computed via a
discontinuous Galerkin (DG) method. We have used third-order polynomials in
space, and a third-order strong-stability-preserving Runge-Kutta method in
time Gottlieb2001 . The reference computation was executed on $10240$ cells
with a CFL number of $\sigma=0.15$.

Convergence plots have been generated in Fig. 7 using CFL $\sigma=0.15$. All
expected orders were obtained. For smaller $\Delta x$, the $l_{1}$-error
converges towards approximately $2\cdot 10^{-11}$, which is the accuracy of
the reference DG solution.

Akin to the earlier cases, a CFL value of $\sigma=0.5$ was attempted for the
construction of the convergence plots. However not all simulations were
stable, more specifically the $3\text{DRKCAT}7\text{-}3$ scheme using $p=4$
diverged for $M=8,16,32$ and $64$ cells.

Figure 6: MDRKCAT 1D Euler solution of $\rho_{0}(x)=0.75+0.5\sin(\pi x)$,
$(\rho u)_{0}(x)=0.25+0.5\sin(\pi x)$ and $E_{0}(x)=0.75+0.5\sin(\pi x)$ up to
$T_{\text{end}}=0.2$ with $\sigma=0.5$. The $3\text{DRKCAT}7\text{-}3$ scheme
has been used on $M=4096$ cells; less nodal points are shown for better visual
distinctness. Figure 7: Convergence order of the explicit MDRKCAT methods
applied to $\rho_{0}(x)=0.75+0.5\sin(\pi x)$, $(\rho
u)_{0}(x)=0.25+0.5\sin(\pi x)$ and $E_{0}(x)=0.75+0.5\sin(\pi x)$ up to
$T_{\text{end}}=0.2$ with $\sigma=0.15$. Very similar behavior can be seen
between the $2\text{DRKCAT}3\text{-}2$ and the $2\text{DRKCAT}4\text{-}2$
schemes. The same can be said for the 5th order MDRKCAT schemes and
$4\text{DRKCAT}6\text{-}2$. For comparison, $l_{1}$-errors of a DG code with
basis functions in $\mathbb{P}_{3}$ and $\mathbb{P}_{4}$ and CFL $\sigma=0.1$
are visualized.

In order to better grasp the efficiency of the MDRKCAT methods, in the same
Fig. 7 convergence plots have been generated by means of a DG code that uses
polynomial basis functions in $\mathbb{P}_{3}$ and $\mathbb{P}_{4}$ for each
cell respectively. A fourth order SSP-RK scheme (2002_Spiteri_Ruuth, , order
4, p.21) has been used as time integrator. The same amount of cells
$M=8,\dots,1024$ has been used as for the MDRKCAT runs, the CFL $\sigma=0.1$
with a maximum eigenvalue of $1.5$ so that $\Delta t=\Delta
x\frac{\sigma}{1.5}$. The $l_{1}$-errors, computed at cell-midpoints, have
been generated relative to the earlier mentioned order three SSP-DG reference
solution.

Overall, the MDRKCAT methods compare well with the DG solutions. A large
discrepancy can be noticed in the manner at which the expected convergence
order is achieved; the MDRKCAT methods gradually head toward order
$\min(2p,q)$ with each refinement, whereas the DG solvers achieve convergence
already going from 32 to 64 cells. This is to be expected: By definition the
MDRKCAT methods only approximate the time derivatives of the flux
$\partial^{k-1}_{t}{f}({w})$. Hence the achieved accuracy is intertwined with
the mesh resolution of the problem at hand. For a lower amount of cells $M$
the numerical flux $F^{n}_{i+1/2}$ at the faces can thus not be an accurate
representation, whereas the DG solvers calculate the fluxes at the half-way
points $i+1/2$ on the basis of the exact flux ${f}({w})$. As soon as enough
cells $M$ are used to finely represent the initial data, full advantage can be
taken of the CAT method.

Moreover, the difference between the methods should be brought into
perspective by studying the amount of effective spatial degrees of freedom
(DOF) and the effective spatial size that influences the order of accuracy. DG
methods make use of numerical integration points on each cell for the
integration of the solution variable multiplied with the chosen basis
functions Cockburn2000 . This illustrates why the DG solvers more quickly
capture the expected convergence order and why a direct comparison of the DG
schemes and the MDRKCAT is difficult in Fig. 7. The actual amount of spatial
DOF used by $\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$ schemes is $\lceil
q/2\rceil M$; each node uses it’s own local stencil in the calculations.
However, as explained in CarrilloPares2019 , the local stencils are merely a
manner to assure that the CAT methods linearly reduce back to Lax-Wendroff
schemes. The same accuracy is achieved by the approximate Taylor methods in
ZorioEtAl  of which the CAT procedure is established. Summing up, we can
conclude that the novel $\mathtt{r}\text{DRKCAT}q\text{-}\mathtt{s}$ schemes
compare well with a state-of-the-art DG solver in terms of accuracy.

##  7 Conclusion and outlook

In this paper we have formulated a family of Jacobian-free multistage
multiderivative solvers for hyperbolic conservation laws, so-called MDRKCAT
methods. Following the Compact Approximate Taylor (CAT) method in
CarrilloPares2019 , instead of computing the exact flux derivative
expressions, local approximations for the time derivatives of the fluxes are
obtained recursively. There are many advantages by virtue of this procedure:
no costly symbolic computations are needed; and we hope that many
multiderivative Runge-Kutta (MDRK) schemes can now actually be of practical
use.

Both theoretically and numerically it is proven that the desired convergence
order $\min(2p,q)$ is achieved, $2p$ being the spatial order and $q$ the
temporal order. Universally among the different test cases the spatial
accuracy is seen to be dominant. A comparison with SSP-DG methods for the
Euler equations shows that MDRKCAT methods compare well with state-of-the-art
schemes in terms of accuracy.

A von Neumann analysis revealed that the stability of the MDRKCAT methods
depends heavily on the number of stages and the underlying high-order Lax-
Wendroff method. The latter one solely utilizes centered differences for the
spatial discretization. Consequently, odd-derivative Runge-Kutta schemes seem
less adequate in conjuction with the CAT algorithm.

In the future, there are two main routes to follow: extend and apply the
scheme to more challenging settings and to further examine the stability
properties of the novel scheme. Concerning more challenging settings the
investigation of multidimensional hyperbolic conservation laws with (possibly)
unstructured meshes and parabolic PDEs with viscous effects are attractive. In
order to accomplish such extensions it might be interesting to combine MDRKCAT
methods with DG techniques SSJ2017 . Presumably, also implicit MDRK schemes
need to be considered to take care of the diffusive effects. A possible
starting point could be the implicit variant of the approximate Taylor
methods, which have been recently developed for ODEs in 2020_Baeza_EtAl .
Concerning the stability properties of the scheme one could think of exploring
more types of MDRK schemes, possibly with SSP properties Seal2015b ;
2021_Gottlieb_EtAl . Moreover, at the same time, it will be possible to
identify more efficient schemes.

## Declarations

Funding J. Zeifang was funded by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the project no. 457811052. The HPC-
resources and services used in this work were provided by the VSC (Flemish
Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and
the Flemish Government.

Conflicts of interest The authors declare that they have no known competing
financial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Availability of data and material The datasets generated and/or analyzed
during the current study are available from the corresponding author on
reasonable request<EMAIL_ADDRESS>

Code availability The code used to generate the results in this work is
available upon reasonable request from the corresponding author
<EMAIL_ADDRESS>
##  Appendix A Butcher tableaux

In this section, we show the multiderivative Runge-Kutta methods used in this
work through their Butcher tableaux. We use three two-derivative methods taken
from TC10 , see Tbl. 3–5; two three-derivative methods taken from TurTur2017 ,
see Tbl. 6–7; and one four-derivative method, constructed for this paper, see
Tbl. 8. This last scheme has been derived from the idea that it should be of
form

$y^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}}=y^{n}+\sum\limits_{k=1}^{\mathtt{r}-1}\frac{\left(c_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color<EMAIL_ADDRESS>t}\right)^{k}}{k!}\Phi^{(k-1)}\left(y^{n}\right)+{\Delta
t}^{\mathtt{r}}\sum\limits_{{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}=1}^{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}-1}a_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}^{(\mathtt{r})}\Phi^{(\mathtt{r}-1)}\left(y^{n,{\color[rgb]{1,0.19,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.19,1}\pgfsys@color@cmyk@stroke{0}{0.81}{0}{0}\pgfsys@color@cmyk@fill{0}{0.81}{0}{0}\nu}}\right)\,,$
(24) for
${\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1,\dots,\mathtt{s}$,
with update $y^{n+1}=y^{n}+\sum\limits_{k=1}^{\mathtt{r}-1}\frac{{\Delta
t}^{k}}{k!}\Phi^{(k-1)}\left(y^{n}\right)+{\Delta
t}^{\mathtt{r}}\sum\limits_{{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}=1}^{\mathtt{s}}b^{(\mathtt{r})}_{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}\Phi^{(\mathtt{r}-1)}(y^{n,{\color[rgb]{0.15,1,0.67}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,1,0.67}\pgfsys@color@cmyk@stroke{0.85}{0}{0.33}{0}\pgfsys@color@cmyk@fill{0.85}{0}{0.33}{0}l}})\,.\hskip
20.00003pt$

These forms have also been used in TC10  and TurTur2017 .

Table 3: $2\text{DRK}3\text{-}2$ \- Third order two-derivative Runge-Kutta scheme using two stages TC10 . $0$ | $0$ | $0$ | $0$ | $0$
---|---|---|---|---
$1$ | $1$ | $0$ | $1/2$ | $0$
| $2/3$ | $1/3$ | $1/6$ | $0$
Table 4: $2\text{DRK}4\text{-}2$ \- Fourth order two-derivative Runge-Kutta scheme using two stages TC10 . $0$ | $0$ | $0$ | $0$ | $0$
---|---|---|---|---
$1/2$ | $1/2$ | $0$ | $1/8$ | $0$
| $1$ | $0$ | $1/6$ | $1/3$
Table 5: $2\text{DRK}5\text{-}3$ \- Fifth order two-derivative Runge-Kutta scheme using three stages TC10 . $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
---|---|---|---|---|---|---
$2/5$ | $2/5$ | $0$ | $0$ | $2/25$ | $0$ | $0$
$1$ | $1$ | $0$ | $0$ | $-1/4$ | $3/4$ | $0$
| $1$ | $0$ | $0$ | $1/8$ | $25/72$ | $1/36$
Table 6: $3\text{DRK}5\text{-}2$ \- Fifth order three-derivative Runge-Kutta scheme using two stages TurTur2017 . $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
---|---|---|---|---|---|---
$2/5$ | $2/5$ | $0$ | $2/25$ | $0$ | $4/375$ | $0$
| $1$ | $0$ | $1/2$ | $0$ | $1/16$ | $5/48$
Table 7: $3\text{DRK}7\text{-}3$ \- Seventh order three-derivative Runge-Kutta scheme using three stages TurTur2017 . The coefficients are given by $c_{2}=\frac{3-\sqrt{2}}{7}$, $c_{3}=\frac{3+\sqrt{2}}{7}$, $a_{32}^{(3)}=\frac{122+71\sqrt{2}}{7203}$, $b_{1}^{(3)}=\frac{1}{30}$, $b_{2}^{(3)}=\frac{1}{15}+\frac{13\sqrt{2}}{480}$, $b_{3}^{(3)}=\frac{1}{15}-\frac{13\sqrt{2}}{480}$.  $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
---|---|---|---|---|---|---|---|---|---
$c_{2}$ | $c_{2}$ | $0$ | $0$ | $c_{2}^{2}/2$ | $0$ | $0$ | $c_{2}^{3}/6$ | $0$ | $0$
$c_{3}$ | $c_{3}$ | $0$ | $0$ | $c_{3}^{2}/2$ | $0$ | $0$ | $c_{2}^{3}/6-a_{32}^{(3)}$ | $a_{32}^{(3)}$ | $0$
| $1$ | $0$ | $0$ | $1/2$ | $0$ | $0$ | $b_{1}^{(3)}$ | $b_{2}^{(3)}$ | $b_{3}^{(3)}$
Table 8: $4\text{DRK}6\text{-}2$ \- Sixth order four-derivative Runge-Kutta scheme using two stages. $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
---|---|---|---|---|---|---|---|---
$1/3$ | $1/3$ | $0$ | $1/18$ | $0$ | $1/162$ | $0$ | $1/1944$ | $0$
| $1$ | $0$ | $1/2$ | $0$ | $1/6$ | $0$ | $1/60$ | $1/40$

## References

  * [1] M. Abramowitz and I. A. Stegun, editors.  Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables.  Dover Publications, Inc., New York, 1964.
  * [2] A. Baeza, R. Bürger, M. d. C. Martí, P. Mulet, and D. Zorío.  On approximate implicit Taylor methods for ordinary differential equations.  Computational and Applied Mathematics, 39(4):304, 2020.
  * [3] H. Carrillo and C. Parés.  Compact approximate Taylor methods for systems of conservation laws.  Journal of Scientific Computing, 80(3):1832–1866, 2019.
  * [4] H. Carrillo, C. Parés, and D. Zorío.  Lax-Wendroff approximate Taylor methods with fast and optimized weighted essentially non-oscillatory reconstructions.  Journal of Scientific Computing, 86(1), 2021.
  * [5] R. Chan and A. Tsai.  On explicit two-derivative Runge-Kutta methods.  Numerical Algorithms, 53:171–194, 2010.
  * [6] A. J. Christlieb, S. Gottlieb, Z. J. Grant, and D. C. Seal.  Explicit strong stability preserving multistage two-derivative time-stepping schemes.  Journal of Scientific Computing, 68:914–942, 2016.
  * [7] B. Cockburn, G. E. Karniadakis, and C.-W. Shu.  The development of discontinuous Galerkin methods.  In B. Cockburn, G. E. Karniadakis, and C.-W. Shu, editors, discontinuous Galerkin Methods: Theory, Computation and Applications, volume 11 of Lecture Notes in Computational Science and Engineering, pages 3–50. Springer, 2000.
  * [8] M. Dumbser, D. S. Balsara, E. F. Toro, and C.-D. Munz.  A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes.  Journal of Computational Physics, 227(18):8209–8253, 2008.
  * [9] M. Dumbser, C. Enaux, and E. F. Toro.  Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws.  Journal of Computational Physics, 227(8):3971–4001, 2008.
  * [10] M. Dumbser, F. Fambri, M. Tavelli, M. Bader, and T. Weinzierl.  Efficient implementation of ADER discontinuous Galerkin schemes for a scalable hyperbolic PDE engine.  Axioms, 7(3), 2018.
  * [11] S. Gottlieb, Z. J. Grant, J. Hu, and R. Shu.  High order unconditionally strong stability preserving multi-derivative implicit and IMEX Runge–Kutta methods with asymptotic preserving properties.  arXiv preprint arXiv:2102.11939, 2021.
  * [12] S. Gottlieb, C.-W. Shu, and E. Tadmor.  Strong stability-preserving high-order time discretization methods.  SIAM Review, 43(1):89–112, 2001.
  * [13] W. Guo, J.-M. Qiu, and J. Qiu.  A new Lax–Wendroff discontinuous Galerkin method with superconvergence.  Journal of Scientific Computing, 65(1):299–326, 2015.
  * [14] E. Hairer and G. Wanner.  Solving ordinary differential equations II.  Springer Series in Computational Mathematics, 1991.
  * [15] A. Jaust, J. Schütz, and D. C. Seal.  Implicit multistage two-derivative discontinuous Galerkin schemes for viscous conservation laws.  Journal of Scientific Computing, 69:866–891, 2016.
  * [16] P. Lax and B. Wendroff.  Systems of conservation laws.  Communications on Pure and Applied Mathematics, 13(2):217–237, 1960\.
  * [17] R. J. LeVeque.  Numerical Methods for Conservation Laws.  Birkhäuser Basel, 1990.
  * [18] J. Li and Z. Du.  A two-stage fourth order time-accurate discretization for Lax–Wendroff type flow solvers I. Hyperbolic conservation laws.  SIAM Journal on Scientific Computing, 38(5):A3046–A3069, 2016.
  * [19] C. Lu and J. Qiu.  Simulations of shallow water equations with finite difference Lax-Wendroff weighted essentially non-oscillatory schemes.  Journal of Scientific Computing, 47(3):281–302, Jun 2011.
  * [20] MathWorks.  Symbolic Math Toolbox.  Natick, Massachusetts, United States, 2020.  URL https://www.mathworks.com/help/symbolic/.
  * [21] M. Ökten Turacı and T. Öziş.  Derivation of three-derivative Runge-Kutta methods.  Numerical Algorithms, 74(1):247–265, 2017.
  * [22] J. Qiu.  Development and comparison of numerical fluxes for LWDG methods.  Numerical Mathematics: Theory, Methods and Applications, 1(4):435–459, 2008.
  * [23] J. Qiu, M. Dumbser, and C.-W. Shu.  The discontinuous Galerkin method with Lax–Wendroff type time discretizations.  Computer Methods in Applied Mechanics and Engineering, 194(42-44):4528–4543, 2005.
  * [24] J. Qiu and C.-W. Shu.  Finite difference WENO schemes with Lax-Wendroff-type time discretizations.  SIAM Journal of Scientific Computing, 24(6):2185–2198, 2003.
  * [25] A. Quarteroni, R. Sacco, and F. Saleri.  Numerical Mathematics.  Springer, 2007.
  * [26] A. Ralston and P. Rabinowitz.  A First Course in Numerical Analysis.  Dover books on mathematics. Dover Publications, 2001.
  * [27] J. Schütz, D. Seal, and A. Jaust.  Implicit multiderivative collocation solvers for linear partial differential equations with discontinuous Galerkin spatial discretizations.  Journal of Scientific Computing, 73:1145–1163, 2017.
  * [28] J. Schütz, D. Seal, and J. Zeifang.  Parallel-in-time high-order multiderivative IMEX methods.  CMAT Preprint UP-21-01, http://www.uhasselt.be/Documents/CMAT/Preprints/2021/UP2101.pdf, 2021.
  * [29] T. Schwartzkopff, M. Dumbser, and C.-D. Munz.  ADER: A high-order approach for linear hyperbolic systems in 2d.  Journal of Scientific Computing, 17(1-4):231–240, 2002.
  * [30] D. Seal, Y. Güçlü, and A. Christlieb.  High-order multiderivative time integrators for hyperbolic conservation laws.  Journal of Scientific Computing, 60:101–140, 2014.
  * [31] R. Spiteri and S. Ruuth.  A new class of optimal high-order strong-stability-preserving time discretization methods.  SIAM Journal on Numerical Analysis, 40, 02 2002.
  * [32] J. C. Strikwerda.  Finite Difference Schemes and Partial Differential Equations, Second Edition.  Society for Industrial and Applied Mathematics, 2004.
  * [33] V. A. Titarev and E. F. Toro.  ADER: Arbitrary high order Godunov approach.  Journal of Scientific Computing, 17(1):609–618, 2002.
  * [34] V. A. Titarev and E. F. Toro.  ADER schemes for three-dimensional non-linear hyperbolic systems.  Journal of Computational Physics, 204(2):715–736, 2005.
  * [35] G. Whitham.  Linear and Nonlinear Waves.  Pure and Applied Mathematics: A Wiley Series of Texts, Monographs and Tracts. Wiley, 2011.
  * [36] D. Zorío, A. Baeza, and P. Mulet.  An approximate Lax-Wendroff-type procedure for high order accurate schemes for hyperbolic conservation laws.  Journal of Scientific Computing, 71:246–273, 2017.