# Well-posedness and numerical approximation of nonlinear conservation laws with hysteresis

Paola Goatin, Stefan Moreti

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

## To cite this version:

Paola Goatin, Stefan Moreti. Well-posedness and numerical approximation of nonlinear conservation laws with hysteresis. 2026. hal-05477008

HAL Id: hal-05477008

https://inria.hal.science/hal-05477008v1

Preprint submitted on 26 Jan 2026

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 CC BY 4.0 - Attribution - International License
# Well-posedness and numerical approximation of nonlinear conservation laws with hysteresis

Paola Goatin *

Stefan Moreti†

## Abstract

This article studies the Cauchy problem for the scalar conservation law

$$ \partial_t u + \partial_t w + \partial_x f(u) = 0, $$

where $w(x, t) = [\mathcal{F}(u)(x, t)]$ is the output of a specific hysteresis operator, namely the Play hysteresis operator, and $f$ is a $C^2$ convex flux function. The hysteresis operator models a rate-independent memory effect, introducing a specific non-local feature into the partial differential equation. We define a suitable notion of entropy weak solution and analyse in detail the Riemann problem. Furthermore, a Godunov-type finite volume numerical scheme is developed to compute approximate solutions. The convergence of the scheme for BV initial data provides the existence of an entropy weak solution. Finally, a stability estimate is established, implying the uniqueness and overall well-posedness of the entropy weak solution.

## 1 Introduction

In this work, we deal with the Cauchy problem for nonlinear conservation laws with hysteresis as follows:

$$ \left\{ \begin{array}{ll} \partial_t u + \partial_t w + \partial_x f(u) = 0 & \text{in } \mathbb{R} \times [0, T[, \\ w = [\mathcal{F}(u, w_0)] & \text{in } \mathbb{R} \times [0, T[, \\ u(x, 0) = u_0(x) & \text{in } \mathbb{R}, \\ w(x, 0) = w_0(x) & \text{in } \mathbb{R}, \end{array} \right. \qquad (1.1) $$

where $\mathcal{F}$ is the so-called Play hysteresis operator, $u_0 \in \text{BV}(\mathbb{R}) \cap \mathbf{L}^1(\mathbb{R})$ is the initial datum for the solution $u$, $w_0 \in \text{BV}(\mathbb{R}) \cap \mathbf{L}^1(\mathbb{R})$ is a suitable space-dependent function for the initial values of the output $w$, $f \in C^2(\mathbb{R}; \mathbb{R})$ is a nonlinear strictly convex function and $T > 0$ is fixed. As for classical conservation laws, the results for $f$ concave can be deduced from the convex case, so we will not describe it in detail.

Hysteresis is a phenomenon commonly observed in various natural and engineered systems, typically characterized by a lag or delay in the system's response to changes in

<footnote>* Université Côte d'Azur, Inria, CNRS, LJAD, 2004 route des Lucioles - BP 93, 06902 Sophia Antipolis Cedex, France. E-mail: paola.goatin@inria.fr</footnote>

<footnote>† Department of Mathematics, University of Trento, Italy. E-mail: stefan.moreti@unitn.it</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: The Play hysteresis operator</figcaption></figure>

the input. For comprehensive accounts of mathematical models for hysteresis and their use in connection with PDEs, we refer the reader to [18] and [27]. Among the various mathematical models used to describe such behaviour, one classical choice is the Play hysteresis operator $\mathcal{F}$

$$ \mathcal{F}: \mathbf{B}([0, T]) \times \mathbb{R} \to \mathbf{B}([0, T]) $$

which represents a memory dependent input-output relationship between a pair of time-dependent scalar functions belonging to a certain functional space $\mathbf{B}([0, T])$. As the unknown $u$ in the PDE is a function both depending on a space and on a time variable, we consider this relationship between the pair $(t \mapsto u(x, t), t \mapsto w(x, t))$ for every $x$, so formally we define

$$ w(x, t) := [\mathcal{F}(u(x, \cdot), w_0(x))](t), \quad \text{a.e. } (x, t). \qquad (1.2) $$

For each fixed $x \in \mathbb{R}$, we briefly describe the input-output relationship of the Play operator referring to Figure 1 (for a more detailed description see [27, Section III.1]). Given an amplitude $a > 0$, we denote

$$ \mathcal{L} := \{(u, w) \in \mathbb{R}^2, |u - w| \le a\} $$

the strip of the feasible states of the system. If the pair $(u(x, t), w(x, t))$ satisfies $|u(x, t) - w(x, t)| < a$ for some $t > 0$, that is if it belongs in the interior of $\mathcal{L}$, and if the input $u(x, \cdot)$ changes in time, then the output $w(x, \cdot)$ will not change until the pair $(u, w)$ will possibly reach one of the two boundary lines of $\mathcal{L}$. If $w(x, t) = u(x, t) - a$, that is the pair $(u, w)$ is on the lower boundary of $\mathcal{L}$ and if the input $u$ increases, then the output $w$ will increase together with $u$; if instead $u$ decreases then $w$ stays constant and the pair $(u, w)$ enters the interior of $\mathcal{L}$. If $(u, w)$ belongs to the upper boundary of $\mathcal{L}$, then the behaviour is symmetric, reversing the role of the monotonicity of $u$.

Given an initial state $(u_0(x), w_0(x)) \in \mathcal{L}$ and the evolution of $u(x, \cdot)$, we can then trace the evolution of $w(x, \cdot)$. We can notice that the value of $w(x, t)$ is not determined pointwisely by $u(x, t)$, indeed for any fixed value of $u$ we have more than one possible value $w$ in the strip $\mathcal{L}$. Hence $w(x, t)$, besides its initial value $w(x, 0) = w_0(x)$, is determined by the whole history $u(x, \cdot)|_{[0,t]}$, i.e. there is a memory effect involved.
In particular, this memory effect is rate-independent, which means that the relation between $u$ and $w$ does not depend on the time-derivative of $u$. Notice that this requirement is essential if we want to draw hysteresis relations as in Figure 1 and it is a general feature of hysteresis phenomena.

As shown in [27] (see also [18]), if we assume $\mathbf{W}^{1,1}$ regularity in time for the input $u$, this heuristic description of the relationship $w = [\mathcal{F}(u, w_0)]$ can be rigorously characterized by the following variational inequality:

$$|u - w| \le a, \quad (u - w - v)\partial_t w \ge 0 \quad \forall v \text{ such that } |v| \le a. \qquad (1.3)$$

In such case also $w(x, \cdot) \in \mathbf{W}^{1,1}([0, T])$ and its time-derivative is described for almost every $t$ as follows

$$\partial_t w(x, t) = \begin{cases} \partial_t u(x, t) & \text{if } w(x, t) = u(x, t) - a \text{ and } \partial_t u(x, t) \ge 0, \\ \partial_t u(x, t) & \text{if } w(x, t) = u(x, t) + a \text{ and } \partial_t u(x, t) \le 0, \\ 0 & \text{otherwise.} \end{cases} \qquad (1.4)$$

In the same work [27], it is also shown how the Play operator can be uniquely extended to an operator that maps continuous inputs to continuous outputs, preserving the same heuristic description.

In this paper, we consider $\mathcal{F}$ as applied to the entropy weak solution of a conservation law, which in general is neither in $\mathbf{W}^{1,1}([0, T])$ or in $\mathbf{C}^0([0, T])$ with respect to $t$, for any fixed $x$. The extension of hysteresis operators to non-regular inputs is a well studied problem, see [8], [19] and [25]. For our intent, we only extend $\mathcal{F}$ to functions with jump discontinuities, performing the following approximation and limit procedure: if the input $u$ has a jump discontinuity at time $t$, we consider a sequence of continuous functions $u_\epsilon$ that fill the discontinuity in a monotone way and such that $u_\epsilon \to u$ pointwisely; then we define $w := \lim_{\epsilon \to 0}[\mathcal{F}(u_\epsilon, w_0)]$ (see Figure 2 for a specific example). In particular, such a construction of the output $w$ as a function in $\mathbf{L}^1([0, T])$ is independent of how we monotonically fill the jump of the input $u$, because of the rate-independence property of the Play operator (see [4] for all the details). We then have the following characterization, which is proven in [4] and will inspire (1.8) in the definition of entropy weak solution.

**Proposition 1.1.** Fix $x \in \mathbb{R}$ and suppose $w_0(x) \in \mathbb{R}$, $u(x, \cdot), w(x, \cdot) \in \text{BV}([0, T; \mathbb{R})$ with a finite number of jump discontinuities, $u(x, 0) = u(x, 0+)$ and $w(x, 0) = w(x, 0+) = w_0(x)$. Then the following are equivalent:

1. $w(x, t) = \mathcal{F}[u(x, \cdot), w_0(x)](t)$ for almost every $t \in [0, T]$;

2. for almost every $t \in [0, T[$, $|w(x, t+) - u(x, t+)| \le a$ and

$$\int_0^t (u(x, s+) - w(x, s+))d(\partial_t w(x, s)) \ge a|\partial_t w|([0, t]), \qquad (1.5)$$

where $\partial_t w$ is interpreted as the measure associated to the distributional time-derivative of $w$, $|\partial_t w|$ its total variation and $u(x, s+), w(x, s+)$ are the right-continuous in time representative of $u$ and $w$ for fixed $x \in \mathbb{R}$.

Given the description of the Play operator $\mathcal{F}$, we introduce the definition of entropy weak solution, which also includes a weaker notion of the relationship $w = [\mathcal{F}(u, w_0)]$.
<figure><img src="image_4.png" /><figcaption>Figure 2: An explicit example of the operator $\mathcal{F}$ applied to $u$ with a jump discontinuity. Here $a = 1$, $u = u_{-} = 0$ for $t < t^{*}$, $u = u_{+} = 2$ for $t > t^{*}$ and $w_{0} = 0$; consequently $w = w_{0}$ for $t < t^{*}$ and $w = w^{*} = 1$ for $t > t^{*}$. In red we highlight the path followed by the couple $(u_{\varepsilon}, w_{\varepsilon})$ which, after the limiting procedure, collapses to $(u_{-}, w_{0})$ for $t < t^{*}$ and $(u_{+}, w^{*})$ for $t > t^{*}$.</figcaption></figure>

**Definition 1.2.** A couple of functions $u, w \in C^{0}([0, T[; L_{\text{loc}}^{1}(\mathbb{R}; \mathbb{R})])$ is an entropy weak solution to (1.1) if:

i) it holds

$$ \int_{0}^{T} \int_{-\infty}^{+\infty} \left[ \left( |u-k| + |w-\hat{k}| \right) \partial_t \phi + \operatorname{sign}(u-k) (f(u)-f(k)) \partial_x \phi \right] dx dt \\ + \int_{-\infty}^{+\infty} (|u_0(x)-k| + |w_0(x)-\hat{k}|) \phi(x,0) dx \geq 0, \quad (1.6) $$

for every $\phi \in C_{c}^{1}(\mathbb{R} \times [0, T[; \mathbb{R}^{+}))$ and for every couple $(k, \hat{k}) \in \mathcal{L}$;

ii) for almost every $(x, t) \in \mathbb{R} \times [0, T[$ it holds

$$ |u(x,t) - w(x,t)| \leq a; \quad (1.7) $$

iii) for almost every $x, u(x, \cdot), w(x, \cdot) \in L^{2}(\mathbb{R}; \mathbb{R})$ and the distributional derivative $\partial_{t}w$ is a measure on $\mathbb{R} \times [0, T[$ that satisfies

$$ \frac{1}{2} \int_{\mathbb{R}} (u(x,t)^2 - u_0(x)^2) dx + \frac{1}{2} \int_{\mathbb{R}} (w(x,t)^2 - w_0(x)^2) dx \leq -a |\partial_t w| (\mathbb{R} \times ]0, t[), \quad (1.8) $$

for almost every $t \in ]0, T[$.

Condition (1.6) is the adaptation to the conservation equation in (1.1) of the classical Kružkov entropy condition [20] and it is inspired by the one introduced in [28]. Also (1.8) was introduced in [28] and it can be interpreted as a weak hysteresis relationship since,
by supposing more regularity on $u$ and $w$, such as space-time Sobolev regularity, and by using the strong form of the PDE: $\partial_t u + \partial_t w + \partial_x f(u) = 0$, (1.8) is equivalent to

$$ \int_0^t \int_{\mathbb{R}} (u(x, s) - w(x, s)) \partial_t w(x, s) \, dx \, ds \geq a \int_0^t \int_{\mathbb{R}} |\partial_t w(x, s)| \, dx \, ds, \quad (1.9) $$

(see [4, 28] for details). Note that (1.9) can be interpreted as the extension to space-time dependent functions of the hysteresis relationship (1.5), when the input and the output are $\mathbf{W}^{1,1}$ in time functions.

Besides [4, 28], scalar hyperbolic conservation laws with hysteresis have been investigated in various applied contexts such as [1, 2, 5, 9–11, 14, 16, 17, 22–24, 29]. In [22, 23] the authors study an initial-boundary value problem where the PDE is similar to ours, with either linear or monotone flux, motivated by application to transport in porous media, where hysteresis is a common feature. There, the hysteresis operator is respectively either a rather general convex-sided Play model or a finite sum of $K$-nonlinear Play operators. They prove well-posedness for their initial-boundary value problem in the class of functions $C^0([0, T]; L^2)0, l$, where $[0, T]$ is the time interval and $]0, l[$, $l > 0$, is the spatial domain, and also existence of time-differentiable solutions in the case with linear flux, differentiable boundary data and 0 initial condition. Their approach relies on the abstract theory of Cauchy problems involving multivalued operators defined on Banach spaces. [28] provides the well-posedness of a Cauchy problem similar to (1.1), but with the presence of a different hysteresis operator, namely the completed delayed Relay. These results can be extended to the Preisach operator [27]. The proof is based on the time-discretization of the partial differential equation, the construction of the corresponding approximate solutions, and, thanks to a priori estimates, their convergence to an exact solution as the discretisation is refined. We also refer the reader to [17], where the previously cited results are well summarized, and to [16] which focuses on entropy conditions for such and similar equations.

Unlike [23, 28] and inspired by the recent works [4, 5], we adopt a more constructive approach focused on the study of characteristics via the resolution of Riemann problems. Further contributions in this direction include [1, 2, 9–11, 14], which study wave propagation in hyperbolic models with hysteresis in the context of traffic flow models (see also [29]). The Cauchy problem (1.1) was studied in the case of linear flux $f(u) = u$ with the Play operator in [4], providing a full analysis of the associated Riemann problems. The same full analysis was also done in [5], where the case of linear flux is still studied but with a more general, versatile and complex hysteretic relationship, namely the one given by the Preisach operator (see [27, Section IV.1] for the description of such hysteresis operator). In particular, the present paper can be seen as an extension of [4] to the case of a general nonlinear convex flux $f$. While in [4, 5] the well-posedness of the Cauchy problem is based on the wave-front tracking approximations, which is a standard tool for proving existence and stability of hyperbolic systems of conservation laws [6], here we develop a finite volume numerical scheme, which is shown to converge to the entropy weak solution. Such scheme can be directly applied to the case of linear flux treated in [4], see Remark 3.13.

Actually, none of the above references uses a Godunov-type numerical approximation and its limit procedure to prove the existence of solutions. This indeed seems to be a novelty of our analysis, in addition to the treatment of the passage to the limit in the
hysteresis relationship (see Section 3.4), similarly to what is done in [4,5]. To the best of the authors' knowledge, the most closely related work remains [22], where also an explicit upwind finite volume numerical scheme is analysed.

We refer the reader to the books [6, 12, 15] for general theory on scalar conservation laws and to [13, 21, 26] for more specific results about finite volume numerical schemes.

This article is then structured as follows: in Section 2 we deal with the Riemann problem, in which the solutions are combinations of shock and rarefaction waves as the flux is nonlinear. In Section 3 the Godunov-type numerical scheme is developed; thanks to BV estimates, it is shown that the approximate solutions generated by the scheme converge as the mesh size tends to 0 to the entropy weak solution of the Cauchy problem (1.1) establishing existence; some numerical examples and simulations are highlighted at the end of this section; finally, in Section 4 it is shown that the solutions to the Riemann problem constructed in Section 2 are entropy admissible, and a stability theorem is stated, which implies uniqueness of the entropy solution.

## 2 The Riemann problem

The Riemann problem associated to (1.1) is the Cauchy problem with the following initial data

$$u_0(x) = \begin{cases} u_l & x < 0, \\ u_r & x \ge 0, \end{cases} \quad \text{and} \quad w_0(x) = \begin{cases} w_l & x < 0, \\ w_r & x \ge 0, \end{cases} \qquad (2.1)$$

such that $|u_0(x) - w_0(x)| \le a$.

By solutions of the Riemann problem, we intend a couple of functions $(u, w) \in C^0([0, T[; L^1_{\text{loc}}(\mathbb{R}; \mathbb{R}^2))$ satisfying the entropy weak formulation of the PDE (1.6) and such that $w(x, t) = [\mathcal{F}(u(x, \cdot), w_0(x))(t)$ for almost every $x \in \mathbb{R}$ and $t \in [0, T[$. It will turn out that, $u$ being the solution of a Riemann problem, $u(x, \cdot)$ is a piecewise continuous function with only jump discontinuities for almost every $x \in \mathbb{R}$. Hence $\mathcal{F}$ has to be interpreted as applied to functions with jump discontinuities, see Proposition 1.1 and the paragraph above.

First of all, we notice that, if an entropy solution $(u, w)$ is discontinuous along a curve $(\sigma(t), t)$, with $u_{\pm} = u_{\pm}(t) := u(\sigma(t)\pm, t)$ and $w_{\pm} = w_{\pm}(t) := w(\sigma(t)\pm, t)$, then the Rankine-Hugoniot condition

$$f(u_{-}) - f(u_{+}) = \sigma'(t)(u_{-} - u_{+} + w_{-} - w_{+}). \qquad (2.2)$$

follows from (1.6). Indeed, it is enough to choose $k \ge \max(u_{-}, u_{+})$ and $\hat{k} \ge \max(w_{-}, w_{+})$ in (1.6) and integrate it by parts to infer

$$\sigma'(t)[u_{-} + w_{-} - u_{+} - w_{+}] - [f(u_{-}) - f(u_{+})] \ge 0; \qquad (2.3)$$

similarly, with $k \le \min(u_{-}, u_{+})$ and $\hat{k} \le \min(w_{-}, w_{+})$ we deduce

$$-\sigma'(t)[u_{-} + w_{-} - u_{+} - w_{+}] + [f(u_{-}) - f(u_{+})] \ge 0, \qquad (2.4)$$

which, together with the previous inequality, implies (2.2). The inequalities (2.3) and (2.4) are derived as in the case of no hysteresis (see e.g. [6, Theorem 4.3 and proof of Theorem 4.4]).

For the explicit construction of solutions, we divide the analysis in the following cases: $u_l = u_r$, $u_l < u_r$ and $u_l > u_r$.
## 2.1 $u_l = u_r$:

In this case the weak solution is $(u(x, t), w(x, t)) = (u_0(x), w_0(x))$ for each $(x, t)$. Notice that if $w_l \neq w_r$ then this solution has a stationary discontinuity for $w$, which satisfies (2.2). Moreover, $w = [\mathcal{F}(u, w_0)]$ holds trivially.

## 2.2 $u_l < u_r$: rarefaction waves

If there was no hysteresis term in the equation, since the flux is convex, we would expect a rarefaction wave type solution for the unknown $u$. Based on the possible different directions of propagation of the waves, we consider the following subcases.

### 2.2.1 $f'(u_l) \geq 0$:

By convexity, in this case $f'(u) \geq 0$ for $u \in [u_l, u_r]$ and the rarefaction has positive speed. Therefore, $u(x, t) \equiv u_l$ and $w(x, t) \equiv w_l$ for $x < 0$. For $x > 0$ fixed, we expect $t \mapsto u(x, t)$ to decrease in a monotone way from $u_r$ to $u_l$, reaching the value $u_l$ in finite time. If we impose the relationship $w = [\mathcal{F}(u, w_0)]$, then it should hold that for each $t$ such that $w_r - a \leq u(x, t) \leq u_r$, as $u(x, \cdot)$ is monotone decreasing, at least formally $w_t(x, t) = 0$, hence $w(x, t) = w_r$. Whereas, if $u(x, t) \leq w_r - a$ and $u(x, \cdot)$ keeps decreasing, the couple $(u, w)$ should follow the upper boundary $u = w + a$ of the hysteresis region, so $w_t(x, t) = u_t(x, t)$.

By this analysis we deduce that, at least formally, $u$ should satisfy the conservation law

$$ \partial_t u + \partial_x \bar{f}_{w_r}(u) = 0, \qquad (2.5) $$

where

$$ \bar{f}_{w_r}(u) := \begin{cases} \frac{1}{2}f(u) + \frac{1}{2}f(w_r - a) & u \leq w_r - a, \\ f(u) & w_r - a \leq u \leq w_r, \end{cases} \qquad (2.6) $$

is still a convex, piecewise $C^1$, flux, see Figure 3, left. Therefore, the classical Riemann problem for (2.5) is solved in a standard way.

In particular, if $w_r - a$ is such that $u_l < w_r - a < u_r$, then $u$ consists of two rarefaction waves, the first from $u_r$ to $w_r - a$, with $\partial_t w = 0$, and the second from $w_r - a$ to $u_l$, with $\partial_t u = \partial_t w$, separated by the intermediate state $u \equiv w_l - a$. Then also $w$ can be computed as $w = [\mathcal{F}(u, w_0)]$ and it will consist of a rarefaction wave that overlaps with the second rarefaction wave of $u$, connecting $w_r$ to the state $u_l + a$, see Figure 4.

If instead we suppose that $w_r - a \leq u_l < u_r$ then $\bar{f}_w(u) = f(u)$ for $u \in [u_l, u_r]$, so $u$ consists of only one rarefaction and $w$ remains constant as $\partial_t w = 0$.

The opposite extreme is when $u_l < u_r = w_r - a$, which corresponds to the case when the couple $(u, w)$ at time $t = 0$ for $x > 0$ belongs already to the upper boundary of the hysteresis region. Then $\bar{f}_w(u) = \frac{1}{2}f(u) + \frac{1}{2}f(w_r - a)$ for $u \in [u_l, u_r]$, so we have a rarefaction wave for both $u$ and $w$ with $w(x, t) = u(x, t) + a$ for each $x > 0$ and $t > 0$.

It is easy to see that by construction $(u, w)$ is a weak solution of the PDE (1.6), as whenever $\partial_t w = 0$, $u$ solves $\partial_t u + \partial_x f(u) = 0$, and when $\partial_t w = \partial_t u$, then $\partial_t u + \frac{1}{2} \partial_x f(u) = 0$ and $\partial_t w + \frac{1}{2} \partial_x f(u) = 0$. Moreover, the relationship $w = [\mathcal{F}(u, w_0)]$ holds in the strong classical sense as $u(x, \cdot) \in \mathbf{W}^{1,1}(0, T)$ for almost every $x$ and we may develop stationary shocks for $w$.
<figure><img src="image_5.png" /><figcaption>Figure 3: An example of $\bar{f}_{w_r}$ (left) and $\hat{f}_{w_l}$ (right) compared with $f$ (gray); the additive constants $\frac{1}{2}f(w_r-a)$ in (2.6) and $\frac{1}{2}f(w_r+a)$ in (2.7) are needed for the continuity of $\bar{f}_{w_r}$ and $\hat{f}_{w_l}$ respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: On the left, the couple $(u, w)$ is shown in the hysteresis plane: the arrow represents the path followed by $t \mapsto (u(x,t), w(x,t))$ for $x > 0$; the dashed line instead represents the stationary shock. On the right, the solution $(u, w)$ in the $(x, t)$ plane: it consists of two rarefaction waves and a stationary shock for $w$ along $x = 0$. In this example, the flux function is $f(u) = \frac{1}{2}u^2$, the parameter $a = 1$, and the initial data $(u_l, w_l) = (1, 0.5)$ and $(u_r, w_r) = (3, 3)$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: An example of the case when $f'(u_r) \le 0$. In this example, the flux function is $f(u) = \frac{1}{2}u^2$, the parameter $a = 1$, and the initial data $(u_l, w_l) = (-3, -3)$ and $(u_r, w_r) = (-1.5, -1)$.</figcaption></figure>

### 2.2.2 $f'(u_r) \le 0$:

In this case we still deal with a rarefaction wave, but moving to the left. As a result, for $x > 0$, $u(x, t) = u_r$ and $w(x, t) = w_r$ for each $t$. Instead, for $x < 0$, $t \mapsto u(x, t)$ increases in a monotone way from $u_l$ to $u_r$. Reasoning as in Subcase 2.2.1, we can check that $u$ formally solves

$$ \partial_t u + \partial_x \hat{f}_{w_l}(u) = 0, $$

with

$$ \hat{f}_{w_l}(u) = \begin{cases} f(u) & w_l - a \le u \le w_l + a, \\ \frac{1}{2}f(u) + \frac{1}{2}f(w_l + a) & w_l + a \le u, \end{cases} \qquad (2.7) $$

see Figure 3, right. Hence, if $u_l < w_l + a < u_r$ the solution $(u, w)$ for $x < 0$ consists of two rarefaction waves, one only for $u$ and the other for both $u$ and $w$, see Figure 5; if instead $u_r \in [u_l, w_l + a]$, then $w(x, t) = w_0(x)$ for every $t$ and there is only one rarefaction wave for $u$; instead if $u_l = w_l + a$, then we have only one rarefaction wave for both $u$ and $w$ and $w(x, t) = u(x, t) - a$ for $x > 0$ and for every $t$.

### 2.2.3 $f'(u_l) \le 0 \le f'(u_r)$:

This subcase is the combination of the previous two Subcases 2.2.1 and 2.2.2. In particular, looking at the direction of the waves, for $x > 0$ we expect $u(x, \cdot)$ to decrease from $u_r$ to $u_*$ where $u_* \in [u_l, u_r]$ is the point of minimum of $f$ on $[u_l, u_r]$. Instead, for $x < 0$, $t \mapsto u(x, t)$ increases from $u_l$ to $u_*$. The idea is then to consider the two Riemann problems with data

$$ u_0^1(x) = \begin{cases} u_* & x < 0, \\ u_r & x \ge 0, \end{cases} \quad \text{and} \quad u_0^2(x) = \begin{cases} u_l & x < 0, \\ u_* & x \ge 0, \end{cases} $$

solve the first one as in Subcase 2.2.1 and consider its solution $(u_1, w_1)$ restricted to $x > 0$, solve the second one as in Subcase 2.2.2 and consider the solution $(u_2, w_2)$ on $x < 0$. We
<figure><img src="image_8.png" /><figcaption>Figure 6: This is an example of solution when $f'(u_l) \le 0 \le f'(u_r)$; in this example, the flux function is $f(u) = \frac{1}{2}u^2$, the parameters $a = 1$, and the initial data $(u_l, w_l) = (-2, -1.5)$ and $(u_r, w_r) = (1, 1.5)$.</figcaption></figure>

finally obtain the solution $(u, w)$ to the original Riemann problem by gluing together these two restrictions.

Notice that, since $f'(u_*) = 0$, then $u(0-, t) = u_* = u(0+, t)$ so $u(x, t)$ is continuous at $x = 0$ for each $t > 0$. Our solution may have a discontinuity in $w$ at $x = 0$, however, as we already pointed out, stationary discontinuities for $w$ are admissible as long as $u$ is continuous, see again (2.2) and Figure 6.

## 2.3 $u_l > u_r$: shock waves

Now we expect $u$ to develop shock discontinuities. Again, in order to study $u(x, \cdot)$, and hence understand $w(x, \cdot)$, we need to distinguish the cases when the shock waves propagate with null, positive or negative speed.

### 2.3.1 $f(u_l) = f(u_r)$:

We expect the formation of a shock for $u$ with 0 speed, so, for each $x$, $u(x, \cdot) = u_0(x)$, and, by imposing $w = [\mathcal{F}(u, w_0)]$, $w(x, \cdot) = w_0(x)$. Then the couple $(u(x, t), w(x, t))(x, t) = (u_0(x), w_0(x))$, satisfying the Rankine-Hugoniot condition (2.2), is also an entropy weak solution to our problem.

**Remark 2.1.** Even if the case $f(u_l) = f(u_r)$ might seem trivial, we wanted to highlight it. Indeed, it may happen that $(u_l - u_r + w_l - w_r) = 0$, with $u_l \ne u_r$, so the Rankine-Hugoniot condition (2.2) does not give any information about the speed of the shock. Of course, we must choose $\sigma'(t) = 0$ to have $w = [\mathcal{F}(u, w_0)]$, but, even without imposing $w = [\mathcal{F}(u, w_0)]$, we will show that the only entropy solution is the one with $\sigma'(t) = 0$, see Proposition 4.1.
### 2.3.2 $f(u_l) > f(u_r)$:

In this subcase, we expect the discontinuity to propagate with positive speed hence for $x < 0$, $u(x, \cdot) \equiv u_l$ so also $w(x, \cdot) \equiv w_l$. Instead, for each $x > 0$, we expect $t \mapsto u(x, t)$ to increase from $u_r$ to $u_l$. Similarly to Subcase 2.2.1 we consider the following flux

$$ \hat{f}_{w_r}(u) = \begin{cases} f(u) & u_r \le u \le w_r + a, \\ \frac{1}{2}f(u) + \frac{1}{2}f(w_r + a) & w_r + a \le u, \end{cases} $$

and again, at least formally, $u$ should satisfy

$$ \partial_t u + \partial_x \hat{f}_{w_r}(u) = 0. $$

When $u_r \le u_l \le w_r + a$, then $\hat{f}_{w_r} = f(u)$ in $[u_r, u_l]$ and the solution for $u$ consists of only one shock and again $w_t = 0$. The other extreme case, that is when $u_r = w_r + a$, generates a solution with a single shock for both $u, w$, connecting $(u_r, u_r - a)$ and $(u_l, u_l - a)$. In all these cases, the speed of the shock is given by the Rankine-Hugoniot condition (2.2).

Instead, when $u_r < w_r + a < u_l$, in order to solve the Riemann problem we analyse the upper concave envelope of $\hat{f}_{w_r}$ between $u_r$ and $u_l$, since $u_l > u_r$ and $\hat{f}_{w_r}$ is certainly not concave. We then denote the following two quantities

$$ \mu_r := \frac{f(w_r + a) - f(u_r)}{(w_r + a) - u_r}, \quad \text{and} \quad \mu_l := \frac{1}{2} \frac{f(u_l) - f(w_r + a)}{u_l - (w_r + a)}, $$

used to construct the upper concave envelope. Notice first of all that $\mu_l > 0$, as $f$ is convex, $u_r \le w_r + a \le u_l$ and $f(u_l) > f(u_r)$.

**$\mu_r \ge \mu_l$**: Under this assumption, the upper concave envelope between $[u_r, u_l]$ of the graph of $\hat{f}_{w_r}$ consists of the two segments connecting the points $(u_r, \hat{f}_{w_r}(u_r))$, $(w_r + a, \hat{f}_{w_r}(w_r + a))$ and $(w_r + a, \hat{f}_{w_r}(w_r + a))$, $(u_l, \hat{f}_{w_r}(u_l))$, the first segment with slope $\mu_r$ and the second with slope $\mu_l$, see Figure 7, left. As a result, when solving the Riemann problem we end up with two shock waves for $u$, one with velocity $\mu_r$ and the other $\mu_l$. By imposing the relationship $w = [\mathcal{F}(u, w_0)]$, we then observe that $w$ has only one discontinuity with speed $\mu_l$ connecting the state $w_r$ to $u_l - a$, see Figure 8. Notice that the Rankine-Hugoniot condition (2.2) holds for both shocks, consequently (1.6) also holds. In fact, for the fastest one, only $u$ is discontinuous, and Rankine-Hugoniot gives $\mu_r$. For the second, instead, we have a discontinuity with the same amplitude for both $u$ and $w$, and Rankine-Hugoniot gives $\mu_l$.

**$\mu_r < \mu_l$**: In this case, the upper concave envelope of $\hat{f}_{w_r}$ consist in just one segment with slope

$$ \frac{\frac{1}{2}f(u_l) + \frac{1}{2}f(w_r + a) - f(u_r)}{u_l - u_r} = \frac{\mu_l I_l + \mu_r I_r}{I_r + I_l} = \tilde{\mu}, $$

where we denoted $I_r = (w_r + a) - u_r$ and $I_l = u_l - (w_r + a)$, see Figure 7, right. So the solution $u$ should consists of one jump discontinuity connecting $u_r$ and $u_l$ with speed given by the slope of that segment. By imposing the relationship $w = \mathcal{F}(u)$, it should hold that $w$ has a shock discontinuity connecting the states $w_r$ and $u_l - a$ with the same speed. However, the couple $(u, w)$ constructed in such a way is not a weak solution to the PDE. Indeed the slope $\tilde{\mu}$ does not satisfy the Rankine-Hugoniot condition which gives

$$ \mu = \frac{f(u_l) - f(u_r)}{u_l - u_r + (u_l - a) - w_r} = \frac{2I_l\mu_l + I_r\mu_r}{I_r + 2I_l}. \qquad (2.8) $$
<figure><img src="image_9.png" /><figcaption>Figure 7: Examples of $\hat{f}_{w_r}$. Left: if $\mu_r \ge \mu_l$, the upper concave envelope between $[u_r, u_l]$ of the graph of $\hat{f}_{w_r}$ consists of the two segments drawn in the picture. Right: if $\mu_r < \mu_l$, the upper concave envelope consists only of one segment.</figcaption></figure>

Necessarily, since we only have one discontinuity both in $u$ and $w$, we must impose its velocity to be equal to $\mu$, so that $(u, w)$ is indeed a weak solution, see Figure 9. Notice that, by doing so, we still ensure that $w = [\mathcal{F}(u, w_0)]$. We will call this kind of shock a “fast shock”.

**Remark 2.2.** In the case $\mu_r > \mu_l$ one could also considered the couple consisting of one “fast shock” with slope given by $\mu$. Such couple $(u, w)$ still satisfies the Rankine-Hugoniot condition with also $w = [\mathcal{F}(u, w_0)]$ but, as we will see in Section 4, it does not satisfy the entropy weak formulation (1.6).

**Remark 2.3.** In the case when $\mu_r < \mu_l$, the heuristic idea is that, since the shock wave is too fast, the couple $(u, w)$ has no time to reach the intermediate state $(u_r, w_r + a)$ but it directly jumps to the state $(u_l, u_l - a)$. Moreover, the slope $\mu$, as defined in (2.8), can be interpreted as a weighted average of the velocities of the two shocks, $\mu_l$ and $\mu_r$, with respective weights $I_r$ and $2I_l$. The weight $I_r$ corresponds to the length of the interval in $u$ where $w$ remains constant, while $I_l$ represents the length of the interval in $u$ where both $u$ and $w$ change values, which is equivalent to the length of the corresponding interval in $w$. This justifies the factor 2 in front of $I_l$.

**Remark 2.4.** It can be easily seen that if $\mu_r = \mu_l$ then also $\mu = \mu_r = \mu_l$. Hence the two approaches coincide.

### 2.3.3 $f(u_l) < f(u_r)$:

In this subcase instead we expect the discontinuity waves to propagate to the left. We then consider the flux

$$ \bar{f}_{w_l}(u) = \begin{cases} \frac{1}{2}f(u) + \frac{1}{2}f(w_l - a) & u < w_l - a, \\ f(u) & w_l - a \le u \le w_l, \end{cases} $$
<figure><img src="image_10.png" /><figcaption>Figure 8: Example of solution for $u_l > u_r$ and $\mu_r > \mu_l$ with $f(u) = \frac{1}{2}u^2$, $a = 1$, and initial data $(u_l, w_l) = (1.5, 2)$ and $(u_r, w_r) = (0.5, 0)$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: Example of solution for $u_l > u_r$ and formation of a “fast shock” with $f(u) = \frac{1}{2}u^2$, $a = 1$ and initial data $(u_l, w_l) = (3.5, 2)$ and $(u_r, w_r) = (-1, -1)$.</figcaption></figure>
and define

$$ \nu_l := \frac{f(u_l) - f(w_l - a)}{u_l - (w_l - a)}, \quad \text{and} \quad \nu_r := \frac{1}{2} \frac{f(w_l - a) - f(u_r)}{(w_l - a) - u_r}. $$

Reasoning as in Subcase 2.3.2, if $u_l < w_l + a < u_r$ and $\nu_l < \nu_r$ then we have two jump discontinuities, one only for $u$ with velocity $\nu_l$ and one for both $u$ and $w$ with velocity $\nu_r$. If instead $u_l < w_l + a < u_r$ and $\nu_l > \nu_r$, then the weak solution develops a “fast shock” for both $u$ and $w$ with slope

$$ \nu = \frac{f(u_l) - f(u_r)}{u_l - u_r + w_l - (u_r + a)} = \frac{J_l \nu_l + 2 J_r \nu_r}{J_l + 2 J_r}, $$

given by the Rankine-Hugoniot condition, where $J_l := u_l - (w_l - a)$ and $J_r = (w_l - a) - u_r$. Notice that this time necessary $\nu_r < 0$, and the same previous remarks, referring to this subcase, hold. We can also reason as before to deal with the cases when $u_r \le w_l + a$ or $u_l = w_l + a$.

**Remark 2.5.** In this section, when constructing the solution, we dealt with all the possible cases to highlight the influence of the hysteresis term in the equation. Nevertheless, defining the flux

$$ \tilde{f}_w(u) = \begin{cases} \frac{1}{2}f(u) + \frac{1}{2}f(w-a) & u < w-a, \\ f(u) & w-a \le u \le w+a, \\ \frac{1}{2}f(u) + \frac{1}{2}f(w+a) & u \ge w+a, \end{cases} \quad (2.9) $$

we can characterize the solution for $u$ in a more general way. Indeed, except for “fast shocks”, the restriction of $u$ on $x > 0$ satisfies the PDE

$$ \partial_t u + \partial_x \tilde{f}_{w_r}(u) = 0, $$

and its restriction on $x < 0$ satisfies

$$ \partial_t u + \partial_x \tilde{f}_{w_l}(u) = 0. $$

This general property will be exploited in the next section to define the numerical scheme by only distinguishing whether a “fast shock” is formed or not, and not by going through all the above cases.

# 3 Existence of solutions

## 3.1 A Godunov-type finite volume scheme

Let us fix $T > 0$ and consider a mesh $\mathcal{T}$ of $\mathbb{R}$ given by an increasing sequence of points $(x_{i+1/2})_{i \in \mathbb{Z}}$, such that $\mathbb{R} = \cup_{i \in \mathbb{Z}}[x_{i-1/2}, x_{i+1/2}] := \cup_{i \in \mathbb{Z}} K_i$. For simplicity we suppose that $|x_{i-1/2} - x_{i+1/2}| = \Delta x$, for some $\Delta x > 0$, for all $i \in \mathbb{Z}$, and we denote by $x_i$ the middle point of $K_i$. We introduce the time step $\Delta t > 0$ and we denote by $t^n := n\Delta t$ and $K_i^n := K_i \times [t^n, t^{n+1}]$. Given an initial datum $u_0, w_0$ we set $U_m := \inf_{\mathbb{R}} u_0$ and $U_M := \sup_{\mathbb{R}} u_0$, so
that $u_0(x) \in [U_m, U_M]$ for each $x$, and we also suppose the following Courant-Friedrichs-Lewy (CFL) condition, see e.g. [26, Section 6.2], to hold

$$ \Delta t \leq \frac{\Delta x}{2L}, \qquad (3.1) $$

where $L := \max_{u \in [U_m, U_M]} |f'(u)|$ is the Lipschitz constant of $f$ on the interval $[U_m, U_M]$. Given piecewise constant data at time $t^n$,

$$ u(x, t^n) = \sum_{i \in \mathbb{Z}} u_i^n 1_{[x_{i-1/2}, x_{i+1/2}]} \quad \text{and} \quad w(x, t^n) = \sum_{i \in \mathbb{Z}} w_i^n 1_{[x_{i-1/2}, x_{i+1/2}]}, $$

we also denote by $\tilde{u}^n(x, t)$ and $\tilde{w}^n(x, t)$ the exact solution to (1.1) defined on $\mathbb{R} \times [t^n, t^{n+1}]$ via the resolution of the Riemann problems centred at $x_{i+1/2}$, $i \in \mathbb{Z}$. Notice that $\tilde{u}^n(x, t)$ and $\tilde{w}^n(x, t)$ are well defined on the whole time interval as the CFL condition (3.1) ensures not only that no waves cross the cell interfaces $\{x_{i+1/2}\} \times [t^n, t^{n+1}]$, but also that there are no wave interactions in $K_i^n$ for any $n$ and $i$.

A finite volume scheme for the approximation of problem (1.1), with mesh $\mathcal{T}$ and time step $\Delta t$, writes as follows

$$ \left\{ \begin{array}{ll} u_i^{n+1} + w_i^{n+1} = u_i^n + w_i^n - \frac{\Delta t}{\Delta x} (f_{i+1/2}^n - f_{i-1/2}^n), & \forall i \in \mathbb{Z}, \\ u_i^0 + w_i^0 = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} (u_0(x) + w_0(x)) \, dx, & \forall i \in \mathbb{Z}, \end{array} \right. \qquad (3.2) $$

where $u_i^n, w_i^n$ are the averages of $u$ and $w$ respectively on $K_i \times \{t^n\}$ and $f_{i+1/2}^n$ denotes the numerical flux at the interface $x_{i+1/2}$, given by the Godunov flux $f_{i+1/2}^n = f(\tilde{u}^n(x_{i+1/2}, t))$, see e.g. [26, Section 6.2]. Notice that $\tilde{u}^n(x_{i+1/2}, t)$ is constant in time as it is the solution to a Riemann problem centred at $x_{i+1/2}$, and that if $\tilde{u}^n(x_{i-1/2}^-, t) \neq \tilde{u}^n(x_{i-1/2}^+, t)$ then $f(\tilde{u}^n(x_{i-1/2}^-, t)) = f(\tilde{u}^n(x_{i-1/2}^+, t))$ by the Rankine-Hugoniot condition, so $f_{i+1/2}^n$ is well defined for every $i \in \mathbb{Z}$. Such numerical flux depends only on the Riemann data of $u_i^n, u_{i+1}^n$ (see Section 2) and, as noticed in [21, Section 13.5], it can be written as $f_{i+1/2}^n = g(u_i^n, u_{i+1}^n)$ where

$$ g(u_l, u_r) = \begin{cases} \min_{u \in [u_l, u_r]} f(u) & u_l \le u_r, \\ \max_{u \in [u_r, u_l]} f(u) & u_l \ge u_r. \end{cases} \qquad (3.3) $$

The recursive scheme (3.2) allows us to compute the approximate sum of the two variables $u + w$, but not the single values of the unknowns. To complete the algorithm, we need to compute $u_i^n, u_{i+1}^n$. To this end, we analyse what happens to the exact solutions $\tilde{u}^n$ in a cell $K_i^n$. This solution is given by the union of (the restrictions to the corresponding quarters of plane of) the solution to the Riemann problem centered at $(x_{i-1/2}, t^n)$ with data $(u_{i-1}^n, w_{i-1}^n)$ and $(u_i^n, w_i^n)$, and the solution to the Riemann problem centred at $(x_{i+1/2}, t^n)$ with data $(u_i^n, w_i^n)$ and $(u_{i+1}^n, w_{i+1}^n)$.

Suppose for now that in $K_i^n$ there are no “fast shocks”. Then, by Remark 2.5, $\tilde{u}^n$ satisfies weakly the following conservation law

$$ \partial_t u + \partial_x \tilde{f}_{w_i^n}(u) = 0, \qquad (3.4) $$

in the whole cell $K_i^n$, where $\tilde{f}_{w_i^n}$ is defined by (2.9). Hence, using the integral formulation of (3.4), we get the scheme for $u$

$$ u_i^{n+1} = u_i^n - \frac{\Delta t}{\Delta x} (\tilde{f}_{w_i^n}(\tilde{u}^n(x_{i+1/2}^-, t)) - \tilde{f}_{w_i^n}(\tilde{u}^n(x_{i-1/2}^+, t))). \qquad (3.5) $$
<figure><img src="image_12.png" /><figcaption>Figure 10: The cell $K_i^n$ split into $K_i^{n,l}$ and $K_i^{n,r}$</figcaption></figure>

If instead there is a “fast shock” in $K_i^n$, then the identification of the flux is not straightforward in general. We then subdivide $K_i^n$ as

$$K_i^n = [x_{i-1/2}, x_i] \times [t^n, t^{n+1}] \cup [x_i, x_{i+1/2}] \times [t^n, t^{n+1}] := K_i^{n,l} \cup K_i^{n,r},$$

and we consider separately the averages of $\tilde{u}^n(x, t^{n+1})$ on $[x_{i-1/2}, x_i]$ and $[x_i, x_{i+1/2}]$, which we denote respectively by $u_i^{n+1,l}$ and $u_i^{n+1,r}$. Indeed, thanks to the CFL condition (3.1), $u_i^{n+1,l}$ and $u_i^{n+1,r}$ depend respectively only on the solution to the Riemann problems centered at $(x_{i-1/2}, t^n)$ and $(x_{i+1/2}, t^n)$.

If there is no “fast shock” in $K_i^{n,l}$, then (3.4) still holds in $K_i^{n,l}$. Hence, applying the discrete formulation of (3.4) in $K_i^{n,l}$,

$$ \begin{aligned} u_i^{n+1,l} &= u_i^n - 2 \frac{\Delta t}{\Delta x} \left( \tilde{f}_{w_i^n}(\tilde{u}^n(x_i^n, t)) - \tilde{f}_{w_i^n}(\tilde{u}^n(x_{i-1/2}^n, t)) \right) \\ &= u_i^n - 2 \frac{\Delta t}{\Delta x} \left( f(u_i^n) - \tilde{f}_{w_i^n}(\tilde{u}^n(x_{i-1/2}^n, t)) \right). \end{aligned} $$

Notice that the space interval here has length $\Delta x/2$ and that $\tilde{u}^n(x_i^n, t) = u_i^n$ as ensured by the CFL condition (3.1).

Otherwise, if a “fast shock” is generated, then we compute explicitly

$$u_i^{n+1,l} = u_i^n + 2 \frac{\Delta t}{\Delta x} \mu (u_{i-1}^n - u_i^n) \quad (3.6)$$

as the jump between $u_i^n$ and $u_{i-1}^n$ travels with speed $\mu$ given by (recall Subcase 2.3.2)

$$\mu = \mu(u_{i-1}^n, u_i^n, w_i^n) = \frac{f(u_{i-1}^n) - f(u_i^n)}{u_{i-1}^n - u_i^n + (u_{i-1}^n - a) - w_i^n}.$$

Since $\tilde{u}^n(x_i, t) \equiv u_i^n$, then at least formally (3.6) means that the total flux crossing the interface $x_{i-1/2} \times [t^n, t^{n+1}]$ is $\mu(u_{i-1}^n - u_i^n) + f(u_i^n)$.

Similarly, if there is no “fast shock” in $K_i^{n,r}$, then

$$u_i^{n+1,r} = u_i^n - 2 \frac{\Delta t}{\Delta x} \left( \tilde{f}_{w_i^n}(\tilde{u}^n(x_{i+1/2}^n, t)) - f(u_i^n) \right),$$
while, if a “fast shock” is present, then

$$u_i^{n+1,r} = u_i^n - 2 \frac{\Delta x}{\Delta t} \nu(u_i^n - u_{i-1}^n),$$

where (recall Subcase 2.3.3)

$$\nu = \nu(u_i^n, u_{i+1}^n, w_i^n) = \frac{f(u_i^n) - f(u_{i+1}^n)}{u_i^n - u_{i+1}^n + w_i^n - (u_{i+1}^n + a)}.$$

We then finally deduce that in general

$$u_i^{n+1} = u_i^n - \frac{\Delta t}{\Delta x} \left( h_1^-(u_i^n, u_{i+1}^n, w_i^n) - h_1^+(u_{i-1}^n, u_i^n, w_i^n) \right), \quad (3.7)$$

where

$$h_1^+(u_{i-1}^n, u_i^n, w_i^n) = \begin{cases} \tilde{f}_{w_i^n}(\tilde{u}^n(x_{i-1/2}+, t)) & \text{if no "fast shock" in } K_i^{n,l}, \\ \mu(u_{i-1}^n - u_i^n) + f(u_i^n) & \text{if "fast shock" in } K_i^{n,l}, \end{cases} \quad (3.8)$$

and

$$h_1^-(u_{i-1}^n, u_i^n, w_i^n) = \begin{cases} \tilde{f}_{w_i^n}(\tilde{u}^n(x_{i+1/2}-, t)) & \text{if no "fast shock" in } K_i^{n,r}, \\ \nu(u_{i+1}^n - u_i^n) + f(u_i^n) & \text{if "fast shock" in } K_i^{n,r}. \end{cases} \quad (3.9)$$

Once $u_i^{n+1}$ is computed by (3.7), we use (3.2) to also compute $w_i^{n+1}$.

Finally, the complete scheme reads as

$$\left\{
\begin{array}{ll}
u_i^{n+1} = u_i^n - \frac{\Delta t}{\Delta x} \left( h_1^-(u_i^n, u_{i+1}^n, w_i^n) - h_1^+(u_{i-1}^n, u_i^n, w_i^n) \right), & \forall i \in \mathbb{Z}, \\
w_i^{n+1} = u_i^n + w_i^n - \frac{\Delta t}{\Delta x} \left( f_{i+1/2}^n - f_{i-1/2}^n \right) - u_i^{n+1}, & \forall i \in \mathbb{Z}, \\
u_i^0 = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} u_0(x) \, dx, \quad w_i^0 = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} w_0(x) \, dx, & \forall i \in \mathbb{Z}, \\
|u_0(x) - w_0(x)| \le a, & \forall x \in \mathbb{R}.
\end{array}
\right.
\quad (3.10)$$

We point out that the scheme is well-posed. Indeed, in order to be able to solve the Riemann problems to define $\tilde{u}^n(x, t)$ and $\tilde{w}^n(x, t)$, due to the hysteresis relationship, we always supposed that $|u_i^n - w_i^n| \le a$. However, starting with compatible initial data, that is $|u_0(x) - w_0(x)| \le a$, then

$$|u_i^0 - w_i^0| \le \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} |u_0(x) - w_0(x)| \, dx \le a, \quad \forall i \in \mathbb{Z}.$$

Consequently, since the solution to the Riemann problem preserves the inequality $|u-w| \le a$ for every $x$ and $t > 0$, (see Section 2), it holds

$$|u_i^{n+1} - w_i^{n+1}| \le \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} |\tilde{u}^n(x, t^{n+1}) - \tilde{w}^n(x, t^{n+1})| \, dx \le a, \quad \forall i \in \mathbb{Z}, \forall n \in \mathbb{N}. \quad (3.11)$$
To study the properties of the scheme, it will be convenient to rewrite (3.10) as

$$ \left\{
\begin{array}{ll}
u_i^{n+1} = u_i^n - \frac{\Delta t}{\Delta x}(h_1^-(u_i^n, u_{i+1}^n, w_i^n) - h_1^+(u_{i-1}^n, u_i^n, w_i^n)), & \forall i \in \mathbb{Z}, \\
w_i^{n+1} = w_i^n - \frac{\Delta t}{\Delta x}(h_2^-(u_i^n, u_{i+1}^n, w_i^n) - h_2^+(u_{i-1}^n, u_i^n, w_i^n)), & \forall i \in \mathbb{Z}, \\
u_i^0 = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} u_0(x) dx, \quad w_i^0 = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} w_0(x) dx, & \forall i \in \mathbb{Z}, \\
|u_0(x) - w_0(x)| \le a, & \forall x \in \mathbb{R}.
\end{array}
\right.
\quad (3.12) $$

where

$$ h_2^{\pm} = f_{i \pm 1/2}^{n} - h_1^{\pm}. $$

We have the following monotonicity properties for $h_j^{\pm}$, $j = 1, 2$.

**Proposition 3.1.** Given $\alpha, \beta, \gamma \in \mathbb{R}$, the following monotonicity properties hold:

* if $|\beta - \gamma| \le a$:

- $h_1^+(\alpha, \beta, \gamma)$ is non-decreasing in $\alpha$ and $\gamma$ and non-increasing in $\beta$;

- $h_2^+(\alpha, \beta, \gamma)$ is non-decreasing in $\alpha$ and $\beta$ and non-increasing in $\gamma$;

* if $|\alpha - \gamma| \le a$:

- $h_1^-(\alpha, \beta, \gamma)$ is non-decreasing in $\alpha$ and non-increasing in $\beta$ and $\gamma$;

- $h_2^-(\alpha, \beta, \gamma)$ is non-decreasing in $\gamma$ and non-increasing in $\alpha$ and $\beta$.

*Proof.* We first suppose $|\beta - \gamma| \le a$ and focus on $h_1^+$. If $\alpha \le \beta$, see Section 2.2, then no shock waves are formed so, recalling its definition (3.8), $h_1^+(\alpha, \beta, \gamma)$ is the minimum in $[\alpha, \beta]$ of $\tilde{f}_\gamma$ as for the standard Godunov flux. More explicitly,

$$ h_1^+(\alpha, \beta, \gamma) = \begin{cases} \min_{u \in [\alpha, \beta]} f(u) & \gamma - a \le \alpha, \\ \min_{u \in [\alpha, \beta]} f(u) & \gamma - a \ge \alpha \text{ and } f'(\gamma - a) \le 0, \\ \min_{u \in [\alpha, \beta]} \frac{1}{2} f(u) + \frac{1}{2} f(\gamma - a) & \gamma - a \ge \alpha \text{ and } f'(\gamma - a) \ge 0. \end{cases} \quad (3.13) $$

Indeed, if in the semi-cell $K_i^{n,l}$, $t \mapsto \tilde{u}^n(x, t)$ possibly decreases from $\beta$ to $\alpha$, hence when $\gamma - a \le \alpha$ the couple $(\tilde{u}^n(x, \cdot), \tilde{w}^n(x, \cdot))$ will always be in the interior of the hysteresis strip, so $\partial_t \tilde{w}^n = 0$ and the problem is solved as in the classical way; instead, when $\gamma - a \ge \alpha$ and $f'(\gamma - a) \le 0$, we consider $\tilde{f}_\gamma(u)$ and its minimum is in $[\gamma - a, \beta]$, hence when $\tilde{f}_\gamma(u) = f(u)$; otherwise, if $\gamma - a \ge \alpha$ and $f'(\gamma - a) \ge 0$, then the minimum is in $[\alpha, \gamma - a]$.

By looking directly at (3.13), it can be easily checked that $h_1^+$ is continuous in every variable when the other two are fixed, hence to prove the monotonicity we just need to prove it in every different possible case of (3.13), ignoring the conditions on the right-hand side of (3.13). It is straightforward then that $\alpha \mapsto h_1^+(\alpha, \beta, \gamma)$ is non-decreasing and that $\beta \mapsto h_1^+(\alpha, \beta, \gamma)$ is non-increasing. Also $\gamma \mapsto h_1^+(\alpha, \beta, \gamma)$ is non-decreasing as either $h_1^+$ is constant in $\gamma$ or non-decreasing under the condition $f(\gamma - a)' \ge 0$.
When $\alpha \ge \beta$, see Section 2.3, paying attention to the case of “fast shocks”, we have

$$h_1^+(\alpha, \beta, \gamma) = \begin{cases} \max_{u \in [\beta, \alpha]} f(u) & \gamma + a \ge \alpha, \\ \max_{u \in [\beta, \alpha]} f(u) & \gamma + a \le \alpha \text{ and } f(\alpha) \le f(\beta), \\ \mu(\alpha - \beta) + f(\beta) & \gamma + a \le \alpha, f(\alpha) \ge f(\beta) \text{ and } \mu_r \le \mu_l, \\ \max_{u \in [\beta, \alpha]} \frac{1}{2}f(u) + \frac{1}{2}f(\gamma + a) & \gamma + a \le \alpha, f(\alpha) > f(\beta) \text{ and } \mu_r \ge \mu_l, \end{cases} \quad (3.14)$$

where $\mu = \mu(\alpha, \beta, \gamma)$. Again it can be checked that $h_1^+$ is continuous in each variable when we fix the another two. In particular, notice that, when $f(\alpha) = f(\beta)$, then second and third definitions coincide and, when $\mu_r = \mu_l$, then the last two coincide (see Remark 2.4). Then again we check the monotonicity for every possible case in (3.14). The first two are constant in $\gamma$ and trivially non-decreasing in $\alpha$ and non-increasing in $\beta$. For the third one instead, using the notation introduced in Subcase 2.3.2 with $u_l = \alpha, u_r = \beta$ and $w_r = \gamma$, we have that since $\partial_\alpha I_r = 0, \partial_\alpha I_l = 1, \partial_\beta I_r = -1, \partial_\beta I_l = 0, \partial_\gamma I_r = 1, \partial_\gamma I_l = -1$, and $\alpha - \beta = I_r + I_l$ then

$$\frac{\partial (\mu(\alpha - \beta) + f(\beta))}{\partial \alpha} = \frac{f'(\alpha)(I_r + 2I_l) - 2(f(\alpha) - f(\beta))}{(I_r + 2I_l)^2} (\alpha - \beta) + \mu \quad (3.15) \\ = \frac{I_r + I_l}{I_r + 2I_l} (f'(\alpha) - 2\mu) + \mu,$$

$$\frac{\partial (\mu(\alpha - \beta) + f(\beta))}{\partial \beta} = \frac{-f'(\beta)(I_r + 2I_l) + (f(\alpha) - f(\beta))}{(I_r + 2I_l)^2} (\alpha - \beta) - \mu + f'(\beta) \\ = \frac{I_r + I_l}{I_r + 2I_l} (\mu - f'(\beta)) - (\mu - f'(\beta)) \quad (3.16) \\ = \frac{I_l}{I_r + 2I_l} (f'(\beta) - \mu),$$

$$\frac{\partial (\mu(\alpha - \beta) + f(\beta))}{\partial \gamma} = \frac{f(\alpha) - f(\beta)}{(I_r + 2I_l)^2} (\alpha - \beta) \quad (3.17) \\ = \frac{I_r + I_l}{I_r + 2I_l} \mu.$$

We observe that (3.15) is positive because $f'(\alpha) - 2\mu \ge f'(\alpha) - 2\mu_l \ge 0$ by $\mu_r \le \mu \le \mu_l$ and the convexity of $f$, and because $\mu \ge 0$ since $f(\alpha) \ge f(\beta)$; (3.16) is negative since for the same previous reasons $f'(\beta) - \mu \le f'(\beta) - \mu_r \le 0$; (3.17) positive since again $\mu \ge 0$. Regarding the fourth case in (3.14), it is non-decreasing in $\alpha$ and non-increasing in $\beta$; it is also non-decreasing in $\gamma$ since $\mu_r \ge \mu_l$ implies $\mu_r \ge 0$, hence $f'(\gamma + a) \ge 0$ by convexity of $f$.

Regarding, $h_2^+$, as $h_2^+ = f_{i-1/2}^n - h_1^+$ and by recalling that $f_{i-1/2}^n = g(\alpha, \beta)$, see (3.3), we have that when $\alpha \le \beta$

$$h_2^+(\alpha, \beta, \gamma) = \begin{cases} 0 & \gamma - a \le \alpha, \\ 0 & \gamma - a \ge \alpha \text{ and } f'(\gamma - a) \le 0, \\ \min_{u \in [\alpha, \beta]} \frac{1}{2}f(u) - \frac{1}{2}f(\gamma - a) & \gamma - a \ge \alpha \text{ and } f'(\gamma - a) \ge 0. \end{cases} \quad (3.18)$$

which, when the other two variables are fixed, is again non-decreasing in $\alpha$, non-increasing in $\gamma$ and constant in $\beta$. Indeed, notice that when $\gamma - a \ge \alpha$ and $f'(\gamma - a) \ge 0$, then the
minimum of $1/2f$ is obtained in $[\alpha, \gamma - a]$ since $f$ is convex and $\beta \ge \gamma - a$, so in this case $h_2^+$ is constant in $\beta$. Similarly, when $\alpha \ge \beta$,

$$h_2^+(\alpha, \beta, \gamma) = \begin{cases} 0 & \gamma + a \ge \alpha, \\ 0 & \gamma + a \le \alpha \text{ and } f(\alpha) \le f(\beta), \\ I_l \mu & \gamma + a \le \alpha, f(\alpha) \ge f(\beta) \text{ and } \mu_r \le \mu_l, \\ \max_{u \in [\beta, \alpha]} \frac{1}{2} f(u) - \frac{1}{2} f(\gamma + a) & \gamma + a \le \alpha, f(\alpha) > f(\beta) \text{ and } \mu_r \ge \mu_l, \end{cases} \quad (3.19)$$

as for the case $\gamma + a \le \alpha, f(\alpha) \ge f(\beta)$ and $\mu_r \le \mu_l$ we have $f_{i-1/2}^n = f(\alpha)$, hence

$$f(\alpha) - (\mu(\alpha - \beta) + f(\beta)) = f(\alpha) - f(\beta) - \frac{I_r + I_l}{I_r + 2I_l} (f(\alpha) - f(\beta)) = I_l \mu.$$

We infer the same conclusion as

$$\begin{aligned} \frac{\partial (I_l \mu)}{\partial \alpha} &= \mu + \frac{f'(\alpha)(I_r + 2I_l) - 2(f(\alpha) - f(\beta))}{(I_r + 2I_l)^2} I_l \\ &= \mu + \frac{I_l}{I_r + 2I_l} (f'(\alpha) - 2\mu) \ge 0, \end{aligned} \quad (3.20)$$

$$\begin{aligned} \frac{\partial (I_l \mu)}{\partial \beta} &= \frac{-f'(\beta)(I_r + 2I_l) + (f(\alpha) - f(\beta))}{(I_r + 2I_l)^2} I_l \\ &= -\frac{I_l}{I_r + 2I_l} (f'(\beta) - \mu) \ge 0, \end{aligned} \quad (3.21)$$

$$\begin{aligned} \frac{\partial (I_l \mu)}{\partial \gamma} &= -\mu + \frac{f(\alpha) - f(\beta)}{(I_r + 2I_l)^2} I_l = -\mu + \frac{I_l}{I_r + 2I_l} \mu \quad (3.22) \\ &= -\frac{I_r + I_l}{I_r + 2I_l} \mu \le 0, \end{aligned}$$

and in the case when $\gamma + a \le \alpha, f(\alpha) > f(\beta)$ and $\mu_r \ge \mu_l$, the maximum of $1/2f(u)$ equals to $1/2f(\alpha)$ which is constant in $\beta$.

Briefly, if we now consider $h_1^-(\alpha, \beta, \gamma)$ with $|\alpha - \gamma| \le a$ then if $\alpha \le \beta$

$$h_1^-(\alpha, \beta, \gamma) = \begin{cases} \min_{u \in [\alpha, \beta]} f(u) & \gamma + a \ge \beta, \\ \min_{u \in [\alpha, \beta]} f(u) & \gamma + a \ge \beta \text{ and } f'(\gamma + a) \ge 0, \\ \min_{u \in [\alpha, \beta]} \frac{1}{2} f(u) + \frac{1}{2} f(\gamma + a) & \gamma + a \ge \beta \text{ and } f'(\gamma + a) \le 0. \end{cases} \quad (3.23)$$

which again is non-decreasing in $\alpha$, non-increasing in $\beta$ and this time also non-increasing in $\gamma$. The same conclusion holds, for $\alpha \ge \beta$ when

$$h_1^+(\alpha, \beta, \gamma) = \begin{cases} \max_{u \in [\beta, \alpha]} f(u) & \gamma - a \le \beta, \\ \max_{u \in [\beta, \alpha]} f(u) & \gamma - a \ge \beta \text{ and } f(\alpha) \ge f(\beta), \\ -\nu(\alpha - \beta) + f(\alpha) & \gamma - a \ge \beta, f(\alpha) \le f(\beta) \text{ and } \nu_l \ge \nu_r, \\ \max_{u \in [\beta, \alpha]} \frac{1}{2} f(u) + \frac{1}{2} f(\gamma - a) & \gamma - a \ge \beta, f(\alpha) < f(\beta) \text{ and } \nu_l \le \nu_r. \end{cases} \quad (3.24)$$
We just point out that

$$ \frac{\partial (-\nu(\alpha - \beta) + f(\alpha))}{\partial \alpha} = \frac{J_r}{J_l + 2J_r} (f'(\alpha) - \nu) \geq 0, \qquad (3.25) $$

$$ \frac{\partial (-\nu(\alpha - \beta) + f(\alpha))}{\partial \beta} = \frac{J_l + J_r}{J_l + 2J_r} (f'(\beta) - 2\nu) + \nu \leq 0, \qquad (3.26) $$

$$ \frac{\partial (-\nu(\alpha - \beta) + f(\alpha))}{\partial \gamma} = \frac{J_l + J_r}{J_l + 2J_r} \nu \leq 0. \qquad (3.27) $$

Regarding $h_2^-$ a similar analysis as for $h_2^+$ can be done, implying the desired conclusions. □

**Proposition 3.2.** $h_j^\pm$ are L-Lipschitz continuous with respect to the first two arguments for $j = 1, 2$. That is

$$ |h_j^\pm(\alpha_1, \beta, \gamma) - h_j^\pm(\alpha_1, \beta, \gamma)| \leq L|\alpha_1 - \alpha_2| $$

and

$$ |h_j^\pm(\alpha, \beta_1, \gamma) - h_j^\pm(\alpha, \beta_2, \gamma)| \leq L|\beta_1 - \beta_2| $$

for $j = 1, 2$ and for all $\alpha, \alpha_1, \alpha_2$ and $\beta, \beta_1, \beta_2$ compatible with the fixed $\gamma$.

*Proof.* Since $h_j^\pm(\alpha, \beta, \gamma)$ are continuous in $\alpha$ and $\beta$, then it is sufficient to check the Lipschitz continuity in each different case of definition. In particular, as $\gamma$ is fixed, it is straightforward to see that $\min_{[\alpha, \beta]} f(u)$ and $\max_{[\beta, \alpha]} f(u)$ are L-Lipschitz continuous with respect to $\alpha$ and $\beta$, so we just prove explicitly the statement in the case of “fast shocks”. Let us now focus on this case for $h_1^+$; then by looking at (3.15) and (3.16) we have that

$$ 0 \leq \frac{\partial (\mu(\alpha - \beta) + f(\beta))}{\partial \alpha} \leq \sup_{[U_m, U_M]} |f'(u)| = L, $$

$$ 0 \geq \frac{\partial (\mu(\alpha - \beta) + f(\beta))}{\partial \beta} \geq -\frac{1}{2} \left( \sup_{[U_m, U_M]} |f'(u)| + \mu \right) \geq -L, $$

showing the Lipschitz continuity. Similar estimates also holds in the other cases for all $h_j^\pm$.

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

A direct consequence of Proposition 3.1 and Proposition 3.2 is the following corollary.

**Corollary 3.3.** Under the CFL condition (3.1), we can write

$$ u_i^{n+1} = H_1(u_{i-1}^n, u_i^n, u_{i+1}^n, w_i^n), \qquad (3.28) $$

$$ w_i^{n+1} = H_2(u_{i-1}^n, u_i^n, u_{i+1}^n, w_i^n), \qquad (3.29) $$

with $H_1, H_2$ non-decreasing in each of their arguments and such that

$$ \begin{aligned} H_1(u_{i-1}^n, u_i^n, u_{i+1}^n, w_i^n) + H_2(u_{i-1}^n, u_i^n, u_{i+1}^n, w_i^n) &= u_i^n + w_i^n - \frac{\Delta t}{\Delta x} (f_{i+1/2}^n - f_{i-1/2}^n) \\ &= u_i^n + w_i^n - \frac{\Delta t}{\Delta x} (g(u_i^n, u_{i+1}^n) - g(u_i^n, u_{i+1}^n)). \end{aligned} \qquad (3.30) $$
Proof. By Proposition 3.2, it is clear that

$$H_1(u_{i-1}^n, u_i^n, u_{i+1}^n, w_i^n) = u_i^n - \frac{\Delta t}{\Delta x} (h_1^-(u_i^n, u_{i+1}^n, w_i^n) - h_1^+(u_{i-1}^n, u_i^n, w_i^n)),$$

is non-decreasing in $u_{i-1}^n, u_{i+1}^n, w_i^n$. Fix now $u_{i-1}^n, u_{i+1}^n, w_i^n$ and consider $\alpha_1 \le \alpha_2$, then, by Proposition 3.2 and under the CFL condition (3.1)

$$ 
\begin{aligned}
H_1(u_{i-1}^n, \alpha_2, u_{i+1}^n, w_i^n) - H_1(u_{i-1}^n, \alpha_1, u_{i+1}^n, w_i^n) &= \\
&= \alpha_2 - \alpha_1 - \frac{\Delta t}{\Delta x} (h_1^-(\alpha_2, u_{i+1}^n, w_i^n) - h_1^-(\alpha_1, u_{i+1}^n, w_i^n)) \\
&\quad + \frac{\Delta t}{\Delta x} (h_1^+(u_{i-1}^n, \alpha_2, w_i^n) - h_1^+(u_{i-1}^n, \alpha_1, w_i^n)) \\
&\ge \alpha_2 - \alpha_1 - 2\frac{\Delta t}{\Delta x} L(\alpha_2 - \alpha_1) \ge 0.
\end{aligned}
$$

Similarly,

$$H_2(u_{i-1}^n, u_i^n, u_{i+1}^n, w_i^n) = w_i^n - \frac{\Delta x}{\Delta t} (h_2^-(u_i^n, u_{i+1}^n, w_i^n) - h_2^+(u_{i-1}^n, u_i^n, w_i^n)),$$

is non-decreasing $u_{i-1}^n, u_i^n, u_{i+1}^n$ as a direct consequence of Proposition 3.1, and, reasoning in the same way as for $u$, it can be proved that it is non-decreasing also in $w_i^n$ by Proposition 3.2 and the CFL condition (3.1). $\square$

## 3.2 Compactness estimates

Given $u_{i-1}^n, u_i^n, u_{i+1}^n$ and $w_i^n$ we define

$$a_{i-1/2}^n := \begin{cases} \frac{\Delta t h_1^+(u_{i-1}^n, u_i^n, w_i^n) - f(u_i^n)}{\Delta x (u_{i-1}^n - u_i^n)} & \text{if } u_i^n \neq u_{i-1}^n, \\ 0 & \text{if } u_i^n = u_{i-1}^n, \end{cases} \quad (3.31)$$

and

$$b_{i+1/2}^n := \begin{cases} \frac{\Delta t h_1^-(u_i^n, u_{i+1}^n, w_i^n) - f(u_i^n)}{\Delta x (u_i^n - u_{i+1}^n)} & \text{if } u_i^n \neq u_{i+1}^n, \\ 0 & \text{if } u_i^n = u_{i+1}^n. \end{cases} \quad (3.32)$$

So, we can rewrite the scheme for $u_i^{n+1}$ in the two following ways

$$u_i^{n+1} = u_i^n + b_{i+1/2}^n(u_{i+1}^n - u_i^n) + a_{i-1/2}^n(u_{i-1}^n - u_i^n) \quad (3.33a)$$

$$= u_i^n (1 - b_{i+1/2}^n - a_{i-1/2}^n) + b_{i+1/2}^n u_{i+1}^n + a_{i-1/2}^n u_{i-1}^n, \quad (3.33b)$$

as $h_1^+(u_{i-1}^n, u_i^n, w_i^n) = f(u_i^n)$ when $u_{i-1}^n = u_i^n$ and $h_1^-(u_i^n, u_{i+1}^n, w_i^n) = f(u_i^n)$ when $u_i^n = u_{i+1}^n$. Similarly, define

$$c_{i-1/2}^n := \begin{cases} \frac{\Delta t h_2^+(u_{i-1}^n, u_i^n, w_i^n)}{\Delta x (w_{i-1}^n - w_i^n)} & \text{if } w_i^n \neq w_{i-1}^n, \\ 0 & \text{if } w_i^n = w_{i-1}^n, \end{cases} \quad (3.34)$$

and

$$d_{i+1/2}^n := \begin{cases} \frac{\Delta t h_2^-(u_i^n, u_{i+1}^n, w_i^n)}{\Delta x (w_i^n - w_{i+1}^n)} & \text{if } w_i^n \neq w_{i+1}^n, \\ 0 & \text{if } w_i^n = w_{i+1}^n. \end{cases} \quad (3.35)$$
Notice now that if $w_{i-1}^n = w_i^n$ then $|w_i^n - u_{i-1}^n| < a$, so when solving the Riemann problem at $x_{i-1/2}$, it can be checked that $w$ remains constant in $K_i^{n,l}$, since the couple $(\tilde{u}^n, \tilde{w}^n)$ remains into the hysteresis region $\mathcal{L}$. This implies then that $h_1^+(u_{i-1}^n, u_i^n, w_i^n) = f_{i-1/2}^n$ and $h_2^+(u_{i-1}^n, u_i^n, w_i^n) = 0$. Similarly, when $w_i^n = w_{i+1}^n$ then $h_2^-(u_i^n, u_{i+1}^n, w_i^n) = 0$. With this in mind, the scheme for $w$ can be rewritten as

$$w_i^{n+1} = w_i^n + d_{i+1/2}^n(w_{i+1}^n - w_i^n) + c_{i-1/2}^n(w_{i-1}^n - w_i^n) \quad (3.36a)$$

$$= w_i^n (1 - d_{i+1/2}^n + c_{i-1/2}^n) + d_{i+1/2}^n w_{i+1}^n + c_{i-1/2}^n w_{i-1}^n. \quad (3.36b)$$

**Lemma 3.4.** Under the CFL condition (3.1), for all $i \in \mathbb{Z}$ and $n \in \mathbb{N}$, it holds that $0 \le a_{i-1/2}^n, b_{i+1/2}^n \le 1/2$ and $0 \le c_{i-1/2}^n, d_{i+1/2}^n \le 1/2$.

*Proof.* Positivity of $a_{i-1/2}^n$ and $b_{i+1/2}^n$ follows from the monotonicity properties of $h_1^\pm$, see Proposition 3.1, and from the fact that $h_1^\pm(u_i^n, u_i^n, w_i^n) = f(u_i^n)$. The estimate $a_{i-1/2}^n, b_{i+1/2}^n \le 1/2$ follows from the Lipschitz continuity of $h_1^\pm$, see Proposition 3.2, together with the CFL condition (3.1).

Regarding $c_{i-1/2}^n$ it is easy to see that if $|u_{i-1}^n - w_i^n| \le a$ then $h_2^+ = 0$, see indeed (3.18) and (3.19). Suppose then for example that $u_{i-1}^n < w_i^n - a$, which in particular implies that $0 \le w_i^n - w_{i-1}^n \le u_i^n - u_{i-1}^n$ as $w_{i-1}^n - a \le u_{i-1}^n < w_i^n - a \le u_i^n$. Hence, as $h_2^+(u_i^n, u_i^n, w_i^n) = 0$, then we can use the monotonicity properties of Proposition 3.1, the Lipschitz continuity stated in Proposition 3.2 and the CFL condition (3.1) to conclude that

$$0 \le c_{i-1}^n = \frac{\Delta t h_2^+(u_{i-1}^n, u_i^n, w_i^n) - h_2^+(u_i^n, u_i^n, w_i^n)}{\Delta x} \frac{u_{i-1}^n - u_i^n}{u_{i-1}^n - u_i^n} \le \frac{1}{2}.$$

With similar reasoning we can prove the estimate when $u_{i-1}^n > w_i^n + a$, noticing that is such case $0 \le w_i^n - w_{i-1}^n \le u_{i-1}^n - u_i^n$. Also the estimates for $d_{i+1/2}^n$ can be deduced in the same way as the ones for $c_{i-1/2}^n$. $\square$

The discrete BV estimate in space then follows.

**Lemma 3.5.** Assume, $u_0, w_0 \in \text{BV}(\mathbb{R})$, $|u_0(x) - w_0(x)| \le a$ and let $\Delta x, \Delta t > 0$ such that the CFL condition (3.1) holds. Then the following inequalities hold

$$\sum_{i \in \mathbb{Z}} |u_{i+1}^{n+1} - u_i^{n+1}| \le \sum_{i \in \mathbb{Z}} |u_{i+1}^n - u_i^n|, \quad \forall n \in \mathbb{N}, \quad (3.37)$$

and

$$\sum_{i \in \mathbb{Z}} |w_{i+1}^{n+1} - w_i^{n+1}| \le \sum_{i \in \mathbb{Z}} |w_{i+1}^n - w_i^n|, \quad \forall n \in \mathbb{N}. \quad (3.38)$$

*Proof.* By using (3.33a) we can write

$$u_{i+1}^{n+1} - u_i^{n+1} = (u_{i+1}^n - u_i^n) (1 - b_{i+1/2}^n - a_{i+1/2}^n) + b_{i+3/2}^n (u_{i+2}^n - u_{i+1}^n) + a_{i-1/2}^n (u_i^n - u_{i-1}^n)$$

Then by Lemma 3.4

$$b_{i+3/2}^n, a_{i-1/2}^n, (1 - b_{i+1/2}^n - a_{i+1/2}^n) \ge 0$$

so

$$|u_{i+1}^{n+1} - u_i^{n+1}| \le |u_{i+1}^n - u_i^n| (1 - b_{i+1/2}^n - a_{i+1/2}^n) + b_{i+3/2}^n |u_{i+1}^n - u_{i+1}^n| + a_{i-1/2}^n |u_i^n - u_{i-1}^n|.$$
Summing this last inequality over $i \in \mathbb{Z}$ we get (3.37).

Starting instead from (3.36a) and by following the same previous reasoning we infer (3.38). □

Notice that the assumption $u_0, w_0 \in BV(\mathbb{R})$ is not really needed, however it ensures that the sums in the inequalities (3.37) and (3.38) are finite.

Let us now denote as $u^\Delta$ and $w^\Delta$ the approximate finite volume solution defined by

$$u^{\Delta}(x,t) := u_{i}^{n}, \quad w^{\Delta}(x,t) := w_{i}^{n} \quad \text{for } x \in [x_{i-1/2}, x_{i+1/2}], \ t \in [t^{n}, t^{n+1}]. \qquad (3.39)$$

We then recover the following bounds on the $L^\infty$ norm and on the total variation in space-time for both $u^\Delta$ and $w^\Delta$, as a consequence of the previous lemmas. The proof of the following results are quite standard for the case of a single unknown $u$, see e.g. [13, Section 5.3], and can be directly adapted to the variable $w$ as it shares the same properties of $u$.

**Proposition 3.6.** Let $u_0, w_0 \in L^\infty(\mathbb{R})$, such that $u_0(x) \in [U_m, U_M]$, $w_0(x) \in [W_m, W_M]$ and $|u_0(x) - w_0(x)| \le a$ for almost every $x \in \mathbb{R}$. Fix also $\Delta x, \Delta t$ such that the CFL condition (3.1) holds and consider the approximate solutions $u^\Delta$ and $w^\Delta$ defined by (3.12), (3.39). Then $u^\Delta$ and $w^\Delta$ satisfy

$$u^{\Delta}(x,t) \in [U_m, U_M] \quad \text{and} \quad w^{\Delta}(x,t) \in [W_m, W_M], \qquad (3.40)$$

for a.e. $(x, t) \in \mathbb{R} \times [0, T]$.

**Proposition 3.7.** Let $u_0, w_0 \in BV(\mathbb{R})$ with $|u_0(x) - w_0(x)| \le a$ for almost all $x \in \mathbb{R}$, and fix $\Delta x, \Delta t > 0$ such that the CFL condition (3.1) holds. Consider the approximate solutions $u^\Delta$ and $w^\Delta$ defined by (3.12) and (3.39). Then, for any $T > 0$, there exist $C_1 = C_1(u_0, T) > 0$ and $C_2 = C_2(w_0, T) > 0$ such that

$$|u^{\Delta}|_{BV(\mathbb{R} \times [0, T])} \le C_1 \quad \text{and} \quad |w^{\Delta}|_{BV(\mathbb{R} \times [0, T])} \le C_2. \qquad (3.41)$$

In particular, the proof of Proposition 3.7 relies on the following $L^1$ continuity in time.

**Lemma 3.8.** Given $w_0, u_0 \in BV(\mathbb{R})$, we have that

$$\int_{\mathbb{R}} |u^{\Delta}(x, s) - u^{\Delta}(x, t)| dx \le 2L |u_0|_{BV(\mathbb{R})} (|t - s| + \Delta t) \qquad (3.42)$$

and

$$\int_{\mathbb{R}} |w^{\Delta}(x, s) - w^{\Delta}(x, t)| dx \le 2L |w_0|_{BV(\mathbb{R})} (|t - s| + \Delta t) \qquad (3.43)$$

for any $s, t \in [0, T[$.

*Proof.* If $s, t \in [t^n, t^{n+1}]$ for some $n \in \mathbb{N}$, then $u^\Delta(x, s) \equiv u^\Delta(x, t)$ and $w^\Delta(x, s) \equiv w^\Delta(x, t)$ so the two inequalities are satisfied trivially.

Let us assume that $s \in [t^{n_1}, t^{n_1+1}]$, $t \in [t^{n_2}, t^{n_2+1}]$ with e.g. $n_1 < n_2$. From (3.33a) and the estimates of Lemma 3.4, under the CFL condition (3.1), it can be shown that

$$\Delta x \sum_{i \in \mathbb{Z}} |u_i^{n+1} - u_i^n| \le 2L \Delta t |u_0|_{BV(\mathbb{R})},$$
for any $n \in \mathbb{N}$ (for the details see [13, proof of Corollary 5.1]). Consequently,

$$
\begin{align*}
\int_{\mathbb{R}} |u^{\Delta}(x, s) - u^{\Delta}(x, t)| dx &\le \sum_{n=n_1}^{n_2-1} \int_{\mathbb{R}} |u^{\Delta}(x, t^{n+1}) - u^{\Delta}(x, t^n)| dx \\
&= \Delta x \sum_{n=n_1}^{n_2-1} \sum_{i \in \mathbb{Z}} |u_i^{n+1} - u_i^n| \le 2L \sum_{n=n_1}^{n_2-1} \Delta t |u_0|_{\text{BV}(\mathbb{R})} \\
&= 2L (n_2 - n_1) \Delta t |u_0|_{\text{BV}(\mathbb{R})} \le 2L (|s-t| + \Delta t) |u_0|_{\text{BV}(\mathbb{R})}.
\end{align*}
$$

By the same reasoning, we get (3.43). □

3.3 Discrete entropy condition

Recalling the definition of g, that is (3.3) and using the monotonicity properties of Corol-
lary 3.3, we can prove the following discrete entropy condition.

**Proposition 3.9.** For each $i \in \mathbb{Z}$, $n \in \mathbb{N}$, $(k, \hat{k}) \in \mathcal{L}$, it holds

$$
\begin{align*}
& |u_i^{n+1} - k| - |u_i^n - k| + |w_i^{n+1} - \hat{k}| - |w_i^n - \hat{k}| + \\
& \frac{\Delta t}{\Delta x} \biggl( g(u_i^n \top k, u_{i+1}^n \top k) - g(u_i^n \perp k, u_{i+1}^n \perp k) \\
& \phantom{\frac{\Delta t}{\Delta x} \biggl(} -g(u_{i-1}^n \top k, u_i^n \top k) + g(u_{i-1}^n \perp k, u_i^n \perp k) \biggr) \le 0, \quad (3.44)
\end{align*}
$$

where $a \top b := \max\{a, b\}$ and $a \perp b := \min\{a, b\}$.

Proof. It is easy to check that, since $|u_i^n - w_i^n| \le a$ and $|k - \hat{k}| \le a$, then $|u_i^{n\top}k - w_i^{n\top}\hat{k}| \le a$ and $|u_i^n \perp k - w_i^n \perp \hat{k}| \le a$. Then by Corollary 3.3 we have that

$$
u_i^{n+1} \le H_1(u_{i-1}^n \top k, u_i^n \top k, u_{i+1}^n \top k, w_i^n \top \hat{k}), \quad w_i^{n+1} \le H_2(u_{i-1}^n \top k, u_i^n \top k, u_{i+1}^n \top k, w_i^n \top \hat{k}), \tag{3.45}
$$

and

$$
k \le H_1(u_{i-1}^n \top k, u_i^n \top k, u_{i+1}^n \top k, w_i^n \top \hat{k}), \quad \hat{k} \le H_2(u_{i-1}^n \top k, u_i^n \top k, u_{i+1}^n \top k, w_i^n \top \hat{k}), \quad (3.46)
$$

as $k = H_1(k, k, k, \hat{k})$ and $\hat{k} = H_2(k, k, k, \hat{k})$. Hence, adding $u_i^{n+1} \top k$ and $w_i^{n+1} \top \hat{k}$, thanks
to (3.30), (3.45) and (3.46), we get the following estimate

$$
u_i^{n+1} \top k + w_i^{n+1} \top \hat{k} \le u_i^n \top k + w_i^n \top \hat{k} - \frac{\Delta t}{\Delta x} \left( g(u_i^n \top k, u_{i+1}^n \top k) - g(u_{i-1}^n \top k, u_i^n \top k) \right). \quad (3.47)
$$

Similarly,

$$
u_i^{n+1} \perp k + w_i^{n+1} \perp \hat{k} \ge u_i^n \perp k + w_i^n \perp \hat{k} - \frac{\Delta t}{\Delta x} \left( g(u_i^n \perp k, u_{i+1}^n \perp k) - g(u_{i-1}^n \perp k, u_i^n \perp k) \right), \quad (3.48)
$$

so subtracting (3.48) to (3.47) we get (3.44). □
## 3.4 Weak Hysteresis Relationship

Our goal is then to show that the approximate solutions constructed via the numerical scheme presented in Section 3.1 satisfy the weak hysteresis relationship, so by starting from (1.9) we would like to prove (1.8). As shock waves solution lack in regularity, we cannot do this directly. Instead, we will exploit the following property.

**Lemma 3.10.** Consider the scheme (3.12), with $\Delta x, \Delta t > 0$ such that the CFL condition (3.1) holds. Moreover, define

$$G(u) := \int_{0}^{u} \xi f'(\xi) d\xi = uf(u) - \int_{0}^{u} f(\xi) d\xi$$

and denote by $\tilde{u}^n$ and $\tilde{w}^n$ the exact solutions defined on $\mathbb{R} \times [t^n, t^{n+1}]$ such that $\tilde{u}^n(x, t^n) = u_i^n$ and $\tilde{w}^n(x, t^n) = w_i^n$ for $x \in K_i$. Then for $t \in [t^n, t^{n+1}]$, if there is an entropic shock in the half cell $K_i^{n,l}$, it holds

$$-\Delta t (G(u_i^n) - G(\tilde{u}(x_{i-1/2}^+, t))) - \frac{1}{2} \int_{x_{i-1/2}}^{x_i} (\tilde{u}^n(x, t^{n+1})^2 - \tilde{u}^n(x, t^n)^2) dx \\ - \frac{1}{2} \int_{x_{i-1/2}}^{x_i} (\tilde{w}^n(x, t^{n+1})^2 - \tilde{w}^n(x, t^n)^2) dx \geq a |\partial_t \tilde{w}^n|(K_i^{n,l}); \quad (3.49)$$

instead, if an entropic shock is present in $K_i^{n,r}$, it holds

$$-\Delta t (G(\tilde{u}(x_{i+1/2}^-, t)) - G(u_i^n)) - \frac{1}{2} \int_{x_i}^{x_{i+1/2}} (\tilde{u}^n(x, t^{n+1})^2 - \tilde{u}^n(x, t^n)^2) dx \\ - \frac{1}{2} \int_{x_i}^{x_{i+1/2}} (\tilde{w}^n(x, t^{n+1})^2 - \tilde{w}^n(x, t^n)^2) dx \geq a |\partial_t \tilde{w}^n|(K_i^{n,r}). \quad (3.50)$$

*Proof.* We detail the proof for (3.49), the other case being similar. Three configuration can occur: a single shock in $\tilde{u}^n$ with $\tilde{w}^n$ constant; two shocks, one both in $\tilde{u}^n$ and $\tilde{w}^n$ and one only in $\tilde{u}^n$; a “fast shock” in both $\tilde{u}^n$ and $\tilde{w}^n$; see Case 2.3 and in particular Subcase 2.3.2 in Section 2.

* If there is only a shock in $u$, then the right hand side of (3.49) is 0. Since the shock joins the states $u_{i-1}^n$ and $u_i^n$ with speed given by the Rankine-Hugoniot condition, see Figure 11, the left hand side reads

$$-\Delta t (G(u_i^n) - G(\tilde{u}(x_{i-1/2}^+, t))) - \frac{1}{2} \int_{x_{i-1/2}}^{x_i} (\tilde{u}^n(x, t^{n+1})^2 - \tilde{u}^n(x, t^n)^2) dx \\ - \frac{1}{2} \int_{x_{i-1/2}}^{x_i} (\tilde{w}^n(x, t^{n+1})^2 - \tilde{w}^n(x, t^n)^2) dx \\ = - \Delta t (G(u_i^n) - G(u_{i-1}^n)) - \frac{1}{2} \frac{f(u_{i-1}^n) - f(u_i^n)}{u_{i-1}^n - u_i^n} \Delta t ((u_{i-1}^n)^2 - (u_i^n)^2) \\ = - \Delta t \left[ (G(u_i^n) - G(u_{i-1}^n)) + \frac{1}{2} (f(u_{i-1}^n) - f(u_i^n)) (u_{i-1}^n + u_i^n) \right]. \quad (3.51)$$
<figure><img src="image_14.png" /><figcaption>Figure 11: The exact solution in the half-cell $K_i^{n,l}$ in the case only a shock in $\tilde{u}^n$ is present (here represented by the dashed line) and it has speed $s$ given by the Rankine-Hugoniot condition.</figcaption></figure>

Replacing $G$ by its definition and doing some computations we get

$$ [(3.51)] = -\Delta t \left[ u_i^n f(u_i^n) - u_{i-1}^n f(u_{i-1}^n) + \int_{u_i^n}^{u_{i-1}^n} f(\xi) \, d\xi + \frac{1}{2} (f(u_{i-1}^n) - f(u_i^n)) (u_{i-1}^n + u_i^n) \right] \\ = \Delta t \left[ \frac{1}{2} (f(u_{i-1}^n) + f(u_i^n)) (u_{i-1}^n - u_i^n) - \int_{u_i^n}^{u_{i-1}^n} f(\xi) \, d\xi \right] \geq 0 $$

since $u_{i-1}^n > u_i^n$ (entropic shock) and $f$ is convex.

* If there are two shocks, one in $\tilde{u}^n$ only, connecting $u_i^n$ to $w_i^n + a$, and the other in both $\tilde{u}^n$ and $\tilde{w}^n$, which connects respectively $w_i^n + a$ to $u_i^n$ and $w_i^n$ to $u_{i-1}^n - a$, both with slope given by the Rankine-Hugoniot condition (see Figure 12), then

$$ \begin{aligned} & \int_{x_{i-1/2}}^{x_i} (\tilde{u}^n(x, t^{n+1})^2 - \tilde{u}^n(x, t^n)^2) \, dx \\ &= \Delta x_1 ((u_{i-1}^n)^2 - (u_i^n)^2) + (\Delta x_2 - \Delta x_1) ((w_i^n + a)^2 - (u_i^n)^2) \\ &= \Delta x_1 ((u_{i-1}^n)^2 - (w_i^n + a)^2) + \Delta x_2 ((w_i^n + a)^2 - (u_i^n)^2) \\ &= \frac{1}{2} \Delta t \frac{f(u_{i-1}^n) - f(w_i^n + a)}{u_{i-1}^n - (w_i^n + a)} ((u_{i-1}^n)^2 - (w_i^n + a)^2) \quad (3.52) \\ & \quad + \Delta t \frac{f(w_i^n + a) - f(u_i^n)}{(w_i^n + a) - u_i^n} ((w_i^n + a)^2 - (u_i^n)^2) \\ &= \frac{1}{2} \Delta t (f(u_{i-1}^n) - f(w_i^n + a)) (u_{i-1}^n + w_i^n + a) \\ & \quad + \Delta t (f(w_i^n + a) - f(u_i^n)) (w_i^n + a + u_i^n) \end{aligned} $$
<figure><img src="image_15.png" /><figcaption>Figure 12: The exact solution in the semi-cell $K_i^{n,l}$ in the case when we have two shocks.</figcaption></figure>

and

$$ 
\begin{aligned}
& \int_{x_{i-1/2}}^{x_i} (\tilde{w}^n(x, t^{n+1})^2 - \tilde{w}^n(x, t^n)^2) \, dx \\
&= \frac{1}{2} \Delta t \frac{f(u_{i-1}^n) - f(w_i^n + a)}{(u_{i-1}^n) - (w_i^n + a)} ((u_{i-1}^n - a)^2 - (w_i^n)^2) \\
&= \frac{1}{2} \Delta t (f(u_{i-1}^n) - f(w_i^n + a)) (u_{i-1}^n - a + w_i^n) \\
&= \frac{1}{2} \Delta t (f(u_{i-1}^n) - f(w_i^n + a)) (u_{i-1}^n + a + w_i^n) - a \Delta t (f(u_{i-1}^n) - f(w_i^n + a))
\end{aligned} 
\quad (3.53) $$

and

$$ 
\begin{aligned}
a |\partial_t \tilde{w}^n| (K_i^{n,l}) &= a \Delta t \frac{1}{2} \frac{f(u_{i-1}^n) - f(w_i^n + a)}{u_{i-1}^n - (w_i^n + a)} |(u_{i-1}^n - a) - w_i^n| \\
&= \frac{1}{2} a \Delta t (f(u_{i-1}^n) - f(w_i^n + a)).
\end{aligned} 
\quad (3.54) $$

Replacing (3.52), (3.53) and (3.54) in (3.49), it becomes

$$ -\Delta t \left[ (G(u_i^n) - G(u_{i-1}^n)) + \frac{1}{2} (f(u_{i-1}^n) - f(w_i^n + a)) (u_{i-1}^n + w_i^n + a) + \frac{1}{2} (f(w_i^n + a) - f(u_i^n)) (w_i^n + a + u_i^n) \right] \geq 0, $$

which can be rewritten as

$$ -\Delta t \left[ \int_{w_i^n+a}^{u_{i-1}^n} f(\xi) \, d\xi - \frac{1}{2} (f(u_{i-1}^n) + f(w_i^n + a)) (u_{i-1}^n - (w_i^n + a)) + \int_{u_i^n}^{w_i^n+a} f(\xi) \, d\xi - \frac{1}{2} (f(w_i^n + a) + f(u_i^n)) (w_i^n + a - u_i^n) \right] \geq 0. $$
As before, since $u_{i-1}^n \ge w_i^n + a \ge u_i^n$ and $f$ is convex, this last inequality is satisfied, inferring the desired conclusion.

* If a “fast shock” is present, connecting $(u_i^n, w_i^n)$ to $(u_{i-1}^n, u_{i-1}^n - a)$ with slope given by Rankine-Hugoniot condition, we get

$$ 
\begin{aligned}
& \int_{x_{i-1/2}}^{x_i} (\tilde{u}^n(x, t^{n+1})^2 - \tilde{u}^n(x, t^n)^2) dx \\
&= \Delta t \frac{f(u_{i-1}^n) - f(u_i^n)}{u_{i-1}^n - u_i^n + (u_{i-1}^n - a) - w_i^n} ((u_{i-1}^n)^2 - (u_i^n)^2) \\
&= \Delta t \frac{f(u_{i-1}^n) - f(u_i^n)}{u_{i-1}^n - u_i^n + (u_{i-1}^n - a) - w_i^n} ((u_{i-1}^n)^2 - (w_i^n + a)^2 + (w_i^n + a)^2 - (u_i^n)^2) \\
&= \Delta t (f(u_{i-1}^n) - f(u_i^n)) (u_{i-1}^n + w_i^n + a) \frac{u_{i-1}^n - (w_i^n + a)}{u_{i-1}^n - u_i^n + (u_{i-1}^n - a) - w_i^n} \\
&\quad + \Delta t (f(u_{i-1}^n) - f(u_i^n)) (w_i^n + a + u_i^n) \frac{(w_i^n + a) - w_i^n}{u_{i-1}^n - u_i^n + (u_{i-1}^n - a) - w_i^n} \\
&= \Delta t (f(u_{i-1}^n) - f(u_i^n)) (u_{i-1}^n + w_i^n + a) \frac{I_l}{I_r + 2I_l} \\
&\quad + \Delta t (f(u_{i-1}^n) - f(u_i^n)) (w_i^n + a + u_i^n) \frac{I_r}{I_r + 2I_l}
\end{aligned}
\quad (3.55) $$

and

$$ 
\begin{aligned}
& \int_{x_{i-1/2}}^{x_i} (\tilde{w}^n(x, t^{n+1})^2 - \tilde{w}^n(x, t^n)^2) dx \\
&= \Delta t \frac{f(u_{i-1}^n) - f(u_i^n)}{u_{i-1}^n - u_i^n + (u_{i-1}^n - a) - w_i^n} ((u_{i-1}^n - a)^2 - (w_i^n)^2) \\
&= \Delta t (f(u_{i-1}^n) - f(u_i^n)) (u_{i-1}^n - a + w_i^n) \frac{(u_{i-1}^n - a) - w_i^n}{u_{i-1}^n - u_i^n + (u_{i-1}^n - a) - w_i^n} \\
&= \Delta t (f(u_{i-1}^n) - f(u_i^n)) (u_{i-1}^n + w_i^n + a) \frac{I_l}{I_r + 2I_l} \\
&\quad - a \Delta t (f(u_{i-1}^n) - f(u_i^n)) \frac{(u_{i-1}^n - a) - w_i^n}{u_{i-1}^n - u_i^n + (u_{i-1}^n - a) - w_i^n}
\end{aligned}
\quad (3.56) $$

and

$$ 
\begin{aligned}
a |\partial_t \tilde{w}^n| (K_i^{n,l}) &= a \Delta t \frac{f(u_{i-1}^n) - f(u_i^n)}{u_{i-1}^n - u_i^n + (u_{i-1}^n - a) - w_i^n} |u_{i-1}^n - a - w_i^n| \\
&= a \Delta t (f(u_{i-1}^n) - f(u_i^n)) \frac{(u_{i-1}^n - a) - w_i^n}{u_{i-1}^n - u_i^n + (u_{i-1}^n - a) - w_i^n}.
\end{aligned}
\quad (3.57) $$

Above, we use the same notation $I_r = (w_r + a) - u_r$, $I_l = u_l - (w_r + a)$ introduced in Subcase 2.3.2. We also define

$$ c_1 = \frac{I_r}{I_r + 2I_l} \quad c_2 = \frac{2I_l}{I_r + 2I_l} $$
and notice that $c_1 + c_2 = 1$.

By(3.55), (3.56) and (3.57), we get

$$
\begin{align}
& -\frac{1}{2} \int_{x_{i-1/2}}^{x_i} (\tilde{u}^n(x, t^{n+1})^2 - \tilde{u}^n(x, t^n)^2) \, dx \nonumber \\
& -\frac{1}{2} \int_{x_{i-1/2}}^{x_i} (\tilde{w}^n(x, t^{n+1})^2 - \tilde{w}^n(x, t^n)^2) \, dx - a |\partial_t \tilde{w}^n| (K_i^{n,l}) \nonumber \\
&= -\frac{1}{2} \Delta t (f(u_{i-1}^n) - f(u_i^n)) (u_{i-1}^n + w_i^n + a) c_2 \nonumber \\
& -\frac{1}{2} \Delta t (f(u_{i-1}^n) - f(u_i^n)) (w_i^n + a + u_i^n) c_1 \nonumber \\
&= -\frac{1}{2} \Delta t (f(u_{i-1}^n) - f(w_i^n + a)) (u_{i-1}^n + w_i^n + a) c_2 \nonumber \\
& -\frac{1}{2} \Delta t (f(w_i^n + a) - f(u_i^n)) (u_{i-1}^n + w_i^n + a) c_2 \nonumber \\
& -\frac{1}{2} \Delta t (f(u_{i-1}^n) - f(w_i^n + a)) (w_i^n + a + u_i^n) c_1 \nonumber \\
& -\frac{1}{2} \Delta t (f(w_i^n + a) - f(u_i^n)) (w_i^n + a + u_i^n) c_1. \tag{3.58}
\end{align}
$$

Besides, we can develop

$$
\begin{align}
& -\Delta t [G(u_i^n) - G(u_{i-1}^n)] \nonumber \\
&= -\Delta t \left[ f(u_i^n)u_i^n - f(u_{i-1}^n)u_{i-1}^n + \int_{u_i^n}^{u_{i-1}^n} f(\xi) \, d\xi \right] \nonumber \\
&= -\Delta t \left[ c_1 f(u_i^n)u_i^n + c_2 f(u_i^n)u_i^n - c_1 f(u_{i-1}^n)u_{i-1}^n - c_2 f(u_{i-1}^n)u_{i-1}^n + \int_{u_i^n}^{u_{i-1}^n} f(\xi) \, d\xi \right] \tag{3.59}
\end{align}
$$

Now, summing the terms involving $c_2$ in the expressions above, we get

$$
\begin{align*}
& -\frac{c_2}{2} (f(u_{i-1}^n) - f(w_i^n + a)) (u_{i-1}^n + w_i^n + a) \\
& -\frac{c_2}{2} (f(w_i^n + a) - f(u_i^n)) (u_{i-1}^n + w_i^n + a) \\
& -c_2 f(u_i^n) u_i^n + c_2 f(u_{i-1}^n) u_{i-1}^n \\
&= \frac{c_2}{2} (f(u_{i-1}^n) + f(w_i^n + a)) (u_{i-1}^n - (w_i^n + a)) \\
& + \frac{c_2}{2} (f(u_i^n) + f(w_i^n + a)) ((w_i^n + a) - u_i^n) \\
& -\frac{c_2}{2} (f(w_i^n + a) - f(u_i^n)) (u_{i-1}^n - u_i^n)
\end{align*}
$$
and, summing the terms involving $c_1$,

$$ 
\begin{aligned}
& -\frac{c_1}{2} (f(u_{i-1}^n) - f(w_i^n + a)) (w_i^n + a + u_i^n) \\
& -\frac{c_1}{2} (f(w_i^n + a) - f(u_i^n)) (w_i^n + a + u_i^n) \\
& -c_1 f(u_i^n) u_i^n + c_1 f(u_{i-1}^n) u_{i-1}^n \\
& = \frac{c_1}{2} (f(u_{i-1}^n) + f(w_i^n + a)) (u_{i-1}^n - (w_i^n + a)) \\
& + \frac{c_1}{2} (f(u_i^n) + f(w_i^n + a)) ((w_i^n + a) - u_i^n) \\
& + \frac{c_1}{2} (f(u_{i-1}^n) - f(w_i^n + a)) (u_{i-1}^n - u_i^n).
\end{aligned}
$$

Therefore, summing (3.58) with (3.59) and recalling that $c_1 + c_2 = 1$, we obtain $\mathcal{I}_1 + \mathcal{I}_2 + \mathcal{I}_3$ with

$$ \mathcal{I}_1 := \Delta t \left[ \frac{1}{2} (f(u_{i-1}^n) + f(w_i^n + a)) (u_{i-1}^n - (w_i^n + a)) - \int_{w_i^n + a}^{u_{i-1}^n} f(\xi) \, d\xi \right] \geq 0 $$

and

$$ \mathcal{I}_2 := \Delta t \left[ \frac{1}{2} (f(u_i^n) + f(w_i^n + a)) ((w_i^n + a) - u_i^n) - \int_{u_i^n}^{w_i^n + a} f(\xi) \, d\xi \right] \geq 0 $$

by convexity of $f$. The third term is

$$ 
\begin{aligned}
\mathcal{I}_3 &:= -\frac{c_2}{2} (f(w_i^n + a) - f(u_i^n)) (u_{i-1}^n - u_i^n) + \frac{c_1}{2} (f(u_{i-1}^n) - f(w_i^n + a)) (u_{i-1}^n - u_i^n) \\
&= \frac{I_r I_l}{I_r + 2I_l} (u_{i-1}^n - u_i^n) (\mu_l - \mu_r) \geq 0
\end{aligned}
$$

since $u_{i-1}^n > u_i^n$ and $\mu_l > \mu_r$ as required by the entropy condition.

This concludes the proof of (3.49). □

The previous Lemma allows to prove the following result.

**Proposition 3.11.** Let $u_0, w_0 \in \mathbf{L}^1(\mathbb{R}) \cap \text{BV}(\mathbb{R})$, such that $|u_0(x) - w_0(x)| \leq a$, and $\Delta x, \Delta t > 0$ satisfying the CFL condition (3.1). Let $(u^\Delta, w^\Delta)$ be the approximate solutions defined by (3.12), (3.39), then this couple satisfies the weak hysteresis relationship (1.8) for almost every $t \in [0, T[$.

*Proof.* Consider a cell $K_i^n$ and the exact solutions $\tilde{u}^n(x, t), \tilde{w}^n(x, t)$ defined on $K_i^n$ with data at time $t^n$ equal to $u_i^n, w_i^n$. We know that, for each $x \in ]x_{i-1/2}, x_{i+1/2}[$, the relationship $\tilde{w}^n(x, t) = \mathcal{F}[\tilde{u}(x, \cdot), w_i^n](t)$ holds, so, by integrating it on $K_i$, we get, by Proposition 1.1

$$ \int_{x_{i-1/2}}^{x_{i+1/2}} \int_{t^n}^{t^{n+1}} (\tilde{u}^n(x, t+) - \tilde{w}^n(x, t+)) d(\partial_t \tilde{w}^n(x, \cdot)) dx \geq a |\partial_t \tilde{w}^n|(K_i^n). \quad (3.60) $$

As usual, we split $K_i^n$ into the two subcells $K_i^{n,l}$ and $K_i^{n,r}$. If in $K_i^{n,l}$ there are no jump discontinuities, then $\tilde{u}^n$ and $\tilde{w}^n$ are Lipschitz continuous and they satisfy the PDE strongly
almost everywhere in $K_i^{n,l}$. In particular, $\partial_t \tilde{w}^n(x, \cdot)$ is not only a measure but a function for a. e. $x$, so using the strong formulation of the PDE, the left hand side of (3.60) restricted to $K_i^{n,l}$ writes as follow

$$
\begin{align}
& \int_{x_{i-1/2}}^{x_i} \int_{t^n}^{t^{n+1}} (\tilde{u}^n(x, t+) - \tilde{w}^n(x, t+)) d(\partial_t \tilde{w}^n(x, \cdot)) dx \nonumber \\
&= \int_{x_{i-1/2}}^{x_i} \int_{t^n}^{t^{n+1}} (\tilde{u}^n(x, t) - \tilde{w}^n(x, t)) \partial_t \tilde{w}^n dt dx \nonumber \\
&= \int_{x_{i-1/2}}^{x_i} \int_{t^n}^{t^{n+1}} (-\tilde{u}^n \partial_t \tilde{u}^n - \tilde{u}^n \partial_x f(\tilde{u}^n)) dt dx - \frac{1}{2} \int_{x_{i-1/2}}^{x_i} (\tilde{w}^n(x, t^{n+1}-)^2 - \tilde{w}^n(x, t^n+)^2) dx \nonumber \\
&= -\Delta t (G(u_i^n) - G(\tilde{u}^n(x_{i-1/2}+, t))) - \frac{1}{2} \int_{x_{i-1/2}}^{x_i} (\tilde{u}^n(x, t^{n+1}-)^2 - \tilde{u}^n(x, t^n+)^2) dx \nonumber \\
&\quad - \frac{1}{2} \int_{x_{i-1/2}}^{x_i} ((\tilde{w}^n(x, t^{n+1}-)^2 - \tilde{w}^n(x, t^n+)^2) dx, \tag{3.61}
\end{align}
$$

where we recall $G(u) := \int_0^u \xi f'(\xi) d\xi = uf(u) - \int_0^u f(\xi) d\xi$.
Similarly, if there are no shocks in $K_i^{n,r}$,

$$
\begin{align}
& \int_{x_i}^{x_{i+1/2}} \int_{t^n}^{t^{n+1}} (\tilde{u}^n(x, t+) - \tilde{w}^n(x, t+)) d(\partial_t \tilde{w}^n(x, \cdot)) dx \nonumber \\
&= -\Delta t (G(\tilde{u}^n(x_{i+1/2}-, t)) - G(u_i^n)) - \frac{1}{2} \int_{x_i}^{x_{i+1/2}} (\tilde{u}^n(x, t^{n+1}-)^2 - \tilde{u}^n(x, t^n+)^2) dx \nonumber \\
&\quad - \frac{1}{2} \int_{x_i}^{x_{i+1/2}} (\tilde{w}^n(x, t^{n+1}-)^2 - \tilde{w}^n(x, t^n+)^2) dx. \tag{3.62}
\end{align}
$$

Combining (3.60), (3.61) and (3.62) to handle rarefactions, together with the results of Lemma 3.10 to account for shocks, we generally obtain in $K_i^n$ the following inequality

$$
\begin{equation}
\begin{split}
& -\Delta t \left( G(\tilde{u}^n(x_{i+1/2}-, t)) - G(\tilde{u}^n(x_{i-1/2}+, t)) \right) - \frac{1}{2} \int_{x_{i-1/2}}^{x_{i+1/2}} (\tilde{u}^n(x, t^{n+1}-)^2 - \tilde{u}^n(x, t^n+)^2) dx \\
& \quad - \frac{1}{2} \int_{x_{i-1/2}}^{x_{i+1/2}} (\tilde{w}^n(x, t^{n+1}-)^2 - \tilde{w}^n(x, t^n+)^2) dx \geq a |\partial_t \tilde{w}^n|(K_i^n)
\end{split}
\tag{3.63}
\end{equation}
$$

As usual, recalling that $\tilde{u}^n(x, t^n) \equiv u_i^n$, $\tilde{w}^n(x, t^n) \equiv w_i^n$ and by Jensen's inequality, we can infer that

$$
\begin{equation}
\begin{split}
& -\Delta t \left( G(\tilde{u}^n(x_{i+1/2}-, t)) - G(\tilde{u}^n(x_{i-1/2}+, t)) \right) \\
& \quad - \frac{1}{2} \int_{x_{i-1/2}}^{x_{i+1/2}} \left( (u_i^{n+1})^2 - (u_i^n)^2 + (w_i^{n+1})^2 - (w_i^n)^2 \right) dx \geq a |\partial_t \tilde{w}^n|(K_i^n).
\end{split}
\tag{3.64}
\end{equation}
$$

Moreover, notice that $\tilde{w}_i^n(x, \cdot)$ is monotone in $t$ for $x \in K_i$ and it may only change monotonicity across $\{x_i\} \times [t^n, t^{n+1}[$, which has measure 0 with respect to $\partial_t \tilde{w}^n$. So, by
using Jensen's inequality once again, it holds that

$$
\begin{align*}
a |\partial_t \tilde{w}^n|(K_i) &= a \int_{x_{i-1/2}}^{x_{i+1/2}} |\tilde{w}^n(x, t^{n+1}) - \tilde{w}^n(x, t^n)| dx \\
&= a \int_{x_{i-1/2}}^{x_{i+1/2}} |\tilde{w}^n(x, t^{n+1}) - w_i^n| dx \\
&\geq a \int_{x_{i-1/2}}^{x_{i+1/2}} |w_i^{n+1} - w_i^n| dx.
\end{align*}
$$

Using this last inequality and summing (3.64) over $i \in \mathbb{Z}$ and recalling (3.39), we finally get

$$
\begin{equation}
\begin{split}
& -\Delta t \sum_{i \in \mathbb{Z}} (G(\tilde{u}^n(x_{i+1/2}-, t)) - G(\tilde{u}^n(x_{i-1/2}+, t))) \\
& \quad - \frac{1}{2} \int_{\mathbb{R}} (u^\Delta(x, t^{n+1})^2 - u^\Delta(x, t^n)^2) dx \\
& \quad - \frac{1}{2} \int_{\mathbb{R}} (w^\Delta(x, t^{n+1})^2 - w^\Delta(x, t^n)^2) dx \\
& \geq a \int_{\mathbb{R}} |w^\Delta(x, t^{n+1}) - w^\Delta(x, t^n)| dx.
\end{split}
\tag{3.65}
\end{equation}
$$

Now, if $\tilde{u}^n(x_{i+1/2}-, t) = \tilde{u}^n(x_{i+1/2}+, t)$ for some $i \in \mathbb{Z}$, then they cancel out in the sum. If instead $\tilde{u}^n(x_{i+1/2}-, t) \neq \tilde{u}^n(x_{i+1/2}+, t)$, then there is a stationary shock in $u$ at $x_{i+1/2}$, i.e. $\tilde{u}^n(x_{i+1/2}-, t) = u_i^n$ and $\tilde{u}^n(x_{i+1/2}+, t) = u_{i+1}^n$ with $f(u_i^n) = f(u_{i+1}^n)$ and $u_i^n > u_{i+1}^n$. (see Subcase 2.3.1). In such case

$$
\begin{align*}
& - (G(\tilde{u}^n(x_{i+1/2}-, t)) - G(\tilde{u}^n(x_{i+1/2}+, t))) \\
&= -u_i^n f(u_i^n) + u_{i+1}^n f(u_{i+1}^n) + \int_{u_{i+1}^n}^{u_i^n} f(\xi) \, d\xi \\
&= \int_{u_{i+1}^n}^{u_i^n} f(\xi) \, d\xi - (u_i^n - u_{i+1}^n) f(u_i^n) \le 0
\end{align*}
$$

where the last inequality is due to the fact that $f$ is convex and $u_i^n > u_{i+1}^n$.

Hence, $-\Delta t \sum_{i \in \mathbb{Z}} (G(\tilde{u}^n(x_{i+1/2}-, t)) - G(\tilde{u}^n(x_{i-1/2}+, t))) \le 0$, so from (3.65) we finally get

$$
\begin{align}
& -\frac{1}{2} \int_{\mathbb{R}} (u^{\Delta}(x, t^{n+1})^2 - u^{\Delta}(x, t^n)^2) dx - \frac{1}{2} \int_{\mathbb{R}} (w^{\Delta}(x, t^{n+1})^2 - w^{\Delta}(x, t^n)^2) dx \tag{3.66} \\
& \geq a \int_{\mathbb{R}} |w^{\Delta}(x, t^{n+1}) - w^{\Delta}(x, t^n)| dx. \nonumber
\end{align}
$$

Let now consider $t \in ]0, T[$ such that $t \in [t^N, t^{N+1}[$ for some $N \in \mathbb{N}$. Then, summing (3.66) over $n = 0, \dots, N-1$, we get

$$
\begin{align}
& -\frac{1}{2} \int_{\mathbb{R}} (u^{\Delta}(x,t)^2 - u_0(x)^2 + w^{\Delta}(x,t)^2 - w_0(x)^2) dx \tag{3.67} \\
& \geq a \sum_{n=0}^{N-1} \int_{\mathbb{R}} |w^{\Delta}(x,t^{n+1}) - w^{\Delta}(x,t^n)| dx = a |\partial_t w^{\Delta}| (\mathbb{R} \times ]0,t[). \nonumber
\end{align}
$$

□
## 3.5 Convergence of approximate solutions: the existence theorem

**Theorem 3.12.** Let $T > 0$, $u_0, w_0 \in BV(\mathbb{R})$, $|u_0(x) - w_0(x)| \le a$, and consider a sequence of meshes $T_m$ and time steps $\Delta t_m$ satisfying the CFL condition (3.1) for every $m \in \mathbb{N}$ and such that $\Delta x_m, \Delta t_m \to 0$ as $m \to +\infty$. Then the sequence of approximate solutions $\{u_m\}_{m \in \mathbb{N}} := \{u^{\Delta m}\}_{m \in \mathbb{N}}$ and $\{w_m\}_{m \in \mathbb{N}} := \{w^{\Delta m}\}_{m \in \mathbb{N}}$ defined by the scheme (3.12) and by (3.39) converges in $\mathbf{L}_{\text{loc}}^1(\mathbb{R} \times [0, T])$ to $u$ and $w$ respectively, where $(u, w)$ is an entropy weak solution to (1.1) according to Definition 1.2.

*Proof.* It can be shown that $u_m, w_m$ converge in $\mathbf{L}_{\text{loc}}^1(\mathbb{R} \times [0, T])$ to a couple $u, w \in C^0([0, T; \mathbf{L}_{\text{loc}}^1(\mathbb{R})])$ and that this couple satisfies (1.6), by following classical arguments; which can be found in Appendix A. To conclude, we then have to show inequalities (1.7) and (1.8).

Condition (1.7) is true since it holds for $u_m, w_m$, see (3.11), and the $\mathbf{L}_{\text{loc}}^1(\mathbb{R} \times [0, T])$ convergence implies (up to a subsequence) almost everywhere pointwise convergence.

Finally, also the weak hysteresis relationship (1.8), which holds for each approximate solution (see Proposition 3.11), passes to the limit, for almost every $t \in [0, T]$, due to the $\mathbf{L}_{\text{loc}}^1(\mathbb{R} \times [0, T])$ convergence and to the lower semicontinuity of the mass of measures (see e.g. [3]). Moreover, since $u_0, w_0 \in \mathbf{L}^1(\mathbb{R}) \cap BV(\mathbb{R})$, it can be checked that $u, w \in \mathbf{L}^1(\mathbb{R} \times [0, T])$, and not only locally. In addition, as $u, w \in \mathbf{L}^{\infty}(\mathbb{R} \times [0, T])$, we also have $u, w \in \mathbf{L}^2(\mathbb{R} \times [0, T])$, which is indeed required by (1.8). □

**Remark 3.13. (The linear case)** The work [4] considers problem (1.1) with linear flux $f(u) = u$. The main result of such paper was the existence and uniqueness of a solution for the Cauchy problem with BV initial data constructed via the wave-front tracking method. However, no investigation about numerical schemes was carried out.

Reasoning as we did in this paper, it can be shown that in a cell $K_i^n$, the solution $\tilde{u}^n$ of the linear equation satisfies the conservation law

$$ \partial_t u + \partial_x \bar{f}_{w_i^n}(u) = 0, \qquad (3.68) $$

where

$$ \bar{f}_{w_i^n} = \begin{cases} \frac{1}{2}u + \frac{1}{2}(w_i^n - a) & u \in (-\infty, w_i^n - a], \\ u & u \in [w_i^n - a, w_i^n + a], \\ \frac{1}{2}u + \frac{1}{2}(w_i^n + a) & u \in [w_i^n + a, +\infty). \end{cases} \qquad (3.69) $$

Then, since waves have positive speed, the solution to the Riemann problem at $x = 0$ coincides with the left datum $u_l$, so

$$ \left\{ \begin{array}{ll} u_i^{n+1} = u_i^n - \frac{\Delta t}{\Delta x} (f_{w_i^n}(u_i^n) - f_{w_i^n}(u_{i-1}^n)), & \forall i \in \mathbb{Z}, \\ w_i^{n+1} = u_i^n + w_i^n - \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) - u_i^{n+1}, & \forall i \in \mathbb{Z}, \\ u_i^0 = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} u_0(x) dx, \quad w_i^0 = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} w_0(x) dx, & \forall i \in \mathbb{Z}, \\ |u_0(x) - w_0(x)| \le a, & \forall x \in \mathbb{R}. \end{array} \right. \qquad (3.70) $$

The numerical solution associated to scheme (3.70) for the linear flux has the same properties of the more general one considered in this paper. In particular, the corresponding
<figure><img src="image_16.png" /><figcaption>Figure 13: Left: double rarefaction wave for $u$ propagating to the right, see Subcase 2.2.1. Middle: double rarefaction wave propagating to the left, see Subcase 2.2.2. Right: centered rarefaction, see 2.2.3.</figcaption></figure>

approximate solution converges to the theoretical solution constructed in [4] via a wave-front tracking algorithm, as the size of the mesh goes to 0, due to the uniqueness of the entropy weak solution.

## 3.6 Some numerical examples

We present some examples of numerical solutions generated via the scheme (3.12) to the Riemann problems (1.1), (2.1). In particular, Figure 13 shows the three possible configurations of rarefaction wave solutions described in Section 2.2, while Figure 14 shows cases involving shock waves, see Section 2.3. In all cases, we consider the hysteresis parameter $a = 1$, the flux function $f(u) = \frac{1}{2}u^2$, whereas the other data and parameters are summarized in Table 1 and $\Delta t$ is such that the CFL condition (3.1) holds as an equality. The analytical solutions of some of these problems are illustrated in Figures 4–9, Section 2. In all cases, the numerical scheme captures well the exact profiles for both $u$ and $w$.

<table><caption>Table 1: Parameters used for the simulations.</caption><thead><tr><th>Figure</th><th>(u<sub>l</sub>, w<sub>l</sub>)</th><th>(u<sub>r</sub>, w<sub>r</sub>)</th><th>Δx</th><th>t</th></tr></thead><tbody><tr><td>Figure 13 (left)</td><td>(1,0.5)</td><td>(3,3)</td><td>10<sup>-3</sup></td><td>0.25</td></tr><tr><td>Figure 13 (middle)</td><td>(-3,-3)</td><td>(-1.5,-1)</td><td>10<sup>-3</sup></td><td>0.25</td></tr><tr><td>Figure 13 (right)</td><td>(-2,-1.5)</td><td>(1,1.5)</td><td>10<sup>-3</sup></td><td>0.4</td></tr><tr><td>Figure 14 (left)</td><td>(1.5,2)</td><td>(0.5,0)</td><td>10<sup>-3</sup></td><td>0.5</td></tr><tr><td>Figure 14 (middle)</td><td>(-0.5,0)</td><td>(-1.5,-2)</td><td>10<sup>-3</sup></td><td>0.5</td></tr><tr><td>Figure 14 (right)</td><td>(1.5,2)</td><td>(-1,-1)</td><td>10<sup>-3</sup></td><td>0.5</td></tr></tbody></table>
<figure><img src="image_17.png" /><figcaption>Figure 14: Left: double shock wave for $u$ propagating to the right with $\mu_r \ge \mu_l$, see Subcase 2.3.2. Middle: double shock wave propagating to the left with $\nu_r \ge \nu_l$, see Subcase 2.3.3. Right: “fast shock” with $\mu_r < \mu_l$, see again Subcase 2.3.2.</figcaption></figure>

We also compute the numerical solution $(u^\Delta, w^\Delta)$ of the Cauchy problem with a more general initial datum, namely $u_0 = e^{-x^2/2}$, with $u_0 = w_0$, the flux function $f(u) = \frac{1}{2}u^2$, $a = 1$, $\Delta x = 10^{-3}$ and $\Delta t$ such that the CFL condition (3.1) holds as an equality. Figure 15 shows the evolution of $u^\Delta(\cdot, t)$ and $w^\Delta(\cdot, t)$ at increasing time steps $t$.

In the same plots, $u^\Delta$ is compared with the solutions to the two equations without hysteresis

$$ \partial_t u + \partial_x f(u) = 0 \quad \text{and} \quad \partial_t u + \frac{1}{2} \partial_x f(u) = 0 $$

which are depicted by the dashed red line and by the dot-dashed blue line, respectively. To compute such solutions we still used the scheme (3.12), with flux either $f(u) = u^2/2$ or $f(u) = u^2/4$, where, by setting the hysteresis parameter $a$ large enough and by also setting $u_0 = w_0$, we ensure $\partial_t w = 0$ for all $(x, t)$. In this way, by considering only the unknown $u$, we are approximating the solutions of the two equations without hysteresis.

Referring again to Figure 15, we observe that the speed with which $u^\Delta$ propagates to the right is intermediate between those of the solutions of the two equations without hysteresis. This behaviour is expected as the equation with hysteresis is a combination of the non-hysteretic two, recall Remark 2.5. Moreover, the presence of hysteresis lowers the peak of $u^\Delta$, as large variations in $u(x, \cdot)$ lead to $(u(x, \cdot), w(x, \cdot))$ following the hysteresis boundary, where $\partial_t w \ne 0$. From the relationship $\partial_t u + \partial_x f(u) = -\partial_t w$, we see that the term $-\partial_t w$ acts as a source term damping large variations of $u$, since when $(u, w)$ is on the boundary of hysteresis regions $\text{sign}(\partial_t u) = -\text{sign}(-\partial_t w)$.

In Table 2, we also highlight the sum of the $L^2$ norms squared in space of $u^\Delta$ and $w^\Delta$, which decreases in time accordingly with (1.8).
<figure><img src="image_18.png" /><figcaption>Figure 15: Graphs of $u^\Delta(\cdot, t)$ (top row) and $u^\Delta(\cdot, t)$ (bottom row) computed at increasing times. The top plots include also the solutions to the equations without hysteresis $\partial_t u + \partial_x f(u) = 0$ (dashed red) and $\partial_t u + \frac{1}{2} \partial_x f(u) = 0$ (dot-dashed blue).</figcaption></figure>

<table><caption>Table 2: The $L^2(\mathbb{R})$ norm square of $u^\Delta(\cdot, t)$, $w^\Delta(\cdot, t)$ and their decreasing in time sum.</caption><thead><tr><th>L<sup>2</sup> norm square</th><th>t = 0</th><th>t = 0.2</th><th>t = 0.4</th><th>t = 0.6</th></tr></thead><tbody><tr><td>||u<sup>Δ</sup>(·, t)||<sup>2</sup><sub>L<sup>2</sup>(R)</sub></td><td>22.1557</td><td>19.2978</td><td>17.0486</td><td>14.9098</td></tr><tr><td>||w<sup>Δ</sup>(·, t)||<sup>2</sup><sub>L<sup>2</sup>(R)</sub></td><td>22.1557</td><td>22.5789</td><td>21.8581</td><td>21.0941</td></tr><tr><td>||u<sup>Δ</sup>(·, t)||<sup>2</sup><sub>L<sup>2</sup>(R)</sub> + ||w<sup>Δ</sup>(·, t)||<sup>2</sup><sub>L<sup>2</sup>(R)</sub></td><td>44.3114</td><td>41.8767</td><td>38.9067</td><td>36.0039</td></tr></tbody></table>

# 4 Stability

Now we show how the condition (1.6) in the definition of entropy weak solution characterizes the entropy-allowed shock discontinuities.

**Proposition 4.1.** Suppose a couple $(u, w) \in C^0([0, T[, L^1_{\text{loc}}(\mathbb{R})]$ to have a shock discontinuity on a curve $(\sigma(t), t)$ between two constant states $(u_-, w_-) \neq (u_+, w_+)$ and to be an entropy weak solution away from that curve. Then this couple satisfies (1.6) and is an entropy weak solution on the whole domain, i.e. this shock is entropy admissible, if and only if $\sigma'$ satisfies the Rankine-Hugoniot condition (2.2), $u_- \ge u_+$ and one of the following holds

$$i) w_{-} = w_{+};$$

$$ii) f(u_{-}) = f(u_{+}) \text{ and } \sigma' = 0;$$

$$iii) w_{-} > w_{+}, f(u_{-}) > f(u_{+}), w_{-} = u_{-} - a \text{ and } w_{+} = u_{+} - a;$$

$$iv) w_{-} > w_{+}, f(u_{-}) > f(u_{+}), w_{-} = u_{-} - a \text{ and } \mu_{+} \le \mu_{-} \text{ where}$$

$$\mu_{+} = \frac{f(w_{+} + a) - f(u_{+})}{(w_{+} + a) - u_{+}}, \quad \text{and} \quad \mu_{-} = \frac{1}{2} \frac{f(u_{-}) - f(w_{+} + a)}{u_{-} - (w_{+} + a)};$$
$$v) w_{-} > w_{+}, f(u_{-}) < f(u_{+}), w_{+} = u_{+} + a \text{ and } w_{-} = u_{-} + a;$$

$$vi) w_{-} > w_{+}, f(u_{-}) < f(u_{+}), w_{+} = u_{+} + a \text{ and } \nu_{+} \le \nu_{-} \text{ where}$$

$$\nu_{-} = \frac{f(u_{-}) - f(w_{-} - a)}{u_{-} - (w_{-} - a)}, \quad \text{and} \quad \nu_{+} = \frac{1}{2} \frac{f(w_{-} - a) - f(u_{+})}{w_{-} - a - u_{+}}.$$

**Remark 4.2.** Before proving the above proposition, let us underline that, when solving the Riemann problem in Case 2.3, Section 2, we only used entropy admissible shocks, see Remarks 2.1 and 2.2.

*Proof.* Integrating by parts (1.6) and following the same argument as in [6, Theorem 4.3] we can rewrite the entropy condition along the discontinuity as

$$\sigma' \left[ \left( |u_+ - k| + |w_+ - \hat{k}| \right) - \left( |u_- - k| + |w_- - \hat{k}| \right) \right] - [\operatorname{sign}(u_+ - k)(f(u_+) - f(k)) - \operatorname{sign}(u_- - k)(f(u_-) - f(k))] \ge 0. \quad (4.1)$$

So the couple $(u, w)$ is an entropy solution if and only if (4.1) holds for every $(k, \hat{k}) \in \mathcal{L}$. In particular, the Rankine-Hugoniot condition (2.2) can be deduced as a necessary condition for weak solutions.

First of all, let us notice that, if $w_{-} = w_{+}$, then (4.1) reduces to the classical entropy condition on shocks and also the Rankine-Hugoniot condition (2.2) becomes the classical one. In such case, it is well known that the necessary and sufficient condition for entropic shocks is $\sigma'$ to satisfy the Rankine-Hugoniot condition with $u_{-} > u_{+}$; this deals with case i).

From now on, we suppose $w_{-} \ne w_{+}$. We distinguish the following cases according to the relationship between $f(u_{-})$, $f(u_{+})$, $u_{-}$ and $u_{+}$.

**Case f(u_{-}) = f(u_{+})**: From (2.2), either $\sigma' = 0$ or $u_{-} - u_{+} + w_{-} - w_{+} = 0$. But if the latter is true, it is sufficient to first choose in (4.1), $(k, \hat{k}) \in \mathcal{L}$ with $k \ge \max(u_{-}, u_{+})$ and $\hat{k} \le \min(w_{-}, w_{+})$, then choose instead $(k, \hat{k}) \in \mathcal{L}$ with $k \le \min(u_{-}, u_{+})$ and $\hat{k} \ge \max(w_{-}, w_{+})$, to infer that

$$\sigma'[u_{-} - u_{+} + w_{+} - w_{-}] = 0,$$

implying $\sigma' = 0$. Notice that the choices of $(k, \hat{k}) \in \mathcal{L}$ are possible as $(u_{-}, w_{-})$, $(u_{+}, w_{+}) \in \mathcal{L}$ with $u_{-} - u_{+} + w_{-} - w_{+} = 0$, moreover $u_{-} - u_{+} + w_{+} - w_{-} \ne 0$ since $u_{-} - u_{+} + w_{-} - w_{+} = 0$ and $(u_{-}, w_{-}) \ne (u_{+}, w_{+})$. So, in any case when $f(u_{-}) = f(u_{+})$, necessarily $\sigma' = 0$, thus (4.1) becomes

$$-\operatorname{sign}(u_{+} - k)(f(u_{+}) - f(k)) + \operatorname{sign}(u_{-} - k)(f(u_{-}) - f(k)) \ge 0. \quad (4.2)$$

Now it is easy to show, thanks to convexity of $f$ and since $f(u_{-}) = f(u_{+})$, that (4.2) is true if and only if $u_{-} \ge u_{+}$. So we just proved that under the Rankine-Hugoniot condition, when $f(u_{-}) = f(u_{+})$, $(u, w)$ is an entropy solution if and only if $\sigma' = 0$ and $u_{-} \ge u_{+}$; thus if and only if $u_{-} \ge u_{+}$ and case ii) holds.

Before proceeding with the other cases, we notice the latter case includes also the case when $u_{-} = u_{+}$ and does not pose any restriction on $w_{-}$ and $w_{+}$. Moreover, we can now
also assume $u_{-} - u_{+} + w_{-} - w_{+} \neq 0$, as the opposite would imply $\sigma' = 0$, so by the Rankine-Hugoniot condition (2.2) $f(u_{-}) = f(u_{+})$. This also allows us, by using again (2.2), to write $\sigma'$ as a function of $(u_{-}, w_{-})$ and $(u_{+}, w_{+})$ and substitute it in (4.1). Thus, we get the equivalent condition

$$ \frac{h(k, \hat{k})}{u_{-} - u_{+} + w_{-} - w_{+}} \geq 0, \qquad (4.3) $$

where

$$ h(k, \hat{k}) := (f(u_{-}) - f(u_{+})) \left[ (|u_{+} - k| + |w_{+} - \hat{k}|) - (|u_{-} - k| + |w_{-} - \hat{k}|) \right] \\ - (u_{-} - u_{+} + w_{-} - w_{+}) [\operatorname{sign}(u_{+} - k)(f(u_{+}) - f(k)) \\ - \operatorname{sign}(u_{-} - k)(f(u_{-}) - f(k))] \qquad (4.4) $$

that must hold for every $(k, \hat{k}) \in \mathcal{L}$. (For sake of precision, $h = h(k, \hat{k}; u_{-}, w_{-}, u_{+}, w_{+})$, but for simplicity we will omit the explicit dependence on $(u_{-}, w_{-})$ and $(u_{+}, w_{+})$ in the notation.) We also notice that

$$ h(u_{-}, w_{-}) = (f(u_{-}) - f(u_{+})) (w_{+} - w_{-}) (\operatorname{sign}(w_{+} - w_{-}) - \operatorname{sign}(u_{+} - u_{-})), $$

and

$$ h(u_{+}, w_{+}) = -(f(u_{-}) - f(u_{+})) (w_{+} - w_{-}) (\operatorname{sign}(w_{+} - w_{-}) - \operatorname{sign}(u_{+} - u_{-})), $$

so, if $\operatorname{sign}(w_{+} - w_{-}) \neq \operatorname{sign}(u_{+} - u_{-})$, then we would have $h(u_{-}, w_{-})$ with opposite sign with respect to $h(u_{+}, w_{+})$. As $h$ is continuous with respect to $(k, \hat{k})$ and $(u_{-}, w_{-})$, $(u_{+}, w_{+}) \in \mathcal{L}$, then we could find multiple couples $(k, \hat{k}) \in \mathcal{L}$ either in neighbourhood of $(u_{-}, w_{-})$ or $(u_{+}, w_{+})$ so that (4.3) does not hold. Consequently (4.3) implies that either $u_{-} > u_{+}$ and $w_{+} > w_{-}$ or $u_{-} < u_{+}$ and $w_{-} < w_{+}$. This allows us to reduce the study to only the following remaining 4 cases:

**Case $u_{-} > u_{+}$ and $f(u_{-}) > f(u_{+})$**: by the previous observation $u_{-} - u_{+} + w_{-} - w_{+} > 0$, hence (4.3) holds if and only if $h(k, \hat{k}) \geq 0$ for every $(k, \hat{k}) \in \mathcal{L}$. To understand the sign of $h$, we split the $(k, \hat{k})$-plane in 9 regions, representing the cases when $k \leq u_{+}$, $u_{+} \leq k \leq u_{-}$ or $u_{-} \leq k$ combined with the cases $\hat{k} \leq w_{+}$, $w_{+} \leq \hat{k} \leq w_{-}$ or $w_{-} \leq \hat{k}$, see Figure 16. We also denote the interior of each of these regions by $R_i$, numbered as shown in the figure.

It can be checked that

$$ h(k, \hat{k}) = \begin{cases} 2(f(u_{-}) - f(u_{+}))(w_{-} - w_{+}) & \text{in } R_{1}, \\ 2[(f(u_{-}) - f(u_{+}))(k - u_{-}) + (f(u_{-}) - f(k))(u_{-} - u_{+} + w_{-} - w_{+})] & \text{in } R_{2}, \\ 0 & \text{in } R_{3}, \\ 2(f(u_{-}) - f(u_{+}))(\hat{k} - w_{+}) & \text{in } R_{4}, \\ h_{5}(k, \hat{k}) & \text{in } R_{5}, \\ -2(f(u_{-}) - f(u_{+}))(w_{-} - \hat{k}) & \text{in } R_{6}, \\ 0 & \text{in } R_{7}, \\ -2[(f(u_{-}) - f(u_{+}))(u_{+} - k) + (f(k) - f(u_{+}))(u_{-} - u_{+} + w_{-} - w_{+})] & \text{in } R_{8}, \\ -2(f(u_{-}) - f(u_{+}))(w_{-} - w_{+}) & \text{in } R_{9}. \end{cases} \qquad (4.5) $$
<figure><img src="image_19.png" /><figcaption>Figure 16: Subdivision of $\mathbb{R}^2$ into the regions $R_i$. The points are defined as follows: $A = (\min(u_-, u_+), \max(w_-, w_+))$, $B = (\max(u_-, u_+), \max(w_-, w_+))$, $C = (\min(u_-, u_+), \min(w_-, w_+))$, and $D = (\max(u_-, u_+), \min(w_-, w_+))$.</figcaption></figure>

The explicit computation of $h$ in regions such as $R_2$, $R_4$, $R_5$, $R_6$, and $R_8$ is not really necessary; instead, we can deduce the properties of $h$ in these regions by analysing its derivatives, as we will see later, which are easier to compute.

We see now that $h > 0$ in $R_1$ and $h < 0$ in $R_9$. Moreover, in $R_4$ and $R_6$ it is constant in $k$ and affine in $\hat{k}$, hence necessarily $h > 0$ in $R_4$ and $h < 0$ in $R_6$. In particular, if we want $h \ge 0$ for every $(k, \hat{k}) \in \mathcal{L}$, it follows that $(\max(u_-, u_+), \max(w_-, w_+)) = (u_-, w_+)$ must lie on the lower boundary of $\mathcal{L}$, see Figure 17, left. Therefore, a necessary condition is $w_- = u_- - a$. In $R_2$ instead, $h$ is constant in $\hat{k}$ and it holds

$$ \frac{\partial^2 h}{\partial^2 \hat{k}} = -2(u_{-} - u_{+} + w_{-} - w_{+})f''(k) < 0, $$

so $k \mapsto h(k, \hat{k})$ is concave on $R_2$. Hence, since $h$ is positive on $R_1$ and zero on $R_3$, then we can conclude that $h \ge 0$ on $R_2$. Notice that as $w_- = u_- - a$, then $R_2 \cap \mathcal{L} \neq \emptyset$, see again Figure 17, left, so the positivity of $h$ on $R_2$ is needed in order to not contradict (4.4).

We focus now on $R_5$ where we have

$$ h(k, \hat{k}) = h_5(k, \hat{k}) = (f(u_{-}) - f(u_{+}))[-u_{+} - u_{-} + 2k - w_{+} - w_{-} + 2\hat{k}] \\ -(u_{-} - u_{+} + w_{-} - w_{+})[-f(u_{+}) - f(u_{-}) + 2f(k)], $$

and we study the level set $h(k, \hat{k}) = 0$. Now as before

$$ \frac{\partial^2 h}{\partial^2 \hat{k}} = -(u_{-} - u_{+} + w_{-} - w_{+})2f''(k) < 0, $$

so the function $k \mapsto h(k, \hat{k})$ is concave, and since for $\hat{k} \in (w_+, w_-)$, $h(u_+, \hat{k}) > 0$ and $h(u_-, \hat{k}) < 0$, then for every $\hat{k} \in (w_+, w_-)$ fixed there exists one and only one $k$ such that $(k, \hat{k}) \in R_5$ and $h(k, \hat{k}) = 0$. Let us consider the set $J = \{k \in ]u_-, u_+[, | \exists \hat{k} \text{ s.t. } h(k, \hat{k}) = 0\}$ and denote by $k_1 := \inf J$ and by $k_2 := \sup J$. As

$$ \frac{\partial h}{\partial \hat{k}} = 2(f(u_{-}) - f(u_{+})) > 0, $$
<figure><img src="image_20.png" /><figcaption>Figure 17: Left: focus on the regions $R_2, R_3$ and $R_6$ in the case $u_- > u_+$ with $f(u_-) > f(u_+)$; since $(u_-, w_-) \in \mathcal{L}$ and $h < 0$ in $R_6$, then necessarily $(u_-, w_-)$ must lie on the lower boundary of $\mathcal{L}$, so that $\mathcal{L} \cap R_6 = \emptyset$. Right: focus on $R_5$, where we can see that, to have $\{h < 0\} \cap R_5 \cap \mathcal{L} = \emptyset$, we need to ensure $k_1 \ge (w_+ + a)$. The gray area represents the region $w < u - a$, which is not included in $\mathcal{L}$.</figcaption></figure>

then also for each $k \in ]k_1, k_2[$ there exists exactly one $\hat{k}$ such that $h(k, \hat{k}) = 0$. We then can apply the implicit function theorem at each of such points, and get the existence of $\varphi:]k_1, k_2[ \to ]w_+, w_-[$ of class $C^1$ such that $h(k, \varphi(k)) = 0$ for each $k \in ]k_1, k_2[$. The idea is to apply such theorem locally on each point, then use the uniqueness of the zeros, to extend it globally on $]k_1, k_2[$. Moreover, the implicit function theorem also states that

$$ \varphi'(k) = -\frac{h_k(k, \varphi(k))}{h_{\hat{k}}(k, \varphi(k))} = -\frac{2(f(u_-) - f(u_+)) - 2(u_- - u_+ + w_- - w_+)f'(k)}{2(f(u_-) - f(u_+))}, $$

consequently

$$ \varphi''(k) = \frac{(u_- - u_+ + w_- - w_+)f''(k)}{f(u_-) - f(u_+)} > 0, $$

meaning that $\varphi$ is convex. We can do the same reasoning for the point $(u_-, w_-)$ as $h(u_-, w_-) = 0$ and get that there exists a function $\phi$ defined in a left neighbourhood of $u_-$ such that $h(k, \phi(k)) = 0$ and

$$ \phi'(u_{-}) = -\frac{2((f(u_{-}) - f(u_{+}))-2(u_{-} - u_{+} + w_{-} + w_{+})f'(u_{+}))}{2(f(u_{-}) - f(u_{+}))}. $$

But since $(u_- - u_+ + w_- - w_+) > u_- - u_+$ and $f'(u_+) \le (f(u_-) - f(u_+))/(u_- - u_+)$, by convexity of $f$, it can be shown that $\phi'(u_{-}) > 0$ so actually, by a uniqueness arguments, $\phi$ coincides with $\varphi$ in a left neighbourhood of $u_-$. So we conclude that $k_2 = u_-$, i.e. $\varphi:]k_1, u_-] \to ]w_+, w_-]$, with $\varphi(u_-) = w_-$. Regarding $k_1$ instead, using the notations $I_+ = (w_+ + a) - u_+$, $I_- = u_- - (w_+ + a)$ and

$$ \mu = \frac{f(u_{-}) - f(u_{+})}{u_{-} - u_{+} + w_{-} - w_{+}} = \frac{\mu_{+}I_{+} + 2\mu_{-}I_{-}}{I_{+} + 2I_{-}}, $$

it can be checked that

$$ h(w_{+} + a, w_{+}) = (I_{+} + 2I_{-})[(\mu - \mu_{+})I_{+} - (\mu - \mu_{-})2I_{-}]. $$
<figure><img src="image_21.png" /><figcaption>Figure 18: Sign of the function $h$ when $u_- > u_+$ and $f(u_-) > f(u_+)$. Left: case $\mu_+ \le \mu_-$, where choosing $(u_-, w_+) = (u_-, u_- - a)$ ensures that $h \ge 0$ for all $(k, \hat{k}) \in \mathcal{L}$. Right: case $\mu_+ > \mu_-$ and $w_- \ne u_+-a$, where even after imposing the necessary condition $w_+ = u_+-a$, it does not hold $h \ge 0$ on $\mathcal{L}$. The gray region represents the region $w < u-a$ which is not in $\mathcal{L}$. From these pictures, we also deduce that $w_- = u_- - a$ together with $w_+ = u_+ - a$ would be enough to get an entropy solution.</figcaption></figure>

Since $I_+ \ge 0, I_- > 0$, then $\mu$ is between $\mu_+$ and $\mu_-$ with equality if and only if $\mu_+ = \mu_- = \mu$, and with $\mu = \mu_-$ when $I_+ = 0$. If now $\mu_+ > \mu_-$ and $I_+ \ne 0$, then $\mu_+ > \mu > \mu_-$ hence $h(w_++a, w_+) < 0$. If instead $\mu_+ \le \mu_-$, then $\mu_+ \le \mu \le \mu_-$, which implies that $h(w_++a, w_+) \ge 0$. Also if $I_+ = 0$ i.e. $u_+ = w_+-a$, then $\mu_- = \mu$ so $h(w_++a, w_+) = 0$. From this and the properties of $h(k, w_+)$, we can conclude that $k_1 < w_+ + a$ when $\mu_+ > \mu_- \ne \mu$ and instead $k_1 \ge w_+ + a$ when either $\mu_+ \le \mu_-$ or $I_+ = 0$, that is $u_+ = w_+-a$. By convexity of $\varphi$, $k_1 \ge w_+ + a$ is sufficient and necessary to have $h \ge 0$ in $R_5 \cap \mathcal{L}$, see Figure 17, right. So, in the current case, $\mu_+ \le \mu_-$ or $u_+ = w_+-a$ is required by the positivity of $h$ in $R_5$. The only region left to study is $R_8$, which is constant in $\hat{k}$ and so in that region $h(k, \hat{k}) = h(k, w_+)$, so we can deduce again that conditions $\mu_+ \le \mu_-$ or $u_+ = w_+-a$ are implied. Finally, by looking at Figure 18, where we picture the sign of $h$ in $\mathcal{L}$, we can conclude that in this case we get an entropy solution if and only if either $u_+ = w_+-a$ or $\mu_+ \le \mu_-$, that is if and only if iii) or iv) hold.

**Case $u_- > u_+$ and $f(u_-) < f(u_+)$:** we can treat this case as the previous one with the explicit form (4.5) of $h$ remaining the same. The main difference is in the sign of $h$, since $\text{sign}(f(u_-) - f(u_+))$ is opposite with respect to the previous case, see Figure 19. Under these assumptions, (4.3) implies that $(u_+, w_+)$ must be on the upper boundary of $\mathcal{L}$ and, moreover, in region $R_5$ we get that $h(k, \hat{k}) = 0$, seen as a graph over $k$, is concave. Then conditions $w_- = u_- + a$ or $\nu_- \ge \nu_+$ are deduced, so that $\{h < 0\} \cap \mathcal{L} = \emptyset$, see again Figure 19. Therefore, in this case we get an entropy solution if and only if either $v$) or $vi$) holds.

**Case $u_- < u_+$ and $f(u_-) > f(u_+)$:** now $u_- - u_+ + w_- - w_+ < 0$, so (4.3) implies $h \le 0$ in $\mathcal{L}$. Reasoning as before, we can check that in $R_3$ and $R_7$, $h \equiv 0$, $h < 0$ in $R_1$ and $h > 0$ in $R_9$. Again, $h$ is affine in $R_4$ and $R_6$, respectively negative and positive. As a result, we can still conclude that the point that is now $(u_+, w_+)$ must be on the lower
<figure><img src="image_22.png" /><figcaption>Figure 19: Sign of $h$ when $u_- > u_+$ and $f(u_-) < f(u_+)$. Left: admissible case when $\nu_- \ge \nu_+$. Right: non admissible case when $\nu_- < \nu_+$.</figcaption></figure>

boundary of $\mathcal{L}$, that is $w_+ = u_+ - a$, otherwise $R_6 \cap \mathcal{L} \neq \emptyset$, contradicting $h \ge 0$. However, in $R_5$ we have

$$h(k, \hat{k}) = (f(u_{-}) - f(u_{+}))[u_{+} + u_{-} - 2k + w_{+} + w_{-} - 2\hat{k}] \\ - (u_{-} - u_{+} + w_{-} - w_{+})[f(u_{+}) + f(u_{-}) - 2f(k)],$$

and we can prove that there exists a $C^2$ function $\varphi : [u_{-}, k_2[ \to [w_{-}, w_{+}[ \text{ such that } \varphi(u_{-}) = w_{-}, h(k, \varphi(k)) = 0 \text{ and}$

$$\varphi''(k) = \frac{(u_{-} - u_{+} + w_{-} - w_{+})f''(k)}{f(u_{-}) - f(u_{+})} < 0.$$

So the graph of $\varphi$ is concave, meaning that, even in the best scenario, that is if $k_2 = u_+$ and $\varphi(u_+) = w_+$, we would have $R_5 \cap \mathcal{L} \cap \{h > 0\} \neq \emptyset$, see Figure 20, left. So this case is not admitted by (4.3).

**Case $u_- < u_+$ and $f(u_-) < f(u_+$):** also this case is incompatible with (4.3). Indeed, we infer again that $(u_+, w_+)$ must lie on the upper boundary of $\mathcal{L}$ and that instead in $R_5$ we have a convex level set 0 of the function $h$. This makes it impossible to have $\{h \ge 0\} \cap \mathcal{L} \neq \emptyset$, contradicting (4.3).

To conclude we notice that if $u_- < u_+$, no entropy solution is admissible, so $u_- \ge u_+$ is a necessary condition. $\square$

Then the following stability result holds, which implies uniqueness of entropy weak solutions.

**Theorem 4.3.** Consider the Cauchy problems with initial conditions $(u_0^1, w_0^1)$ and $(u_0^2, w_0^2)$ respectively, where $u_0^i, w_0^i \in \mathbf{L}^1(\mathbb{R}) \cap \text{BV}(\mathbb{R})$. Let us denote by $(u_1, w_1)$ and $(u_2, w_2)$ two entropy weak solutions of the corresponding Cauchy problems. Then, it holds

$$\int_{-\infty}^{+\infty} (|u_1 - u_2|(x, t) + |w_1 - w_2|(x, t)) dx \le \int_{-\infty}^{+\infty} (|u_0^1 - u_0^2| + |w_0^1 - w_0^2|) dx, \quad (4.6)$$

for almost every $t \in [0, T[$.
<figure><img src="image_23.png" /><figcaption>Figure 20: Left: sign of the function $h$ when $u_{-} < u_{+}$ and $f(u_{-}) > f(u_{+})$; independently from $(u_{-}, w_{-})$, $(u_{+}, w_{+}) \in \mathcal{L}$, $\mathcal{L} \cap \{h > 0\} \neq \emptyset$. Right: $u_{-} < u_{+}$ and $f(u_{-}) < f(u_{+})$; $\mathcal{L} \cap \{h > 0\} \neq \emptyset$, so when $u_{-} < u_{+}$ no shock wave solutions are admissible.</figcaption></figure>

The proof is an adaptation of the standard doubling of variables method by Kruzkov [20]. In [28] an analogous result for the case of delayed-relay hysteresis is proven, see also [27, Section IV.1] for the description of such operator. The adjustment to Play hysteresis is straightforward.

## A Details of the proof of the Theorem 3.12

We show that the sequences $u_m, w_m$ defined in Theorem 3.12 converge in $L_{\text{loc}}^1(\mathbb{R} \times [0, T])$ to a couple of functions $u, w \in C^0([0, T[; L_{\text{loc}}^1(\mathbb{R}))$ which satisfy (1.6).

*Proof.* Using a classical argument (see e.g. [6, Theorem 2.6]), by Propositions 3.6 and 3.7, Lemma 3.8 and as a consequence of Helly's compactness theorem, there exists $u, w \in C^0([0, T[; L_{\text{loc}}^1(\mathbb{R}))$ such that (up to a subsequence)

$$u_m, w_m \to u, w \quad \text{in } L_{\text{loc}}^1(\mathbb{R} \times [0, T]), \text{ as } m \to +\infty.$$

Now, following the proof of [13, Theorem 5.3], we show that $u, w$ satisfies (1.6). We consider $\{(u_i^n, w_i^n) | i \in \mathbb{Z} \text{ and } n \in \mathbb{N}\}$ given by the scheme (3.12) associated to $\Delta x_m$ and $\Delta t_m$, and we fix $\phi \in C_c^1(\mathbb{R} \times [0, T])$ non negative. Notice that $\forall i \in \mathbb{Z}, n \in \mathbb{N}, u_i^n$ and $w_i^n$ depend also on $m$, but to easy the notation, we will omit this index. From Proposition 3.9, the discrete entropy condition (3.44) holds for every $m$. Multiplying then this inequality by $\phi(x, t^n)$, integrating it over $x \in K_i$ and summing over $i$ and $n$, we get

$$A'_m + A''_m + B_m \le 0, \qquad (A.1)$$
where

$$
\begin{align*}
A'_m &:= \sum_{i \in \mathbb{Z}} \sum_{n \in \mathbb{N}} (|u_i^{n+1} - k| - |u_i^n - k|) \int_{K_i} \phi(x, t^n) dx, \\
A''_m &:= \sum_{i \in \mathbb{Z}} \sum_{n \in \mathbb{N}} (|w_i^{n+1} - \hat{k}| - |w_i^n - \hat{k}|) \int_{K_i} \phi(x, t^n) dx, \\
B_m &:= \sum_{i \in \mathbb{Z}} \sum_{n \in \mathbb{N}} (G_k(u_i^n, u_{i+1}^n) - G_k(u_{i-1}^n, u_i^n)) \frac{\Delta t_m}{\Delta x_m} \int_{K_i} \phi(x, t^n) dx,
\end{align*}
$$

with $G_k(\alpha, \beta) := g(\alpha \top k, \beta \top k) - g(\alpha \perp k, \beta \perp k)$, (recall (3.3)).

Regarding $A'_m$, we can see that

$$
\begin{align*}
& - \int_{\mathbb{R}} \int_{\Delta t_m}^{T} |u_m(x,t) - k| \partial_t \phi(x, t - \Delta t_m) dt dx - \sum_{i \in \mathbb{Z}} |u_i^0 - k| \int_{K_i} \phi(x, 0) dx \\
&= - \sum_{n=1}^{+\infty} \sum_{i \in \mathbb{Z}} \iint_{K_i^n} |u_i^n - k| \partial_t \phi(x, t - \Delta t_m) dt dx - \sum_{i \in \mathbb{Z}} |u_i^0 - k| \int_{K_i} \phi(x, 0) dx \\
&= - \sum_{n=1}^{+\infty} \sum_{i \in \mathbb{Z}} \int_{K_i} |u_i^n - k| (\phi(x, t^n) - \phi(x, (n-1)\Delta t_m)) dx - \sum_{i \in \mathbb{Z}} |u_i^0 - k| \int_{K_i} \phi(x, 0) dx \\
&= - \sum_{n \in \mathbb{N}} \sum_{i \in \mathbb{Z}} (|u_i^n - k| - |u_i^{n+1} - k|) \int_{K_i} \phi(x, t^n) dx = A'_m.
\end{align*}
$$

Hence, from this equality, as $u_m \to u$ in $\mathbf{L}_{\text{loc}}^{1}(\mathbb{R} \times [0, T])$, $\sum u_i^0 1_{K_i} \to u_0$ in $\mathbf{L}_{\text{loc}}^{1}(\mathbb{R})$, $\phi$ has compact support and $\partial_t \phi(\cdot, \cdot - \Delta t_m) \to \partial_t \phi$ in $\mathbf{L}^\infty(\mathbb{R} \times [0, T])$ by uniform continuity, we get

$$
\lim_{m \to \infty} A'_m = - \int_0^T \int_{-\infty}^{+\infty} |u - k| \partial_t \phi \, dx \, dt - \int_{-\infty}^{+\infty} |u_0(x) - k| \phi(x, 0) \, dx. \quad (A.2)
$$

By the same reasoning, it also holds

$$
\lim_{m \to \infty} A_m'' = - \int_0^T \int_{-\infty}^{+\infty} |w - \hat{k}| \partial_t \phi \, dx \, dt - \int_{-\infty}^{+\infty} |w_0(x) - \hat{k}| \phi(x, 0) \, dx. \quad (A.3)
$$

To compute instead the limit of $B_m$, we introduce the following quantity $B_{1,m}$

$$
B_{1,m} := - \sum_{n \in \mathbb{N}} \int_{t^n}^{t^{n+1}} \int_{-\infty}^{+\infty} G_k(u_m, u_m) \partial_x \phi(x, t^n) dx dt.
$$

By the strong $\mathbf{L}_{\text{loc}}^1(\mathbb{R} \times [0, T])$ convergence for $u_m$, the Lipschitz continuity of $G_k$ (this is a consequence of the Lipschitz continuity of $g$ defined by (3.3)) and the convergence of $\partial_x \phi(\cdot, n\Delta t) 1_{[n\Delta t_m, (n+1)\Delta t_m]}(\cdot)$ in $\mathbf{L}^\infty(\mathbb{R} \times [0, T])$, we can immediately notice that

$$
\lim_{m \to \infty} B_{1,m} = - \int_0^T \int_{-\infty}^{+\infty} G_k(u, u) \partial_x \phi \, dx \, dt = - \int_0^T \int_{-\infty}^{+\infty} \text{sign}(u-k)(f(u)-f(k)) \partial_x \phi \, dx \, dt.
$$
We then have just to compare $B_m$ to $B_{1,m}$, so we rewrite the latter as follows

$$
\begin{align*}
B_{1,m} &= -\Delta t_m \sum_{i \in \mathbb{N}} \sum_{i \in \mathbb{Z}} G_k(u_i^n, u_i^n)(\phi(x_{i+1/2}, t^n) - \phi(x_{i-1/2}, t^n)) \\
&= -\Delta t_m \sum_{i \in \mathbb{N}} \sum_{i \in \mathbb{Z}} (G_k(u_{i-1}^n, u_{i-1}^n) - G_k(u_i^n, u_i^n))\phi(x_{i-1/2}, t^n) \\
&= \Delta t_m \sum_{i \in \mathbb{N}} \sum_{i \in \mathbb{Z}} (G_k(u_i^n, u_i^n) - G_k(u_{i-1}^n, u_i^n))\phi(x_{i-1/2}, t^n) + \\
&\quad + \Delta t_m \sum_{i \in \mathbb{N}} \sum_{i \in \mathbb{Z}} (G_k(u_{i-1}^n, u_i^n) - G_k(u_i^n, u_{i-1}^n))\phi(x_{i-1/2}, t^n)
\end{align*}
$$

Similarly, we rewrite

$$
\begin{align*}
B_m &= \sum_{i \in \mathbb{Z}} \sum_{n \in \mathbb{N}} (G_k(u_i^n, u_{i+1}^n) - G_k(u_i^n, u_i^n) + G_k(u_i^n, u_i^n) - G_k(u_{i-1}^n, u_i^n)) \frac{\Delta t_m}{\Delta x_m} \int_{K_i} \phi(x, t^n) dx \\
&= \sum_{i \in \mathbb{Z}} \sum_{n \in \mathbb{N}} (G_k(u_i^n, u_i^n) - G_k(u_{i-1}^n, u_i^n)) \frac{\Delta t_m}{\Delta x_m} \int_{K_i} \phi(x, t^n) dx + \\
&\quad + \sum_{i \in \mathbb{Z}} \sum_{n \in \mathbb{N}} (G_k(u_{i-1}^n, u_i^n) - G_k(u_i^n, u_i^n)) \frac{\Delta t_m}{\Delta x_m} \int_{K_i} \phi(x, t^n) dx.
\end{align*}
$$

Then

$$
|B_{1,m} - B_m| \le C_1 C_2 \Delta t_m \Delta x_m \sum_{i=i_0}^{i_1} \sum_{n=0}^{N} |u_i^n - u_{i-1}^n|,
$$

where $i_0, i_1$ and $N$ (all dependent on $m$) are the indices for which the support of $\phi$ is
contained in $[i_0 \Delta x_m, i_1 \Delta x_m] \times [0, (N+1) \Delta t_m[$, $C_1$ is a constant given by the Lipschitz
continuity of $G_k$ on $[U_m, U_M]$ and $C_2$ is a constant given as a consequence of the $C^1$
regularity of $\phi$. Indeed, $C_2$ uniform in $i, n$ and $m$, can be found such that

$$
\phi(x_{i-1/2}, t^n) - C_2 \Delta x_m \leq \frac{1}{\Delta x_m} \int_{K_i} \phi(x, t^n) dx \leq \phi(x_{i-1/2}, t^n) + C_2 \Delta x_m.
$$

In particular, by Frechet-Kolmogorov's compactness theorem for $\mathbf{L}^p$ spaces (see e.g. [7]),
taking as translation parameter $h_m := \Delta x_m$, we have, as $\Delta x_m \to 0$, that

$$
0 = \lim_{m \to \infty} \iint_{\text{spt}(\phi)} |u_m(x, t) - u_m(x - \Delta x_m, t)| \, dx \, dt = \lim_{m \to \infty} \Delta t_m \Delta x_m \sum_{i=i_0}^{i_1} \sum_{n=0}^{N} |u_i^n - u_{i-1}^n|,
$$

meaning that

$$
\lim_{m \to \infty} B_m = \lim_{m \to \infty} B_{1,m} = - \int_0^T \int_{-\infty}^{+\infty} \text{sign}(u-k)(f(u)-f(k)) \partial_x \phi \, dx dt. \quad (A.4)
$$

Finally, taking the limit as $m$ goes to infinity in (A.1), from (A.2), (A.3) and (A.4) we
obtain (1.6).

<figure><img src="image_24.png" /></figure>
# Acknowledgments

This work was mainly written while Stefan Moreti was visiting the ACUMES team at Inria Centre at Université Côte d'Azur in Sophia Antipolis, France. SM was partially supported by the INdAM - GNAMPA Project, CUP E53C25002010001: "Analisi e controllo per alcuni problemi di evoluzione". He also gratefully acknowledges the Inria financial support during his research stay, as well as the hospitality of the ACUMES team.

# References

[1] D. Amadori, A. Bressan, and W. Shen. Conservation Laws with Discontinuous Gradient-Dependent Flux: the Unstable Case. Preprint, arXiv:2411.13444, 2024.

[2] D. Amadori, A. Bressan, and W. Shen. Conservation laws with discontinuous gradient-dependent flux: the stable case. *Math. Models Methods Appl. Sci.*, 35(6):1421–1469, 2025.

[3] L. Ambrosio, N. Fusco, and D. Pallara. *Functions of bounded variation and free discontinuity problems*. Oxford Math. Monogr. Oxford: Clarendon Press, 2000.

[4] F. Bagagiolo and S. Moreti. Wave-front tracking for a quasi-linear scalar conservation law with hysteresis. *J. Math. Anal. Appl.*, 543(2):Paper No. 128900, 25, 2025.

[5] F. Bagagiolo and S. Moreti. Wave-front tracking for a quasi-linear scalar conservation law with hysteresis II: the case of Preisach. *J. Math. Anal. Appl.*, 558(1):130345, 2026.

[6] A. Bressan. *Hyperbolic systems of conservation laws, volume 20 of Oxford Lecture Series in Mathematics and its Applications*. Oxford University Press, Oxford, 2000. The one-dimensional Cauchy problem.

[7] H. Brezis. *Functional analysis, Sobolev spaces and partial differential equations*. New York, NY: Springer, 2011.

[8] M. Brokate and J. Sprekels. *Hysteresis and phase transitions*. New York Springer, 1996.

[9] A. Corli and H. Fan. Hysteresis and stop-and-go waves in traffic flows. *Mathematical Models and Methods in Applied Sciences*, 29:2637–2678, 2019.

[10] A. Corli and H. Fan. String stability in traffic flows. *Applied Mathematics and Computation*, 443:24, 2023.

[11] A. Corli and H. Fan. The hysteretic Aw-Rascle-Zhang model. *Stud. Appl. Math.*, 153(4):38, 2024.

[12] L. C. Evans. *Partial differential equations*. American Mathematical Society, Providence, 2nd ed. edition, 2010.

[13] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In J. L. Lions and P. Ciarlet, editors, *Solution of equations in $\mathbb{R}^n$ (Part 3). Techniques of scientific computing (Part 3)*, volume 7 of *Handbook of Numerical Analysis*, pages 713–1020. Elsevier, 2000.

[14] H. Fan. Conservation laws with hysteretic fluxes. *Journal of Differential Equations*, 394:1–30, 2024.
[15] H. Holden and N. H. Risebro. Front tracking for hyperbolic conservation laws. Berlin: Springer, 2nd edition edition, 2015.

[16] J. Kopfová. Entropy condition for a quasilinear hyperbolic equation with hysteresis. *Differential and Integral Equations*, 18:451 – 467, 2005.

[17] P. Kordulová. Quasilinear hyperbolic equations with hysteresis. *Journal of Physics: Conference Series*, 55:135, 2006.

[18] M. Krasnosel'skii and A. Pokrovskii. *Systems with Hysteresis*. Berlin Springer, 1989.

[19] P. Krejci and P. Laurencot. Generalized variational inequalities. *J. Convex Anal.*, 9:159–183, 2002.

[20] S. N. Kruˇ zkov. First order quasilinear equations with several independent variables. *Mat. Sb. (N.S.)*, 81(123):228–255, 1970.

[21] R. J. LeVeque. *Numerical methods for conservation laws*. Basel: Birkhäuser, 2nd ed. edition, 1992.

[22] M. Peszynska and R. Showalter. Approximation of scalar conservation law with hysteresis. *SIAM Journal on Numerical Analysis*, 58:962–987, 2020.

[23] M. Peszyńska and S. R.E. A transport model with adsorption hysteresis. *Differential and Integral Equations*, 11:327 – 340, 1998.

[24] B. Plohr, D. Marchesin, P. Bedrikovetsky, and P. Krause. Modeling hysteresis in porous media flow via a relaxation. *Comput. Geosci.*, 5:225 – 256, 2001.

[25] V. Recupero. BV-extension of rate independent operators. *Math. Nachr.*, 282(1):86–98, 2009.

[26] E. F. Toro. *Riemann solvers and numerical methods for fluid dynamics. A practical introduction*. Berlin: Springer, 3rd ed. edition, 2009.

[27] A. Visintin. *Differential models of hysteresis*. Berlin: Springer, 1994.

[28] A. Visintin. Quasilinear first-order PDEs with hysteresis. *J. Math. Anal. Appl.*, 312:401–419, 2005.

[29] H. M. Zhang. A mathematical theory of traffic hysteresis. *Transportation Research Part B: Methodological*, 33(1):1–23, 1999.