Nicolas Seguin
email: nicolas.seguin@univ-rennes1.fr

Magali Tournus
email: magali.tournus@centrale-marseille.fr.

Asymptotic preserving discretisation of a Jin-Xin model with implicit equilibrium manifold on a bounded domain

Keywords: Asymptotic Preserving scheme, Hyperbolic Relaxation, Boundary layer. Subject Classifications: 65N08, 65N12, 35L10, 35L65

In this paper, we design and analyze a numerical scheme which approximates a Jin-Xin linear system with implicit equilibrium on a bounded domain. This scheme relaxes toward the asymptotic limit of the linear system. The main properties of the limiting scheme are that it does no require to invert the implicit function defining the manifold, and that it provides an accurate discretization of the boundary conditions.

Introduction

The Jin-Xin model, introduced in [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF], is a 2×2 linear hyperbolic system with a nonlinear dissipative source term which writes    ∂ t s ε + ∂ x w ε = 0,

∂ t w ε + ∂ x s ε = 1 ε (f (s ε ) -w ε ).
When the source term becomes infinitely sharp, i.e. when ε → 0, the conservation law

∂ t ρ + ∂ x f (ρ) = 0.
is obtained, under the subcharacteristic condition |f | 1. Then, the system can be viewed as a dissipative approximation of entropy weak solutions of conservation laws. A large literature is dedicated to this convergence, see for instance [START_REF] Natalini | Convergence to equilibrium for the relaxation approximations of conservation laws[END_REF], [START_REF] Natalini | Convergence of a relaxation approximation to a boundary value problem for conservation laws[END_REF], [START_REF] Aregba-Driollet | Kinetic approximation of a boundary value problem for conservation laws[END_REF], [START_REF] Serre | Relaxations semi-linéaire et cinétique des systèmes de lois de conservation[END_REF]. . . In its usual formulation, the equilibrium manifold is explicit, and is given by {w = f (s)}.

In [START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF][START_REF] Tournus | A model of calcium transport along the rat nephron[END_REF], a model is introduced and analyzed for the evolution of the concentration of chemical species dissolved in a fluid moving along the loop of Henle in the human kidney. It corresponds to a countercurrent exchanger, i.e. a U-shaped circuit, made of two parallel tubes in which a fluid is flowing in opposite directions, connected at one of their ends. In the first tube, fluid moves with positive velocity 1 and has a concentration denoted by u ε (x, t), whereas in the other tube, fluid moves with negative velocity -1 and has a concentration denoted by v ε (x, t). The positive constant ε is the characteristic time associated with the chemical exchanges between the two tubes through the medullary interstitial region. A nonlinear function h encodes the dynamics of the exchange. The governing equations finally are

     ∂ t u ε + ∂ x u ε = 1 ε (h(v ε ) -u ε ), ∂ t v ε -∂ x v ε = 1 ε (u ε -h(v ε )).
(1)

This system also admits the alternative form, by defining

s ε = u ε + v ε and w ε = u ε -v ε ,    ∂ t s ε + ∂ x w ε = 0, ∂ t w ε + ∂ x s ε = 2 ε h(s ε -w ε ) 2 - (s ε + w ε ) 2 . ( 2
)
Formally, this system converges when ε → 0 towards the conservation law

∂ t (h(v) + v) + ∂ x (h(v) -v) = 0. ( 3
)
This implicit equation is well-posed as soon as the flux h(v)-v can be uniquely defined as a function of the unknown h(v)+v, following the Kruzhkov's theory [START_REF] Kružkov | First order quasilinear equations with several independent variables[END_REF]. This will be the case in our study, see assumptions below. However, from the numerical point of view, it is not straightforward to obtain a convergent and conservative numerical scheme for equation (3) starting with a classical scheme for [START_REF] Aregba-Driollet | Kinetic approximation of a boundary value problem for conservation laws[END_REF] and letting ε go to 0, without inverting the flux h(v) -v (with respect to v or h(v) + v). This will be the first goal of our study. As mentioned above, the countercurrent exchanger model is completed by specific boundary conditions. Denoting the domain by [0, L], with L > 0, the initial-boundary value problem (IBVP) writes

               ∂ t u ε + ∂ x u ε = 1 ε (h(v ε ) -u ε ), t > 0, x ∈ [0, L], ∂ t v ε -∂ x v ε = 1 ε (u ε -h(v ε )), t > 0, x ∈ [0, L],
u ε (0, t) = u b , v ε (L, t) = αu ε (L, t), t > 0, (u ε , v ε )(x, 0) = (u 0 , v 0 )(x),

x ∈ [0, L],

(S ε )
where the reflection capacity α is assumed to be in (0, 1), u b ∈ R and initial conditions (u 0 , v 0 ) are of bounded variations

u 0 ∈ BV ([0, L]), v 0 ∈ BV ([0, L]). (4)
Some results of this paper are stated under the additional technical assumption that the initial conditions are at equilibrium

u 0 (x) = h(v 0 (x)), x ∈ [0, L]. (5)
The IBVP for the Jin-Xin model has been studied by several authors, see [START_REF] Aregba-Driollet | Kinetic approximation of a boundary value problem for conservation laws[END_REF], [START_REF] Carbou | Semilinear behavior for totally linearly degenerate hyperbolic systems with relaxation[END_REF], [START_REF] Xin | Stiff well-posedness and asymptotic convergence for a class of linear relaxation systems in a quarter plane[END_REF], [START_REF] Yong | Basic structures of hyperbolic relaxation systems[END_REF]. . . More specifically, the well-posedness and the asymptotic analysis of the IBVP (S ε ) are given in [START_REF] Perthame | A simple derivation of BV bounds for inhomogeneous relaxation systems[END_REF].

In order to understand the IBVP when ε → 0, let us provide the assumptions on function h: there exists two positive constant β µ such that 1 < β h (v) µ, and h(0) = 0.

(

) 6
As a consequence, the function

f : h(v) + v → h(v) -v (7)
is increasing. Following the classical theory of IBVP for conservation laws provided in [START_REF] Bardos | First order quasilinear equations with boundary conditions[END_REF], only the boundary condition at x = 0 persists, and the limit IBVP is thus

     ∂ t (h(v) + v) + ∂ x (h(v) -v) = 0, t > 0, x ∈ [0, L], h(v(0, t)) = u b , t > 0, (h(v) + v)(x, 0) = (h(v 0 ) + v 0 )(x), x ∈ [0, L]. (S 0 )
The convergence of solutions of (S ε ) to solutions of (S 0 ) is provided in [START_REF] Perthame | A simple derivation of BV bounds for inhomogeneous relaxation systems[END_REF]. We complement in the present paper the analysis of [START_REF] Perthame | A simple derivation of BV bounds for inhomogeneous relaxation systems[END_REF] showing the existence of a relaxation boundary layer at x = L if the intersection between the equilibrium manifold

M eq = {(u, v) ∈ R 2 | u = h(v)}
and the boundary manifold

M b = {(u, v) ∈ R 2 | v = αu}
is empty. In this paper, the goal is to obtain and analyze a numerical scheme which fits with the limit (S 0 ) and which is a accurate discretization of the boundary conditions of (S 0 ), using threepoint schemes. At x = 0, the approximation of the boundary condition is classical, using a ghost cell and imposing inside the Dirichlet value. At x = L, in order to avoid any numerical boundary layer, the easiest way is to obtain the upwind scheme when ε → 0, which does not depends on any ghost cell since f > 0. We also show by numerical tests that this numerical treatment also provides an accurate approximation of the relaxation boundary layer on coarse mesh. Let us sum up the requirements we presented on the approximation of the IBVP (S ε ):

1. Provide a first-order approximation of the relaxation system (S ε ) without any nonlinear inversion of the function h.

2. When ε → 0, obtain a first-order approximation of the IBVP (S 0 ),

• without the use of any nonlinear inversion of the function f defined by ( 7),

• with an upwind discretisation of the flux f .

An asymptotic preserving scheme (S ε,∆ ) is usually defined as a convergent scheme for the system (S ε ), which tends to become a convergent scheme (S ∆ ) for the limiting equation as ε goes to zero. In other words, an asymptotic preserving scheme is a scheme (S ε,∆ ) such that the following diagram is commutative.

(S ε ) (S 0 ) (S ε,∆ ) (S ∆ ) 1 ∆ → 0 2 ε → 0 3 ε → 0 4 ∆ → 0 3
In the specific context of time-explicit numerical schemes, a necessary condition is that the CFL condition for (S ε,∆ ) is uniform in ε.

A first idea to build an asymptotic preserving scheme is to use a splitting method [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]. The scheme (S ∆ ) we obtain at the limit is highly diffusive, and generates a numerical boundary layer at x = L. The alternative method we use in the present paper is based on the use of well-balanced schemes, introduced by [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF] and developed in [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF] for the sake of asympptotic preserving scheme. The main idea is to cleverly approximate the source term in order to end up with the wanted discretized version of the flux at the limit. In our context, let us point out that the term ∆t ε + ∆x behaves like

         ∆t ε as ∆x goes to zero, ∆t ∆x as ε goes to zero. ( 8
)
This work fits into the more general problem of building AP schemes with constraint on the limiting scheme (S). After the pionneering work [START_REF] Shi | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF], the authors of [START_REF] Berthon | Asymptotic-preserving Godunov-type numerical schemes for hyperbolic systems with stiff and nonstiff relaxation terms[END_REF] developed a somewhat generic method to make optional the choice of the numerical scheme in the asymptotic regime ε = 0. Properties of stability and convergence are automatically given by the construction provided in [START_REF] Berthon | Asymptotic-preserving Godunov-type numerical schemes for hyperbolic systems with stiff and nonstiff relaxation terms[END_REF], as the scheme they obtain at the limit can be seen as a convex combination of well-known schemes. Our specific problem cannot be directly solved using their method since it provides us with a scheme that requires to invert h. Since the scheme is built by hand and does not correspond to any classical scheme at the limit, we are left with analyzing its basic properties by hand as well.

The oultline of the paper is the following. In Section 2, we provide the definitions of the solutions of the IBVP's (S ε ) and (S 0 ) and the associated well-posedness and asymptotic results. We also describe the relaxation boundary layer which appears as soon as M eq ∩M b = ∅. In Section 3, we design a numerical scheme which fulfills all the above-mentioned requirements, and state the main results of convergence, showing the asymptotic compatibility of the approach (the so-called asymptotic preserving property [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF]). Sections 4 and 5 are dedicated to the proofs of convergence of the scheme, respectively when ε > 0 and when ε = 0. The last section contains numerical results including comparison with the classical splitting method.

2 Well-posedness, zero-relaxation limit, and relaxation boundary layer

In all the following, we assume that assumption (6) is fulfilled, so that function f defined by [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF] exists and is increasing.

Definitions and existing results

Let us provide the definition of weak solutions of the relaxation IBVP (S ε ), regardless of their smoothness.

Definition 1. Consider any initial data (u 0 , v 0 ) satisfying (4), and ε > 0. A weak solution of the 2 satisfying Φ(x, T ) = Ψ(x, T ) = 0 and Ψ(0, t) = 0, the following equality holds

IBVP (S ε ) is a couple of functions (u ε , v ε ) ∈ C (0, T ); L 1 [0, L] ∩ L ∞ ([0, T ]; BV [0, L]) such that for all (Φ, Ψ) ∈ C 1 ([0, T ] × [0, L])
T 0 L 0 u ε ∂ t Φ + v ε ∂ t Ψ + u ε ∂ x Φ -v ε ∂ x Ψ dx dt = - 1 ε T 0 L 0 (h(v ε ) -u ε ) (Φ -Ψ) dx dt - T 0 u b Φ(0, t) dt - T 0 u ε (L, t) αΨ(L, t) -Φ(L, t) dt - L 0 Φ(x, 0)u 0 (x) + Ψ(x, 0)v 0 (x) dx. (9)
The first result is the well-posedness of the relaxation IBVP (S ε ).

Theorem 1 (Well-posedness of the relaxation IBVP [START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF]). Under assumption [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF], there is a unique weak solution to the relaxation IBVP (S ε ), in the sense of Definition 1.

Now, let us define the entropy weak solutions of the zero-relaxation limit, following the theories provided in [START_REF] Kružkov | First order quasilinear equations with several independent variables[END_REF] and [START_REF] Bardos | First order quasilinear equations with boundary conditions[END_REF]. Definition 2. Consider u 0 , v 0 satisfying (4) and [START_REF] Carbou | Semilinear behavior for totally linearly degenerate hyperbolic systems with relaxation[END_REF]

and v b ∈ R. An entropy weak solution to (S 0 ) is a function v ∈ C (0, T ); L 1 [0, L] ∩ L ∞ ([0, T ]; BV [0, L]) such that
1. for all non negative Φ ∈ C 1 ([0, T ) × (0, L)), and for all k ∈ R,

T 0 L 0 |h(v) + v -(h(k) + k)|∂ t Φ + |h(v) -v -(h(k) -k)|∂ x Φ dx dt + L 0 |h(v 0 (x)) + v 0 (x) -(h(k) + k)|Φ(x, 0) dx 0, (10) 2. for all k in the interval I(v(0, t), u b ) sign h(v(0, t)) + v(0, t) -(h(u b ) + u b ) h(v(0, t)) -v(0, t) -(h(k) -k) 0. ( 11
)
Existence and uniqueness of such entropy solution follows from the theory developed in [START_REF] Bardos | First order quasilinear equations with boundary conditions[END_REF]. We state the following result of convergence partially proved in [START_REF] Perthame | A simple derivation of BV bounds for inhomogeneous relaxation systems[END_REF].

Theorem 2 (Convergence [START_REF] Perthame | A simple derivation of BV bounds for inhomogeneous relaxation systems[END_REF]). We assume (4), ( 5) and [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]. Consider a family of solutions (u ε , v ε ) ε>0 to the relaxation IBVP (S ε ). Then there exists a function v which is an entropy weak solution to (S 0 ) in the sense of Definition 2 such that

u ε -→ ε→0 h(v), v ε -→ ε→0 v, L 1 ([0, L] × [0, T ]).
The outline of the proof is as follows. First, a dissipative formulation for (S ε ) is obtained. Combined with L ∞ estimates for u ε and v ε , this proves that (u ε -h(v ε )) goes to zero in L 1 ([0, T ] × [0, L]) (see [START_REF] James | Convergence results for some conservation laws with a reflux boundary condition and a relaxation term arising in chemical engineering[END_REF]). The dissipative formulation also implies that the weak solution (u ε , v ε ) to the linear system (S ε ) satisfies the following entropy formulation: for all non negative Φ ∈ C 1 ([0, T ] × [0, L]) such that Φ(., T ) = 0, and for all k ∈ R,

T 0 L 0 (|u ε -h(k)| + |v ε -k|)∂ t Φ + (|u ε -h(k)| -|v ε -k|)∂ x Φ dx dt + 2 ε T 0 L 0 |u ε -h(v ε )|Φ(x, t) dx dt + L 0 |u 0 (x) -h(k)| + |v 0 (x) -k| Φ(x, 0) dx + T 0 |u b -h(k)| -|v ε (0, t) -k| Φ(0, t)dt - T 0 |u ε (L, t) -h(k)| -|αu ε (L, t) -k| Φ(L, t)dt 0. ( 12
)
The non-linear formulation ( 12) is then passed to the limit: non linear-quantities |u ε -h(k)| and |v ε -k| converges towards |h(v) -h(k)| and |v -k| for some v ∈ L ∞ ([0, T ] × [0, L]) using BV estimates obtained in [START_REF] Perthame | A simple derivation of BV bounds for inhomogeneous relaxation systems[END_REF]. By considering test functions Φ satisfying Φ(0, t) = Φ(L, t) = 0, and letting ε go to zero in [START_REF] Shi | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF], we obtain that v satisfies the first item of Definition 2. Then, we use the same method as in [START_REF] Natalini | Recent results on hyperbolic relaxation problems[END_REF], i.e. we consider for any g ∈ C 1 ([0, T ]) sequences of test functions Φ m such that ∂ x Φ m (x, t) converges toward g(t)δ(x = 0) and we let m go to infinity, which proves that v satisfies the second item of Definition 2.

Study of the relaxation boundary layer

This section is devoted to the existence of the boundary layer in the framework of continuous solutions. We first state that the solution (u ε , v ε ) to (S ε ) is uniformly bounded from above and below.

Proposition 1 (Uniform L ∞ bounds). We assume ( 4) and ( 6). Then, there exists u min , v min , u max and v max which depend on u b , u 0 , v 0 , α, β and µ such that the solution (u ε , v ε ) to (S ε ) satisfies the following estimates

0 < u min u ε (t, x) u max , 0 < v min v ε (t, x) v max , a.e.(x, t) ∈ [0, L] × [0, T ].
We prove here Proposition 1. For U b > 0, we introduce the stationary system

             dU ε dx (x) = 1 ε h(V ε (x)) -U ε (x) , - dV ε dx (x) = 1 ε U ε (x) -h(V ε (x)) , U ε (0) = U b , V ε (L) = αU ε (L). (13)
It was proved in [START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF] that (13) admits a unique solution and that this solution is continuously differentiable and non-negative. The proof of Proposition 1 is based on the following Lemma.

Lemma 1. We assume [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]. Then, there are four scalar numbers u min , u max , v min and v max depending on U b , α, β, µ such that the solution (U ε , V ε ) to (13) satisfies

0 < u min U ε (x) u max , 0 < v min V ε (x) v max , x ∈ [0, L], ε > 0. ( 14
)
Let us prove Lemma 1.

Proof. First

Step. A bound from below for U ε -V ε .

We add the two lines of ( 13) and obtain

d dx (U ε -V ε )(x) = 0 which implies that U ε (x)-V ε (x) does not depend on x. Then, since U ε (0)-V ε (0) = U b -V ε (0) U b and U ε (L)-V ε (L) = (1-α)U ε (L) 0, we have 0 U ε -V ε U b . ( 15
)
We now here that U ε -V ε is uniformly bounded from below by some k min > 0. Let us assume by contradiction that

∀ε 0 > 0, ∀δ > 0, ∃ε < ε 0 such that U ε -V ε < δ. We pick 0 < δ < min β -1 µ U b , (1 -α)U b and ε 0 > 0. Consider ε < ε 0 such that U ε -V ε < δ. Then, for ε < ε 0 , we have (1 -α)U ε (L) = U ε (L) -V ε (L) < δ.
We also have

dU ε dx (0) = 1 ε {h (V ε (0)) -U ε (0)} 1 ε {h (V ε (0)) -V ε (0) -δ} using U ε (0) -V ε (0) < δ 1 ε {(β -1)V ε (0) -δ} , using (6)
1 ε {(β -1)(U b -δ) -δ} 0 since δ < β -1 β U b . The function U ε is continuous, U ε (L) < δ 1 -α U b , U ε (0) = U b and dU ε dx (
0) > 0, then, there exists

x ε ∈ (0, L) such that U ε (x ε ) = max x∈[0,L] {U ε (x)}. Then dU ε dx (x ε ) = 0,
and the first line of ( 13)

implies h (V ε (x ε )) = U ε (x ε ). Then we have U ε (x ε ) -V ε (x ε ) = h (V ε (x ε )) -V ε (x ε ) (β -1)V ε (x ε ) = (β -1)h -1 (U ε (x ε )) β -1 µ U b , which contradicts U ε (x ε )-V ε (x ε ) < δ for δ < β -1 µ U b . Then, by contradiction, there is k min > 0 such that k min < U ε -V ε , ε > 0. ( 16
)
We denote

K ε := U ε -V ε . Second Step. Uniform bounds for U ε and V ε .
Existence of u max . We have U ε (0) = U b , and from [START_REF] Natalini | Recent results on hyperbolic relaxation problems[END_REF] we deduce that

U ε (L) U b 1 -α . If we assume that U ε reaches its maximal value at x ε ∈ (0, L), then, 0 = dU ε dx (x ε ) = 1 ε h(U ε (x ε ) -K ε ) -U ε (x ε ) ,
and thus, using [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF] and

K ε U b , U ε (x ε ) = h(U ε (x ε ) -K ε ) β (U ε (x ε ) -K ε ) βU ε (x ε ) -βU b , which is U ε (x ε ) β β -1 U b .
Then in any case we can set

u max = max 1 1 -α , β β -1 U b .
Existence of u min . Using ( 16) we have

U ε (L) k min 1 -α . If the minimum of U ε is reached at x ε ∈ (0, L), equation (13) gives again directly U ε (x ε ) = h(U ε (x ε ) -K ε ) µU ε (x ε ) -µk min , and U ε (x ε ) µ µ -1 k min . We also recall that U ε (0) = U b .
Then, in any case, we can then set

u min = min k min 1 -α , k min µ µ -1 , U b . Existence of v max . V ε (x) = U ε (x) -K ε u max and we can set v max = u max . Existence of v min . We have V ε (L) = αU ε (L) αu min .
If we assume that V ε reaches its minimum at x = 0, then, V ε is increasing at x = 0. Equation ( 13) then implies h(V ε (0)) -U b 0, and then

V ε (0) U b µ . Now if we assume that V ε reaches its minimum at x ε ∈ (0, L), we have d dx V ε (x ε ) = 0 and then h(V ε (x ε )) = U ε (x ε ) which implies V ε (x ε ) u min µ .
In any case, we can set v min = min u min µ , U b µ , αu min . This ends the proof of Lemma 1.

The comparison principle in [START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF] gives that 0 u 0 (x) U ε (x) and 0 v 0 (x)

V ε (x) implies 0 u ε (t, x) U ε (x) and 0 v ε (t, x) V ε (x). The choice U b = max{u b , u 0 L ∞ , v 0 L ∞ } ends the proof of Proposition 1. The solution u ε , v ε is then uniformly contained in the rectangle [u min , u max ] × [v min , v max ].
As depicted on Figure 2.2, depending on u b , u 0 , v 0 , α, β, µ but not on ε, either M b ∩ M eq = ∅, or there exists (u

I , v I ) ∈ [u min , u max ] × [v min , v max ] such that M b ∩ M eq = {(u I , v I )}.
In the first case, a boundary layer appears. Proposition 2 (Existence of the boundary layer). We assume [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]. Assuming that M b ∩ M eq = ∅, and that the solution (u ε , v ε ) to (S ε ) is continuous with respect to x, then there exists D > 0, independent of ε, and η(ε) > 0 such that

M b ∩ M eq = ∅ M b ∩ M eq = {(u I , v I )}
|h(v ε (x, t)) -u ε (x, t)| > D, t ∈ [0, T ], x ∈ [L -η(ε), L].
Proof. Proposition 1 shows that M b and M eq are the graphs of two continuous functions, respectively v = αu and v = h -1 (u), defined on the compact set [u min , u max ] and which do not intersect. Then there exists m > 0 which does not depend on ε such that

∀(x 0 , y b ) ∈ M b , ∀(x 0 , y eq ) ∈ M eq , |y eq -y b | > m.
Thus for all t ∈ [0, T ], we have

|u ε (L, t) -h(v ε (L, t))| = |u ε (L, t) -h(αu ε (L, t))| β|h -1 (u ε (L, t)) -αu ε (L, t)| βm, since (u ε(L,t) , h -1 (u ε (L, t))) ∈ M eq and (u ε (L, t), αu ε (L, t)) ∈ M b . Since u ε , v ε and

Construction of an Asymptotic Preserving scheme and main results

In the context of the finite volume schemes framework, we consider a mesh of N disjoint cells

C k , k ∈ [[1, N ]].
Let ∆x be the size of each cell and let ∆t be the time step. The final time is denoted by T , and the number of iterations is denoted by n f , so that n f ∆t = T . The approximated value of the function ρ(x, t) for x ∈ C k and t ∈ [(n -1)∆t, n∆t] is denoted by ρ n k .

Construction of the scheme

We detail here the requirements we impose on the numerical schemes. R-1. For simplicity, the scheme (S ε,∆ ) is explicit, and its stencil contains 3 points. R-2. The scheme (S ε,∆ ) is upwind in the sense that the fluxes u and -v are computed using only a one-sided approximation to the derivative. R-3. The scheme (S ∆ ) is upwind in the sense that at each time step, the updated value of the conservative quantity ρ n+1 k only depends on {ρ n , k}.

Based on remark (8), we consider a class of schemes of the form

     u n+1 ε,k = u n ε,k - ∆t ∆x u n ε,k -u n ε,k-1 + ∆t ε + ∆x S u (u n , v n ), v n+1 ε,k = v n ε,k - ∆t ∆x -v n ε,k+1 + v n ε,k - ∆t ε + ∆x S v (u n , v n ), (17)
where S u and S v are two ways to discretize the source term. Both S u and S v should be consistent with h(v) -u. The sum of the equations of ( 17) for ε = 0 is

u n+1 ε,k + v n+1 ε,k = u n ε,k + v n ε,k - ∆t ∆x u n ε,k -u n ε,k-1 + v n ε,k -v n ε,k+1 -S u + S v . ( 18
)
The condition R-3 then leads us to impose that the discretization of the flux is

u n ε,k -u n ε,k-1 + v n ε,k -v n ε,k+1 -S u + S v = h(v n ε,k ) -v n ε,k -(h(v n ε,k-1 ) -v n ε,k-1
). Since we restricted ourselves to linear upwind schemes with a three-point stencil, there exists 3 real numbers a, b, c such that

S u = h(v n ε,k-1 ) -u n ε,k-1 -av n ε,k+1 + bv n ε,k -cv n ε,k-1 , S v = h(v n ε,k ) -u n ε,k + (1 -a)v n ε,k+1 + (b -2)v n ε,k + (1 -c)v n ε,k-1 . (19)
Among the class of numerical schemes [START_REF] Perthame | A simple derivation of BV bounds for inhomogeneous relaxation systems[END_REF] which satisfy [START_REF] Tournus | Analysis of a simplified model of the urine concentration mechanism[END_REF], one can check that the ones which are stable in L ∞ ∩ BV are those where a, b and c satisfy a = 0, b = c, 1 b β. The schemes we select are then written for 1 b β

             u n+1 ε,k -u n ε,k ∆t + u n ε,k -u n ε,k-1 ∆x = 1 ε + ∆x h(v n ε,k-1 ) -u n ε,k-1 + bv n ε,k -bv n ε,k-1 , v n+1 ε,k -v n ε,k ∆t + v n ε,k -v n ε,k+1 ∆x = - 1 ε + ∆x h(v n ε,k ) -u n ε,k + v n ε,k+1 + (b -2)v n ε,k + (1 -b)v n ε,k-1 , u n 0 = u b , v n 0 = h -1 (u b ), v n N +1 = αu n N . (20)
For simplicity, we focus on the case b = 1 and define the scheme for k

∈ [[1, N ]] and n ∈ N              u n+1 ε,k -u n ε,k ∆t + u n ε,k -u n ε,k-1 ∆x = 1 ε + ∆x h(v n ε,k-1 ) -u n ε,k-1 + v n ε,k -v n ε,k-1 , v n+1 ε,k -v n ε,k ∆t + v n ε,k -v n ε,k+1 ∆x = - 1 ε + ∆x h(v n ε,k ) -u n ε,k + v n ε,k+1 -v n ε,k , u n 0 = u b , v n 0 = h -1 (u b ), v n N +1 = αu n N .
(S ε,∆ )

The scheme (S ε,∆ ) satisfy R-1, R-2 and R-3. The sequence of following results states that the AP diagram is commutative.

Convergence results

For all ε > 0, let us define the following functions

u ε,∆ (x, t) = n∈N k∈[1,N ] u n ε,k 1 [n∆t,(n+1)∆t)×C k (x, t), v ε,∆ (x, t) = n∈N k∈[1,N ] v n ε,k 1 [n∆t,(n+1)∆t)×C k (x, t), (21)
where the sequence (u ε,∆ , v ε,∆ ) is given by the scheme (S ε,∆ ).

Theorem 3. Assuming ( 4), [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF], and the CFL condition µ∆t ∆x, the approximate solution (u ε,∆ , v ε,∆ ) defined in [START_REF] Vovelle | Convergence of finite volume monotone schemes for scalar conservation laws on bounded domains[END_REF] satisfies

u ε,∆ -u ε L 1 ((0,T )×[0,L]) → ∆→0 0, v ε,∆ -v ε L 1 ((0,T )×[0,L]) → ∆→0 0.
where (u ε , v ε ) is the unique solution of (S ε ).

The scheme (S ∆ ) is obtained by setting ε = 0 in (S ε,∆ ), and enables us to build a sequence

(u n k , v n k ) , k ∈ [0, N ], n 0. Let us define u ∆ (x, t) = n∈N k∈[1,N ] u n k 1 [n∆t,(n+1)∆t)×C k (x, t), v ∆ (x, t) = n∈N k∈[1,N ] v n k 1 [n∆t,(n+1)∆t)×C k (x, t). (22)
Theorem 4. Assuming (4), ( 5), [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF], and the CFL condition µ∆t ∆x, the approximate solution (u ∆ , v ∆ ) defined in [START_REF] Xin | Stiff well-posedness and asymptotic convergence for a class of linear relaxation systems in a quarter plane[END_REF] satisfies

u ∆ -h(v) L 1 ((0,T )×[0,L]) → ∆→0 0, v ∆ -v L 1 ((0,T )×[0,L]) → ∆→0 0,
where v is the unique solution of (S 0 ).

We prove in the next section that the schemes (S ε,∆ ) are convergent for all ε > 0, and that they relax toward an upwind convergent scheme when ε goes to zero. In Section 2, we stated the results that justify the arrow 2 of the diagram. We focus here on arrows 1 , 3 and 4 .

4 Convergence of the relaxation scheme (S ε,∆ ) as ∆ → 0

This section is devoted to the proof of Theorem 3. To avoid cumbersome notations, we drop the indices ε in the quantities u n ε,k and v n ε,k . Throughout Propositions 3, 4, and 5, we prove uniform estimates on the functions u ε,∆ and v ε,∆ that enables us to pass to the limit using strong compactness.

Proposition 3 (Conservation, monotonicity and positivity). We assume (4), [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]. Then the sequence scheme (S ε,∆ ) satisfies the following properties:

i) The quantity u n k + v n k is preserved, i.e. there exists a numerical flux (G n k+ 1

2

) i,n such that

u n+1 k + v n+1 k = u n k + v n k - ∆t ∆x G n k+ 1 2 -G n k-1 2 , k = k ∈ [[1, N ]], n ∈ N.
ii) Under the Courant-Friedrichs-Levy condition

∆t ∆x µ , (23)
the scheme (S ε,∆ ) is monotone in the sense that we can write

u n+1 k = G(u n k-1 , u n k , v n k-1 , v n k ) v n+1 k = H(u n k , v n k-1 , v n k , v n k+1 ),
where G and H are non-decreasing functions with respect to each of their variables.

iii) The scheme (S ε,∆ ) preserves positivity:

if ∀k ∈ [1, N ], u 0 k 0, v 0 k 0, then ∀n 0, ∀k ∈ [[1, N ]], u n k 0, v n k 0.
Proof. We first write the scheme in a conservative form:

u n+1 k = u n k - ∆t ∆x u n k -u n k-1 + ∆x 2(∆x + ε) (h(v n k ) -u n k + v n k+1 -v n k ) -(h(v n k-1 ) -u n k-1 + v n k -v n k-1 ) + ∆x 2(∆x + ε) h(v n k ) + h(v n k-1 ) -(u n k + u n k-1 ) + v n k+1 -v n k-1 , v n+1 k = v n k - ∆t ∆x v n k -v n k+1 + ∆x 2(∆x + ε) (h(v n k-1 ) -u n k-1 + v n k -v n k-1 ) -(h(v n k ) -u n k + v n k+1 -v n k ) - ∆x 2(∆x + ε) h(v n k ) + h(v n k-1 ) -(u n k + u n k-1 ) + v n k+1 -v n k-1 , (24)
which proves (i). To prove the monotonicity property (ii), let us write the scheme under the form

u n+1 k = 1 - ∆t ∆x u n k + ∆t ∆x - ∆t ∆x + ε u n k-1 + ∆t ∆x + ε h(v n k-1 ) -v n k-1 + ∆t ∆x + ε v n k := G(u n k-1 , u n k , v n k-1 , v n k ), v n+1 k = 1 - ∆t ∆x + ∆t ∆x + ε v n k + ∆t ∆x - ∆t ∆x + ε v n k+1 - ∆t ∆x + ε h(v n k ) + ∆t ∆x + ε u n k := H(u n k , v n k-1 , v n k , v n k+1 ). (25)
For any ∆x > 0, ∆t > 0, it is clear from the assumptions on h that G is non-decreasing with respect to

u n k-1 , v n k-1 , v n k , and that H is non-decreasing with respect to v n k-1 , v n k+1 , v n k .
By the CFL condition [START_REF] Yong | Basic structures of hyperbolic relaxation systems[END_REF] and since µ > 1, we have ∆t < ∆x, which also implies that G is non-decreasing with

u n k and ∂H ∂v n k (u n k , v n k-1 , v n k , v n k+1 ) = ∂ ∂v n k 1 - ∆t ∆x + ∆t ∆x + ε v n k - ∆t ∆x + ε h(v n k ) 1 - ∆t ∆x + ∆t ∆x + ε -µ ∆t ∆x + ε . > 1 -µ ∆t ∆x 0.
In conclusion the scheme is monotone provided that the stability condition ( 23) is satisfied. Finally, we notice that G(0, 0, 0, 0, 0, 0) = H(0, 0, 0, 0, 0, 0) = 0, and the positivity (iii) follows directly from the monotonicity (ii).

We first prove L ∞ estimates that are useful to prove BV estimates.

Proposition 4 (L ∞ estimate). We assume (4), ( 6), the CFL condition [START_REF] Yong | Basic structures of hyperbolic relaxation systems[END_REF] and ∆x < (1 -α)ε. Then there exists a function M such that the solution (u n k , v n k ) to the scheme (S ε,∆ ) satisfies

∀n 0, ∀k ∈ [[1, N ]], u n k M (ε, u 0 L ∞ , v 0 L ∞ ), v n k M (ε, u 0 L ∞ , v 0 L ∞ ).
We start with the following lemma which states the existence of a super-solution.

Lemma 2 (Existence of a super-solution). We assume (6) and we fix ∆x (1 -α)ε. For any δ > 0, there exists U b u b depending on δ, a pair of vectors (U, V ) ∈ R N +1 × R N +2 which may depend on ε, and a function M :

R + × R + → R + such that            U k -U k-1 = ∆x ∆x + ε h(V k-1 ) -U k-1 + V k -V k-1 , k ∈ [[1, N ]], V k -V k+1 = ∆x ∆x + ε -h(V k ) + U k -V k+1 + V k , k ∈ [[0, N ]], U 0 = U b , V N +1 αU N , (26)
and δ U k M (δ, ε), δ V k M (δ, ε), k ∈ [[1, N ]]. (27)
Proof of Lemma 2. We denote by r = ∆x/(∆x + ε). We have 0 < r < 1 -α. To prove Lemma 2, we are decoupling the difficulties. We first prove the existence of solutions for the system

         U k -U k-1 = r h(V k-1 ) -U k-1 + V k -V k-1 , k ∈ [[1, N ]], V k -V k+1 = r -h(V k ) + U k -V k+1 + V k , k ∈ [[0, N ]], U 0 = U b , V N +1 = V b , (28)
where V b ∈ R + and U b ∈ R + are given, using a fixed point argument. Then, we use a shooting method to prove that there is at least one value for V b for which the solution to (28) satisfies V N +1 αU N , which makes it a solution to (26) as well. In the last step, we prove that we can always find U b (δ) the estimates (27).

Step 1. Existence of a solution for (28). We fix U b u b and V b ∈ R + . We build here a solution (U, V ) to (28) by defining U as a fixed point of the following operator Φ : R N +1 → R N +1 (once the vector U is defined, the vector V is directly deduced from the second line of (28)). Given Ū ∈ R N +1 , we define U := Φ[ Ū ] as the unique solution of the system

         U k -U k-1 = r h(V k-1 ) -U k-1 + V k -V k-1 , k ∈ [[1, N ]], V k -V k+1 = r -h(V k ) + Ūk -V k+1 + V k , k ∈ [[0, N ]], U 0 = U b , Ū0 = U b , V N +1 = V b . (29)
We prove here that for any M such that U b h(M ) and

V b M, we have Φ [0, h(M )] N +1 ⊂ [0, h(M )] N +1
. Indeed, let us assume Ū ∈ [0, h(M )] N +1 . We define (V 0 , . . . , V N ) as the unique vector satisfying the second line of (29) and V N +1 = V b . Let us assume that for some k ∈ [[0, N ]], we have 0 V k+1 M (true for k = N ), then the second line of (29) gives us

(rh + (1 -r)Id)V k = r Ūk + (1 -r)V k+1 , k ∈ [[0, N ]],
Since -r 1 -r h(V N ) < 0, and since α -1 + ∆x ε U N 0 for ∆x ε(1 -α), we have P (V b ) < 0 as soon as U 0 -V 0 < 0, i.e. as soon as We prove now that if we choose U b large enough so that

V b > U b exp µ ε , using (36)
α ε L + ε 2 exp - µ ε U b = δ,
then ( 27) is guaranteed. Indeed, using the first line of (26), we obtain by direct induction

U k (1 -r) k U b (1 -r) N +1 U b , k ∈ [[1, N ]], (37)
and using (36) we have

V N +1 αU N +1 , V k exp - µ ε V N +1 , k ∈ [[0, N ]]. (38)
Combining ( 37) and (38), we deduce, since α < 1 and exp -

µ ε < 1, min min k∈[[1,N ]] U k , min k∈[[1,N +1]] V k α exp - µ ε (1 -r) N +1 U b .
Using N ∆x = L, we have

(1 -r) N +1 = 1 - ∆x ∆x + ε N +1 = N ε L + N ε N +1 , and since N → N ε L + N ε N +1 is increasing, we have (1 -r) N +1 ε L + ε 2 ,
and thus min min

k∈[[1,N ]] U k , min k∈[[1,N +1]] V k δ. According to Step 1, for M such that V b M and U b h(M ), we have 0 U k h(M ), 0 V k M, k ∈ [[1, N ]]. ( 39
)
Then, estimate (27) holds for

M (δ, ε) = 1 + L + ε ε 2 exp 2µ ε δ.
This ends the proof of Lemma 2.

Proof of Proposition 4. Consider the approximations (u n k ) k,n and (v n k ) k,n given by (S ε,∆ ). We denote by (U k ) k∈[0,N ] and (V k ) k∈[0,N +1] the vectors of R N given by Lemma 2 and corresponding to δ :

= max{ u 0 L ∞ , v 0 L ∞ }.
Let us assume that for some n 0, we have for all

k ∈ [[1, N ]], u n k U k and v n k V k . This is true for n = 0 since u 0 k = u 0 (k∆x) δ U k and v 0 k = v 0 (k∆x) δ V k . Since the scheme (S ε,∆ ) is monotone, we have u n+1 k = G(u n k-1 , u n k , v n k-1 , v n k ) G(U k-1 , U k , V k-1 , V k ) v n+1 k = H(u n k , v n k-1 , v n k , v n k+1 ) H(U k , V k-1 , V k , V k+1 ), (40)
where G and H are defined in (24). Since U k , V k satisfy (26), we have

G(U k-1 , U k , V k-1 , V k ) U k H(U k , V k-1 , V k , V k+1 ) V k . ( 41
)
The combination of ( 41) and (40) leads to

u n+1 k U k , v n+1 k V k , k ∈ [[1, N ]].
By induction on n, we have then

u n k U k , v n k V k , k ∈ [[1, N ]], n 0,
and thus the result follows from Lemma 2 with M (ε,

u 0 L ∞ , v 0 L ∞ ) = M (max{ u 0 L ∞ , v 0 L ∞ }, ε).
We define

T V (u n ) = N -1 k=0 |u n k+1 -u n k |, T V (v n ) = N -1 k=0 |v n k+1 -v n k |.
Proposition 5 (spatial BV estimate). We assume ( 4) and [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]. Under the CFL condition ( 23) and assuming ∆x < (1 -α)ε, for u n k and v n k given by the scheme (S ε,∆ ), there exists K such that

T V (u n ) + T V (v n ) T V (u 0 ) + T V (v 0 ) + T K(ε, u 0 L ∞ , v 0 L ∞ ).
Proof. We first write

N -1 k=0 |u n+1 k+1 -u n+1 k | + N k=0 |v n+1 k+1 -v n+1 k | = N -1 k=1 u n+1 k+1 -u n+1 k | + N -1 k=1 v n+1 k+1 -v n+1 k | M k + u n+1 1 -u b + v n+1 1 -v 0 M 0
We consider separately the terms M k and M 0 .

Step 1: M k . Using the numerical scheme (S ε,∆ ) for 1 k N -1, we have

|u n+1 k+1 -u n+1 k | 1 - ∆t ∆x (u n k+1 -u n k ) + ∆t ∆x - ∆t ∆x + ε |u n k -u n k-1 | + ∆t ∆x + ε |v n k+1 -v n k | + ∆t ∆x + ε h(v n k ) -h(v n k-1 ) -(v n k -v n k-1 ) |v n+1 k+1 -v n+1 k | 1 - ∆t ∆x - ∆t ∆x + ε (v n k+1 -v n k ) - ∆t ∆x + ε (h(v n k+1 ) -h(v n k )) + ∆t ∆x - ∆t ∆x + ε |v n k+2 -v n k+1 | + ∆t ∆x + ε |u n k+1 -u n k |.
Under the CFL condition [START_REF] Yong | Basic structures of hyperbolic relaxation systems[END_REF] which guarantees the positivity of the coefficients, the terms of M k can be reorganized the following way

M k N -1 k=1 1 - ∆t ∆x + ∆t ∆x + ε |u n k+1 -u n k | + N -2 k=0 ∆t ∆x - ∆t ∆x + ε |u n k+1 -u n k | + N -1 k=1 1 - ∆t ∆x + (2 -µ) ∆t ∆x + ε |v n k+1 -v n k | + N -2 k=0 (µ -1) ∆t ∆x + ε |v n k+1 -v n k | + N k=2 ∆t ∆x - ∆t ∆x + ε |v n k+1 -v n k |.
which is

M k N -2 k=1 |u n k+1 -u n k | + N -2 k=2 |v n k+1 -v n k | + ∆t ∆x - ∆t ∆x + ε |u n 1 -u b | + 1 - ∆t ∆x + ∆t ∆x + ε |u n N -u n N -1 | + 1 - ∆t ∆x + ∆t ∆x + ε |v n 2 -v n 1 | + 1 + (1 -µ) ∆t ∆x + ε |v n N -v n N -1 | + (µ -1) ∆t ∆x + ε |v n 1 -v 0 | + ∆t ∆x - ∆t ∆x + ε |αu n N -v n N |. (42)
Step 2: M 0 . The term corresponding to k = 0 is treated the following way

M 0 = 1 - ∆t ∆x u n 1 + ∆t ∆x - ∆t ∆x + ε u b + ∆t ∆x + ε [h(v 0 ) -v 0 ] + ∆t ∆x + ε v n 1 -u b + 1 - ∆t ∆x - ∆t ∆x + ε v n 1 + ∆t ∆x - ∆t ∆x + ε v n 2 - ∆t ∆x + ε h(v n 1 ) + ∆t ∆x + ε u n 1 -v 0 ,
we rearrange the terms and plug the equality

u b = h(v 0 ) M 0 = 1 - ∆t ∆x |u n 1 -u b | + ∆t ∆x + ε |v n 1 -v 0 | + ∆t ∆x - ∆t ∆x + ε |v n 2 -v n 1 | + ∆t ∆x + ε |u n 1 -u b | + 1 - ∆t ∆x (v n 1 -v 0 ) - ∆t ∆x + ε (h(v n 1 ) -h(v 0 )) .
and thus

M 0 1 - ∆t ∆x + ∆t ∆x + ε |u n 1 -u b | + ∆t ∆x - ∆t ∆x + ε |v n 2 -v n 1 | + 1 - ∆t ∆x + (1 -µ) ∆t ∆x + ε |v n 1 -v 0 |. (43)
We combine now (42)and ( 43) to obtain

N -1 k=0 |u n+1 k+1 -u n+1 k | + |v n+1 k+1 -v n+1 k | N -1 k=0 |u n k+1 -u n k | + |v n k+1 -v n k | + K(ε, u 0 L ∞ , v 0 L ∞ )∆t, (44) where
K(ε, u 0 L ∞ , v 0 L ∞ ) = (α+1)M (ε, u 0 L ∞ , v 0 L ∞ ), since Proposition 4 implies ε ε + ∆x |αu n N - v n N | (α + 1)M (ε, u 0 L ∞ , v 0 L ∞ ).
Proposition 5 follows by immediate induction.

Proof. (of Theorem 3) We use the following lemma.

Lemma 3 (Theorem 2.4. [START_REF] Bressan | Hyperbolic systems of conservation laws[END_REF]). Consider a sequence of functions u ∆ : [0, +∞) × [0, 1] → R with the following properties:

1. TV(u ∆ (t, .)) C, |u ∆ (t, x)| M, ∀(t, x) ∈ [0, T ] × [0, L] 2. [0,L] |u ∆ (t, x) -u ∆ (s, x)|dx L 1 |t -s| + L 2 ∆, ∀s, t 0, then, there exists a function u ∈ L ∞ ([0, T ] × [0, L]) ∩ C((0, T ); L 1 ([0, L])), such that T V (u(t, .)) C and lim ∆→0 u ∆ = u, L 1 ([0, T ] × [0, L]).
In [START_REF] Bressan | Hyperbolic systems of conservation laws[END_REF], Theorem 2.4 is stated for L 2 = 0, but the proof can be easily adapted for L 2 = 0. We now prove Theorem 3. Consider the sequence of functions (u ∆ , v ∆ ) defined in [START_REF] Vovelle | Convergence of finite volume monotone schemes for scalar conservation laws on bounded domains[END_REF]. Using the first item of Lemma 3 is satisfied. Fix t < s. There exists n and m natural numbers such that t ∈ [n∆t, (n + 1)∆t) and s ∈ [m∆t, (m + 1)∆t). Then, using (S ε,∆ ) and Proposition 5, we have

[0,L] |u ∆ (t, x) -u ∆ (s, x)|dx = N k=1 ∆x|u n k -u m k | m-1 =n N k=1 ∆x|u +1 k -u k | ∆t m-1 =n N k=1 |u k -u k-1 | + ∆x ε h(v k-1 ) -u k-1 + ∆x ε v k -v k-1 ∆t m-1 =n T V (u 0 ) + T V (v 0 ) + T K(ε, u 0 L ∞ , v 0 L ∞ ) + 3 + µ ε M (ε, u 0 L ∞ , v 0 L ∞ ) L(|t -s| + ∆t), with L = T V (u 0 ) + T V (v 0 ) + T K(ε, u 0 L ∞ , v 0 L ∞ ) + 3 + µ ε M (ε, u 0 L ∞ , v 0 L ∞ ).
Thus, the second item of Lemma 3 is satisfied as well, which guarantees the existence of u ε and

v ε in L ∞ ([0, T ] × [0, L]) ∩ C((0, T ); L 1 ([0, L])), such that T V (u ε (t, .
)) and T V (v ε (t, .)) are finite, and such that lim

∆→0 u ε,∆ = u ε , lim ∆→0 v ε,∆ = v ε , L 1 ([0, T ] × [0, L]).
We do not detail the fact that the limit (u ε , v ε ) is a weak solution to the linear system (S ε ), but similar (and easier) arguments developed in next section when proving that (u, v) is the unique entropy solution to (S 0 ) can be applied.

Convergence of the equilibrium scheme (S ∆ )

In this section, we focus on the equilibrium scheme. The goal is to prove Theorem 4. It is allowed to allocate the value 0 to the parameter ε in the scheme (S ε,∆ ). The scheme obtained for ε = 0 gives us two sequences u n k and v n k . We prove here that the sequences obtained are a good discretization of the solution to the system (S 0 ). We denote by λ = ∆t ∆x .

The scheme we obtain at the limit can be written, for k ∈ [[1, N ]] and n 0, as

     v n+1 k = v n k -λ h(v n k ) -u n k , u n+1 k = u n k -λ u n k -v n k -(h(v n k-1 ) -v n k-1 ) , h(v 0 k ) = u 0 k , k ∈ [1, N ], u n 0 = u b , h(v n 0 ) = u b , (S ∆ )
or equivalently, using the intermediate variable

s n k := u n k + v n k , as s n+1 k = s n k -λ h(v n k ) -v n k -(h(v n k-1 ) -v n k-1 ) , (45)
v n+1 k = v n k -λ(h(v n k ) -u n k ), (46)
u n+1 k = s n+1 k -v n+1 k , (47)
h(v 0 k ) = u 0 k , k ∈ [1, N ], u n 0 = u b , h(v n 0 ) = u n 0 , n 0. ( 48
)
The scheme (S ∆ ) provides a solver for the system S which does not require to invert the nonlinear function h at each time step. The scheme can be interpreted the following way. Let us assume that the quantities u n k , v n k are given for a time step n and recall that s n k = u n k + v n k . The variable s n+1 k is updated using (45), which discretizes the scalar equation

∂ t (u + v) + ∂ x (h(v) -v) = 0. Then, we aim to define v n+1 k := h -1 (u n k ) without inverting h.
This is approximately done with step (46), that may be seen as an approximation of the first iteration of the Newton iteration scheme solving h(x) = u n k with the inital guess x = v n k . Note that this would be exactly the first iteration of the Newton iteration scheme in the case where λ = (h (v n k )) -1 . Let us now prove that the solution to the scheme (S ∆ ) converges toward the entropy solution of S. To do so, we first prove monotonicity and a priori BV bounds on u ∆ , v ∆ to guarantee the convergence of the sequence toward a couple (u, v) using Helly theorem. To make sure that the limit is the entropy solution to S, we then adapt the proof of the Lax-Wendroff theorem to our case.

The property is then true for

n = 1 since u 1 k = U λ (u 0 k , v 0 k-1 , v 0 k ), v 1 k = V λ (u 0 k , v 0 k )
, and is easily generalized by induction. To obtain the BV bounds we notice that we have directly from the definition of the total variation and using a triangle inequality and (23) that

T V (u n+1 ) (1 -λ) T V (u n ) + λµ T V (v n ), T V (v n+1 ) λ T V (u n ) + (1 -λµ) T V (v n ).
Summing the two lines, we get

T V (u n+1 ) + T V (v n+1 ) T V (u n ) + T V (v n ),
which proves Lemma 5.

Lemma 6 (Discrete entropy inequalities). We assume (4), [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF] and the CFL condition [START_REF] Yong | Basic structures of hyperbolic relaxation systems[END_REF]. Then there exists C > 0 such that the numerical approximations u n k and v n k obtained by the scheme (S ∆ ) satisfiy the following discrete entropy inequality for n 0 and k

∈ [[1, N ]] |u n+1 k + v n+1 k -(h(κ) + κ)| -|u n k + v n k -(h(κ) + κ)| + λ G n k+1/2 -G n k-1/2 0 ( 58
)
where

G k+1/2 = h(v n k κ) -v n k κ -(h(v n k ⊥ κ) -v n k ⊥ κ) (59)
Proof. We recall that

a b = max(a, b), a ⊥ b = min(a, b), |a -b| = a b -a ⊥ b.
We have

(u n+1 k + v n+1 k ) (h(κ) + κ) u n+1 k h(κ) + v n+1 k κ
and the monotonicity of U λ and V λ implies that

u n+1 k h(κ) + v n+1 k κ U λ (u n k h(κ), v n k-1 κ, v n k κ) + V λ (u n k h(κ), v n k κ). (60)
Rearranging the terms of (60) leads to

(u n+1 k + v n+1 k ) (h(κ) + κ) u n k h(κ) + v n k κ -λ{h(v n k κ) -v n k κ -(h(v n k-1 κ) -v n k-1 κ)} (61)
and for the same reason we have

(u n+1 k + v n+1 k ) ⊥ (h(κ) + κ) u n k ⊥ h(κ) + v n k ⊥ κ -λ{h(v n k ⊥ κ) -v n k ⊥ κ -(h(v n k-1 ⊥ κ) -v n k-1 ⊥ κ)} (62)
The subtraction of (62) to (61) gives

|u n+1 k + v n+1 k -(h(κ) + κ)| -|u n k + v n k -(h(κ) + κ)| + λ G n k+1/2 -G n k-1/2 0,
The result of Lemma 6 thus holds.

with L s = (µ -1) T V (u 0 ) + T V (v 0 ) . We also have using (45) and Lemma 7

[0,L] |v ∆ (t, x) -v ∆ (τ, x)|dx = N k=1 ∆x|v n k -v m k | m-1 =n N k=1 ∆x|v +1 k -v k | m-1 =n ∆t N k=1 h(v n ) -u k m-1 =n ∆t N k=1 1 µ 2 1 -C T V (u 0 ) + T V (v 0 ) L v (|t -τ | + ∆t), with L v = 1 µ 2 1 -C T V (u 0 ) + T V (v 0
) . Thus, since u ∆ = s ∆ -v ∆ , we have the same type of estimate for u ∆ and Lemma 3 guarantees the existence of u and v belonging to

L ∞ ([0, T ] × [0, L]) ∩ C((0, T ); L 1 ([0, L]
)), such that T V (u(t, .)) C, T V (v(t, .)) C and such that lim

∆→0 u ∆ = u, lim ∆→0 v ∆ = v, L 1 ([0, T ] × [0, L]).
We prove now that the limit is an entropy solution to (S 0 ) and that u = h(v) almost everywhere. To do so, we assume first that the problem is posed on the whole space line, i.e. that x ∈ R. We begin by rewriting the discrete entropy inequality we proved on Lemma 6

∆x |u n+1 k + v n+1 k -(h(κ) + κ)| -|u n k + v n k -(h(κ) + κ)| + ∆t G n k+1/2 -G n k-1/2 0. (65) Let us consider a nonnegative test function ϕ ∈ C ∞ c ([0, T ] × R), with T > 0. We introduce ϕ n k = 1 ∆x x i+1/2
x i-1/2 ϕ(n∆t, x) dx.

We multiply (65) by ϕ n k and sum over n ∈ N and k ∈ [1, N ]. Therefore, if we define

A ∆t = ∆x n∈N k∈Z |u n+1 k + v n+1 k -(h(κ) + κ)| -|u n k + v n k -(h(κ) + κ)| ϕ n k , (66)
B ∆t = ∆t n∈N k∈Z G n k+1/2 -G n k-1/2 ϕ n k , (67)
the discrete inequality entropy becomes

A ∆t + B ∆t 0. ( 68
)
The goal is to pass inequality (68) to the limit as ∆t goes to zero and to check that v satisfies the classical continuous entropy inequality associated with (S 0 ). Let us begin with proving that

A ∆t → A 0 =: - R + R |h(v)+v-(h(κ)+κ)|∂ t ϕ dx dt- R |u 0 +v 0 -(h(κ)+κ)|ϕ(0, x) dx when ∆t → 0.
First, we split A ∆t into two parts:

A ∆t = Ā∆t + Ã∆t , with Ā∆t = ∆x n∈N k∈Z |h(v n+1 k ) + v n+1 k -(h(κ) + κ)| -|h(v n k ) + v n k -(h(κ) + κ)| ϕ n k , Ã∆t = ∆x n∈N k∈Z R(u n+1 k , v n+1 k , κ) -R(u n k , v n k , κ) ϕ n k , with R(u, v, κ) = |u + v -(h(κ) + κ)| -|h(v) + v -(h(κ) + κ)|.
The convergence of Ā∆t is classical. Indeed, using the Abel rule (discrete integration by part),

Ā∆t = -∆x k∈Z |h(v 0 k ) + v 0 k -(h(κ) + κ)|ϕ 0 k + n∈N |h(v n+1 k ) + v n+1 k -(h(κ) + κ)|(ϕ n+1 k -ϕ n k ) = - R |h(v ∆t (0, x)) + v ∆t (0, x) -(h(κ) + κ)|ϕ(0, x) dx - R R + |h(v(t + ∆t, x)) + v(t + ∆t, x) -(h(κ) + κ)| ϕ(t + ∆t, x) -ϕ(t, x) ∆t dt dx
Since the two terms converge strongly, we obtain Ā∆t → -

R + R |h(v) + v -(h(κ) + κ)|∂ t ϕ dx dt - R |h(v 0 ) + v 0 -(h(κ) + κ)|ϕ(0, x) dx
when ∆t → 0. We obtain for Ã∆t by similar calculations

Ã∆t = -∆x k∈Z R(u 0 k , v 0 k , κ)ϕ 0 k -∆x k∈Z n∈N R(u n+1 k , v n+1 k , κ)(ϕ n+1 k -ϕ n k ) which is equal to - R R(u (0, x), v ∆t (0, x), κ)ϕ(0, x)dx - R R + R(u ∆t (t + ∆t), v ∆t (t + ∆t), κ) ϕ(t + ∆t, x) -ϕ(t, x) ∆t dt dx. It is straightforward that - R R(u ∆t (0, x), v ∆t (0, x), κ)ϕ(0, x)dx → - R |u 0 +v 0 -(h(κ)+κ)|-|h(v 0 )+v 0 -(h(κ)+κ)| ϕ(0, x) dx
while for the last term, we have

R R + R(u ∆t (t + ∆t), v ∆t (t + ∆t), κ) ϕ(t + ∆t, x) -ϕ(t, x) ∆t dt dx ∂ t ϕ L ∞ (R + ×R) (µ + 1)T sup t 0 u ∆t -h(v ∆t ) L 1 (R) (t),
which tends to 0 when ∆t → 0 using Lemma 7. Gathering all these results gives us A ∆t → A 0 .

Let us now focus on B ∆t . One may first remark that

G n k+1/2 = h(v n k κ) -v n k κ -(h(v n k ⊥ κ) -v n k ⊥ κ)
since h is increasing. Therefore, we have

B ∆t = -∆t n∈N k∈Z G k-1/2 (ϕ n k -ϕ n k-1 )
which is equal to

- R + R (h(v ∆t (n∆t, x) κ) -v ∆t (n∆t, x) κ-(h(v ∆t (n∆t, x) ⊥ κ) -v ∆t (n∆t, x) ⊥ κ) ϕ(n∆t, x) -ϕ(n∆t, x -∆x) ∆x dx dt
and tends toward -

R + R (h(v ∆t κ) -v ∆t κ -(h(v ∆t ⊥ κ) -v ∆t ⊥ κ)(t, x) ∂ x ϕ(t, x) dx dt.
As this stage, we have proved that any limit of the numerical scheme (S ∆ ) satisfies the entropy inequalities [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF], restricting the support of ϕ to [0, T ) × (0, L). Concerning the second point of Definition 2, one can use the Otto's formalism and invoke a more general result of convergence, see for instance [START_REF] Vovelle | Convergence of finite volume monotone schemes for scalar conservation laws on bounded domains[END_REF].

This ends the proof to Theorem 4.

Numerical illustrations

We now present some numerical results which illustrate the good behavior of our numerical scheme. We only focus on the case of a linear source term h(v) = µv, with µ = 3, and the boundary conditions are u ε (0, t) = 1 and v ε (L, t) = αu ε (L, t), with α = 0.1 and L = 1. We choose the same initial data for all experiments:

∀x ∈ (0, L), u ε (x, 0) = u ε (x, 0) = 1.
Let us note that this initial data is neither compatible with the left boundary condition, nor with the equilibrium u = h(v). As a result, we expect to see a right-going wave initiated by the left boundary condition and also a boundary layer at the right boundary.

We compare our numerical scheme (called the AP scheme in the sequel) with the classical splitting method, using an implicit Euler method for the source term:

u n+1/2 k = u n k -λ(u n k -u n k-1 ) v n+1/2 k = v n k + λ(v n k+1 -v n k ) u n+1 k = u n+1/2 k + ∆t ε (h(v n+1 k ) -u n+1 k ) v n+1 k = v n+1/2 k + ∆t ε (u n+1 k -h(v n+1 k )) (69)
where the second part can be explicitly solved since the source term is linear. In all the numerical tests, we used ∆t = ∆x/3.

Figure 2: Comparison of the AP scheme and of the splitting method with a reference solution for several mesh sizes: 50 (up), 200 (center), 800 (bottom) -u + v vs. space.

Different mesh sizes for the relaxation model

The first test illustrates the accuracy of both schemes according to the number of cells: we use successively 50, 200 and 800 cells, while ε = 10 -2 and T = 1.

In Figure 2, one can check that both schemes seem to converge towards the same profile, the reference solution, which has been computed using the AP scheme with 3000 cells. The splitting method is clearly more diffusive than the AP scheme, in particular when the number of cells is small. The wave and the boundary layer are better approximated by the AP scheme.

Behavior of the relaxation boundary layer

In this test, we study how the numerical schemes approximate the boundary layer at the right boundary. To do so, we use a final time T equal to 5 which corresponds to a stationary solution. Three values of ε are used: 10 -1 , 10 -2 and 10 -5 . The number of cells is 100 for all the tests, so that when ε = 10 -5 , one may expect an under-resolved boundary layer.

In Figure 3, we have plotted the results provided by the schemes, and a reference solution computed by the AP scheme with 1000 cells. Let us mention that we represent only the right part of the domain in order to better see the differences. For ε = 10 -1 , the profiles provided by the two numerical schemes are similar but one can note that the point at the right boundary given by the AP scheme is significantly greater than the points obtained by the splitting method and the reference solution. When ε is equal to 10 -2 , this difference increases. However, the shape of the boundary layer is much better approximated by the AP scheme than by the splitting method. The case of ε = 10 -5 leads to much larger discrepancies. The boundary layer of the reference solution is so tiny that it cannot be seen, so that in the figure, one can only see a constant state. The AP scheme provides the same constant state, this is due to its upwind nature when ε is very small. On the contrary, the splitting method leads to a large numerical boundary layer, which would remain even for ε = 0. The only way to make it disappear would be to let the number of cells tend to infinity.

Numerical results for the equilibrium case

We now investigate the behavior of both schemes when ε = 0, using different numbers of cells. Results are plotted at time T = 1. Up and at the center of Figure 4, the results of the splitting method and of the AP scheme are shown, for the unknown u + v. One can check that no boundary layer is present with the AP scheme. The splitting method, which still suffer from a dependence of the right boundary condition, provides a numerical boundary layer at the right, which reduces when the number of cells increases. Moreover, since the AP scheme is an upwind scheme when ε = 0, it is less diffusive and more accurate that the splitting method.

The figure at the bottom represents |h(v) -u| in order to understand of the results are far from the equilibrium. Since the AP scheme is fully explicit and do not use the inverse of function h, one cannot expect to be exactly at the equilibrium. The gap from the equilibrium appears near the discontinuity and disappears when the number of cells increases.

Figure 1 :

 1 Figure 1: Plot of the two equilibrium manifolds. Either M eq (blue) intersects M b (red) inside [u min , u max ] × [v min , v max ](right), either it does not (left).

  h are continuous, this implies that there exists η that may depend on ε such that |h(v ε -u ε )| > βm/2 for |x -L| < η, and Proposition 2 holds for D = βm/2.

  . As a conclusion, for any U b u b and V b satisfying V b > U b exp µ ε , any solution to (28) is a solution to (26). Step 4. Estimates from below and from above for the super-solution. We fix U b u b and V b = U b exp µ ε + 1 and denote by (U, V ) the solution to (26) we selected in Step 1.

Figure 3 :Figure 4 :

 34 Figure3: Comparison of the AP scheme and of the splitting method with a reference solution for several values of ε: 10 -1 (up), 10 -2 (center), 10 -5 (bottom) -u + v vs. space. 29

since (rh + (1 -r)Id) is invertible and monotone. By induction, estimate (30) holds for any k ∈ [[0, N ]]. We now define U by the first line of (29) and U 0 = U b . For k = 0, we have the estimate 0 U k h(M ). The first line of (29) gives for any k ∈ [[0, N ]]

By induction, this directly implies, for all k ∈ [[0, N -1]], 0 U k+1 h(M ).

We can now apply the Brouwer fixed point theorem to conclude that Φ admits at least one fixed point U . We notice that for any of these fixed points U , the couple (U, V ) where V is defined by the second line of (28), is a solution to (28), which guarantees the existence of a solution to (28). This solution satisfies 0

for any M such as U b h(M ) and V b M , and, moreover,

From now on, for each U b u b , V b ∈ R + , we choose a solution to (29) obtained through the process described in Step 1 and denote it by (U, V ).

Step 2. An intermediate estimate on the solution. We prove here an estimate on the solution to (28) defined in Step 1. We have

By direct induction,

we conclude combining (34) and (35) that

Step 3. The shooting method. We fix U b u b . We define the shooting function P as P : V b → αU N -V b , where U N is defined as the N -th component of the vector U selected in Step 1. We prove here that there exists V b ∈ R + such that P (V b ) 0, i.e. such that the (selected) solution to (28) is a solution to (26) as well. We have

Using (33) and (28), we obtain

Let us define the functions

so that the numerical scheme (S ∆ ) can be rewritten

Lemma 4 (Monotonicity of the numerical scheme). We assume ( 4) and [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]. Under the CFL condition (23), the numerical scheme (S ∆ ) is monotone, in the sense that 1. It preserves constant solutions at equilibrium: for all v ∈ R + ,

2. The functions U λ and V λ are nondecreasing with respect to each variable.

Proof. The first item of Lemma 4 is clear. To prove the second item, we rewrite U λ and V λ as

Clearly λ > 0, (h -Id)by assumption and the CFL condition (23) can be written as 1 -λ > 0. This implies that U λ is increasing. Besides, assuming [START_REF] Yong | Basic structures of hyperbolic relaxation systems[END_REF], one has (Id-λh) = 1-λh 1-λµ > 0, then w(Id -λh) is increasing. This implies that V λ is increasing.

Lemma 5 (A priori bounds). We assume (4), [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF] and the CFL condition [START_REF] Yong | Basic structures of hyperbolic relaxation systems[END_REF]. Then the following estimates are satisfied for the sequences (u n k ) and (v n k ) defined by (S ∆ ).

1. L ∞ bounds. For n 0 and k

2. BV bounds. For n 0,

Proof. The L ∞ bounds follow directly from the monotonicity of the scheme. Indeed, assuming

We define for all vector (u k ) k∈Z

Lemma 7 (Control of the deviation with respect to equilibrium). We assume (4), ( 6) and the CFL condition [START_REF] Yong | Basic structures of hyperbolic relaxation systems[END_REF]. Then there is C < 1 such that the discrepancy to equilibrium for the numerical approximations u n k and v n k obtained by the scheme (S ∆ ) is controlled in the L 1 norm by

In particular, since λ 1 µ and µ > 1, we have -1 < C < 1.

Proof. The Rolle's theorem gives us the existence of

Combinig ( 64) with (S ∆ ) gives us

We use a triangle inequality, multiply by ∆x and sum from 1 to N to get

which implies the estimate (63) by induction on n.

Proposition 6 (Convergence of the numerical scheme). We assume (4), ( 5), ( 6) and the CFL condition [START_REF] Yong | Basic structures of hyperbolic relaxation systems[END_REF], then the sequence of numerical approximations

loc when ∆t → 0, up to a subsequence. Proof. of Proposition 6. We apply the same method used for Theorem 3. The sequences (u ∆ ) ∆ and (v ∆ ) ∆ satisfy the estimate 1 required for applying Lemma 3 (see Lemma 5) and we have, using (45) and Lemma 5,