C H Chalons

P Engel

C H Rohde
email: crohde@math.uni-stuttgart.de

A CONSERVATIVE AND CONVERGENT SCHEME FOR UNDERCOMPRESSIVE SHOCK WAVES

Undercompressive shock waves arise in numerous physical applications. We propose a class of conservative finite-volume type schemes to approximate weak solutions of conservation laws that contain undercompressive shock waves. We prove the convergence of a subsequence of approximate solutions towards a generalized entropy solution if the mesh width tends to zero. The proof relies on a refined BV compactness analysis, which accounts for the effect of the kinetic relation that drives the undercompressive wave. At the same time we establish a new proof for the existence of solutions to the underlying model. Numerical experiments supplement the analytical results.

1. Introduction. Let f ∈ C 2 (R), u 0 ∈ L ∞ (R), and T > 0 be given. We consider solutions u : R × [0, T ) → R of the scalar conservation law

u t + f (u) x = 0 in R × (0, T ), (1.1)
with initial condition

u(•, 0) = u 0 in R. (1.2)
If the flux is nonlinear, the Cauchy problem (1.1), (1.2) cannot be solved within the class of smooth functions. Therefore, weak solutions, i.e., distributional solutions in L ∞ (R × (0, T )) have to be considered. However, (1.1), (1.2) can have multiple weak solutions. Uniqueness is restored in the class of entropy solutions. That are weak solutions which satisfy the entropy inequality η(u) t + q(u) x ≤ 0 for all entropies η in D ′ (R × [0, T )).

(

Thereby, an entropy η ∈ C 2 (R) is a convex function and the entropy flux q ∈ C 2 (R) is given through q ′ = η ′ f ′ . Weak solutions of shock wave type that satisfy the entropy condition (1.3) are necessarily compressive. That means that characteristics impinge on the discontinuity from both sides. We refer to [START_REF] Dafermos | Hyperbolic conservation laws in continuum physics[END_REF] for a comprehensive treatment of the classical scalar well-posedness theory.

In a wide range of applications (cf. [START_REF] Abeyaratne | Kinetic relations and the propagation of phase boundaries in solids[END_REF][START_REF] Bertozzi | Undercompressive shocks in thin film flows[END_REF][START_REF] Colombo | Sonic and kinetic phase transitions with applications to Chapman-Jouguet deflagrations[END_REF][START_REF] Colombo | Pedestrian flows and non-classical shocks[END_REF][START_REF] Kluwick | Non-classical kinematic shocks in suspensions of particles in fluids[END_REF][START_REF] Van Duijn | A new class of entropy solutions of the Buckley-Leverett equation[END_REF] for an incomplete list of examples), discontinuous waves appear that take the form of undercompressive waves. In the scalar case of problem (1.1), undercompressive shock waves are characterized by the fact that characteristics pass through the shock wave. For the precise definition see Section 2 below. Note that this is only possible in the case that the flux in (1.1) has inflection points. As a matter of fact, undercompressive waves are ruled out by (1.3). In any case, physically relevant weak solutions with undercompressive waves should still be consistent with the second law of thermodynamics. Thus, they are supposed to satisfy η(u) t + q(u) x ≤ 0 for (at least) one entropy η in D ′ (R × [0, T )).

(1.4)

Obviously condition (1.4) does no ensure unique solvability of the Cauchy problem (1.1), (1.2). Other criteria are required.

An approach to guarantee unique solvability is embedding the first order equation (1.1) into an enlarged model of type

u ε t + f (u ε ) x = R[ε; u ε ] in R × (0, T ). (1.5)
Here, ε > 0 is a small parameter and R[ε; •] some operator with R[0; •] ≡ 0. In the limit ε → 0 the sequence of solutions {u ε } ε>0 to (1.5), (1.2) selects a unique weak solution for (1.1), (1.2) (provided solutions u ε exist and converge at all in an appropriate function space). To recover undercompressive waves that satisfy (1.4) in the limit, the operator R has to be chosen such that solutions of (1.5), (1.2) dissipate the entropy η. Examples typically take the form of higher-order differential or integral operators, and can be found in e.g. [START_REF] Bertozzi | Undercompressive shocks in thin film flows[END_REF][START_REF] Corli | Singular limits for a parabolic-elliptic regularization of scalar conservation laws[END_REF][START_REF] Hayes | Non-classical shocks and kinetic relations: scalar conservation laws[END_REF][START_REF] Rohde | Scalar conservation laws with mixed local and nonlocal diffusion-dispersion terms[END_REF][START_REF] Van Duijn | A new class of entropy solutions of the Buckley-Leverett equation[END_REF].

Another approach is motivated by the study of configurational forces for phase transformation in solids [START_REF] Abeyaratne | Kinetic relations and the propagation of phase boundaries in solids[END_REF][START_REF] Truskinovsky | Kinks versus shocks, in Shock induced transitions and phase structures in general media[END_REF]. This ansatz enforces uniqueness by restrictions on undercompressive shock waves through kinetic relations and nucleation criteria. A kinetic relation is an algebraic condition on the traces across an undercompressive shock wave. In view of (1.4) and the application background, kinetic functions are usually chosen to determine exactly the amount of entropy dissipation with respect to η. For a discussion of both methods and its analytical implications we refer to the book [START_REF] Lefloch | The theory of classical and nonclassical shock waves[END_REF].

The numerical approximation of undercompressive waves is a nontrivial issue. On the one hand the exact solutions are neither monotone nor total-variation diminishing [START_REF] Lefloch | The theory of classical and nonclassical shock waves[END_REF]. As a consequence, the highly successful finite volume schemes with these properties are ruled out from the beginning. On the other hand, one of the selection mechanisms outlined above has to be accounted for on the discrete level. Following the idea of embedding, one discretizes first (1.5) and then identifies the role of ε with the spatial mesh parameter. Thus, one obtains a consistent discretization for (1.1) that mimics the effect of the operator R in (1.5) if the mesh width vanishes. The schemes from the papers [START_REF]Nonclassical shocks and kinetic relations: finite difference schemes[END_REF][START_REF] Kissling | The computation of nonclassical shock waves with a heterogeneous multiscale method[END_REF][START_REF] Lefloch | Fully discrete, entropy conservative schemes of arbitrary order[END_REF][START_REF] Lefloch | Why many theories of shock waves are necessary: kinetic functions, equivalent equations, and fourth-order models[END_REF][START_REF] Lefloch | High-order schemes, entropy inequalities, and nonclassical shocks[END_REF] fall into this class. Up to our knowledge, the proof of convergence for these methods is an open issue. The other class of schemes relies on the kinetic relation directly. The first of such schemes is the Glimm-type ansatz in [START_REF] Lefloch | Propagating phase boundaries: formulation of the problem and existence via the Glimm method[END_REF], see also [START_REF] Chalons | Computing undercompressive waves with the random choice scheme[END_REF]. Deterministic versions that use an extra tracking of the undercompressive waves have been introduced in [START_REF] Chalons | Fast relaxation solvers for hyperbolic-elliptic phase transition problems[END_REF][START_REF] Thomas | A level-set approach to the computation of twinning and phase-transition dynamics[END_REF][START_REF] Merkle | Computation of dynamical phase transitions in solids[END_REF][START_REF] Merkle | The sharp-interface approach for fluids with phase change: Riemann problems and ghost fluid techniques[END_REF][START_REF] Zhong | Computational methods for propagating phase boundaries[END_REF]. The drawback of all these schemes is the fact that the discrete solution does not conserve the integral of the solution. To overcome this problem, conservative finite-volume type methods have been suggested in [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF][START_REF] Perrier | A conservative method for the simulation of the isothermal Euler system with the van-der-Waals equation of state[END_REF]. These schemes use a special numerical flux close to undercompressive waves (and standard monotone fluxes away from the wave). From the point of view of error analysis these schemes can usually be proven to converge if the initial data corresponds to a single undercompressive wave that satisfies the given kinetic relation.

In this contribution we introduce a new conservative method. It imposes the kinetic relation directly, and thus belongs to the second class of schemes described above. We prove the convergence towards a generalized entropy solution. The proof applies to solutions which contain a finite number of undercompressive waves but is not restricted to well-prepared initial data as in many previous works. Our proof does not only yield the weak solution, but does also identify the location of the phase boundary, the associated entropy dissipation, and its characteristic structure.

In Section 2 we give a rigorous definition of the generalized entropy solutions we are searching for. We introduce the novel method in Section 3. In Section 4 we state the convergence theorem (Theorem 4.4). Section 5 is devoted completely to the proof. The last section contains numerical experiments that illustrate the analytical findings and the future potential of the method.

2. The Mathematical Model. We assume that the flux f ∈ C 2 (R) has a concave-convex shape with inflection point at u = 0, i.e.,

uf ′′ (u) > 0 for all u ∈ R \ {0} and f ′′ (0) = 0. (2.1)
With this assumption and having in mind the application scenario of phase transformations we define two phases by

Z + = {u > 0} and Z -= {u < 0}.
By a phase boundary we mean a function u = u(x, t) of shock wave type, i.e.,

u(x, t) =    u -: x -āt < 0 u + : x -āt > 0 ā = ā(u -, u + ) := f (u -) -f (u + ) u --u + , (2.2)
such that either u ± ∈ Z ± or u ± ∈ Z ∓ holds. The choice of ā guarantees the validity of the Rankine-Hugoniot conditions. Phase boundaries can be of compressive/Laxian type if

f ′ (u + ) ≤ ā ≤ f ′ (u -) (2.3)
holds, or of undercompressive type if either

f ′ (u + ) < ā, ā > f ′ (u -) or f ′ (u + ) > ā, ā < f ′ (u -) (2.4)
hold. Note that by (2.3) also the (one-sided) characteristic case is included for the sake of simplicity. We are in particular interested in weak solutions of (1.1) that contain undercompressive phase boundaries. For the initial data we assume that u 0 ∈ BV ∩L 1 (R) and

xu 0 (x) < 0 (2.5)
holds for almost all x ∈ R. The property (2.5) implies that exactly one phase change at x = 0 is prescribed by u 0 . The extension of the analysis of this paper to initial data with (finitely) many phase boundaries is possible but would lead to a much more heavy notational overhead. In view of phase transformation applications the assumption (2.1) on the flux function appears to be a generic situation [START_REF] Abeyaratne | Kinetic relations and the propagation of phase boundaries in solids[END_REF]. The treatment of the interesting case of flux functions with multiple inflection points and thus problems with multiple phases is out of the scope of this paper.

As the next step, we introduce a generalized notion of entropy solution for weak solutions that contain undercompressive waves. Let us recall that a weak solution of (1.1) is a function u ∈ L ∞ (R × (0, T )) such that

R×R + uψ t + f (u) ψ x dx dt = - R ψ(x, 0)u 0 (x) dx (2.6)
holds for all ψ ∈ C 1 0 (R × [0, T )). Rewriting (1.4) leads to

R×R + η(u)ψ t + q(u)ψ x dx dt ≥ - R ψ(x, 0)η u 0 (x) dx (2.7)
for all non-negative functions ψ ∈ C 1 0 (R × [0, T )). In recent years several well-posedness theories for weak solutions have been developed that allow for undercompressive waves [START_REF] Lefloch | The theory of classical and nonclassical shock waves[END_REF]. We detail here the concept of a kinetic relation.

Definition 2.1 (Kinetic Relation). A function ϕ :

Z + → Z -is called kinetic relation if ∃ L > 0 : ϕ(u) -ϕ(u ′ ) ≤ L u -u ′ ∀ u, u ′ ∈ Z + , (2.8a) ϕ is decreasing, (2.8b) ϕ ♯ (u) ≥ ϕ(u) and ϕ(0) = 0 (2.8c) hold. In (2.8c) the function ϕ ♯ : Z + → Z -is defined by ā u, ϕ ♯ (u) = ā u, ϕ(u) and ϕ ♯ (u) = ϕ(u).
(2.9)

According to (2.1) the function ϕ ♯ is well-defined.

A kinetic relation ϕ :

Z + → Z -is called compatible with an entropy η ∈ C 2 (R) if ā(u, ϕ(u)) η(u) -η(ϕ(u)) -q(u) -q(ϕ(u)) ≤ 0 (2.10)
holds for all u ∈ Z + , with ā as in (2.2). To illustrate the interplay of kinetic relations and the notion of entropy compatibility consider the following example that will also be used later in Section 6. Example 2.2. Let f (u) = 1 3 u 3 + cu with c ∈ R and the kinetic relation ϕ(u) = -κu with κ ∈ [ 1 2 , 1] be given. Consider the entropy pairs (η, q) given by

η(u) = 1 r u r and q(u) = 1 r + 2 u r+2 + c r u r , r 2 ∈ N. (2.

11)

The left-hand side of (2.10) is computed as

2(r -1) κ r+2 -1 + (r + 2) κ r+1 -κ r + κ 2 -κ 3r(r + 2) u r+2 ≤ 0,
thus ϕ is compatible with all entropies from (2.11).

In the following we will construct an entropy which is not compatible with ϕ. We choose for k ∈ R the Kruzkov entropy η(u) = |u -k| with entropy flux q(u) = sign(u-

k) f (u) -f (k) . Fix u ∈ Z + and take k ∈ [ϕ(u), u], then (2.10) is equivalent to t(k) = ā(u, ϕ(u)) u + ϕ(u) -2k -f (u) + f (ϕ(u)) -2f (k) ≤ 0.
If we maximize t(k) for k ∈ [ϕ(u), u] we find the maximum at k * = -

1 3 √ 3κ 2 -3κ + 3u and t(k * ) > 0 for κ ∈ [ 1 2 , 1].
In the sequel we search for functions in the space

V := v ∈ C 0 [0, T ]; L 1 (R) v(•, t) ∈ BV(R), ∀t ∈ [0, T ] .
(2.12)

We can now present our generalized notion of entropy solutions. Definition 2.3 ((η, ϕ)-solution). Let an entropy η ∈ C 2 (R) and a kinetic relation ϕ be given such that ϕ is compatible with η. A weak solution u ∈ V of (1.1) together with p ∈ C 0 (0, T ) is called an (η, ϕ)-solution of (1.1) if it satisfies the following properties:

1. The function u fulfills the entropy condition (2.7) for η. 2. p is differentiable for almost all t ∈ (0, T ) and fulfills there

p t (t) = ā u -(t), u + (t) .
(2.13)

3.

For almost all t ∈ (0, T ) we have waves. As we will see later in Example 6.4 a phase boundary can change its type which has to be reflected in the formulation of the weak solution: the first condition in (2.14) identifies the discontinuity across p as an undercompressive case, the second one as Laxian wave (see [START_REF] Lefloch | The theory of classical and nonclassical shock waves[END_REF]). (iv) The limits in (2.15) are well-defined, since u(•, t) ∈ BV(R) for all t ∈ [0, T ]. (v) In view of assumption (2.5) we will have u -(t) ∈ Z + in (2.14). In the general case the inequality in (2.14) should be written as

u + (t) = ϕ(u -(t)) or u + (t) ≥ ϕ ♯ (u -(t)). ( 2
u -(t)u + (t) ≥ u -(t)ϕ ♯ (u -(t)).
We believe that the notion of weak solution from Definition 2.3 can be used to prove the well-posedness of the initial value problem (1.1). Rigorous proofs for existence and uniqueness can be found in [START_REF] Lefloch | The theory of classical and nonclassical shock waves[END_REF]. Let us discuss here the Riemann problem for the sake of illustration. Consider for u L > 0 and u R ∈ R the Riemann initial data

u 0 (x) = u L for x < 0, u R for x > 0. (2.16)
The cases u R ≥ u L and u R ∈ [0, u L ) correspond to one-phase initial data. They are solved by a single rarefaction wave and a single Laxian shock wave, respectively. We concentrate on two-phase initial datum (cf. (2.5)).

Example 2.5. Recalling the results from [START_REF] Lefloch | The theory of classical and nonclassical shock waves[END_REF] we know that solutions of the Riemann problem (1.1), (2.16) have the following structure for a given kinetic relation ϕ :

Z + → Z -: (i) u R ∈ [ϕ ♯ (u L ), 0): The solution u(x, t) = u L for x t ≤ ā(u L , u R ) u R for ā(u L , u R ) ≤ x t is a Laxian wave. (ii) u R ∈ [ϕ(u L ), ϕ ♯ (u L )): The solution u(x, t) =        u L for x t ≤ ā(u L , ϕ(u L )) ϕ(u L ) for ā(u L , ϕ(u L )) ≤ x t ≤ ā(ϕ(u L ), u R ) u R for ā(ϕ(u L ), u R ) ≤ x t
is an undercompressive phase boundary followed by a shock wave.

(iii) u R < ϕ(u L ): The solution u(x, t) =              u L for x t ≤ ā(u L , ϕ(u L )) ϕ(u L ) for ā(u L , ϕ(u L )) ≤ x t ≤ f ′ (ϕ(u L )) ( f ′ Z-) -1 ( x t ) for f ′ (ϕ(u L )) ≤ x t ≤ f ′ (u R ) u R for f ′ (u R ) ≤ x t
is an undercompressive phase boundary followed a rarefaction wave. Note that the solution from (ii) is neither monotone nor total variation diminishing. Both properties hold for classical entropy solutions [START_REF] Dafermos | Hyperbolic conservation laws in continuum physics[END_REF]. It is easy to check that the function u in (i), (ii), and (iii) is an (η, ϕ)-solution for all entropies (2.11).

The Numerical Method.

In this section we propose a numerical method to approximate (η, ϕ)-solutions in the sense of Definition 2.3. The approximation of (η, ϕ)-solutions requires a very accurate flux balancing close to the undercompressive waves. In particular the numerical method must detect two-wave structures as in cases (ii) and (iii) in Example 2.5. Standard finite volume schemes fail in these situations. Therefore, the new explicit time-stepping method uses a background grid in the regions away from the phase boundary. At the phase boundary we cut the space-time cells into two smaller cells according to its exact position. As this might lead to extremely small cells and corresponding vanishing time steps, the cutted cells are merged with the closest spatial background grid cell in the same phase. Then we use this time-space calculation grid to perform an explicit time stepping with exact flux balancing at the interface and standard fluxes away from the phase boundary. Finally, for the new time step, the results are projected to the spatial part of the background grid. The approach is in the spirit of moving mesh methods [START_REF] Sanders | The moving grid method for nonlinear hyperbolic conservation laws[END_REF]. We give some illustration in Figure 3.1. The background grid is plotted with dotted lines, the space calculation grid with solid lines. Note that the cell edges of the calculation grid might vary during a single time step.

For the formal presentation of the method in Algorithm 3.1 below, consider a uniform (for the sake of simplicity) background grid on R×[0, T ) with mesh parameter A possible fifth case is described as the mirrored case to h 4 .

x m(n)-2 x m(n)-1 x m(n) x m(n)+1 x m(n)+2 x m(n)+3 t n t n+1 t n+2 t n+3 t n+4 p n p n+1 p n+2 p n+3 p n+4 j 1 j 2 j 3 j 4
h > 0. It consists for spatial cells

I j := [x j-1 2 , x j+ 1 2 ) of elements I j × [t n , t n+1 ) with x j+ 1 2 := j + 1 2 h (j ∈ Z) and 0 = t 0 < t 1 < • • • < t N := T for some N ∈ N.
The time steps will be related to the spatial mesh parameter by a CFL-like condition.

The spatial grid midpoints are given by x j := jh.

The discrete position of the phase boundary at time t n is denoted as p n ∈ R. We represent p n uniquely with numbers m(n) ∈ Z and α n ∈ [0, 1) by the formula

p n = m(n) + α n - 1 2 h. (3.1)
Thus m(n) denotes the index of the spatial cell, where the discrete phase boundary is located, i.e., p n ∈ I m(n) . At t = 0 we have p 0 = 0 and therefore m(0) = 0 and α 0 = 1 2 . Furthermore, we use the following notations: We define ∆m(n) := m(n + 1)m(n) and the set J(n) by

J(n) := Z \ m(n) -1, m(n), m(n) + 1 .
We write ∆m + (n) := max(0, ∆m(n)) and ∆m -(n) :=min(0, ∆m(n)). Last, we introduce the indices m

-(n) = m(n) -∆m -(n) and m + (n) = m(n) + ∆m + (n).
We define the initial values for j ∈ J(0) by integral averages, i.e.,

u 0 j := 1 h x j+ 1 2 x j-1 2 u 0 (x) dx (3.2a) and u 0 1 , u 0 -1 by u 0 1 := 2 3h 3 2 h 0 u 0 (x) dx , u 0 -1 := 2 3h 0 -3 2 h u 0 (x) dx . (3.2b)
Note that u 0 0 is not initialized. Furthermore, we ensure that u 0 j ∈ Z ± ⇔ j ∈ Z ∓ . Let F : R 2 → R be a numerical flux, i.e., F (u, u) = f (u) for u ∈ R. Finally we select an entropy η ∈ C 2 and a kinetic relation ϕ : Z + → Z -which is compatible with η. Algorithm 3.1. Let {u n j } j and p n be given for n ∈ N 0 . Then {u n+1 j } j and p n+1 are calculated by the following procedure. 1. Consider the (η, ϕ)-solution of (1.1), (1.2), (2.16)

with u L = u n m(n)-1 and u R = u n m(n)+1
. Define the velocity s n of the phase boundary by

s n := ā(ϕ(u n m(n)-1 ), u n m(n)-1 ) if u n m(n)+1 < ϕ ♯ (u n m(n)-1 ), ā(u n m(n)+1 , u n m(n)-1 ) if u n m(n)+1 ≥ ϕ ♯ (u n m(n)-1 ), (3.3)
and the numerical phase boundary flux h n by

h n := f (u n m(n)-1 ) -s n u n m(n)-1 . (3.4)
2. Define the time step ∆t through

∆t := h 8 Ĉ . (3.5)
For the definition of the constant Ĉ > 0 see Assumption 4.2 below. 3. Calculate m(n + 1) ∈ Z and α n+1 ∈ [0, 1) from α n and s n such that

m(n + 1) + α n+1 = m(n) + α n + ∆t h s n (3.6)
holds. 4. Calculate {u n+1 j } j by the following formulas. For j ∈ J(n) ∩ J(n + 1) (type j 1 in Figure 3.1) we use

u n+1 j := u n j - ∆t h f n j+ 1 2 -f n j-1 2 , f n j+ 1 2 = F (u n j , u n j+1 ). (3.7)
We give the formulas for the other cells for the different cases ∆m(n) ∈ {-1, 0, 1} separately: For ∆m(n) = 0 (type j 2 and j 3 in Figure 3.1) we have

u n+1 m(n+1)-1 := (1 + α n ) u n m(n)-1 -∆t h h n -f n m(n)-3 2 1 + α n+1 , u n+1 m(n+1)+1 := (2 -α n ) u n m(n)+1 -∆t h f n m(n)+ 3 2 -h n 2 -α n+1 . (3.8)
For ∆m(n) = 1 (type j 2 and j 4 in Figure 3.1) we have

u n+1 m(n+1)-2 := u n+1 m(n+1)-1 := (1 + α n ) u n m(n)-1 -∆t h h n -f n m(n)-3 2 2 + α n+1 , u n+1 m(n+1)+1 := (2 -α n ) u n m(n)+1 + u n m(n)+2 -∆t h f n m(n)+ 5 2 -h n 2 -α n+1 .
(3.9)

For ∆m(n) = -1 (type j 5 and j 3 in Figure 3.1) we get

u n+1 m(n+1)-1 := (1 + α n ) u n m(n)-1 + u n m(n)-2 -∆t h h n -f n m(n)-5 2 1 + α n+1 , u n+1 m(n+1)+1 := u n+1 m(n+1)+2 := (2 -α n ) u n m(n)+1 -∆t h f n m(n)+ 3 2 -h n 3 -α n+1 .
(3.10)

The three expressions (3.8), (3.9), and (3.10) can be summarized in the compact form

u n+1 m(n+1)-1 u n+1 m -(n)-1    := (1 + α n )u n m(n)-1 + ∆m -(n)u n m(n)-2 -∆t h h n -f n m -(n)-3 2 1 + ∆m + (n) + α n+1 (3.11a)
and

u n+1 m(n+1)+1 u n+1 m + (n)+1    := (2 -α n )u n m(n)+1 + ∆m + (n)u n m(n)+2 -∆t h f n m + (n)+ 3 2 -h n 2 + ∆m -(n) -α n+1 . (3.11b)
Finally, let us define the approximate (η, ϕ)-solution u h : R × [0, T ] → R for x ∈ I j and t ∈ [t n , t n+1 ) through

u h (x, t) =      u n j for j ∈ J(n), u n m(n)-1 for j / ∈ J(n) and x < p n + s n (t -t n ), u n m(n)+1 for j / ∈ J(n) and x ≥ p n + s n (t -t n ).
(3.12)

The symbol u n h denotes u h (•, t n ). As the first simple but striking property of Algorithm 3.1 we get Proposition 3.2. Algorithm 3.1 is conservative, i.e., for all n

∈ N R u h (x, t n ) dx = R u 0 (x) dx .
Proof. From (3.1), (3.6), (3.7), (3.11), and (3.12) we get

R u h (x, t n+1 ) dx = h   j∈J(n+1) u n+1 j + 1 + α n+1 u n+1 m(n+1)-1 + 2 -α n+1 u n+1 m(n+1)+1   = h   j∈J(n)∩J(n+1) u n+1 j + 1 + ∆m + (n) + α n+1 u n+1 m(n+1)-1 + 2 + ∆m -(n) -α n+1 u n+1 m(n+1)+1   = h   j∈J(n)∩J(n+1) u n j + ∆m -(n)u n m(n)-2 + (1 + α n ) u n m(n)-1 + (2 -α n ) u n m(n)+1 + ∆m + (n)u n m(n)+2   = h   j∈J(n) u n j + (1 + α n ) u n m(n)-1 + (2 -α n ) u n m(n)+1   = R u h (x, t n ) dx ,
which implies the statement by induction. Remark 3.3. Thanks to the time step restriction (3.5) we have ∆m(n) ∈ {-1, 0, 1}. Therefore, ∆m

+ (n), ∆m -(n) ∈ {0, 1}. We interpret ∆m + (n) (∆m -(n))
as a trigger telling whether the phase boundary moves one cell to the right (left) or not. Thus the following equalities hold:

∆m + (n)u n m(n)+2 = ∆m + (n)u n m + (n)+1 , (1 -∆m + (n))u n m(n)+1 = (1 -∆m + (n))u n m + (n)+1 .
4. The Main Theorem. In this section we state the main result. The proof will be given in the subsequent Section 5. We summarize the assumptions on the model and the scheme.

Assumption 4.1 (Assumptions on the Model). Let f ∈ C 2 (R) with (2.1) and η ∈ C 2 (R) be a fixed entropy. Let ϕ : Z + → Z -be a kinetic relation in the sense of Definition 2.1 that is compatible with η. Let u 0 ∈ L 1 ∩ BV (R) with (2.5).

Assumption 4.2 (Assumptions on the Scheme).

Let M := max(1, L) u 0 L ∞ (R) . Let F ∈ C 1 (R 2 ,
R) be a numerical flux consistent with f , which is monotone and locally Lipschitz-continuous, i.e.,

F (u, u) = f (u), (4.1a)
F (•, u) is increasing, F (u, •) is decreasing, (4.1b)
F (u, v) -F (u ′ , v ′ ) ≤ Ĉ u -u ′ + v -v ′ (4.1c) for all u, u ′ , v, v ′ ∈ [-M, M ]. Here Ĉ is the local Lipschitz constant of F on the domain [-M, M ] 2 . Let Q ∈ C 1 (R 2 ,
R) be a numerical entropy flux function, which is consistent to q, local Lipschitz-continuous, and compatible to F , i.e., Q(u, u) = q(u), (4.2a)

∃ C = C(M ) : F (u, v) -F (u ′ , v ′ ) ≤ C(M ) u -u ′ + v -v ′ , (4.2b)
∂ 2 Q(u, v) = η′ (v)∂ 2 F (u, v) (4.2c) for all u, v, u ′ , v ′ ∈ [-M, M ]. Remark 4.3.
Let us mention that there are consistent, monotone, and locally Lipschitz-continuous numerical fluxes that allow the construction of a corresponding numerical entropy flux with (4.2). Examples are provided by the Lax-Friedrichs and the Enquist-Osher flux [START_REF] Kröner | Numerical schemes for conservation laws[END_REF].

Theorem 4.4 (Convergence and Existence). Let Assumptions 4.1 and 4.2 be fulfilled. Consider for h > 0 the sequence {u h } h>0 of approximate solutions given by Algorithm 3.1. Then we have:

• There exists a subsequence {u hn } n∈N of {u h } h>0 and u ∈ V such that lim n→∞ u hnu L 1 = 0.

• The limit u is a (η, ϕ)-solution of (1.1) in the sense of Definition 2.3.

If we additionally assume that there exists a unique solution, then the statement is also true for the complete sequence.

Proof. The proof is the combination of the Lemmata 5.3, 5.4, and 5.6, Helly's Theorem, Theorem 5.8, and Theorem 5.13 (see below).

∆m(n) = 1 =⇒ α n+1 ≤ 1 4 , (5.1a
)

∆m(n) = -1 =⇒ α n+1 ≥ 3 4 . (5.1b) Furthermore, if u n m(n)+1 < ϕ ♯ (u n m(n)-1 ), then we have for all ξ ∈ (-∞, ϕ ♯ (u n m(n)-1 )) f ′ (ξ) -s n ≥ 0. (5.2)
Proof. From (3.5), (3.6), and ∆m(n) = m(n + 1)m(n) we get

α n+1 = α n + ∆t h s n -∆m(n)    ≤ ∆t h Ĉ ≤ 1 4 for ∆m(n) = 1, ≥ 1 -∆t h Ĉ ≥ 3 4 for ∆m(n) = -1.
Note that α n ∈ [0, 1) for all n ≥ 0. The inequality (5.2) follows directly from (2.1).

In this section we will prove Theorem 4.4 step by step. Note that we will extract several times subsequence of {u h } h>0 . All subsequences are also denoted by {u h } h>0 to avoid an unnecessary notation overload.

Reformulation of the Scheme in Harten

Form. As a first step, we rewrite Algorithm 3.1 in a Harten-like form (see [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF]), which enables us to apply BV compactness theory. We get for j ∈ J(n) ∩ J(n + 1) the formula

u n+1 j = u n j + ∆t h C n j+ 1 2 u n j+1 -u n j - ∆t h D n j-1 2 u n j -u n j-1 (5.3) with C n j+ 1 2 := - F (u n j , u n j+1 ) -F (u n j , u n j ) u n j+1 -u n j and D n j-1 2 := F (u n j , u n j ) -F (u n j-1 , u n j ) u n j -u n j-1
.

Next, we rewrite the cell values near the phase boundary. With (3.3), (3.4), and (4.1a) we get

u n+1 m(n+1)-1 = u n+1 m -(n)-1 = u n m -(n)-1 - ∆t h D n m -(n)-3 2 u n m -(n)-1 -u n m -(n)-2 + ∆m -(n)S n -u n m(n)-1 -u n m(n)-2 ,
(5.4a)

u n+1 m(n+1)+1 = u n+1 m + (n)+1 = u n m + (n)+1 + ∆t h C n m + (n)+ 3 2 u n m + (n)+2 -u n m + (n)+1 - ∆t h D n m(n)+ 1 2 u n m(n)+1 -ϕ(u n m(n)-1 ) -∆m + (n)S n + u n m(n)+2 -u n m(n)+1
(5.4b) with

C n m + (n)+ 3 2 := - 1 2 + ∆m -(n) -α n+1 F (u n m + (n)+1 , u n m + (n)+2 ) -F (u n m + (n)+1 , u n m + (n)+1 ) u n m + (n)+2 -u n m + (n)+1 , D n m(n)+ 1 2 :=      ā" u n m(n)+1 ,ϕ(u n m(n)-1 ) " -s n 2+∆m -(n)-α n+1 for u n m(n)+1 < ϕ ♯ (u n m(n)-1 ), 0 for u n m(n)+1 ≥ ϕ ♯ (u n m(n)-1 ), D n m -(n)-3 2 := 1 1 + ∆m + (n) + α n+1 F (u n m -(n)-1 , u n m -(n)-1 ) -F (u n m -(n)-2 , u n m -(n)-1 ) u n m -(n)-1 -u n m -(n)-2 , S n -:= 1 - 1 1 + α n+1 1 + ∆t h ā(u n m(n)-2 , u n m(n)-1 ) ,
and

S n + := 1 - 1 2 -α n+1 1 - ∆t h ā(u n m(n)+1 , u n m(n)+2 ) .
We use the local Lipschitz-continuity of the numerical flux to bound the generalized Harten coefficients provided we have a uniform L ∞ -bound on the approximate solution. To ease the notation we set

C n m(n)-1 2 = C n m(n)+ 1 2 = D n m(n)-1 2
= 0. Lemma 5.2. Let the assumptions of Theorem 4.4 and

-M ≤ u h (x, t) ≤ M , (x, t) ∈ R × [0, T ]
be satisfied. Then we have

∆t h C n j+ 1 2 , ∆t h D n j+ 1 2 , ∆m -(n)S n -, ∆m + (n)S n + ∈ 0, 1 2
for all h > 0, j ∈ Z, and n ∈ {0, . . . , N }. Proof. By using (3.5), (4.1b) and (4.1c) we get the result for all values except

D n m(n)+ 1 2 . For the exceptional value D n m(n)+ 1 2
we make a distinction of cases: The

case u n m(n)+1 ≥ ϕ ♯ (u n m(n)-1 ) is trivial. In the case u n m(n)+1 < ϕ ♯ (u n m(n)-1 ) we can deduce from (2.1) the positivity of D n m(n)+ 1 2 . We get the inequality ∆t h D n m(n)+ 1 2 ≤ 1 2
, since s n and ā(•, •) are bounded by Ĉ from (4.1c).

5.3.

A Uniform L ∞ -Bound for the Approximate Solution. Lemma 5.3. Let the assumptions of Theorem 4.4 be fulfilled. Then the following estimate holds for all h > 0 and t ∈ [0, T ]:

u h (•, t) L ∞ (R) ≤ M := max(1, L) u 0 L ∞ (R) .
(5.5)

Proof. In view of (3.2) it is obvious that u 0 h L ∞ (R) ≤ u 0 L ∞ (R) ≤ M holds. From (5.3) we get for all j ∈ J(n) ∩ J(n + 1)

u n+1 j ≤ 1 - ∆t h C n j+ 1 2 - ∆t h D n j-1 2 u n j + ∆t h C n j+ 1 2 u n j+1 + ∆t h D n j-1 2 u n j-1
≤ max u n j-1 , u n j , u n j+1 .

(5.6)

With formula (5.4a) and Lemma 5.2 we get also

u n+1 m(n+1)-1 = u n+1 m -(n)-1 ≤ 1 - ∆t h D n m -(n)-3 2 -∆m -(n)S n - u n m -(n)-1 + ∆t h D n m -(n)-3 2 u n m -(n)-2 + ∆m -(n)S n -u n m(n)-1 ≤ max u n m -(n)-1 , u n m -(n)-2 , u n m(n)-1
.

(5.7)

From (5.6) and (5.7) we get by induction

u n j ≤ u 0 h L ∞ (-∞,p 0 ) ≤ u 0 L ∞ (-∞,p 0 ) for all j ≤ m(n) -1 and n ∈ {0, . . . , N }.
We proceed with the cell values right to the phase boundary, i.e., j ≥ m(n) + 1. With the last formula and (2.8a) we immediately get

ϕ(u n m(n)-1 ) ≤ ϕ(0) + L u n m(n)-1 -0 ≤ L u 0 L ∞ (-∞,p 0 ) .
(5.8)

From (5.4b) and in view of Remark 3.3 we get

u n+1 m(n+1)+1 = u n+1 m + (n)+1 ≤ 1 - ∆t h C n m + (n)+ 3 2 -(1 -∆m + (n)) ∆t h D n m(n)+ 1 2 -∆m + (n)S n + u n m + (n)+1 + ∆m + (n) S n + - ∆t h D n m(n)+ 1 2 u n m(n)+1 + ∆t h C n m + (n)+ 3 2 u n m + (n)+2 + ∆t h D n m(n)+ 1 2 ϕ(u n m(n)-1 ) ≤ max u n m + (n)+1 , u n m + (n)+2 , u n m(n)+1 , ϕ(u n m(n)-1 ) .
From this and again (5.6), we obtain by induction

u n j ≤ max(1, L) u 0 L ∞ for all j ≥ m(n) + 1 and n ∈ {0, . . . , N }.
(5.9)

Combining the results (5.8) and (5.9) implies (5.5).

A Uniform Total Variation Bound for the Approximate Solution.

We derive a uniform estimate on the total variation (TV) of u h . Similarly as in Section 5.3 it is possible to estimate the total variation left to the phase boundary only using values left to it.

Lemma 5.4. Let the assumptions of Theorem 4.4 be fulfilled. Then the following formula holds for all N ∈ {1, . . . , N }

j<m(N )-1 u N j+1 -u N j + N -1 n=0 u n+1 m(n+1)-1 -u n m(n)-1 ≤ j<m(0)-1 u 0 j+1 -u 0 j .
Proof. Let us first define w n j+ 1 2 = u n j+1u n j . We get

w n+1 j+ 1 2 ≤ 1 - ∆t h D n j+ 1 2 - ∆t h C n j+ 1 2 w n j+ 1 2 + ∆t h D n j-1 2 w n j-1 2 + ∆t h C n j+ 3 2 w n j+ 3 2 .
Then we sum over j < m(n) -1 and get

j<m(n+1)-1 w n+1 j+ 1 2 = j<m -(n)-1 w n+1 j+ 1 2 ≤ j<m -(n)-1 w n j+ 1 2 - ∆t h D n m -(n)-3 2 w n m -(n)-3 2 - ∆t h C n m -(n)-1 2 =0 w n m -(n)-1 2 = j<m(n)-1 w n j+ 1 2 - ∆t h D n m -(n)-3 2 w n m -(n)-3 2 -∆m -(n)w n m(n)-3 2 ≤ j<m(n)-1 w n j+ 1 2 -u n+1 m(n+1)-1 -u n m(n)-1 .
The last inequality follows directly from (5.4a), since

u n+1 m(n+1)-1 -u n m(n)-1 ≤ ∆t h D n m -(n)-3 2 u n m -(n)-1 -u n m -(n)-2 + ∆m -(n)S n -u n m(n)-1 -u n m(n)-2 ≤ ∆t h D n m -(n)-3 2 w n m -(n)-3 2 + ∆m -(n)w n m(n)-3 2 .
Finally, we sum over 0 ≤ n ≤ N -1 and get the result by rearranging the terms

j<m(N )-1 w N j+ 1 2 + N -1 n=0 u n+1 m(n+1)-1 -u n m(n)-1 ≤ j<m(0)-1 w 0 j+ 1 2 .
Let us now recall the definition of the total variation. For v ∈ L 1 (R) we define

TV(v) := lim h→0 R 1 h v(x + h) -v(x) dx .
A function v is total variation bounded if TV(v) < ∞. For a piecewise constant function as u n h , we get

TV(u n h ) = j / ∈{m(n)-1,m(n)} u n j+1 -u n j + u n m(n)+1 -u n m(n)-1 .
We further introduce the operator TV * as

TV * (u n h ) := j / ∈{m(n)-1,m(n)} u n j+1 -u n j + u n m(n)+1 -ϕ(u n m(n)-1 ) + ϕ(u n m(n)-1 ) -u n m(n)-1 .
We have the trivial relation

TV(u n h ) ≤ TV * (u n h ).
However, the advantage of TV * is that the difference TV * (u n+1 h ) -TV * (u n h ) can be estimated independently of u n m(n)+1 , see proof of Lemma 5.5 below. Note that the maximally possible increase in total variation is not realized within one time step but only in the asymptotic regime. That makes it hard to control TV(u n h ) directly. Lemma 5.5. Let the assumptions from Theorem 4.4 be fulfilled. Then we have for all h > 0 and n ∈ {0, . . . , N }

TV(u n h ) ≤ (1 + 2L) TV * (u 0 ). (5.10)
Here, TV * (u 0 ) is given by

TV * (u 0 ) := TV (u 0 ) -|u + -u -| + u + -ϕ(u -) + ϕ(u -) -u -, u ± = u 0 (0±).
Proof. The following proof exploits again the fact that the scheme can be written in Harten form. For n ≥ 0 we have

TV * (u n+1 h ) -TV * (u n h ) = j∈Z\[m -(n)-1,m + (n)] u n+1 j+1 -u n+1 j -u n j+1 -u n j + u n+1 m(n+1)+1 -ϕ(u n+1 m(n+1)-1 ) -u n m(n)+1 -ϕ(u n m(n)-1 ) -∆m + (n) u n m(n)+2 -u n m(n)+1 + ϕ(u n+1 m(n+1)-1 ) -u n+1 m(n+1)-1 -ϕ(u n m(n)-1 ) -u n m(n)-1 -∆m -(n) u n m(n)-1 -u n m(n)-2 + ∆m + (n) u n+1 m(n+1)-1 -u n+1 m(n+1)-2 + ∆m -(n) u n+1 m(n+1)+2 -u n+1 m(n+1)+1 =: L 1 + L 2 + L 3 + L 4 .
For L 1 we get from (5.3)

L 1 ≤ - ∆t h C n m + (n)+ 3 2 u n m + (n)+2 -u n m + (n)+1 + ∆t h D n m(n)+ 1 2 u n m(n)+1 -ϕ(u n m(n)-1 ) + ∆m + (n)S n + u n m(n)+2 -u n m(n)+1 - ∆t h D n m -(n)-3 2 u n m -(n)-1 -u n m -(n)-2 + ∆m -(n)S n -u n m(n)-1 -u n m(n)-2 =: R 1 + R 2 + R 3 . (5.11)
For the terms L 2 and L 3 we use the following representation of ϕ(u n+1 m(n+1)-1 )

ϕ(u n+1 m(n+1)-1 ) = ϕ(u n m(n)-1 ) -T n ϕ u n+1 m(n+1)-1 -u n m(n)-1
(5.12) with

T n ϕ := ϕ(u n+1 m(n+1)-1 ) -ϕ(u n m(n)-1 ) u n+1 m(n+1)-1 -u n m(n)-1 . From (2.8b) it is clear that T n ϕ ∈ [-L, 0].
To further examine L 2 , we use (5.4b), (5.12), and get

L 2 ≤ -R 1 -R 2 + T n ϕ u n+1 m(n+1)-1 -u n m(n)-1 . (5.13)
With (5.4a) and (5.12) we get

L 3 ≤ -R 3 + T n ϕ u n+1 m(n+1)-1 -u n m(n)-1 . (5.14)
By the definition of the scheme, see (5.4a) and (5.4b), we have L 4 = 0. Finally, we sum all terms from (5.11), (5.13), (5.14) and use Lemma 5.2 to get

TV * (u n+1 h ) -TV * (u n h ) ≤ 2L u n+1 m(n+1)-1 -u n m(n)-1 .
By summing all time steps up to N ≥ 0 we get the final result:

TV * (u N +1 h ) = TV * (u 0 h ) + N n=0 TV * (u n+1 h ) -TV * (u n h ) ≤ TV * (u 0 h ) + 2L N n=0 u n+1 m(n+1)-1 -u n m(n)-1 ≤ (1 + 2L) TV * (u 0 h ).
For the last inequality we used Lemma 5.4. Finally, we get (5.10).

5.5. Continuity in Time. Lemma 5.6 (Continuity in Time). Under the assumptions of Theorem 4.4, there is a C > 0 such that for all

0 ≤ t 1 < t 2 ≤ T u h (•, t 2 ) -u h (•, t 1 ) L 1 (R) ≤ C(t 2 -t 1 ).
Proof. Thanks to the Harten form of the scheme the proof follows exactly the same lines as for monotone schemes, see e.g. [START_REF] Kröner | Numerical schemes for conservation laws[END_REF], Theorem 2.3.9. From Lemmas 5.3, 5.5, and 5.6 and Helly's theorem imply that there is a subsequence of {u h } h>0 and a function u ∈ V such that the first assertion of Theorem 4.4 holds. In the following we will show that the limit function u is a (η, ϕ)-solution of (1.1).

A Cell Entropy Inequality.

In this section we will prove cell entropy inequalities, which will be used in Section 5.7 to prove that u h converges to a function that obeys (2.7).

Lemma 5.7 (Cell Entropy Inequality). Under the assumptions of Theorem 4.4, the entropy inequality

η(u n+1 j ) ≤ η(u n j ) - ∆t h q n j+ 1 2 -q n j-1 2 (5.15)
holds for all j ∈ J(n) ∩ J(n + 1), n ∈ {0, . . . , N }, and h > 0. Furthermore, we have

1 + ∆m + (n) + α n+1 η(u n+1 m(n)-1 ) ≤ (1 + α n ) η(u n m(n)-1 ) + ∆m -(n)η(u n m(n)-2 ) - ∆t h q n m(n) -q n m -(n)- 3 2 (5.16)
and

2 + ∆m -(n) -α n+1 η(u n+1 m(n)+1 ) ≤ (2 -α n ) η(u n m(n)+1 ) + ∆m + (n)η(u n m(n)+2 ) - ∆t h q n m + (n)+ 3 2 -q n m(n) .
(5.17)

Here, the entropy fluxes are defined by

q n j+ 1 2 = Q u n j , u n j+1 and q n m(n) = q(u n m(n)-1 ) -s n η(u n m(n)-1 ).
Proof. The proof for (5.15) can be done as for monotone schemes, see e.g. [START_REF] Kröner | Numerical schemes for conservation laws[END_REF], Proposition 3.3.26.

We show only (5.17) for ∆m(n) = 1, which is the most complicated case. The case ∆m(n) = 1 follows the same lines. Using (3.11b) we get

u n+1 m(n+1)+1 = T 1 + T 2 + (1 -α n+1 ) (T 3 + T 4 ) 2(2 -α n+1 ) (5.18)
with

T 1 := u n m(n)+2 -2 ∆t h f n m(n)+ 5 2 -f (u n m(n)+2 ) , T 2 := u n m(n)+2 -2 ∆t h f (u n m(n)+2 ) -f n m(n)+ 3 2 , T 3 := u n m(n)+1 - 2 1 -α n+1 ∆t h f (u n m(n)+1 ) -s n u n m(n)+1 -h n , T 4 := u n m(n)+1 - 2 1 -α n+1 ∆t h f n m(n)+ 3 2 -f (u n m(n)+1 ) .
As the next step, we state the inequalities

η(T 1 ) ≤ η(u n m(n)+2 ) -2 ∆t h q n m(n)+ 5 2 -q(u n m(n)+2 ) , (5.19a)
η(T 2 ) ≤ η(u n m(n)+2 ) -2 ∆t h q(u n m(n)+2 ) -q n m(n)+ 3 2 , ( 5.19b
)

η(T 3 ) ≤ η(u n m(n)+1 ) - 2 1 -α n+1 ∆t h q(u n m(n)+1 ) -s n η(u n m(n)+1 ) -q n m(n) , (5.19c) η(T 4 ) ≤ η(u n m(n)+1 ) - 2 1 -α n+1 ∆t h q n m(n)+ 3 2 -q(u n m(n)+1 ) . (5.19d)
We give the proof of (5.19c), all others inequalities are similar as in the standard situation [START_REF] Kröner | Numerical schemes for conservation laws[END_REF]. In the Laxian case, where s n = ā u n m(n)+1 , u n m(n)-1 , the inequality is trivially fulfilled, since T 3 = u n m(n)+1 and the wave connecting u n m(n)-1 and u n m(n)+1

fulfils the entropy inequality

s n η(u n m(n)-1 ) -η(u n m(n)+1 ) -q(u n m(n)-1 ) -q(u n m(n)+1 ) ≤ 0.
In the undercompressive case, we set a = u n m(n)+1 , λ = We have to show the inequality g(ϕ(u n m(n)-1 )) ≤ 0, which together with (2.10) implies (5.19c). We will prove g(b) ≤ 0 for the whole interval b ∈ -M, ϕ ♯ (u n m(n)-1 ) . Here M is defined as before by M := max(1, L) u 0 L ∞ .

By differentiating g we get

g ′ (b) = η ′ a -λ f (a) -f (b) -s n (a -b) λ f ′ (b) -s n + λ q ′ (b) -s n η ′ (b) = λ f ′ (b) -s n η ′ a -λ f (a) -f (b) -s n (a -b) -η ′ (b) = λ f ′ (b) -s n η ′′ (ξ 1 ) a -b -λ f (a) -f (b) -s n (a -b) = λ f ′ (b) -s n η ′′ (ξ 1 ) 1 -λ(f ′ (ξ 2 ) -s n ) (a -b),
where

ξ 1 ∈ I a -λ f (a) -f (b) -s n
λ(f ′ (ξ 2 ) -s n ) ≤ 2 8 Ĉ(1 -α n+1 ) f ′ (ξ 2 ) + |s n | ≤ 2 • 2 Ĉ 8 Ĉ 3 4 < 1 we get g ′ (b)(a -b) ≥ 0 for b ∈ -M, ϕ ♯ (u n m(n)-1
) . As we have g(a) = 0, we get g(b) ≤ 0 for b ∈ -M, ϕ ♯ (u n m(n)-1 ) . Putting (5.19a)- (5.19d) together, we get

2 -α n+1 η(u n+1 m(n)+1 ) ≤ 1 2 η(T 1 ) + η(T 2 ) + (1 -α n+1 ) η(T 3 ) + η(T 4 ) = (2 -α n ) η(u n m(n)+1 ) + η(u n m(n)+2 ) - ∆t h q n m(n)+ 5 2 -q n m(n) .
The proof for (5.16) is done in the same way.

5.7. Lax-Wendroff Theorem. In this section we state a Lax-Wendroff theorem for the scheme.

Theorem 5.8 (Lax-Wendroff). Let the assumptions of Theorem 4.4 be fulfilled. If the sequence {u h } h>0 converges to a limit function u ∈ V , then the limit function u is a weak solution in the sense of (2.6) and an entropy solution in the sense of (2.7).

Proof. For the proof we refer to [START_REF] Volker | A Lax-Wendroff type theorem for unstructured quasi-uniform grids[END_REF], which states a result for general quasiuniform grid. We point out that the computational grid introduced in Section 3 fulfils the requirements of a quasi-uniform grid.

5.8. Recovery of the Phase Boundary Position. To prove Theorem 4.4 completely it remains to verify the pointwise properties 2. and 3. in Definition 2.3. In this section we construct the function p : [0, T ] → R from the Definition 2.3. We will verify (2.13)

p h -→ p in C 0 ([0, T ]) for h → 0.
5.9. Trace Values. After the construction of p the trace values u -and u + at the phase boundary are well-defined. The next definition introduces more quantities defined on the phase boundary.

Definition 5.9. Let u ∈ V be the limit of {u h } h>0 and p ∈ Lip([0, T ]) the limit of {p h } h>0 . We define the phase boundary speed s ∈ L 1 ([0, T ]), the phase boundary flux H ∈ L 1 ([0, T ]) and the entropy dissipation D ∈ L 1 ([0, T ]) by the following formulas:

s(t) = ā u + (t), u -(t) , H(t) = f (u -(t)) -s(t)u -(t), D(t) = -s(t) η(u + (t)) -η(u -(t)) + q(u + (t)) -q(u -(t)) .
All quantities defined in Definition 5.9 have discrete analogons, which we define in the next definition.

Definition 5.10. We define the discrete traces values u * h,-, u * h,+ ∈ L 1 ([0, T ]), the discrete phase boundary speed s * h ∈ L 1 ([0, T ]), the discrete phase boundary flux H * h ∈ L 1 ([0, T ]) and the discrete phase boundary dissipation

D * h ∈ L 1 ([0, T ]) by u * h,-(t) = u n m(n)-1 , u * h,+ (t) = u n m(n)+1 , s * h (t) = s n , H * h (t) = h n = f (u n m(n)-1 ) -s n u n m(n)-1 , D * h (t) = h ∆t α n+1 η(u n+1 m(n)-1 ) + (1 -α n+1 )η(u n+1 m(n)+1 ) -α n η(u n m(n)-1 ) -(1 -α n )η(u n m(n)+1 ) + q(u n m(n)+1 ) -q(u n m(n)-1 )
for t ∈ [t n , t n+1 ). The next lemma provides convergence properties of these quantities: Lemma 5.11. Under the assumptions of Theorem 4.4, there are functions u * -, u * + , H * , s * , D * ∈ L 1 ([0, T ]), such that for a subsequence it is true that

u * h,-→ u * -, u * h,+ → u * + , s * h → s * , H * h → H * , and D * h → D * for h → 0 in L 1 ([0, T ]). Furthermore, we have for almost all t ∈ [0, T ] D * (t) = -ā u * -(t), u * + (t) η(u * + (t)) -η(u * -(t)) + q(u * + (t)) -q(u * -(t)) . (5.20)
Proof. Note that u * h,-, u * h,+ , s * h , H * h , and D * h are TV-bound by using Lemma 5.4. Helly's theorem selects a converging subsequence. For the convergence of D * h and the limit representation (5.20) we need (3.5) to control the term h ∆t and the iteration formula (3.6) for the representation.

In the next step we analyze the relationship between the * -quantities in Lemma 5.11 and corresponding trace quantities of the limit function u in Definition 5.9.

Lemma 5.12. Under the assumptions of Theorem 4.4, we have

s(t) = s * (t), H(t) = H * (t), and
D(t) ≤ D * (t) for almost all t ∈ [0, T ]. Proof. s(t) = s * (t): From the facts that s * h = ∂p h ∂t → s * in L 1 ([0, T ]), p h → p in C 0 ([0, T ]
), the existence of s = ∂p ∂t and the uniform bound of s * h L ∞ , we immediately get the result. H(t) = H * (t): We get for almost all 0 ≤ t 1 < t 2 ≤ T the mass balance

p h (t2) -∞ u h (x, t 2 ) dx - p h (t1) -∞ u h (x, t 1 ) dx = t 2 t 1 H * h (t) dt + O(∆t).
In the limit h → 0 we get by We see that D

p h → p in C 0 ([0, T ]) and u h (•, t) → u(•, t) p(t2) -∞ u(x, t 2 ) dx - p(t1) -∞ u(x, t 1 ) dx = t 2

x=p(t)

= Dµ 1 {x=p(t)} . Here µ 1 is the one-dimensional Hausdorff measure and λ the Lebesgue measure. We have to be more careful in the discrete setting to obtain a similar result. We define

D n j := 1 ∆t η(u n j+1 ) -η(u n j ) + 1 h q n j+ 1 2 -q n j-1 2 , D * h := 0≤n<N j / ∈{m(n)-1,m(n),m(n)+1} D n j λ [x j-1 2 ,x j+ 1 2 )×[t n ,t n+1 ) + D * h µ 1 {x=p h (t)} . Again we have D * h x=p h (t) = D * h µ 1 {x=p h (t)} . From D * h ≤ 0 (see Lemma 5.

7) and D *

h → D in the sense of measures the result follows directly.

Let us mention that Lemma 5.12 shows that p fulfils (2.13), which is equivalent to s(t) = s * (t). With the help of the just introduced definitions and lemmas, we can state the main result of this section. We note that the relation s(t) = s * (t) from Lemma 5.12 is equivalent to (2.13) in the definition of the (η, ϕ)-solutions. It remains to verify condition (2.14) In the following we will drop the index h for better reading.

Case 1: There are infinitely many h > 0 such that u n m(n)+1 < ϕ ♯ (u n m(n)-1 ). In this case we can restrict ourselves to a subsequence of {u h } h>0 which fulfils this condition in every step. Recall H * (t) from Definition 5. We measure the L 1 ([-1, 1])-error at time t = 1. The measured errors and resulting EOCs are listed in Table 6 Example 6.2 (Experimental Order of Convergence II). In contrast to Example 6.1 we use here initial data such that the phase boundary is not reproduced exactly. We achieve this by putting two discontinuities into the initial conditions, which evolve as single Laxian shock wave and phase boundary, and collide after some time. After the interaction the phase boundary contributes to the total error. The initial conditions are

u 0 (x) =        3 2 for x < -1 2 , 1
for -1 2 < x < 0, -1 2 for x > 0.

The solution of this problem is illustrated in Figure 6.1. We measure the error in the L 1 ([-1, 1])-norm at time t = 0.45 (see Table 6.1). As in Example 6.1 we observe an EOC of 1.

Fig. 3 . 1 .

 31 Fig. 3.1. Background and calculation grid for four different geometric configurations h 1 to h 4 .

5 .

 5 Proof of Theorem 4.4. 5.1. Preliminary Note. Lemma 5.1 (Properties of α n+1 and s n ). Under the assumptions of Theorem 4.4, we have for n = 0, . . . , N -1 the inequalities

  g(b) := η aλ f (a)f (b)s n (ab)η(a) + λ q(a)q(b)s n (η(a)η(b)) .

  (ab) , b and ξ 2 ∈ I (a, b) are values arising in the mean value theorem. Here, I(a, b) = [a, b] ∪ [b, a]. By Lemma 5.1, (3.5), and

t 1 H

 1 * (t) dt . (5.21) Since u is a weak solution, we get (5.21) also for H instead of H * , what finishes the proof. D(t) ≤ D * (t): To prove this, we introduce a dissipation measure D ∈ M(R × [0, T ]) by D = ∂ t η(u) + ∂ x q(u).

Case 2 : 1 u = -1 2 u = -3 4 u = - 9 8Fig. 6 . 1 .

 2124961 Fig. 6.1. Solution to initial conditions given in Example 6.2. Solid lines indicate phase boundaries and dashed line shock waves.

  in the next section. We define a continuous version p h : [0, T ] → R of the position of the discrete phase boundary (p n ) n≥0 from (3.1) byp h (t) = p n + s n (tt n ) for all t ∈ t n , t n+1 .The set {p h : [0, T ] → R} h>0 is uniformly Lipschitz-continuous, since s n defined

	in (3.3) is bounded, since Lemma 5.3 states that u n m(n)-1 n Arzela-Ascoli theorem we can select a subsequence with	is bounded. By the

  . Theorem 5.13. Under the assumptions of Theorem 4.4, we have for almost every t ∈ [0, T ] that eitheru + (t) = ϕ(u -(t)) or u + (t) ≥ ϕ ♯ (u -(t)) (5.22) is fulfilled. Proof. Fix t ∈ [0, T ]and select a subsequence of {u h } h>0 such that u * h,-(t), u * h,+ (t), s * h (t), and D * h (t) converge. This is possible for almost every t ∈ [0, T ] and we will show (5.22) under these assumptions. Define a sequence {n h } h>0 such that n h = max n ∈ N t n h ≤ t , t n h = n∆t h =

	nh 8 Ĉ .
	For this sequence we have
	u * h,-(t) = u n h m(n

h )-1 → u * -(t), . . . for h → 0.

Table 6 . 1

 61 .1. We observe an approximately linear rate, which we can expect for Riemann initial data, since the undercompressive wave is exactly reproduced in this case. Errors and EOC for the calculations done in Example 6.1 and 6.2 with different numbers of grid cells N .

		Example 6.1	Example 6.2
	N	L 1 -Error EOC	L 1 -Error EOC
	40	0.01567	0.02791
	80	0.00858 0.868	0.01592 0.810
	160	0.00449 0.934	0.00709 1.167
	320	0.00226 0.993	0.00389 0.864
	640	0.00113 1.000	0.00240 0.696
	1280	0.00056 1.000	0.00143 0.754
	2560	0.00028 1.000	0.00072 0.992
	5120	0.00014 1.000	0.00036 0.985

Acknowledgments. The authors P. Engel and C. Rohde would like to thank the German Research Foundation (DFG) for financial support of the project within the Cluster of Excellence in Simulation Technology (EXC 310/1) at the University of Stuttgart.

). The value Ψ(v, w) ∈ R\{v, w} is defined as the solution of ā(v, w) = ā v, Ψ(v, w) .

By the same argument as in Case 1 we get that the limit equation (5.24) has also three solutions given by v, w, Ψ(v, w). The situation u -(t) = v, u + (t) = w fulfils the second part of (5.22), but u -(t) = v, u + (t) = Ψ(v, w) in general not. However, the latter situation does not show up: we prove that for u -(t) = v, u + (t) = Ψ(v, w) we get the contradiction D(t) ≤ D * (t) < D(t).

(5.25)

The first inequality in (5.25) was already shown in Lemma 5.12. To show the second one, we use (5.20

(5.26)

The inequality (5.26) is finally correct, since the shock wave connecting w to Ψ(v, w) is not entropy consistent.

6. Numerical Experiments. In this section we will demonstrate the efficiency of Algorithm 3.1 numerically. We do this by calculating the so-called experimental order of convergence (EOC). For a given sequence of grid widths h i and corresponding numerical solutions u hi we define the EOC by

By u we denote here the (η, ϕ)-solution with η(u) = 1 2 u 2 and ϕ to be fixed. Furthermore, we will present an experiment showing the convergence of the trace values for non-Riemann initial data. The test was proposed in [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF] and we compare our results with the results given there. We give also an example of an (η, ϕ)-solution switching between a Laxian and an undercompressive phase transition, compare to Remark 2.4 (iii).

In the following experiments are done on equi-distant meshes with N grid cells. For the numerical flux we use an upwind numerical flux except for Example 6.4, where a Lax-Friedrichs numerical flux is used.

Example 6.1 (Experimental Order of Convergence I). For Algorithm 3.1 we proved the convergence to an (η, ϕ)-solution. The question of the rate remains open. Therefore we study the experimental order of convergence (EOC) for the scheme. We use the simple flux f (u) = 1 3 u 3 with ϕ(u) = -3 4 u as a kinetic relation. The entropy plays no role in this example. First we study the Riemann problem (1.1), (2.16) with u L = 1 and u R = -1 2 on the domain [-1, 1] for different grids with N cells. The (η, ϕ)-solution is given by

(c) Trace value history from t = 0 to t = 0.8. -

in time for the same calculation. This confirms the statement of Theorem 5.8. In Figure 6.2(c) the trace values

are plotted. We observe convergence to the kinetic relation for N → ∞, as already proved in Theorem 5.13. Calculating the EOC for the traces values gives: 0.1537, 0.5030 and 0.8765. We expect a convergence rate of 0.5. Compared to [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF] we get the same qualitative behavior.

- Example 6.4 (Laxian and Undercompressive Phase Boundaries). In this example we give initial conditions for which the phase boundary switches from a Laxian to an undercompressive phase boundary. We used f (u) = u 3 -1 2 u and ϕ(u) = -3 4 u. The initial datum is given by

For the discretization we used 5000 grid cells. The initial conditions are constructed such that the phase boundary propagates as a Laxian wave to the right. At approximately t = 0.1 a rarefaction wave hits the phase boundary and makes it undercompressive for t > 0.2. The undercompressive phase boundary now travels to the left. In Figure 6.3(c) the evolution of the traces values ϕ ♯ (u n m(n)-1 ) and u n m(n)+1 in time is plotted. The intersection point marks the switch from a Laxian to an undercompressive phase boundary. We refer to Remark 2.4 for the definition and (3.3) for the realization in the scheme.