Ivan Nourdin
email: ivan.nourdin@upmc.fr

Giovanni Peccati
email: giovanni.peccati@gmail.com

Anthony Réveillac

Multivariate normal approximation using Stein's method and Malliavin calculus

Keywords: Breuer-Major CLT, fractional Brownian motion, Gaussian processes, Malliavin calculus, Normal approximation, Stein's method, Wasserstein distance théorème de la limite centrale de Breuer et Major, mouvement brownien fractionnaire, processus gaussiens, calcul de Malliavin, approximation normale, méthode de Stein, distance de Wasserstein

, Nualart and Ortiz-Latorre (2007), Peccati (2007) and Nourdin and Peccati (2007b, 2008); in particular, they apply to approximations by means of Gaussian vectors with an arbitrary, positive definite covariance matrix. Among several examples, we provide an application to a functional version of the Breuer-Major CLT for fields subordinated to a fractional Brownian motion.

Introduction

Let Z ∼ N (0, 1) be a standard Gaussian random variable on some probability space (Ω, F, P ), and let F be a real-valued functional of an infinite-dimensional Gaussian field. In the papers [START_REF] Nourdin | Stein's method on Wiener chaos[END_REF][START_REF] Nourdin | Stein's method and exact Berry-Esséen asymptotics for functionals of Gaussian fields[END_REF] it is shown that one can combine Stein's method (see e.g. [START_REF] Chen | Stein's method for normal approximation[END_REF], [START_REF] Ch | A bound for the error in the normal approximation to the distribution of a sum of dependent random variables[END_REF] or [START_REF] Ch | Approximate computation of expectations[END_REF]) with Malliavin calculus (see e.g. [START_REF] Nualart | The Malliavin calculus and related topics (2 nd Edition)[END_REF]), in order to deduce explicit (and, sometimes, optimal) bounds for quantities of the type d(F, Z), where d stands for some distance between the law of F and the law of Z (e.g., d can be the Kolmogorov or the Wasserstein distance). The aim of this paper is to extend the results of [START_REF] Nourdin | Stein's method on Wiener chaos[END_REF][START_REF] Nourdin | Stein's method and exact Berry-Esséen asymptotics for functionals of Gaussian fields[END_REF] to the framework of the multidimensional Gaussian approximation in the Wasserstein distance. Once again, our techniques hinge upon the use of infinite-dimensional operators on Gaussian spaces (like the divergence operator or the Ornstein-Uhlenbeck generator) and upon an appropriate multidimensional version of Stein's method (in a form close to Chatterjee and Meckes [START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF], but see also Reinert and Röllin [20]). As a result, we will obtain explicit bounds, both in terms of Malliavin derivatives and contraction operators, thus providing a substantial refinement of the main findings by Nualart and Ortiz-Latorre [START_REF] Nualart | Central limit theorem for multiple stochastic integrals and Malliavin calculus[END_REF] and Peccati and Tudor [START_REF] Peccati | Gaussian limits for vector-valued multiple stochastic integrals[END_REF]. Note that an important part of our computations (see e.g. Lemma 3.7) are directly inspired by those contained in [START_REF] Nualart | Central limit theorem for multiple stochastic integrals and Malliavin calculus[END_REF]: we shall indeed stress that this last reference contains a fundamental methodological breakthrough, showing that one can deal with (possibly multidimensional) weak convergence on a Gaussian space, by means of Malliavin-type operators and "characterizing" differential equations. See [START_REF] Nourdin | Non-central convergence of multiple integrals[END_REF] for an application of these techniques to non-central limit theorems. Incidentally, observe that the paper [START_REF] Peccati | Gaussian approximations of multiple integrals[END_REF], which is mainly based on martingale-type techniques, also uses distances between probability measures (such as the Prokhorov distance) to deal with multidimensional Gaussian approximations on Wiener space, but without giving explicit bounds.

The rationale behind Stein's method is better understood in dimension one. In this framework, the starting point is the following crucial result, proved e.g. in [START_REF] Ch | A bound for the error in the normal approximation to the distribution of a sum of dependent random variables[END_REF].

Lemma 1.1 (Stein's Lemma).

A random variable Y is such that Y Law = Z ∼ N (0, 1) if and only if, for every continuous and piecewise continuously differentiable function f : R → R such that E f ′ (Z) < ∞, one has

E[f ′ (Y ) -Y f (Y )] = 0.
(1.1)

The fact that a random variable Y satisfying (1.1) is necessarily Gaussian can be proved by several routes: for instance, by taking f to be a complex exponential, one can show that the characteristic function of Y , say ψ(t), is necessarily a solution to the differential equation ψ ′ (t) + tψ(t) = 0, and therefore ψ(t) = exp(-t 2 /2); alternatively, one can set f (x) = x n , n = 1, 2, ..., and observe that (1.1) implies that, for every n, one must have E(Y n ) = E(Z n ), where Z ∼ N (0, 1) (note that the law of Z is determined by its moments).

Heuristically, Lemma 1.1 suggests that the distance d(Y, Z), between the law of a random variable Y and that of Z ∼ N (0, 1), must be "small" whenever

E[f ′ (Y ) -Y f (Y )] ≃ 0,
for a sufficiently large class of functions f . In the seminal works [START_REF] Ch | A bound for the error in the normal approximation to the distribution of a sum of dependent random variables[END_REF][START_REF] Ch | Approximate computation of expectations[END_REF], Stein proved that this somewhat imprecise argument can be made rigorous by means of the use of differential equations. To see this, for a given function g : R → R, define the Stein equation associated with g as

g(x) -E[g(Z)] = h ′ (x) -xh(x), ∀x ∈ R, (1.2)
(we recall that Z ∼ N (0, 1)). A solution to (1.2) is a function h which is Lebesgue-almost everywhere differentiable, and such that there exists a version of h ′ satisfying (1.2) for every x ∈ R. If one assumes that g ∈ Lip(1) (that is, if g Lip ≤ 1, where • Lip stands for the usual Lipschitz seminorm), then a standard result (see e.g. [START_REF] Ch | Approximate computation of expectations[END_REF]) yields that (1.2) admits a solution h such that h ′ ∞ ≤ 1 and h ′′ ∞ ≤ 2. Now recall that the Wasserstein distance between the laws of two real-valued random variables Y and X is defined as

d W (Y, X) = sup g∈Lip(1) |E[g(Y )] -E[g(X)]| ,
and introduce the notation

F W = {f : f ′ ∞ ≤ 1, f ′′ ∞ ≤ 2}
. By taking expectations on the two sides of (1.2), one obtains finally that, for Z ∼ N (0, 1) and for a generic random variable Y ,

d W (Y, Z) ≤ sup f ∈F W E[f ′ (Y ) -Y f (Y )] , (1.3)
thus giving a precise meaning to the heuristic argument sketched above (note that an analogous conclusion can be obtained for other distances, such as the total variation distance or the Kolmogorov distance -see e.g. [START_REF] Chen | Stein's method for normal approximation[END_REF] for a discussion of this point). We stress that the topology induced by d W , on probability measures on R, is stronger than the topology induced by weak convergence.

The starting point of [START_REF] Nourdin | Stein's method on Wiener chaos[END_REF][START_REF] Nourdin | Stein's method and exact Berry-Esséen asymptotics for functionals of Gaussian fields[END_REF] is that a relation such as (1.3) can be very effectively combined with Malliavin calculus, whenever Y is a centered regular functional of some infinite dimensional Gaussian field. To see this, denote by DY the Malliavin derivative of Y (observe that DY is a random element with values in some adequate Hilbert space H), and write L to indicate the (infinite-dimensional) Ornstein-Uhlenbeck generator (see Section 2 below for precise definitions). One crucial relation proved in [START_REF] Nourdin | Stein's method on Wiener chaos[END_REF], and then further exploited in [START_REF] Nourdin | Stein's method and exact Berry-Esséen asymptotics for functionals of Gaussian fields[END_REF], is the upper bound

d W (Y, Z) ≤ E|1 -DY, -DL -1 Y H |. (1.4)
As shown in [START_REF] Nourdin | Stein's method on Wiener chaos[END_REF], when specialized to the case of Y being equal to a multiple Wiener-Itô integral, relation (1.4) yields bounds that are intimately related with the CLTs proved in [START_REF] Nualart | Central limit theorem for multiple stochastic integrals and Malliavin calculus[END_REF] and [START_REF] Nualart | Central limit theorems for sequences of multiple stochastic integrals[END_REF]. See [START_REF] Nourdin | Stein's method and exact Berry-Esséen asymptotics for functionals of Gaussian fields[END_REF] for a characterization of the optimality of these bounds; see again [START_REF] Nourdin | Stein's method on Wiener chaos[END_REF] for extensions to non-Gaussian approximations and for applications to the Breuer-Major CLT (stated and proved in [START_REF] Breuer | Central limit theorems for nonlinear functionals of Gaussian fields[END_REF]) for functionals of a fractional Brownian motion.

The principal contribution of the present paper (see e.g. the statement of Theorem 3.5 below) consists in showing that a relation similar to (1.4) continues to hold when Z is replaced by a d-dimensional (d ≥ 2) Gaussian vector F = (F 1 , ..., F d ) of smooth functionals of a Gaussian field, and d W is the Wasserstein distance between probability laws on R d (see Definition 3.1 below). Our results apply to Gaussian approximations by means of Gaussian vectors with arbitrary positive definite covariance matrices. The proofs rely on a multidimensional version of the Stein equation (1.2), that we combine with standard integration by parts formulae on an infinite-dimensional Gaussian space. Our approach bears some connections with the paper by Hsu [START_REF] Hsu | Characterization of Brownian motion on manifolds through integration by parts[END_REF], where the author proves an hybrid Stein/semimartingale characterization of Brownian motions on manifolds, via Malliavintype operators.

The paper is organized as follows. In Section 2 we provide some preliminaries on Malliavin calculus. Section 3 contains our main results, concerning Gaussian approximations by means of vectors of Gaussian random variables with positive definite covariance matrices. Finally, Section 4 deals with two applications: (i) to a functional version of the Breuer-Major CLT (see [START_REF] Breuer | Central limit theorems for nonlinear functionals of Gaussian fields[END_REF]), and (ii) to Gaussian approximations of functionals of finite normal vectors, providing a generalization of a technical result proved by Chatterjee in [START_REF] Chatterjee | Fluctuation of eigenvalues and second order Poincaré inequalities[END_REF].

Preliminaries and notation

In this section, we recall some basic elements of Malliavin calculus for Gaussian processes. The reader is referred to [START_REF] Nualart | The Malliavin calculus and related topics (2 nd Edition)[END_REF] for a complete discussion of this subject. Let X = {X(h), h ∈ H} be an isonormal Gaussian process on a probability space (Ω, F, P ). This means that X is a centered Gaussian family indexed by the elements of an Hilbert space H, such that, for every pair h, g ∈ H one has that E[X(h)X(g)] = h, g H .

We let L 2 (X) be shorthand for the space L 2 (Ω, σ(X), P ). It is well known that every random variable F ∈ L 2 (X) admits the chaotic expansion F = E(F ) + ∞ n=1 I n (f n ) where the deterministic kernels f n , n ≥ 1, belong to H ⊙n and the convergence of the series holds in L 2 (X). One sometimes uses the notation I 0 (f 0 ) = E[F ]. In the particular case where H := L 2 (T, A, µ), with (T, A) a measurable space and µ is a σ-finite measure without atoms, the random variable I n (f n ) coincides with the multiple Wiener-Itô integral (of order n) of f n with respect to X (see [START_REF] Nualart | The Malliavin calculus and related topics (2 nd Edition)[END_REF]Section 1.1.2.]).

Let f ∈ H ⊙p , g ∈ H ⊙q and 0 ≤ r ≤ p ∧ q. We define the rth contraction f ⊗ r g of f and g as the element of H ⊗(p+q-2r) given by

f ⊗ r g := ∞ i 1 ,...,ir=1
f, e i 1 ⊗ . . . ⊗ e ir H ⊗r ⊗ g, e i 1 ⊗ . . . ⊗ e ir H ⊗r , where {e k , k ≥ 1} is a complete orthonormal system in H. Note that f ⊗ 0 g = f ⊗ g; also, if p = q, then f ⊗ p g = f, g H ⊗p . Note that, in general, f ⊗ r g is not a symmetric element of H ⊗(p+q-2r) ; the canonical symmetrization of f ⊗ r g is denoted by f ⊗ r g. We recall the product formula for multiple stochastic integrals:

I p (f )I q (g) = p∧q r=0 r! p r p q I p+q-2r (f ⊗ r g).
Now, let S be the set of cylindrical functionals F of the form

F = ϕ(X(h 1 ), . . . , X(h n )), (2.1)
where n ≥ 1, h i ∈ H and the function ϕ ∈ C ∞ (R n ) is such that its partial derivatives have polynomial growth. The Malliavin derivative DF of a functional F of the form (2.1) is the square integrable H-valued random variable defined as

DF = n i=1 ∂ i ϕ(X(h 1 ), . . . , X(h n ))h i , (2.2)
where ∂ i ϕ denotes the ith partial derivative of ϕ. In particular, one has that DX(h) = h for every h in H. By iteration, one can define the mth derivative D m F of F ∈ S , which is an element of L 2 (Ω; H ⊙m ), for m ≥ 2. As usual D m,2 denotes the closure of S with respect to the norm • m,2 defined by the relation

F 2 m,2 = E[F 2 ] + m i=1 E[ D i F 2 H ⊗i ].
Note that every finite sum of Wiener-Itô integrals always belongs to D m,2 (∀m ≥ 1). The Malliavin derivative D satisfies the following chain rule formula:

if ϕ : R n → R is in C 1 b
(defined as the set of continuously differentiable functions with bounded partial derivatives) and if (F 1 , . . . , F n ) is a random vector such that each component belongs to D 1,2 , then ϕ(F 1 , . . . , F n ) is itself an element of D 1,2 , and moreover

Dϕ(F 1 , . . . , F n ) = n i=1 ∂ i ϕ(F 1 , . . . , F n )DF i .
(2.

3)

The divergence operator δ is defined as the dual operator of D. Precisely, a random element u of L 2 (Ω; H) belongs to the domain of δ (denoted by Domδ) if there exists a constant

c u satisfying |E[ DF, u H ]| ≤ c u F L 2 (Ω)
for every F ∈ S ; in this case, the divergence of u, written δ(u), is defined by the following duality property:

E[F δ(u)] = E[ DF, u H ], ∀F ∈ D 1,2 .
(2.4)

The crucial relation (2.4) is customarily called the (Malliavin) integration by parts formula.

In what follows, we shall denote by T = {T t : t ≥ 0} the Ornstein-Uhlenbeck semigroup. We recall that, for every t ≥ 0 and every F ∈ L 2 (X),

T t (F ) = ∞ n=0 e -nt J n (F ), (2.5)
where, for every n ≥ 0 and for the rest of the paper, the symbol J n denotes the projection operator onto the nth Wiener chaos, that is onto the closed linear subspace of L 2 (X) generated by the random variables of the form H n (X(h)) with h ∈ H such that h H = 1, and H n the nth Hermite polynomial defined by (4.1). Note that T is indeed the semigroup associated with an infinite-dimensional stationary Gaussian process with values in R H , having the law of X as an invariant distribution (see e.g. [14, Section 1.4] for a more detailed discussion of the Ornstein-Uhlenbeck semigroup in the context of Malliavin calculus; see Barbour [START_REF] Barbour | Stein's method for diffusion approximations[END_REF] for a version of Stein's method involving Ornstein-Uhlenbeck semigroups on infinite-dimensional spaces; see Götze [START_REF] Götze | On the rate of convergence in the multivariate CLT[END_REF] for a version of Stein's method based on multi-dimensional Ornstein-Uhlenbeck semigroups). The infinitesimal generator of the Ornstein-Uhlenbeck semigroup is noted L. A square integrable random variable F is in the domain of L (noted DomL) if F belongs to the domain of δD (that is, if F is in D 1,2 and DF ∈ Domδ) and, in this case, LF = -δDF. One can prove that LF is such that LF = -∞ n=0 nJ n (F ). As an example, if F = I q (f q ), with f q ∈ H ⊙q , then LF = -qF . Note that, for every F ∈ DomL, one has E(LF ) = 0. The inverse L -1 of the operator L acts on zero-mean random variables

F ∈ L 2 (X) as L -1 F = -∞ n=1 1 n J n (F ).
In particular, for every q ≥ 1 and every F = I q (f q ) with f q ∈ H ⊙q , one has that L -1 F = -1 q F . We conclude this section by recalling two important characterizations of the Ornstein-Uhlenbeck semigroup and its generator. i) Mehler's formula. Let F be an element of L 2 (X), so that F can be represented as an application from R H into R. Then, an alternative representation (due to Mehler) of the action of the Ornstein-Uhlenbeck semigroup T (as defined in (2.5)) on F , is the following:

T t (F ) = E[F (e -t a + 1 -e -2t X)] | a=X , t ≥ 0, (2.6)
where a designs a generic element of R H . See Nualart [14, Section 1.4.1] for more details on this and other characterizations of T .

ii) Differential characterization of L. Let F ∈ L 2 (X) have the form F = f (X(h 1 ), ..., X(h d )), where f ∈ C 2 (R d ) has bounded first and second derivatives, and h i ∈ H, i = 1, ..., d. Then,

LF = d i,j=1 ∂ 2 f ∂x i ∂x j (X(h 1 ), ..., X(h d )) h i , h j H - d i=1 ∂f ∂x i (X(h 1 ), ..., X(h d ))X(h i ). (2.7)
See Propositions 1.4.4 and 1.4.5 in [START_REF] Nualart | The Malliavin calculus and related topics (2 nd Edition)[END_REF] for a proof and some generalizations of (2.7).

Stein's method and Gaussian vectors

We start by giving a definition of the Wasserstein distance, as well as by introducing some useful norms over classes of real-valued matrices.

Definition 3.1. (i) The Wasserstein distance between the laws of two R d -valued random vectors X and Y , noted d W (X, Y ), is given by

d W (X, Y ) := sup g∈H ; g Lip ≤1 E[g(X)] -E[g(Y )] ,
where H indicates the class of Lipschitz functions, that is, the collection of all func-

tions g : R d → R such that g Lip := sup x =y |g(x) -g(y)| x -y R d < ∞ (with • R d the usual Euclidian norm on R d ).
A op := sup x R d =1 Ax R d . Remark 3.2.
1. For every d ≥ 1 the topology induced by d W , on the class of all probability measures on R d , is strictly stronger than the topology induced by weak convergence (see e.g. Dudley [START_REF] Dudley | Real Analysis and Probability (2 nd Edition)[END_REF]Chapter 11]).

2. The reason why we focus on the Wasserstein distance is nested in the statement of the forthcoming Lemma 3.3. Indeed, according to relation (3.4), in order to control the second derivatives of the solution of the Stein equation (3.3) associated with g, one must use the fact that g is Lipschitz.

3. According to the notation introduced in Definition 3.1(ii), relation (2.7) can be rewritten as

LF = C, Hessf (Z) H.S. -Z, ∇f (Z) R d , (3.1)
where Z = (X(h 1 ), ..., X(h d )), and C = {C(i, j) : i, j = 1, ..., d} is the

d × d covariance matrix such that C(i, j) = E(X(h i )X(h j )) = h i , h j H .
Given a d × d positive definite symmetric matrix C, we use the notation N d (0, C) to indicate the law of a d-dimensional Gaussian vector with zero mean and covariance C. The following result, which is basically known (see e.g. [START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF] or [START_REF] Reinert | Multivariate normal approximation with Stein's method of exchangeable pairs under a general linearity condition[END_REF]), is the d-dimensional counterpart of Stein's Lemma 1. 1. In what follows, we provide a new proof which is almost exclusively based on the use of Malliavin operators. for every twice differentiable function f :

R d → R such that E| C, Hessf (Y ) H.S. | + E| Y, ∇f (Y ) R d | < ∞, it holds that E[ Y, ∇f (Y ) R d -C, Hessf (Y ) H.S. ] = 0. (3.2) (ii) Consider a Gaussian random vector Z ∼ N d (0, C). Let g : R d → R belong to C 2 (R d )
with first and second bounded derivatives. Then, the function U 0 (g) defined by

U 0 g(x) := 1 0 1 2t E[g( √ tx + √ 1 -tZ) -g(Z)]dt
is a solution to the following differential equation (with unknown function f ):

g(x) -E[g(Z)] = x, ∇f (x) R d -C, Hessf (x) H.S. , x ∈ R d . (3.3)
Moreover, one has that

sup x∈R d Hess U 0 g(x) H.S. ≤ C -1 op C 1/2 op g Lip . (3.4)
Remark 3.4.

1. If C = σ 2 I d for some σ > 0 (that is, if Z is composed of i.i.d. centered Gaussian random variables with common variance equal to σ 2 ), then Proof of Lemma 3.3. We start by proving Point (ii). First observe that, without loss of generality, we can suppose that Z = (Z 1 , ..., Z d ) := (X(h 1 ), ...X(h d )), where X is an isonormal Gaussian process over H = R d , the kernels h i belong to H (i = 1, ..., d), and

C -1 op C 1/2 op = σ -2 I d op σ 2 I d 1/2 op = σ -1 .
h i , h j H = E(X(h i )X(h j )) = E(Z i Z j ) = C(i, j)
. By using the change of variable 2u = -log t, one can rewrite U 0 g(x) as follows

U 0 g(x) = ∞ 0 {E[g(e -u x + 1 -e -2u Z)] -E[g(Z)]}du. Now define g(Z) := g(Z) -E[g(Z)
], and observe that g(Z) is by assumption a centered element of L 2 (X). For q ≥ 0, denote by J q ( g(Z)) the projection of g(Z) on the qth Wiener chaos, so that J 0 ( g(Z)) = 0. According to Mehler's formula (2.6),

E[g(e -u x + 1 -e -2u Z)]| x=Z -E[g(Z)] = E[ g(e -u x + 1 -e -2u Z)]| x=Z = T u g(Z),
where x denotes a generic element of R d . In view of (2.5), it follows that

U 0 g(Z) = ∞ 0 T u g(Z)du = ∞ 0 q≥1 e -qu J q ( g(Z))du = q≥1 1 q J q ( g(Z)) = -L -1 g(Z).
Since g belongs to C 2 (R d ) with bounded first and second derivatives, it is easily seen that the same holds for U 0 g. By exploiting the differential representation (3.1), one deduces that

Z, ∇U 0 g(Z) R d -C, HessU 0 g(Z) H.S. = -LU 0 g(Z) = LL -1 g(Z) = g(Z) -E[g(Z)].
Since the matrix C is positive definite, we infer that the support of the law of Z coincides with R d , and therefore (e.g. by a continuity argument) we obtain that

x, ∇U 0 g(x) R d -C, Hess U 0 g(x) H.S. = g(x) -E[g(Z)],
for every x ∈ R d . This yields that the function U 0 g solves the Stein's equation (3.3).

To prove the estimate (3.4), we first recall that there exists a unique non-singular symmetric matrix A such that A 2 = C, and that one has that On the other hand, by noting h A -1 (x) = h(A -1 x), one obtains by standard computations (recall that A is symmetric) that Hess U

A -1 Z ∼ N d (0, I d ). Now write U 0 g(x) = h(A -1 x), where h(x) = 1 0 1 2t E[g A ( √ tx + √ 1 -tA -1 Z) -g A (A -1 Z)]dt,
0 g(x) = Hess h A -1 (x) = A -1 Hess h(A -1 x)A -1 , yielding sup x∈R d
Hess U 0 g(x) H.S. = sup

x∈R d A -1 Hess h(A -1 x)A -1 H.S.
= sup

x∈R d A -1 Hess h(x)A -1 H.S. ≤ A -1 2 op sup x∈R d Hess h(x) H.S. (3.6) ≤ A -1 2 op A op g Lip (3.7) ≤ C -1 op C 1/2 op g Lip . (3.8)
The chain of inequalities appearing in formulae (3.6)-(3.8) are mainly a consequence of the usual properties of the Hilbert-Schmidt and operator norms. Indeed, to prove inequality (3.6) we used the relations

A -1 Hess h(x)A -1 H.S. ≤ A -1 op Hess h(x)A -1 H.S. ≤ A -1 op Hess h(x) H.S. A -1 op ; relation (3.7
) is a consequence of (3.5); finally, to show the inequality (3.8), one uses the fact that

A -1 op ≤ A -1 A -1 op = C -1 op and A op ≤ AA op = C op .
We are now left with the proof of Point (i) in the statement. The fact that a vector Y ∼ N d (0, C) necessarily verifies (3.2) can be proved by standard integration by parts.

On the other hand, suppose that Y verifies (3.2). Then, according to Point (ii), for every g ∈ C 2 (R d ) with bounded first and second derivatives,

E(g(Y )) -E(g(Z)) = E( Y, ∇U 0 g(Y ) R d -C, Hess U 0 g(Y ) H.S. ) = 0,
where Z ∼ N d (0, C). Since the collection of all such functions g generates the Borel σ-field on R d , this implies that Y Law = Z, thus yielding the desired conclusion.

The following statement is the main result of this paper. Its proof makes a crucial use of the integration by parts formula (2.4) discussed in Section 2.

d W (F, Z) ≤ C -1 op C 1/2 op E C -Φ(DF ) 2 H.S (3.9) = C -1 op C 1/2 op d i,j=1 E[(C(i, j) -DF i , -DL -1 F j H ) 2 ], (3.10)
where we write Φ(DF ) to indicate the matrix Φ(DF ) := { DF i , -DL -1 F j H : 1 ≤ i, j ≤ d}.

Proof. We start by proving that, for every g ∈ C 2 (R d ) with bounded first and second derivatives,

|E[g(F )] -E[g(Z)]| ≤ C -1 op C 1/2 op g Lip E C -Φ(DF ) 2 H.S .
To prove such a claim, observe that, according to Point (ii

) in Lemma 3.3, E[g(F )] - E[g(Z)] = E[ F, ∇U 0 g(F ) R d -C, HessU 0 g(F ) H.S. ]. Moreover, E[ C, HessU 0 g(F ) H.S. -F, ∇U 0 g(F ) R d ] = E   d i,j=1 C(i, j)∂ 2 ij U 0 g(F ) - d i=1 F i ∂ i U 0 g(F )   = d i,j=1 E C(i, j)∂ 2 ij U 0 g(F ) - d i=1 E LL -1 F i ∂ i U 0 g(F ) (since E(F i ) = 0) = d i,j=1 E C(i, j)∂ 2 ij U 0 g(F ) + d i=1 E δ(DL -1 F i )∂ i U 0 g(F ) (since δD = -L) = d i,j=1 E C(i, j)∂ 2 ij U 0 g(F ) - d i=1 E D(∂ i U 0 g(F )), -DL -1 F i H (by (2.4)) = d i,j=1 E C(i, j)∂ 2 ij U 0 g(F ) - d i,j=1 E ∂ 2 ji U 0 g(F ) DF j , -DL -1 F i H (by (2.3)) = d i,j=1 E ∂ 2 ij U 0 g(F ) C(i, j) -DF i , -DL -1 F j H = E Hess U 0 g(F ), C -Φ(DF ) H.S. ≤ E Hess U 0 g(F ) 2 H.S E C -Φ(DF ) 2
H.S (by the Cauchy-Schwarz inequality)

≤ C -1 op C 1/2 op g Lip E C -Φ(DF ) 2 H.S (by (3.4)).
To prove the Wasserstein estimate (3.9), it is sufficient to observe that, for every globally Lipschitz function g such that g Lip ≤ 1, there exists a family {g ε : ε > 0} such that:

(i) for each ε > 0, the first and second derivatives of g ε are bounded;

(ii) for each ε > 0, one has that g ε Lip ≤ g Lip ;

(iii) as ε → 0, g ε -g ∞ ↓ 0.

For instance, we can choose

g ε (x) = E g(x + √ εN ) with N ∼ N d (0, I d ).
Observe that Theorem 3.5 generalizes relation (1.4) (that was proved in [12, Theorem 3.1]). We now aim at applying Theorem 3.5 to vectors of multiple stochastic integrals. Corollary 3.6. Fix d ≥ 2 and 1 ≤ q 1 ≤ . . . ≤ q d . Consider a vector

F := (F 1 , . . . , F d ) = (I q 1 (f 1 ), . . . , I q d (f d )) with f i ∈ H ⊙q i for any i = 1 . . . , d. Let Z ∼ N d (0, C), with C positive definite. Then, d W (F, Z) ≤ C -1 op C 1/2 op 1≤i,j≤d E C(i, j)- 1 q j DF i , DF j H 2 . ( 3

.11)

Proof. We have -L -1 F j = 1 q j F j so that the desired conclusion follows from (3.10).

When one applies Corollary 3.6 in concrete situations (see e.g. Section 4 below), one can use the following result in order to evaluate the right-hand side of (3.11). Lemma 3.7. Let F = I p (f ) and G = I q (g), with f ∈ H ⊙p and g ∈ H ⊙q (p, q ≥ 1). Let a be a real constant. If p = q, one has the estimate:

E a - 1 p DF, DG H 2 ≤ (a -p! f, g H ⊗p ) 2 + p 2 2 p-1 r=1 (r -1)! 2 p -1 r -1 4 (2p -2r)! f ⊗ p-r f 2 H ⊗2r + g ⊗ p-r g 2 H ⊗2r .
On the other hand, if p < q, one has that

E a - 1 q DF, DG H 2 ≤ a 2 + p! 2 q -1 p -1 2 (q -p)! f 2 H ⊗p g ⊗ q-p g H ⊗2p + p 2 2 p-1 r=1 (r -1)! 2 p -1 r -1 2 q -1 r -1 2 (p + q -2r)! f ⊗ p-r f 2 H ⊗2r + g ⊗ q-r g 2 H ⊗2r . Remark 3.8. 1. Recall that E I p (f )I q (g) = p! f, g H ⊗p if p = q, 0 otherwise.
2. In order to estimate the right-hand side of (3.11), we see that it is sufficient to asses the quantity f i ⊗ r f i H ⊗2(q i -r) for any i ∈ {1, . . . , d} and r ∈ {1, . . . , q i -1} on the one hand, and f i , f j H ⊗q i for any 1 ≤ i, j ≤ d such that q i = q j on the other hand.

Proof of Lemma 3.7 (see also [START_REF] Nualart | Central limit theorem for multiple stochastic integrals and Malliavin calculus[END_REF]Lemma 2]). Without loss of generality, we can assume that H = L 2 (A, A , µ), where (A, A ) is a measurable space, and µ is a σ-finite and non-atomic measure. Thus, we can write DF, DG H = p q I p-1 (f ), I q-1 (g) H = p q

A I p-1 f (•, t) I q-1 g(•, t) µ(dt) = p q A p∧q-1 r=0 r! p -1 r q -1 r I p+q-2-2r f (•, t) ⊗ r g(•, t) µ(dt) = p q p∧q-1 r=0 r! p -1 r q -1 r I p+q-2-2r (f ⊗ r+1 g) = p q p∧q r=1 (r -1)! p -1 r -1 q -1 r -1 I p+q-2r (f ⊗ r g).
It follows that

E a - 1 q DF, DG H 2 (3.12) =      a 2 + p 2 p r=1 (r -1)! 2 p-1 r-1 2 q-1 r-1 2 (p + q -2r)! f ⊗ r g 2 H ⊗(p+q-2r) if p < q, (a -p! f, g H ⊗p ) 2 + p 2 p-1 r=1 (r -1)! 2 p-1 r-1 4 (2p -2r)! f ⊗ r g 2 H ⊗(2p-2r) if p = q. If r < p ≤ q then f ⊗ r g 2 H ⊗(p+q-2r) ≤ f ⊗ r g 2 H ⊗(p+q-2r) = f ⊗ p-r f, g ⊗ q-r g H ⊗2r ≤ f ⊗ p-r f H ⊗2r g ⊗ q-r g H ⊗2r ≤ 1 2 f ⊗ p-r f 2 H ⊗2r + g ⊗ q-r g 2 H ⊗2r . If r = p < q, then f ⊗ p g 2 H ⊗(q-p) ≤ f ⊗ p g 2 H ⊗(q-p) ≤ f 2 H ⊗p g ⊗ q-p g H ⊗2p . If r = p = q, then f ⊗ p g = f, g H ⊗p .
By plugging these last expressions into (3.12), we deduce immediately the desired conclusion.

Let us now recall the following result, which is a collection of some of the findings contained in the papers by Peccati and Tudor [START_REF] Peccati | Gaussian limits for vector-valued multiple stochastic integrals[END_REF] and Nualart and Ortiz-Latorre [START_REF] Nualart | Central limit theorem for multiple stochastic integrals and Malliavin calculus[END_REF]. Theorem 3.9 (See [START_REF] Nualart | Central limit theorem for multiple stochastic integrals and Malliavin calculus[END_REF][START_REF] Peccati | Gaussian limits for vector-valued multiple stochastic integrals[END_REF]). Fix d ≥ 2 and let C = {C(i, j) : i, j = 1, ..., d} be a d × d positive definite matrix. Fix integers 1 ≤ q 1 ≤ . . . ≤ q d . For any n ≥ 1 and i = 1, . . . , d, let f

(n) i belong to H ⊙q i . Assume that F (n) = (F (n) 1 , . . . , F (n) d ) := (I q 1 (f (n) 1 ), . . . , I q d (f (n) d )) n ≥ 1, is such that lim n→∞ E[F (n) i F (n) j ] = C(i, j), 1 ≤ i, j ≤ d. (3.13)
Then, as n → ∞, the following four assertions are equivalent:

(i) For every 1 ≤ i ≤ d, F (n) i
converges in distribution to a centered Gaussian random variable with variance C(i, i).

(ii) For every

1 ≤ i ≤ d, E (F (n) i ) 4 → 3C(i, i) 2 .
(iii) For every 1 ≤ i ≤ d and every

1 ≤ r ≤ q i -1, f (n) i ⊗ r f (n) i H ⊗2(q i -r) → 0.
(iv) The vector F (n) converges in distribution to a d-dimensional Gaussian vector N d (0, C). Moreover, if C(i, j) = δ ij , where δ ij is the Kronecker symbol, then either one of conditions (i)-(iv) above is equivalent to the following:

(v) For every 1 ≤ i ≤ d, DF (n) i 2 H L 2 -→ q i .
We conclude this section by pointing out the remarkable fact that, for vectors of multiple Wiener-Itô integrals of arbitrary length, the Wasserstein distance metrizes the weak convergence towards a Gaussian vector with positive definite covariance. Note that the next statement also contains a generalization of the multidimensional results proved in [START_REF] Nualart | Central limit theorem for multiple stochastic integrals and Malliavin calculus[END_REF] to the case of an arbitrary covariance. Proposition 3.10. Fix d ≥ 2, let C be a positive definite d × d symmetric matrix, and let 1 ≤ q 1 ≤ . . . ≤ q d . Consider vectors

F (n) := (F (n) 1 , . . . , F (n) d ) = (I q 1 (f (n) 1 ), . . . , I q d (f (n) d )), n ≥ 1, with f (n) i
∈ H ⊙q i for every i = 1 . . . , d. Assume moreover that F (n) satisfies condition (3.13). Then, as n → ∞, the following three conditions are equivalent:

(a) d W (F (n) , Z) → 0. (b) For every 1 ≤ i ≤ d, q -1 i DF (n) i 2 H L 2
-→ C(i, i) and, for every

1 ≤ i = j ≤ d, DF (n) i , -DL -1 F (n) j H = q -1 j DF (n) i , DF (n) j H L 2 -→ C(i, j). (c) F (n) converges in distribution to Z ∼ N d (0, C).
Proof. Since convergence in the Wasserstein distance implies convergence in distribution, the implication (a) → (c) is trivial. The implication (b) → (a) is a consequence of relation (3.11). Now assume that (c) is verified, that is, F (n) converges in law to Z ∼ N d (0, C) as n goes to infinity. By Theorem 3.9 we have that, for any i ∈ {1, . . . , d} and r ∈ {1, . . . , q i -1},

f (n) i ⊗ r f (n) i H ⊗2(q i -r) -→ n→∞ 0.
By combining Corollary 3.6 with Lemma 3.7 (see also Remark 3.8(2)), one therefore easily deduces that, since (3.13) is in order, condition (b) must necessarily be satisfied.

Applications

Convergence of marginal distributions in the functional Breuer-Major CLT

In this section, we use our main results in order to derive an explicit bound for the celebrated Breuer-Major CLT for fractional Brownian motion (fBm). We recall that a fBm B = {B t : t ≥ 0}, with Hurst index H ∈ (0, 1), is a centered Gaussian process, started from zero and with covariance function E(B s B t ) = R(s, t), where

R(s, t) = 1 2 t 2H + s 2H -|t -s| 2H ; s, t ≥ 0.
For any choice of the Hurst parameter H ∈ (0, 1), the Gaussian space generated by B can be identified with an isonormal Gaussian process of the type X = {X(h) : h ∈ H}, where the real and separable Hilbert space H is defined as follows: (i) denote by E the set of all R-valued step functions on [0, ∞), (ii) define H as the Hilbert space obtained by closing E with respect to the scalar product 1 [0,t] , 1 [0,s] H = R(t, s). In particular, with such a notation, one has that B t = X(1 [0,t] ). The reader is referred e.g. to [START_REF] Nualart | The Malliavin calculus and related topics (2 nd Edition)[END_REF] for more details on fBm, including crucial connections with fractional operators. We also define ρ(•) to be the covariance function associated with the stationary process

x → B x+1 -B x (x ∈ R), that is ρ(x) := 1 2 |x + 1| 2H + |x -1| 2H -2|x| 2H ∼ |x|→∞ H|2H -1| |x| 2H-2 .
Now fix an integer q ≥ 2, assume that H < 1 -1 2q and set

S n (t) = 1 σ √ n ⌊nt⌋-1 k=0 H q (B k+1 -B k ), t ≥ 0,
where H q is the qth Hermite polynomial defined as H q (x) = (-1) q e x 2 /2 d q dx q e -x 2 /2 , x ∈ R,

and where σ = q! r∈Z ρ 2 (r). According e.g. to the main results in [START_REF] Breuer | Central limit theorems for nonlinear functionals of Gaussian fields[END_REF] or [START_REF] Giraitis | CLT and other limit theorems for functionals of Gaussian processes[END_REF], one has the following CLT:

{S n (t), t ≥ 0} f.d.d.
-

d W S n (t i ) -S n (t i-1 ) √ t i -t i-1 1≤i≤d ; N d (0, I d ) ≤ c ×              n -1 2 if H ∈ (0, 1 2 ] n H-1 if H ∈ ( 1 2 , 2q-3 2q-2 ] n qH-q+ 1 2 if H ∈ ( 2q-3 2q-2 , 2q-1 2q )
.

Proof. Fix d ≥ 1 and t 0 = 0 < t 1 < . . . < t d . In the sequel, c will denote a constant independent of n, which can differ from one line to another. First, observe that

S n (t i ) -S n (t i-1 ) √ t i -t i-1 = I q (f (n) i ) with f (n) i = 1 σ √ n √ t i -t i-1 ⌊nt i ⌋-1 k=⌊nt i-1 ⌋ 1 ⊗q [k,k+1] .
In [START_REF] Nourdin | Stein's method on Wiener chaos[END_REF], proof of Theorem 4.1, it is shown that, for any i ∈ {1, . . . , d} and r ∈ {1, . . . , q i -1}:

f (n) i ⊗ r f (n) i H ⊗2(q i -r) ≤ c ×              n -1 2 if H ∈ (0, 1 2 ] n H-1 if H ∈ ( 1 2 , 2q-3 2q-2 ] n qH-q+ 1 2 if H ∈ ( 2q-3 2q-2 , 2q-1 2q ) . ( 4.2)
Moreover, when 1 ≤ i < j ≤ d, we have:

f (n) i , f (n) j H ⊗q = 1 σ 2 n √ t i -t i-1 √ t j -t j-1 ⌊nt i ⌋-1 k=⌊nt i-1 ⌋ ⌊nt j ⌋-1 l=⌊nt j-1 ⌋ ρ q (l -k) = c n ⌊nt j ⌋-⌊nt i-1 ⌋-1 |r|=⌊nt j-1 ⌋-⌊nt i ⌋+1 (⌊nt j ⌋ -1 -r) ∧ (⌊nt i ⌋ -1) -(⌊nt j-1 ⌋ -r) ∨ (⌊nt i-1 ⌋) ρ q (r) ≤ c ⌊nt i ⌋ -⌊nt i-1 ⌋ -1 n |r|≥⌊nt j-1 ⌋-⌊nt i ⌋+1 ρ(r) q = O n 2qH-2q+1 , as n → ∞, (4.3)
the last equality coming from

|r|≥N ρ(r) q = O( |r|≥N |r| 2qH-2q ) = O(N 2qH-2q+1 ), as N → ∞.
Finally, by combining (4.2), (4.3), Corollary 3.6 and Lemma 3.7, we obtain the desired conclusion.

Vector-valued functionals of finite Gaussian sequences

Let Y = (Y 1 , ..., Y n ) ∼ N n (0, I n ), and let f : R n → R be an absolutely continuous function such that f and its partial derivatives have subexponential growth at infinity. The following result has been proved by Chatterjee in [START_REF] Chatterjee | Fluctuation of eigenvalues and second order Poincaré inequalities[END_REF], in the context of limit theorems for linear statistics of eigenvalues of random matrices. We use the notation d TV to indicate the total variation distance between laws of real valued random variables. In what follows, we shall use Theorem 3.5 in order to deduce a multidimensional generalization of Proposition 4.2 (with the Wasserstein distance replacing total variation). where the functions T ab (•) are defined as T ab (y) = Proof. Without loss of generality, we can assume that Y = (Y 1 , ..., Y n ) = (X(h 1 ), ..., X(h n )), where X is an isonormal Gaussian process over some Hilbert space H, and h i , h l H = E(X(h i )X(h l )) = K(i, l). According to Theorem 3. By specializing the previous statement to the case n = d and f j (y) = y j , j = 1, ..., d, one obtains the following simple bound on the Wasserstein distance between Gaussian vectors of the same dimension (the proof is straightforward and omitted).

(

  ii) The Hilbert-Schmidt inner product and the Hilbert-Schmidt norm on the class of d × d real matrices, denoted respectively by •, • H.S. and • H.S. , are defined as follows: for every pair of matrices A and B, A, B H.S. := Tr(AB T ) and A H.S. := A, A H.S. . (iii) The operator norm of a d × d matrix A over R is given by

Lemma 3 . 3 .

 33 Fix an integer d ≥ 2 and let C = {C(i, j) : i, j = 1, ..., d} be a d × d positive definite symmetric real matrix. (i) Let Y be a random variable with values in R d . Then Y ∼ N d (0, C) if and only if,

2 .

 2 Unlike formulae (1.1) and (1.2) (associated with one-dimensional Gaussian approximations) the relation (3.2) and the Stein equation (3.3) involve second-order differential operators. A discussion of this fact is detailed e.g. in [4, Theorem 4].

  and g A (x) = g(Ax). Note that, sinceA -1 Z ∼ N d (0, I d ), the function h solves the Stein's equation x, ∇h(x) R d -∆h(x) = g A (x) -E[g A (Y )],where Y ∼ N d (0, I d ). We can now use the same arguments as in the proof of Lemma 3 in[START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF] to deduce that sup x∈R d Hess h(x) H.S. ≤ g A Lip ≤ A op g Lip .(3.5)

Theorem 3 . 5 .

 35 Fix d ≥ 2 and let C = {C(i, j) : i, j = 1, ..., d} be a d × d positive definite matrix. Suppose that Z ∼ N d (0, C) and that F = (F 1 , . . . , F d ) is a R d -valued random vector such that E[F i ] = 0 and F i ∈ D 1,2 for every i = 1, . . . , d. Then,

Proposition 4 . 2 (√ 1 -

 421 Lemma 5.3 in [3]). Assume that the random variable W = f (Y ) has zero mean and unit variance, and denote by Z ∼ N (0, 1) a standard Gaussian random variable. Then d TV (W, Z) ≤ 2Var(T (Y )) 1/2 , where the function T (•) is defined asT (y) = t Y ) dt.

Proposition 4 . 3 . 1

 431 Let Y ∼ N n (0, K), where K = {K(i, l) : i, l = 1, ..., n} is a n × n positive definite matrix. Consider absolutely continuous functions f j : R n → R, j = 1, ..., d. Assume that each random variable f j (Y ) has zero mean, and also that each function f j and its partial derivatives have subexponential growth at infinity. Denote by Z ∼ N d (0, C) a Gaussian vector with values in R d and with positive definite covariance matrix C = {C(a, b) : a, b = 1, ..., d}. Finally, write W = (W 1 , ..., W d ) = (f 1 (Y ), ..., f d (Y )). Then, d W (W, Z) ≤ C -1 op C C(a, b) -T ab (Y )) 2 ]

5 ,√ 1 -

 51 it is therefore sufficient to show that, for every a, b = 1, ..., d, T ab (Y ) = DW a , -DL -1 W b H . To prove this last claim, introduce the two H-valued functions Θ a (y) and Θ b (y), defined for y ∈ R d as follows: t Y ) h j dt.By using (2.2), it is easily seen that Θ a (Y ) = DW a . Moreover, by using e.g. formula(3.46) in[START_REF] Nourdin | Stein's method on Wiener chaos[END_REF], one deduces that Θ b (Y ) = -DL -1 W b . Since T ab (Y ) = Θ a (Y ), Θ b (Y ) H , the conclusion is immediately obtained.

Corollary 4 . 4 .

 44 Let Y ∼ N d (0, K) and Z ∼ N d (0, C), where K and C are two positive definite covariance matrices. Then dW (Y, Z) ≤ Q(C, K) × C -K H.S., whereQ(C, K) := min{ C -1 op C 1/2 op , K -1 op K 1/2 op }.

Acknowledgments. We thank an anonymous referee for insightful remarks. Part of this paper has been written while I. Nourdin and G. Peccati were visiting the Institute of Mathematical Sciences of the National University of Singapore, in the occasion of the "Workshop on Stein's Method", from March 31 until April 4, 2008. These authors heartily thank A. Barbour, L. Chen, K.-P. Choi and A. Xia for their kind hospitality and support.


Dus Mathias
email: mathias.dus@math.univ-toulouse.fr

Exponential stability of a general slope limiter scheme for scalar conservation laws subject to a dissipative boundary condition

Keywords: Exponential stability, bounded variation, scalar conservation laws, flux limiter scheme

In this paper, we establish the exponential BV stability of general systems of discretized scalar conservation laws with positive speed. The focus is on numerical approximation of such systems using a wide class of slope limiter schemes built from the upwind monotone flux. The proof is based on a Lyapunov analysis taken from the continuous theory [11] and a careful use of Harten formalism.

1 Introduction

Literature Review

In this work, we are interested in the exponential stability of systems of scalar conservation laws using boundary feedback laws. The system under consideration is a set of d nonlinear scalar conservation laws coupled at the boundary by a square matrix H of size d:

∀1 ≤ i ≤ d,    ∂ t R i + ∂x[f i (R i )] = 0 R i (t, 0) = [HR(t, 1)] i R i (0, x) = R 0 i (x) ∈ BV ([0, 1]) (1)
where R : R + × [0, 1] → R d (d ∈ N * ), f i : R → R. For coherence, it is assumed that all characteristic velocities are positive and consequently, the boundary condition in (1) is adapted. In [START_REF] Dus | Bv exponential stability for systems of scalar conservation laws using saturated feedback[END_REF], the author proves the well-posedness of such system using front-tracking techniques and a classical Kruzkov analysis. The case where the sign of characteristic velocities is not fixed is out of the scope of this article. This corresponds to a problem of traffic junction treated in [START_REF] Andreianov | Well-posedness for vanishing viscosity solutions of scalar conservation laws on a network[END_REF][2] just to mention a few.

More specifically, the focus is on the exponential stability of [START_REF] Alabau-Boussouira | Nonlinear damped partial differential equations and their uniform discretizations[END_REF]. The problem is equivalent to find sufficient conditions on H such that for any R 0 initial data, the solution to (1) converges exponentially fast toward zero in the sense that ∀t ≥ 0, ||R(t, .)|| X ≤ Ce -γt ||R 0 || X [START_REF] Andreianov | A theory of L 1 -dissipative solvers for scalar conservation laws with discontinuous flux[END_REF] where C, γ > 0 are constants independent on t and || • || X is a norm on a Banach space X.

In [7, Chapter 1], typical examples of systems modeled by hyperbolic PDEs with such feedback boundary conditions are cited; the telegrapher equations for electrical lines, the shallow water (Saint-Venant) equations for open channels [START_REF] Hayat | On boundary stability of inhomogeneous 2 × 2 1-D hyperbolic systems for the C 1 norm[END_REF], the isothermal Euler equations for gas flow in pipelines or even the Aw-Rascle equations [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] for road traffic. However it should be noted that, in previous examples, there is often an in-domain coupling which is not present in our model. Moreover, fluxes are sometimes not scalar which renders the analysis far more complicated. Hence, the family of systems we study, constitutes a simplified model for more realistic systems cited above. It introduces techniques helping in the complex studies on general systems of conservation and balance laws from [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]Chapter 1].

Additionally, the stabilization of similar systems with non-local terms receive more and more attention. We can cite [START_REF] Bernard | Adaptive output-feedback stabilization of non-local hyperbolic PDEs[END_REF] where the authors added a nonlocal zeroth order term to be stabilized. In this article, uncertainties on parameters and on the state of the system are allowed and an adaptive command built from an observer is designed. In [START_REF] Coron | Output feedback stabilization for a scalar conservation law with a nonlocal velocity[END_REF], authors propose a spectral analysis to stabilize a scalar linear transport equation with a non-local velocity. The control is exerted at the boundary. Then, by a Lyapunov analysis they prove a local stability result for the nonlinear version of the system. Spaces L p on [0, 1] are embedded with their canonical norms ||.|| L p . The space Lip loc (R + ; E) corresponds to functions from R + to a Banach space E which are locally Lipschitz. Let f be a C 1 function from R to R. Let a, b ∈ R, then if a = b, the notation:

f (b) -f (a) b -a := f (a) (3)
is imposed all along the article.

Stability of the PDE system without discretization

For the continuous PDE system, there are many results giving sufficient conditions on H to ensure exponential stability. One can mention the introduction of [START_REF] Dus | Bv exponential stability for systems of scalar conservation laws using saturated feedback[END_REF] for a review of some of these results. The criterion which is relevant in this paper is given here:

ρ∞(H) := inf ∆∈D + d (R) |∆H∆ -1 |∞ < 1. (4)
Owing (4), we proved the exponential stability in space X = BV ([0, 1]) of systems of scalar conservation laws with positive velocities of propagation in [START_REF] Dus | Bv exponential stability for systems of scalar conservation laws using saturated feedback[END_REF]. In this work, the Banach space X = BV ([0, 1]) is also considered and the focus is on the exponential stability of the corresponding numerical solutions.

The numerical problem

In this paper, the focus is on finite volume approximations of system (1):

R n+1 j -R n j dt + f n j+1/2 -f n j-1/2 dx = 0
where n, j are respectively the time and space index. More precisely, we find sufficient condition on H such that the discretized version of ( 1) is exponentially stable.

The upwind flux f n j+1/2 = f (R n j ) is known for its simplicity and its good properties. It is consistent and under a classic hypothesis of CF L, it is Total Variation Decreasing (TVD) and monotone [START_REF] Eymard | Solution of Equation in R n (Part 3)[END_REF] for scalar equations. These characteristics allow to prove easily that when the parameter of discretization tends towards zero, the numerical solution tends towards the unique entropy solution of the problem. For an introduction on the notion of entropy solution, we refer to [START_REF] Eymard | Solution of Equation in R n (Part 3)[END_REF]. However, it is over diffusive and only first order accurate. In order to correct such behavior, one can use the additional precision of a second order scheme taking:

f n j+1/2 = f (R n j + Rn j )
where for example Rn

j = R n j+1 -R n j dx dx 2
. However, a second order scheme cannot be TVD by Godunov Theorem and have oscillatory behaviors around discontinuities. Moreover, it can also be unstable as it is the case for the example cited just above. This is why one introduces a slope limiter φ [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] whose role is to damp the effect of the second order flux around discontinuities. The slope limiter scheme is defined taking:

Rj := φ R n j -R n j-1 dx , R n j+1 -R n j dx dx 2
where φ is called the slope limiter that will be defined carefully later. The main contribution of the paper is to find sufficient condition on φ and H to prove the exponential BV ([0, 1]) stability of the numerical solution built from the slope limiter scheme giving an explicit formula for the dissipation rate. There are several results for the stabilization of discretized 1D wave equation. A very instructive survey from Trélat [START_REF] Trélat | Stabilization of Semilinear PDEs, and Uniform Decay under Discretization[END_REF] gives an insight of techniques for the stabilization of discretized hyperbolic PDEs. One can also cite [START_REF] Tebou | Uniform boundary stabilization of the finite difference space discretization of the 1 -d wave equation[END_REF] where the authors prove the boundary stabilization of a semi-discretized 1D wave equation. Adding a numerical viscosity term damping high frequencies, they prove a uniform exponential stability decay of the energy uniformly with respect to the parameter of discretization. The article [START_REF] Tébou | Uniform exponential long time decay for the space semidiscretization of a locally damped wave equation via an artificial numerical viscosity[END_REF] deals with the case of an in-domain damping.

The main concern is to obtain a uniform observability inequality whose constant is independent on the parameter of discretization. This is not an easy task since the observability of the continuous system does not necessarily imply the same for the fully discrete or the semi discrete system. The survey [START_REF] Zuazua | Propagation, observation, and control of waves approximated by finite difference methods[END_REF] summarizes most of the main results about this question. The main problem comes from high frequency modes and their group velocity. The group velocity of two consecutive high frequency modes tends to zero as the discretization parameter goes to zero which prevent from using Ingham's inequality. This makes the observation (at a finite fixed time) very difficult as the discretization get finer [START_REF] Marica | Localized solutions for the finite difference semi-discretization of the wave equation[END_REF]. This is why it is primordial to add a non physical numerical viscosity in order to damp these high frequencies.

For completeness, we cite [START_REF] Ervedoza | Observability properties of a semi-discrete 1D wave equation derived from a mixed finite element method on nonuniform meshes[END_REF] where the author prove the uniform exponential decay for another semi-discretized scheme simulating an in-domain damped wave equation. The case of semilinear wave equation is treated in [START_REF] Alabau-Boussouira | Nonlinear damped partial differential equations and their uniform discretizations[END_REF]. The method to deduce exponential stabilization for the full discretized case is given in [START_REF] Ervedoza | Uniformly exponentially stable approximations for a class of damped systems[END_REF].

However to the author's knowledge, few results are known about the stabilization of discretized transport equations of the form (1). The method presented here is very different from observation inequalities techniques. We rather use discrete Lyapunov functions to prove exponential decay of our solution. In a L 2 framework, the articles [START_REF] Gerster | Discretized feedback control for systems of linearized hyperbolic balance laws[END_REF][START_REF] Banda | Numerical discretization of stabilization problems with boundary controls for systems of hyperbolic conservation laws[END_REF] prove a stabilization result for discretized linear balance laws coupled in the domain and at the boundary using such techniques. However, it does not seem to be generalizable to a BV context for nonlinear scalar conservation laws. In section 3, we address this problem introducing a discrete Lyapunov functional in BV.

Problem setting

For wellposedness and scheme convergence results, we consider a more general boundary condition:

   ∂ t R + ∂x[f (R)] = 0 R(t, 0) = g(R(t, 1)) R(0, x) = R 0 (x) ∈ BV ([0, 1]).
(

) 5
where the following hypothesis on f and g hold:

f is diagonal in the sense that:

∀1 ≤ i ≤ d, R ∈ R d , f i (R) = f i (R i ).
Moreover, all the f i seen as functions from R to R are non decreasing. Doing so, the boundary condition of (5) has sense. -There exist two positive constants 0 < v min < vmax such that for all 1 ≤ i ≤ d,

∀(u, v) ∈ R 2 , v min |u -v| ≤ |f i (u) -f i (v)| ≤ vmax|u -v|. ( 6
)
-The boundary feedback function g : R d → R d is Lipschitz with Lipschitz constant Lg > 0. Moreover, we impose g(0) = 0.

Remark 1 The components of a solution R to (5) interact with each others only at the boundary.

In this paper, the framework is BV ([0, 1]). This space is embedded with its classical norm || • || BV ([0,1]) defined by

∀R ∈ BV ([0, 1]), ||R|| BV ([0,1]) = T V [0,1] (R) + ||R|| L 1 ([0,1]) (7)
where:

T V (R) = sup n, (x1,...,xn) n-1 i=1 |R(x i+1 ) -R(x i )| .
The reason why the BV space is considered is because any function with bounded variations has a left limit and a right limit at each point x of [0, 1]. Hence, it is easy to define the trace operator and impose a boundary condition. Moreover, BV ([0, 1]) has a very interesting property of compactness which will be very useful when we will pass to the limit when the parameter of discretization tends towards zero. These properties are summed up in a lemma and a theorem (see references [9, Lemma 2.1, Theorem 2.4]) :

Lemma 1 Let R : [0, 1] → R d with bounded variations. Then for all x ∈ (0, 1), the left and right limit

R(x -) = lim y→x - R(y), R(x + ) = lim y→x + R(y)
exist. Moreover, R(0 + ) and R(1 -) are also well defined and R has at most countably many point of discontinuities.

The following theorem is from Helly and states the compactness of BV ([0, 1]) in L 1 loc (R + , L 1 ([0, 1])).

Theorem 1 Let (Rν )ν be a sequence of functions from R + × [0, 1] to R d such that there exists constants C, M and L satisfying

∀ν > 1, x ∈ [0, 1], t ≥ 0, T V [0,1] (Rν (t, .)) ≤ C, |Rν (t, x)| ≤ M (8)
and

∀0 ≤ t, s ≤ T, ||Rν (t, .) -Rν (s, .)|| L 1 ([0,1]) ≤ L|t -s|. (9)
Then there exists a subsequence (Rµ)µ converging toward a certain R in

L 1 loc (R + , L 1 ([0, 1]
)) and this limit satisfies (8)- [START_REF] Bressan | Hyperbolic systems of conservation laws[END_REF] with Rν replaced by R.

When the initial data is supposed to be L ∞ only, an entropic solution exists and is unique [START_REF] Kružkov | First order quasilinear equations with several independent variables[END_REF]. Moreover, it is possible to define the trace in a certain sense [START_REF] Yu | Existence of strong traces for generalized solutions of multidimensional scalar conservation laws[END_REF][START_REF] Vasseur | Strong traces for solutions of multidimensional scalar conservation laws[END_REF]. However, the method we use to prove exponential stability crucially uses the BV hypothesis and this is why the L ∞ framework is out of the scope of this paper.

Outline

The article is organized as follows. In section 2, we precisely define the slope limiter scheme and prove the convergence in BV ([0, 1]) of the corresponding numerical solution. As a by-product, a well-posedness result for system (5) is proved. In section 3, the framework is restricted to linear boundary feedback operator g(R) = HR. Then, the exponential stability of the scheme in BV ([0, 1]) is proved giving an explicit formula for the dissipation rate. Passing to the limit, we give another proof of the exponential stability of the continuous system which was already stated in [START_REF] Dus | Bv exponential stability for systems of scalar conservation laws using saturated feedback[END_REF]. In section 4, simulations are given to illustrate results from section 3. A study of saturated controls is also given. Finally, the conclusions and perspectives are exposed in the last part of this work.

Well-posedness and numerical approximation results

In this section, we introduce rigorously the slope limiter scheme corresponding to [START_REF] Banda | Numerical discretization of stabilization problems with boundary controls for systems of hyperbolic conservation laws[END_REF] and prove the convergence of such a scheme. Firstly, let us consider a space step dx = 1/N (N ∈ N * ) and a time step dt > 0 such that the following CFL condition holds (recall the definition of vmax in ( 6)):

vmaxν ≤ 1 -ξ (10)
with

ν := dt dx
and where ξ is a real number such that:

0 < ξ < 1.
Doing so, the space-time mesh is given by:

∀n ∈ N, 1 ≤ j ≤ N, x j := (j -1/2)dx t n := ndt , C j := (x j -dx/2, x j + dx/2). - - - - - - - - - - - - - - - - - C 1 C 2 • • • • • • C N x 1 x 2 • • • • • • x N Fig. 1:
The space grid All along this paper, the superscript n is exclusively allocated to designate time indices. Whereas, subscripts i and j are the indices corresponding respectively to the vector component and space grid. The numerical approximation R

0 ∆x of R 0 ∈ BV ([0, 1]) is given below ∀1 ≤ j ≤ N, x ∈ C j , R 0 ∆x (x) = 1 dx Cj R 0 (x)dx. ( 11
)
In our context, the definition of a four point finite volume scheme is given below:

Definition 1 A four-point finite volume scheme is defined a follows:

-The scheme is initialized defining (R 0 j ) 1≤j≤N as in ( 11). -Then the numerical solution is built by induction. For n ∈ N, (R n j ) 1≤j≤N is supposed to be given. We define (R n+1 j ) 1≤j≤N using the following relation

∀ 2 ≤ j ≤ N -1, R n+1 j -R n j dt + f n j+1/2 -f n j-1/2 dx = 0 ( 12
)
where f n j+1/2 is the approximation of the flux at the right interface of cell j:

∀1 ≤ j ≤ N -1, f n j+1/2 = f (R n j-1 , R n j , R n j+1 ). ( 13
)
The function f is called the three point approximation of the flux f . The choice of such a function will be given later with careful justifications. -At the boundary the scheme degenerates into the upwind one:

           R n+1 N -R n N dt + f (R n N ) -f (R n N -1 ) dx = 0 R n+1 1 -R n 1 dt + f (R n 1 ) -f (R n 0 ) dx = 0 R n 0 := g(R n N ). (14)
There are many possibilities for the choice of the function f . For example, one can take:

f (R n j-1 , R n j , R n j+1 ) = f (R n j
). It corresponds to the classic upwind scheme when the transport velocity is positive. Recalling what was said in the introduction, the first main advantage of such scheme is that it is TVD and consistent allowing convergence towards a weak solution to [START_REF] Banda | Numerical discretization of stabilization problems with boundary controls for systems of hyperbolic conservation laws[END_REF]. The other advantage is that it is monotone forcing the scheme to converge towards the unique entropy solution of (5). However, it is known for its diffusivity giving poor results in simulation. This is why one introduces the slope limiter flux [31] [23]:

f (R n j-1 , R n j , R n j+1 ) = f (R n j + Rn j ) (15)
where the term Rn j is given by a piecewise linear approximation of the solution:

∀1 ≤ i ≤ d, Rn i,j = φ R n i,j -R n i,j-1 dx , R n i,j+1 -R n i,j dx dx 2
where φ is assumed to be of the form: [START_REF] Ervedoza | Uniformly exponentially stable approximations for a class of damped systems[END_REF] such that there exists 0 ≤ β < 1/2 such that:

φ(u, v) = 0 if v = 0 φr(u/v)ψ dx (v) otherwise
|ψ dx (v)| ≤ v and |ψ dx (v)| ≤ dx -β (17)
and φr : R → R and for all r ≤ 0, φr(r) = 0. The natural choice for ψ dx is:

ψ dx (v) = v if |v| ≤ dx -β sign(v)dx -β otherwise. Hence, if R n i,j+1 -R n i,j = 0, then: Rn i,j := φr r n i,j+1/2 ψ dx R n j+1 -R n j dx dx 2 .
The variable r n i,j+1/2 ∈ R d is defined by:

∀1 ≤ i ≤ d, r n i,j+1/2 := R n i,j -R n i,j-1 R n i,j+1 -R n i,j
.

Remark 2 When the flux is linear, slope limiter schemes and flux limiter schemes coincide [23, p. 114].

The slope limiter φ is designed such that when the solution is regular i.e r j+1/2 ≈ 1, the flux f is closed to a second order flux. On the contrary, when the solution is not regular (r j+1/2

1) or admits a local extrema (r j+1/2 ≤ 0), the second order scheme gives poor oscillatory results. In this case, we prefer to give more weight on the upwind TVD scheme imposing φ ≈ 0. We refer to [23, Section 6.6-6.9] for a more complete description of how slope limiter schemes were designed.

We define the following coefficient for all 1 ≤ i ≤ d, 1 ≤ j ≤ N -1, which will be useful all along this article:

a n i,j-1 := ν f i (R n i,j + Rn i,j ) -f i (R n i,j-1 + Rn i,j-1 ) R n i,j + Rn i,j -R n i,j-1 -Rn i,j-1 × 1 + φr(r n i,j+1/2 ) 2r n i,j+1/2 ψ dx ( R n i,j+1 -R n i,j dx ) R n i,j+1 -R n i,j dx - φr(r n i,j-1/2 ) 2 ψ dx ( R n i,j -R n i,j-1 dx ) R n i,j -R n i,j-1 dx . (18)
For the pathologic cases, one defines a n i,j-1 as follows:

-If R n i,j -R n i,j-1 = 0 (ie r n i,j+1/2 = 0 and r n i,j-1/2 ill-defined), then a n i,j-1 := νf i (R n i,j ). -Otherwise, if R n i,j -R n i,j-1 = 0 and R n i,j+1 -R n i,j = 0 (r n i,j+1/2 is ill-defined), the coefficient φr(r n i,j+1/2 ) 2r n i,j+1/2 ψ dx ( R n i,j+1 -R n i,j dx ) R n i,j+1 -R n i,j dx
is replaced by zero in [START_REF] Gerster | Discretized feedback control for systems of linearized hyperbolic balance laws[END_REF].

-If R n i,j + Rn i,j -R n i,j-1 -Rn i,j-1 = 0, then recall that fi(R n i,j + Rn i,j )-fi(R n i,j-1 + Rn i,j-1 ) R n i,j + Rn i,j -R n i,j-1 -Rn i,j-1 = f i (R n i,j + Rn i,j
) by the convention of notation (3). For j = N , we set:

a n i,N -1 := ν f i (R n i,N ) -f i (R n i,N -1 ) R n i,N -R n i,N -1
.

The diagonal matrix with a n i,j as entries is denoted A n j := diag{a n i,j | 1 ≤ i ≤ d}. Hence, the scheme can be reformulated under an upwind form:

R n+1 j = R n j -A n j-1 (R n j -R n j-1 ).
In order to ensure that our scheme is TVD, we impose the following condition of the slope limiter φ:

Hypothesis 1 There exists 0 ≤ vnum ≤ v min such that:

0 ≤ φr(r) ≤ min 2 1 νvmax -1 r, 2 1 - vnum v min
This hypothesis is fundamental to ensure Harten's condition [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF]:

Lemma 2 If Hypothesis 1 is satisfied, then Harten's condition is satisfied:

νvnum ≤ a i n,j ≤ 1. ( 19
)
Proof Looking at [START_REF] Gerster | Discretized feedback control for systems of linearized hyperbolic balance laws[END_REF], one has:

a n i,j ≤ νvmax 1 + φr(r n i,j+1/2 ) 2r n i,j+1/2 , a n i,j ≥ νv min 1 - φr(r n i,j-1/2 ) 2 .
By Hypothesis 1:

νvmax 1 + φr(r n i,j+1/2 ) 2r n i,j+1/2 ≤ 1, νv min 1 - φr(r n i,j-1/2 ) 2 ≥ νvnum
which immediately gives the result of the Lemma.

Before going into the main result of this section, we introduce a new notation. Let R ∆x be the piecewise constant function equal to R n j on each cell [ndt, (n + 1)dt] × C j . This numerical approximation helps us proving the existence of a solution to [START_REF] Banda | Numerical discretization of stabilization problems with boundary controls for systems of hyperbolic conservation laws[END_REF]. Meanwhile, Kruzhkov theory allows to prove the uniqueness. To prove the following Theorem 2, techniques are quite classical and plenty of works deal with the convergence of slope limiter schemes for scalar conservation laws in several dimensions and when the domain is the whole real line [START_REF] Chainais-Hillairet | Second-order finite-volume schemes for a non-linear hyperbolic equation: error estimate[END_REF] [21] [START_REF] Osher | On the convergence of difference approximations to scalar conservation laws[END_REF]. In our case, this is a bit more complex since we have to take into account the boundary condition.

Theorem 2 System (5) admits a unique entropy solution R ∈ L ∞ loc (R + , BV ([0, 1])) ∩ Lip loc (R + , L 1 ([0, 1]
)) in the sense that for all T > 0:

-The entropy decay estimate holds:

∀k ∈ R d , d i=1 T 0 1 0 |R i -k i |∂ t ϕ i +(f i (R i )-f i (k i ))sign(R i -k i )∂xϕ i dxdt ≥ 0 (20) for all ϕ i ≥ 0 and ϕ i ∈ C 1 c (]0, T [×]0, 1[).
-R(0, .) = R 0 in the almost everywhere sense.

-R(., 0 + ) = g(R(., 1 -)) in the almost everywhere sense.

Moreover, the following convergence properties hold

∀t ≥ 0, dt > 0, T V [0,1] (R(t, .)) ≤ lim sup N →+∞ sup s∈[t,t+dt] T V [0,1] (R ∆x (s, .)). ( 21
)
∀t ≥ 0, lim

N →∞ ||R(t, .) -R ∆x (t, .)|| L 1 ([0,1]) = 0. ( 22
)
Proof The proof is given in Appendix A.

3 Exponential BV stability

In this section, we focus on a particular case of system of (5).

   ∂ t R + ∂x[f (R)] = 0 R(t, 0) = HR(t, 1) R(0, x) = R 0 (x) ∈ BV ([0, 1]). ( 23
)
where H ∈ M d (R). We introduce the following stability hypothesis.

Hypothesis 2

The feedback matrix H satisfies:

ρ 1 (H) < 1. Remark 3 By [11, Remark 1.4] ∀M ∈ M d (R), ρ 1 (H) = ρ∞(H)
so that Hypothesis 2 can also be written as:

ρ∞(H) < 1.
The following Lyapunov functional is proposed for functions which are constant on the space mesh.

Definition 2 Let N be in N * , γ > 0 and P ∈ D + d (R). For any function R piecewise constant on cells C j taking its values in R d , the BV Lyapunov functional L is given by

L(R) = d i=1 P i N -1 j=0 |R i,j+1 -R i,j |e -γxj
where

x 0 := -dx/2, R 0 := HR N .
Next lemma ensures the equivalence between L and || • || BV ([0,1]) defined in [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF].

Lemma 3 Suppose Hypothesis 2. For all P ∈ D + d (R + ) such that |P HP -1 | 1 ≤ 1 and for all γ > 0, there exists a constant C(P, γ, H) > 1 such that for all N ∈ N * and all functions R piecewise constant on cells (C j ) j :

L(R) C(P, γ, H) ≤ ||R|| BV ([0,1]) ≤ C(P, γ, H)L(R). ( 24
)
Proof We define the extension R of R at C 0 := (-dx, 0) setting R(x) = HR(1) on C 0 . In [13, Lemma 2], it is proved that:

L(R) C(P, γ, H) ≤ T V ( R) + R L 1 ([0,1]) ≤ C(P, γ, H)L(R). (25)
Additionally, as R is an extension on the left boundary of R and by the definition of the BV norm, it holds:

R BV ([0,1]) ≤ T V ( R) + R L 1 ([0,1]) = R BV ([0,1]) + |R 1 -R0 | + | R0 |dx.
As R0 = HR N by definition and as

• L ∞ ([0,1]) ≤ C • BV ([0,1]) (C >
0 is a constant depending on the parameters of the problem only), one has:

R BV ([0,1]) ≤ T V ( R) + R L 1 ([0,1]) ≤ C(H) R BV ([0,1]) .
Hence, using [START_REF] Osher | On the convergence of difference approximations to scalar conservation laws[END_REF], it holds:

L(R) C(P, γ, H) ≤ ||R|| BV ([0,1]) ≤ C(P, γ, H)L(R).
where we may have changed the constant C(P, γ, H). The lemma is proved.

The reason why we introduce a discrete Lyapunov functional is because we cannot deduce the exponential stability of the numerical solution from the stability of the "continuous" one proven in [START_REF] Dus | Bv exponential stability for systems of scalar conservation laws using saturated feedback[END_REF]. Indeed, the convergence of the numerical solution towards the continuous solution is only proven to be in [START_REF] Eymard | Solution of Equation in R n (Part 3)[END_REF]Theorem 5.2]. This is not enough to prove an exponential stability result even in L 1 ([0, 1]). This is why, we need to study the stability property of the scheme in itself introducing a discrete Lyapunov functional. The main result of this paper is given in the following theorem whose proof is the object of the section.

L 1 loc (R + ; L 1 ([0, 1])) in general
Theorem 3 Under Hypothesis 1-2 and for all γ > 0, 0 ≤ c < vnum, P ∈ D + d (R d )

such that |P HP -1 | 1 < e -γ , there exists ε such that if γdt < ε, the discrete Lyapunov estimate holds:

∀n ≥ 0, L(R n ) ≤ e -cγndt L(R 0 ).
To understand the consequences of Theorem 3, we consider a limiter φr satisfying Hypothesis 1 with vnum = αv min and α ∈ [0, 1]. Then, a Sweby diagram [START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF] is drawn in Figure 2 where the colored region represents the zone corresponding to Hypothesis 1. Theorem 3 states that for such φr the dissipation rate of the numerical solution is bounded from below by αv min (= vnum). Remark that the best estimation vnum = v min ⇐⇒ α = 1 is obtained for only one scheme φr = 0 which is the upwind scheme.

y = 2( 1 νvmax -1)r r y y = 2(1 -α) φr(r)

Fig. 2: The Sweby diagram

Proof Let N, n be positive integers. By the definition of L, we have:

L(R n+1 ) = d i=1 P i N -1 j=0 |R n+1 i,j+1 -R n+1 i,j |e -γxj .
The ith component of the Lyapunov functional L is given here:

L i (R n+1 ) := N -1 j=0 |R n+1 i,j+1 -R n+1 i,j |e -γxj .
For 1 ≤ i ≤ d and 1 ≤ j ≤ N -1, we estimate the intensity of the discontinuity

R n+1 i,j+1 -R n+1 i,j . R n+1 i,j+1 -R n+1 i,j = (1 -a n i,j )(R n i,j+1 -R n i,j ) + a n i,j-1 (R n i,j -R n i,j-1
) where a n j was defined in [START_REF] Gerster | Discretized feedback control for systems of linearized hyperbolic balance laws[END_REF]. Because of the Harten's condition [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF] we have for 1 ≤ i ≤ d and 1 ≤ j ≤ N -1:

|R n+1 i,j+1 -R n+1 i,j | ≤ (1 -a n i,j )|R n i,j+1 -R n i,j | + a n i,j-1 |R n i,j -R n i,j-1 |. (26)
Multiplying by e -γxj and summing over 1 ≤ j ≤ N -1, we get

N -1 j=1 |R n+1 i,j+1 -R n+1 i,j |e -γxj ≤ N -2 j=1 (1 -a n i,j (1 -e -γdx ))|R n i,j+1 -R n i,j |e -γxj +(1 -a n i,N -1 )|R n i,N -R n i,N -1 |e -γx N -1 +a n i,0 |R n i,1 -R n i,0 |e -γdx/2 .
Coming back to the Lyapunov functional L i , it holds:

L i (R n+1 ) ≤ N -2 j=1 (1 -a n i,j (1 -e -γdx ))|R n i,j+1 -R n i,j |e -γxj +(1 -a n i,N -1 )|R n i,N -R n i,N -1 |e -γx N -1 +a n i,0 |R n i,1 -R n i,0 |e -γdx/2 +|R n+1 i,1 -R n+1 i,0 |e γdx/2 .

Now, we estimate

|R n+1 i,1 -R n+1 i,0 |: R n+1 i,1 -R n+1 i,0 = R n i,1 + a n i,0 (R n i,0 -R n i,1 ) -R n+1 i,0 = (R n i,1 -R n i,0 )(1 -a n i,0 ) + R n i,0 -R n+1 i,0 . (27)
Taking the absolute value and by triangle inequality, one obtains:

|R n+1 i,1 -R n+1 i,0 | ≤ (1 -a n i,0 )|R n i,1 -R n i,0 | + |R n i,0 -R n+1 i,0 |.
This gives:

L i (R n+1 ) ≤ N -2 j=1 (1 -a n i,j (1 -e -γdx ))|R n i,j+1 -R n i,j |e -γxj +(1 -a n i,N -1 )|R n i,N -R n i,N -1 |e -γx N -1 +(1 -a n i,0 (1 -e -γdx ))|R n i,1 -R n i,0 |e γdx/2 +|R n+1 i,0 -R n i,0 |e γdx/2
which is no more than:

L i (R n+1 ) ≤ N -2 j=0 (1 -a n i,j (1 -e -γdx ))|R n i,j+1 -R n i,j |e -γxj +(1 -a n i,N -1 )|R n i,N -R n i,N -1 |e -γx N -1 +|R n+1 i,0 -R n i,0 |e γdx/2 .
Using the fact that there exists ε such that for γdx ≤ ε, 1 -e -γdx ≥ c vnum γdx (because c vnum < 1 by assumption), we get:

L i (R n+1 ) ≤ N -2 j=0 (1 -a n i,j cγdx vnum )|R n i,j+1 -R n i,j |e -γxj +(1 -a n i,N -1 )|R n i,N -R n i,N -1 |e -γx N -1 +|R n+1 i,0 -R n i,0 |e γdx/2 .
Then we add and subtract (1

-a n i,N -1 cγdx vnum )|R n i,N -R n i,N -1 |e -γ(N -1
)dx in the first and second line respectively to get:

L i (R n+1 ) ≤ N -1 j=0 (1 -a n i,j cγdx vnum )|R n i,j+1 -R n i,j |e -γxj -a n i,N -1 (1 - cγdx vnum )|R n i,N -R n i,N -1 |e -γx N -1 +|R n+1 i,0 -R n i,0 |e γdx/2 .
Finally, using the fact that

a n N -1 (R n N -R n N -1 ) = -(R n+1 N -R n N
) and the fact that 1 -a n i,j cγdx vnum ≤ 1 -cγdt (see [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF]), the estimate of L i writes:

L i (R n+1 ) ≤ (1 -cγdt)L i (R n ) -(1 - cγdx vnum )|R n+1 i,N -R n i,N |e -γx N -1 +|R n+1 i,0 -R n i,0 |e γdx/2 .
Using the discrete boundary condition ( 14), one changes the last boundary term:

L i (R n+1 ) ≤ (1 -cγdt)L i (R n ) -(1 - cγdx vnum )|R n+1 i,N -R n i,N |e -γx N -1 +|[HR n+1 N ] i -[HR n N ] i |e γdx/2 .
Multiplying by P i and summing over all the 1 ≤ i ≤ d, an estimate on L writes:

L(R n+1 ) ≤ (1 -cγdt)L(R n ) - d i=1 P i |R n+1 i,N -R n i,N |(1 - cγdx vnum )e -γx N -1 + d i=1 P i |[HR n+1 N ] i -[HR n N ] i |e γdx/2 .
Using the definition of |•| 1 , the boundary terms can be compared with each others:

L(R n+1 ) ≤ (1 -cγdt)L(R n ) +(|P HP -1 | 1 e γdx/2 -(1 - cγdx vnum
)e -γx N -1 )

d i=1 P i |R n+1 i,N -R n i,N |. As P ∈ D + d (R) is such that |P HP -1 | 1 < e -γ
, we get for γdx sufficiently small (that is to say γdx < ε(c, vnum, H)):

L(R n+1 ) ≤ 1 -cγdt L(R n ).
Summing with respect to time, one obtains:

L(R n+1 ) ≤ 1 -cγdt n+1 L(R 0 ).
Using the fact that 1 -x ≤ e -x , the result of Lemma 3 is proved.

Remark 4 The rate of exponential convergence γc given by Theorem 3 is weaker than the one given in [START_REF] Dus | Bv exponential stability for systems of scalar conservation laws using saturated feedback[END_REF] for the continuous setting γv min . This is mainly due to the choice of the limiter which guaranties only [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF]. To recover the rate of convergence from [START_REF] Dus | Bv exponential stability for systems of scalar conservation laws using saturated feedback[END_REF], one should impose:

νv min ≤ a i n,j ≤ 1
which is true for the upwind scheme but not for other limiters.

The following Theorem was already proved in [START_REF] Dus | Bv exponential stability for systems of scalar conservation laws using saturated feedback[END_REF]. The main difference here is that the proof comes from the convergence of our finite volume scheme.

Theorem 4 Under Hypothesis 2 and for all

P ∈ D + d (R), γ > 0 such that |P HP -1 | 1 < e -γ ,
there exists a constant C(P, γ, H) depending only on P, γ, H such that the entropy solution to [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]

verifies ∀t ≥ 0, ||R(t, •)|| BV ([0,1]) ≤ C(P, γ, H) exp -γv min t ||R 0 || BV ([0,1]) .
Proof For this proof, we take φ ≡ 0 which corresponds to the upwind scheme. As a consequence, it is possible to take vnum = v min in Hypothesis 1. Then, it suffices to pass to the limit as N goes to infinity in Lemma 3 using estimates [START_REF] Kružkov | First order quasilinear equations with several independent variables[END_REF], [START_REF] Kröner | Convergence of higher order upwind finite volume schemes on unstructured grids for scalar conservation laws in several space dimensions[END_REF] and the equivalence between norms L, || • || BV ([0,1]) (see [START_REF] Marica | Localized solutions for the finite difference semi-discretization of the wave equation[END_REF]).

Simulations

In this section, we illustrate numerically that slope limiter schemes are clearly less dissipative than the upwind scheme and capture better the behavior of the continuous solution. For computations, we consider the case where d = 4, dt dx = 0.4 and the flux is given by:

∀1 ≤ i ≤ 4, x ∈ R, f i (R i ) = 0.5R i + 0.5 2 i-1 arctan(R i ).
In such a case, min and max velocities are v min = 0.5 and vmax = 1. Concerning the limiter, it is chosen as:

φr(r) = max{min{r, 1}, 0}, ψ dx (v) = v if |v| ≤ dx -β sign(v)dx -β otherwise
with β = 0.49. This limiter is called a minmod limiter.

Dissipativity of the scheme

If the focus is on the dissipativity of the scheme alone, a relevant boundary condition is:

H =     0 0 0 1 1 0 0 0 0 1 0 0 0 0 1 0     .
In such a case, the boundary condition does not dissipate energy. For a smooth initial data:

R 1 (t = 0, .) = • • • = R 4 (t = 0, .) = cos(2πx),
we compare the BV ([0, 1]) norm of the solution given by the upwind scheme with the one given by the limiter one. It is expected that the limiter scheme is less dissipative than the upwind one when the solution is regular. Fig. 3: The BV norm of the solution for the non dissipative case One sees in Figure 3a-3b that this is indeed the case. Obviously when the grid is fine enough, the dissipation created by the upwind scheme is negligible and results for both schemes are similar. Next, we can do the same analysis imposing this time, a dissipative feedback boundary condition:

H =     0.5 0.1 0 0.1 0.1 0.5 0.1 0 0 0.1 0.5 0.1 0.1 0 0.1 0.5     .
An optimization routine gives ρ 1 (H) = 0.7. Moreover, a less regular initial data is examined:

       R 0 1 = x ≥ 0.5 R 0 2 = -(x < 0.5) + (x ≥ 0.5) R 0 3 = x ≥ 0.5 R 0 4 = -(x < 0.5) + (x ≥ 0.5).
One obtains Figures 4a-4b for the evolution of the BV norm:

(a) dt = 10 -2 (b) dt = 5 × 10 -4
Fig. 4: The BV norm of the solution for the dissipative case

The black dashed line corresponds to the estimation of the rate of convergence γvnum (vnum = 0.5v min , γ = -log(0.7) here) given by Theorem 3. We see this is not optimal since the rate of convergence is underestimated. Then, when the solution is small enough, it should be wise to take v min = min i f i (0) in Theorem 3 in order to estimate the rate of convergence when the solution is small (in L ∞ ([0, 1])). This rate is represented by the dashed dotted line for the upwind scheme (vnum = v min ) in figures 4a-5b. This estimation fits better with the numerical experiments. With these results, one sees that Theorem 3 gives only a lower bound on the exponential convergence rate and not a precise estimation of it. Another remark to make is that the limiter scheme captures better the dissipation rate, the upwind scheme dissipating too much when the grid is coarse.

If we change the limiter, we have similar numerical results. For example, one can take the Van Leer limiter:

φr(r) = r + |r| 1 + |r| . (a) dt = 10 -2 (b) dt = 5 × 10 -4
Fig. 5: The Van-Leer case

For completeness, we plot a part of the left trace of the solution for t < 15 (dt = 5 × 10 -4 ). One clearly sees that the solution is rapidly damped when time evolves.

Comparison with the exact solution

An interesting analysis would be to compare convergence rates when it is possible to calculate the exact solution using the characteristics method. This is done by taking d = 2 and a linear flux given by:

f 1 (R 1 ) = R 1 f 2 (R 2 ) = 1.6 × R 2 .
The initial data is taken as follows:

R 0 1 = cos(2πx) R 0 2 = sin(2πx).
For the non dissipative boundary condition, one takes:

H = 0 1 1 0 .
One obtains Figure 7 where the energy dynamics are compared for both slope limiter and upwind schemes ( dt dx = 0.4) and also for the exact solution: For the dissipative boundary condition, one takes:

H = 0.2 0.5 0.5 0.2
for which ρ 1 (H) = 0.7. In this case, the evolution of the BV norm is represented in Figure 8 in the same way as in Figure 7. For both dissipative and non dissipative cases, we can make the same conclusions as in section 4.1 in the sense that the upwind scheme dissipates too much when the grid is coarse whereas the slope limiter scheme fits better with the exact solution.

Saturated control

In this section, we analyze an example of system of scalar conservation laws for d = 2 with saturated feedback control law. The same example as [START_REF] Dus | Bv exponential stability for systems of scalar conservation laws using saturated feedback[END_REF] is considered:

   ∂ t R + ∂x[f (R)] = 0 R(t, 0) = HR(t, 1) + Bu(t) R(0, x) = R 0 (x) ∈ BV ([0, 1]) (28)
where:

u(t) = σ(KR(t, 1)).
The operator σ is the saturation by component with level of saturation equal to one:

∀1 ≤ i ≤ 2, x ∈ R, σ i (x) = x if |x| ≤ 1 σ i (x) = sign(x) otherwise.
Matrices are defined as follows.

H = 0 1.1 1 0 , B = I 2 , K = 0 -0.1050 -0.1045 0 .
As nonlinear flux, one takes f (R) = ΛR + 0.2(arctan(R 1 ), arctan(R 2 )) with

Λ = 1 0 0 √ 2 .
For the scheme parameters, it is imposed that dt = 10 -2 and ν = 0.4. Then, we take R 0 constant with value in (-40, 40) 2 and look if the solution does not blow up at infinite time in BV norm. In this way, it is possible to estimate the projection of the basin of attraction onto the space of constant initial data. Results are given below for the upwind and the minmod slope limiter scheme: We see that the upwind scheme overestimates the basin of attraction whereas it is not the case for the minmod slope limiter scheme.

21

The exponential BV stability of general slope limiter scheme for system [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] has been established with an estimate of the dissipation rate. The main idea was to use Harten's formalism. Then, using new discrete Lyapunov techniques, we concluded on the exponential stability of the numerical solution in the spirit of the continuous case. The dissipation rate is estimated from below uniformly with respect to the discretization. Finally, the convergence of the scheme towards the entropy solution was established using the boundedness of ψ dx .

There are some questions which remain open. The case of time dependent flux may not change the philosophy of the proof of exponential stability. However if we add a source term coupling the equations inside the domain the analysis is far from being obvious. Then the case where there are positive and negative characteristic velocities is not treated but this is only a question of change of variable. Finally combining both difficulties of velocities with different signs and the one of source terms, we know that a Lyapunov function does not exist in the continuous case [START_REF] Bastin | On boundary feedback stabilization of non-uniform linear 2 × 2 hyperbolic systems over a bounded interval[END_REF]. It does not seem extravagant to suppose that such Lyapunov function does not exist in the discrete case. Backstepping techniques were designed to solve this difficulty when the flux is linear but it is not obvious that such methods can be directly applied to the corresponding numerical system.

A Proof of Theorem 2

The uniqueness can be proved by classical Kruzhkov techniques [START_REF] Dus | Bv exponential stability for systems of scalar conservation laws using saturated feedback[END_REF]Appendix B]. For the existence and the convergence result, it relies on classical techniques of BV estimates based on the CFL condition [START_REF] Chainais-Hillairet | Second-order finite-volume schemes for a non-linear hyperbolic equation: error estimate[END_REF]. This is why we do not give any proof of these estimates and focus only on the effect of the boundary condition which is less classical. For more details, we refers to [START_REF] Dus | Exponential stability of a general slope limiter scheme for scalar conservation laws subject to a dissipative boundary condition[END_REF].

Let N be a positive integer and denote R ∆x the numerical approximation given in Definition 1. We recall that we defined R n 0 as:

R n 0 = g(R n N ).
A.1 The L ∞ estimate Now, define for all ndt < ν 2 := τ , the index jν (n) such that 2 ν (n -1)dt ≤ jν (n)dx < 2 ν ndt. From now on, we suppose that the time index n verifies ndt < τ . Under a CFL like condition [START_REF] Chainais-Hillairet | Second-order finite-volume schemes for a non-linear hyperbolic equation: error estimate[END_REF], it is easy to prove that:

max ndt<τ jν (n)≤j≤N |R n j | ≤ max 0≤j≤N |R 0 j | (29)
and that:

max ndt<τ 0≤j≤jν (n) |R n j | ≤ max ndt<τ |R n 0 | + max 0≤j≤N |R 0 j |.
Then, the Lipschitz character of the boundary condition allows to transform last inequality

max ndt<τ 0≤j≤jν (n) |R n j | ≤ C(Lg) max ndt<τ |R n N | + max 0≤j≤N |R 0 j |
where C(Lg) is a constant depending only on Lg. Finally, one can easily bound the left boundary term by the L ∞ norm of the initial data to get:

max ndt<τ 0≤j≤jν (n) |R n j | ≤ (1 + C(Lg)) max 0≤j≤N |R 0 j |. (30)
Gathering ( 30) and ( 29), it holds:

R ∆x L ∞ ([0,τ ]×[0,1]) ≤ (1 + C(Lg)) R 0 ∆x L ∞ ([0,1]) .
Repeating this process, we obtain the following L ∞ estimate for any T > 0:

R ∆x L ∞ ([0,T ]×[0,1]) ≤ C(Lg, T ) R 0 ∆x L ∞ ([0,1]) .
where C depends only on T , Lg, ξ and vmax. By the fact that

(R 0 ∆x ) N is bounded in L ∞ (recall that R 0 is in BV ([0, 1])), the L ∞ estimate holds: R ∆x L ∞ ([0,T ]×[0,1]) ≤ C R 0 L ∞ ([0,1]) (31)
where C depends only on T , Lg, ξ and vmax.

A.2 The TV estimate

We have,

T V (R n • ) = T V [0,jν ] (R n • ) + T V [jν ,N ] (R n • ). (32)
Under the CFL condition [START_REF] Chainais-Hillairet | Second-order finite-volume schemes for a non-linear hyperbolic equation: error estimate[END_REF], it is easy to prove that the first term on the left hand side of (32) can be bounded by a boundary term and an initial data term:

T V [0,jν ] (R n • ) ≤ T V [0,n] (R • 0 ) + T V (R 0 • ).
Moreover, by the Lipschitz character of the boundary condition, we get:

T V [0,jν ] (R n • ) ≤ LgT V [0,n] (R • N ) + T V (R 0 • ).
Finally under [START_REF] Chainais-Hillairet | Second-order finite-volume schemes for a non-linear hyperbolic equation: error estimate[END_REF], We can bound T V

[0,n] (R • N ) by T V (R 0 • ) to obtain: T V [0,jν ] (R n • ) ≤ (1 + Lg)T V (R 0 • ). (33)
Similarly, it is possible to prove that:

T V [jν ,N ] (R n • ) ≤ T V (R 0 • ). (34)
and with [START_REF] Zuazua | Propagation, observation, and control of waves approximated by finite difference methods[END_REF], the following T V estimate holds:

∀ndt < τ, T V (R n • ) ≤ C(Lg, T )T V (R 0 • ).
Repeating the process, one gets the following bound in T V for ndt < T with T > 0:

T V (R n • ) ≤ C(Lg, T )T V (R 0 • )
where C(Lg, T ) is a constant that depends solely on Lg and T . By the fact that (R 0 • ) N converges towards R 0 in BV ([0, 1]), one bounds the right hand side of last equation:

∀ndt ≤ T, T V (R n • ) ≤ C(Lg, T, R 0 ). ( 35
) A.3 The L 1 ([0, 1]) continuity estimate For all integers n > 1, N j=0 |R n+1 j -R n j |dx = N j=1 |A n j-1 (R n j -R n j-1 )|dx + |R n+1 0 -R n 0 |dx ≤ N j=1 |A n j-1 ||R n j -R n j-1 |dx + Lg|R n+1 N -R n N |dx = N -1 j=0 |A n j ||R n j+1 -R n j |dx + Lg|A n N -1 (R n N -R n N -1 )|dx ≤ (1 + Lg)T V (R n • )dx.
By the triangle inequality, it holds for n, m > 0

N j=0 |R n+m j -R n j |dx ≤ dx n+m-1 l=n (1 + Lg)T V (R l • ).
And by the boundedness of the total variation (35), we get

N j=0 |R n+m j -R n j |dx ≤ C(Lg, (n + m)dt, R 0 )mdx.
which is no more than:

R ∆x ((n + m)dt, •) -R ∆x (ndt, •) L 1 ([0,1]) ≤ C(Lg, f, (n + m)dt, R 0 )mdt. (36)
A. [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] The L 1 loc (R + ) continuity estimate

We define for T > 0, n T := E(T /dt). For all integers 0 ≤ j ≤ N -1, we have the following estimate:

n T n=0 |R n j+1 -R n j |dt = n T n=0 dt|(a n+1 j ) -1 (R n+1 j+1 -R n j+1 )| ≤ Cdt n T n=0 |R n+1 j+1 -R n j+1 |.
where we have used the fact that the a n i,j s are bounded from below because of the Harten's condition [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF]. By the triangular inequality, it holds for j 0 < j 1 :

n T n=0 |R n j 1 -R n j 0 |dt ≤ Cdt j 1 j=j 0 T V [1,n T +1] (R • j+1 ). ( 37
)
Similarly to proof of the boundedness of the space total variation (see ( 35)), we can prove the corresponding result for the time total variation. That is to say:

∀T > 0, 0 ≤ j ≤ N, T V [0,n T ] (R • j ) ≤ C(T, R 0 , Lg).
Injecting this in (37), the following continuity estimate writes: R ∆x (•, j 1 dx) -R ∆x (•, j 0 dx) L 1 ([0,T ]) ≤ C(T, R 0 , Lg, f )(j 1 -j 0 )dx. (38)

A.5 Conclusion

By Helly's Theorem, the sequence (R ∆x ) N converges (up to a subsequence) in L 1 loc (R + ×[0, 1]) to a function R ∈ L 1 loc (R + , BV ([0, 1])) ∩ Lip loc (R + , L 1 ([0, 1])). It remains to prove that R is an entropy solution to [START_REF] Banda | Numerical discretization of stabilization problems with boundary controls for systems of hyperbolic conservation laws[END_REF].

Let k ∈ R d and define q k (u) := f (u ⊥ k) -f (k ⊥ u), η k (u) = |u -k|. The proof of the following result crucially depends on [START_REF] Eymard | Solution of Equation in R n (Part 3)[END_REF] and is given in a more general context in [START_REF] Chainais-Hillairet | Second-order finite-volume schemes for a non-linear hyperbolic equation: error estimate[END_REF]. We give it here for completeness.

Lemma 4 (Discrete entropy estimate)

For all k ∈ R d , N ∈ N * , we have for all n ∈ N, 1 ≤ j ≤ N, 1 ≤ i ≤ d:

1 dt (|R n+1 i,j -k i | -|R n i,j -k i |) + 1 dx (q k (R n i,j + Rn i,j ) -q k (R n i,j-1 + Rn i,j-1 )) ≤ Cdx 1-2β . ( 39
)
In order to ease the reading, the i index is dropped until the end of the proof. We take ϕ ∈ Cc(]0, T [×]0, 1[). Multiplying (39) by I n j (ϕ) :=

t n+1
t n

x j +dx/2

x j -dx/2 ϕ(t, x)dtdx and summing over n, j, one gets:

n T n=0 N j=1 η k (R n+1 j ) -η k (R n j ) dt I n j (ϕ)+ n T n=0 N j=1
q k (R n j + Rn j ) -q k (R n j-1 + Rn j-1 ) dx I n j (ϕ) ≤ Cdx 1-2β .

A numerical integration by parts (owing the fact that ϕ is zero at the boundary of [0, T ]×[0, 1]), we get:

n T n=0 N j=1 η k (R n j ) I n-1 j (ϕ) -I n j (ϕ) dt + n T n=0 N j=1
q k (R n j + Rn j )

I n j (ϕ) -I n j+1 (ϕ) dx ≥ Cdx 1-2β .
Passing to the limit owing that β < 1/2, one gets T 0 1 0 η k (R)∂tϕ + q k (R)∂xϕdxdt ≥ 0.

A.5.1 Boundary conditions and convergence estimates

To finish the proof of the existence of a solution to [START_REF] Banda | Numerical discretization of stabilization problems with boundary controls for systems of hyperbolic conservation laws[END_REF], it remains to show that the initial and boundary conditions are verified in the almost everywhere sense. Let ε > 0 and 0 < s < 1 ||R(0, .) -R ∆x (0, . By the fact that R 0 ∆x converges towards R 0 in L 1 ([0, 1]), we deduce that R(0, .) = R 0 in a L 1 sense and R(0, .) = R 0 almost everywhere.

Remark 5

The same procedure can be repeated for any t ≥ 0 and prove [START_REF] Kružkov | First order quasilinear equations with several independent variables[END_REF].

Remark 6 We can use the same technique to prove that the boundary condition R(t, 0) = g(R(t, 1)) is satisfied in the almost everywhere sense. In this case, the continuity estimate (38) is the key tool (see [13, p.24] for details).

The final step is to prove the convergence estimate [START_REF] Kröner | Convergence of higher order upwind finite volume schemes on unstructured grids for scalar conservation laws in several space dimensions[END_REF]. We give the sketch of the proof (see [13, p.24-25] for more details). In fact ( 21) is a consequence of Helly's Theorem applied to the sequence (R ∆x ) N > Ñ where Ñ > 1, on a time interval [t, t + dt] (dt > 0). Then it suffices to take the limit when Ñ goes to infinity and when dt goes to zero.

Notation:

  For all R ∈ R d , |R| designates the canonical euclidean norm of R. For matrices M ∈ M d (R), |M | = sup |M R| |R|=1, R∈R d . For 1 ≤ p ≤ ∞ and all matrices M ∈ M d (R), |M |p is subordinate norm of l p . D + d (R)is the set of diagonal strictly positive matrices. The real ρp(M ) for matrices M ∈ M d (R) is defined by: ρp(M ) := inf ∆∈D + d (R) |∆M ∆ -1 |p.

  (a) dt = 10 -2 (b) dt = 5 × 10 -4

Fig. 6 :

 6 Fig. 6: Trace of the solution

3 Fig. 7 :

 37 Fig. 7: Non dissipative case: comparison with the exact solution

(a) dt = 10 - 2 (b) dt = 10 - 3 Fig. 8 :

 1021038 Fig. 8: Dissipative case: comparison with the exact solution

Fig. 9 :

 9 Fig. 9: Comparison of basins of attraction (dt = 10 -2 )

Fig. 10 :

 10 Fig. 10: Basin of attraction for dt = 10 -3

t t 0 ||Rt 1 0 1 0

 011 )|| L 1 ([0,1]) ≤ ||R(0, .) -R(s, .)|| L 1 ([0,1]) + ||R(s, .) -R ∆x (s, .)|| L 1 ([0,1]) +||R ∆x (s, .) -R ∆x (0, .)|| L 1 ([0,1]) ≤ 2C(g, R 0 )s + ||R(s, .) -R ∆x (s, .)|| L 1 ([0,1]) .where we have used (36) to get the last inequality and the fact thatR ∈ Lip loc (R + , L 1 ([0, 1])).Integrating with respect to s on an interval [0, t] for 0 ≤ t ≤ 1, one gets||R(0, .) -R ∆x (0, .)|| L 1 ([0,1]) ≤ C(g, R 0 )t + 1 (s, .) -R ∆x (s, .)|| L 1 ([0,1]) ds ≤ C(g, R 0 )t + 1 ||R(s, .) -R ∆x (s, .)|| L 1 ([0,1]) ds Taking t = ε 2C(g,R 0 ) and N sufficiently large such that ||R(s, .) -R(s, .)|| L 1 ([0,1]) ds ≤ ε 2 4C(g, R 0 ), one finally obtains:||R(0, .) -R ∆x (0, .)|| L 1 ([0,1]) ≤ ε.