Olivier Bokanowski

Nicolas Forcadel
email: nicolas.forcadel@ensta.fr

Hasnaa Zidani
email: hasnaa.zidani@ensta.fr

L 1 -error estimates for numerical approximations of Hamilton-Jacobi-Bellman equations in dimension 1

Keywords: Hamilton-Jacobi-Bellman equations, Fast Marching Method, Ultra-Bee scheme, L 1 -error estimate, anti-diffusive scheme Mathematics Subject Classification: 49L99, 65M15

The goal of this paper is to study some numerical approximations of particular Hamilton-Jacobi-Bellman equations in dimension 1 and with possibly discontinuous initial data. We investigate two anti-diffusive numerical schemes, the first one is based on the Ultra-Bee scheme and the second one is based on the Fast Marching Method. We prove the convergence and derive L 1 -error estimates for both schemes. We also provide numerical examples to validate their accuracy in solving smooth and discontinuous solutions.

Introduction

This paper discuss two explicit numerical approximations of the following one-space dimensional Hamilton-Jacobi-Bellman (HJB) equation

ϑ t + max (f + (x)ϑ x , f -(x)ϑ x ) = 0 in R × (0, T ) ϑ(•, 0) = v 0 in R (1.1)
where v 0 ∈ L ∞ (R). In particular, v 0 can be discontinuous.

In optimal control theory, the solution ϑ of equation (1.1) corresponds to the value function of an optimization problem [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF][START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]. It often happens that this function, as well as the "final" cost v 0 , is discontinuous (for instance for target or Rendez-Vous problems). The discontinuities of ϑ will represent, for example, the interface between the domain of admissible trajectories and the one of prohibited trajectories and then it is very important to localize the discontinuities. This is the reason why, in the discontinuous case, the classical monotone schemes for HJB equations ([8, 11, 1, 13]) are no more adapted. Indeed, if we attempt to use these schemes, we observe an increasing numerical diffusion around discontinuities, and this is due to the fact that monotone schemes use at some level finite differences and/or interpolation technics.

In this work, we investigate two different schemes to solve (1.1) for discontinuous initial data. The first one is a Fast Marching Method type scheme. This method, introduced by Sethian [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF], is a very efficient scheme to solve numerically the eikonal equation for given positive velocity c(x). This scheme has been improved by Carlini et al in [START_REF] Carlini | Convergence of a generalized fast marching method for a non-convex eikonal equation[END_REF][START_REF] Forcadel | Comparison principle for a generalized fast marching method[END_REF] where the case of changing sign velocity is considered.

Recall that the FMM method was built [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Forcadel | Comparison principle for a generalized fast marching method[END_REF] to deal with eikonal equations with initial condition taking values in {0, 1}, this schemes allows to concentrate numerical effort only in a Narrow band around the interface separating the zone of 0-value from the zone of 1-value.

Here we consider the case where the initial condition v 0 is any bounded lower semi continuous (l.s.c) function. We first define a level-set approximation w 0 of v 0 in the following form: given p ≥ 1 and given (h k ) k=1,••• ,p ⊂ R * + , we set

w 0 (x) = p k=1
h k w 0,k (x), with w 0,k (x) := 1 if v 0 (x) > p i=1 h i 0 otherwise. Now for each level k = 1, • • • , p, the function w 0,k takes values only in {0, 1}. Therefore, we propose an algorithm based on the FMM which make evolve each level-set function w 0,k , for k = 1, • • • , p. Hence, we obtain for every k an approximation ϑ ρ k (with ρ := (∆x, (h k ) k )) of the solution of (1.1) associated to the initial condition w 0,k . Thanks to a comparison principle, we prove that an approximation of the solution ϑ of (1.1) is obtained by:

ϑ ρ (t, x) = p k=1 h k ϑ ρ k (t, x).
We derive an L 1 -error estimate (in finite time t) in order of ∆x, that is the mesh step size.

Let us mention that Y. Brenier [START_REF] Brenier | Averaged multivalued solutions for scalar conservation laws[END_REF] has used a similar level decomposition in the case of conservation laws.

The second scheme is a modified version of the Ultra-Bee scheme for HJB equations proposed in [START_REF] Bokanowski | Anti-diffusive schemes for linear advection and application to Hamilton-Jacobi-Bellman equations[END_REF] and which convergence has been proved in [4]. Let us mention that this scheme was first studied in [START_REF] Després | A non-linear anti-diffusive scheme for the linear advection equation[END_REF][START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF] for linear advection equations with constant velocity. In this case, it is proved that the scheme is exact whenever the initial function takes values in {0, 1} and the discontinuities are separated by 3∆x (∆x being the mesh step size). The scheme keeps nice anti-diffusive properties when we deal with advection or HJB equations with changing sign velocities and initial condition taking only values 0 -1. A generalisation of the Ultrabee scheme is also proposed in [4] for HJB equation with l.s.c bounded initial condition v 0 . This generalisation use additional steps of truncation and prediction when two discontinuities get close (closer than 3∆x).

In this paper, we use the level-set decomposition of v 0 (as explained above). We lead back to an HJB equation in the form of (1.1) with initial condition w 0,k which takes only values 0 and 1. The evolution of each level-set function can be accurately approximated by Ultrabee scheme, and the resulting approximation of the solution ϑ of (1.1) is very satisfactory. The Ultrabee scheme combined with level-set decomposition has almost the same L 1 -error bound than the Ultrabee scheme studied in [4], but numerically it seems that the method proposed in this paper give more accurate results (see Section 5, for a numerical comparison).

The aim of future work is to use the level set decomposition idea in order to treat twodimensional (or more) HJB evolution equations of the form u t + max α (f α (x) • ∇u) = 0 by related schemes.

The paper is organized as follow. In Section 2 we present our main results: a scheme based on the Fast Marching Method (FMM), an Ultra-Bee scheme (UB), and main convergence results for both schemes in an L 1 -error approximation bound. Section 3 is devoted to some preliminary results. Next Section deals with the convergence proof for the FMM. Numerical simulations are finally presented in Section 5. Some technical proofs are postponed to the Appendix.

Main results

In this section, we present the convergence results for the FMM scheme and for the Ultra-Bee scheme.

The principal assumptions on the dynamics are the following ones:

(H1) f + and f -are L-Lipschitz continuous. (H2) ∃ε > 0 ∀x ∈ R, f -(x) + ε ≤ f + (x).
Remark 2.1. This last assumption will allow us to compare the velocities f + and f -on different nearby points. It can be replaced by

(H2') f -(x) ≤ f + (x)
, ∀x ∈ R, and f + and f -are non-decreasing functions on R.

or by

(H2") f -(x) ≤ 0 ≤ f + (x), ∀x ∈ R.
On the initial condition v 0 , we assume that (H3) v 0 ∈ L ∞ (R), v 0 is lower semi-continuous, and has a finite number of extrema, in the following sense:

      
There exist real numbers A 1 , . . . , A q+1 and B 1 , . . . , B q with

A 1 = -∞ ≤ B 1 < A 2 < • • • < B q ≤ A q+1 = +∞, (with possibly B 1 = -∞ or B q = +∞), such that v 0 ր on each [A i , B i [, v 0 ց on each ]B i , A i+1 ], and v 0 (B i ) = min(v 0 (B - i ), v 0 (B + i )).
(A i are local minima of v 0 , and B i are local maxima of v 0 ). We also assume that v 0 ≥ 0 and is compactly supported: ∃α, β ≥ 0,

supp(v 0 ) ⊂ [α, β]. (2.2)
We finally assume that v 0 is locally Lipschitz continuous in a neighbourhood of its local minima:

(H4) ∃δ 0 > 0, ∀i = 2, . . . q, v 0 is Lipschitz continuous in [A i -δ 0 , A i + δ 0 ]
Let us recall the definition of total variation of a real-valued function. Let ∆x > 0 be a step size of a spatial grid, and let x j := j ∆x denote a uniform mesh, with j ∈ Z. Consider also

Level-set decomposition

x j+ 1 2 := (j + 1 2
) ∆x, and

I j :=]x j-1 2 , x j+ 1 2 [.
We define w 0 , a level set decomposition of v 0 , and the function w 0,k of level k, by:

w 0,k (x i ) := 1 hk <v 0 (x i ) = 1 if hk < v 0 (x i ) 0 otherwise, (2.4a) w 0,k (x) = w 0,k (x i ) for x ∈ I i . (2.4b)
We also set

w 0 (x) := k=1,...,p h k w 0,k (x), x ∈ R. (2.4c) Remark 2.2. Definition (2.4) clearly implies that w k (t, x) ∈ {0, 1}. Moreover, if 1 ≤ k 1 ≤ k 2 ≤
p, then from the comparison principle [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF] and the fact that

w 0,k 1 (x) ≤ w 0,k 2 (x), ∀x ∈ R, we obtain w k 1 (t, x) ≤ w k 2 (t, x), ∀t ≥ 0 and x ∈ R.
Now, the idea is to propose two algorithms to compute numerically the approximation ϑ ρ k (where ρ = (∆x, h) represents the space discretization) of the solution w k of (1.1) with initial data w 0,k . This two scheme are based on the Fast Marching Method and on the Ultra-Bee Method.

As soon as we have computed the numerical solution ϑ ρ k , a natural approximation of the solution ϑ of (1.1) is simply given by ϑ ρ = p k=1 h k ϑ ρ k (see Proposition 3.3). We now describe in details the two algorithms as well as the convergence result we obtain.

First, we give an error approximation estimate between v 0 and its projection w 0 . The proof is left to the reader.

Lemma 2.2. (Error at initial time)

We have the following estimate

w 0 -v 0 L 1 (R) ≤ (β -α)h + T V (v 0 )∆x, ∀x ∈ R.
(2.5)

Next, we will compare the evolution of the viscosity solutions of (1.1) associated to initial data v 0 and w 0 . Proposition 2.3. Assume (H1)-(H2)-(H3) with ∆x ≤ δ 0 . Let w k (resp. w) be the viscosity solution of (1.1) with initial data w 0,k (resp. w 0 ). Then

(i) w(t, x) = k=1,...,p h k w k (t, x), ∀t > 0, x ∈ R, (ii) w(t, .) -ϑ(t, .) L 1 (R) ≤ e Lt (β -α + M 0 t) h + e Lt T V (v 0 ) + M 1 ∆x.
where

M 0 = M 0 (v 0 , f ) := i=2,...,q (|f + (A i )| + |f -(A i )|) and M 1 = max j,∃i, A i ∈ Īj v ′ 0 L ∞ (I j )
are constant.

Remark 2.3. In this proposition indeed only f -(x) ≤ f + (x) is needed, not (H2).

Fast Marching Method

The idea is to make evolve each level set w 0,k using an adaptation of the Generalized Fast Marching Method introduced in [START_REF] Carlini | Convergence of a generalized fast marching method for a non-convex eikonal equation[END_REF] (see also [START_REF] Forcadel | Comparison principle for a generalized fast marching method[END_REF]).

Notations and algorithm

Let ∆x > 0 be a mesh step size of a uniform grid. For k = 1, • • • , p and i ∈ Z, we consider:

θ 0,k i := 2w 0,k (x i ) -1 = 1 if v 0 (x i ) > hk -1 if v 0 (x i ) ≤ hk (2.6)
(the introduction of θ is just useful to formulate the algorithm in a simple way and in particular to have some symmetry properties of the algorithm).

As in [START_REF] Carlini | Convergence of a generalized fast marching method for a non-convex eikonal equation[END_REF], we also define approximated piece-wise constant velocity functions f + and f -, for

x ∈]x j-1 2 , x j+ 1 2
[ and for j ∈ Z:

f α (x) := 0 if ∃ i ∈ {j ± 1} s.t.f α (x i )f α (x j ) ≤ 0 and |f α (x j )| ≤ |f α (x i )|, f α (x j ) otherwise.
(

This "regularization" allows us to introduce a numerical band of zero to separate the region where the velocity is positive from the one where the velocity is negative. This separation is needed to avoid the duplication of the front (see [START_REF] Carlini | Convergence of a generalized fast marching method for a non-convex eikonal equation[END_REF]).

Let us remark that, from (1.1), when θ 0,k i = -1 = -θ 0,k i+1 , then the discontinuity evolves with the velocity f + (to the right if f + > 0 and to the left if f + < 0), while when θ 0,k i = 1 = -θ 0,k i+1 , the discontinuity evolves with the velocity f -(see Figure 1).

We now define, for each control α ∈ {-, +}, the stencil of grid points useful to compute the value at point

x i , for i ∈ Z U n,k α (i) :=              i + 1 if θ n,k i = -θ n,k i+1 = -α1 and f α (x i ) < 0 i -1 if θ n,k i-1 = -θ n,k i = -α1 and f α (x i ) > 0 ∅ otherwise and U n,k α = i U n,k α (i).
The set U n,k will play the role of the frozen points of the classical Fast Marching Method. We point out that the set U n,k α (i) is either empty or a singleton. We also define a set of Narrow Bands by:

NB n,k α := i, U n,k α (i) = ∅ , NB n,k := NB n,k + ∪ NB n,k - and NB n := k=1,..,p NB n,k . i -1 ∈ U + (i) i ∈ NB + i ∈ NB + i + 1 ∈ U + (i)
Figure 1: Representation of the useful points and of the Narrow Band for the control α = +. Left: f + (x i ) > 0 and the discontinuity moves to the right. Right: f + (x i ) < 0 and the discontinuity moves to the left.

We describe now our FMM for Hamilton-Jacobi-Bellman equation (1.1). This is an adaptation of the one proposed in Carlini et al. [START_REF] Carlini | Convergence of a generalized fast marching method for a non-convex eikonal equation[END_REF][START_REF] Forcadel | Comparison principle for a generalized fast marching method[END_REF]. In order to track correctly the evolution we need to introduce a discrete function τ n,k i,α ∈ R + , defined only for the points i ∈ U n,k α , to represent the approximated physical time for the front propagation at the nodes i for the level set k, the control α and at the n-th iteration of the algorithm (FMM).

The idea of the algorithm is then very simple. For each point i of the Narrow Band N B, we compute a tentative value τi of the arrival time of the front, using the time of the useful points. We then find the minimum of the τi and we accept the nodes that realize the minimum (i.e. we change the value of the θ) and we iterate the process. Let us now give our algorithm in details.

Initialization: for n = 0, initialize the field θ 0,k as in (2.6), and set

τ 0,k i,α := 0 if i ∈ U 0,k α , +∞ otherwise, for k = 1, • • • , p and α = ±.
Loop: for n ≥ 1, 1. Compute τ n-1,k on NB n-1,k as follows: for α = ±, define

τ n-1,k i,α :=        τ n-1,k ı,α + ∆x | f α (x i )| if i ∈ NB n-1,k α , +∞ otherwise,
where ı ∈ U n-1,k α (i), and set

τ n-1,k i := min α∈{+,-} τ n-1,k i,α . 2. Set t n := min τ n-1,k i , i ∈ NB n-1,k , k ∈ {1, ..., p} .
3. Define the new accepted point

NA n,k = {i ∈ NB n-1,k , τ n-1,k i = t n }.
4. Update the values of θ n,k :

θ n,k i = -θ n-1,k i if i ∈ NA n,k θ n-1,k i otherwise 5. Reinitialize τ n,k i,α on U n,k α : τ n,k i,α = min(t n , τ n-1,k i,α ) if i ∈ U n,k α +∞ otherwise 6. If t n ≥ T then stop. Else, set n := n + 1.
Remark 2.4. Let us remark that in our algorithm, the minimum time t n is taken on all the level sets. This allows us in particular to have a comparison principle between the level sets (see Corollary 2.5).

Main results for the FMM scheme

We extend the function θ n,k in the following way

θ ρ,k (t, x) := θ n,k i if x ∈ [x i-1 2 , x i+ 1 2 + ∆x) and t ∈ [t n , t n+1 ) (2.8)
where ρ denotes (∆x, h). Hence, we define a function ϑ ρ by:

ϑ ρ (t, x) := p k=1 θ ρ,k (t, x) + 1 2 h k .
(2.9)

First we shall check that ϑ ρ well define a numerical approximation of the solution ϑ of (1.1). This claim is a consequence of a comparison principles for the numerical level-set functions (θ ρ,k ).

Theorem 2.4. (Discrete comparison principle)

Let 1 < k 1 < k 2 ≤ p.
For all n ∈ N and all i ∈ Z, we have either

θ n,k 1 i > θ n,k 2 i or θ n,k 1 i = θ n,k 2 i =: σ i = ±1 and if i ∈ U n,k 1 α ∩ U n,k 2 α , then τ n,k 1 i,α ≤ τ n,k 2 i,α if σ i = +1 τ n,k 1 i,α ≥ τ n,k 2 i,α if σ i = -1
The proof of this theorem is technical and is given in Appendix B. A first straightforward consequence of this discrete comparison principle can be formulated as follows:

Corollary 2.5. (Comparison principle for the level-set functions)

Let 1 ≤ k 1 < k 2 ≤ p. Then θ ρ,k 2 ≤ θ ρ,k 1 .
Now, we can give the statement of the main result of this section. (The proof will be done in Section 4.) Theorem 2.6. (Convergence of the FMM scheme) Assume (H1)-(H4), and let ρ = (∆x, h) with ∆x < min( ε L , δ 0 ) and h as in (2.3). The numerical solution ϑ ρ given by the FMM scheme, defined as in (2.9), converges to the viscosity solution ϑ of (1.1), and for t ≥ 0, the following error estimate holds:

ϑ ρ (t, .) -ϑ(t, .) L 1 (R) ≤ 5 2 e Lt + 3Lte Lt T V (v 0 ) + M 1 ∆x + e Lt (β -α + 2M 0 t)h, (2.10)
where

M 0 = M 0 (v 0 , f ) := i=2,...,q (|f + (A i )| + |f -(A i )|) and M 1 = max j,∃i, A i ∈ Īj v ′ 0 L ∞ (I j )
are constant.

Remark 2.5. Furthermore if h is chosen to be of the order of ∆x (for instance using h k ≡ h := ∆x, ∀k), we deduce a global estimate of order ∆x in L 1 -error.

Remark 2.6. In the level set decomposition, we can choose (h j ) j such that v 0 (A i ) = w 0 (A i ) for i = 2, . . . , q. In this case, assumption (H4) is not needed (see the proof of Proposition 2.3 in the appendix).

Remark 2.7. When the velocities f + and f -depend on time, it is possible to adapt the algorithm as in [START_REF] Forcadel | Comparison principle for a generalized fast marching method[END_REF] and to obtain the comparison principle and the convergence result (in the same way as [START_REF] Carlini | Convergence of a generalized fast marching method for a non-convex eikonal equation[END_REF][START_REF] Forcadel | Comparison principle for a generalized fast marching method[END_REF]). Nevertheless, we are not able, in this case, to prove the L 1 error estimate.

Ultra-Bee scheme

Algorithm (UB)

Let ∆t > 0 be a constant time step, and t n := n∆t for n ≥ 0. Let us notice that in the FMM approach, each iteration takes into account the evolution of all levels-set functions w k . On the contrary, the UB scheme should be performed starting from each w 0,k independently of the others. This scheme aims to compute, for every k = 1, • • • , p, a numerical approximation of the averages w n,k j := 1 ∆x I j w k (t n , x) dx, for j ∈ Z. Since the function w k (t n , •) takes only values in {0, 1}, their averages w n,k j contain the information of the discontinuities localization. The UB scheme gives an accurate approximation of (w n,k j ) j as long as the discontinuities are separated by more than 2∆x. Otherwise, when two discontinuities are sufficiently close, a truncation step is made in order to avoid numerical diffusion around these discontinuities. The scheme takes the following form.

V n+1,1 j -V n j ∆t + max α=±   f α (x j ) V n,L j+ 1 2 ,α -V n,R j-1 2 ,α ∆x   = 0, (2.11a) V n+1 j = [Trunc(V n+1,1 )] j , (2.11b)
with the initialization:

V 0 j := 1 ∆x I j w 0,k (x)dx.
(2.12)

Here V n,L j+ 1 2 ,α and V n,R j+ 1
We first set, for j ∈ Z and α = ±,

ν j,α := ∆t ∆x f α (x j ),
the "local CFL" number. We assume that |ν j,α | ≤ 1, ∀j ∈ Z, and for α = ±.

(2.13)

We also consider

if ν j,α > 0,      b + j,α := max(V n j , V n j-1 ) + 1 ν j,α (V n j -max(V n j , V n j-1 )), B + j,α := min(V n j , V n j-1 ) + 1 ν j,α (V n j -min(V n j , V n j-1 )), (2.14a) if ν j,α < 0,      b - j,α := max(V n j , V n j+1 ) + 1 |ν j,α | (V n j -max(V n j , V n j+1 )), B - j,α := min(V n j , V n j+1 ) + 1 |ν j,α | (V n j -min(V n j , V n j+1 )).
(2.14b)

Under condition (2.13), these numbers satisfy b

+ j,α ≤ B + j,α , b - j,α ≤ B - j,α
, and correspond to flux limiters that ensure stability properties. Now, we define the UB scheme as follows (see [START_REF] Bokanowski | Anti-diffusive schemes for linear advection and application to Hamilton-Jacobi-Bellman equations[END_REF]4]).

UB Algorithm: For each level-set function w 0,k , k = 1, • • • , p, we consider the following evolution algorithm.

Initialization: We compute the initial averages (V 0 j ) j∈Z as in (2.12).

Loop:

For n ≥ 0, We compute V n+1 = (V n+1 j
) j∈Z by: A) Evolution. Define "fluxes" V n j+ 1 2 ,α for α ∈ {-, +} as follows:

• If ν j,α ≥ 0, set V n,L j+ 1 2 ,α :=    min(max(V n j+1 , b + j,α ), B + j,α ) if ν j,α > 0 V n j+1 if ν j,α = 0 and V n j = V n j-1 V n j if ν j,α = 0 and V n j = V n j-1 , • If ν j,α ≤ 0, set V n,R j-1 2 ,α :=    min(max(V n j-1 , b - j,α ), B - j,α ) if ν j < 0 V n j-1 if ν j,α = 0 and V n j = V n j+1 V n j if ν j,α = 0 and V n j = V n j+1 , (where b + j,α , b - j,α , B + j,α and B - j,α are defined by (2.14a)-(2.14b)). • If ν j,α ≥ 0 and ν j+1,α > 0, set V n,R j+ 1 2 ,α := V n,L j+ 1 2 ,α . • If ν j+1,α ≤ 0 and ν j,α < 0, set V n,L j+ 1 2 ,α := V n,R j+ 1 2 ,α . • If ν j,α < 0 and ν j+1,α > 0, then set V n,R j+ 1 2 ,α := V n j+1 if V n j+1 = V n j+2 V n j otherwise and V n,L j+ 1 2 ,α := V n j if V n j = V n j-1 V n j+1 otherwise. (2.15) Set V n+1,1 j := min {α=±} V n j -ν j,α V n,L j+ 1 2 ,α -V n,R j-1 2 ,α .
B) Truncation: Set V n+1 := Trunc(V n+1,1 ) as follows.

• For all indexes j such that

     max(V n+1,1 j-1 , V n+1,1 j+1 ) < 1, and V n+1,1 j = 1 or 0 < max(V n+1,1 j , V n+1,1 j+1 ) < 1, and V n+1,1 j-1 = V n+1,1 j+2 = 0, set V n+1 j-1 = V n+1 j = V n+1 j+1 := 0. • Otherwise set V n+1 j := V n+1,1 j .
For every k = 1, • • • , p, we associate to the scheme values (V n j ) j , the l.s.c. step function ϑ ρ k defined for every t ≥ 0, x ∈ R by

ϑ ρ k (t, x) := V n j if x ∈]x j-1 2 , x j+ 1 2 [, t ∈ [t n , t n+1 [, min(V n j , V n j+1 ) if x = x j+ 1 2 , t ∈ [t n , t n+1 [. (2.16)
The UB scheme approximation of the solution ϑ of (1.1) is finally determined by

ϑ ρ (t, x) := p k=1 h k ϑ ρ k (t, x).
(2.17) Remark 2.8. A general version of the Ultra-bee scheme is given in [4], for any l.s.c. initial condition in L 1 loc (R). Here, the algorithm (and specially the truncation step) is specified to the case of an initial condition taking values only in {0, 1}. Also, in the algorithm of [4] there is a prediction step that is unnecessary in our context.

Several stability and convergence properties of this scheme can be found in [START_REF] Bokanowski | Anti-diffusive schemes for linear advection and application to Hamilton-Jacobi-Bellman equations[END_REF]4].

Convergence Result for the Ultra-Bee scheme

For every k = 1, • • • , p, the function ϑ ρ k (given by (2.16)) corresponds to a step wise function of level k. We have the following L 1 error estimate between the solution ϑ of (1.1) and its numerical approximation ϑ ρ given by the Ultra-Bee scheme. (Proof is postponed to the end of Section 3.) Theorem 2.7. (Error estimate) Assume (H1)-(H2)-(H3). We assume that ∆x ≤ min( ε 2L , δ 0 ). Let ϑ ρ be as in (2.17) defined by the algorithm UB. Then for all t n ≥ 0, we have the following estimate:

ϑ ρ (t n , .) -ϑ(t n , .) L 1 (R) ≤ (1 + e Ltn + Lt n e Ltn )T V (v 0 ) + M 1 ∆x + e Ltn (β -α + 2M 0 t n )h,
where M 0 and M 1 are the same constant as in Theorem 2.6.

Proof of Theorem 2.7: ¿From [4, Theorem 4], for every k = 1, • • • , p we have

ϑ ρ k (t n , .) -w k (t n , .) L 1 (R) ≤ (Lt n e Ltn + 1)T V (w 0,k )∆x.
Summing for k = 1, . . . , p, we obtain

ϑ ρ (t n , .) -w(t n , .) L 1 (R) ≤ (Lt n e Ltn + 1)T V (v 0 )∆x
where we have used that

p k=1 h k T V (w 0,k ) = T V (w 0 ) ≤ T V (v 0 ) (see for instance [4, Lemma B.3]).
Together with Proposition 2.3, we get the desired estimate.

Preliminary results

From now on, consider an initial condition v 0 satisfying (H3), ϑ solution of (1.1) with the initial condition v 0 , and the level set decomposition w 0 of v 0 (w 0 defined as in (2.4)).

Mesh approximation. First, we define exact and approximated characteristics that will be useful throughout the paper. As the dynamics f -and f + are Lipschitz continuous, then for any a ∈ R we can define characteristics X a,+ and X a,-as the solutions of the following Cauchy problems:

Ẋa,+ (t) = f + (X a,+ (t)), X a,+ (0) = a, and Ẋa,-(t) = f -(X a,-(t)). X a,-(0) = a. (3.18)
In general, the differential equation

χa (t) = f + (χ a (t)), a.e. t ≥ 0, χ a (0) = a, (3.19)
may have more than one absolutely continuous solution. The non-uniqueness comes from the behaviour on boundary points (x j+ 1 2

) in the case when the velocity vanishes (or changes sign).

Throughout this paper, we shall denote by X S a,+ the function defined by: X S a,+ is an absolutely continuous solution of (3.19), and

(3.20) if ∃t * ≥ 0, ∃j ∈ Z s.t. X S a,+ (t * ) = x j+ 1 2 , f + (x j )f + (x j+1 ) ≤ 0 then X S a,+ (t) = x j+ 1 2 ∀t ≥ t * .
We have uniqueness of such solution (see [4, Appendix A]). We construct X S a,-in a similar way. By using the same arguments as in [4, Lemma 1 and Lemma 9], we get: Lemma 3.1. Assume that (H1) and (H2) hold. Let a, b be in R. The following assertions are satisfied:

(i) Let s ≤ t and assume that X S b,-(θ) ≥ X S a,+ (θ), for every θ ∈ [s, t]. Then |X S b,-(t) -X S a,+ (t)| + 2∆x ≤ e L(t-s) (|X S b,-(s) -X S a,+ (s)| + 2∆x). (ii) If a ≥ b + ∆x and ∆x ≤ ε L , then X S a,+ (t) ≥ X S b,-(t) + ∆x for every t ≥ 0.
Also, we have the following representation of the solutions

w k , k = 1, • • • , p. Lemma 3.2. (Lemma 2 of [4])
Assume that (H1)-(H2) hold.Then the unique viscosity solution of (1.1) with initial condition w 0,k is given by:

w k (t, x) = min y∈[X x,+ (-t),X x,-(-t)] w 0,k (y), ∀t > 0, x ∈ R. (3.21)
We also consider the function w S k which is defined in an analogous way as in (3.21), but with the approximated characteristics X S x,+ , X S x,-instead of X x,+ and X x,-:

w S k (t, x) := min y∈[X S x,+ (-t),X S x,-(-t)] w 0,k (y), ∀t > 0, x ∈ R, (3.22)
The approximate function w S k will play an important role throughout the paper: the two studied schemes will give approximation of the function w S k . By using the same arguments as in [4, Proposition 1], we have the following L 1 -error estimate:

w S k (t, .) -w k (t, .) L 1 (R) ≤ 3Lte Lt T V (w 0,k ) ∆x. (3.23)
The aim now is to approximate w S k by some discrete scheme.

k = 1, • • • , p, an approximation ϑ ρ k of w S k such that ϑ ρ k (t, .) -w S k (t, .) L 1 (R) ≤ C t T V (w 0,k ) ∆x, (3.24)
for some constant C t ≥ 0. Define a global approximation by

ϑ ρ := k=1,...,p h k ϑ ρ k .
Then we have the estimate

ϑ ρ (t, .) -ϑ(t, .) L 1 (R) ≤ (C t + e Lt + 3Lte Lt )T V (v 0 ) + M 1 ∆x + e Lt (β -α + 2M 0 t)h, (3.25)
where M 0 and M 1 are defined as in Theorem 2.6.

Proof: Set w S (t, .) := k=1,...,p h k w S k (t, .). By summing the estimate (3.24) for k = 1, . . . , p, we obtain

ϑ ρ (t, .) -w S (t, .) L 1 (R) ≤ C t p k=1 h k T V (w 0,k )∆x ≤ C t T V (v 0 )∆x,
where we have used that

p k=1 h k T V (w 0,k ) = T V (w 0 ) ≤ T V (v 0 ) (see for instance [4, Lemma B.3]).
On the other hand, by using (3.23), we obtain

w S (t, .) -w(t, .) L 1 (R) ≤ 3Lte Lt T V (v 0 ) ∆x.
We conclude the proof by combining Proposition 2.3 and the previous bounds.

Fast Marching Method 4.1 General remarks on the algorithm

Proposition 4.1. The following properties hold:

(i) For all i ∈ U n,k α , we have τ n,k i,α ≤ t n (ii) If i ∈ NA n,k , then τ n,k i,α = t n if i ∈ U n,k α +∞ otherwise (iii) If i ∈ U n-1,k α ∩ U n,k α , then τ n,k i,α = τ n-1,k i,α . (iv) If i ∈ U n,k α \U n-1,k α , then τ n,k i,α = t n .
Proof (i) This is a straightforward consequence of Step 5 of the algorithm.

(ii) By Step 5 of the algorithm, we just have to prove

that if i ∈ NA n ∩U n,k α , then min(τ n-1,k α , t n ) = t n .
Let us consider the case when α = + and i ∈ U n,k + (i + 1) (the other cases being similar). Thus, f + (x i+1 ) > 0, θ n,k i = -1 and θ n,k i+1 = 1. Also, since i ∈ N A n , then we have:

θ n-1,k i = 1, and i ∈ N B n-1,k . (4.26) Assume that min(τ n-1,k + , t n ) = τ n-1,k

+

. This implies in particular that τ n-1,k + < ∞ and so i ∈ U n-1,k + . We claim that: i ∈ NB n-1,k -, and i + 1 ∈ NA n . Indeed, (4.26) implies that i ∈ U n-1,k + (i -1) and f + (x i-1 ) < 0 which gives that i ∈ NB n-1,k

+ . Since i ∈ N B n-1,k , we get that i ∈ N B n-1,k - . Moreover, the fact that θ n-1,k i = 1 and i ∈ NB n-1,k - implies that i + 1 ∈ U n-1,k - (i). It yields that f -(x i ) < 0 and θ n-1,k i+1 = -1. Since θ n-1,k i+1 = -1 = -θ n,k
i+1 , we deduce that i + 1 ∈ NA n . On the other hand, the fact that

θ n-1,k i+1 = -1 and f -(x i ) implies that i + 1 ∈ NB n-1 -. But i+1 ∈ NA n and so i+1 ∈ NB n-1,k + . This implies that f + (x i+1 ) < 0, which leads to a contradiction. (iii) If i ∈ U n-1,k α ∩ U n,k
α , then by Step 5 of the algorithm, we have

τ n,k i,α = min(τ n-1,k i,α , t n ) = τ n-1,k i,α . (iv) If i ∈ U n,k α \U n-1,k α
, then by Step 5 of the algorithm, we have

τ n,k i,α = min(τ n-1,k i,α , t n ) = t n
where we have used that

τ n-1,k i,α = +∞ since i ∈ U n-1,k α .
4.2 Proof of Theorem 2.6

In view of Proposition 3.3 we mainly have to deal with the case p = 1 (one level set).

For simplicity of notation we denote

θ n ≡ θ n,1 . Let (a i ) i=1,...,d , (b i ) i=1,...,d ⊂ (Z + 1 2 )∆x be such that a 1 < b 1 < a 2 < b 2 < ... < a d < b d .
We consider an initial data of the form

w 0 (x) := d i=1 1 ]a i ,b i [ (x). (4.27)
We have the following convergence result for one level set which proof is given in the following subsection:

Proposition 4.2. (Convergence of the FMM scheme for one level) Assume (H1)-(H2) and p = 1. Let ρ = ∆x ≤ ε L . We have the following error estimate between the numerical solution ϑ ρ given by the FMM scheme, defined as in (2.9), and the approximate solution w S defined in (3.22)

ϑ ρ (t, .) -w S (t, .) L 1 (R) ≤ 3 2 e Lt T V (v 0 )∆x. (4.28)
The proof of this proposition is given in the following subsection.

Proof of Theorem 2.6 Theorem 2.6 is a straightforward consequence of Propositions 3.3 and 4.2.

Proof of Proposition 4.2

Before to give the proof of Proposition 4.2, we have to give some definitions and preliminary results.

Notations. We define the numerical discontinuities in the following way: Let a ∈ (Z + 1 2 )∆x and i

0 a ∈ Z such that a = x i 0 a -1 2 . Assume that θ 0 i 0 a = α1 = -θ 0 i 0
a -1 for some α ∈ {+, -}. In particular, the numerical discontinuity starting from a will move with the velocity f α . We denote by (t n,∆ a ) n the sequence of time at which the numerical discontinuity starting from a moves and by X S,∆ a,α (t n,∆ a ) the position of the discontinuity at time t n,∆ a . More precisely, we define:

t 0,∆ a = 0 and X S,∆ a,α = a = x i 0 a -1 2 and for n ≥ 1,                          t n,∆ a = inf t m > t n-1,∆ a s.t. i n-1 a ∈ NA m or i n-1 a -1 ∈ NA m , i n a =      i n-1 a + 1 if f α (x i n-1 a ) > 0 and i n-1 a ∈ NA m , i n-1 a -1 if f α (x i n-1 a -∆x) < 0 and i n-1 a -1 ∈ NA m , i n-1 a otherwise X S,∆ a,α (t n,∆ a ) = x i n a -1 2 .
(4.29)

We now define the extinction time T ∆ a of the numerical discontinuity starting from a by

T ∆ a = inf{t n,∆ a , θ m i n a = θ m i n a -1 for m such that t m = t n,∆ a } Remark 4.1.
In the definition of i n a (4.29), we have i n a = i n-1 a only if the discontinuity X S,∆ a,α do not move and disappear.

Lemma 4.3. (Profile of the discontinuity) Let i 0 a ∈ N be such that θ 0

i 0 a -1 = α1 = -θ 0 i 0 a with α ∈ {+, -}. Let n ∈ N be such that t n,∆ a < T ∆ a . Then, for all m ∈ N such that t n,∆ n ≤ t m < t n+1,∆ a , we have θ m i n a = -θ m i n a -1 = α1.
This lemma claims in fact that the numerical discontinuity starting from a and evolving with the velocity f α is located at a discontinuity of θ m and always keep the same profile. Proof. Let us assume that α = + (the case α =being similar). By recurrence, let us assume that

θ l i n-1 a = -θ l i n-1 a -1 = α1 for all l ∈ N such that t n-1,∆ a ≤ t l < t n,∆ a . Let us define m such that t m = t n,∆
a . We have

i n-1 a ∈ NA m or i n-1 a -1 ∈ NA m .
Step 1:

i n a = i n-1 a
By contradiction, assume that i n a = i n-1 a . Let us assume that i n-1 a ∈ NA m (the other case being similar). This implies in particular that f

+ (x i n-1 a ) ≤ 0. Moreover, we have θ m-1 i n-1 a = 1 = -θ m-1 i n-1 a -1 . This implies that i n-1 a ∈ NB m-1 + . But i n-1 a ∈ NA m and so i m-1 a ∈ NB m-1 - which implies that f -(x i n-1 a ) < 0.
(4.30)

Since i n-1 a ∈ NA m , we also deduce that θ m i n a = θ m i n-1 a = -1. But t n,∆ a < T ∆ a and so θ m i n-1 a -1 = θ m i n a -1 = 1. This implies that i n-1 a -1 ∈ NA m . Since i n a = i n-1 a , we then deduce that f + (x i n-1 a -1 ) ≥ 0 and so i n-1 a -1 ∈ NB m-1 + . This implies that i n-1 a -1 ∈ NB m-1 - and so f -(x i n-1 a -1 ) > 0. This contradicts (4.30).
Step 2:

θ m i n a = -θ m i n a -1 = 1 Let us assume that i n a = i n-1 a + 1 (the case i n a = i n-1 a -1 being similar). This implies that i n-1 a ∈ NA m . But θ m-1 i n-1 a = 1 and so θ m i n a -1 = θ m i n-1 a = -1. Moreover, since t n,∆ a < T ∆ a , we deduce that θ m i n a = -θ m i n a -1 = 1.
Step 3: conclusion By definition of t n+1,∆ a , for all l such that t n,∆ a ≤ t l < t n+1,∆ a , we have i n a ∈ NA l and i n a -1 ∈ NA l and so θ l )

Let i 0 a ∈ N be such that θ 0 i 0 a = α1 = -θ 0 i 0 a -1 with α ∈ {+, -}. Let n ∈ N. Assume that t n,∆ a < T ∆ a . Then, for all m ∈ N such that t n,∆ n ≤ t m < t n+1,∆ a , we have τ m i n a ,α = t n,∆ a if f α (x i n a -∆x) < 0 τ m i n a -1,α = t n,∆ a if f α (x i n a ) > 0 Remark 4.2. If f α (x i n α -∆x) < 0, since by Lemma 4.3 we have θ m i n a = -θ m i n a -1 , we deduce that i n a ∈ U m α (i n a -1)
and so τ m i n a ,α is well defined. In the same way, if f α (x i n a ) > 0, then i n a -1 ∈ U m α (i n a ) and τ m i n a -1,α is well defined. Proof of Lemma 4.4 We assume that α = + and that f α (x i n a -∆x) < 0 (the others cases being similes). Let m * be such that t m * = t n,∆ a . The proof is decomposed into two steps:

Step 1 τ m * i n a ,+ = t n,∆ a . If n = 0, then m * = 0 and τ 0 i n a = 0 = t 0 = t n,∆ a . Let us treat the case n ≥ 1. We claim that i n a = i n-1 a -1. Indeed, if i n a = i n-1 a + 1, then, by (4.29), we have f + (x i n a -∆x) = f + (x i n-1 a
) > 0, which is absurd. By (4.29), we then deduce that

i n a = i n-1 a -1 ∈ NA m * . This implies that τ m * i n a ,+ = t m * = t n,∆ a .
Step 2

τ m i n a ,+ = t n,∆ a for all m such that t n,∆ a < t m < t n+1,∆ a . By Lemma 4.3, for all m such that t n,∆ a < t m < t n+1,∆ a , we have θ m i n a = -θ m i n a -1 = 1. Since f + (x i n a -1 ) < 0, this implies that i n a ∈ U m + (i n a -1) for all m such that t n,∆ a < t m < t n+1,∆ a
. By Proposition 4.1 (iii), we then get that

τ m i n a ,+ = τ m * i n a ,+ = t n,∆
a . This ends the proof of the lemma.

For a ∈ (Z + 1 2 )∆x, we also define the time when the discontinuity X S a,+ changes of mesh by

t 0 a = 0 ∀n ≥ 0, t n+1 a = inf{t > t n a , |X S a,+ (t) -X S a,+ (t n a )| = ∆x}. For b ∈ (Z + 1
2 )∆x, we define the time when the discontinuity X S b,-changes of mesh in the same way. We define the extinction time of the discontinuity a i and b i by

T a i = T b i := inf{t, X S a i ,+ (t) = X S b i ,-(t)}.
The following proposition claims that the FMM scheme computes exactly the position of the discontinuity before it disappear. Proposition 4.5. (Exact computation of the discontinuity before the meeting of discontinuity) Let a ∈ (Z + 1 2 )∆x and i 0 a ∈ N be such that a = x i 0 a -1

2

. We assume that θ 0

i 0 a = α1 = -θ 0 i 0 a -1 with α ∈ {+, -}. For all n ∈ N, if t n,∆ a < T ∆ a , then t n a = t n,∆ a X S a,α (t n a ) = X S,∆ a,α (t n,∆ a ).
Proof of Proposition 4.5

Let us assume that α = +, the case α =being similar. By recurrence, assume that

t n-1 a = t n-1,∆ a X S a,+ (t n-1 a ) = X S,∆ a,+ (t n-1,∆ a
).

Since t n,∆ a < ∞, we have either

f + (X S a,+ (t n-1 a ) -∆x) < 0 or f + (X S a,+ (t n-1 a )) > 0.
Let us assume that f + (X S a,+ (t n-1 a ) -∆x) < 0 (the other case being similar). This implies that the discontinuities X S a,+ (t n-1 a ) and X S,∆ a,+ (t n-1 a ) will move to the left. For the approached discontinuity, we have

t n a = t n-1 a + ∆x | f + (X S a,+ (t n-1 a ) -∆x)| and X S a,+ (t n a ) = X S a,+ (t n-1 a ) -∆x.
We now turn to the numerical discontinuity. We define m such that t m = t n,∆ a . In particular, since

f + (x i n-1 a -1 ) < 0, we have i n a = i n-1 a -1 ∈ NA m .
Let us now compute t m . For simplicity of notation, let us denote i n-1 a -1 by i. By the algorithm, we have

τ m-1 i,+ = τ m-1 i n a ,+ + ∆x | f + (x i )| .
We now claim that i ∈ NB n-1 -. By contradiction, assume that i ∈ NB n-1 -. By Lemma 4.3, we have θ m-1 i = -1. This implies that θ m-1 i-1 = 1 and f -(x i ) > 0. This gives in particular that i -1 ∈ NB m-1 -. Moreover, since f + (x i-1 ) ≤ 0, we also have deduce that i -1 ∈ NB m-1 + and so i -1 ∈ NA m which gives

θ m i-1 = θ m-1 i-1 = 1.
Since i ∈ NA m , we also deduce that θ m i = 1. This contradicts the fact that t n,∆ a < T ∆ a and proves that i ∈ NB n-1 -. We then deduce that τ n-1 i,-= +∞ and so

t m = τ m-1 i = τ m-1 i,+ .
By Lemma 4.4, we have τ m-1

i n a ,+ = t n-1,∆ a . So we recover that t n,∆ a = t m = t n a . Moreover, we have i n a = i n-1 a -1 and X S,∆ a (t n,∆ a ) = x i n a -1 2 = x i n-1 a -1 2 -∆x = X S,∆ a,+ (t n-1,∆ a ) -∆x = X S a,+ (t n-1 a ) -∆x = X S a,+ (t n a ).
This ends the proof of the proposition.

The following proposition claims that the discontinuities X S,∆ a i ,+ (t) and X S,∆ b i-1 ,-(t) cannot meet. This essentially comes from our assumption (H2). )∆x be such that

a 1 < b 1 < a 2 < b 2 < ... < a d < b d .
Assume that

θ 0 i = 1 if ∃ j s.t a j < x i < b j -1 if ∃ j s.t b j < x i < a j+1
Then, for all t ≥ 0, i ∈ {1, ..., d -1}, we have

X S,∆ a i ,+ (t) -X S,∆ b i-1 ,-(t) ≥ ∆x.
Proof of Proposition 4.6 By contradiction, let us define

t * = inf{t ≥ 0, ∃i ∈ {1, ..., d} s.t. X S,∆ a i ,+ (t) -X S,∆ b i-1 ,-(t) < ∆x}.
We denote by ı the index such that X S,∆ a ı ,+ (t * ) -X S,∆ b ı-1 ,-(t * ) < ∆x and by a = a ı , b = b ı-1 . Let us define n and m such that

t n+1,∆ a ≤ t * < t n+2,∆ a t m+1,∆ b ≤ t * < t m+2,∆ b
In particular, we have either

t * = t n+1,∆ a or t * = t m+1,∆ b . Finally we define i ∈ Z such that x i = X S,∆ a,+ (t * ) + 1 2 .
The proof is decomposed into three cases:

Case 1: t n+1,∆ a = t m+1,∆ b = t *
In this case, the two discontinuities have moved at time t * and we have (since we have X S,∆ a, In the first case, we then have f -

+ (t n,∆ a )- X S,∆ b,-(t m,∆ b ) ≥ ∆x) i n+1 a = i m b = i m+1 b -1 = i n a -1 or i m b + 1 = i m+1 b = i n+1 a = i n a -1. X S,∆ a,+ (t n,∆ a ) = X S,∆ b,-(t m+1,∆ b ) = X S,∆ b,-(t m,∆ b ) X S,∆ a,+ (t n+1,∆ a ) x i ∆x x i-1 = X S,∆ b,-(t m+1,∆ b ) X S,∆ a,+ (t n,∆ a ) x i X S,∆ a,+ (t n+1,∆ a ) X S,∆ b,-(t m,∆ b ) ∆x
(x i ) = f -(x i m b ) > 0 and f + (x i ) = f + (x i n a -∆x) < 0 which contradicts the fact that f + ≥ f -.
In the other case, we have

0 < f -(x i m b ) = f -(x i-1 ) = f -(x i-1 ) ≤ f + (x i-1 ) -ε ≤ f + (x i ) = f + (x i ) = f + (x i n a ) < 0
where we have used assumption (H2) and the fact that f α (x j ) = 0 implies that f α (x j ) = f α (x j ). This is absurd.

Case 2: t m+1,∆ b < t n+1,∆ a = t * In this case, only the discontinuity X S,∆ a,+ has moved at time t * and we have (since

X S,∆ a,+ (t n,∆ a ) - X S,∆ b,-(t m,∆ b ) ≥ ∆x) i m+1 b = i n+1 a = i n a -1.
We then deduce that

f + (x i ) = f + (x i ) = f + (x i n+1 a
) < 0 and so f -(x i ) < 0. We define k such that

X S,∆ a,+ (t n,∆ a ) = X S,∆ b,-(t m+1,∆ b ) X S,∆ a,+ (t n+1,∆ a ) x i ∆x
Figure 3: Representation of the discontinuity before the meeting.

t k = t n+1,∆ a = t * .
Step 1:

Ordering τ k-1 i,-≤ τ k-1 i+1,+ . By Lemma 4.4, since f + (x i ) = f + (x i n+1 a ) = f + (x i n a -1 ) < 0, we get τ k-1 i+1,+ = τ k-1 i n a ,+ = t n,∆ a .
Moreover, using assumption (H2), we deduce that

f -(x i m+1 b -1 ) = f -(x i-1 ) ≤ f + (x i-1 ) -ε ≤ f + (x i ) < 0.
We now claim that f -(x i-1 ) = f -(x i-1 ). By contradiction, if it is not the case then (since

f -(x i-1 ) < 0 and f -(x i ) < 0) f -(x i-2 ) > 0 and |f -(x i-1 )| ≤ |f -(x i-2 )|.
But, by assumption (H2), we also have

f -(x i-1 ) ≤ f + (x i ) < 0 and 0 < f -(x i-2 ) ≤ f + (x i-1 )
and so

|f + (x i-1 )| ≥ |f -(x i-2 )| ≥ |f -(x i-1 )| ≥ |f + (x i )|.
This implies that f + (x i ) = 0, which is absurd. This implies that f -(x i-1 ) = f -(x i-1 ) < 0. By Lemma 4.4, we then get

τ k-1 i,-= τ k-1 i m+1 b ,-= t m+1,∆ b .
To get the result, we then have to prove that t m+1,∆ b ≤ t n,∆ a . By contradiction, assume that t m+1,∆ b > t n,∆ a . We then have

t n,∆ a < t m+1,∆ b < t n+1,∆ a < t m+2,∆ b . Moreover, at time t = max(t m,∆ b , t n,∆ a ) < t * , we have X S,∆ a,+ ( t) = X S,∆ a,+ (t n,∆ a ) = x i n a -1 2 = x i+ 1 2 , X S,∆ b,+ ( t) = X S,∆ a,+ (t m,∆ b ) = x i m b -1 2 = x i+ 1 2
which contradicts the definition of t * .

Step 2: contradiction We have

t n+1,∆ a = t k = τ k-1 i,+ =τ k-1 i+1,+ + ∆x | f + (x i )| ≥ τ k-1 i,-+ ∆x | f + (x i-1 )| ≥τ k-1 i,-+ ∆x | f -(x i-1 )| = τ k-1 i-1,- ≥t k .
We then deduce that τ k-1 i-1,-= t k and so i m+1

b

-1 = i-1 ∈ NA k This implies that t m+2,∆ b = t k = t * . This is absurd. Case 3: t n+1,∆ a < t m+1,∆ b = t *
This case can be treated in the same way as Case 2.

The following lemma claims that when the discontinuities X S,∆ a,+ (t), X S a,+ (t), X S,∆ b,-(t) and X S b,-move in the same direction, then X S,∆ a,+ (t) and X S,∆ b,-(t) meet before X S a,+ (t) and X S b,-. This comes from the fact that the numerical discontinuities cannot be in the same mesh.

(Z + 1 2 )∆x be such that a < b. Let θ 0 = 1 [a,b[ . Assume that f + (a)f -(b) ≥ 0. Then T ∆ a = T ∆ b ≤ T a = T b .
Proof of Lemma 4.7 By contradiction, assume that T ∆ a = T ∆ b > T a . We assume that f + (a) ≥ 0 and f -(b) ≥ 0 (the other case is similar). This implies that the two discontinuities will move to the right. Let us define n, m such that

t n a < T a ≤ t n+1 a t m b < T a ≤ t m+1 b
In particular, since the discontinuities move to the right, we have for all t

∈ [t n,∆ a , t n+1,∆ a [ X S,∆ a,+ (t n a ) = X S a,+ (t n a ) ≤ X S a,+ (t) < X S a,+ (t n a ) + ∆x and X S,∆ b,-(t m b ) = X S b,-(t m b ) ≤ X S b,-(t) < X S b,-(t m b ) + ∆x Let us define i = i n a and j = i n b .
In particular, we have

X S a,+ (t) ∈ I i and X S b,-(t) ∈ I j for t ∈ [t n a , t n+1 a [. Since X S,∆ a,+ (t * ) = X S,∆ b,-(t * ) for t * = sup(t n a , t n b ), we deduce that i < j. This implies in particular that t n+1 a = T a = T b . Moreover, since X S a,+ (t) ∈ I i , X S b,-(t) ∈ I j for t ∈ [t n a , t n+1 a [, i < j and X S a,+ (t n+1 a ) = X S b,-(t n+1 a
), we deduce that i = j -1 and

X S a,+ (t n a ) + ∆x = X S a,+ (T a ) = X S b,-(t n b ) = X S b,-(t * ) = x i+ 1 2 .
This implies that the discontinuity X S b,-(t n b ) do not move and then f -(x i+1 ) = 0. We then deduce that X S,∆ b,-(t) = x i+ 1 2 for all t ≥ t n b . Moreover, by definition of t n+1 a , we have

T a = t n+1 a = t n a + ∆x b f + (x i ) = t n+1,∆ a . So at time t n+1,∆ a , we have X S,∆ a,+ (t n+1,∆ a ) = X S,∆ b,-(t n+1,∆ a ) = x i+ 1 2 and so T ∆ a = T ∆ b = t n+1,∆ a
= T a which is absurd.

Proposition 4.8. (Error estimate after the meeting of discontinuities)

Let a, b ∈ (Z + 1 2 )∆x be such that a < b. Let θ 0 = 1 [a,b[ . For all t ∈ [inf(T a , T ∆ a ), sup(T a , T ∆ a )[, we have R 1 [X S a i ,+ (t),X S b i ,-(t)] (x) -1 [X S,∆ a i ,+ (t),X S,∆ b i ,-(t)] (x) dx ≤ 3∆xe Lt . (4.31) Proof of Proposition 4.8 If f + (a)f -(b) ≥ 0, then (4.31
) is a consequence of Lemma 4.7 and Lemma 3.1 (i).

If f + (a) ≤ 0 and f -(b) ≥ 0, then the discontinuities X S a,+ and X S,∆ a,+ will move to the left while the discontinuities X S b,-and X S,∆ b,-will move to the right. We then deduce that T a = T ∆ a = ∞ and so the result is trivial.

Let us then assume that f + (a) ≥ 0 and f -(b) ≤ 0. Then the discontinuities X S a,+ and X S,∆ a,+ will move to the right while the discontinuities X S b,-and X S,∆ b,-will move to the left. Let us assume that T a < T ∆ a (the other case being similar) and let us define n, m such that

t n a < T a ≤ t n+1 a t m b < T a ≤ t m+1 b .
We recall that, by Proposition 4.5, we have

X S a,+ (t n a ) = X S,∆ a,+ (t n a ) = x i n a = x i and X S b,-(t m b ) = X S,∆ b,-(t m b ) = x i m b = x j . Moreover, we have ∀t ∈ [sup(t n a , t m b ), T a ) X S a,+ (t) ∈ I i and X S b,- (t)
∈ I j and j = i + 1. We then deduce that

X S,∆ b,-(t * ) -X S,∆ a,+ (t * ) ≤ ∆x for t * = sup(t n a , t m b
). Since X S,∆ b,-moves to the left and X S,∆ a,+ moves to the right, we deduce (4.31).

We are now able to give the proof of Proposition 4.2:

Proof of Proposition 4.2 (continued). Recall that we assumed w 0 of the form (4.27). By definition of w S and ϑ ρ , we have

w S (x, t) = d i=1 1 [X S a i ,+ (t),X S b i ,-(t)] (x) and ϑ ρ (x, t) = d i=1 1 [X S,∆ a i ,+ (t),X S,∆ b i ,-(t)] (x)
We then deduce that

w S (•, t) -ϑ ρ (•, t) L 1 (R) ≤ d i=1 R 1 [X S a i ,+ (t),X S b i ,-(t)] (x) -1 [X S,∆ a i ,+ (t),X S,∆ b i ,-(t)] (x) dx
For each i ∈ {1, ..d}, we then have to estimate

I i = R 1 [X S a i ,+ (t),X S b i ,-(t)] (x) -1 [X S,∆ a i ,+ (t),X S,∆ b i ,-(t)] (x) dx
We distinguish three cases:

Case 1: t ≥ sup(T a i , T ∆ a i ) In this case, we have [X S a i ,+ (t), X S b i ,-(t)] = ∅ and [X S,∆ a i ,+ (t), X S,∆ b i ,-(t)] = ∅
and so I 1 = 0.

Case 2: inf(T a i , T ∆ a i ) ≤ t < sup(T a i , T ∆ a i )
In this case, by Proposition 4.8, we have I 1 ≤ 3∆xe Lt .

Case 3: t < inf(T a i , T ∆ a i ) In this case, we have

I 1 = R 1 [X S a i ,+ (t),X S b i ,-(t)]∆[X S,∆ a i ,+ (t),X S,∆ b i ,-(t)] (x)dx ≤|X S a i ,+ (t) -X S,∆ a i ,+ (t)| + |X S b i ,-(t) -X S,∆ b i ,-(t)| ≤2∆x
where A∆B is the symmetric difference of the sets A and B and where we have used Proposition 4.5 for the last line.

We then deduce that we always have

I i ≤ 3e Lt ∆x and so w S (t, •) -ϑ ρ (t, •) L 1 (R) ≤ 3de Lt ∆x = 3 2 e Lt T V (v 0 )∆x since T V (v 0 ) = 2d.

Numerical Simulations

In all the following examples, we consider an equation on (x min , x max ) in the form of:

ϑ t + max(f -(x)ϑ x , f + ϑ x ) = 0 t ≥ 0, x ∈ (x min , x max ), ϑ(0, x) = v 0 (x); x ∈ (x min , x max ), (5.32)
with periodic boundary conditions. We will denote by N x the number of mesh points considered in (x min , x max ), and by p the number of levels used in the level-set decomposition of v 0 (see

(2.3)-(2.4a)).
Example 1: We first consider an advection equation with constant velocity, f -= f + ≡ -1, on (x min , x max ) = (-2, 2) with periodic boundary conditions, with v 0 defined as follows:

v 0 (x) := 0.64 if x ∈ [0.2, 0.6], max(1 -x 2 , 0) otherwise
We show in figure 4, the numerical solution compared to the exact solution, with parameters N x = p = 50 (and CFL number 0.75 for the Ultra-bee Scheme). The solution is periodic of period T = 4 and we show the solution at time t = 12 (3 periods, thus the exact solution recovers its initial position). For this example, we remark a very good behaviour of the both schemes even for long time. The L 1 error produced by FMM comes only from the level-set decomposition of v 0 . Then the advection of each level-set function is exact (constant velocity). The Ultra-Bee scheme provides also a very good solution, and the L 1 error corresponding to this algorithm comes from the decomposition of v 0 and also from the truncation made around the maxima.

In Table 1 we show the L 1 -error for the two schemes and see that they are both of first order. We consider now the case of f -= -1 and f + = -0.9. The domain is (x min , x max ) = (-2, 2), and the initial data is given by v 0 (x) := max(max(0, 1 -|x|), max(0, 0.7 -|x -0.2|)).

N x = p FMM
We compare the Ultra-Bee scheme, the FMM method and the exact solution in Fig. In this example, we consider an advection equation on the domain (x min , x max ) = (0, 1) with periodic boundary conditions, and with variable velocity

f -(x) = f + (x) = -1 -0.5 sin(2πx),
The initial data is given by v 0 (x) := max(1 -16(x -0.5). The solution can be shown to be periodic of period T = 2/ √ 3.

Example 4: In this example, we consider the case of variable velocity functions:

f -(x) = -sin(πx) -0.5, f + (x) = -sin(πx),
on the domain (x min , x max ) = (-1, 1) with periodic boundary conditions, and with the initial data: v 0 (x) := | sin(π(x -0.5)/2)|.

The classical Ultra-Bee scheme (see [4]) tends to project the solution on a class of step functions. The Ultrabee scheme combined with a level-set decomposition does not show diffusion. It furthermore gives a good approximation that is not amplified for longer times.

A Proof of Proposition 2.3

We first give a stability result by deriving an L 1 -error estimate between the solution ϑ of (1.1), corresponding to the initial data v 0 , and the solution u of (1.1) associated to another initial data u 0 .

Proposition A.1. Assume (H1) and (H2). Let u 0 and v 0 be two real valued l.s.c. functions such that v 0u 0 ∈ L 1 (R), v 0 satisfies (H3), and let u (resp. ϑ) be the l.s.c. solutions of (1.1) with initial data u 0 (resp. v 0 ). We suppose that for all interval I ⊂ R, v 0 ր on I ⇒ u 0 ր on I, v 0 ց on I ⇒ u 0 ց on I;

(1.33)

Then, for all t ≥ 0, we have ϑ(t, .)u(t, .) L 1 (R) ≤ e Lt v 0u 0 L 1 (R) + M 0 t e Lt max i=2,...,q |v 0 (A i )u 0 (A i )|. (1.34) where

M 0 = M 0 (v 0 , f ) := i=2,...,q |f + (A i )| + |f -(A i )|.
Remark A.1. Indeed in Proposition A.1, instead of (H2) it suffices that f -(x) ≤ f + (x).

Proof of Proposition A.1 We recall that the case when ∀i = 2, . . . , q, v 0 (A i ) = u 0 (A i ) has already been proved in [4] (in view of assumption (1.33), this case means that the values of v 0 and u 0 are the same at local minima of v 0 ). In order to treat the general case, we start with the following lemma.

Lemma A.2. Assume (H1)-(H2)-(H3). Let v0 (x) := v 0 (x) if x / ∈ (A i ) i=2,...,q and v0 (A i ) := v 0 (A i )δ i , with δ i ≥ 0, for i = 2, . . . , q. Let θ be the viscosity solution associated to (1.1) with initial data v0 . Then

|| θ(t, .) -ϑ(t, .)|| L 1 (R) ≤ M 0 t e Lt max i δ i . Proof. Let I x (t) := [X x,+ (-t), X x,-(-t)]. We first have θ(t, x) := min y∈Ix(t) v0 (y) = min min y∈Ix(t) v 0 (y), min i,A i ∈Ix(t) (v 0 (A i ) -δ i ) and ϑ(t, x) := min y∈Ix(t) v 0 (y) = min min y∈Ix(t) v 0 (y), min i,A i ∈Ix(t) v 0 (A i ) .
Hence by difference,

| θ(t, x) -ϑ(t, x)| ≤ max i, A i ∈Ix(t) |(v 0 (A i ) -δ i ) -v 0 (A i )| ≤ max i, A i ∈Ix(t) δ i .
In the case ∀i,

A i / ∈ I x (t) then θ(t, x) = ϑ(t, x). Hence || θ(t, .) -ϑ(t, .)|| L 1 (R) ≤ (max i δ i ) i L {x, A i ∈ I x (t)} ≤ (max i δ i ) i |X A i ,+ (t) -X A i ,- (t)|.
(1.35) where L is Lebesgue's measure. On the other hand, if

X a,+ (0) = X a,-(0) = A i , then |X a,+ (t) -A i | = | t 0 f + (X a,+ (s) ds| ≤ t 0 (|f + (X a,+ (s)) -f + (A i )| + |f + (A i )|) ds ≤ t|f + (A i )| + t 0 L|X a,+ (s) -A i | ds.
Using a Gronwall Lemma, we obtain

|X a,+ (t) -A i | ≤ |f + (A i )| e Lt -1 L ≤ |f + (A i )| te Lt .
(where we have used that (e x -1)/x ≤ e x with x = Lt). In the same way,

|X a,-(t) -A i | ≤ |f -(A i )| te Lt .
Hence the desired result.

We come back to the proof of Proposition A.1. We consider now δ u i := (u 0 (A i )v 0 (A i )) + , and δ v i := (v 0 (A i )u 0 (A i )) + . As in the previous Lemma, we define ū0 by ū0 (x) = u 0 (x) if x / ∈ (A i ), and ū0 (A i ) := u 0 (A i )δ u i . In the same way we consider v0 such that v0 (x) = v 0 (x) if x / ∈ (A i ) and v0 (A i ) := v 0 (A i )δ v i . We finally denote by ū and θ the solution of (1.1) with initial condition ū0 and v0 . Then we have

ϑ -u L 1 (R) ≤ θ -ϑ L 1 (R) + θ -ū L 1 (R) + ū -u L 1 (R) . Then θ(t, .) -ū(t, .) L 1 (R) ≤ e Lt v0 -ū0 L 1 (R) ≤ e Lt v 0 -u 0 L 1 (R) since v0 (A i ) = ū0 (A i ) ∀i = 2, .
. . , q and using the result from [4, Proposition 3]. Also both terms θϑ L 1 (R) and ūu L 1 (R) are controlled by Lemma A.2:

ū(t, .) -u(t, .) L 1 (R) + θ(t, .) -ϑ(t, .) L 1 (R) ≤ M 0 t e Lt max i (δ u i , δ v i ) Since max(δ u i , δ v i ) ≤ |v 0 (A i ) -u 0 (A i )| we conclude to the desired result. Proof of Proposition 2.3: (i)
We proceed by recursion on the number p in the level set decomposition of w 0 .

First we notice that using Lemma 3.2 the viscosity solution of (1.1) with initial data λw 0,1 (for a given λ > 0) is given by λw 1 . This proves the result when p = 1. Now we assume that p ≥ 2 and that the result is true for up to p -1 levels. Let w

(1) 0 (x) := p-1 k=1 h k w 0,k (x) and w

(2) 0 (x) := h p w 0,p (x). We denote by w (1) (resp w (2) ) the viscosity solution of (1.1) with initial data w 0 . We want to prove that w ≡ w (1) + w (2) . Using the representation of Lemma 3.2 we have for a given t ≥ 0 and a given x, w(t, x) = inf y∈I w  Let a ∈ I be such that w(t, x) = w

0 (a). We shall prove that w

(A i ) -w 0 (A i )|. Let j be such that A i ∈ I j . We then have |v 0 (A i ) -w 0 (A i )| ≤ |v 0 (A i ) -v 0 (x j )| + |v 0 (x j ) -w 0 (x j )| ≤ ||v ′ 0 || L ∞ (I j ) ∆x + h (1.36)
The result is then a direct consequence of Proposition A.1 and estimates (2.5) and (1.36).

B Proof of the discrete comparison principle

Proposition B.1. (Comparison principle on the times) Let 1 < k 1 < k 2 ≤ p. For all n ∈ N and all i ∈ Z, we have either

θ n,k 1 i > θ n,k 2 i or θ n,k 1 i = θ n,k 2 i =: σ i = ±1 and if i ∈ U n,k 1 α ∩ U n,k 2 α , then τ n,k 1 i,α ≤ τ n,k 2 i,α if σ i = +1 τ n,k 1 i,α ≥ τ n,k 2 i,α if σ i = -1
Proof of Proposition B.1 By contradiction, let n be the first index such that the condition does not hold and we denote by i a node where the condition is not fullfilled. In particular, we have

     θ n,k 1 i < θ n,k 2 i or θ n,k 1 i = θ n,k 2 i = σ i , i ∈ U n,k 1 α ∩ U n,k 2 α
and σ i τ n,k 1 i,α > σ i τ n,k 2 i,α

(2.37)

The proof is decomposed in several cases: Case 1: i ∈ NA n,k 1 \ NA n,k 2 Since the condition is true at step n -1, we have in particular θ n-1,k 1 i ≥ θ n-1,k 2 i . Sub-case 1.1: θ n-1,k 1 i > θ n-1,k 2 i Then θ n-1,k 1 i = 1 and θ n-1,k 2 i = -1. Since i ∈ NA n,k 1 \ NA n,k 2 , we deduce that θ n,k 1 i = θ n,k 2 i = -1 = σ i . We then have τ n,k 2 i,α ≤ t n = τ n,k 1 i,α

where we have used the fact that i ∈ U n,k 2 α and Proposition 4.1(i) for the first inequality and the fact that i ∈ NA n,k 1 ∩ U n,k 1 α joint to Proposition 4.1(ii) for the last equality. This contradicts (2.37).

Sub-case 1.2:

θ n-1,k 1 i = θ n-1,k 2 i Sub-case 1.2.1: θ n-1,k 1 i = θ n-1,k 2 i = -1 Since i ∈ NA n,k 1 \ NA n,k 2 ,
we deduce that θ n,k 1 i = 1 and θ n,k 2 i = -1. This contradicts (2.37).

Sub-case 1.2.1:

θ n-1,k 1 i = θ n-1,k 2 i = 1 Since i ∈ NA n,k 1 \ NA n,k 2 , we have θ n,k 1 i = -1 and θ n-1,k 1 i = θ n-1,k 2 i = θ n,k 2 i = 1. Let ᾱ be such that τ n-1,k 1 i, ᾱ = τ n-1,k 1 i = t n . Let ı ∈ U n-1,k 1 ᾱ (i)
. This implies that θ n-1,k 1 ı = -1. Since θ n-1,k 1 ≥ θ n-1,k 2 , we get that θ n-1,k 2 ı = -1 and so ı ∈ U n-1,k 2 ᾱ (i). We then deduce that τ n-1,k 1 ı,ᾱ ≥ τ n-1,k 2 ı,ᾱ .

Using Step 1 of the algorithm, we deduce that

t n = τ n-1,k 1 i = τ n-1,k 1 i, ᾱ ≥ τ n-1,k 2 i, ᾱ ≥ τ n-1,k 2 i .
Contradiction since i ∈ NA n,k 2 .

Case 2: i ∈ NA n,k 2 \ NA n,k 1 This case can be treated in the same way of Case 1.

Case 3: i ∈ NA n,k 1 ∩ NA n,k 2 Sub-case 3.1:

θ n-1,k 1 i = θ n-1,k 2 i
Since i ∈ NA n,k 1 ∩ NA n,k 2 , we deduce that θ n,k

1 i = θ n,k 2 i
, and so by (2.37), i ∈ U n,k 1 α ∩ U n,k 2 α . By Proposition 4.1(ii), we deduce that τ n,k 1 i,α = τ n,k 2 i,α = t n .

Contradiction.

Sub-case 3.2:

θ n-1,k 1 i > θ n-1,k 2 i
In this case, we have θ n-1,k 1 = -1. This is contradicts the fact that θ n-1,k 2 i ≤ θ n-1,k 1 i .

Case 4: i ∈ NA n,k 1 ∪ NA n,k 2 Sub-case 4.1: θ n,k 1 i < θ n,k 2 i Since i ∈ NA n,k 1 ∪ NA n,k 2 , we have

θ n-1,k 1 i = θ n,k 1 i < θ n,k 2 i = θ n-1,k 2 i .

This is absurd.

Sub-case 4.2: θ n,k 1 i = θ n,k 2 i Since i ∈ NA n,k 1 ∪ NA n,k 2 , we have

θ n-1,k 1 i = θ n,k 1 i = θ n,k 2 i = θ n-1,k 2 i = σ i = ±1. From (2.37), i ∈ U n,k 1 α ∩ U n,k 2 α
and σ i τ n,k 1 i,α > τ n,k 2 i,α .

(2.38)

Sub-case 4.2.1:

i ∈ U n-1,k 1 α ∩ U n-1,k 2 α
In this case, we have, since the condition holds at step n -1

σ i τ n-1,k 1 i,α ≤ σ i τ n-1,k 2 i,α
.

Using the fact that i ∈ U n-1,k

1 α ∩ U n,k 1 α , i ∈ U n-1,k 2 α ∩ U n,k 2 α
joint to Proposition 4.1(iii), we get that τ n-1,k 1 i,α = τ n,k 1 i,α and τ n-1,k 2 i,α = τ n,k 2 i,α . This implies that

σ i τ n,k 1 i,α ≤ σ i τ n,k 2 i,α .
This contradicts (2.38).

Sub-case 4.2.2: i ∈ U n- where we have used the fact that i ∈ U n,k 2 α \ U n-1,k 2 α joint to Proposition 4.1(iv) for the first equality and the fact that i ∈ U n-1,k 1 α ∩ U n,k 1 α joint to Proposition 4.1(iii) for the last inequality. This contradicts (2.38).

Sub-case 4.2.2.2: σ i = -1 Let us treat the case α = +, the case α =being similar. We have that i ∈ U n,k 1 + (i + 1) ∩ U n,k 2 + (i + 1) This implies that f + (x i+1 ) > 0 and θ n,k 1 i+1 = θ n,k 2 i+1 = 1. Moreover, since θ n-1,k 1 i = -1 and i ∈ U n-1,k 1 + , we get that i ∈ U n-1,k 1 + (i + 1). This implies that θ n-1,k 1 i+1 = 1. Since i ∈ U n-1,k 2 + (i + 1), θ n-1,k 2 i = -1 and f + (x i+1 ) > 0, we deduce that θ n-1,k 2 = -1, we also deduce that U n-1,k 2 -(i + 1) = ∅. This implies that j ∈ NB n-1,k 2 . This contradicts the fact that θ n-1,k 2

Definition 2 . 1 .

 21 Let w be a real-valued function. the total variation of w is defined by: T V (w) := sup    j=1,...,k |w(y j+1 )w(y j )|; k ∈ N * , and (y j ) 1≤j≤k+1 non-decreasing sequence    .

1 .

 1 By Step 2, we deduce the result.

  Lemma 4.4. (τ i n a ,α is the time when the discontinuity a reached the point x i n a -1 2

Proposition 4. 6 .

 6 (No meeting for minima) Assume (H2) and ∆x ≤ ε L . Let (a i ) i=1,...,d , (b i ) i=1,...,d ⊂ (Z + 1 2

Figure 2 :

 2 Figure 2: Representation of the discontinuity before the meeting.

Lemma 4. 7 .

 7 (Numerical discontinuities meet before approached discontinuities) Let a, b ∈

Figure 4 :Figure 5 :

 45 Figure 4: at t = 12, with P = 50 points and N h = 50.

Figure 6 :Figure 7 :

 67 Figure 6: TEST2: UB scheme at t = 0, , 1.2, 1.6, with N x = 50 and p = 30, CFL=0.625

  and also by w the viscosity solution of (1.1) with initial data w

  ), where I = [X x,+ (-t), X x,-(-t)]. We can assume that w(1) 0 is non constant on I otherwise the result is obvious.

  ) is not minimal on I, i.e., w ) for some ā in I. If v 0 (ā) > p i=1 h i , then by definition of w 0,k we have w 0,k (ā) = 1 ∀k = 1, . . . , p, and thus w (1) 0 (ā) = p-1 k=1 h k . This contradicts the fact that w

  h k in all cases). Hence v 0 (ā) < p i=1 h i , and thus w(2) 0 (ā) = 0, from which we obtain that w(t, x) = w

  ). Since a ∈ I, this contradicts the fact that a is a minimum of w ) = inf y∈I w(1) 0 (y). Now we have also w

0

  is non constant on I (and in a we have v 0 (a) ≤ p k=1 h k ). Hence w (2) 0 (a) is minimal on I. (ii) In view of Proposition A.1, we have to estimate |v 0

i = 1 and θ n-1,k 2 i= - 1 . 1 +∩ 2 - 1 -∩ 2 + 1 +∩ N B n- 1 ,k 2 + 1 -∩ N B n- 1 ,k 2 - because θ n- 1 ,k 1 i= θ n- 1 ,k 2 i 2 - 2 i= - 1 2 i- 1 = 1 . 1 +

 21121211211211122212111 Since i ∈ NA n,k 1 ∩ NA n,k 2 , we have i ∈ NB n-1,k NB n-1,k or i ∈ NB n-1,k NB n-1,k (we cannot have i ∈ N B n-1,k or i ∈ N B n-1,kand the velocity is the same). Let us treat the first case, the other being similar. The fact that i ∈ NB n-1,k and θ n-1,k implies that θ n-1,k Similarly, i ∈ NB n-1,k and θ n-1,k 1 i = 1 implies that θ n-1,k 1 i-1

i+1= - 1 . 2 i+1= - 1 , 2 +

 1212 Since f + (x i+1 ) > 0 and θ n-1,k we deduce that U n-1,k (i + 1) = ∅. Since θ n-1,k 2 i = θ n-1,k 2

  i+1

j= - 1 = 1 α 2 α 1 i= t n = τ n,k 2 i 1 α\U n- 1 ,k 1 α, i ∈ U n,k 2 α\U n- 1 ,k 2 α

 11212111212 -θ n,k 2 j . Sub-case 4.2.3: i ∈ U n-1,k 2 α \U n-1,kThis case can be treated in the same way of Sub-case 4.2.2.Sub-case 4.2.4: i ∈ U n-1,k 1 α ∪ U n-1,kIn this case, we have τ n,kwhere we have used the fact that i ∈ U n,k joint to Proposition 4.1(iv) This contradicts (2.38).

  Let us consider steps (h k ) k=1,...,p such that h k > 0 and

			p
			k=1	h k > ||v 0 || ∞ , and let
		k
	hk :=	i=1	h i , for k ≥ 1, and h := sup 1≤k≤p	h k .	(2.3)

Table 1 :

 1 Advection with constant velocity

			Ulta-Bee
	25	3.74e-2	4.38e-2
	50	1.52e-2	1.67e-2
	100	7.41e-3	7.80e-3
	200	4.14e-3	4.24e-3
	Example 2:

Table 2 :

 2 Example 2, CFL= 0.5 Example 3:

	2. Discreti-

  1,k 1

α \U n-1,k 2 α Sub-case 4.2.2.1: σ i = 1

In this case, we have

τ n,k 2 i,α = t n ≥ τ n,k 1 i,α

,α are numerical fluxes that will be defined below, while Trunc denotes a truncation operator that will also be made precise below.

Classical Ultra-Bee scheme UB algorithm (combined to level-set decomposition)

Figure 9: TEST4


Laurent S

Anne R

Marie-Franc ¸oise M

Jean-Claude S

Mathieu G a , Sead T a , Didier B c André E

Edmond P C

A first generation bovine BAC-based physical map

A first generation clone-based physical map for the bovine genome was constructed combining, fluorescent double digestion fingerprinting and sequence tagged site (STS) marker screening. The BAC clones were selected from an Inra BAC library (105 984 clones) and a part of the CHORI-240 BAC library (26 500 clones). The contigs were anchored using the screening information for a total of 1303 markers (451 microsatellites, 471 genes, 127 EST, and 254 BAC ends). The final map, which consists of 6615 contigs assembled from 100 923 clones, will be a valuable tool for genomic research in ruminants, including targeted marker production, positional cloning or targeted sequencing of regions of specific interest.

INTRODUCTION

The detection of loci affecting economically important traits represents a major objective in livestock genomics. It should ultimately lead to more efficient breeding schemes (marker-assisted selection or MAS) and improve the accuracy and intensity of selection programs [START_REF] Georges | Livestock genomics comes of age[END_REF][START_REF] Haley | Livestock QTLs -bringing home the bacon?[END_REF]. In this perspective, genetic maps constructed in various livestock species [START_REF] Barendse | A medium-density genetic linkage map of the bovine genome[END_REF][START_REF] De Gortari | A second-generation linkage map of the sheep genome[END_REF][START_REF] Kappes | A second-generation linkage map of the bovine genome[END_REF][START_REF] Maddox | An enhanced linkage map of the sheep genome comprising more than 1000 loci[END_REF][START_REF] Vaiman | A genetic linkage map of the male goat genome[END_REF] are sufficient to detect regions containing genes and QTL. The identification and cloning of the corresponding genes may be achieved by standard positional cloning, taking advantage of the existence of large insert libraries and searching for transcribed sequences in these regions. Indeed, large insert libraries are available for livestock species, such as YAC libraries [START_REF] Alexander | Construction and characterization of a large insert porcine YAC library[END_REF][START_REF] Broom | Construction of a large-insert yeast artificial chromosome library from sheep DNA[END_REF][START_REF] Hills | A bovine YAC library containing four-to five-fold genome equivalents[END_REF][START_REF] Libert | Construction of a bovine genomic library of large yeast artificial chromosome clones[END_REF][START_REF] Rogel-Gaillard | Construction of a swine YAC library allowing an efficient recovery of unique and centromeric repeated sequences[END_REF] and BAC libraries [START_REF] Buitkamp | Construction and characterization of a gridded cattle BAC library[END_REF][START_REF] Cai | Construction and characterization of a bovine bacterial artificial chromosome library[END_REF][START_REF] Eggen | Construction and characterization of a bovine BAC library with four genome-equivalent coverage[END_REF][START_REF] Rogel-Gaillard | Construction of a swine BAC library: application to the characterization and mapping of porcine type C endoviral elements[END_REF][START_REF] Schibler | Construction and extensive characterization of a goat bacterial artificial chromosome library with threefold genome coverage[END_REF][START_REF] Vaiman | Construction and characterization of a sheep BAC library of three genome equivalents[END_REF][START_REF] Zhu | A 5× genome coverage bovine BAC library: production, characterization, and distribution[END_REF]. In addition, high-resolution comparative mapping projects have been developed in ruminants [START_REF] Band | An ordered comparative map of the cattle and human genomes[END_REF][START_REF] Schibler | Comparative gene mapping: A fine-scale survey of chromosome rearrangements between ruminants and humans[END_REF] and pigs [START_REF] Pinton | Localization of 113 anchor loci in pigs: improvement of the comparative map for humans, pigs, and goats[END_REF] in order to benefit from the spectacular progress of human and mouse genome projects, which will facilitate the comparative positional candidate gene approach. Likewise, systematic EST mapping projects have been initiated [START_REF] Smith | Sequence evaluation of four pooled-tissue normalized bovine cDNA libraries and construction of a gene index for cattle[END_REF] to efficiently address the crucial step of finding candidate genes for a trait of interest.

To date, genes responsible for economic trait loci have been identified in ruminants [START_REF] Grobet | A deletion in the bovine myostatin gene causes the double-muscled phenotype in cattle[END_REF][START_REF] Maddox | An enhanced linkage map of the sheep genome comprising more than 1000 loci[END_REF][START_REF] Mulsant | Mutation in bone morphogenetic protein receptor-IB is associated with increased ovulation rate in Booroola Merino ewes[END_REF]. In all cases, the authors used a contig assembly approach based on a walking strategy starting from only two entry points. Comparative mapping could be an efficient way to find new entry points and thus speed up the contig construction [START_REF] Schibler | Fine mapping suggests that the goat Polled Intersex Syndrome and the human Blepharophimosis Ptosis Epicanthus Syndrome map to a 100-kb homologous region[END_REF]. This sequential bi-directional walking strategy is however labor intensive and complicated by repetitive sequences (SINE or LINE) frequently encountered in either BAC end sequences or sequence tagged sites (STS). As in humans [START_REF] Chumakov | A YAC contig map of the human genome[END_REF][START_REF] Hudson | An STS-based map of the human genome [see comments[END_REF], the construction of genome wide physical maps in cattle will provide essential tools for efficient positional cloning, target DNA marker development or large-scale sequencing.

Several approaches have been developed to build physical maps, including hybridization-based methods [START_REF] Evans | Physical mapping of complex genomes by cosmid multiplex analysis[END_REF][START_REF] Mozo | Use of the IGF BAC library for physical mapping of the Arabidopsis thaliana genome[END_REF], restriction-based fingerprinting methods [START_REF] Coulson | Genome linking with yeast artificial chromosomes[END_REF][START_REF] Olson | Random-clone strategy for genomic restriction mapping in yeast[END_REF] or sequencing of a large number of BAC ends that serve as sequence-tagged connectors [START_REF] Venter | A new strategy for genome sequencing[END_REF]. Since the restriction fingerprinting approach is less hindered by the presence of repeated sequences than the hybridization methods, is faster and more automation friendly, it has been considered a reasonable compromise to detect and measure clone overlaps. Two different fingerprinting approaches were initially developed, using either single [START_REF] Olson | Random-clone strategy for genomic restriction mapping in yeast[END_REF] or double complete digestion [START_REF] Coulson | Toward a physical map of the genome of the nematode Caenorhabditis elegans[END_REF] methods. Fingerprinting techniques were later improved in order to increase data quality and to support the generation of fingerprints from a number of clones representative of large genomes [START_REF] Gregory | Genome mapping by fluorescent fingerprinting[END_REF][START_REF] Marra | High throughput fingerprint analysis of large-insert clones[END_REF]. Moreover, data analysis was greatly improved by the development of the FPC software [START_REF] Soderlund | Contigs built with fingerprints, markers, and FPC V4.7[END_REF], which combines automation and interactive graphics.

The availability of stable large-size clones and FPC software has renewed the fingerprinting method. Currently, BAC libraries have been extensively used to build numerous chromosome specific or whole genome sequence ready maps by fingerprinting [START_REF] Cao | A 12-Mb complete coverage BAC contig map in human chromosome 16 p13.1[END_REF][START_REF] Chen | An integrated physical and genetic map of the rice genome[END_REF][START_REF] Gregory | A physical map of the mouse genome[END_REF][START_REF] Hoskins | A BAC-based physical map of the major autosomes of Drosophila melanogaster[END_REF][START_REF] Mcpherson | A physical map of the human genome[END_REF][START_REF] Mozo | A complete BACbased physical map of the Arabidopsis thaliana genome[END_REF][START_REF] Tao | Bacterial artificial chromosome-based physical map of the rice genome constructed by restriction fingerprint analysis[END_REF]. Whole genome maps are being constructed for a number of organisms including the rat, cow, zebrafish, sorghum, maize and tomato (see www.genome.clemson.edu/fpc and www.bcgsc.edu for links to the corresponding web sites).

In this paper, we report the assembly of a first generation clone-based physical map, for the bovine genome, combining fluorescent double digestion fingerprinting and STS marker screening. This map represents the first report of genome wide physical mapping using the fluorescent double digestion fingerprinting method.

MATERIALS AND METHODS

BAC DNA preparation

BAC clones were grown in 2 mL 96-well blocks for 18 h 30 min at 37 • C with shaking, using 900 µL 2YT containing 12.5 µg•mL -1 chloramphenicol. Bacteria were pelleted by centrifugation at 3200 ×g for 8 min. The blocks were inverted to discard the supernatant and stored at -20 • C for one week until DNA preparation.

DNA preparation was performed using a modified alkaline lysis procedure. Cell pellets were resuspended by addition of 100 µL of TGE + 30 µg•mL -1 RNAse and vigorous vortexing. Lysis was achieved by a 3 min incubation step after addition of 100 µL of NaOH 0.2 M + SDS 1%. One hundred microliters of ice-cold potassium acetate (1.32 M, pH 4.8) were then added to each well, followed by 100 µL of lithium chloride (8 M). After 20 min at -20 • C, the lysis products were purified using Millipore filter plates (Multiscreen MANLY). The samples were then precipitated with 250 µL of isopropanol and centrifuged at 3200 ×g for 30 min with isopropanol. DNA pellets were washed with 250 µL of 70% ethanol per well and centrifuged for 10 min. The blocks were then inverted on paper towel to drain excess ethanol from the pellet and placed in a vacuum hybridization oven for 5 min to dry the DNA. Resuspension was achieved by addition of 10 µL of 10 mM Tris-HCl (pH 8.0) + 10 mg•mL -1 RNAse and vortexing after a 30 min incubation step at 37 • C.

The DNA concentration of 24 samples was estimated by fluorimetry for each plate. The plates were then calibrated at a mean concentration of 75 ng•µL -1 by addition of 10 mM Tris-HCl.

Fingerprinting

Restriction enzyme digestions and dye labeling were carried out simultaneously in a 10 µL reaction volume using 300-400 ng of BAC DNA, 2 units HindIII (Promega), 3 units HaeIII (Promega), 0.5 units Sequenase II (Amersham Pharmacia Biotech) and 3 pmoles ddA (R110, R6G or TAMRA from NEN) in restriction buffer C (Promega). The reactions were carried out for 90 min at 37 • C. In order to normalize signal intensities, R6G-labeled samples were diluted twice by adding 10 µL of sterile water.

Sample pooling and reaction cleanup

Multiscreen filter plates (Millipore, MAHVN45), filled with 360 µL of preswollen Sephadex G50 (G50 superfine 80 g•L -1 , Amersham Pharmacia Biotech) were used to achieve dye removal and sample desalting. Samples (10 µL each) from three plates labeled with different dyes were then pooled on the center of G50 mini-columns using a Hydra 96 dispenser. After a fiveminute incubation step, the samples were recovered by spinning down for 5 min at 910 ×g and stored at -20 • C until injection.

Sample injection and capillary electrophoresis

Seven microliters of the purified fingerprinting reaction were transferred to the injection plate and three microliters of a size-standard mix (2.8 µL formamide + 0.2 µL ET 900 ROX standard) were added. The samples were denatured in a heater block for 3 min at 90 • C prior to injection in a MegaBace 1000 automated 96 capillary DNA sequencer. Injection parameters and run conditions were respectively 3 kV for 90 s and 10 kV for 100 min using dye set 2 filters.

Data analysis and cleanup

The runs were analyzed with the Genetic Profiler software developed to perform the genotyping analysis on the MegaBace. After spectral matrix correction and peak identification, this software creates a "pks" file for each capillary that contains the scan value, the peak height and the peak width for all peaks detected in each channel. Since the Genetic Profiler software takes all peaks into account, it was necessary to eliminate background and artifactual peaks before exporting data to an FPC formatted file. Therefore, we developed a Visual Basic software (available upon request to the author) taking a "pks" file as input and applying three filters. The first one removes all peaks with a signal height lower than 200 units or not in the range of the analysis (55 to 750 bp); after an iterative estimation of the mean and standard deviations of the most likely distribution of peak height, a second filter eliminates external bands, the height of which is deviated by more than 2.5 standard deviations from the mean; the third filter removes peaks in the size range (+/-1 bp) of known artifacts (76 and 129 bp additional bands generated by FAM, 80 bp by TAMRA and 102 and 160 bp by R6G). This software rewrites "pks" files, so that modifications can be visualized by the Genetic Profiler. Good results were obtained, since about 70% of the sample were perfectly cleaned, except for data with low signals or when multiple clones (cross contaminations) were present in the same well. Therefore, to provide reliable data, all samples were manually checked and edited if necessary before exporting to FPC. As described in [START_REF] Ding | Contig assembly of bacterial artificial chromosome clones through multiplexed fluorescence-labeled fingerprinting[END_REF], we submited fragment sizes as ×10 in order to deal with some limitations of FPC: we could thus set the TOL value at 5 to take into account the fact that the standard deviation of fragment size measurements was less than one base.

BAC library screening

A PCR-based screening was performed on the Inra bovine BAC library as previously described [START_REF] Eggen | Construction and characterization of a bovine BAC library with four genome-equivalent coverage[END_REF]. Briefly, primer pairs for markers were obtained from the BOVMAP database (http://locus.jouy.inra.fr). PCR reactions were performed on PTC-100 thermocyclers (MJ Research) in a 15 µL reaction volume with 1× standard buffer supplemented with 125 µM dNTP, 1.5 mM MgCl 2 , 0.5 µM of each primer and 0.035 U•µL -1 Taq polymerase (Promega). The samples were preheated for 5 min at 94 • C, subjected to 35 cycles (94 • C for 20 s, 50-60 • C for 30 s and 72 • C for 30 s), and a final extension step of 5 min at 72 • C.

BAC end sequencing

The Nucleobond AX100 kit (Macherey-Nagel) was used to prepare a BAC DNA suitable for end sequencing according to the manufacturer's recommendations. Sequencing reactions were performed as previously described [START_REF] Schibler | Fine mapping suggests that the goat Polled Intersex Syndrome and the human Blepharophimosis Ptosis Epicanthus Syndrome map to a 100-kb homologous region[END_REF].

Contig assembly

The map was constructed with FPC starting from an initial stringent build and using an incremental process, which consisted in joining contigs together based on end-end comparisons. This time consuming step requires generally to manually fuse and reanalyze contigs at higher cutoff values. In order to automate this work and to trace weak joins introduced by manual editing, we developed a new strategy using virtual markers. In fact, FPC allows a less stringent cutoff for clones sharing one or more markers and the "Incremental Build Contig" (IBC) option takes into consideration new marker data to automatically merge and reanalyse contigs. We thus transcribed end-end comparisons results in term of marker data. For example, end-end comparisons showed that clone bI0053B12 matched clone bI0133E09 at a 2 × 10 -11 value. We thus defined a Virtual End (VE) Marker (called V I0053B12 2E-11) hitting these two clones. Marker data were incorporated in FPC through *.ace files using the replace marker command. The same strategy was used to incorporate singleton at higher cutoff.

In practice, the first map was obtained at a 10 -13 cutoff value using the pure Sulston method, with the Tolerance parameter set to 5 (0.5 base in fact). Contigs with more than 5 Q (questionable clones that FPC failed to place properly in the map) were automatically reanalyzed at a lower cutoff value using the DQER (from 10 -13 to 10 -25 ) from FPC software. Three large contigs with more than 10 Q remained unchanged after reanalysis. These contigs encompassed 1059 clones, which displayed more than 75% shared bands with an unusually high number of related clones. These clones were discarded as they may represent clones with large repeated or duplicated sequences. Moreover, all contigs encompassing more than 15 clones were manually edited and split when they appeared doubtful. Singletons were then incorporated up to a 10 -10 cutoff value using virtual markers (called S clone name CutOff in FPC).

The next building step was an end-end comparison at a 10 -12 cutoff value. We retained as valuable, only reciprocal and unique fusions. Moreover, contigs split by the last DQER step and manual editing, could not be merged. Then, we automatically fused contigs based on end-end comparisons using Virtual Ends (VE) Markers (called V Clone Name CutOff in FPC) and the IBC option (see above). Fusions, which produced contigs with more than 5 Q, were rejected by deleting the corresponding virtual marker and reanalyzing with IBC. Remaining VE markers are denoted validated VE markers in Table II. All steps required to create or delete VE markers, were done automatically using an Access Database, managing files containing the detailed end-end results (by printing in a file the FPC standard output), the summarized end-end results (by saving the FPC results window), the contig results (by saving the FPC By Contigs window) and the IBC results (by saving the FPC results window). Four additional steps were performed as described, making it possible to merge ends sequentially from 10 -11 to 10 -08 . In a last step, we merged contigs based on our framework markers and checked contigs with more than 50 clones to split manually those which seemed to be doubtful.

RESULTS

Fingerprinting of BAC clones

Two libraries were used: the Inra BAC library [START_REF] Eggen | Construction and characterization of a bovine BAC library with four genome-equivalent coverage[END_REF] with 105 984 clones and part of the CHORI-240 (www.chori.org/bacpac) BAC library (26 500 clones). About 9% of the clones were lost during the fingerprinting process, due to poor DNA preparation or low DNA yield (7%), injection failure or damaged capillaries (2%). Cross contamination was identified in about 2% of the wells.

Empty or small clones showing a low number of bands were rejected at the data cleanup step, in order to retain only clones with more than 7 bands for Inra clones and 11 bands for CHORI-240 clones. Moreover, since a clone with a high number of bands may correspond to two different clones in a same well, INRA clones with more than 63 bands and CHORI-240 clones with more than 75 bands were marked as cancelled in FPC. Finally, 102 725 clones were transferred to FPC, 40 of which were fingerprinted twice and 1802 were cancelled. Table I summarizes all these results and Figure 1 shows the distribution of band numbers after data cleanup.

Contig assembly

Preliminary maps built using different cutoff values show that cutoff values ranging from 10 -12 to 10 -08 produced contigs consistent with either the expected results based on statistical analysis [START_REF] Lander | Genomic mapping by fingerprinting random clones: a mathematical analysis[END_REF] or on our contigs built by chromosome walking (unpublished data). We used two criteria to evaluate the quality of the construct: the number of clones poorly incorporated (Q clones) and the number of contigs exhibiting at least one Q clone (Q contigs). As shown in Figure 2, higher cutoff values (10 -7 ) generate many false positive contigs (an increase in Q clones), with one contig containing up to 70 000 clones. Moreover, Figure 2 shows clearly that map buildings are more stable in terms of Q clones and Q contig number for cutoff values in the 10 -12 to 10 -09 range.

Moreover, we analyzed more precisely the differences between 10 -10 and 10 -12 cutoff value maps after multiple DQER steps. We observed that in most cases, contigs built at a cutoff value of 10 -10 represent end fusions of 2 or more contigs built at 10 -12 . Thus, contig assembly seemed to be reliable in the 10 -12 to 10 -09 cutoff range. However, to provide reliable data, we used an incremental process to build our final map, as described in the Materials and methods starting from a stringent build obtained at a 10 -13 cutoff. We stopped this incremental process at 10 -08 , since higher cutoff values generated too many non reciprocal and non unique end fusions (Fig. 3).

Table II shows the results obtained at each fusion step, the map encompassing 6890 contigs obtained using 4604 virtual markers and 4157 fusions.

Map validation and anchoring by PCR screening

A total of 1390 markers was screened on the Inra BAC library, 87 of them being absent from our library. Thus the screening information for a total of 1303 markers (451 microsatellites, 471 genes, 127 EST, and 254 BAC-ends) was used to anchor the contigs and validate our building strategy. These markers were derived from existing bovine genetic and radiation hybrid maps or according to their position on the human genome. In particular, we used 80 type I and II markers placed on a BTA26 comprehensive radiation hybrid map [START_REF] Gautier | A comprehensive radiation hybrid map of bovine chromosome 26 (BTA26): comparative chromosomal organization between HSA10q and BTA26 and BTA28[END_REF].

Fifty-five contigs contained at least one VE marker flanked by screened markers. In all cases, the locations of flanking markers were consistent, suggesting that VE markers introduced no error.

About thirty contigs were built in the course of this work by chromosome walking using framework preliminary maps with VE markers. BAC ends were sequenced and screened on the BAC library, making it possible to confirm 152 VE markers (data not shown). No erroneous fusion could be detected. These screening data made it possible to fuse some more contigs, resulting in the first generation bovine physical map.

The bovine physical map

The current release (June 2003) encompasses 6615 contigs, 747 of them being anchored. The contigs contain 15 clones on average (from two to 164 clones). Table III shows the number of contigs, the Q contigs and number of Q clones for 7 contig size classes.

Our 6615 contigs cover a total of 791 214 bands, each of them representing about 3.5 kb based on 82 000 Inra clones (∼120 kb) showing 35 bands and 19 000 CHORI clones (∼180 kb) showing 48 bands. Thus, our contigs cover about 400 kb on average and the map covers about 2769 Mb, i.e., more than 90% of the bovine genome. Singletons and cancelled clones containing large repetitive sequences represent about 117 000 additional bands (i.e., ∼14%). Therefore, the BAC map covers virtually all the bovine genome. However, this may represent a slight overestimation because of undetected overlaps.

In terms of anchored contigs, the map covers 171 500 bands representing about 600 Mb, i.e. 20% of the bovine genome.

The bovine physical map will be continuously updated, based on new marker screening data. It will be publicly available through the BovMap database and webFPC at http://locus.jouy.inra.fr/cgi-bin/lgbc/mapping/common/ intro2.pl?BASE=cattle.

Inra BAC clones can be obtained at the Biological Resources Center (http://www-crb.jouy.inra.fr).

A first draft of a BTA 26 physical map

In order to test the genome coverage of our contigs and the usefulness of our map, we screened a first set of 34 evenly spaced markers from BTA 26. It was thus possible to anchor 23 contigs and to detect two regions not represented in our library (Fig. 4). These 22 contigs covered about 30 Mb, i.e. 50% of the estimated 60 Mb total length of BTA26 [START_REF] Gautier | A comprehensive radiation hybrid map of bovine chromosome 26 (BTA26): comparative chromosomal organization between HSA10q and BTA26 and BTA28[END_REF]. Thus, screening of about 30 evenly spaced markers made it possible to increase from the chromosome coverage with anchored contigs 20% to 50%.

The remaining 44 markers from the BTA26 RH map [START_REF] Gautier | A comprehensive radiation hybrid map of bovine chromosome 26 (BTA26): comparative chromosomal organization between HSA10q and BTA26 and BTA28[END_REF] were then screened: five new contigs were identified and six previously identified contigs were enlarged or fused.

Finally this preliminary BTA26 physical map covers about 35 Mb, which correspond to about 60% of BTA26.

The results are concordant in terms of gene order with the radiation hybrid map data except for five contigs (Ctg520, Ctg146, Ctg525, Ctg530 and Ctg150). These discrepancies concern very closely related markers separated by a distance under the resolution limit of our panel [START_REF] Gautier | An extensive and comprehensive radiation hybrid map of bovine Chromosome 15: comparison with human Chromosome 11[END_REF].

In addition, these results are consistent in terms of contig length, except for Ctg146 and Ctg525 for which the physical sizes (1700 kb and 2050 kb) represent about 170 and 90 cR 3000 instead of 60 and 70 cR 3000 respectively, based on the converting ratio previously defined from human physical data [START_REF] Gautier | An extensive and comprehensive radiation hybrid map of bovine Chromosome 15: comparison with human Chromosome 11[END_REF].

DISCUSSION

We developed a first generation genome-wide BAC-based physical map of the bovine genome. This map consists of 6756 contigs assembled from 100 923 clones selected in two libraries. Restriction profiling was achieved using a fluorescent fingerprinting technique, based on capillary electrophoresis to analyze samples. This strategy has proven to be efficient in terms of speed and data quality. Four persons were sufficient to generate about 2000 fingerprint profiles per day on a MegaBace 1000, using little automation (one Hydra 96). Capillary electrophoresis reduces the variations of fragment sizes to less than 1 bp in a range of 50 to 750 bp, making it possible to more accurately dedect overlaps. However, efficient fragment analysis software is necessary to simplify and speed up the data cleanup process, which required two persons for six months.

Our iterative strategy makes it possible to not only automate the manual merging step but also to precisely monitor the contig assembly process and thus to stop it before data is corrupted. No error could be detected based on 50 contigs containing at least one VE marker flanked by screened markers. Since these 50 contigs represent 118 virtual fusions, we can expect, with a risk lower than 5%, that less than 2.5% of the contigs could be falsely branched by VE markers.

Moreover, the use of VE markers allows any user to trace all 'weak joints' introduced by manual editing. Indeed, the only way to detect a merge at a high cutoff value is to reanalyze this contig. VE markers clearly indicate which clones and which cutoff values were used to merge contigs. This could be of great interest for users accessing contigs through webFPC, for example.

The clones analyzed in this study were submitted to a double digestion, which generated about 40 bands on average. It should be possible to achieve a three to fourfold increase of band number since the MegaBace can accurately resolve fragments in the 50-750 bp range. A higher number of bands could be of great interest to detect smaller overlaps. For example, two clones with 30% overlap could share 11 bands out of 35 or 22 bands out of 70. Their coincidence scores would be 2E -08 and 9E -10 , respectively. Many restriction enzymes could be used simultaneously to generate more bands and thus detect 20% overlaps. This strategy could be an alternative to the strategy of Ding et al. [START_REF] Ding | Contig assembly of bacterial artificial chromosome clones through multiplexed fluorescence-labeled fingerprinting[END_REF]. However, these authors suggest that too many bands (∼100 bands) could be troublesome for FPC to handle. A three-digestion strategy generating about 70 bands and permitting 25-30% could thus be a good compromise.

Our BAC-based map will be a valuable tool for genomic research in ruminants, including targeted marker production, positional cloning or targeted sequencing of regions of specific interest. Even if our map provides only a four-genome equivalent coverage, it may not be worthwhile to spend more time, fingerprinting additional clones. Three elements support this statement:

-Firstly, STS content mapping may be more effective than fingerprinting for achieving gap closure and contig joining. About 3000 additional markers evenly distributed on the whole genome should be sufficient to achieve a 60% coverage with anchored contigs based on our results on BTA26. About 6000 markers would allow most gap closure. At present, comparative mapping data combined with the human sequence makes it possible to quickly identify these 6000 markers and to develop bovine specific primers from the numerous bovine EST available in databases. This is encouraging since the main applications of this physical map could be dedicated to large contig construction to assist positional cloning in the ETL regions.

-Secondly, our map provides a good framework to initiate a strategy similar to that of Gregory et al. [START_REF] Gregory | A physical map of the mouse genome[END_REF] and establish high-resolution syntenies among ruminant, human and mouse genomes. About 60 000 bovine BES from the CHORI-240 library have been submitted to GenBank, 10 000 of them corresponding to clones integrated in our contigs. End sequencing of singletons and clones at the end of our contigs will provide about 60 000 additional BES. BLASTN comparison with the human genome should thus provide at least 5000 significant "hits", making it possible to align most bovine contigs along the human genome. The deduced bovine contig juxtaposition could be helpful to identify potential overlapping contigs, which could be easily joined by PCR screening.

-Thirdly, an international physical map is under development (www.bcgsc.bc.ca/projects/bovine mapping, http://www.livestockgenomics. csiro.au/cattle.shtml) by analyzing single digest fingerprints obtained from 280 000 BAC. Amongst these, 18 913 clones from the CHORI-240 library were also incorporated in our map. Even if two different fingerprinting strategies were used, clones common to both maps could serve as anchors to identify news fusions between contigs from the two maps. Cross validation between these two independently constructed maps should provide a reliable framework to start whole genome sequencing projects.

Figure 1 .

 1 Figure 1. Band number distribution after data cleanup.

Figure 2 .

 2 Figure 2. CutOff effect on the reliability of the map.

Figure 3 .

 3 Figure 3. Possible and valuable fusions at each fusion step.

Figure 4 .

 4 Figure 4. First draft of the BTA26 physical map. Correspondence of the physical map of BTA26 (right) with the radiation hybrid map (left) are indicated by dotted lines. Markers and contigs from the first set are painted in black. The results from the second set are in gray. Contig numbers are indicated as Ctg. X correspond to clones absent in the library for a specific marker and arrows represent discrepancies between the radiation and the physical map.

Table I .

 I Library statistics.

	Library	Number Number Empty Small Number of Number of Mean nb Mean FPC Mean
		of clones of lost clones clones cancelled	clones of bands sizes (bp) size (bp)
			clones			clones	in FPC
	Inra	105 984	7%	5%	6%	1225	82 010	35	10 658 114 000
	CHORI-240 26 500	17%	6%	2%	577	18 913	48	15 451 180 000
	Total	132 484	9% 5.2% 5.2%	1802	100 923	37	11 609

Table II .

 II Results of the bovine map iterative automated assembly. clones poorly incorporated into contigs. Q contigs: contigs containing at least one Q clone. VE markers: virtual end marker obtained by ends-ends comparison in FPC. Validated VE markers: VE markers that do not increase the number of Q after merging the contigs and thus retained in subsequent analysis.

	Step	Contigs Q contigs Q Contigs	VE	Validated Q clones
				>5 Qs	markers VE Markers
	Build 10 -13 DQER 10 944	1741	7			2794
	Build VE 10 -12	9994	1785	12	1167	1064	3058
	Build VE 10 -11	9099	1835	12	1124	1051	3263
	Build VE 10 -10	8364	1855	12	949	820	3433
	Build VE 10 -09	7648	1894	12	877	811	3608
	Build VE 10 -08	6890	1942	12	1121	858	3848
	Q clones:

Table III .

 III The June 2003 bovine map.

		Number of clones in the contig >200 200-100 100-50 50-25 25-10 9-3	2	Total
	Contigs	0	10	196	874	2032 2815 688	6615
	Q contigs (>5 Q)	0	8	6	5	0	0	0	19
	Q clones		117	909	1752	1140	0	0	3918
	Anchored Contigs	0	10	103	223	264	132	15	747

ACKNOWLEDGEMENTS

We are grateful to Dr. H. Hayes for comments on this manuscript and thank colleagues of the Inra Animal Genetics Department for support during this work, especially the CTIG team and T. Coudert. This work was supported by grants from AGENAE.