# A CWENO large time-step scheme for Hamilton–Jacobi equations

E. Carlini, R. Ferretti, S. Preda, and M. Semplice

We propose a high order numerical scheme for time-dependent first order
Hamilton–Jacobi–Bellman equations. In particular we propose to combine a semi-
Lagrangian scheme with a Central Weighted Non-Oscillatory reconstruction. We
prove a convergence result in the case of state- and time-independent
Hamiltonians.

Numerical simulations are presented in space dimensions one and two, also for
more general state- and time-dependent Hamiltonians, demonstrating superior
performance in terms of CPU time gain compared with a semi-Lagrangian scheme
coupled with Weighted Non-Oscillatory reconstructions.

Keywords: semi-Lagrangian schemes, Hamilton–Jacobi equations, $\mathsf{CWENO}$
methods
AMS Subject Classification: Primary, 65N12, 65M10; Secondary, 49L25

##  1 Introduction

In this paper, we address the numerical approximation of the following
hyperbolic Hamilton–Jacobi–Bellman (HJB) equation:

$\begin{cases}v_{t}(t,x)+H(t,x,Dv(t,x))=0,&\text{for
}t,x\in(0,T)\times\mathbb{R}^{d},\\\ v(0,x)=v_{0}(x),&\text{for
}x\in\mathbb{R}^{d},\end{cases}$ (1.1)

where $v:(0,T)\times\mathbb{R}^{d}\to\mathbb{R}$, $Dv$ stands for the spatial
gradient, and $H:(0,T)\times\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$.

Numerous schemes have been proposed to approximate (1.1), but only a small
number are aimed at high order accuracy. High order finite difference schemes
based on $\mathsf{ENO}$ (Essentially Non-Oscillatory) reconstruction, defined
in [HEOC87], were proposed in [OS91] and extended to unstructured grids in
[Abg94]. Second order Godunov-type schemes based on global projection
operators are discussed in [LT00]. Weighted $\mathsf{ENO}$ ($\mathsf{WENO}$)
schemes were introduced in [JS96, LOC94] and combined with finite difference
schemes for Hamilton–Jacobi (HJ) equations in [JP00]. A fifth order central
scheme based on $\mathsf{WENO}$ reconstruction has been developed in [HT03].

Over the past two decades, semi-Lagrangian (SL) schemes have been employed to
discretize the HJ equation. In [FF94], a high order semi-discrete SL scheme is
proposed to discretize stationary HJB equations, while in [FF02], SL schemes
are applied to evolutionary HJB equations. A high order SL scheme for HJB
equations is presented in [CFR06] by combining the SL technique with
$\mathsf{WENO}$ reconstructions.

A different approach to achieve high order accuracy involves filtered schemes.
In [BFS16, FPT20a], filtered schemes are combined with monotone finite
difference schemes to discretize first order evolutionary HJ equations.

For an overview of numerical methods for first order Hamilton–Jacobi
equations, we also refer to the book [FF14].

We will develop our theory for the HJB equation related to a finite horizon
optimal control. More precisely, given a compact $A\in\mathbb{R}^{m}$, let the
set $\mathcal{A}=\\{\alpha:[0,T]\to A,\text{ measurable}\\}$ denote the
admissible controls. Given the running cost
$f_{C}:(0,T)\times\mathbb{R}^{d}\times A\to\mathbb{R}$ and the dynamics
$f_{D}:(0,T)\times\mathbb{R}^{d}\times A\to\mathbb{R}^{d}$ for the system

$\dot{y}(s)=f_{D}(t-s,y(s),\alpha(s)),\quad s\in(0,t],\quad y(0)=x,$ (1.2)

we consider the following Hamiltonian

$H(t,x,p)=\max_{a\in A}(-f_{D}(t,x,a)\cdot p-f_{C}(t,x,a)).$ (1.3)

Let us assume that

  * (H1)

$f_{C},f_{D}$ are continuous and bounded. Moreover, for every $a\in A$, the
functions $f_{C}(\cdot,\cdot,a)$, $f_{D}(\cdot,\cdot,a)$ are Lipschitz
continous, with Lipschitz constants independent of $a\in A$.

  * (H2)

$v_{0}$ is Lipschitz continous and bounded.

Under Assumptions (H1)–(H2), problem (1.1) admits a unique viscosity solution
$v$, which is Lipschitz continuous and bounded. Moreover, the Dynamic
Programming Principle holds (see [BCD97, Chapter 3], i.e., for any $h>0$

$v(t,x)=\underset{\alpha\in\mathcal{A}}{\inf}\left\\{\int_{0}^{h}f_{C}(t-s,y_{x,t}(s),\alpha(s)){\rm{d}}s+v(y_{x,t}(h),t-h)\right\\},$
(1.4)

where $y_{x,t}$ denotes the solution to (1.2).

A numerical scheme based on the previous description involves a minimization
procedure over the set of controls; in turn, each function evaluation requires
to approximate the cost integral and to interpolate the solution at the
previous time step at the foot $y_{x,t}(h)$ of the characteristics. In this
respect, $\mathsf{WENO}$ reconstructions, especially in higher space
dimensions, are not efficient, due to fact that the linear weights of
$\mathsf{WENO}$ depend on the reconstruction point.

A novel paradigm for non-oscillatory reconstruction operators has been
introduced in [LPR00], where the authors suggested to blend, in a
$\mathsf{WENO}$-like fashion, polynomials of different degrees, allowing to
overcome some difficulties of non-existence, non-positivity and dependence on
the reconstruction point of the $\mathsf{WENO}$ linear weights. The idea has
been further developed into the so-called $\mathsf{CWENO}$ reconstruction and
exploited in more spatial dimensions, also in the case of adaptive mesh
refinement and non-uniform grids [BGFB20, ZQ16, BBMZ19, ZCZ08, DBSR17]. The
technique has also been exploited in finite difference schemes for Hamilton-
Jacobi equations on Cartesian meshes via dimensional splitting [ZQ17, ZSQ19]
and on general meshes [ZQ20]. General results for establishing the convergence
order of a $\mathsf{CWENO}$ reconstruction have been presented in [CPSV18,
CSV19, SV20].

One of the main advantages of $\mathsf{CWENO}$ over the traditional
$\mathsf{WENO}$ is that the central approach provides a reconstruction
polynomial that is defined everywhere in the reconstruction cell and that can
be evaluated later, with no essential extra cost, at many different
reconstruction points. This is guaranteed by the independence of the linear
weights from the reconstruction point. Furthermore, as shown in [SCR16,
DBSR17, CS19, CPSV18], the $\mathsf{CWENO}$ approach allows to avoid the
dimensional splitting procedure and this is advantageous also on Cartesian
meshes when the number of reconstruction points per cell is high.

In this paper we want to exploit the positive features of $\mathsf{CWENO}$ to
obtain a high order SL scheme that is more efficient than the one of [CFR06],
which are based on $\mathsf{WENO}$. The rest of the paper is organized as
follows. In §2 the general principles of the SL scheme are recalled, whereas
in §3 a $\mathsf{CWENO}$ interpolation from point values is constructed in one
and two space dimensions. §4 proves a convergence result in the framework of
[Fer02], and finally §5 provides an extensive number of numerical examples to
validate the scheme.

##  2 Numerical scheme

To obtain an approximate version of (1.1), we need first to discretize the
control problem in time. Let $\mathrm{\Delta t}>0$ be a time step,
$t_{n}=n\mathrm{\Delta t}$ a uniform time grid with $n=0,\ldots,N_{T}$, where
$N_{T}=\lceil\frac{T}{\mathrm{\Delta t}}\rceil$, and consider the Dynamic
Programming Principle on a single time step $[t_{n},t_{n+1}]$, by choosing in
(1.4) $h=\mathrm{\Delta t}$ and $t=t_{n+1}$.

When using a $\nu$-stages scheme for (1.2), we discretize an admissible
control $\alpha\in\mathcal{A}$ via a sequence
$\underline{a}=(a_{1},a_{2},\ldots,a_{\nu})\in A^{\nu}$. Then the foot of the
characteristic corresponding to the control $\underline{a}$ is

$y^{n}(x,\underline{a})\simeq y_{x,t_{n+1}}(\mathrm{\Delta t})$

for any $n=0,\ldots,N_{T}-1$, $x\in{\mathbb{R}}^{d}$ and $\underline{a}\in
A^{\nu}$.

Moreover, a suitable quadrature formula based on $\nu$ quadrature nodes is
introduced to discretize the integral term related to the running cost,

$C^{n}(x,\underline{a})\simeq\int_{0}^{\mathrm{\Delta
t}}f_{C}(t_{n+1}-s,y_{x,t_{n+1}}(s),\alpha(s))\,{\rm d}s,$

for any $n=0,\ldots,N_{T}-1$, $x\in{\mathbb{R}}^{d}$ and $\underline{a}\in
A^{\nu}$ .

Let us now introduce the space grid with space step $\mathrm{\Delta x}$ so
that $x_{i}=i\mathrm{\Delta x}$ with $i\in\mathbb{Z}^{d}$ and denote
$y_{i}^{n}(\underline{a})=y^{n}(x_{i},\underline{a})$ and
$C^{n}_{i}(\underline{a})=C^{n}(x_{i},\underline{a})$. We approximate the
solution $v(x,t)$ of (1.1) by a discrete function $u^{n}_{i}\simeq
v(x_{i},t_{n})$ defined on the space-time grid, computed by the following
iterative scheme, for $n=0,\dots,N_{T}-1$

$\begin{cases}u^{n+1}_{i}=\underset{\underline{a}\in
A^{\nu}}{\min}\\{R[u^{n}](y^{n}_{i}(\underline{a}))+C^{n}_{i}(\underline{a})\\}&i\in\mathbb{Z}^{d},\\\
u^{0}_{i}=v_{0}(x_{i}),&i\in\mathbb{Z}^{d},\end{cases}$ (2.1)

where $R[u^{n}](x)$ denotes the spatial reconstruction at
$x\in{\mathbb{R}}^{d}$ of the numerical solution
$u^{n}=(u^{n}_{i})_{i\in\mathbb{Z}^{d}}$.

For a given $\underline{a}\in A^{\nu}$, $y^{n}_{i}(\underline{a})$ is computed
via a $\nu$-stage Runge–Kutta (RK) methods as follows:

$\displaystyle y^{n}_{i}(\underline{a})$ $\displaystyle=x_{i}+\mathrm{\Delta
t}\sum_{k=1}^{\nu}b_{k}K_{k}(\underline{a}),$ (2.2) $\displaystyle
K_{k}(\underline{a})$ $\displaystyle=f_{D}\,(t_{n}+(1-c_{k})\mathrm{\Delta
t},X_{k}(\underline{a}),a_{k}),$ $\displaystyle X_{k}(\underline{a})$
$\displaystyle=x_{i}+\mathrm{\Delta
t}\sum_{j=1}^{k-1}A_{kj}K_{j}(\underline{a}),$

where $b_{k}$,$c_{k}$,$A_{kj}$ are the coefficients of the Butcher tableau
defining the RK method used to solve (1.2). For our purposes, we resort to the
Forward Euler method for first order schemes, to Heun’s method

$\begin{array}[]{c|cc}0\\\ 1&1\\\
\hline\cr&\nicefrac{{1}}{{2}}&\nicefrac{{1}}{{2}}\end{array}$ (2.3)

to get second order accuracy, and to RK method with tableau

$\begin{array}[]{c|ccc}0\\\ \nicefrac{{1}}{{2}}&\nicefrac{{1}}{{2}}&\\\
1&-1&2&\\\
\hline\cr&\nicefrac{{1}}{{6}}&\nicefrac{{2}}{{3}}&\nicefrac{{1}}{{6}}\end{array}$
(2.4)

to get third order accuracy.

In order to compute $C^{n}_{i}(\underline{a})$, we employ a suitable
quadrature rule and replace the function evaluations at the nodes with the
numerical solution of (1.2):

$\displaystyle\int_{0}^{\mathrm{\Delta
t}}f_{C}(t_{n+1}-s,y_{x_{i},t_{n+1}}(s),\alpha(s))ds$ (2.5) $\displaystyle=$
$\displaystyle\mathrm{\Delta
t}\sum_{k}w_{k}f_{C}(t_{n+1}-\xi_{k}\mathrm{\Delta
t},y_{x_{i},t_{n+1}}(\xi_{k}\mathrm{\Delta t}),\alpha(\xi_{k}\mathrm{\Delta
t}))+\mathcal{O}(\mathrm{\Delta t}^{Q+1})$ $\displaystyle=$
$\displaystyle\mathrm{\Delta
t}\sum_{k}w_{k}f_{C}(t_{n+1}-\xi_{k}\mathrm{\Delta
t},\tilde{y}_{x_{i},t_{n+1}}(\xi_{k}\mathrm{\Delta
t}),\alpha(\xi_{k}\mathrm{\Delta t}))+\mathcal{O}(\mathrm{\Delta
t}^{Q+1}+\mathrm{\Delta t}^{R+1})$

where $\xi_{k},w_{k}$ are the nodes and the weights of the quadrature rule and
$\tilde{y}_{x_{i},t_{n+1}}(\xi_{k}\mathrm{\Delta t})$ is a suitable
approximation of the numerical solution of $y_{x_{i},t_{n+1}}(s)$ for
$s\in[0,\mathrm{\Delta t}]$; $Q$ is the accuracy of the quadrature rule and
$R$ describes the extra error introduced by the approximation of
$y_{x_{i},t_{n+1}}(\xi_{k}\mathrm{\Delta t})$ with
$\tilde{y}_{x_{i},t_{n+1}}(\xi_{k}\mathrm{\Delta t})$ within the quadrature
rule. It is tempting to choose a quadrature rule whose nodes coincide with the
abscissae $\\{c_{k}\\}$ of the RK scheme and replace, in the last
approximation of (2.5), $\tilde{y}_{x_{i},t_{n+1}}(\xi_{k}\mathrm{\Delta t})$
with the stage values of the RK scheme $X_{k}$ and
$\alpha(\xi_{k}\mathrm{\Delta t})$ with $a_{k}$, for $k=1,\ldots,\nu$.

For the Heun method, the numerical computation of $C^{n}_{i}(\underline{a})$
would then be performed with the trapezoidal rule, whose nodes are at
$\xi=0,1$ and for which $\tilde{y}_{x_{i},t_{n+1}}(0)=x_{i}$ and
$\tilde{y}_{x_{i},t_{n+1}}(\mathrm{\Delta t})=X_{2}$ are computed with the
correct accuracy by the RK scheme (2.3).

For the third order RK scheme (2.4), whose nodes are at
$\xi=~{}0,\nicefrac{{1}}{{2}},1$, one would then choose the approximations
$\tilde{y}_{x_{i},t_{n+1}}(0)=x_{i}$,
$\tilde{y}_{x_{i},t_{n+1}}(\frac{1}{2}\mathrm{\Delta t})=X_{2}$ and
$\tilde{y}_{x_{i},t_{n+1}}(\mathrm{\Delta t})=X_{3}$. In this special case,
even if both $y_{x_{i},t_{n+1}}(\xi_{2}\mathrm{\Delta t})-X_{2}$ and
$y_{x_{i},t_{n+1}}(\xi_{3}\mathrm{\Delta t})-X_{3}$ are
$\mathcal{O}(\mathrm{\Delta x}^{2})$, these errors cancel each other out when
computing the linear combination of the Simpson’s quadrature rule, resulting
in a third order accuracy for the approximation of $C_{i}^{n}(\underline{a})$.
For more general RK processes, this issue can be overcome by resorting to the
continuous extension [Zen86] to compute $\tilde{y}_{x_{i},t_{n+1}}(\xi)$ with
the proper accuracy.

To compute the solution of (2.1), we employ tabulation on a coarse grid in
$A$, followed by a Nelder-Mead algorithm, adjusted so that no vertex of the
simplices can exit the compact set $A$. At each step of the minimization we
need to evaluate the reconstruction at the feet of the characteristics
$y^{n}_{i}(\underline{a})$ for each node $x_{i}$ and for each discrete control
$\underline{a}$. Thus an efficient and accurate reconstruction operator is
crucial to obtain a robust and fast numerical scheme. Since (2.1) has in
general nonsmooth solutions, standard high order interpolation give rise to
oscillations in the numerical solution. [CFR06], a $\mathsf{WENO}$
reconstruction has been applied to overcome this problem. The $\mathsf{CWENO}$
approach, instead, is more efficient when many evaluations of the
reconstruction in each given cell are required. For this reason, we propose to
apply $\mathsf{CWENO}$ interpolation within our semi-Lagrangian scheme (2.1).

##  3 $\mathsf{CWENO}$ reconstruction

We recall here the definition of the $\mathsf{CWENO}$ and $\mathsf{CWENOZ}$
operators, following the presentation of [CPSV18, CSV19]; in particular,
despite the present finite-difference setting in place of the finite-volume
one, we will still be able to exploit the general theorems for the accuracy of
the reconstructions on smooth solutions proven in the aforementioned papers.

For a point $x\in{\mathbb{R}}^{d}$, let $\Omega$ be the grid cell containing
it and $\mathcal{S}_{\Omega}$ be the set of its vertices. In order to achieve
better than first order accuracy, we need to consider stencils
$\mathcal{S}\supset\mathcal{S}_{\Omega}$. We associate to any $\mathcal{S}$ a
polynomial $P^{(r)}_{\mathcal{S}}(x)\in\mathbb{P}^{r}_{d}$ which interpolates
the data in $\mathcal{S}$, i.e. such that $P^{(r)}_{\mathcal{S}}(x_{j})=u_{j}$
$\forall j\in\mathcal{S}$. Larger symmetric stencils will define highly
accurate interpolators, which however would be very oscillatory when the data
in $\mathcal{S}$ represent a non smooth function. Instead, polynomials
associated to smaller stencils, biased in a specific direction, could avoid
the discontinuities in the data.

The general idea of $\mathsf{CWENO}$ and $\mathsf{CWENOZ}$ reconstruction is
to blend, in a nonlinear and data-dependent fashion, a high order accurate
interpolating polynomial with a set of lower order ones to produce an
essentially non-oscillatory reconstruction polynomial for the cell $\Omega$.
This can be later evaluated at any point in $x\in\Omega$ with negligible cost.

The nonlinear selection or blending of polynomials relies on oscillation
indicators $\mathrm{I}[P]$, which are in general scalar quantities associated
to a polynomial $P$, designed in such a way that $\mathrm{I}[P]\rightarrow 0$
under grid refinement, if $P$ is associated to smooth data, and
$\mathrm{I}[P]\asymp 1$, in presence of a discontinuity in the stencil
$\mathcal{S}$. In this work we rely on the classical Jiang–Shu oscillation
indicators [JS96], suitably modified to accomodate for the regularity of the
solution of HJB problems under consideration [JP00, FPT20b].

We recall now the definition of the $\mathsf{CWENO}$ and $\mathsf{CWENOZ}$
reconstructions given in [CSV19].

###### Definition 3.1

Given a stencil $\mathcal{S}_{\text{opt}}$, including $\mathcal{S}_{\Omega}$,
let $P_{\text{opt}}\in\mathbb{P}^{G}_{n}$ (optimal polynomial) be the
polynomial of degree $G$, associated to $\mathcal{S}_{\text{opt}}$. Further,
let $P_{1},P_{2},...,P_{m}$ be a set of $m\geq 1$ polynomials of degree $g$
with $g<G$, associated to substencil $\mathcal{S}_{k}$ such that
$\mathcal{S}_{\Omega}\subset\mathcal{S}_{k}\subset\mathcal{S}_{\text{opt}}$
$\forall k=1,...,m$. Let also $\left\\{d_{k}\right\\}_{k=0}^{m}$ be a set of
strictly positive real coefficients such that $\sum_{k=0}^{m}d_{k}=1$.

The $\mathsf{CWENO}$ and $\mathsf{CWENOZ}$ operators compute a reconstruction
polynomial

$\displaystyle P_{\text{rec}}^{CW}$
$\displaystyle=\mathsf{CWENO}(P_{\text{opt}},P_{1},...,P_{m})\in\mathbb{P}^{G}_{n},$
(3.1) $\displaystyle P_{\text{rec}}^{CWZ}$
$\displaystyle=\mathsf{CWENOZ}(P_{\text{opt}},P_{1},...,P_{m})\in\mathbb{P}^{G}_{n},$

as follows:

  1. 1.

First, introduce the polynomial $P_{0}$ defined as

$P_{0}(x)=\frac{1}{d_{0}}\left(P_{\text{opt}}(x)-\sum_{k=1}^{m}d_{k}P_{k}(x)\right)\in\mathbb{P}^{G}_{n};$
(3.2)

  2. 2.

compute suitable regularity indicators

$\mathrm{I}_{0}=\mathrm{I}[P_{\text{opt}}],\quad\mathrm{I}_{k}=\mathrm{I}[P_{k}],\;k\geq
1;$ (3.3)

  3. 3.

compute the nonlinear coefficients $\left\\{\omega_{k}\right\\}_{k=0}^{m}$ or
$\left\\{\omega^{Z}_{k}\right\\}_{k=0}^{m}$ as

    1. (a)

$\mathsf{CWENO}$ operator: for $k=0,...,m$,

$\alpha_{k}=\frac{d_{k}}{(\mathrm{I}_{k}+\epsilon)^{l}},\quad\omega_{k}^{CW}=\frac{\alpha_{k}}{\sum_{i=0}^{m}\alpha_{i}},$
(3.4)

    2. (b)

$\mathsf{CWENOZ}$ operator: for $k=0,...,m$,

$\alpha_{k}^{Z}=d_{k}\left(1+\left(\frac{\tau}{\mathrm{I}_{k}+\epsilon}\right)^{l}\right),\quad\omega_{k}^{CWZ}=\frac{\alpha_{k}^{Z}}{\sum_{i=0}^{m}\alpha_{i}^{Z}},$
(3.5)

where $\epsilon$ is a small positive quantity, $l\geq 1$, and, in the case of
$\mathsf{CWENOZ}$, $\tau$ is a global smoothness indicator;

  4. 4.

finally, define the reconstruction polynomial as

$P_{\text{rec}}^{CW}(x)=\sum_{k=0}^{m}\omega_{k}^{CW}P_{k}(x)\in\mathbb{P}^{G}_{n},$
(3.6)
$P_{\text{rec}}^{CWZ}(x)=\sum_{k=0}^{m}\omega_{k}^{CWZ}P_{k}(x)\in\mathbb{P}^{G}_{n}.$
(3.7)

Note that the reconstruction polynomial defined in (3.6) and in (3.7) can be
evaluated at any point in the computational cell $\Omega$ at a very low
computational cost and this happens because the reconstruction procedure, in
both cases, starts with the definition of the linear coefficients
$\left\\{d_{k}\right\\}_{k=0}^{m}$ that do not depend on the reconstruction
point. Therefore, the nonlinear coefficients given by (3.4) and (3.5) are
computed once per cell and not once per reconstruction point, as in the
standard $\mathsf{WENO}$. We also remark that, since every polynomial taking
place in the reconstruction procedures is required to satisfy the
interpolation constraint on $\mathcal{S}_{\Omega}$, then also
$P_{\text{rec}}^{CW}$ and $P_{\text{rec}}^{CWZ}$ satisfy the same constraint
on the vertices of $\Omega$.

The accuracy and non-oscillatory properties of $\mathsf{CWENO}$ and
$\mathsf{CWENOZ}$ schemes are guaranteed by the dependence of their nonlinear
weights (3.4) and (3.5) on the regularity indicators $\mathrm{I}_{k}$. On
smooth data, the nonlinear weights are driven sufficiently close to the
optimal ones, so that $P_{\text{rec}}\approx P_{\text{opt}}$ and the
reconstruction reaches the optimal order of accuracy $G+1$. On the other hand,
when a discontinuity is present in $\mathcal{S}_{\text{opt}}$, both
$\mathrm{I}_{0}\asymp 1$ and at least one $\mathrm{I}_{\hat{k}}\asymp 1$ for
some $\hat{k}\in\\{1,...,m\\}$. Then the formulas (3.4) and (3.5) for the
nonlinear weights will ensure that $\omega_{0}\approx 0$ and
$\omega_{\hat{k}}\approx 0$ for all $\hat{k}$ such that $P_{\hat{k}}$ would
bring oscillations in the reconstruction. The reconstruction polynomial will
then be a linear combination of all polynomials of degree $g$ that are not
affected by the discontinuity; the accuracy of the reconstruction thus reduces
to $g+1$, but spurious oscillations would be tamed.

###  3.1 One spatial dimension

Let us describe the reconstruction for any point $x$ in the cell
$\Omega=[x_{j},x_{j+1}]$, so that $\mathcal{S}_{\Omega}=\\{j,j+1\\}$. We
consider $\mathcal{S}_{\text{opt}}=\\{j-1,j,j+1,j+2\\}$ and introduce the
optimal cubic polynomial
$P_{\text{opt}}(x)=Q(x)=\sum_{i=0}^{3}z_{i}[(x-x_{j})/\mathrm{\Delta x}]^{i}$
that interpolates the data $u_{j-1},u_{j},u_{j+1},u_{j+2}$ at the nodes
$x_{j-1},x_{j},x_{j+1},x_{j+2}$. Next, we consider two parabolas $P_{L}(x)$
and $P_{R}(x)$ that interpolate only the data $u_{j-1},u_{j},u_{j+1}$ and,
respectively, $u_{j},u_{j+1},u_{j+2}$.

The reconstruction operators thus compute

$\displaystyle P_{\text{rec}}^{CW}$
$\displaystyle=\mathsf{CWENO}(Q,P_{L},P_{R})\in\mathbb{P}^{3}_{1},$ (3.8)
$\displaystyle P_{\text{rec}}^{CWZ}$
$\displaystyle=\mathsf{CWENOZ}(Q,P_{L},P_{R})\in\mathbb{P}^{3}_{1}.$

For any polynomial, its oscillation indicator is defined as

$\mathrm{I}[P]=\sum_{\alpha\geq 2}\mathrm{\Delta
x}^{2\alpha-3}\int_{x_{j}}^{x_{j+1}}\left(\frac{d^{(\alpha)}P}{dx^{\alpha}}\right)^{2}dx.$
(3.9)

The above definition is similar to the classical definition of the oscillation
indicators for the $\mathsf{WENO}$ reconstruction as given in [JS96], except
that the first derivative is not included in the sum. This choice is the
appropriate one for Hamilton–Jacobi equations, whose solution can be at worst
continuous with kinks, see [JP00], and ensures that at worst
$\mathrm{I}[P]=\mathcal{O}(1)$.

Let $P(x)=\sum_{i=0}^{3}z_{i}[(x-x_{j})/\mathrm{\Delta x}]^{i}$ be a
polynomial of degree up to $3$. Its indicator (3.9) can be written as a
quadratic form of its coefficients given by

$\mathrm{I}[P]=\frac{1}{\mathrm{\Delta
x}^{2}}(4z_{2}^{2}+12z_{2}z_{3}+48z_{3}^{2})$

or equivalently, denoting with $\vec{z}$ the vector of coefficients,

$\mathrm{I}[P]=\frac{1}{\mathrm{\Delta x}^{2}}\vec{z}^{T}M\vec{z}\text{,
with}\quad M=\begin{pmatrix}0&0&0&0\\\ 0&0&0&0\\\ 0&0&4&6\\\ 0&0&6&48\\\
\end{pmatrix}.$ (3.10)

Denoting by $U$ the vector of data $(u_{j-1},u_{j},u_{j+1},u_{j+2})^{T}$,
since the coefficients of the polynomial linearly depends on $U$, we can
express the regularity indicators also as a quadratic form on the data being
interpolated. More precisely, the vector of coefficients $\vec{z}$ can be
expressed as $\vec{z}=V^{-1}U$, where $V^{-1}$ is the inverse of the
Vandermonde matrix. Thus, in the one-dimensional case, one obtains that the
matrices of the quadratic forms expressed in terms of $U$ scale again globally
as $1/\mathrm{\Delta x}^{2}$:

$\mathrm{I}[Q]=\frac{1}{\mathrm{\Delta
x}^{2}}U^{T}A_{Q}U,\quad\mathrm{I}[P_{k}]=\frac{1}{\mathrm{\Delta
x}^{2}}U^{T}A_{k}U\quad(k=L,R).$ (3.11)

In our case, we get

$A_{Q}=\begin{pmatrix}4/3&-7/2&3&-5/6\\\ -7/2&10&-19/2&3\\\ 3&-19/2&10&-7/2\\\
-5/6&3&-7/2&4/3\\\ \end{pmatrix},$ $A_{L}=\begin{pmatrix}1&-2&1&0\\\
-2&4&-2&0\\\ 1&-2&1&0\\\ 0&0&0&0\\\ \end{pmatrix},\hskip
5.69046ptA_{R}=\begin{pmatrix}0&0&0&0\\\ 0&1&-2&1\\\ 0&-2&4&-2\\\ 0&1&-2&1\\\
\end{pmatrix}.$

The $\mathsf{CWENO}$ reconstruction applied in the numerical tests of this
paper is then defined by choosing linear coefficients
$d_{L}=d_{R}=\nicefrac{{1}}{{8}}$ and $d_{0}=\nicefrac{{3}}{{4}}$, $l=2$ and
$\epsilon=\mathrm{\Delta x}^{2}$ in the above construction.

For the optimal definition of the $\mathsf{CWENOZ}$ reconstruction, in order
to exploit the results of [CSV19], we need to study the Taylor expansions of
the indicators centered at the point $x_{0}=(x_{j}+x_{j+1})/2$, namely the
center of the reconstruction cell $\Omega$. They are given by

$\displaystyle\mathrm{I}[P_{L}]$
$\displaystyle=B-u^{\prime\prime}(x_{0})u^{\prime\prime\prime}(x_{0})\mathrm{\Delta
x}^{3}+\frac{1}{4}u^{\prime\prime\prime}(x_{0})^{2}\mathrm{\Delta
x}^{4}+\mathcal{O}(\mathrm{\Delta x}^{5}),$ (3.12)
$\displaystyle\mathrm{I}[P_{R}]$
$\displaystyle=B+u^{\prime\prime}(x_{0})u^{\prime\prime\prime}(x_{0})\mathrm{\Delta
x}^{3}+\frac{1}{4}u^{\prime\prime\prime}(x_{0})^{2}\mathrm{\Delta
x}^{4}+\mathcal{O}(\mathrm{\Delta x}^{5}),$ $\displaystyle\mathrm{I}[P_{L}]$
$\displaystyle=B+\frac{13}{12}u^{\prime\prime\prime}(x_{0})^{2}\mathrm{\Delta
x}^{4}+\mathcal{O}(\mathrm{\Delta x}^{5}),$

where $B=u^{\prime\prime}(x_{0})^{2}\mathrm{\Delta x}^{2}$. Thus, defining

$\tau=\bigg{|}2\mathrm{I}[Q]-\mathrm{I}[P_{L}]-\mathrm{I}[P_{R}]\bigg{|},$
(3.13)

all terms up to $\mathcal{O}(\mathrm{\Delta x}^{3})$ cancel. Note that this is
the optimal definition of $\tau$ since it is never possible to cancel all the
$\mathcal{O}(\mathrm{\Delta x}^{4})$ terms in a convex combination of the
three indicators. The hypotheses of Theorem $24$ in [CSV19] hold and we can
guarantee the optimal order of accuracy of the $\mathsf{CWENOZ}$
reconstruction for $l=2$.

###  3.2 Two spatial dimensions

In two space dimensions, we will not rely on dimensional splitting. Such an
approach would in fact nullify the advantages of a fast evaluation of the
reconstruction polynomial on a large number of arbitrary points in the
reconstruction cell. Our reconstruction will instead generate a polynomial of
two variables, like it is done in [SCR16, CS19] for Cartesian meshes or in
general mesh settings [ZQ20, DBSR17, BGFB20].

Let $\Omega$ be the cell of the Cartesian grid containing the reconstruction
point $x$. The setup of the stencils for the reconstruction is illustrated in
Fig. 1: the cell $\Omega$, with corners $(x_{i},y_{j})$ and
$(x_{i+1},y_{j+1})$, is hatched in red and the reconstruction stencil is
shaded in gray.

Figure 1: Stencils of the two-dimensional $\mathsf{CWENO}$ and
$\mathsf{CWENOZ}$ reconstructions. The red hatched region represents the cell
$\Omega_{i,j}$ in which we compute the reconstruction. The vertices of
$\mathcal{S}_{\text{opt}}$ are enclosed in the grey shaded region, while the
stencils for the low degree polynomials are enclosed in the coloured squares:
blue, green, orange and purple, respectively for the north-east, north-west,
south-west and south-east polynomial.

The reconstruction operators compute

$\displaystyle P_{\text{rec}}^{CW}$
$\displaystyle=\mathsf{CWENO}(Q^{(3)};Q^{(2)}_{\text{ne}},Q^{(2)}_{\text{se}},Q^{(2)}_{\text{sw}},Q^{(2)}_{\text{nw}})\in\mathbb{P}^{3}_{1}\otimes\mathbb{P}^{3}_{1},$
(3.14) $\displaystyle P_{\text{rec}}^{CWZ}$
$\displaystyle=\mathsf{CWENOZ}(Q^{(3)};Q^{(2)}_{\text{ne}},Q^{(2)}_{\text{se}},Q^{(2)}_{\text{sw}},Q^{(2)}_{\text{nw}})\in\mathbb{P}^{3}_{1}\otimes\mathbb{P}^{3}_{1},$

with the optimal polynomial set to the bicubic polynomial interpolating the 16
data points in the square with corners $(x_{i-1},y_{j-1})$ and
$(x_{i+2},y_{j+2})$ enclosed in the grey shaded region of Fig 1. The other
four polynomials are the biquadratic polynomials that interpolate the four
$3\times 3$ substencils in the north-east, south-east, south-west and north-
west directions, as illustrated in Fig. 1, named $\mathcal{S}_{\text{ne}}$,
$\mathcal{S}_{\text{nw}}$, $\mathcal{S}_{\text{sw}}$ and
$\mathcal{S}_{\text{se}}$, respectively. In particular, the blue box with
vertices $(x_{i},y_{j})$ and $(x_{i+2},y_{j+2})$ encloses
$\mathcal{S}_{\text{ne}}$, the green box with vertices $(x_{i-1},y_{j})$ and
$(x_{i+1},y_{j+2})$ encloses $\mathcal{S}_{\text{nw}}$, the orange box with
vertices $(x_{i-1},y_{j-1})$ and $(x_{i+1},y_{j+1})$ encloses
$\mathcal{S}_{\text{sw}}$ and the purple box with vertices $(x_{i},y_{j-1})$
and $(x_{i+2},y_{j+1})$ encloses $\mathcal{S}_{\text{se}}$. The basis for
bicubic and biquadratic polynomials are defined by tensorization and are thus
respectively given by

$\displaystyle\mathcal{B}_{3}$
$\displaystyle=\\{1,\hat{x},\hat{y},\hat{x}^{2},\hat{x}\hat{y},\hat{y}^{2},\hat{x}^{3},\hat{x}^{2}\hat{y},\hat{x}\hat{y}^{2},\hat{y}^{3},\hat{x}^{3}\hat{y},\hat{x}^{2}\hat{y}^{2},\hat{x}\hat{y}^{3},\hat{x}^{3}\hat{y}^{2},\hat{x}^{2}\hat{y}^{3},\hat{x}^{3}\hat{y}^{3}\\},$
$\displaystyle\mathcal{B}_{2}$
$\displaystyle=\\{1,\hat{x},\hat{y},\hat{x}^{2},\hat{x}\hat{y},\hat{y}^{2}\\},$

where $\hat{x}=(x-x_{i})/\mathrm{\Delta x}$ and
$\hat{y}=(y-y_{j})/\mathrm{\Delta y}$.

The oscillation indicators are defined as

$\mathrm{I}[P]=\sum_{|\vec{\alpha}|\geq
2}\Delta^{2|\vec{\alpha}|-4}\int_{\Omega_{i,j}}\left[\frac{\partial
P}{\partial x^{\alpha_{1}}\partial
y^{\alpha_{2}}}\right]^{2}\hbox{d}x\hbox{d}y$ (3.15)

where $\Delta$ is the diameter of the cell $\Omega_{i,j}$. These indicators
are inspired to the classical ones of [HS99], but the powers of $\Delta$
appearing in (3.15) have been adjusted for the expected regularity of the
solution to Hamilton–Jacobi equations. In particular the scaling
$\Delta^{2|\vec{\alpha}|-4}$ is a generalization to two space dimensions of
the choice of [JP00]; it is introduced so that $\mathrm{I}[P]\asymp 1$ for
functions with discontinous first derivative and
$\mathrm{I}[P]=\mathcal{O}(\Delta^{2})$ on regular solutions. A detailed
discussion of this scaling can be found in [FPT20b].

Denoting with $\vec{z}$ the vector of coefficients along the basis
$\mathcal{B}_{3}$ of a generic polynomial in two variables of degree up to 3,
and considering a uniform grid of size $\mathrm{\Delta x}$, its oscillation
indicator can be expressed as

$\mathrm{I}[P]=\frac{1}{\mathrm{\Delta x}^{2}}\vec{z}^{T}M\vec{z}$ (3.16)

for a $16\times 16$ matrix $M$We point out that in the case of Cartesian but
uneven grids, there would still be a global factor inversely proportional to
the local grid size, but matrix $M$ would also depend on the aspect ratios and
on the ratios of sizes of nearby cells.

Consider now the stencil of the reconstruction given by
$\\{(x_{i-1+k},y_{j-1+l})\\}_{\tiny{\begin{subarray}{c}l=0,\ldots,3\\\
k=0,\ldots,3\end{subarray}}}$ and denote by $U$ the corresponding vector of
values $\\{(u_{i-1+k},y_{j-1+l})\\}_{\tiny{\begin{subarray}{c}l=0,\ldots,3\\\
k=0,\ldots,3\end{subarray}}}$ with a prescribed ordering, e.g. lexicographic.
Thanks to the linear relation among the coefficients $\vec{z}$ and the data
vector $U$, the regularity indicators can be expressed in the form

$\mathrm{I}[Q]=\frac{1}{\mathrm{\Delta
x}^{2}}U^{T}A_{Q}U,\quad\mathrm{I}[Q_{k}]=\frac{1}{\mathrm{\Delta
x}^{2}}U^{T}A_{k}U\quad(k=\text{ne},\text{se},\text{sw},\text{nw}).$ (3.17)

The matrices $A_{Q}$ and $A_{k}$ will in general depend on the local aspect
ratio of cells and on the neighbours size ratio.

The $\mathsf{CWENO}$ reconstruction is then defined by choosing linear
coefficients $d_{k}=1/16$ and $d_{0}=1-\sum_{k}d_{k}$ with
$k\in\\{\text{ne},\text{se},\text{sw},\text{nw}\\}$, $l=2$ and
$\epsilon=\mathrm{\Delta x}^{2}$.

The $\mathsf{CWENOZ}$ reconstruction is computed choosing $l=2$ and

$\tau=\left|4\mathrm{I}[Q^{(3)}]-\mathrm{I}[Q^{(2)}_{\text{ne}}]-\mathrm{I}[Q^{(2)}_{\text{se}}]-\mathrm{I}[Q^{(2)}_{\text{sw}}]-\mathrm{I}[Q^{(2)}_{\text{nw}}]\right|$
(3.18)

in order to guarantee the optimal order of accuracy of the $\mathsf{CWENOZ}$
reconstruction, see [CSV19, Theorem 24]. In fact, computing the Taylor
expansions of the indicators centered at the center of the cell $\Omega$, one
gets, for smooth enough data,

$\displaystyle\mathrm{I}[Q^{(2)}_{\text{ne}}]$
$\displaystyle=B+2u_{xx}u_{xxx}\mathrm{\Delta
x}^{3}+2u_{yy}u_{yyy}\mathrm{\Delta x}^{3}+\mathcal{O}(\mathrm{\Delta
x}^{4}),$ (3.19) $\displaystyle\mathrm{I}[Q^{(2)}_{\text{nw}}]$
$\displaystyle=B-2u_{xx}u_{xxx}\mathrm{\Delta
x}^{3}+2u_{yy}u_{yyy}\mathrm{\Delta x}^{3}+\mathcal{O}(\mathrm{\Delta
x}^{4}),$ $\displaystyle\mathrm{I}[Q^{(2)}_{\text{sw}}]$
$\displaystyle=B-2u_{xx}u_{xxx}\mathrm{\Delta
x}^{3}-2u_{yy}u_{yyy}\mathrm{\Delta x}^{3}+\mathcal{O}(\mathrm{\Delta
x}^{4}),$ $\displaystyle\mathrm{I}[Q^{(2)}_{\text{se}}]$
$\displaystyle=B+2u_{xx}u_{xxx}\mathrm{\Delta
x}^{3}-2u_{yy}u_{yyy}\mathrm{\Delta x}^{3}+\mathcal{O}(\mathrm{\Delta
x}^{4}),$ $\displaystyle\mathrm{I}[Q^{(3)}]$
$\displaystyle=B+\mathcal{O}(\mathrm{\Delta x}^{4})$

where $B=(u_{xx}^{2}+u_{xy}^{2}+u_{yy}^{2})\mathrm{\Delta x}^{2}$. Therefore,
the coefficients chosen for the definition of $\tau$ in (3.18), cancel all
terms up to $\mathcal{O}(\mathrm{\Delta x}^{3})$, thus ensuring that the
hypotheses of [CSV19, Theorem 24] hold.

##  4 Convergence

The convergence analysis will be carried out in one space dimension, and
follow the guidelines of [Fer02, CFR06], which will be briefly reviewed here.
Assume that equation (1.1) is posed on ${\mathbb{R}}$ in the simplified form

$\begin{cases}v_{t}(t,x)+H(Dv(t,x))=0&t,x\in(0,T)\times{\mathbb{R}},\\\
v(0,x)=v_{0}(x),&x\in{\mathbb{R}}.\end{cases}$ (4.1)

Also assume that the Hamiltonian function $H$ is a $W^{2,\infty}$ function and
that

$H^{\prime\prime}(p)\geq m_{H}>0.$ (4.2)

The convexity assumption (4.2) implies, by the Fenchel duality formula, that
$H(\cdot)$ can always be written in the form

$H(Dv(t,x))=\sup_{a\in\mathbb{R}}\\{a\,Dv(t,x)-H^{*}(a)\\},$

where $H^{*}$ denotes the Legendre transform of $H$. Moreover, the restriction
of the supremum can be limited to a suitable compact set $A$, so that problem
(4.1) can be recast in the case of a Hamiltonian of the form (1.3), by setting
$f_{D}(t,x,a)=-a$ and $f_{C}(t,x,a)=H^{*}(a)$.

The key assumption is that for any Lipschitz continuous function $v(x)$, once
defined the sequence $V=\\{v_{j}\\}_{j}=\\{v(x_{j})\\}_{j}$, the interpolation
operator $R[V]$ satisfies $R[V](x_{k})=v(x_{k})$ and, for some constant $C<1$:

$|R[V](x)-R_{1}[V](x)|\leq C\max_{x_{k}\in U(x)}|v_{k+1}-2v_{k}+v_{k-1}|$
(4.3)

where by $R_{1}[V]$ we denote the $\mathbb{P}_{1}$ (i.e., piecewise linear)
interpolation on the sequence $V$, and by $U(x)=(x-h_{-}\mathrm{\Delta
x},x+h_{+}\mathrm{\Delta x})$ the interval containing the stencil of the
reconstruction $R[V](x)$. For example, in our case the reconstruction is
performed taking two nodes on the left and two nodes on the right of the point
$x$, so that $h_{-}=h_{+}=2$.

Under such assumptions it is possible (see [Fer02]) to prove the following

###### Theorem 4.1

Consider the scheme (2.1) applied to equation (4.1), and assume that (4.2) and
(4.3) hold, that $\mathrm{\Delta x}=\mathcal{O}(\mathrm{\Delta t}^{2})$ and
that $v_{0}$ is Lipschitz continuous. Then, the numerical solution
$U^{n}=\\{u_{j}^{n}\\}_{j}$ (with $u_{j}^{n}$ defined by (2.1)) satisfies

$\|R[U^{n}]-v(n\mathrm{\Delta t})\|_{\infty}\to 0$

(where $v$ is the solution of (4.1)) for $0\leq n\leq T/\mathrm{\Delta t}$, as
$\mathrm{\Delta t}\to 0$.

It is proved in [Fer02] that condition (4.3) is satisfied for Lagrange
reconstructions up to the fifth order, if the reconstruction stencil includes
the interval $[x_{j},x_{j+1}]$. Since $\mathsf{WENO}$ interpolation is
performed by taking a convex combination of polynomials which satisfy (4.3) up
the degree 5 for the partial polynomials (and therefore, up to the degree 9
for the global interpolant), this fact is used in [CFR06] to obtain
convergence via Theorem 4.1. This essentially boils down to proving that the
all the function appearing in the convex combination defining $R$ satisfy
(4.3), and we plan to apply the same principle here.

In the case of $\mathsf{CWENO}$, once set $d_{L}=d_{R}=d$ and recast $P_{0}$
as

$P_{0}(x)=\frac{1}{1-2d}(Q(x)-dP_{L}(x)-dP_{R}(x)),$

we estimate the left-hand side of (4.3) as

$\displaystyle|R-R_{1}|$ $\displaystyle=$
$\displaystyle|\omega_{0}P_{0}+\omega_{L}P_{L}+\omega_{R}P_{R}-R_{1}|$
$\displaystyle=$
$\displaystyle|\omega_{0}(P_{0}-R_{1})+\omega_{L}(P_{L}-R_{1})+\omega_{R}(P_{R}-R_{1})|$
$\displaystyle\leq$
$\displaystyle\max\\{|P_{0}-R_{1}|,|P_{L}-R_{1}|,|P_{R}-R_{1}|\\},$

where the identity $\omega_{0}+\omega_{L}+\omega_{R}=1$ has been used twice.
As it has been proved in [Fer02], both $P_{L}$ and $P_{R}$ satisfy (4.3). It
remains then to check that the same property holds for $P_{0}$. We have:

$\displaystyle|P_{0}-R_{1}|$ $\displaystyle=$
$\displaystyle\left|\frac{1}{1-2d}(Q-dP_{L}-dP_{R})-R_{1}\right|$
$\displaystyle=$
$\displaystyle\left|\frac{1}{1-2d}(Q-dP_{L}-dP_{R}+(2d-1)R_{1})\right|$
$\displaystyle=$
$\displaystyle\left|\frac{1}{1-2d}\big{(}(Q-R_{1})-d(P_{L}-R_{1})-d(P_{R}-R_{1})\big{)}\right|$
$\displaystyle\leq$
$\displaystyle\frac{1}{1-2d}|Q-R_{1}|+\frac{d}{1-2d}|P_{L}-R_{1}|+\frac{d}{1-2d}|P_{R}-R_{1}|.$

Denoting now by $C_{r}$ the constant appearing in (4.3) for the specific case
of an interpolation of degree $r$, we finally obtain an estimate in the form

$|P_{0}-R_{1}|\leq\left(\frac{1}{1-2d}C_{3}+\frac{2d}{1-2d}C_{2}\right)\max_{x_{k}\in
U(x)}|v_{k+1}-2v_{k}+v_{k-1}|.$

As proved in [Fer02], $C_{2}=1/8$ and $C_{3}\approx 0.2533$; then, it turns
out that (4.3) is satisfied for $d\lesssim 0.332$, and therefore, with this
additional condition, all the assumptions of Theorem 4.1 are satisfied, and
the scheme (2.1) converges to the viscosity solution of (4.1).

We point out that our choice of $d_{0}=3/4$, i.e. $d=d_{L}=d_{R}=1/8$ in the
$\mathsf{CWENO}$ and $\mathsf{CWENOZ}$ reconstructions employed in this paper
fully satisfies the above requirement.

##  5 Numerical tests

Here, we present some numerical tests in order to assess the performance of
the schemes proposed in this paper. In particular, we focus on the expected
accuracy of the schemes, considering the $L^{1}$-norm of the error, and the
computational times, reported in seconds. In particular we will compare the SL
schemes using the $\mathsf{CWENO}$ and $\mathsf{CWENOZ}$ reconstructions of §3
and those employing the dimensionally-split $\mathsf{WENO}$ approach of
[CFR06].

In order to account for bounded domains, we consider either periodic boundary
conditions or extrapolation technique. It is noteworthy that in all the
numerical tests where extrapolation is applied, the solution has a compact
support, or the boundary of the numerical domain is entirely of outflow type.
In both cases, these conditions ensure no loss in accuracy in the numerical
treatment of the boundary. In all tests, unless otherwise stated, we employ an
extrapolation technique. For more general boundary conditions, we refer to
[CCDS23] for approximations of Neumann-type boundary conditions and to
[BCCF21] for a second order accurate treatment of Dirichlet boundary
conditions.

Although the convergence analysis requires a relation $\mathrm{\Delta
x}=\mathcal{O}(\mathrm{\Delta t}^{2})$ between the time and space steps, in
practice, smaller time steps are also allowed. Heuristically, assuming exact
minimization, and supposing a third order time discretization is used, then
the consistency error is given by

${\cal{T}}(\mathrm{\Delta t},\mathrm{\Delta
x})=\mathcal{O}\left(\mathrm{\Delta t}^{3}+\frac{\mathrm{\Delta
x}^{q}}{\mathrm{\Delta t}}\right),$

where $q$ is the order of the interpolation in space. For smooth solutions,
the choice $\mathrm{\Delta t}=\mathcal{O}\left(\mathrm{\Delta x}^{q/4}\right)$
optimizes this error and leads to a numerical convergence rate of $3q/4$.

In particular, choosing $q=4$ (i.e., a cubic reconstruction), we obtain

${\cal{T}}(\mathrm{\Delta t},\mathrm{\Delta
x})=\mathcal{O}\left(\mathrm{\Delta t}^{3}+\frac{\mathrm{\Delta
x}^{4}}{\mathrm{\Delta t}}\right),$

which results in an overall order of $3$ for the choice $\mathrm{\Delta
t}=\mathcal{O}(\mathrm{\Delta x})$. In the particular case where the
Hamiltonian is independent of $(t,x)$, the characteristics solving problem
(1.2) are affine, and even a first order time discretization computes them
exactly, reducing the consistency error to

${\cal{T}}(\mathrm{\Delta t},\mathrm{\Delta
x})=\mathcal{O}\left(\frac{\mathrm{\Delta x}^{q}}{\mathrm{\Delta t}}\right).$

In this case, the larger the time step, the smaller the consistency errors.

In what follows $N$ indicates the number of grid points per direction, so that
$\mathrm{\Delta x}\sim 1/N$. The tests have been performed on the cluster
Galileo 100 hosted at
CINECA111https://www.hpc.cineca.it/systems/hardware/galileo100/.

#### Test 1: Passive advection.

Figure 2: Initial condition (left) and numerical solution at $T=1$ (right) for Test 1 computed with $\mathsf{CWENOZ}$ reconstruction and $N=81$. Figure 3: Reconstructions count for Test 1 in the final time step on a grid of $81\times 81$ points. | $\mathsf{WENO}$ | $\mathsf{CWENO}$ | $\mathsf{CWENOZ}$
---|---|---|---
$N$ |  $L^{1}$-err |  $L^{1}$-ord |  $L^{1}$-err |  $L^{1}$-ord |  $L^{1}$-err |  $L^{1}$-ord
$21\times 21$ | $2.97e-03$ |  | $2.98e-03$ |  | $2.98e-03$ |
$41\times 41$ | $5.46e-04$ | $2.44$ | $4.88e-04$ | $2.61$ | $4.84e-04$ | $2.62$
$81\times 81$ | $1.24e-04$ | $2.13$ | $8.18e-05$ | $2.57$ | $7.85e-05$ | $2.62$
$161\times 161$ | $2.34e-05$ | $2.41$ | $1.22e-05$ | $2.75$ | $1.12e-05$ | $2.80$
$321\times 321$ | $3.92e-06$ | $2.58$ | $1.72e-06$ | $2.82$ | $1.49e-06$ | $2.92$
$641\times 641$ | $6.11e-07$ | $2.68$ | $2.27e-07$ | $2.92$ | $1.79e-07$ | $3.05$
$1281\times 1281$ | $8.87e-08$ | $2.78$ | $2.92e-08$ | $2.96$ | $2.06e-08$ | $3.12$
Table 1:  Errors at time $T=1$ for Test 1, $\mathsf{WENO}$, $\mathsf{CWENO}$ and $\mathsf{CWENOZ}$ schemes.  | $\mathsf{WENO}$ | $\mathsf{CWENO}$ | $\mathsf{CWENOZ}$
---|---|---|---
grid | CPU time | CPU time |  $\%$ gain | CPU time |  $\%$ gain
$21\times 21$ | $1.05e+00$ | $7.22e-01$ | $31.21$ | $1.01e-01$ | $90.36$
$41\times 41$ | $1.10e+00$ | $1.06e+00$ | $3.39$ | $7.20e-01$ | $34.25$
$81\times 81$ | $7.84e+00$ | $5.42e+00$ | $30.87$ | $5.46e+00$ | $30.39$
$161\times 161$ | $6.19e+01$ | $4.36e+01$ | $29.51$ | $4.37e+01$ | $29.35$
$321\times 321$ | $4.86e+02$ | $3.38e+02$ | $30.47$ | $3.36e+02$ | $30.95$
$641\times 641$ | $3.87e+03$ | $2.70e+03$ | $30.31$ | $2.67e+03$ | $30.99$
$1281\times 1281$ | $3.09e+04$ | $2.12e+04$ | $31.25$ | $2.12e+04$ | $31.32$
Table 2: CPU times for Test 1, $\mathsf{WENO}$ and $\mathsf{CWENO}$ schemes.

As a first benchmark problem, we consider the uniform rotation of a scalar
field in two space dimensions, expressed by the linear transport equation

$v_{t}-f_{D}\cdot Dv=0,$

which is included in our setting by choosing in (1.3) $A=\emptyset$ and
$f_{C}=0$. The speed of propagation is given by the vector field
$f_{D}=(-2\pi(x_{2}-0.5),2\pi(x_{1}-0.5))$ in the square domain
$[0,1]\times[0,1]$. We choose as initial condition the globally
$\mathcal{C}^{2}$ function

$v(0,x)=v_{0}(r)=M\left(1+\frac{r^{3}}{R^{3}}\left(-1+3\,\frac{r-R}{R}\left(1-2\,\frac{r-R}{R}\right)\right)\right)$

with $r=|x-(0.3,0.7)|$ and parameters $M=0.15$, $R=0.15$. Numerical solutions
are computed at the final time $T=1$, using a relationship of $\mathrm{\Delta
t}=3\mathrm{\Delta x}$ among the discretization steps. Since the advection is
rigid, the exact solution at final time coincides with the initial condition.
We point out that the characteristic lines are not affine, therefore we use
the third order RK scheme and a third order reconstruction scheme as described
in section §2 and §3, respectively. Plots of the initial data and of the
numerical solution are shown in Fig. 2.

As expected, since the transport is rigid and the initial data is globally
$\mathcal{C}^{2}$, all the schemes ($\mathsf{WENO}$, $\mathsf{CWENO}$ and
$\mathsf{CWENOZ}$) yield solutions converging with the same order, see Table
1. However, central reconstructions are more accurate with respect to the
traditional $\mathsf{WENO}$ one; in particular, $\mathsf{CWENOZ}$ achieves the
best results: its errors are about one-half than those of $\mathsf{WENO}$. The
fundamental difference among the schemes is in their computational costs,
which are shown in Table 2. The high number of reconstructions performed in
each cell (see Fig. 3) allows the $\mathsf{CWENO}$ and $\mathsf{CWENOZ}$ to
save about the $30\%$ of time. This happens because Central $\mathsf{WENO}$
reconstructions compute their nonlinear coefficients only once per cell
instead of once per reconstruction point. Looking at Fig. 3, note in
particular that the high number of reconstructions performed near the boundary
is related to the boundedness of the domain: in fact, except for the case of
periodic boundary conditions, each time a foot of a characteristic falls
outside the domain, it is projected back in the first cell near the boundary
in order to compute the reconstruction.

#### Test 2: One–dimensional semi-concave data.

$-2$$-1$$0$$1$$2$$-1$$-0.5$$0$exactnumerical

$-2$$-1$$0$$1$$2$$0$$100$$200$$300$

Figure 4: Left: the exact solution of Test 2 in red and the numerical one in blue computed with a third order $\mathsf{CWENOZ}$ reconstruction and $N=81$. Right: reconstructions count in the last time step. | $\mathsf{WENO}$ | $\mathsf{CWENO}$ | $\mathsf{CWENOZ}$
---|---|---|---
$N$ |  $L^{1}$-err |  $L^{1}$-ord |  $L^{1}$-err |  $L^{1}$-ord |  $L^{1}$-err |  $L^{1}$-ord
$81$ | $3.56e-06$ |  | $2.24e-06$ |  | $1.78e-06$ |
$161$ | $2.83e-07$ | $3.65$ | $1.80e-07$ | $3.64$ | $1.44e-07$ | $3.63$
$321$ | $2.45e-08$ | $3.53$ | $1.59e-08$ | $3.50$ | $1.30e-08$ | $3.47$
$641$ | $2.61e-09$ | $3.23$ | $1.96e-09$ | $3.02$ | $1.75e-09$ | $2.90$
Table 3:  Errors at time $T=1$ for Test 2, $\mathsf{WENO}$, $\mathsf{CWENO}$ and $\mathsf{CWENOZ}$ schemes.  | $\mathsf{WENO}$ | $\mathsf{CWENO}$ | $\mathsf{CWENOZ}$
---|---|---|---
grid | CPU time | CPU time |  $\%$ gain | CPU time |  $\%$ gain
$81$ | $6.19e-03$ | $4.83e-03$ | $21.95$ | $4.86e-03$ | $21.48$
$161$ | $9.79e-03$ | $8.88e-03$ | $9.31$ | $8.64e-03$ | $11.76$
$321$ | $2.74e-02$ | $2.39e-02$ | $13.03$ | $2.39e-02$ | $12.84$
$641$ | $9.42e-02$ | $7.65e-02$ | $18.81$ | $7.68e-02$ | $18.41$
$1281$ | $3.48e-01$ | $2.83e-01$ | $18.65$ | $2.90e-01$ | $16.72$
$2561$ | $1.33e+00$ | $1.10e+00$ | $17.15$ | $1.10e+00$ | $16.83$
$5121$ | $5.28e+00$ | $4.30e+00$ | $18.51$ | $4.29e+00$ | $18.77$
$10241$ | $2.09e+01$ | $1.68e+01$ | $19.42$ | $1.68e+01$ | $19.46$
Table 4: CPU times for Test 2, $\mathsf{WENO}$ and $\mathsf{CWENO}$ schemes.

This test deals with the HJ equation

$\begin{cases}v_{t}(t,x)+\frac{1}{2}|\partial_{x}v(t,x)|^{2}=0\\\
v(0,x)=v_{0}(x)=\min(-\cos(\pi x/2),0),\end{cases}$ (5.1)

in the domain $[-2,2]$, with homogeneous boundary conditions, which are
treated using the extrapolation technique. The approximate solution is
computed, at $T=1$, with $\mathrm{\Delta t}=10\mathrm{\Delta x}$, while the
exact solution $v(t,x)$ is defined as follows:

$v(t,x)=\min(0,\frac{1}{2}[a^{*}(t,x)]^{2}+v_{0}(x-ta^{*}(t,x))),$ (5.2)

where, for a given $(t,x)$, $a^{*}(t,x)$ represents the optimal control, which
is constant along the characteristics, and it can be computed as the solution
of

$a^{*}(t,x)=g(t,x,a^{*}(t,x))$ (5.3)

with the function $g$ defined by

$g(t,x,a^{*}(t,x))=\begin{cases}-\frac{2}{\pi
t}\arcsin{\frac{2a^{*}(t,x)}{\pi}}+\frac{x}{t}\quad&\text{if}\;|a^{*}(t,x)|\leq\frac{\pi}{2},\\\
\frac{\pi}{2}\quad&\text{if}\;a^{*}(t,x)>\frac{\pi}{2},\\\
-\frac{\pi}{2}\quad&\text{if}\;a^{*}(t,x)<-\frac{\pi}{2}.\\\ \end{cases}$
(5.4)

For our purposes, since the exact solution of (5.3) is not available, we
compute an approximate solution by using a fixed point algorithm with initial
guess $a^{*}_{0}(t,x)=\frac{1}{2}x^{2}t$.

Since the characteristics are affine, in this test we have combined the third
order spatial reconstructions with RK1 time stepping, which computes these
curves exactly. Plots of the numerical and exact solutions are shown in the
left panel of Figure 4 and $L^{1}$ errors are listed in Table 3. In this case
too, all schemes compute approximations affected by errors of the same order
but the central reconstructions turn out to be more accurate (errors are about
$2/3$ than the $\mathsf{WENO}$ case). In fact, the feet of the characteristics
fall in regular regions of $v$ at each time step, and therefore all the
reconstructions are close to the optimal cubic one. Table 4 reports the
computational costs of each scheme, showing an average gain of $16-19\%$ for
both the Central $\mathsf{WENO}$ reconstructions. This is again the effect of
the multiple reconstructions performed in each cell (see the right panel of
Figure 4, in which the number of reconstructions performed in each cell during
the last time step is depicted). We point out that, as in Test 1, the high
number of reconstructions computed in the first and in the last cell is
related to the boundary conditions.

#### Test 3: One–dimensional eikonal equation.

$0$$1$$2$$3$$4$$5$$6$$-1$$0$$1$exactnumerical

$0$$1$$2$$3$$4$$5$$6$$100$$200$$300$

Figure 5: Left: the exact solution of Test 3 in red and the numerical one in blue computed with a third order $\mathsf{CWENOZ}$ scheme and $N=63$. Right: reconstructions count in the last time step. | $\mathsf{WENO}$ | $\mathsf{CWENO}$ | $\mathsf{CWENOZ}$
---|---|---|---
$N$ |  $L^{1}$-err |  $L^{1}$-ord |  $L^{1}$-err |  $L^{1}$-ord |  $L^{1}$-err |  $L^{1}$-ord
$126$ | $1.28e-05$ |  | $1.38e-05$ |  | $1.43e-05$ |
$252$ | $1.54e-06$ | $3.05$ | $1.65e-06$ | $3.07$ | $1.70e-06$ | $3.07$
$503$ | $1.93e-07$ | $3.00$ | $2.03e-07$ | $3.02$ | $2.10e-07$ | $3.02$
$1006$ | $2.56e-08$ | $2.91$ | $2.63e-08$ | $2.95$ | $2.71e-08$ | $2.95$
Table 5:  Errors at time $T=1$ for Test 3, $\mathsf{WENO}$, $\mathsf{CWENO}$ and $\mathsf{CWENOZ}$ schemes.  | $\mathsf{WENO}$ | $\mathsf{CWENO}$ | $\mathsf{CWENOZ}$
---|---|---|---
grid | CPU time | CPU time |  $\%$ gain | CPU time |  $\%$ gain
$126$ | $5.92e-02$ | $5.42e-02$ | $8.43$ | $5.64e-02$ | $4.81$
$252$ | $2.11e-01$ | $1.94e-01$ | $8.04$ | $1.94e-01$ | $8.05$
$503$ | $7.97e-01$ | $7.37e-01$ | $7.53$ | $7.48e-01$ | $6.05$
$1006$ | $3.13e+00$ | $2.83e+00$ | $9.75$ | $2.87e+00$ | $8.45$
$2011$ | $1.19e+01$ | $1.11e+01$ | $6.67$ | $1.11e+01$ | $6.79$
$4022$ | $4.73e+01$ | $4.27e+01$ | $9.58$ | $4.21e+01$ | $11.03$
$8043$ | $1.80e+02$ | $1.65e+02$ | $8.63$ | $1.62e+02$ | $10.02$
$16085$ | $6.96e+02$ | $6.41e+02$ | $7.96$ | $6.38e+02$ | $8.35$
Table 6: CPU times for Test 3, $\mathsf{WENO}$ and $\mathsf{CWENO}$ schemes.

In this test, we consider the HJ equation

$\begin{cases}v_{t}(t,x)+\frac{1}{2}|\partial_{x}v(t,x)|^{2}-f(t,x)=0\\\
v(0,x)=v_{0}(x),\end{cases}$ (5.5)

with $f(t,x)=-\sin(x)+(\frac{9}{8}+\frac{t^{2}-3t}{2})\cos^{2}(x)$ and
$v_{0}(x)=\frac{3}{2}\sin(x)$. We consider the problem in $[0,2\pi]$ with
periodic boundary conditions and compute the numerical solution at time
$T=0.5$, with $\mathrm{\Delta t}=\mathrm{\Delta x}$. This problem has a known
exact solution given by $v(t,x)=(\frac{3}{2}-t)\sin(x)$.

In this case, characteristics are not affine and both the dynamic and the cost
function explicitly depend on $x$ and $t$. Therefore this is a specific test
for the third order accuracy of our scheme, which depends on the use of the
third order RK scheme (2.4) and on the approximation of the integral of the
cost function.

Plots of the exact and of the numerical solution are shown in the left panel
of Fig.5, while in the right panel the reconstruction count in the last time
step is depicted. We point out that in this case, with periodic boundary
conditions, there is no artificial increase in the first and in the last cell
near the boundary.

Errors and computational times are listed in Tables 5 and 6. As in the
previous tests, due to the regularity of the exact solution, all schemes
perform almost equally when compared in terms of the error. Looking at the
computational costs, the Central $\mathsf{WENO}$ reconstructions allow for a
saving of about $6-10\%$.

#### Test 4: Two–dimensional semi-convex data.

Figure 6: Test 4. Initial condition (left) and numerical solution computed
with $\mathsf{CWENO}$ at $T=1$ (right) with $N=81$.

& $0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.7$$0.8$$0.9$$1$exact$\mathsf{CWENO}$$\mathsf{CWENOZ}$
$0.4$$0.45$$0.5$$0.55$$0.6$$0.2$$0.3$exact$\mathsf{CWENO}$$\mathsf{WENO}$
$0.4$$0.45$$0.5$$0.55$$0.6$$0.2$$0.3$exact$\mathsf{CWENOZ}$$\mathsf{WENO}$
$0.9$$0.95$$1$$1.05$$1.1$$0$$0.5$$1$$\cdot
10^{-2}$exact$\mathsf{CWENO}$$\mathsf{WENO}$
$0.9$$0.95$$1$$1.05$$1.1$$0$$0.5$$1$$\cdot
10^{-2}$exact$\mathsf{CWENOZ}$$\mathsf{WENO}$

Figure 7: Test 4. Comparison of overshoots and undershoots of the numerical
solutions computed with the traditional $\mathsf{WENO}$ scheme and with the
central $\mathsf{WENO}$ ones on a $161\times 161$ grid. In the left panels
$\mathsf{WENO}$ and $\mathsf{CWENO}$ are compared, while in the right panels
similar plots are shown for $\mathsf{WENO}$ and $\mathsf{CWENOZ}$. The exact
solution is always represented with a thick black line. First row: scatter
plots of all data points with $r\in[0,0.2]$ corresponding to the spike of the
solution. Second row: scatter plots of all data points with $r\in[0.4,0.6]$
corresponding to regular region of the solution. Third row: scatter plots of
all data points with $r\in[0.9,1.1]$ corresponding to the singular region of
the initial data.

| $\mathsf{WENO}$ | $\mathsf{CWENO}$ | $\mathsf{CWENOZ}$
---|---|---|---
$N$ |  $L^{1}$-err |  $L^{1}$-ord |  $L^{1}$-err |  $L^{1}$-ord |  $L^{1}$-err |  $L^{1}$-ord
$41\times 41$ | $4.02e-02$ |  | $3.38e-02$ |  | $3.38e-02$ |
$81\times 81$ | $2.25e-02$ | $0.84$ | $1.82e-02$ | $0.90$ | $1.81e-02$ | $0.90$
$161\times 161$ | $1.15e-02$ | $0.96$ | $9.01e-03$ | $1.01$ | $8.99e-03$ | $1.01$
$321\times 321$ | $5.56e-03$ | $1.05$ | $4.10e-03$ | $1.14$ | $4.09e-03$ | $1.14$
$641\times 641$ | $2.85e-03$ | $0.97$ | $2.03e-03$ | $1.02$ | $2.02e-03$ | $1.02$
Table 7:  Errors at time $T=1$ for Test 4, $\mathsf{WENO}$, $\mathsf{CWENO}$
and $\mathsf{CWENOZ}$ schemes.

| $\mathsf{WENO}$ | $\mathsf{CWENO}$ | $\mathsf{CWENOZ}$
---|---|---|---
grid | CPU time | CPU time |  $\%$ gain | CPU time |  $\%$ gain
$41\times 41$ | $2.88e+00$ | $1.98e+00$ | $31.37$ | $1.96e+00$ | $32.07$
$81\times 81$ | $2.19e+01$ | $1.54e+01$ | $29.70$ | $1.53e+01$ | $30.07$
$161\times 161$ | $1.72e+02$ | $1.19e+02$ | $31.10$ | $1.18e+02$ | $31.35$
$321\times 321$ | $1.37e+03$ | $9.38e+02$ | $31.56$ | $9.34e+02$ | $31.86$
$641\times 641$ | $1.09e+04$ | $7.44e+03$ | $31.59$ | $7.51e+03$ | $30.93$
Table 8: CPU times for Test 4, $\mathsf{WENO}$ and $\mathsf{CWENO}$ schemes.

| $\mathsf{WENO}$ | $\mathsf{CWENO}$ | $\mathsf{CWENOZ}$
---|---|---|---
$N$ | min | max | min | max | min | max
$21\times 21$ | $-5.67e-05$ | $9.76e-01$ | $-6.08e-03$ | $9.80e-01$ | $-6.08e-03$ | $9.80e-01$
$41\times 41$ | $-4.62e-05$ | $9.89e-01$ | $-5.19e-03$ | $9.92e-01$ | $-5.19e-03$ | $9.92e-01$
$81\times 81$ | $-1.21e-05$ | $9.95e-01$ | $-3.39e-03$ | $9.96e-01$ | $-3.39e-03$ | $9.96e-01$
$161\times 161$ | $-3.45e-06$ | $9.98e-01$ | $-1.97e-03$ | $9.98e-01$ | $-1.97e-03$ | $9.98e-01$
Table 9: Overshoots and undershoots observed in Test 4, $\mathsf{WENO}$,
$\mathsf{CWENO}$ and $\mathsf{CWENOZ}$ schemes.

In this test, the HJ equation

$\begin{cases}v_{t}(t,x)+\frac{1}{2}|Dv(t,x)|^{2}=0\\\
v(0,x)=v_{0}(x)=\max(1-|x|^{2},0),\end{cases}$ (5.6)

is considered in $[-2,2]^{2}$, using extrapolation techniques to treat the
homogeneous boundary condition. The exact solution of this problem is known
and is given, for $t\geq 1/2$, by

$v(t,x)=\begin{cases}\frac{(|x|-1)^{2}}{2t}\quad&\text{if}\;|x|\leq 1,\\\
0\quad&\text{if}\;|x|>1.\end{cases}$ (5.7)

We set final time $T=0.5$ and compute the approximate solution with
$\mathrm{\Delta t}=\frac{5}{4}\mathrm{\Delta x}$. Initial condition and final
numerical solution are shown in Figure 6, while $L^{1}$ errors and CPU times
are reported respectively in Tables 7 and 8.

Since in this case the feet of characteristics always fall in the singular
region of the solution, we observe, as expected, that all numerical schemes
are degraded to first order accuracy. From Table 7 one can see that the
$L^{1}$-norm of the error is about $30\%$ lower for the central
reconstructions, as in the previous cases. From Table 9 and Figure 5 we can
see that near the kink in the solution, the central reconstructions produce
larger undershoots in the flat region, while being closer to the exact
solution in the nonzero region. Finally, from Table 8 one can observe that
central reconstructions schemes are about $30\%$ faster.

#### Test 5: Front propagation with obstacles.

Figure 8: Reachable sets (5.9) computed at each time step for Test 5 on a grid $101\times 81$ and $T=3$, using $\mathsf{WENO}$ (left) and $\mathsf{CWENO}$ (right) reconstructions. In both panel the black line represents a reference solution computed on a grid $1001\times 801$, the red circle represents the target and the red rectangles represent the obstacles. | $\mathsf{WENO}$ | $\mathsf{CWENO}$ | $\mathsf{CWENOZ}$
---|---|---|---
grid | CPU time | CPU time |  $\%$ gain | CPU time |  $\%$ gain
$126\times 126$ | $6.28e+00$ | $5.34e+00$ | $14.94$ | $5.30e+00$ | $15.66$
$251\times 251$ | $4.56e+01$ | $3.96e+01$ | $13.04$ | $3.96e+01$ | $13.06$
$501\times 501$ | $3.61e+02$ | $3.11e+02$ | $13.87$ | $3.20e+02$ | $11.50$
$1001\times 1001$ | $2.81e+03$ | $2.43e+03$ | $13.64$ | $2.40e+03$ | $14.59$
$2001\times 2001$ | $2.20e+04$ | $1.92e+04$ | $13.01$ | $1.89e+04$ | $14.00$
Table 10: CPU times for Test 5, $\mathsf{WENO}$ and $\mathsf{CWENO}$ schemes.

Last, we consider a state constrained Zermelo problem, proposed in [BFZ10],
for a swimmer that has to reach a circular island at $(1,0)$, facing a non-
constant current flowing in the horizontal direction. Two obstacles are placed
in the domain, near the target.

More precisely, to consider the minimum time problem the running cost
$f_{C}(t,x,a)$ is set to zero, and the nonlinear dynamics is given by
$f_{D}(t,x,a)=(2-0.5x_{2}^{2}+a_{1},a_{2})$, with
$A:=\\{a\in\mathbb{R}^{2}\,|\,\|a\|=(a_{1}^{2}+a^{2}_{2})^{1/2}=1\\}$. The
swimmer has unit speed and can move in any direction $a\in A$; the current has
speed $2-0.5x_{2}^{2}$ and flows in the direction $(1,0)$. The target is the
set $\mathcal{X}=\\{x\;|\;v_{0}(x)\leq 0\\}$, where
$v_{0}(x)=C\,{\rm{min}}\big{(}\|(x_{1}-1,x_{2})\|-r,r\big{)}$ with $r=0.25$
and $C=20$. In order to take into account the presence of the two obstacles,
we follow the level set approach, as proposed in [BFZ10], and represent the
set of constrained states as $\mathcal{K}=\\{x\in\mathbb{R}^{2}\;|\;g(x)\leq
0\\}$ where g(x)=max(-γ, C( γ-max (—x_1-0.3—,—x_2-0.4—)), C(
γ-max(—x_1+1—,—x_2+1.5—) ) ), with $C=20$, $\gamma=0.2$.

For every $t\in[0,T]$, the set of points from which the swimmer can reach the
target before time $t$ (Backward Reachable set) is represented by the level
set
$\mathcal{R}[0,t]=\\{x\in\mathbb{R}^{2}\,|\,\exists\,s\in[0,t],\,v(x,s)\leq
0\\}$ where $v$ is the solution to a constrained version of (1.1)

$\min(v_{t}+H(t,x,Dv),v-g(x))=0$ (5.8)

with initial condition $v(x,0)={\rm{max}}(v_{0}(x),g(x))$, (see [BFZ10, Remark
2]). The space domain is defined as $[-3,2]\times[-2,2]$ and the final time to
$T=3$. By using scheme (2.1), we discetize (5.8) as min(u^n+1_i-a
∈A^νmin{R[u^n](y^n_i(a))} ,u^n+1-g(x_i))=0, resulting in the following time
marching scheme

${u}_{i}^{n+1}=\max\left(\underset{\underline{a}\in
A^{\nu}}{\min}\\{R[u^{n}](y^{n}_{i}(\underline{a}))\\},g(x_{i})\right),$

from which we define the numerical approximation of the seachable set
$\mathcal{R}[0,t_{n}]$ as

$\mathcal{R}^{n}=\\{x_{i}\,|\,\exists\,k\in\\{0,\dots,n\\},u^{k}_{i}\leq
0\\}.$ (5.9)

The solution is computed with RK3 and $\mathrm{\Delta t}=\mathrm{\Delta x}$.
In Figure 8, we show that, even in this very complex situation, no relevant
difference in the accuracy of the computed solution can be observed between
the $\mathsf{WENO}$ and the Central $\mathsf{WENO}$ approaches. On the other
hand, Central reconstruction schemes provide a reduction of about $13-14\%$ of
the computational time.

## Conclusion

In this study, a high order semi-Lagrangian scheme has been developed to
approximate first order Hamilton–Jacobi–Bellman equations. The primary
challenge lies in the nonsmooth nature of the solutions, which would lead to
spurious oscillations when using unlimited high order polynomial
interpolation.

To address this issue, the reconstruction of the solution at the previous time
step at the foot of characteristics has to be performed with a non-oscillatory
scheme. Since one has to perform very many interpolations during the
minimization procedure, a $\mathsf{CWENO}$ technique was proposed to be
coupled with the semi-Lagrangian approach. Our study demonstrates that, in
terms of errors, our new scheme maintains the favorable behavior of
$\mathsf{WENO}$ schemes of [CFR06], producing about $30\%$ more accurate
results in smooth regions at the price of some extra over/undershoots.
However, its computational cost is significantly lower, by a $10-30\%$
depending on the specific test.

Additionally, we have established a convergence result in a simpler case, and
provided several numerical simulations to further validate the effectiveness
of our proposed scheme.

### Future Directions

Future research directions could explore the inclusion of high order
treatments for boundary conditions and extensions of the convergence analysis
to more general Hamiltonians. Investigating the scheme’s performance in more
complex scenarios would also contribute to a comprehensive understanding of
its capabilities.

### Acknowledgments

The first two authors would like to thank the Italian Ministry of Instruction,
University and Research (MIUR) for supporting this research with funds coming
from the PRIN Project 2022 (2022238YY5, entitled ”Optimal control problems:
analysis, approximation”). All authors are members of the GNCS-INdAM Gruppo
Nazionale per il Calcolo Scientifico of Istituto Nazionale di Alta Matematica,
Francesco Severi, P.le Aldo Moro, Roma, Italy.

## References

  * [Abg94] R. Abgrall.  On essentially non-oscillatory schemes on unstructured meshes: analysis and implementation.  J. Comput. Phys., 114(6):45–48, 1994.
  * [BBMZ19] A. Baeza, R. Bürger, P. Mulet, and D. Zorío.  Central WENO schemes through a global average weight.  J. Sci. Comput., 78(1):499–530, 2019.
  * [BCCF21] L. Bonaventura, E. Calzola, E. Carlini, and R. Ferretti.  Second order fully semi-Lagrangian discretizations of advection-diffusion-reaction systems.  J. Sci. Comput., 88(23), 2021.
  * [BCD97] M. Bardi and I. Capuzzo-Dolcetta.  Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations.  Systems & Control: Foundations & Applications. Birkhäuser Boston, Inc., Boston, MA, 1997.  With appendices by Maurizio Falcone and Pierpaolo Soravia.
  * [BFS16] O. Bokanowski, M. Falcone, and S. Sahu.  An efficient filtered scheme for some first order time-dependent Hamilton-Jacobi equations.  SIAM J. Sci. Comput., 38(1):A171–A195, 2016.
  * [BFZ10] O. Bokanowski, N. Forcadel, and H. Zidani.  Reachability and minimal times for state constrained nonlinear problems without any controllability assumption.  SIAM J. Control Optim., 48(7):4292–4316, 2010.
  * [BGFB20] D. S. Balsara, S. Garain, V. Florinski, and W. Boscheri.  An efficient class of WENO schemes with adaptive order for unstructured meshes.  J. Comput. Phys., 404:109062, 2020.
  * [CCDS23] E. Calzola, E. Carlini, X. Dupuis, and F. J. Silva.  A semi-Lagrangian scheme for Hamilton-Jacobi-Bellman equations with oblique derivatives boundary conditions.  Numer. Math., 153(1):49–84, 2023.
  * [CFR06] E. Carlini, R. Ferretti, and G. Russo.  A weighted essentially nonoscillatory, large time-step scheme for Hamilton-Jacobi equations.  SIAM J. Sci. Comput., 27(3):1071–1091, 2006.
  * [CPSV18] I. Cravero, G. Puppo, M. Semplice, and G. Visconti.  CWENO: uniformly accurate reconstructions for balance laws.  Math. Comp., 87(312):1689–1719, 2018.
  * [CS19] M. J. Castro and M. Semplice.  Third- and fourth-order well-balanced schemes for the shallow water equations based on the CWENO reconstruction.  Int. J. Numer. Meth. Fluid, 89(8):304–325, 2019.
  * [CSV19] I. Cravero, M. Semplice, and G. Visconti.  Optimal definition of the nonlinear weights in multidimensional Central WENOZ reconstructions.  SIAM J. Numer. Anal., 57(5):2328–2358, 2019.
  * [DBSR17] M. Dumbser, W. Boscheri, M. Semplice, and G. Russo.  Central weighted ENO schemes for hyperbolic conservation laws on fixed and moving unstructured meshes.  SIAM J. Sci. Comput., 39(6):A2564–A2591, 2017.
  * [Fer02] R. Ferretti.  Convergence of semi-Lagrangian approximations to convex Hamilton-Jacobi equations under (very) large Courant numbers.  SIAM J. Numer. Anal., 40(6):2240–2253, 2002.
  * [FF94] M. Falcone. and R. Ferretti.  Discrete time high-order schemes for viscosity solutions of Hamilton-Jacobi-Bellman equations.  Numer. Math., 67(3):315–344, 1994.
  * [FF02] M. Falcone. and R. Ferretti.  Semi-Lagrangian schemes for Hamilton-Jacobi equations, discrete representation formulae and Godunov methods.  J. Comput. Phys., 175(2):559–575, 2002.
  * [FF14] M. Falcone. and R. Ferretti.  Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2014.
  * [FPT20a] M. Falcone, G. Paolucci, and S. Tozza.  Convergence of adaptive filtered schemes for first order evolutionary Hamilton-Jacobi equations.  Numer. Math., 145(2):271–311, 2020.
  * [FPT20b] M. Falcone, G. Paolucci, and S. Tozza.  Multidimensional smoothness indicators for first-order Hamilton-Jacobi equations.  J. Comput. Phys., 409, 2020.
  * [HEOC87] A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy.  Uniformly high-order accurate essentially nonoscillatory schemes III.  J. Comput. Phys., 71(2):231–303, 1987.
  * [HS99] C. Hu and C. W. Shu.  Weighted essentially non-oscillatory schemes on triangular meshes.  J. Comput. Phys., 150(1):97–127, 1999.
  * [HT03] T. Y. Hou and E. Tadmor, editors.  High-Order Schemes for Multi-Dimensional Hamilton-Jacobi Equations, Berlin, Heidelberg, 2003. Springer Berlin Heidelberg.
  * [JP00] G.-S. Jiang and D. Peng.  Weighted ENO schemes for Hamilton-Jacobi equations.  SIAM J. Sci. Comput., 21(6):2126–2143, 2000.
  * [JS96] G. S. Jiang and C. W. Shu.  Efficient implementation of weighted ENO schemes.  J. Comput. Phys., 126:202–228, 1996.
  * [LOC94] X. D. Liu, S. Osher, and T. Chan.  Weighted essentially non-oscillatory schemes.  J. Comput. Phys., 115(1):200–212, 1994.
  * [LPR00] D. Levy, G. Puppo, and G. Russo.  Compact central WENO schemes for multidimensional conservation laws.  SIAM J. Sci. Comput., 22(2):656–672, 2000.
  * [LT00] C. T. Lin and E. Tadmor.  High-resolution nonoscillatory central schemes for Hamilton-Jacobi equations.  SIAM J. Sci. Comput., 21(6):2163–2186, 2000.
  * [OS91] S. Osher and C. W. Shu.  High–order essentially non–oscillatory schemes for Hamilton-Jacobi equation.  SIAM J. Numer. Anal., 28:907–922, 1991.
  * [SCR16] M. Semplice, A. Coco, and G. Russo.  Adaptive mesh refinement for hyperbolic systems based on third-order Compact WENO reconstruction.  J. Sci. Comput., 66:692–724, 2016.
  * [SV20] M. Semplice and G. Visconti.  Efficient implementation of adaptive order reconstructions.  J. Sci. Comput., 83(1), 2020.
  * [ZCZ08] J. Zhou, L. Cai, and F.-Q. Zhou.  New high-resolution scheme for three-dimensional nonlinear hyperbolic conservation laws.  Appl. Math. Comp., 198(2):770–786, 2008.
  * [Zen86] M. Zennaro.  Natural continuous extensions of Runge-Kutta methods.  Math. Comp., 46(173):119–133, 1986.
  * [ZQ16] J. Zhu and J. Qiu.  A new fifth order finite difference WENO scheme for solving hyperbolic conservation laws.  J. Comput. Phys., 318:110–121, 2016.
  * [ZQ17] J. Zhu and J. Qiu.  A new fifth order finite difference WENO scheme for Hamilton-Jacobi equations.  Numer. Meth. for PDEs, 33(4):1095–1113, 2017.
  * [ZQ20] J. Zhu and J. Qiu.  A new type of high-order WENO schemes for Hamilton-Jacobi equations on triangular meshes.  Comm. Computat. Phys., 27(3):897–920, 2020.
  * [ZSQ19] F. Zheng, C. W. Shu, and J. Qiu.  High order finite difference Hermite WENO schemes for the Hamilton–Jacobi equations on unstructured meshes.  Comp. & Fluids, 183:53–65, 2019.

## Appendix

We report here the matrix involved in the definition (3.16).

$M=\setcounter{MaxMatrixCols}{16}\begin{pmatrix}0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&4&0&0&6&2&0&0&3&\frac{4}{3}&0&2&1&\frac{3}{2}\\\
0&0&0&0&1&0&0&1&1&0&1&1&1&1&1&1\\\
0&0&0&0&0&4&0&0&2&6&0&\frac{4}{3}&3&1&2&\frac{3}{2}\\\
0&0&0&6&0&0&84&3&0&0&42&2&0&28&\frac{3}{2}&21\\\
0&0&0&2&1&0&3&\frac{32}{3}&1&0&\frac{31}{2}&\frac{31}{3}&1&15&\frac{152}{15}&\frac{147}{10}\\\
0&0&0&0&1&2&0&1&\frac{32}{3}&3&1&\frac{31}{3}&\frac{31}{2}&\frac{152}{15}&15&\frac{147}{10}\\\
0&0&0&0&0&6&0&0&3&84&0&2&42&\frac{3}{2}&28&21\\\
0&0&0&3&1&0&42&\frac{31}{2}&1&0&\frac{989}{5}&15&1&\frac{954}{5}&\frac{147}{10}&\frac{933}{5}\\\
0&0&0&\frac{4}{3}&1&\frac{4}{3}&2&\frac{31}{3}&\frac{31}{3}&2&15&\frac{3992}{45}&15&\frac{1918}{15}&\frac{1918}{15}&\frac{737}{4}\\\
0&0&0&0&1&3&0&1&\frac{31}{2}&42&1&15&\frac{989}{5}&\frac{147}{10}&\frac{954}{5}&\frac{933}{5}\\\
0&0&0&2&1&1&28&15&\frac{152}{15}&\frac{3}{2}&\frac{954}{5}&\frac{1918}{15}&\frac{147}{10}&\frac{56076}{35}&\frac{737}{4}&\frac{161501}{70}\\\
0&0&0&1&1&2&\frac{3}{2}&\frac{152}{15}&15&28&\frac{147}{10}&\frac{1918}{15}&\frac{954}{5}&\frac{737}{4}&\frac{56076}{35}&\frac{161501}{70}\\\
0&0&0&\frac{3}{2}&1&\frac{3}{2}&21&\frac{147}{10}&\frac{147}{10}&21&\frac{933}{5}&\frac{737}{4}&\frac{933}{5}&\frac{161501}{70}&\frac{161501}{70}&\frac{721401}{25}\\\
\end{pmatrix}.$

Also, we report the matrices involved in the definitions (3.17). We assume
that the vector $U$ of data in the $4\times 4$ stencil is ordered
lexicographically, i.e. $x$ direction faster and $y$ direction slower.

$A_{\text{ne}}=\setcounter{MaxMatrixCols}{16}\begin{pmatrix}0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&\frac{3497}{720}&\frac{-623}{90}&\frac{1787}{720}&0&\frac{-623}{90}&\frac{1667}{180}&\frac{-571}{180}&0&\frac{1787}{720}&\frac{-571}{180}&\frac{797}{720}\\\
0&0&0&0&0&\frac{-623}{90}&\frac{1141}{90}&\frac{-229}{45}&0&\frac{1667}{180}&\frac{-1547}{90}&\frac{1187}{180}&0&\frac{-571}{180}&\frac{278}{45}&\frac{-421}{180}\\\
0&0&0&0&0&\frac{1787}{720}&\frac{-229}{45}&\frac{1817}{720}&0&\frac{-571}{180}&\frac{1187}{180}&\frac{-293}{90}&0&\frac{797}{720}&\frac{-421}{180}&\frac{827}{720}\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&\frac{-623}{90}&\frac{1667}{180}&\frac{-571}{180}&0&\frac{1141}{90}&\frac{-1547}{90}&\frac{278}{45}&0&\frac{-229}{45}&\frac{1187}{180}&\frac{-421}{180}\\\
0&0&0&0&0&\frac{1667}{180}&\frac{-1547}{90}&\frac{1187}{180}&0&\frac{-1547}{90}&\frac{1472}{45}&\frac{-1157}{90}&0&\frac{1187}{180}&\frac{-1157}{90}&\frac{887}{180}\\\
0&0&0&0&0&\frac{-571}{180}&\frac{1187}{180}&\frac{-293}{90}&0&\frac{278}{45}&\frac{-1157}{90}&\frac{571}{90}&0&\frac{-421}{180}&\frac{887}{180}&\frac{-109}{45}\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&\frac{1787}{720}&\frac{-571}{180}&\frac{797}{720}&0&\frac{-229}{45}&\frac{1187}{180}&\frac{-421}{180}&0&\frac{1817}{720}&\frac{-293}{90}&\frac{827}{720}\\\
0&0&0&0&0&\frac{-571}{180}&\frac{278}{45}&\frac{-421}{180}&0&\frac{1187}{180}&\frac{-1157}{90}&\frac{887}{180}&0&\frac{-293}{90}&\frac{571}{90}&\frac{-109}{45}\\\
0&0&0&0&0&\frac{797}{720}&\frac{-421}{180}&\frac{827}{720}&0&\frac{-421}{180}&\frac{887}{180}&\frac{-109}{45}&0&\frac{827}{720}&\frac{-109}{45}&\frac{857}{720}\\\
\end{pmatrix}$

$A_{\text{nw}}=\setcounter{MaxMatrixCols}{16}\begin{pmatrix}0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&\frac{1817}{720}&\frac{-229}{45}&\frac{1787}{720}&0&\frac{-293}{90}&\frac{1187}{180}&\frac{-571}{180}&0&\frac{827}{720}&\frac{-421}{180}&\frac{797}{720}&0\\\
0&0&0&0&\frac{-229}{45}&\frac{1141}{90}&\frac{-623}{90}&0&\frac{1187}{180}&\frac{-1547}{90}&\frac{1667}{180}&0&\frac{-421}{180}&\frac{278}{45}&\frac{-571}{180}&0\\\
0&0&0&0&\frac{1787}{720}&\frac{-623}{90}&\frac{3497}{720}&0&\frac{-571}{180}&\frac{1667}{180}&\frac{-623}{90}&0&\frac{797}{720}&\frac{-571}{180}&\frac{1787}{720}&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&\frac{-293}{90}&\frac{1187}{180}&\frac{-571}{180}&0&\frac{571}{90}&\frac{-1157}{90}&\frac{278}{45}&0&\frac{-109}{45}&\frac{887}{180}&\frac{-421}{180}&0\\\
0&0&0&0&\frac{1187}{180}&\frac{-1547}{90}&\frac{1667}{180}&0&\frac{-1157}{90}&\frac{1472}{45}&\frac{-1547}{90}&0&\frac{887}{180}&\frac{-1157}{90}&\frac{1187}{180}&0\\\
0&0&0&0&\frac{-571}{180}&\frac{1667}{180}&\frac{-623}{90}&0&\frac{278}{45}&\frac{-1547}{90}&\frac{1141}{90}&0&\frac{-421}{180}&\frac{1187}{180}&\frac{-229}{45}&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&\frac{827}{720}&\frac{-421}{180}&\frac{797}{720}&0&\frac{-109}{45}&\frac{887}{180}&\frac{-421}{180}&0&\frac{857}{720}&\frac{-109}{45}&\frac{827}{720}&0\\\
0&0&0&0&\frac{-421}{180}&\frac{278}{45}&\frac{-571}{180}&0&\frac{887}{180}&\frac{-1157}{90}&\frac{1187}{180}&0&\frac{-109}{45}&\frac{571}{90}&\frac{-293}{90}&0\\\
0&0&0&0&\frac{797}{720}&\frac{-571}{180}&\frac{1787}{720}&0&\frac{-421}{180}&\frac{1187}{180}&\frac{-229}{45}&0&\frac{827}{720}&\frac{-293}{90}&\frac{1817}{720}&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ \end{pmatrix},$

$A_{\text{sw}}=\setcounter{MaxMatrixCols}{16}\begin{pmatrix}\frac{857}{720}&\frac{-109}{45}&\frac{827}{720}&0&\frac{-109}{45}&\frac{887}{180}&\frac{-421}{180}&0&\frac{827}{720}&\frac{-421}{180}&\frac{797}{720}&0&0&0&0&0\\\
\frac{-109}{45}&\frac{571}{90}&\frac{-293}{90}&0&\frac{887}{180}&\frac{-1157}{90}&\frac{1187}{180}&0&\frac{-421}{180}&\frac{278}{45}&\frac{-571}{180}&0&0&0&0&0\\\
\frac{827}{720}&\frac{-293}{90}&\frac{1817}{720}&0&\frac{-421}{180}&\frac{1187}{180}&\frac{-229}{45}&0&\frac{797}{720}&\frac{-571}{180}&\frac{1787}{720}&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
\frac{-109}{45}&\frac{887}{180}&\frac{-421}{180}&0&\frac{571}{90}&\frac{-1157}{90}&\frac{278}{45}&0&\frac{-293}{90}&\frac{1187}{180}&\frac{-571}{180}&0&0&0&0&0\\\
\frac{887}{180}&\frac{-1157}{90}&\frac{1187}{180}&0&\frac{-1157}{90}&\frac{1472}{45}&\frac{-1547}{90}&0&\frac{1187}{180}&\frac{-1547}{90}&\frac{1667}{180}&0&0&0&0&0\\\
\frac{-421}{180}&\frac{1187}{180}&\frac{-229}{45}&0&\frac{278}{45}&\frac{-1547}{90}&\frac{1141}{90}&0&\frac{-571}{180}&\frac{1667}{180}&\frac{-623}{90}&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
\frac{827}{720}&\frac{-421}{180}&\frac{797}{720}&0&\frac{-293}{90}&\frac{1187}{180}&\frac{-571}{180}&0&\frac{1817}{720}&\frac{-229}{45}&\frac{1787}{720}&0&0&0&0&0\\\
\frac{-421}{180}&\frac{278}{45}&\frac{-571}{180}&0&\frac{1187}{180}&\frac{-1547}{90}&\frac{1667}{180}&0&\frac{-229}{45}&\frac{1141}{90}&\frac{-623}{90}&0&0&0&0&0\\\
\frac{797}{720}&\frac{-571}{180}&\frac{1787}{720}&0&\frac{-571}{180}&\frac{1667}{180}&\frac{-623}{90}&0&\frac{1787}{720}&\frac{-623}{90}&\frac{3497}{720}&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ \end{pmatrix},$

$A_{\text{se}}=\setcounter{MaxMatrixCols}{16}\begin{pmatrix}0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&\frac{1817}{720}&\frac{-293}{90}&\frac{827}{720}&0&\frac{-229}{45}&\frac{1187}{180}&\frac{-421}{180}&0&\frac{1787}{720}&\frac{-571}{180}&\frac{797}{720}&0&0&0&0\\\
0&\frac{-293}{90}&\frac{571}{90}&\frac{-109}{45}&0&\frac{1187}{180}&\frac{-1157}{90}&\frac{887}{180}&0&\frac{-571}{180}&\frac{278}{45}&\frac{-421}{180}&0&0&0&0\\\
0&\frac{827}{720}&\frac{-109}{45}&\frac{857}{720}&0&\frac{-421}{180}&\frac{887}{180}&\frac{-109}{45}&0&\frac{797}{720}&\frac{-421}{180}&\frac{827}{720}&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&\frac{-229}{45}&\frac{1187}{180}&\frac{-421}{180}&0&\frac{1141}{90}&\frac{-1547}{90}&\frac{278}{45}&0&\frac{-623}{90}&\frac{1667}{180}&\frac{-571}{180}&0&0&0&0\\\
0&\frac{1187}{180}&\frac{-1157}{90}&\frac{887}{180}&0&\frac{-1547}{90}&\frac{1472}{45}&\frac{-1157}{90}&0&\frac{1667}{180}&\frac{-1547}{90}&\frac{1187}{180}&0&0&0&0\\\
0&\frac{-421}{180}&\frac{887}{180}&\frac{-109}{45}&0&\frac{278}{45}&\frac{-1157}{90}&\frac{571}{90}&0&\frac{-571}{180}&\frac{1187}{180}&\frac{-293}{90}&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&\frac{1787}{720}&\frac{-571}{180}&\frac{797}{720}&0&\frac{-623}{90}&\frac{1667}{180}&\frac{-571}{180}&0&\frac{3497}{720}&\frac{-623}{90}&\frac{1787}{720}&0&0&0&0\\\
0&\frac{-571}{180}&\frac{278}{45}&\frac{-421}{180}&0&\frac{1667}{180}&\frac{-1547}{90}&\frac{1187}{180}&0&\frac{-623}{90}&\frac{1141}{90}&\frac{-229}{45}&0&0&0&0\\\
0&\frac{797}{720}&\frac{-421}{180}&\frac{827}{720}&0&\frac{-571}{180}&\frac{1187}{180}&\frac{-293}{90}&0&\frac{1787}{720}&\frac{-229}{45}&\frac{1817}{720}&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
\end{pmatrix}.$

The matrix $A_{\text{opt}}$ appearing in (3.17) is splitted into two matrices
$A_{\text{opt}}^{(1)}$ and $A_{\text{opt}}^{(2)}$ reporting the first eight
columns and the second eight columns of $A_{\text{opt}}$.

$A_{\text{opt}}^{(1)}=\begin{pmatrix}\frac{2903}{1575}&\frac{-182131}{37800}&\frac{10859}{2700}&\frac{-13889}{12600}&\frac{-182131}{37800}&\frac{60463}{4800}&\frac{-58881}{5600}&\frac{871553}{302400}\\\
\frac{-182131}{37800}&\frac{48799}{3150}&\frac{-58687}{4200}&\frac{10859}{2700}&\frac{60463}{4800}&\frac{-340087}{8400}&\frac{1226693}{33600}&\frac{-58881}{5600}\\\
\frac{10859}{2700}&\frac{-58687}{4200}&\frac{48799}{3150}&\frac{-182131}{37800}&\frac{-58881}{5600}&\frac{1226693}{33600}&\frac{-340087}{8400}&\frac{60463}{4800}\\\
\frac{-13889}{12600}&\frac{10859}{2700}&\frac{-182131}{37800}&\frac{2903}{1575}&\frac{871553}{302400}&\frac{-58881}{5600}&\frac{60463}{4800}&\frac{-182131}{37800}\\\
\frac{-182131}{37800}&\frac{60463}{4800}&\frac{-58881}{5600}&\frac{871553}{302400}&\frac{48799}{3150}&\frac{-340087}{8400}&\frac{426623}{12600}&\frac{-11137}{1200}\\\
\frac{60463}{4800}&\frac{-340087}{8400}&\frac{1226693}{33600}&\frac{-58881}{5600}&\frac{-340087}{8400}&\frac{815527}{6300}&\frac{-2958593}{25200}&\frac{426623}{12600}\\\
\frac{-58881}{5600}&\frac{1226693}{33600}&\frac{-340087}{8400}&\frac{60463}{4800}&\frac{426623}{12600}&\frac{-2958593}{25200}&\frac{815527}{6300}&\frac{-340087}{8400}\\\
\frac{871553}{302400}&\frac{-58881}{5600}&\frac{60463}{4800}&\frac{-182131}{37800}&\frac{-11137}{1200}&\frac{426623}{12600}&\frac{-340087}{8400}&\frac{48799}{3150}\\\
\frac{10859}{2700}&\frac{-58881}{5600}&\frac{221047}{25200}&\frac{-363499}{151200}&\frac{-58687}{4200}&\frac{1226693}{33600}&\frac{-1536757}{50400}&\frac{841823}{100800}\\\
\frac{-58881}{5600}&\frac{426623}{12600}&\frac{-1536757}{50400}&\frac{221047}{25200}&\frac{1226693}{33600}&\frac{-2958593}{25200}&\frac{10709107}{100800}&\frac{-1536757}{50400}\\\
\frac{221047}{25200}&\frac{-1536757}{50400}&\frac{426623}{12600}&\frac{-58881}{5600}&\frac{-1536757}{50400}&\frac{10709107}{100800}&\frac{-2958593}{25200}&\frac{1226693}{33600}\\\
\frac{-363499}{151200}&\frac{221047}{25200}&\frac{-58881}{5600}&\frac{10859}{2700}&\frac{841823}{100800}&\frac{-1536757}{50400}&\frac{1226693}{33600}&\frac{-58687}{4200}\\\
\frac{-13889}{12600}&\frac{871553}{302400}&\frac{-363499}{151200}&\frac{66427}{100800}&\frac{10859}{2700}&\frac{-58881}{5600}&\frac{221047}{25200}&\frac{-363499}{151200}\\\
\frac{871553}{302400}&\frac{-11137}{1200}&\frac{841823}{100800}&\frac{-363499}{151200}&\frac{-58881}{5600}&\frac{426623}{12600}&\frac{-1536757}{50400}&\frac{221047}{25200}\\\
\frac{-363499}{151200}&\frac{841823}{100800}&\frac{-11137}{1200}&\frac{871553}{302400}&\frac{221047}{25200}&\frac{-1536757}{50400}&\frac{426623}{12600}&\frac{-58881}{5600}\\\
\frac{66427}{100800}&\frac{-363499}{151200}&\frac{871553}{302400}&\frac{-13889}{12600}&\frac{-363499}{151200}&\frac{221047}{25200}&\frac{-58881}{5600}&\frac{10859}{2700}\\\
\end{pmatrix},$

$A_{\text{opt}}^{(2)}=\begin{pmatrix}\frac{10859}{2700}&\frac{-58881}{5600}&\frac{221047}{25200}&\frac{-363499}{151200}&\frac{-13889}{12600}&\frac{871553}{302400}&\frac{-363499}{151200}&\frac{66427}{100800}\\\
\frac{-58881}{5600}&\frac{426623}{12600}&\frac{-1536757}{50400}&\frac{221047}{25200}&\frac{871553}{302400}&\frac{-11137}{1200}&\frac{841823}{100800}&\frac{-363499}{151200}\\\
\frac{221047}{25200}&\frac{-1536757}{50400}&\frac{426623}{12600}&\frac{-58881}{5600}&\frac{-363499}{151200}&\frac{841823}{100800}&\frac{-11137}{1200}&\frac{871553}{302400}\\\
\frac{-363499}{151200}&\frac{221047}{25200}&\frac{-58881}{5600}&\frac{10859}{2700}&\frac{66427}{100800}&\frac{-363499}{151200}&\frac{871553}{302400}&\frac{-13889}{12600}\\\
\frac{-58687}{4200}&\frac{1226693}{33600}&\frac{-1536757}{50400}&\frac{841823}{100800}&\frac{10859}{2700}&\frac{-58881}{5600}&\frac{221047}{25200}&\frac{-363499}{151200}\\\
\frac{1226693}{33600}&\frac{-2958593}{25200}&\frac{10709107}{100800}&\frac{-1536757}{50400}&\frac{-58881}{5600}&\frac{426623}{12600}&\frac{-1536757}{50400}&\frac{221047}{25200}\\\
\frac{-1536757}{50400}&\frac{10709107}{100800}&\frac{-2958593}{25200}&\frac{1226693}{33600}&\frac{221047}{25200}&\frac{-1536757}{50400}&\frac{426623}{12600}&\frac{-58881}{5600}\\\
\frac{841823}{100800}&\frac{-1536757}{50400}&\frac{1226693}{33600}&\frac{-58687}{4200}&\frac{-363499}{151200}&\frac{221047}{25200}&\frac{-58881}{5600}&\frac{10859}{2700}\\\
\frac{48799}{3150}&\frac{-340087}{8400}&\frac{426623}{12600}&\frac{-11137}{1200}&\frac{-182131}{37800}&\frac{60463}{4800}&\frac{-58881}{5600}&\frac{871553}{302400}\\\
\frac{-340087}{8400}&\frac{815527}{6300}&\frac{-2958593}{25200}&\frac{426623}{12600}&\frac{60463}{4800}&\frac{-340087}{8400}&\frac{1226693}{33600}&\frac{-58881}{5600}\\\
\frac{426623}{12600}&\frac{-2958593}{25200}&\frac{815527}{6300}&\frac{-340087}{8400}&\frac{-58881}{5600}&\frac{1226693}{33600}&\frac{-340087}{8400}&\frac{60463}{4800}\\\
\frac{-11137}{1200}&\frac{426623}{12600}&\frac{-340087}{8400}&\frac{48799}{3150}&\frac{871553}{302400}&\frac{-58881}{5600}&\frac{60463}{4800}&\frac{-182131}{37800}\\\
\frac{-182131}{37800}&\frac{60463}{4800}&\frac{-58881}{5600}&\frac{871553}{302400}&\frac{2903}{1575}&\frac{-182131}{37800}&\frac{10859}{2700}&\frac{-13889}{12600}\\\
\frac{60463}{4800}&\frac{-340087}{8400}&\frac{1226693}{33600}&\frac{-58881}{5600}&\frac{-182131}{37800}&\frac{48799}{3150}&\frac{-58687}{4200}&\frac{10859}{2700}\\\
\frac{-58881}{5600}&\frac{1226693}{33600}&\frac{-340087}{8400}&\frac{60463}{4800}&\frac{10859}{2700}&\frac{-58687}{4200}&\frac{48799}{3150}&\frac{-182131}{37800}\\\
\frac{871553}{302400}&\frac{-58881}{5600}&\frac{60463}{4800}&\frac{-182131}{37800}&\frac{-13889}{12600}&\frac{10859}{2700}&\frac{-182131}{37800}&\frac{2903}{1575}\\\
\end{pmatrix}.$