# Local fully discrete entropy stability for a genuine second-order scheme for scalar hyperbolic equations

Christophe Berthon, Ludovic Martaud

<figure><img src="image_1.png" /></figure>

## To cite this version:

Christophe Berthon, Ludovic Martaud. Local fully discrete entropy stability for a genuine second-order scheme for scalar hyperbolic equations. 2025. hal-05248101

HAL Id: hal-05248101

https://hal.science/hal-05248101v1

Preprint submitted on 16 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Local fully discrete entropy stability for a genuine second-order scheme for scalar hyperbolic equations

Christophe Berthon,\* Ludovic Martaud†

September 10, 2025

## Abstract

This work concerns the numerical approximations of the weak solutions of scalar hyperbolic conservation laws. After showing how to bypass the barrier theorems for the linear advection, the derivation of a second-order entropy-satisfying scheme is presented for non-linear equations. The fully discrete stability result is established for regular strictly convex entropy and under a parabolic CFL-like condition. Some numerical experiments are done to assess the accuracy and the stability of the proposed scheme.

## 1 Introduction

The present work concerns the numerical approximation of the weak solutions of scalar conservation laws in the one space dimension on a bounded domain $I \subset \mathbb{R}$ given by

$$ \left\{ \begin{array}{l} \partial_t w + \partial_x f(w) = 0, \quad x \in I, t > 0, \\ w(x, t = 0) = w_0(x), \quad x \in I. \end{array} \right. \qquad (1) $$

The unknown is $w(x, t) \in \mathbb{R}$, $f: \mathbb{R} \to \mathbb{R}$ is a smooth given flux function and $w_0: I \to \mathbb{R}$ is a given measurable function in $L^1(I)$. For the sake of simplicity, only periodic boundary conditions are considered in this work.

According to [29, 30, 39] (see also [18, 19, 32, 45]), it is well-known that the solutions of (1) may develop, in a finite time, discontinuities and that the weak solutions are in general non unique. In order to rule out non-admissible weak solutions, (1) must be endowed with entropy inequalities. In this regard, weak solutions of (1) containing discontinuities, verify an entropy inequality (for instance, see [29, 30, 39]) given by

$$ \partial_t \eta(w) + \partial_x G(w) \le 0 \quad \text{in} \quad \mathcal{D}'(I \times (0, +\infty)), \qquad (2) $$

where $\eta: \mathbb{R} \to \mathbb{R}$ is a convex function and $G: \mathbb{R} \to \mathbb{R}$ is the entropy flux function satisfying $\eta'(w)f'(w) = G'(w)$ for all $w \in \mathbb{R}$. A weak solution of (1) is called an entropy-satisfying solution if and only if the entropy inequality (2) holds for any pair entropy-entropy flux $(\eta, G)$.

As for the numerical approximation, the weak solutions of (1) are approximated within the standard finite volume framework. In this regard and considering $N \in \mathbb{N}^*$, the periodic-space domain $I$ and time interval $[0, +\infty)$ are discretized by uniform meshes of size $\Delta x, \Delta t > 0$ such that $(x_{i+\frac{1}{2}}, t^{n+1}) = (x_{i-\frac{1}{2}}, t^n) + (\Delta x, \Delta t)$, for all $(i, n) \in \mathbb{Z}/\mathbb{N}\mathbb{Z} \times \mathbb{N}$. Considering such meshes, the numerical approximation of the weak solutions of (1) is given by the following piecewise constant function:

$$ w_{\Delta}(x, t^n) = w_i^n \approx \frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} w(x, t^n) \, dx, \quad \text{if } x \in (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}), \quad \forall i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}, \qquad (3) $$

<footnote>\*Laboratoire de Mathématiques Jean Leray, CNRS UMR 6629, Université de Nantes, 2 rue de la Houssinière, BP 92208, 44322 Nantes, France, Email address: christophe.berthon@univ-nantes.fr</footnote>

<footnote>†Univ Rennes, Inria Bretagne Atlantique Email address: ludovic.martaud@inria.fr.</footnote>
and the discrete periodic boundary conditions are enforced with the sequence $(w_i^n)_{i \in \mathbb{Z}/N\mathbb{Z}}$. Over the past fifty years, numerous strategies have been proposed to evolve in time the approximation (3) and to define suitable updated states $(w_i^{n+1})_{i \in \mathbb{Z}/N\mathbb{Z}}$ (for instance, see [11, 18, 19, 21, 32, 37, 45] and references therein). In the present work, three-points conservative finite volume schemes are used, so that the updated state reads

$$w_i^{n+1} = w_i^n - \frac{\Delta t}{\Delta x} (\mathcal{F}(w_i^n, w_{i+1}^n) - \mathcal{F}(w_{i-1}^n, w_i^n)) \quad (4)$$

where $\mathcal{F}: (\mathbb{R})^2 \to \mathbb{R}$ is a numerical flux function satisfying the consistency relation $\mathcal{F}(w, w) = f(w)$ (see [19, 21, 45] for the details). The consistency of the numerical flux function is required for the convergence of (4) but the obtained limit solution is not necessarily entropy-satisfying. As a consequence, non-admissible discontinuous waves may appear (for instance, see [9, 26]). To correct such irrelevant solutions, $(w_i^{n+1})_{i \in \mathbb{Z}/N\mathbb{Z}}$ has to satisfy discrete entropy inequalities in the form

$$\eta(w_i^{n+1}) \le \eta(w_i^n) - \frac{\Delta t}{\Delta x} (\mathcal{G}(w_i^n, w_{i+1}^n) - \mathcal{G}(w_{i-1}^n, w_i^n)) \quad (5)$$

where $\mathcal{G}: (\mathbb{R})^2 \to \mathbb{R}$ is a numerical entropy flux function satisfying the consistency relation $\mathcal{G}(w, w) = G(w)$. Discrete entropy inequalities (5) are satisfied by several first-order schemes, among of them the exact Godunov scheme [20, 21], the kinetic schemes [3, 27], the HLL scheme [21, 22], the HLLC scheme [44], some relaxation schemes [2, 4, 10], the numerical strategy introduced by Tadmor [42, 43] or the staggered schemes [16, 23, 24].

As far as high-order numerical approximations are concerned, the situation turns out to be drastically distinct and very challenging. Several works devoted to high-order schemes attempted to exhibit discrete entropy inequalities (5) but, from a general point of view, the time discretization is an important technical obstacle. Hence, semi-discrete entropy estimates $\frac{d}{dt}\eta(w_i(t)) \le -\frac{1}{\Delta x}(\mathcal{G}(w_i(t), w_{i+1}(t)) - \mathcal{G}(w_{i-1}(t), w_i(t)))$, are often considered for finite volumes schemes [7, 8, 17, 43] or within the Discontinuous Galerkin (DG) finite element framework [28, 36]. Unfortunately, such semi-discrete entropy inequalities are not sufficient to rule out non-admissible discontinuities in the converged solutions [6]. In [2, 5, 34, 35], fully discrete entropy estimates are introduced as follows:

$$\eta(w_i^{n+1}) \le \bar{\eta}_i^n - \frac{\Delta t}{\Delta x} (\mathcal{G}(w_i^n, w_{i+1}^n) - \mathcal{G}(w_{i-1}^n, w_i^n)) \quad (6)$$

for a suitable entropy average $\bar{\eta}_i^n \approx \eta(w_i^n)$. Such discrete entropy inequality strategies are not fully relevant since Lax-Wendroff Theorem [31] cannot be successfully applied (see [14] for the details). In addition, in [14], numerical experiments exhibited the capture of non-admissible shock solutions for MUSCL schemes which satisfy (6).

This work concerns the design of a second-order scheme (4) also satisfying a local fully discrete entropy inequality (5). The paper is organized as follows. Section 2 concerns a discussion about the celebrate order barrier theorems by Godunov in [20] and Schonbek in [38]. According to these two crucial papers, it comes (for instance see also [19, 45]): that “if a 3-point finite volume scheme (4) satisfies the discrete entropy inequality (5) then it is at most first-order accurate. Involving to a suitable definition of the order of accuracy, we exhibit the limitations of the well-known order barrier theorems. Moreover, adopting a very easy linear scheme for the usual advection equation, we present a high-order, i.e. larger than 2, in the form (4) which satisfies the discrete entropy inequality (5). In Section 3, a second-order scheme is derived and its discrete entropy stability property is given. Section 4 is devoted to the proof of the second-order accuracy while Section 5 concerns the proof of the entropy stability. Section 6 is finally dedicated to the numerical experiments.

## 2 About the order limitation of entropy-satisfying three-points finite volume schemes

Considering second-order, eventually high-order, finite volume schemes endowed with a discrete entropy inequality, we immediately have to emphasize that an entropy satisfying 3-point scheme (4) is at most
first-order accurate. The reader is easily convinced by such a restriction after the order barrier theorems by Godunov [20] and Schonbek [38]. The believe of this order limitation is enforced, from several decades, by numerous books and articles (for instance, see [19, 45]). Arguing some modern definition of the order of accuracy, in this section we show how bypass the barriers imposed by Godunov [20] and Schonbek [38].

First, let us consider the order restrictions enforced by Godunov. Indeed, in [20], the author considers the discretization of the linear advection equation to claim: Among schemes of second-order accuracy for the (linear advection) equation there is none which satisfies the monotonicity condition. To address such a statement, Godunov adopt second-order schemes that remain exact for initial data given by a polynomial of second order degree. According to such a second-order scheme definition, Godunov proves the impossibility to get a monotone scheme.

Next, we recall that Schonbek states in [38] that a scheme, given by (4), is necessarily first-order to be entropy satisfying according to (5). This time, to establish this result, the author adopts second-order scheme such that $|\mathcal{F}(w_i^n, w_{i+1}^n) - \mathcal{F}^{LW}(w_i^n, w_{i+1}^n)| \le \mathcal{O}(|w_{i+1}^n - w_i^n|^2)$, where $\mathcal{F}^{LW}(w_i^n, w_{i+1}^n)$ stands for the second-order Lax-Wendroff numerical flux function [19, 31].

From now on, we notice that Godunov [20] and Schonbek [38] adopt two distinct definitions of the order of accuracy. But, nowadays the definition of the order of accuracy turns out to be weaker. In the present section, we adopt the certainly more common definitions; namely the truncation errors [19] and the order of accuracy in the weak sense [4]. Both definition are given as follow:

**Definition 2.1** (Truncation errors [19]). The space and time order of accuracy of a scheme (4) are the largest numbers $k \ge 1$ and $\ell \ge 1$ such that we have for any smooth solution $w$ of (1) and for $\frac{\Delta t}{\Delta x}$ kept constant, $\tau = \mathcal{O}(\Delta t^k) + \mathcal{O}(\Delta x^\ell)$ with $\Delta t, \Delta x \to 0$ where we have set

$$ \tau = \frac{1}{\Delta t} (w(x, t + \Delta t) - w(x, t)) + \frac{1}{\Delta x} (\mathcal{F}(w(x, t), w(x + \Delta x, t)) - \mathcal{F}(w(x - \Delta x, t), w(x, t))). $$

**Definition 2.2** (Order of accuracy in the weak sense [4]). For any smooth solution of (1), let us set $w_i^n = \frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} w(x, t^n) dx$. With a constant ratio $\frac{\Delta t}{\Delta x}$, the scheme (4) is said of k-order in time and $\ell$-order in space if

$$ \frac{1}{\Delta t} \left( w_i^{n+1} - \frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} w(x, t^n + \Delta t) dx \right) = \frac{1}{\Delta x} \left( \mathcal{R}_{i+\frac{1}{2}} - \mathcal{R}_{i-\frac{1}{2}} \right), \quad (7) $$

where $\mathcal{R}_{i+\frac{1}{2}} = \mathcal{O}(\Delta t^k) + \mathcal{O}(\Delta x^\ell)$ is a scheme residue.

From now on, let us underline that the numerical evaluation of the order of accuracy is, in general, obtained assessing $\| w^n - w(., t^n) \|_{L^1}$, $\| w^n - w(., t^n) \|_{L^2}$ and $\| w^n - w(., t^n) \|_{L^\infty}$. Such a practical numerical order of accuracy is nothing but a global extension of the weak sense of the order.

Moreover, since the ratio $\frac{\Delta t}{\Delta x}$ stays constant, we may consider a global order given by $\min(k, \ell)$ because of $\mathcal{O}(\Delta t^k) + \mathcal{O}(\Delta x^\ell) = \mathcal{O}(\Delta x^{\min(k, \ell)})$.

In fact, the above definitions are similar and they naturally contain the order of accuracy adopted by Godunov and by Schonbek. However, it is worth noticing that approximations with second-order according to the truncation error Definition 2.1 or to the weak sense Definition 2.2 may be first-order in the sense of Godunov or Schonbek. As a consequence, it is clearly possible to envisage a second-order scheme (4), eventually high-order, in a modern definitions 2.1 or 2.2, which stays Godunov and/or Schonbek first-order, such that a discrete entropy inequality (5) is satisfied. The well-known barrier theorems can be circumvented.

To illustrate our purpose, let us consider the linear advection equation given by $\partial_t w + a \partial_x w = 0$, with $a > 0$. To approximate the weak solutions of this equation, we adopt the usual upwind scheme [19, 45] given by

$$ w_i^{n+1} = w_i^n - \frac{a \Delta t}{\Delta x} (w_i^n - w_{i-1}^n), \quad (8) $$
with the CFL time restriction $\frac{a\Delta t}{\Delta x} \le 1$. This scheme is well-known to be first-order and to satisfy all the expected stability properties : $L^1$-stable, $L^2$-stable, $L^\infty$-stable, TVD, maximum principle preserving, monotonicity, entropy preserving (for instance, see [19, 45]). Now, we establish that (8) can be high-order.

**Lemma 2.1.** Let us set $C > 0$ a generic constant. With a constant ratio $\frac{\Delta t}{\Delta x}$, the truncation error of the scheme (8) reads $\tau = (1 - \frac{a\Delta t}{\Delta x}) C\Delta x$, while the scheme residue (7) in the weak sense order is $R_{i+\frac{1}{2}} = (1 - a\frac{\Delta t}{\Delta x}) C\Delta x$.

The proof of the above statement is classical and it is left to the reader. Of course, this result clearly establishes that the upwind scheme (8) is at least first-order.

However, it is possible to make very easily the upwind scheme (8) high-order. Indeed, let us fix the time increment as follows:

$$a \frac{\Delta t}{\Delta x} = 1 - \varepsilon \Delta x^{k-1}, \qquad (9)$$

with $k \ge 1$ and $\varepsilon > 0$ such that $1 - \varepsilon \Delta x^{k-1} > 0$. We emphasize that, as soon as $\Delta x$ is small enough, it is sufficient to get $\varepsilon = 1$.

Now, concerning the order of accuracy of the upwind scheme (8) with (9), we immediately have $\tau = \varepsilon C \Delta x^k$ and $R_{i+\frac{1}{2}} = \varepsilon C \Delta x^k$, so that the upwind scheme turns out to be $k$-order with $k \ge 1$. But we also remark that the upwind scheme (8) is not modified and the CFL condition $\frac{a\Delta t}{\Delta x} \le 1$ is preserved. As a consequence, we have established the following result.

**Theorem 2.1.** The upwind scheme (8) with a time step given by (9) is $k$-order with $k \ge 1$. Moreover, the scheme is $L^1$-stable, $L^2$-stable, $L^\infty$-stable, TVD, maximum principle preserving, monotonicity, and it satisfies a discrete entropy inequality (5) for all entropy pairs $(\eta, G)$.

Once again, we underline that the upwind scheme, even with a time step given by (9), remains first-order according to both Godunov and Schonbek order definitions.

Let us present some numerical experiments to illustrate our purpose. We adopt a computational domain given by $I = (-1, 1]$ with periodic boundary conditions. The transport velocity is fixed by $a = 2$ and $\varepsilon > 0$ in (9) satisfies $\varepsilon = 1$. First, we approximate a smooth solution obtained by considering an initial data given by $w_0(x) = 0.25 + 0.5 \sin(\pi x)$. The $L^1$, $L^2$- and $L^\infty$-errors are evaluated at the final time $t = 7$. The obtained numerical orders of accuracy are reported on the top of Table 1 for $k = 2$ and $k = 5$. As expected, we get the correct orders of accuracy. In the second simulation, a discontinuous initial data is used under the form $w_0(x) = \chi_{|x|<0.25}(x)$ where $\chi$ denotes the indicator function. Once again the numerical order of accuracy is evaluated by considering $L^1$, $L^2$- and $L^\infty$-errors. As presented on the bottom of Table 1, we obtain the correct order of the scheme which is given by $(k-2) + \frac{1}{p}$ in $L^p$ for discontinuous solutions. Moreover, in Figure 1 we present the performed approximate solutions for both initial data. We underline that usual spurious oscillations never occur since the upwind scheme (8) satisfied all the expected stability properties.

## 3 Local fully discrete entropy stability of a second-order finite volume scheme

This section concerns the definition of a second-order finite volume scheme satisfying (5). For the sake of completeness, the below lemma is given to characterize a formal second-order finite volume scheme (4) according to Definition 2.2.

**Lemma 3.1** (Weak consistency [4]). Let $u(x, t) \in \mathbb{R}$ be a smooth function and consider $u_i^n = \frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} u(x, t^n) \, dx$. Let $\mathfrak{F}: (\mathbb{R})^2 \to \mathbb{R}$ be a continuous function which verifies the consistency relation $\mathfrak{F}(u, u) = f(u)$ for all $u \in \mathbb{R}$ and which defines the following finite volume scheme:

$$u_i^{n+1} = u_i^n - \frac{\Delta t}{\Delta x} (\mathfrak{F}(u_i^n, u_{i+1}^n) - \mathfrak{F}(u_{i-1}^n, u_i^n)). \qquad (10)$$
<table><caption>Table 1: Numerical order of accuracy for the scheme (8), (9) with $k = 2$ (resp. $k = 5$) on the left, (resp. right) and endowed with $w_0(x) = 0.25 + 0.5 \sin(\pi x)$ (resp. $w_0(x) = \chi_{|x|\le 0.25}(x)$) on the top (resp. bottom).</caption><thead><tr><th>Cells</th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th><th>Cells</th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><td>100</td><td>-</td><td>-</td><td>-</td><td>100</td><td>-</td><td>-</td><td>-</td></tr><tr><td>200</td><td>1.9647</td><td>1.9646</td><td>1.9644</td><td>200</td><td>4.9991</td><td>4.9995</td><td>4.9995</td></tr><tr><td>400</td><td>1.9943</td><td>1.9943</td><td>1.9942</td><td>400</td><td>4.9998</td><td>4.9998</td><td>4.9998</td></tr><tr><td>800</td><td>2.0021</td><td>2.0021</td><td>2.0021</td><td>800</td><td>4.9999</td><td>4.9999</td><td>4.9999</td></tr><tr><td>1600</td><td>2.0023</td><td>2.0023</td><td>2.0023</td><td>1600</td><td>4.9999</td><td>4.9999</td><td>4.9999</td></tr><tr><td>3200</td><td>2.0014</td><td>2.0014</td><td>2.0014</td><td>3200</td><td>4.9999</td><td>4.9999</td><td>4.9999</td></tr></tbody></table><table><thead><tr><th>Cells</th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th><th>Cells</th><th>L<sup>1</sup></th><th>L<sup>2</sup></th><th>L<sup>∞</sup></th></tr></thead><tbody><tr><td>102</td><td>-</td><td>-</td><td>-</td><td>102</td><td>-</td><td>-</td><td>-</td></tr><tr><td>202</td><td>0.9972</td><td>0.4990</td><td>-0.0016</td><td>202</td><td>3.9993</td><td>3.4991</td><td>2.9993</td></tr><tr><td>402</td><td>1.0009</td><td>0.5004</td><td>-0.0005</td><td>402</td><td>3.9999</td><td>3.4998</td><td>2.9999</td></tr><tr><td>802</td><td>1.0020</td><td>0.5010</td><td>-0.0003</td><td>802</td><td>3.9999</td><td>3.4999</td><td>2.9999</td></tr><tr><td>1602</td><td>1.0014</td><td>0.5006</td><td>-0.0002</td><td>1602</td><td>3.9999</td><td>3.4999</td><td>2.9999</td></tr><tr><td>3202</td><td>1.0008</td><td>0.5003</td><td>-0.0001</td><td>3202</td><td>3.9999</td><td>3.4999</td><td>2.9999</td></tr></tbody></table>

If $\mathfrak{F}$ satisfies $\mathfrak{F}(u_i^n, u_{i+1}^n) = f(u(x_{i+\frac{1}{2}}, t^n)) + \mathcal{O}(\Delta x^2)$, where the $\mathcal{O}(\Delta x^2)$ term is uniform with respect to $i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}$ then (10) is second-order consistent in the weak sense.

*Proof.* The reader can be refer to [1, 4] for the details of the proof. $\square$

<figure><img src="image_3.png" /><figcaption>Figure 1: Numerical simulations for the scheme (8), (9) with $w_0(x) = 0.25+0.5 \sin(\pi x)$ (resp. $w_0(x) = \chi_{|x|<0.25}(x)$) on the left (resp. center and right).</figcaption></figure>

Equipped with the above lemma, the proposed scheme (4) in this work is defined by the following numerical flux:

$$ \mathcal{F}(w_i^n, w_{i+1}^n) = \frac{1}{2}(f(w_i^n) + f(w_{i+1}^n)) - \frac{1}{2}\Lambda_{i+\frac{1}{2}}^n(w_{i+1}^n - w_i^n), \quad (11a) $$

$$ \Lambda_{i+\frac{1}{2}}^n = \gamma_{i+\frac{1}{2}}^n |w_{i+1}^n - w_i^n| + \sigma_{i+\frac{1}{2}}^n \Delta x, \quad (11b) $$

where $\gamma_{i+\frac{1}{2}}^n, \sigma_{i+\frac{1}{2}}^n > 0$ have to be restricted according to the entropy stability (5). The main originality in (11) is the term $\Lambda_{i+\frac{1}{2}}^n = \gamma_{i+\frac{1}{2}}^n |w_{i+1}^n - w_i^n| + \sigma_{i+\frac{1}{2}}^n \Delta x$ which stands for the artificial viscosity. For smooth solutions, this viscosity has $\mathcal{O}(\Delta x)$ magnitude and thus it ensures the second-order accuracy. For discontinuous solutions, the sequel shows the relevancy of $\Lambda_{i+\frac{1}{2}}^n$ to enforce (5). As a consequence, the scheme (4), (11) does not involve limitations techniques in contrast with other usual approaches as MUSCL technique [46], (W)ENO schemes [33, 40, 41] or DG schemes [15, 25, 41].
In the forthcoming developments, the entropy-entropy flux pair $(\eta, G)$ is assumed to belong to $(C^3(\mathbb{R}, \mathbb{R}))^2$. For the sake of clarity in the notations, sequences $(\delta_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}}$, $(\psi_{i+\frac{1}{2}}^{\pm}(s))_{(i,s) \in \mathbb{Z}/\mathbb{N}\mathbb{Z} \times [0,1]}$, $(\Delta \mathcal{F}_i^n)_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}}$ now are introduced as follows:

$$ \delta_{i+\frac{1}{2}}^n = w_{i+1}^n - w_i^n, \qquad (12a) $$

$$ \begin{aligned} \psi_{i+\frac{1}{2}}^{+}(s) = & \frac{1}{2}(1-s)^2(G'''(w_i^n + s\delta_{i+\frac{1}{2}}^n) - \eta'(w_i^n)f''(w_i^n + s\delta_{i+\frac{1}{2}}^n)) \\ & - \frac{1}{4}(\eta''f')'(w_i^n + s\delta_{i+\frac{1}{2}}^n), \end{aligned} \qquad (12b) $$

$$ \begin{aligned} \psi_{i+\frac{1}{2}}^{-}(s) = & \frac{1}{2}(1-s)^2(G'''(w_{i+1}^n - s\delta_{i+\frac{1}{2}}^n) - \eta'(w_{i+1}^n)f''(w_{i+1}^n - s\delta_{i+\frac{1}{2}}^n)) \\ & - \frac{1}{4}(\eta''f')'(w_{i+1}^n - s\delta_{i+\frac{1}{2}}^n), \end{aligned} \qquad (12c) $$

$$ \Delta \mathcal{F}_i^n = \mathcal{F}(w_i^n, w_{i+1}^n) - \mathcal{F}(w_{i-1}^n, w_i^n). \qquad (12d) $$

The main result about the scheme (4), (11) is the following:

**Theorem 3.1** (Local fully discrete entropy stability of a second-order finite volume scheme). Consider $(\eta, G) \in (C^3(\mathbb{R}, \mathbb{R}))^2$ a strictly convex entropy-entropy flux pair (i.e. $\eta'' > 0$) which endowed (1). Let $(\gamma_{i+\frac{1}{2}}^n, \sigma_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}}$ be bounded sequences as $\Delta x \to 0$ and that define $(\Lambda_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}} \in \mathbb{R}_+^*$ according to (11b). The following hold:

(i) The scheme (4), (11) defines a second-order accuracy scheme in the sense of Lemma 3.1.

(ii) In addition, if $\sigma_{i+\frac{1}{2}}^n > 0$ and if $\gamma_{i+\frac{1}{2}}^n > 0$ is such that

$$ \gamma_{i+\frac{1}{2}}^n \geq \max \left( \frac{|\int_0^1 \psi_{i+\frac{1}{2}}^+(s) \, ds|}{\int_0^1 (1-s)\eta''(w_i^n + s\delta_{i+\frac{1}{2}}^n) \, ds}, \frac{|\int_0^1 \psi_{i+\frac{1}{2}}^-(s) \, ds|}{\int_0^1 (1-s)\eta''(w_{i+1}^n - s\delta_{i+\frac{1}{2}}^n) \, ds} \right), \qquad (13) $$

for all $i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}$, then there exists $\frac{\Delta t}{\Delta x^2} > 0$ small enough such that the updated sequence $(w_i^{n+1})_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}}$ given by (4), (11) satisfies a local fully discrete entropy inequality (5) with the following numerical entropy flux:

$$ \begin{aligned} \mathcal{G}(w_i^n, w_{i+1}^n) = & \frac{1}{2}(G(w_i^n) + G(w_{i+1}^n)) - \frac{1}{2}\Lambda_{i+\frac{1}{2}}^n(\eta(w_{i+1}^n) - \eta(w_i^n)) \\ & - \frac{1}{8}((\eta''f')(w_i^n) + (\eta''f')(w_{i+1}^n))(w_{i+1}^n - w_i^n)^2. \end{aligned} \qquad (14) $$

The constraint (13) is fully explicit and, according to the strictly convexity of $\eta(w)$, $\gamma_{i+\frac{1}{2}}^n$ remains bounded. The formulation (13) may be also considerably simplified depending on $f(w)$ in (1) and depending on the couple $(\eta, G)$. In this regard, Section 6 details several cases. Section 6 also presents practical heuristic choices for the restriction $\frac{\Delta t}{\Delta x^2}$ which is implicitly given in Theorem 3.1. The sequel deals with the proof of Theorem 3.1 and the following section concerns the second-order accuracy.

# 4 Proof of the second-order accuracy

This section concerns the proof of the second-order accuracy of the scheme (4), (11). The proof is given below.

*Proof of Theorem 3.1-(i).* According to Lemma 3.1, the proof will be complete as soon as $\mathcal{F}(u_i^n, u_{i+1}^n) = f((u(x_{i+\frac{1}{2}}, t^n)) + \mathcal{O}(\Delta x^2)$ will be established with $u_i^n = \frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} u(x, t^n) \, dx$ and $\mathcal{F}(u_i^n, u_{i+1}^n)$ given by (11).
Since $u(x, t)$ is regular, a Taylor expansion in the neighborhood of $x_{i+\frac{1}{2}}$ reads

$$ 
\begin{aligned}
u_i^n &= \frac{1}{\Delta x} \int_{-\Delta x}^{0} u(x + x_{i+\frac{1}{2}}, t^n) \, dx = u(x_{i+\frac{1}{2}}, t^n) - \frac{\Delta x}{2} \partial_x u(x_{i+\frac{1}{2}}, t^n) + \mathcal{O}(\Delta x^2), \\
u_{i+1}^n &= \frac{1}{\Delta x} \int_{0}^{\Delta x} u(x + x_{i+\frac{1}{2}}, t^n) \, dx = u(x_{i+\frac{1}{2}}, t^n) + \frac{\Delta x}{2} \partial_x u(x_{i+\frac{1}{2}}, t^n) + \mathcal{O}(\Delta x^2).
\end{aligned}
$$

As a consequence, using the above Taylor expansions in the definition of $\Lambda_{i+\frac{1}{2}}^n$ given in (11b), the following holds:

$$ 
\begin{aligned}
\Lambda_{i+\frac{1}{2}}^n (u_{i+1}^n - u_i^n) &= \left( \gamma_{i+\frac{1}{2}}^n |\Delta x \partial_x u(x_{i+\frac{1}{2}}, t^n)| + \mathcal{O}(\Delta x^2) \right) + \sigma_{i+\frac{1}{2}}^n \Delta x \left( \Delta x \partial_x u(x_{i+\frac{1}{2}}, t^n) + \mathcal{O}(\Delta x^2) \right), \\
&= \Delta x^2 \left( \gamma_{i+\frac{1}{2}}^n |\partial_x u(x_{i+\frac{1}{2}}, t^n)| + \mathcal{O}(\Delta x) \right) + \sigma_{i+\frac{1}{2}}^n \left( \partial_x u(x_{i+\frac{1}{2}}, t^n) + \mathcal{O}(\Delta x) \right).
\end{aligned}
\quad (15) $$

If $\partial_x u(x_{i+\frac{1}{2}}, t^n) = 0$ and since $\gamma_{i+\frac{1}{2}}^n, \sigma_{i+\frac{1}{2}}^n > 0$ are bounded as $\Delta x \to 0$ the above formulation gives $\Lambda_{i+\frac{1}{2}}^n(u_{i+1}^n - u_i^n) = \mathcal{O}(\Delta x^2)$. Now, if $\partial_x u(x_{i+\frac{1}{2}}, t^n) > 0$ then there exists $\Delta x > 0$ small enough such that $|\partial_x u(x_{i+\frac{1}{2}}, t^n) + \mathcal{O}(\Delta x)| = \partial_x u(x_{i+\frac{1}{2}}, t^n) + \mathcal{O}(\Delta x)$. In this case, (15) reads

$$ 
\begin{aligned}
\Lambda_{i+\frac{1}{2}}^n (u_{i+1}^n - u_i^n) &= \Delta x^2 \left( \gamma_{i+\frac{1}{2}}^n (\partial_x u(x_{i+\frac{1}{2}}, t^n))^2 + \sigma_{i+\frac{1}{2}}^n \partial_x u(x_{i+\frac{1}{2}}, t^n) + \mathcal{O}(\Delta x) \right) \\
&= \mathcal{O}(\Delta x^2).
\end{aligned}
$$

Since the case $\partial_x u(x_{i+\frac{1}{2}}, t^n) < 0$ would be derived with similar computations, it results $\Lambda_{i+\frac{1}{2}}^n(u_{i+1}^n - u_i^n) = \mathcal{O}(\Delta x^2)$.

On the other hands, the regularity of $u \mapsto f(u)$ and Taylor expansions in the neighborhood of $x_{i+\frac{1}{2}}$ lead to

$$ 
\begin{aligned}
\frac{1}{2} (f(u_i^n) + f(u_{i+1}^n)) &= \frac{1}{2} f(u(x_{i+\frac{1}{2}}, t^n)) - \frac{\Delta x}{2} \partial_x f(u(x_{i+\frac{1}{2}}, t^n)) + \mathcal{O}(\Delta x^2) \\
&\quad + \frac{1}{2} f(u(x_{i+\frac{1}{2}}, t^n)) + \frac{\Delta x}{2} \partial_x f(u(x_{i+\frac{1}{2}}, t^n)) + \mathcal{O}(\Delta x^2), \\
&= f(u(x_{i+\frac{1}{2}}, t^n)) + \mathcal{O}(\Delta x^2).
\end{aligned}
$$

Using the above estimate and $\Lambda_{i+\frac{1}{2}}^n(u_{i+1}^n - u_i^n) = \mathcal{O}(\Delta x^2)$ in the definition of $\mathcal{F}$ given by (11), the equality $\mathcal{F}(u_i^n, u_{i+1}^n) = f(u(x_{i+\frac{1}{2}}, t^n)) + \mathcal{O}(\Delta x^2)$ follows, that concludes the proof. $\square$

According to the above proof, $\Lambda_{i+\frac{1}{2}}^n(u_{i+1}^n - u_i^n) = \mathcal{O}(\Delta x^2)$ plays a central role for the accuracy of the scheme (4), (11). This quantity is also essential for the discrete entropy stability that is proved in the next section.

# 5 Proof of the local fully discrete entropy stability

This section is dedicated to the proof of the entropy stability (5). Since $\eta \in C^2(\mathbb{R}, \mathbb{R})$, a Taylor expansion of $\eta(w_i^{n+1})$ writes

$$ \eta(w_i^{n+1}) = \eta(w_i^n) - \frac{\Delta t}{\Delta x} \eta'(w_i^n) \Delta \mathcal{F}_i^n + \left(\frac{\Delta t}{\Delta x}\right)^2 \int_0^1 (1-s) \eta''(w_i^n - s \frac{\Delta t}{\Delta x} \Delta \mathcal{F}_i^n) (\Delta \mathcal{F}_i^n)^2 ds, \quad (16) $$

where $(\Delta \mathcal{F}_i^n)_{i \in \mathbb{Z}/\mathbb{N}}$ is defined in (12d). Equipped with this expansion, the proof of Theorem 3.1-(ii) consists at first to reformulate (16) to show the existence of a symmetric matrix $A_i^n(s) \in \mathcal{M}_{2,2}(\mathbb{R})$ such that

$$ \frac{\eta(w_i^{n+1}) - \eta(w_i^n)}{\Delta t} + \frac{\mathcal{G}(w_i^n, w_{i+1}^n) - \mathcal{G}(w_{i-1}^n, w_i^n)}{\Delta x} = \frac{1}{2\Delta x} \int_0^1 (\Delta_i^n)^\top A_i^n(s) \Delta_i^n ds, \quad (17) $$
where $\Delta_i^n = (\delta_{i-\frac{1}{2}}^n, \delta_{i+\frac{1}{2}}^n)^T \in \mathbb{R}^2$ and $\mathcal{G}(w_i^n, w_{i+1}^n)$ given by (14). To address such an issue, several lemma are now derived and the first one given below concerns a reformulation of $(\Delta \mathcal{F}_i^n)^2$.

**Lemma 5.1** (Reformulation for $(\Delta \mathcal{F}_i^n)^2$). Consider $(\delta_{i+\frac{1}{2}}^n, \Delta \mathcal{F}_i^n)_{i \in \mathbb{Z}/\mathbb{N}Z}$ defined in (12a), (12d) respectively and also $(a_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}Z}$, $(B_i^n)_{i \in \mathbb{Z}/\mathbb{N}Z} \in \mathbb{R} \times \mathcal{M}_{2,2}(\mathbb{R})$ such that

$$a_{i+\frac{1}{2}}^n = \int_0^1 f'(w_i^n + s\delta_{i+\frac{1}{2}}^n) \, ds, \qquad (18)$$

$$B_i^n = \begin{pmatrix} (a_{i-\frac{1}{2}}^n + \Lambda_{i-\frac{1}{2}}^n)^2 & (a_{i-\frac{1}{2}}^n + \Lambda_{i-\frac{1}{2}}^n)(a_{i+\frac{1}{2}}^n - \Lambda_{i+\frac{1}{2}}^n) \\ (a_{i-\frac{1}{2}}^n + \Lambda_{i-\frac{1}{2}}^n)(a_{i+\frac{1}{2}}^n - \Lambda_{i+\frac{1}{2}}^n) & (a_{i+\frac{1}{2}}^n - \Lambda_{i+\frac{1}{2}}^n)^2 \end{pmatrix}. \qquad (19)$$

The equality $(\Delta \mathcal{F}_i^n)^2 = \frac{1}{4}(\Delta_i^n)^T B_i^n \Delta_i^n$ is satisfied.

*Proof.* Since $w \mapsto f(w)$ is regular, the definitions of $\Delta \mathcal{F}_i^n, a_{i+\frac{1}{2}}^n$ given in (12d), (18) and a Taylor expansion in the neighborhood of $w_i^n$ give

$$\begin{aligned} (\Delta \mathcal{F}_i^n)^2 &= \frac{1}{4} (f(w_{i+1}^n) - f(w_{i-1}^n) - \Lambda_{i+\frac{1}{2}}^n \delta_{i+\frac{1}{2}}^n + \Lambda_{i-\frac{1}{2}}^n \delta_{i-\frac{1}{2}}^n)^2, \\ &= \frac{1}{4} (f(w_{i+1}^n) - f(w_i^n) + f(w_i^n) - f(w_{i-1}^n) - \Lambda_{i+\frac{1}{2}}^n \delta_{i+\frac{1}{2}}^n + \Lambda_{i-\frac{1}{2}}^n \delta_{i-\frac{1}{2}}^n)^2, \\ &= \frac{1}{4} (a_{i+\frac{1}{2}}^n \delta_{i+\frac{1}{2}}^n + a_{i-\frac{1}{2}}^n \delta_{i-\frac{1}{2}}^n - \Lambda_{i+\frac{1}{2}}^n \delta_{i+\frac{1}{2}}^n + \Lambda_{i-\frac{1}{2}}^n \delta_{i-\frac{1}{2}}^n)^2, \\ &= \frac{1}{4} ((a_{i+\frac{1}{2}}^n - \Lambda_{i+\frac{1}{2}}^n) \delta_{i+\frac{1}{2}}^n + (a_{i-\frac{1}{2}}^n + \Lambda_{i-\frac{1}{2}}^n) \delta_{i-\frac{1}{2}}^n)^2. \end{aligned}$$

Rewriting this last equality with the definition of $B_i^n$ given by (19), the statement of lemma follows. $\square$

Equipped with the above lemma, the goal now is to reformulate $-\eta'(w_i^n)\Delta \mathcal{F}_i^n$. According to the definitions of $\Delta \mathcal{F}_i^n, \mathcal{F}(w_i^n, w_{i+1}^n)$ given by (12d), (11), $-\eta'(w_i^n)\Delta \mathcal{F}_i^n$ also writes

$$-\eta'(w_i^n)\Delta \mathcal{F}_i^n = -\frac{1}{2}\eta'(w_i^n)(f(w_{i+1}^n) - f(w_i^n)) - \frac{1}{2}\eta'(w_i^n)(f(w_i^n) - f(w_{i-1}^n)) \\ + \frac{1}{2}\Lambda_{i+\frac{1}{2}}^n \eta'(w_i^n)\delta_{i+\frac{1}{2}}^n - \frac{1}{2}\Lambda_{i-\frac{1}{2}}^n \eta'(w_i^n)\delta_{i-\frac{1}{2}}^n. \qquad (20)$$

The three following lemma establish reformulations for (20) and the first one given below concerns $\Lambda_{i+\frac{1}{2}}^n \eta'(w_i^n)\delta_{i+\frac{1}{2}}^n/2 - \Lambda_{i-\frac{1}{2}}^n \eta'(w_i^n)\delta_{i-\frac{1}{2}}^n/2$.

**Lemma 5.2.** Let $\eta \in C^2(\mathbb{R}, \mathbb{R})$ and consider the sequence $(\delta_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}Z}$ given by (12a). For all sequences $(\Lambda_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}Z} \in \mathbb{R}$, the following equality holds:

$$\begin{aligned} & \frac{1}{2}\Lambda_{i+\frac{1}{2}}^n \eta'(w_i^n)\delta_{i+\frac{1}{2}}^n - \frac{1}{2}\Lambda_{i-\frac{1}{2}}^n \eta'(w_i^n)\delta_{i-\frac{1}{2}}^n \\ &= \frac{1}{2}\Lambda_{i+\frac{1}{2}}^n (\eta(w_{i+1}^n) - \eta(w_i^n)) - \frac{1}{2}\Lambda_{i-\frac{1}{2}}^n (\eta(w_i^n) - \eta(w_{i-1}^n)) \\ & \quad - \frac{1}{2}\sum_{\pm} \Lambda_{i \pm \frac{1}{2}}^n (\delta_{i \pm \frac{1}{2}}^n)^2 \int_0^1 (1-s)\eta''(w_i^n \pm s\delta_{i \pm \frac{1}{2}}^n) \, ds. \end{aligned}$$

*Proof.* Since $w \mapsto \eta(w)$ is regular, Taylor expansions in the neighborhood of $w_i^n$ read

$$\eta'(w_i^n)\delta_{i+\frac{1}{2}}^n = \eta(w_{i+1}^n) - \eta(w_i^n) - (\delta_{i+\frac{1}{2}}^n)^2 \int_0^1 (1-s)\eta''(w_i^n + s\delta_{i+\frac{1}{2}}^n) \, ds, \qquad (21)$$

$$-\eta'(w_i^n)\delta_{i-\frac{1}{2}}^n = \eta(w_{i-1}^n) - \eta(w_i^n) - (\delta_{i-\frac{1}{2}}^n)^2 \int_0^1 (1-s)\eta''(w_i^n - s\delta_{i-\frac{1}{2}}^n) \, ds. \qquad (22)$$

Computing the quantity $\Lambda_{i+\frac{1}{2}}^n (21)/2 + \Lambda_{i-\frac{1}{2}}^n (22)/2$, the expected reformulation is immediate that concludes the proof. $\square$
Now, it is necessary to rewrite $-\eta'(w_i^n)(f(w_{i+1}^n) - f(w_i^n))/2$ in the right-hand side of (20). In this regard, the aim is to use the relation $G'(w) = (\eta'f')(w)$ in order to couple Taylor expansions of both functions $w \mapsto f(w)$ and $w \mapsto G(w)$. The lemma detailed below illustrates this strategy.

**Lemma 5.3** (Reformulation for $-\eta'(w_i^n)(f(w_{i+1}^n) - f(w_i^n))/2$). Let $(\eta, G) \in (C^3(\mathbb{R}, \mathbb{R}))^2$ be an entropy-entropy flux pair and also consider $(\delta_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}}$ be given by (12a). Considering $s \mapsto \psi_{i+\frac{1}{2}}^+(s)$ defined in (12b), $-\eta'(w_i^n)(f(w_{i+1}^n) - f(w_i^n))/2 \in \mathbb{R}$ also reads

$$-\frac{1}{2}\eta'(w_i^n)(f(w_{i+1}^n) - f(w_i^n)) = -\frac{1}{2}(G(w_{i+1}^n) - G(w_i^n)) + \frac{1}{8}((\eta''f')(w_i^n) + (\eta''f')(w_{i+1}^n))(\delta_{i+\frac{1}{2}}^n)^2 + \frac{1}{2}(\delta_{i+\frac{1}{2}}^n)^3 \int_0^1 \psi_{i+\frac{1}{2}}^+(s) \, ds.$$

*Proof.* Since $w \mapsto (f, G)(w)$ is regular, three-order Taylor expansions in the neighborhood of $w_i^n$ write

$$\frac{1}{2}(f(w_{i+1}^n) - f(w_i^n)) = \frac{1}{2}f'(w_i^n)\delta_{i+\frac{1}{2}}^n + \frac{1}{4}f''(w_i^n)(\delta_{i+\frac{1}{2}}^n)^2 + \frac{1}{4}(\delta_{i+\frac{1}{2}}^n)^3 \int_0^1 (1-s)^2 f'''(w_i^n + s\delta_{i+\frac{1}{2}}^n)ds, \quad (23)$$

$$\frac{1}{2}(G(w_{i+1}^n) - G(w_i^n)) = \frac{1}{2}G'(w_i^n)\delta_{i+\frac{1}{2}}^n + \frac{1}{4}G''(w_i^n)(\delta_{i+\frac{1}{2}}^n)^2 + \frac{1}{4}(\delta_{i+\frac{1}{2}}^n)^3 \int_0^1 (1-s)^2 G'''(w_i^n + s\delta_{i+\frac{1}{2}}^n)ds. \quad (24)$$

As the entropy flux $w \mapsto G(w)$ is defined by $G' = \eta'f'$, it also satisfies $G'' = \eta''f' + \eta'f''$, that leads to express $-\eta'(w_i^n)(f(w_{i+1}^n) - f(w_i^n))/2$ in (23) as follows:

$$-\frac{1}{2}\eta'(w_i^n)(f(w_{i+1}^n) - f(w_i^n)) = -\frac{1}{2}G'(w_i^n)\delta_{i+\frac{1}{2}}^n - \frac{1}{4}(\eta'f'')(w_i^n)(\delta_{i+\frac{1}{2}}^n)^2 - \frac{1}{4}\eta'(w_i^n)(\delta_{i+\frac{1}{2}}^n)^3 \int_0^1 (1-s)^2 f'''(w_i^n + s\delta_{i+\frac{1}{2}}^n)ds, \\ = -\frac{1}{2}G'(w_i^n)\delta_{i+\frac{1}{2}}^n - \frac{1}{4}G''(w_i^n)(\delta_{i+\frac{1}{2}}^n)^2 + \frac{1}{4}(\eta''f')(w_i^n)(\delta_{i+\frac{1}{2}}^n)^2 - \frac{1}{4}\eta'(w_i^n)(\delta_{i+\frac{1}{2}}^n)^3 \int_0^1 (1-s)^2 f'''(w_i^n + s\delta_{i+\frac{1}{2}}^n)ds.$$

Using (24) in the above equality, the following reformulation holds:

$$-\frac{1}{2}\eta'(w_i^n)(f(w_{i+1}^n) - f(w_i^n)) = -\frac{1}{2}(G(w_{i+1}^n) - G(w_i^n)) + \frac{1}{4}(\eta''f')(w_i^n)(\delta_{i+\frac{1}{2}}^n)^2 + \frac{1}{4}(\delta_{i+\frac{1}{2}}^n)^3 \left( \int_0^1 (1-s)^2 \left( G'''(w_i^n + s\delta_{i+\frac{1}{2}}^n) - \eta'(w_i^n)f'''(w_i^n + s\delta_{i+\frac{1}{2}}^n) \right) ds \right). \quad (25)$$

However, the quantity $(\eta''f')(w_i^n)(\delta_{i+\frac{1}{2}}^n)^2/4$ also writes

$$\frac{1}{4}(\eta''f')(w_i^n)(\delta_{i+\frac{1}{2}}^n)^2 = \frac{1}{8}((\eta''f')(w_i^n) + (\eta''f')(w_{i+1}^n))(\delta_{i+\frac{1}{2}}^n)^2 - \frac{1}{8}((\eta''f')(w_{i+1}^n) - (\eta''f')(w_i^n))(\delta_{i+\frac{1}{2}}^n)^2, \\ = \frac{1}{8}((\eta''f')(w_i^n) + (\eta''f')(w_{i+1}^n))(\delta_{i+\frac{1}{2}}^n)^2 - \frac{1}{8}(\delta_{i+\frac{1}{2}}^n)^3 \int_0^1 (\eta''f')'(w_i^n + s\delta_{i+\frac{1}{2}}^n) ds.$$
Hence, associating the above equality to (25) then arguing the definition $s \mapsto \psi_{i+\frac{1}{2}}^{+}(s)$ given by (12b), the statement of lemma follows. $\square$

Following the computations of the above proof, it is now possible to achieve the reformulation of the right-hand side of (20) with the below corollary that concerns the quantity $-\eta'(w_i^n)(f(w_i^n) - f(w_{i-1}^n))/2$.

**Corollary 5.1** (Reformulation for $-\eta'(w_i^n)(f(w_i^n) - f(w_{i-1}^n))/2$). Let $(\eta, G) \in (C^3(\mathbb{R}, \mathbb{R}))^2$ be an entropy-entropy flux pair and also consider $(\delta_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}}$ given by (12a). Considering $s \mapsto \psi_{i+\frac{1}{2}}^-(s)$ given by (12c), $-\eta'(w_i^n)(f(w_i^n) - f(w_{i-1}^n))/2 \in \mathbb{R}$ also reads

$$ -\frac{1}{2}\eta'(w_i^n)(f(w_i^n) - f(w_{i-1}^n)) = -\frac{1}{2}(G(w_i^n) - G(w_{i-1}^n)) - \frac{1}{8}((\eta''f')(w_i^n) + (\eta''f')(w_{i-1}^n))(\delta_{i-\frac{1}{2}}^n)^2 + \frac{1}{2}(\delta_{i-\frac{1}{2}}^n)^3 \int_0^1 \psi_{i-\frac{1}{2}}^-(s) \, ds. $$

*Proof.* The proof is the same than once used for Lemma 5.3 (mutadis mutantdis). $\square$

According to Lemma 5.2, 5.3 and also Corollary 5.1, (20) henceforward reads

$$ -\eta'(w_i^n)\Delta\mathcal{F}_i^n = -(\mathcal{G}(w_i^n, w_{i+1}^n) - \mathcal{G}(w_{i-1}^n, w_i^n)) + \frac{1}{2}\sum_{\pm} (\delta_{i\pm\frac{1}{2}}^n)^2 \int_0^1 \left(\delta_{i\pm\frac{1}{2}}^n \psi_{i\pm\frac{1}{2}}^{\pm}(s) - \Lambda_{i\pm\frac{1}{2}}^n(1-s)\eta''(w_i^n \pm s\delta_{i\pm\frac{1}{2}}^n)\right) ds. \quad (26) $$

Using the above equality and Lemma 5.1, a formulation under the form (17) is now reachable for the second-order scheme (4), (11). This is detailed in the next lemma.

**Lemma 5.4** (Existence of the formulation (17)). Let $(\eta, G) \in (C^3(\mathbb{R}, \mathbb{R}))^2$ be an entropy-entropy flux pair and consider $(\delta_{i+\frac{1}{2}}^n, \Delta\mathcal{F}_i^n, B_i^n)_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}} \in (\mathbb{R})^2 \times \mathcal{M}_{2,2}(\mathbb{R})$ given by (12a) (12d), (19). If $s \mapsto A_i^n(s) \in \mathcal{M}_{2,2}(\mathbb{R})$ is defined as

$$ A_i^n(s) = C_i^n(s) + \frac{\Delta t}{2\Delta x}(1-s)\eta''(w_i^n - s\frac{\Delta t}{\Delta x}\Delta\mathcal{F}_i^n)B_i^n, \quad \forall s \in [0, 1], \quad (27) $$

where $s \mapsto C_i^n(s) \in \mathcal{M}_{2,2}(\mathbb{R})$ is such that

$$ C_i^n(s) = \begin{pmatrix} \delta_{i-\frac{1}{2}}^n \psi_{i-\frac{1}{2}}^-(s) & 0 \\ 0 & \delta_{i+\frac{1}{2}}^n \psi_{i+\frac{1}{2}}^+(s) \end{pmatrix} - (1-s) \begin{pmatrix} \Lambda_{i-\frac{1}{2}}^n \eta''(w_i^n - s\delta_{i-\frac{1}{2}}^n) & 0 \\ 0 & \Lambda_{i+\frac{1}{2}}^n \eta''(w_i^n + s\delta_{i+\frac{1}{2}}^n) \end{pmatrix}, \quad (28) $$

and where $\psi_{i+\frac{1}{2}}^{\pm}(s)$ are defined in (12b), (12c), then the second-order scheme (4), (11) verifies (17).

*Proof.* Since $A_i^n(s) \in \mathcal{M}_{2,2}(\mathbb{R})$ given by (27) obviously is symmetric, the proof will be completed as soon as the formulation (17) will be reached. In this regard, arguing Lemma 5.1, the Taylor expansion (16) is equivalent to

$$ \eta(w_i^{n+1}) = \eta(w_i^n) - \frac{\Delta t}{\Delta x}\eta'(w_i^n)\Delta\mathcal{F}_i^n + \left(\frac{\Delta t}{2\Delta x}\right)^2 \int_0^1 (1-s)\eta''(w_i^n - s\frac{\Delta t}{\Delta x}\Delta\mathcal{F}_i^n)(\Delta_i^n)^T B_i^n \Delta_i^n \, ds, $$
where $\Delta_i^n = (\delta_{i-\frac{1}{2}}^n, \delta_{i+\frac{1}{2}}^n)^\top \in \mathbb{R}^2$. Arguing (26) and the definition of $C_i^n(s)$ given in (28), the above Taylor expansion becomes

$$ \eta(w_i^{n+1}) - \eta(w_i^n) = -\frac{\Delta t}{\Delta x} (\mathcal{G}(w_i^n, w_{i+1}^n) - \mathcal{G}(w_{i-1}^n, w_i^n)) + \frac{\Delta t}{2\Delta x} \int_0^1 (\Delta_i^n)^\top C_i^n(s) \Delta_i^n ds \\ + \left(\frac{\Delta t}{2\Delta x}\right)^2 \int_0^1 (1-s)\eta''(w_i^n - s\frac{\Delta t}{\Delta x}\Delta_i^n)(\Delta_i^n)^\top B_i^n \Delta_i^n ds. $$

Since (17) can be deduced from $A_i^n(s)$ defined in by (27) and also from a reorganization of the terms of this last equality, that achieves the proof. $\square$

Equipped with (17), the proof of the discrete entropy inequality (5) is now equivalent to show that $\int_0^1 A_i^n(s) ds \in \mathcal{M}_{2,2}(\mathbb{R})$ has two negative eigenvalues for all $i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}$. According to (27), $\int_0^1 A_i^n(s) ds$ is symmetric in $\mathcal{M}_{2,2}(\mathbb{R})$ and it depends on $w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t$. As a consequence, it is sufficient to show the existence of $\Delta t > 0$ such that the trace (resp. determinant) of $\int_0^1 A_i^n(s) ds$ is negative (resp. positive) for all $w_i^n, \delta_{i\pm\frac{1}{2}}^n$. In this regard, two additional functions are now introduced as follows:

$$ \Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t) := \mathrm{Tr}(\int_0^1 A_i^n(s) \mathrm{d}s), \quad \Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t) := \det(\int_0^1 A_i^n(s) \mathrm{d}s). \quad (29) $$

Considering the above functions, the proof of Theorem 3.1-(ii) will be achieved as soon as the existence of $\Delta t > 0$ such that

$$ \Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t) < 0, \quad \Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t) > 0, \quad \forall (w_i^n, \delta_{i\pm\frac{1}{2}}^n) \in \mathbb{R}^3, \quad (30) $$

will be established. However, since $w \mapsto (f, \eta, G)(w)$ are assumed belongs to $(C^3(\mathbb{R}, \mathbb{R}))^3$, standard continuity arguments hold in the regimes $\Delta t \to 0$. As a consequence, such arguments are now used for the proof of Theorem 3.1-(ii).

*Proof of Theorem 3.1-(ii).* Considering $i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}$ fixed and arguing the definitions of $\Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t)$, $A_i^n(s)$ given by (29), (27), $\Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, 0)$ writes

$$ \begin{aligned} \Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, 0) &= \mathrm{Tr}(\int_0^1 A_i^n(s) \mathrm{d}s)|_{\Delta t=0}, \\ &= \mathrm{Tr}(\int_0^1 C_i^n(s) \mathrm{d}s), \\ &= \sum_{\pm} \int_0^1 \left( \delta_{i\pm\frac{1}{2}}^n \psi_{i\pm\frac{1}{2}}^{\pm}(s) - \Lambda_{i\pm\frac{1}{2}}^n (1-s)\eta''(w_i^n \pm s\delta_{i\pm\frac{1}{2}}^n) \right) \mathrm{d}s, \end{aligned} $$

where $C_i^n(s)$ is given by (28). Using (11b) to expand $\Lambda_{i\pm\frac{1}{2}}^n$ and since $x = \mathrm{sign}(x)|x|$ for all $x \in \mathbb{R}$, the above formulation becomes

$$ \Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, 0) = \sum_{\pm} |\delta_{i\pm\frac{1}{2}}^n| \int_0^1 \left( \mathrm{sign}(\delta_{i\pm\frac{1}{2}}^n) \psi_{i\pm\frac{1}{2}}^{\pm}(s) - \gamma_{i\pm\frac{1}{2}}^n (1-s) \eta''(w_i^n \pm s\delta_{i\pm\frac{1}{2}}^n) \right) \mathrm{d}s \\ - \Delta x \sum_{\pm} \sigma_{i\pm\frac{1}{2}}^n \int_0^1 (1-s) \eta''(w_i^n \pm s\delta_{i\pm\frac{1}{2}}^n) \mathrm{d}s. $$

Using the constraint (13) that restricts $\gamma_{i+\frac{1}{2}}^n > 0$ and as $\sigma_{i\pm\frac{1}{2}}^n, \eta'' > 0$, the right-hand side of the above equality is upper-bounded as follows:

$$ \Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, 0) \leq -\Delta x \sum_{\pm} \sigma_{i\pm\frac{1}{2}}^n \int_0^1 (1-s) \eta''(w_i^n \pm s\delta_{i\pm\frac{1}{2}}^n) \mathrm{d}s < 0. $$
This last inequality coupled to a standard continuity argument leads to the existence of $\Delta t_i > 0$ such that $\Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t_i) < 0$. Since such an existence holds for all $i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}$, the time step $0 < \Delta t = \min_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}} \Delta t_i$ ensures $\Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t) < 0$. The rest of the proof now deals with $\Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t) > 0$.

For $i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}$ fixed, using the quantity $\Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t)$, $A_i^n(s)$ defined in (29), (27), $\Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, 0)$ writes

$$ \begin{aligned} \Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, 0) &= \det\left(\int_0^1 A_i^n(s) \, ds\right)|_{\Delta t=0}, \\ &= \prod_{\pm} \int_0^1 \left(\delta_{i\pm\frac{1}{2}}^n \psi_{i\pm\frac{1}{2}}^{\pm}(s) - \Lambda_{i\pm\frac{1}{2}}^n (1-s)\eta''(w_i^n \pm s\delta_{i\pm\frac{1}{2}}^n)\right) \, ds. \end{aligned} \quad (31) $$

Now, using (11b) to reformulate $\Lambda_{i\pm\frac{1}{2}}^n$ and denoting $r_{i\pm\frac{1}{2}}^n$ such that $-r_{i\pm\frac{1}{2}}^n = \int_0^1 (\text{sign}(\delta_{i\pm\frac{1}{2}}^n)\psi_{i\pm\frac{1}{2}}^{\pm}(s) - \gamma_{i\pm\frac{1}{2}}^n(1-s)\eta''(w_i^n \pm s\delta_{i\pm\frac{1}{2}}^n)) \, ds \le 0$, which is non-positive in virtue of the constraint (13), (31) also writes

$$ \begin{aligned} \Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, 0) &= \prod_{\pm} \left( -|\delta_{i\pm\frac{1}{2}}^n| r_{i\pm\frac{1}{2}}^n - \Delta x \sigma_{i\pm\frac{1}{2}}^n \int_0^1 (1-s) \eta''(w_i^n \pm s \delta_{i\pm\frac{1}{2}}^n) \, ds \right), \\ &= \prod_{\pm} \left( |\delta_{i\pm\frac{1}{2}}^n| r_{i\pm\frac{1}{2}}^n + \Delta x \sigma_{i\pm\frac{1}{2}}^n \int_0^1 (1-s) \eta''(w_i^n \pm s \delta_{i\pm\frac{1}{2}}^n) \, ds \right), \\ &= \prod_{\pm} |\delta_{i\pm\frac{1}{2}}^n| r_{i\pm\frac{1}{2}}^n + \Delta x^2 \prod_{\pm} \sigma_{i\pm\frac{1}{2}}^n \int_0^1 (1-s) \eta''(w_i^n \pm s \delta_{i\pm\frac{1}{2}}^n) \, ds \\ &\quad + \Delta x \sum_{\pm} r_{i\pm\frac{1}{2}}^n |\delta_{i\pm\frac{1}{2}}^n| \sigma_{i\mp\frac{1}{2}}^n \int_0^1 (1-s) \eta''(w_i^n \mp s \delta_{i\mp\frac{1}{2}}^n) \, ds, \end{aligned} $$

where both symbols $\pm$ and $\mp$ mean $+$ or $-$ but they are such that $\mp = -$ if $\pm = +$ and conversely. As $r_{i\pm\frac{1}{2}}^n \ge 0$ and $\sigma_{i\pm\frac{1}{2}}^n, \eta'' > 0$, $\Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, 0)$ in the above formulation is lower-bounded as follows:

$$ \Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, 0) \ge \Delta x^2 \prod_{\pm} \sigma_{i\pm\frac{1}{2}}^n \int_0^1 (1-s) \eta''(w_i^n \pm s \delta_{i\pm\frac{1}{2}}^n) \, ds > 0. $$

This last inequality leads to the existence of $\Delta t_i > 0$ such that $\Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t_i) > 0$. As a consequence, $0 < \Delta t = \min_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}} \Delta t_i$ enforces $\Psi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t) > 0$.

Therefore, under the constraint (13) and $\sigma_{i\pm\frac{1}{2}}^n, \eta'' > 0$, there exists $\Delta t > 0$ such that the symmetric matrix $\int_0^1 A_i^n(s) \, ds \in \mathcal{M}_{2,2}(\mathbb{R})$ given by (27) is also definite negative. Thus, using the formulation (17), the second-order scheme (4), (11) satisfies (5) for all strictly convex entropy-entropy flux pair $(\eta, G) \in (C^3(\mathbb{R}, \mathbb{R}))^2$, that concludes the proof. $\square$

The above proof stands for a CFL-like restriction for $\Delta t > 0$. This restriction is implicitly given but it is at least under the form $\frac{\Delta t}{\Delta x^2} \sim \mathcal{O}(1)$. Such a CFL-like condition is highlighted by $\Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t)$ in the regimes $\delta_{i\pm\frac{1}{2}}^n \to 0$ that read

$$ \lim_{\delta_{i\pm\frac{1}{2}}^n \to 0} \Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t) = \frac{\Delta x}{2} \eta''(w_i^n) \left( -\sum_{\pm} \sigma_{i\pm\frac{1}{2}}^n + \frac{\Delta t}{2\Delta x^2} \sum_{\pm} (f'(w_i^n) \mp \Delta x \sigma_{i\pm\frac{1}{2}}^n)^2 \right). $$

As a consequence, when $\delta_{i\pm\frac{1}{2}}^n \to 0$, the inequality $\Phi(w_i^n, \delta_{i\pm\frac{1}{2}}^n, \Delta t) < 0$ is enforced if $\frac{\Delta t}{\Delta x^2} \sim \mathcal{O}(1)$. Such parabolic CFL-like conditions have already been exhibited for high-order numerical schemes in [12, 13, 47] and the practical choice of $\frac{\Delta t}{\Delta x^2} > 0$ for (4), (11) is discussed in the next section that concerns the numerical experiments.
# 6 Numerical experiments

## 6.1 Practical implementation

This section concerns the practical implementation of the second-order scheme (4), (11) focusing on the enforcement of a single discrete entropy inequality (5). In this regard, considering a given strictly convex entropy-entropy flux pair $(\eta, G) \in (C^3(\mathbb{R}, \mathbb{R}))^2$, the scheme (4), (11) needs the definition of the triplet $(\sigma_{i+\frac{1}{2}}^n, \gamma_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}}$, $\frac{\Delta t}{\Delta x^2} > 0$ in a such way (5) holds. The selection procedure of this triplet starts from $(\eta, G, f)$ to compute $(\psi_{i+\frac{1}{2}}^{\pm}(s))_{(i,s) \in \mathbb{Z}/\mathbb{N}\mathbb{Z} \times [0,1]}$ according to (12b), (12c). Then, for a given $\sigma_{i+\frac{1}{2}}^n > 0$, $\gamma_{i+\frac{1}{2}}^n > 0$ has to enforces (13). Finally, $\frac{\Delta t}{\Delta x^2} > 0$ has to be selected in a such way $\int_0^1 A_i^n(s) ds \in \mathcal{M}_{2,2}(\mathbb{R})$ given by (27) is definite negative. Theorem 3.1 gives explicit constraints for $(\gamma_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}}$ but the restriction for $\frac{\Delta t}{\Delta x^2} > 0$ remains implicit. As a consequence, $\frac{\Delta t}{\Delta x^2} > 0$ is selected according to the following heuristic:

$$ \Delta t = \text{CFL} \Delta x \min_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}} \left( \gamma_{i+\frac{1}{2}}^n |w_{i+1}^n - w_i^n| + \sigma_{i+\frac{1}{2}}^n \Delta x \right), \quad (32) $$

where CFL $\in (0, 1]$. In the sequel, the above heuristic is also checked with the quantity

$$ \frac{\Delta S}{\Delta t} := \frac{\eta(w_i^{n+1}) - \eta(w_i^n)}{\Delta t} + \frac{\mathcal{G}(w_i^n, w_{i+1}^n) - \mathcal{G}(w_{i-1}^n, w_i^n)}{\Delta x}. \quad (33) $$

The above quantity $\Delta S/\Delta t \in \mathbb{R}$ obviously is equivalent to (5) and the choices of $(\sigma_{i+\frac{1}{2}}^n, \gamma_{i+\frac{1}{2}}^n)_{i \in \mathbb{Z}/\mathbb{N}\mathbb{Z}}$, $\frac{\Delta t}{\Delta x^2} > 0$ have to ensure $\Delta S/\Delta t \le 0$. In order to illustrate the above implementation, the first numerical experiment concerns the linear transport equation.

## 6.2 Linear transport equation

This section is devoted to the advection problem endowed by the quartic entropy defined by

$$ f(w) = aw, \quad a = 1, \quad \eta(w) = w^2/2 + w^4/4, \quad G(w) = a(w^2/2 + w^4/4). \quad (34) $$

According to the above formulation, (13) reads $\gamma_{i+\frac{1}{2}}^n \ge \max \left( \frac{2|w_{i+1}^n|}{(w_i^n + w_{i+1}^n)^2 + (w_i^n)^2 + 2}, \frac{2|w_i^n|}{(w_i^n + w_{i+1}^n)^2 + (w_{i+1}^n)^2 + 2} \right)$. Since a straightforward computation shows that the right hand-side of this inequality is bounded by 1, it is sufficient to satisfy $\gamma_{i+\frac{1}{2}}^n \ge 1$. As a consequence and for the sake of simplicity, $\sigma_{i+\frac{1}{2}}^n, \gamma_{i+\frac{1}{2}}^n$ are selected according to

$$ \sigma_{i+\frac{1}{2}}^n = \gamma_{i+\frac{1}{2}}^n := \alpha \ge 1. \quad (35) $$

The above definition ensures (13) and two values for $\alpha \ge 1$ are tested in the sequel.

At first, the smooth initial data $w_0(x) = \sin(\pi x)$ is considered over the periodic domain $I = [-1, 1]$. The final time is $t = 0.1$ and the time step $\Delta t > 0$ is selected according to (32) with CFL = 0.9. Figure 2 displays the numerical results and Table 2 reports the order of convergence.

<figure><img src="image_4.png" /><figcaption>Figure 2: Second-order approximation and $\Delta S/\Delta t$ defined in (33) for $f(w) = aw$ endowed by $w_0(x) = \sin(\pi x)$ on a mesh made of $N = 400$ cells.</figcaption></figure>
<table><caption>Table 2: Errors and order evaluations for the second-order scheme (4), (11) for $f(w) = aw$ endowed by $w_0(x) = \sin(\pi x)$.</caption><thead><tr><th colspan="7">α = 1</th></tr><tr><th>cells</th><th>L<sup>1</sup></th><th>order</th><th>L<sup>2</sup></th><th>order</th><th>L<sup>∞</sup></th><th>order</th></tr></thead><tbody><tr><td>100</td><td>8.2E-04</td><td>-</td><td>6.6E-04</td><td>-</td><td>7.8E-04</td><td>-</td></tr><tr><td>200</td><td>2.1E-04</td><td>2.0</td><td>1.7E-04</td><td>2.0</td><td>2.0E-04</td><td>2.0</td></tr><tr><td>400</td><td>5.2E-05</td><td>2.0</td><td>4.2E-05</td><td>2.0</td><td>4.9E-05</td><td>2.0</td></tr><tr><td>800</td><td>1.3E-05</td><td>2.0</td><td>1.0E-05</td><td>2.0</td><td>1.2E-05</td><td>2.0</td></tr><tr><td>1600</td><td>3.3E-06</td><td>2.0</td><td>2.6E-06</td><td>2.0</td><td>3.1E-06</td><td>2.0</td></tr></tbody></table><table><thead><tr><th colspan="7">α = 10.0</th></tr><tr><th>cells</th><th>L<sup>1</sup></th><th>order</th><th>L<sup>2</sup></th><th>order</th><th>L<sup>∞</sup></th><th>order</th></tr></thead><tbody><tr><td>100</td><td>8.0E-03</td><td>-</td><td>6.2E-03</td><td>-</td><td>6.3E-03</td><td>-</td></tr><tr><td>200</td><td>2.0E-03</td><td>2.0</td><td>1.6E-03</td><td>2.0</td><td>1.6E-03</td><td>2.0</td></tr><tr><td>400</td><td>5.1E-04</td><td>2.0</td><td>4.0E-04</td><td>2.0</td><td>4.0E-04</td><td>2.0</td></tr><tr><td>800</td><td>1.3E-04</td><td>2.0</td><td>9.9E-05</td><td>2.0</td><td>1.0E-04</td><td>2.0</td></tr><tr><td>1600</td><td>3.2E-05</td><td>2.0</td><td>2.5E-05</td><td>2.0</td><td>2.5E-05</td><td>2.0</td></tr></tbody></table>

A good agreement to the exact solution is observed in Figure 2. The inequality $\Delta S/\Delta t \le 0$ is also observed that means the enforcement of the discrete entropy inequality for the couple $(\eta, G)$ defined in (34). The numerical errors in Table 2 are greater for $\alpha = 10$ but for both values of $\alpha$, the second-order accuracy is exactly observed. Now, the discontinuous initial data $\chi_{|x|\le 0.25}(x)$ is considered over the periodic domain $[-0.5, 0.5]$. The final time is $t = 0.1$ and the time step $\Delta t > 0$ is selected according to (32) with CFL = 0.9. Figure 3 displays the numerical results and Table 3 reports the order of convergence.

<figure><img src="image_5.png" /><figcaption>Figure 3: Second-order approximation and $\Delta S/\Delta t$ defined in (33) for $f(w) = aw$ endowed with $w_0(x) = \chi_{|x|\le 0.25}(x)$ on a mesh made of $N = 400$ cells.</figcaption></figure>

Both values of $\alpha$ generate spurious oscillations but their amplitudes genuinely are reduced when $\alpha = 10$. In this case, the numerical results are very diffusive. For both values $\alpha = 1, \alpha = 10$, the inequality $\Delta S/\Delta t \le 0$ is checked that means the enforcement of the discrete entropy inequality (5) for the couple $(\eta, G)$ defined in (34).

<table><thead><tr><th colspan="7">α = 1</th></tr><tr><th>cells</th><th>L<sup>1</sup></th><th>order</th><th>L<sup>2</sup></th><th>order</th><th>L<sup>∞</sup></th><th>order</th></tr></thead><tbody><tr><td>100</td><td>3.8E-02</td><td>-</td><td>9.7E-02</td><td>-</td><td>4.6E-01</td><td>-</td></tr><tr><td>200</td><td>2.5E-02</td><td>0.6</td><td>7.8E-02</td><td>0.3</td><td>4.9E-01</td><td>0.1</td></tr><tr><td>400</td><td>1.6E-02</td><td>0.6</td><td>6.2E-02</td><td>0.3</td><td>5.0E-01</td><td>0.0</td></tr><tr><td>800</td><td>1.0E-02</td><td>0.6</td><td>4.9E-02</td><td>0.3</td><td>5.2E-01</td><td>0.0</td></tr><tr><td>1600</td><td>6.7E-03</td><td>0.6</td><td>3.9E-02</td><td>0.3</td><td>5.3E-01</td><td>0.0</td></tr></tbody></table>
<table><caption>Table 3: Errors and order evaluations for the second-order scheme (4), (11) for $f(w) = aw$ endowed by $w_0(x) = \chi_{|x|\le 0.25}(x)$.</caption><thead><tr><th colspan="7">α = 10</th></tr><tr><th>cells</th><th>L<sup>1</sup></th><th>order</th><th>L<sup>2</sup></th><th>order</th><th>L<sup>∞</sup></th><th>order</th></tr></thead><tbody><tr><td>100</td><td>5.9E-02</td><td>-</td><td>1.3E-01</td><td>-</td><td>4.6E-01</td><td>-</td></tr><tr><td>200</td><td>3.7E-02</td><td>0.7</td><td>1.1E-01</td><td>0.3</td><td>4.7E-01</td><td>0.0</td></tr><tr><td>400</td><td>2.3E-02</td><td>0.7</td><td>8.5E-02</td><td>0.3</td><td>4.8E-01</td><td>0.0</td></tr><tr><td>800</td><td>1.5E-02</td><td>0.7</td><td>6.7E-02</td><td>0.3</td><td>4.8E-01</td><td>0.0</td></tr><tr><td>1600</td><td>9.3E-03</td><td>0.7</td><td>5.3E-02</td><td>0.3</td><td>4.9E-01</td><td>0.0</td></tr></tbody></table>

## 6.3 The Burgers equation

This section is devoted to the Burgers equation endowed by the quadratic entropy that reads

$$f(w) = w^2/2, \quad \eta(w) = w^2/2, \quad G(w) = w^3/3. \qquad (36)$$

In this case, (13) reads $\gamma_{i+\frac{1}{2}}^n \ge \frac{1}{6} \max(|1 - 2w_i^n|, |1 - 2w_{i+1}^n|)$. Since the right hand-side of this inequality is not bounded with respect to the sequence $(w_i^n)$, the opportunity is taken to test the choice (35) for the non-linear Burgers equation (36). The heuristic choice (35) is obviously checked with the criteria $\Delta S/\Delta t$ defined in (33).

In the first numerical test, the smooth initial data $w_0(x) = \sin(\pi x)$ is considered over the periodic domain $I = [-1, 1]$. With a final time small enough, here given by $t = 0.05$, the exact solution remains smooth so that the order of accuracy can be assessed. The time step $\Delta t > 0$ is selected according to (32) with CFL = 0.9. Figure 4 displays the numerical results and Table 4 reports the order of convergence.

<figure><img src="image_6.png" /><figcaption>Figure 4: Second-order approximation and $\Delta S/\Delta t$ defined in (33) for $f(w) = w^2/2$ endowed by $w_0(x) = \sin(\pi x)$ on a mesh made of $N = 400$ cells.</figcaption></figure>

<table><thead><tr><th colspan="7">α = 1</th></tr><tr><th>cells</th><th>L<sup>1</sup></th><th>order</th><th>L<sup>2</sup></th><th>order</th><th>L<sup>∞</sup></th><th>order</th></tr></thead><tbody><tr><td>100</td><td>5.5E-04</td><td>-</td><td>5.4E-04</td><td>-</td><td>8.2E-04</td><td>-</td></tr><tr><td>200</td><td>1.4E-04</td><td>2.0</td><td>1.4E-04</td><td>2.0</td><td>2.1E-04</td><td>2.0</td></tr><tr><td>400</td><td>3.5E-05</td><td>2.0</td><td>3.4E-05</td><td>2.0</td><td>5.2E-05</td><td>2.0</td></tr><tr><td>800</td><td>8.6E-06</td><td>2.0</td><td>8.5E-06</td><td>2.0</td><td>1.3E-05</td><td>2.0</td></tr><tr><td>1600</td><td>2.2E-06</td><td>2.0</td><td>2.1E-06</td><td>2.0</td><td>3.2E-06</td><td>2.0</td></tr></tbody></table>
<table><caption>Table 4: Errors and order evaluations for the second-order scheme (4), (11) for $f(w) = w^2/2$ endowed by $w_0(x) = \sin(\pi x)$.</caption><thead><tr><th colspan="7">α = 10.0</th></tr><tr><th>cells</th><th>L<sup>1</sup></th><th>order</th><th>L<sup>2</sup></th><th>order</th><th>L<sup>∞</sup></th><th>order</th></tr></thead><tbody><tr><td>100</td><td>5.3E-03</td><td>-</td><td>4.3E-03</td><td>-</td><td>5.6E-03</td><td>-</td></tr><tr><td>200</td><td>1.3E-03</td><td>2.0</td><td>1.1E-03</td><td>2.0</td><td>1.4E-03</td><td>2.0</td></tr><tr><td>400</td><td>3.3E-04</td><td>2.0</td><td>2.8E-04</td><td>2.0</td><td>3.6E-04</td><td>2.0</td></tr><tr><td>800</td><td>8.3E-05</td><td>2.0</td><td>6.9E-05</td><td>2.0</td><td>9.1E-05</td><td>2.0</td></tr><tr><td>1600</td><td>2.1E-05</td><td>2.0</td><td>1.7E-05</td><td>2.0</td><td>2.3E-05</td><td>2.0</td></tr></tbody></table>

A good agreement to the exact solution is observed in Figure 4. The inequality $\Delta S/\Delta t \le 0$ is also satisfied that implies the enforcement of the discrete entropy inequality for the couple $(\eta, G)$ defined in (36). For $\alpha = 1$ and $\alpha = 10$, the second-order accuracy is exactly once again observed. The case of discontinuous solution is now investigated considering the initial data $\chi_{|x|\le 0.25}(x)$ over the periodic domain $[-0.5, 0.5]$. The final time is $t = 0.1$ and the time step $\Delta t > 0$ is selected according to (32) with CFL = 0.9. Figure 5 displays the numerical results and Table 5 reports the order of convergence.

<figure><img src="image_7.png" /><figcaption>Figure 5: Second-order approximation and $\Delta S/\Delta t$ defined in (33) for $f(w) = w^2/2$ endowed with $w_0(x) = \chi_{|x|\le 0.25}(x)$ on a mesh made of $N = 400$ cells.</figcaption></figure>

Both values of $\alpha$ generate spurious oscillations but their amplitudes once again are reduced when $\alpha = 10$. In this case, the numerical results are very diffusive. For both values $\alpha = 1$, $\alpha = 10$, the inequality $\Delta S/\Delta t \le 0$ is checked so the discrete entropy inequality (5) is satisfied for the couple $(\eta, G)$ defined in (36).

<table><caption>Table 5: Errors and order evaluations for the second-order scheme (4), (11) for $f(w) = w^2/2$ endowed by $w_0(x) = \chi_{|x|\le 0.25}(x)$.</caption><thead><tr><th colspan="7">α = 1</th></tr><tr><th>cells</th><th>L<sup>1</sup></th><th>order</th><th>L<sup>2</sup></th><th>order</th><th>L<sup>∞</sup></th><th>order</th></tr></thead><tbody><tr><td>100</td><td>1.7E-02</td><td>-</td><td>4.7E-02</td><td>-</td><td>3.2E-01</td><td>-</td></tr><tr><td>200</td><td>9.4E-03</td><td>0.8</td><td>3.3E-02</td><td>0.5</td><td>3.2E-01</td><td>0.0</td></tr><tr><td>400</td><td>5.2E-03</td><td>0.9</td><td>2.3E-02</td><td>0.5</td><td>3.2E-01</td><td>0.0</td></tr><tr><td>800</td><td>2.8E-03</td><td>0.9</td><td>1.6E-02</td><td>0.5</td><td>3.2E-01</td><td>0.0</td></tr><tr><td>1600</td><td>1.5E-03</td><td>0.9</td><td>1.1E-02</td><td>0.5</td><td>3.1E-01</td><td>0.0</td></tr></tbody></table><table><thead><tr><th colspan="7">α = 10</th></tr><tr><th>cells</th><th>L<sup>1</sup></th><th>order</th><th>L<sup>2</sup></th><th>order</th><th>L<sup>∞</sup></th><th>order</th></tr></thead><tbody><tr><td>100</td><td>4.0E-02</td><td>-</td><td>9.3E-02</td><td>-</td><td>4.3E-01</td><td>-</td></tr><tr><td>200</td><td>2.1E-02</td><td>0.9</td><td>6.3E-02</td><td>0.6</td><td>4.3E-01</td><td>0.0</td></tr><tr><td>400</td><td>1.1E-02</td><td>1.0</td><td>4.2E-02</td><td>0.6</td><td>4.3E-01</td><td>0.0</td></tr><tr><td>800</td><td>5.6E-03</td><td>1.0</td><td>2.9E-02</td><td>0.6</td><td>4.3E-01</td><td>0.0</td></tr><tr><td>1600</td><td>2.8E-03</td><td>1.0</td><td>2.0E-02</td><td>0.5</td><td>4.3E-01</td><td>0.0</td></tr></tbody></table>
References

[1] M. BADSI, C. BERTHON, AND L. MARTAUD, *A family of second-order dissipative finite volume schemes for hyperbolic systems of conservation laws*, The SMAI Journal of computational mathematics, 9 (2023), pp. 31–60.

[2] C. BERTHON, *Stability of the MUSCL schemes for the Euler equations*, Comm. Math. Sci. 3 (2005), pp. 133–158.

[3] F. BOUCHUT, *Construction of BGK models with a family of kinetic entropies for a given system of conservation laws*, Journal of statistical physics, 95 (1999), pp. 113–170.

[4] F. BOUCHUT, *Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources*, Frontiers in Mathematics, Birkhäuser Verlag, Basel, 2004.

[5] F. BOUCHUT, C. BOURDARIAS, AND B. PERTHAME, *A MUSCL method satisfying all the numerical entropy inequalities*, Mathematics of Computation, 65 (1996), pp. 1439–1462.

[6] F. BOUCHUT AND T. M. DE LUNA, *An entropy satisfying scheme for two-layer shallow water equations with uncoupled treatment*, ESAIM: Mathematical Modelling and Numerical Analysis, 42 (2008), pp. 683–698.

[7] P. CHANDRASHEKAR, *Kinetic energy preserving and entropy stable finite volume schemes for compressible Euler and Navier-Stokes equations*, Communications in Computational Physics, 14 (2012).

[8] X. CHENG, *A fourth order entropy stable scheme for hyperbolic conservation laws*, Entropy (Basel, Switzerland), 21 (2019).

[9] B. COCKBURN, F. COQUEL, AND P. G. LEFLOCH, *Convergence of the finite volume method for multidimensional conservation laws*, SIAM Journal on Numerical Analysis, 32 (1995), pp. 687–705.

[10] F. COQUEL, J.-M. HÉRARD, AND K. SALEH, *A splitting method for the isentropic baer-nunziato two-phase flow model*, in ESAIM: Proceedings, vol. 38, EDP Sciences, 2012, pp. 241–256.

[11] F. COQUEL AND B. PERTHAME, *Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics*, SIAM Journal on Numerical Analysis, 35 (1998), pp. 2223–2249.

[12] E. DERIAZ, *Stability conditions for the numerical solution of convection-dominated problems with skew-symmetric discretizations*, SIAM Journal on Numerical Analysis, 50 (2012), pp. 1058–1085.

[13] E. DERIAZ AND P. HALDENWANG, *Non-linear CFL conditions issued from the von Neumann stability analysis for the transport equation*, Journal of Scientific Computing, 85 (2020), p. 5.

[14] V. DESVEAUX AND C. BERTHON, *An entropy preserving mood scheme for the euler equations*, International Journal On Finite Volumes, 11 (2014), pp. 1–39.

[15] A. DURAN AND F. MARCHE, *Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms*, Comput. & Fluids, 101 (2014), pp. 88–104.

[16] A. DURAN, J. VILA, AND R. BARAILLE, *Semi-implicit staggered mesh scheme for the multi-layer shallow water system*, C. R. Acad. Sci. Paris., (2017).

[17] U. FJORDHOLM, S. MISHRA, AND E. TADMOR, *Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws*, SIAM Journal on Numerical Analysis, 50 (2012), pp. 544–573.

[18] E. GODLEWSKI AND P.-A. RAVIART, *Hyperbolic systems of conservation laws*, vol. 3/4 of Mathématiques & Applications (Paris) [Mathematics and Applications], Ellipses, Paris, 1991.
[19] E. GODLEWSKI AND P.-A. RAVIART, *Numerical approximation of hyperbolic systems of conservation laws*, vol. 118 of Applied Mathematical Sciences, Springer-Verlag, New York, 1996.

[20] S. K. GODUNOV, *A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics*, Mat. Sb. (N.S.), 47 (1959), pp. 271–306.

[21] A. HARTEN, P. LAX, AND B. VAN LEER, *On upstream differencing and Godunov-type schemes for hyperbolic conservation laws*, SIAM review, 25 (1983), pp. 35–61.

[22] A. HARTEN AND P. D. LAX, *A random choice finite difference scheme for hyperbolic conservation laws*, SIAM Journal on Numerical Analysis, 18 (1981), pp. 289–315.

[23] R. HERBIN, W. KHERIJI, AND J. LATCHÉ, *On some implicit and semi-implicit staggered schemes for the shallow water and euler equations*, Mathematical Modelling and Numerical Analysis, 48 (2014), pp. 1807–1857.

[24] R. HERBIN, J. LATCHÉ, AND T. NGUYEN, *Consistent explicit staggered schemes for compressible flows part i: the barotropic euler equations*, (2013).

[25] A. HILTEBRAND AND S. MISHRA, *Entropy stable shock capturing space-time discontinuous galerkin schemes for systems of conservation laws*, Numerische Mathematik, 126 (2014), pp. 103–151.

[26] F. ISMAIL AND P. L. ROE, *Affordable, entropy-consistent Euler flux functions II: Entropy production at shocks*, Journal of Computational Physics, 228 (2009), pp. 5410–5436.

[27] B. KHOBALATTE AND B. PERTHAME, *Maximum principle on the entropy and second-order kinetic schemes*, Math. Comp., 62 (1994), pp. 119–131.

[28] D. KUZMIN, H. HAJDUK, AND A. RUPP, *Limiter-based entropy stabilization of semi-discrete and fully discrete schemes for nonlinear hyperbolic problems*, 2021, https://arxiv.org/abs/2107.11283.

[29] P. LAX, *Shock waves and entropy*, in Contributions to nonlinear functional analysis (Proc. Sympos., Math. Res. Center, Univ. Wisconsin, Madison, Wis., 1971), Academic Press, New York, 1971, pp. 603–634.

[30] P. LAX, *Hyperbolic systems of conservation laws and the mathematical theory of shock waves*, Society for Industrial and Applied Mathematics, Philadelphia, Pa., 1973. Conference Board of the Mathematical Sciences Regional Conference Series in Applied Mathematics, No. 11.

[31] P. LAX AND B. WENDROFF, *Systems of conservation laws*, Comm. Pure Appl. Math., 13 (1960), pp. 217–237.

[32] P. G. LEFLOCH, *Hyperbolic systems of conservation laws*, Lectures in Mathematics ETH Zürich, Birkhäuser Verlag, Basel, 2002. The theory of classical and nonclassical shock waves.

[33] S. NOELLE, Y. XING, AND C. W. SHU, *High-order well-balanced finite volume WENO schemes for shallow water equation with moving water*, J. Comput. Phys., 226 (2007), pp. 29–58.

[34] B. PERTHAME AND Y. QIU, *A variant of van Leer's method for multidimensional systems of conservation laws*, Journal of Computational Physics, 112 (1994), pp. 370–381.

[35] B. PERTHAME AND C. W. SHU, *On positivity preserving finite volume schemes for Euler equations*, Numerische Mathematik, 73 (1996), pp. 119–130.

[36] F. RENAC, *Entropy stable DGSEM for nonlinear hyperbolic systems in nonconservative form with application to two-phase flows*, Journal of Computational Physics, 382 (2019).
[37] P. L. ROE, *Approximate Riemann solvers, parameter vectors, and difference schemes*, J. Comput. Phys., 43 (1981), pp. 357-372.

[38] M. E. SCHONBEK, *Second-order conservative schemes and the entropy condition*, Mathematics of Computation, 44 (1985), pp. 31–38.

[39] D. SERRE, *Systems of conservation laws. 1*, Cambridge University Press, Cambridge, 1999. Hyperbolicity, entropies, shock waves, Translated from the 1996 French original by I. N. Sneddon.

[40] J. SHI, Y. ZHANG, AND C. SHU, *Resolution of high order WENO schemes for complicated flow structures*, Journal of Computational Physics, 186 (2003), pp. 690–696.

[41] C. W. SHU, *High-order finite difference and finite volume WENO schemes and discontinuous Galerkin methods for CFD*, International Journal of Computational Fluid Dynamics, 17 (2003), pp. 107–118.

[42] E. TADMOR, *The numerical viscosity of entropy stable schemes for systems of conservation laws. I*, Mathematics of Computation, 49 (1987), pp. 91–103.

[43] E. TADMOR, *Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems*, Acta Numerica, 12 (2003), pp. 451–512.

[44] E. TORO, M. SPRUCE, AND W. SPEARES, *Restoration of the contact surface in the HLL-Riemann solver*, Shock waves, 4 (1994), pp. 25–34.

[45] E. F. TORO, *Riemann solvers and numerical methods for fluid dynamics*, Springer-Verlag, Berlin, third ed., 2009. A practical introduction.

[46] B. VAN LEER, *Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method*, Journal of computational Physics, 32 (1979), pp. 101–136.

[47] H. ZAKERZADEH AND U. FJORDHOLM, *High-order accurate, fully discrete entropy stable schemes for scalar conservation laws*, IMA Journal of Numerical Analysis, 36 (2015), pp. 633–654.