Madalina Deaconu
email: madalina.deaconu@inria.fr

Samuel Herrmann
email: samuel.herrmann@u-bourgogne.fr

Initial-Boundary Value Problem for the heat equation -A stochastic algorithm

Keywords: Initial-Boundary Value Problem, heat equation, random walk, mean-value formula, heat balls, Riesz potential, submartingale, randomized algorithm. 2010 AMS subject classifications: primary 35K20, 65C05, 60G42; secondary 60J22

The Initial-Boundary Value Problem for the heat equation is solved by using a new algorithm based on a random walk on heat balls. Even if it represents a sophisticated generalization of the Walk on Spheres (WOS) algorithm introduced to solve the Dirichlet problem for Laplace's equation, its implementation is rather easy. The definition of the random walk is based on a new mean value formula for the heat equation. The convergence results and numerical examples permit to emphasize the efficiency and accuracy of the algorithm.

Introduction

In this paper, we study the Initial-Boundary Value Problem (IBVP) associated to the heat equation and develop a new method of simulation based on the Walk on Moving Sphere Algorithm (WOMS). The main objective is to construct an efficient approximation to the solution of the IBVP. The solution is a C 1,2 function u satisfying

   ∂ t u(t, x) = ∆ x u(t, x), ∀(t, x) ∈ R + × D, u(t, x) = f (t, x), ∀(t, x) ∈ R + × ∂D, u(0, x) = f 0 (x), ∀ x ∈ D, (1.1)
where f is a continuous function defined on R + × ∂D, f 0 is continuous on D and D denotes a bounded finitely connected domain in R d . For compatibility reasons we have also f (0, x) = f 0 (x).

The foundation stone of our work is the probabilistic representation for the solution of a partial differential equation. Suppose that we are looking for the solution u(t, x) of some PDE defined on the whole space R d . Under suitable hypothesis we can use the classical form u(t, x) = E[f (t, X t )] where (X t ) t∈R + is a stochastic process, satisfying a stochastic differential equation, and f a known function. In order to approximate u(t, x), the Strong Law of large Number allows us to construct Monte Carlo methods once we are able to propose an approximating procedure for the stochastic process (X t ) t∈R + .

The problem is more difficult when considering problems with boundary conditions. Nevertheless if some regularity is provided we can also find a probabilistic approach. A generic representation, for the solution of the Dirichlet problem in a domain D (the solution does not depend on time), is

u(x) = E f (X τ D ) exp - τ D 0 k(X s )ds - τ D 0 g(X t ) exp - t 0 k(X s )ds dt ,
where f, g, k are given functions, X 0 = x and τ D = inf{t ≥ 0; X t ∈ ∂D}. We refer to several classical books for more details [START_REF] Bass | Diffusions and elliptic operators. Probability and its Applications[END_REF][START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF][START_REF] Friedman | Stochastic differential equations and applications[END_REF][START_REF] Øksendal | Stochastic differential equations[END_REF]. The problem is hard to address as, in order to give an approximation, we need to approach the hitting time, the exit position and sometimes even the path of the process X t up to exit the domain D.

In particular situations we need to characterize either the hitting time τ D or the exit position X τ D , and these problems reveal quite difficult. The main goal of our work is to handle a more complex situation by unearthing numerical algorithms for the couple (τ D , X τ D ) itself.

To fix ideas and present a brief history, consider the simple Dirichlet problem for Laplace's equation in a smooth and bounded domain D ⊂ R d : ∆u(x) = 0, ∀x ∈ D u(x) = f (x), ∀x ∈ ∂D.

We recall the associated probabilistic representation: u(x) = E x [f (X τ D )] where (X t , t ≥ 0) here stands for the d-dimensional Brownian motion starting in x. The original idea in order to approximate u(x) by using the walk on spheres algorithm (WOS) goes back to Müller [START_REF] Muller | Some continuous Monte Carlo methods for the Dirichlet problem[END_REF]. The idea consists in constructing a step by step R d -valued Markov chain (x n , n ≥ 0) with initial point x 0 := x which converges towards a limit x ∞ , x ∞ and X τ D being identically distributed. Let us roughly describe (x n ): first, we choose S 0 the largest sphere centered in x 0 and included in D. The first exit point x 1 from the sphere S 0 for the Brownian motion starting from x 0 has an uniform distribution on ∂S 0 and is easy to sample.

The construction is pursued with the new starting point given by x 1 (see Figure 1). The algorithm goes on and stops while reaching the boundary ∂D. In order to avoid an infinite sequence of hitting times the stopping criteria of the algorithm includes a ε test: we stop the Markov chain as soon as δ(x n , ∂D) ≤ ε (δ represents here the Euclidean distance in R d ). Convergence results depending on ε and on the regularity of ∂D can be found in Müller [START_REF] Muller | Some continuous Monte Carlo methods for the Dirichlet problem[END_REF] and Mascagni and Hwuang [START_REF] Mascagni | ǫ-shell error analysis for "walk on spheres[END_REF]. Generalization of this result to a constant drift, by means of convergence theorems for discrete-time martingales, was proposed in the work of Villa-Moralès [START_REF] Villa-Morales | On the Dirichlet problem[END_REF], [START_REF] Villa-Morales | Solution of the Dirichlet problem for a linear second-order equation by the monte carlo method[END_REF]. Binder and Braverman [START_REF] Binder | The rate of convergence of the walk on spheres algorithm[END_REF] gave also the complete characterization of the rate of convergence for the WOS in terms of the local geometry of D. Other elliptic problems have been studied by Gurov, Whitlock and Dimov [START_REF] Gurov | A Grid Free Monte Carlo Algorithm for Solving Elliptic Boundary Value Problems[END_REF].

If needed, we can also approach the boundary hitting time by using the explicit form of its probability distribution function. However, a real difficult leap appears when we want to move from the simulation of X t to the simulation of (t, X t ). For example, if the domain is a sphere then X τ D can be simulated by the uniform random variable on the ∂D while τ D has an explicit pdf function which is not suited for numerical approaches as it depends on the Bessel function.

In previous works [START_REF] Deaconu | Hitting time for Bessel processes-walk on moving spheres algorithm (WoMS)[END_REF], [START_REF] Deaconu | Simulation of hitting times for Bessel processes with non integer dimension[END_REF], [START_REF] Deaconu | The walk on moving spheres: a new tool for simulating Brownian motion's exit time from a domain[END_REF] the authors discussed the connexion between the hitting times of the Bessel process and Brownian ones and introduced a new technique for approximating both the hitting time and the exit position. These previous studies on the hitting time form the foundation of our current work. We propose a new algorithm, involving a random walk on heat balls belonging to the domain [0, t] × D (see [START_REF] Evans | Partial differential equations[END_REF] p.53 for a definition

0 1 0 1 0 1 (T 0 , X 0 ) Time axis (T 1 , X 1 ) (T 2 , X 2 ) D Figure 2: Markov chain (T n , X n ) n≥0
of the heat ball) which approaches (τ D , X τ D ) in general domains. Thus we obtain a method for approximating the solution of the equation (1.1). Let us mention at this stage that Sabelfeld [START_REF] Karlovich | Monte Carlo methods in boundary value problems[END_REF] already described a random walk (τ n , Y n ) n≥0 for solving the Initial-Boundary Value Problem for the heat equation. His approach is essentially different: first of all his random walk is valued only on the boundary [0, t] × ∂D and secondly the main argument is based on solving an integral equation of the second kind rather than using Monte-Carlo techniques. A nicely written description of the method can be find in [START_REF] Sabelfeld | Random walks on boundary for solving PDEs[END_REF]. Let us just note that such algorithm which permits to evaluate u(t, x) is less accurate for large time t or non convex domains D.

Let us now introduce the main results concerning the algorithm Random Walk on Heat Balls which approximates (τ D , X τ D ), X being a d-dimensional Brownian motion. We just introduce first some preliminary notations: we recall that δ(x, ∂D) is the Euclidean distance between the point x and the boundary of the domain and introduce the function α(u, v) = min(u, e 2d δ 2 (v, ∂D)). In the following, (U n ) n≥1 stands for a sequence of independent uniformly distributed random vectors on [0, 1] ⌊d/2⌋+1 , Π U n denote the product of all its coordinates, (G n ) n≥1 is a sequence of independent standard Gaussian r.v. and (V n ) n≥1 a sequence of independent uniformly distributed random vectors on the unit sphere of dimension d, centered on the origin. We assume these three sequences to be independent. Let us define:

R n+1 := Π U n+1 2/d exp -(1 - 2 d ⌊ d 2 ⌋)G 2 n+1
and construct a sequence (T n , X n ) n≥0 by the following procedure (Figure 2).

ALGORITHM

Initialisation: Fix ε > 0. The initial value of the sequence (T n , X n ) is (T 0 , X 0 ) = (t, x).

Step n: The sequence is defined by recurrence as follows: for n ≥ 0,

T n+1 = T n -α(T n , X n )R n+1 , X n+1 = X n + 2 α(T n , X n )ψ d (R n+1 )V n+1 , where ψ d (t) = t log(t -d/2 ). Stop If α(T n , X n ) ≤ ε then N ε = n 1. If δ 2 (X Nε , ∂D) ≤ 2εd e then choose X ε ∈ ∂D such that δ(X Nε , X ε ) = δ(X Nε , ∂D) and define T ε := T Nε . 2. If δ 2 (X Nε , ∂D) > 2εd
e then set T ε = 0 and X ε := X Nε .

Algorithm outcomes: We get thus (T ε , X ε ) and N ε the number of steps.

We propose an approximation of the solution to (1.1) by using the definition:

u ε (t, x) = E (t,x) [f (T ε , X ε )1 {Xε∈∂D} ] + E (t,x) [f 0 (X ε )1 {Xε / ∈∂D} ], for (t, x) ∈ [0, T ] × D.
We will prove the convergence of this approximation in Proposition 4.1:

Convergence result. Let us assume that the Initial-boundary Value Problem (1.1) admits an unique C 1,2 ([0, T ] × D)-solution u, defined by (3.3). We introduce the approximation u ε given by (4.2). Then u ε converges towards u, as ε → 0, uniformly with respect to (t, x). Moreover there exist κ T,D (u) > 0 and ε > 0 such that

|u(t, x) -u ε (t, x)| ≤ κ T,D (u) √ ε, ∀ε ≤ ε 0 , (t, x) ∈ [0, T ] × D.
An important result, based on the construction of a submartingale related to the Riesz potential, completes the convergence of the algorithm: Efficiency result. Let D ⊂ B(0, 1) be a 0-thick domain. The number of steps N ε , of the approximation algorithm, is almost surely finite. Moreover there exist constants C > 0 and ε 0 > 0 both independent of (t, x) such that

E[N ε ] ≤ C| log ε|, for all ε ≤ ε 0 .
The material is organized as follows. In the second section we present mean value properties for the heat equation which plays a central role in the definition of the algorithm. The third section constructs the Random Walk on Heat Balls used to solve the Initial-Boundary Value Problem. In Section 4, we introduce the stopping procedure of the algorithm and prove the convergence result. The rate of the algorithm is also analyzed. We end up the paper with numerical results for two particular domains. These illustrations corroborate the accuracy of the algorithm.

A mean value property associated to the heat equation

In this section we will discuss the link between solutions of the heat equation and a particular version of the mean value property. This link is also an essential tool in the study of the classical Dirichlet problem.

Let us first note that due to the time reversion, the solution of the Initial-Boundary Value Problem for the heat equation is directly related to the solution of the Terminal-Boundary value problem for the backward heat equation (heat equation with negative diffusion). Due to this essential property, we are going to first present a mean value property for the backward heat equation and then deduce a similar property for the heat equation.

Let A be an open non empty set of The mean value formula (2.2) is quite different and more general than the classical formula associated to the heat equation (see, for instance, Theorem 3 on page 53 in [START_REF] Evans | Partial differential equations[END_REF]). Nevertheless, after some transformations on (2.2), it is possible to obtain the classical mean value property. These transformations consist in time reversion and integration with respect to a particular probability distribution function with compact support. These main ideas appear implicit in the proof of Proposition 2.4, the details being left to the reader.

R + × R d . Definition 2.1. A function h : A → R is said to be a reverse temperature in A if h is a C 1,2 -function satisfying ∂ t h(t, x) + ∆ x h(t, x) = 0, ∀(t, x) ∈ A. ( 2
Proof. Let a ∈ R + be defined by a = α d/2 Γ(d/2)2 d/2-1 and let us consider the associated function

ψ a,d (t) := 2t log a Γ(d/2)t d/2 2 d/2-1 .
We introduce (W t , t ≥ 0) a standard d-dimensional Brownian motion and define by τ a,d the following hitting time τ a,d = inf{t ≥ 0 :

W t = ψ a,d (t)}.
Let us just notice that this hitting time is bounded by α

= a Γ(d/2)2 d/2-1 2/d
and its distribution function is given by Proposition 5.1 in [START_REF] Deaconu | Simulation of hitting times for Bessel processes with non integer dimension[END_REF]

p a,d (t) = 1 2at ψ d a,d (t), 0 ≤ t ≤ α. (2.4)
Furthermore the exit location W τ a,d is uniformly distributed on the sphere of radius ψ a,d (τ a,d ).

Let us consider h a reverse temperature on A. By Itô's formula, we obtain

h(t + τ a,d , x + √ 2W τ a,d ) = h(t, x) + τ a,d 0 ∂ t h(t + s, x + √ 2W s ) ds + √ 2 τ a,d 0 ∂ x h(t + s, x + √ 2W s ) dW s + τ a,d 0 ∆ x h(t + s, x + √ 2W s ) ds.
If a is small enough, then (t + τ a,d , x + √ 2W τ a,d ) ∈ A a.s. Using the fact that h is a reverse temperature in A, in particular, the continuity of ∂ x h is known, we can prove that the stochastic integral introduced in the Itô formula is a martingale. Hence the stopping time theorem leads to

h(t, x) = E[h(t + τ a,d , x + √ 2W τ a,d )].
By (2.4), we get

h(t, x) = 1 σ(S d ) α 0 S d h(t + u, x + √ 2ψ a,d (u)y) 1 2au ψ d a,d (u) dσ(y) du.
We introduce the change of variable u = αs, such that s ∈]0,

H( √ α) = H(0) + H ′ (0) √ α + α 2 H ′′ (0) + o(α) (2.5)
where o(α) is uniform with respect to both s and y variables. The derivatives of H can be computed explicitly and we get:

H(0) = h(t, x), H ′ (0) = 2ψ d (s) d j=1 ∂ x j h(t, x)y j , H ′′ (0) = 2s∂ t h(t, x) + 4ψ 2 d (s) 1≤i,j≤d ∂ 2 x i x j h(t, x)y i y j .
Applying the mean value property to both sides of (2.5), we obtain

d j=1 ∂ x j h(t, x)A 0 j + ∂ t h(t, x)A 1 + 1≤i,j≤d ∂ x i x j h(t, x)A i,j = o(α), (2.6)
where

A 0 j = 2 √ α Γ(d/2) (s,y)∈[0,1]×S d 1 s ψ d+1 d (s)y j dσ(y) ds, A 1 = α Γ(d/2) (s,y)∈[0,1]×S d ψ d d (s) dσ(y) ds, A i,j = 2α Γ(d/2) (s,y)∈[0,1]×S d 1 s ψ d+2 d (s)y i y j dσ(y) ds.
By symmetry arguments, we have A 0 j = 0 and A i,j = 0 for i = j. Let X d be a random variable whose probability distribution function is

p d (t) = 1 Γ(d/2)t ψ d d (t)1 [0,1] (t).
Let us just notice that p d (t) = αp a,d (αt), p a,d being defined by (2.4). Then

X d = e -G
where G is a random variable which has the gamma distribution of parameters (d + 2)/2 and 2/d. In particular, X d has the same distribution as (U 1 . . . U (d+2)/2 ) 2/d if d is even (here (U i ) i∈N is a sequence of standard uniform independent random variables) and X d has the same distribution as (U

1 . . . U ⌊d+2⌋/2 ) 2/d e -G 2 /d if d is odd (here G is a standard Gaussian r.v. independent of the sequence (U i ) i ). Therefore if d is even, we deduce A 1 = αE[X d ] = αE[U 2/d 1 ]E[U 2/d 2 ] . . . E[U 2/d (d+2)/2 ] = α d d + 2 (d+2)/d .
For the odd case,

A 1 = αE[X d ] = αE[U 2/d 1 ]E[U 2/d 2 ] . . . E[U 2/d ⌊d+2⌋/2 ]E[e -G 2 /d ] = α d d + 2 ⌊d+2⌋/d E[e -G 2 /d ]. Let us now compute A i,i for 1 ≤ i ≤ n. First we observe that S d y 2 i dσ(y) = 1 d d j=1 S d y 2 j dσ(y) = 1 d .
So using a convenient change of variable, we get

A i,i = 2α dΓ(d/2) 1 0 1 s ψ d+2 d (s) ds = 2α(d + 2) d 2 Γ(d/2) Γ((d + 2)/2) 1 0 1 Γ((d + 2)/2) t (d+2)/d 1 t ψ d+2 d+2 (t) dt = α d + 2 d E[X (d+2)/d d+2 ] = α d d + 2 (d+2)/d if d is even, and A i,i = α d d+2 ⌊d+2⌋/d E[e -G 2 /d ] if d is odd. So we note that for any d ∈ N * , we proved that A 0 j = 0, A i,j = δ ij A 1
, where δ ij is the Kronecker's symbol. Equation (2.6) leads therefore to (2.1).

In order to prove the equivalence between the notion of reverse temperature and the mean value (MV) property defined in (2.2), we prove that the MV formula implies the regularity of the solution. This is the case when the dimension of the space is large enough. Intuitively the regularity increases as the dimension of the space increases.

(2.2). Then h is a C 1,2 (A, R)-function.
In the historical proof of the regularity associated to the Laplace operator (Proposition 2.5 in [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]), the key argument is to introduce a convolution with respect to a C ∞ -function with compact support. For the heat equation, one needs to handle quite differently and will not be able to prove regularity of infinite order in any case.

Proof. Let us assume that h : A → R is a bounded function satisfying the mean value property for α small enough, smaller than some α M > 0. We introduce a C

∞ (R + , R + )- probability distribution function function k whose support is included in [α m , α M ] with 0 < α m < α M . Rewriting the straightforward equality ∞ 0 k(α)h(t, x) dα = h(t, x)
leads to the following mean value property

h(t, x) = 1 2π d/2 (α,s,y)∈R + ×]0,1[×S d k(α)h(t + αs, x + 2 √ αψ d (s)y) ψ d d (s) s dσ(y) ds dα. (2.7)
In order to do calculation, it is more convenient in this situation to use spherical coordinates in S d ; for y ∈ S d we define (θ 1 , . . . , θ d-1 ) by

y 1 = cos(θ 1 ), y k = k-1 j=1 sin(θ j ) cos(θ k ), for 2 ≤ k ≤ d, where θ 1 , . . . , θ d-2 ∈ [0, π], θ d-1 ∈ [0, 2π] and θ d = 0. The change of measure is therefore given by dσ(y) ←→ sin d-2 θ 1 . . . sin 2 (θ d-3 ) sin(θ d-2 )dθ 1 . . . dθ d-1 .
Let us now consider another system of coordinates (u, z 1 , . . . , z d ) replacing (α, s, θ 1 , . . . , θ d-1 ) and given by:

u = αs, z 1 = 2 √ αψ d (s) cos(θ 1 ), z k = 2 √ αψ d (s) k-1 j=1 sin(θ j ) cos(θ k ), for 2 ≤ k ≤ d,
we recall that θ d = 0 by convention. Let us just observe that

z 2 = 4αψ 2 d (s) = 4αs log(s -d/2 ) = 4u log(s -d/2 ). Hence s = exp -z 2 2du
and α = u exp z 2 2du . The Jacobian determinant of the change of variable is equal to:

Jac := ∂(u, z 1 , z 2 , . . . , z d ) ∂(α, s, θ 1 , . . . , θ d-1 ) = (2s √ αψ ′ d (s) - √ αψ d (s))(2 √ αψ d (s)) d-1 sin d-2 θ 1 . . . sin 2 (θ d-3 ) sin(θ d-2 ) = -d 2 d-2 α d/2 sψ d (s) d-2 sin d-2 θ 1 . . . sin 2 (θ d-3 ) sin(θ d-2 ).
The mean value property therefore becomes:

h(t, x) = 1 2π d/2 (u,z)∈R + ×R d k(ue z 2 2du )h(t + u, x + z) z 2 e (2-d) z 2 4du d 2 d u d/2+1 dz du. (2.8)
Let us just note that, for the particular choice k(x) = d 2 r -d/2 x d/2-1 1 { x ≤r 2 /(4π)} with r > 0, we obtain the classical mean value formula presented in the statement of Theorem 3 p. 53 in [START_REF] Evans | Partial differential equations[END_REF]. Of course in this case the condition on the smoothness of k is not satisfied.

Since α > 0 if and only if u > 0 and since the support of k belongs to R * + , we can replace the integration domain R + × R d by R d+1 . The formula (2.8) can be written as the following convolution integral

h(t, x) = R d+1 h(u, z)w(u -t, z -x)dz du, (t, x) ∈ A, (2.9)
where

w(t, x) := k(te x 2 2dt ) 2π d/2 x 2 e (2-d) x 2 4dt d 2 d t d/2+1 , (t, x) ∈ R d+1 .
(2.10)

In fact, the support of the function w is compact due to k whose support belongs to [α m , α M ].

The aim is now to deduce the regularity property of h from that of w.

Step 1. Regularity with respect to the time variable.

For any (t 0 , x 0 ) ∈ A, we choose a small neighborhood of the form ]t 0η, t 0 + η[×B(x 0 , η) which is contained in A and we are going to prove that h is C 1 with respect to the time variable in this neighborhood.

Let us denote by

α := (u -t)e z-x 2 2d(u-t) then ∂α ∂t = -1 - z -x 2 2d(u -t) α (u -t) = -1 - z -x 2 2d(u -t) e z-x 2 2d(u-t) .
(2.11)

Introducing g : R → R + defined by g(x) = (d2 d+1 π d/2 ) -1 k(x)x 1-d/2
which is a C ∞ function with compact support due to the regularity of the function k and since the support of k does not contain a small neighborhood of the origin, we obtain w(u

-t, z -x) = g(α) z-x 2 (u-t) 2 . Hence ∂w ∂t (u -t, z -x) = z -x 2 g ′ (α) (u -t) 2 ∂α ∂t + 2g(α) (u -t) 3 = z -x 2 - g ′ (α) (u -t) 2 1 - z -x 2 2d(u -t) α (u -t) + 2g(α) (u -t) 3 .
(2.12)

Let us fix x ∈ R d . We can observe that both (t, u, z) → w(ut, zx) and (t, u, z) → ∂w ∂t (u -t, z -x) are continuous for t = u that is for (u, z) belonging to R ×R d \ ({t} ×R d ) (the complementary set is negligible for the Lebesgue measure). Moreover, due to the compact support of k, for any ε > 0 small enough, there exists a constant

κ ε > 0 such that z -x 2 ≤ 2d(u -t) log(α M /(u -t)) ≤ κ ε (u -t) 1-ε . (2.13)
Hence (2.12) implies the existence of C 0 , C 1 and C 2 independent of z, u and t such that

∂w ∂t (u -t, z -x) ≤ C 0 1 z -x 6/(1-ε)-2 + 1 z -x 8/(1-ε)-4 1 {0≤u≤C 1 } 1 { z-x ≤C 2 } (2.14)
for any (t, x) in ]t 0η, t 0 + η[×B(x 0 , η) ⊂ A. The right hand side of (2.14) is integrable as soon as the space dimension satisfies d > 4. Since h is a bounded function the Lebesgue theorem permits to apply results involving differentiation under the integral sign. The function h is therefore C 1 with respect to the time variable.

Step 2. Regularity with respect to the space variable.

The computation is quite similar as in the first step. We have:

∂w ∂x i (u -t, z -x) = - (z i -x i ) (u -t) 2 g ′ (α)α z -x 2 d(u -t) + 2g(α) , 1 ≤ i ≤ d, and
∂ 2 w ∂x i ∂x j (u -t, z -x) = α(z i -x i )(z j -x j ) d(u -t) 3 g ′′ (α)α + g ′ (α) d(u -t) z -x 2 + 4g ′ (α) + 1 (u -t) 2 g ′ (α)α d(u -t) z -x 2 + 2g(α) δ ij .
The second derivative is a continuous function on R d+1 \ {u = t}. Since g is C ∞ with a compact support, we know that α is bounded and therefore, using (2.13), we obtain the following bound:

∂ 2 w ∂x i ∂x j (u -t, z -x) ≤ C 0 z -x 4 (u -t) 4 + z -x 2 (u -t) 3 + 1 (u -t) 2 1 {0≤u≤C 1 } 1 { z-x ≤C 2 } ≤ C 0 1 z -x 8/(1-ε)-4 + 1 z -x 6/(1-ε)-2 + 1 z -x 4/(1-ε) 1 {0≤u≤C 1 } 1 { z-x ≤C 2 } .
(Here C 0 is just a generic constant: the values can change from one computation line to the following one). The conclusion is based on the same argument presented in Step 1: the boundedness of h permits to conclude that h is C 2 with respect to the space variable as soon as d > 4.

Let us note that the heat equation and the Laplace equation have some similar properties. In particular, solutions of these equations are automatically smooth. More precisely, if u ∈ C 1,2 (A, R) and is a reverse temperature then u ∈ C ∞ (A × R + , R) (see for instance Theorem 8 page 59 in [START_REF] Evans | Partial differential equations[END_REF]). In fact, as soon as the dimension d is large enough (d > 4), the mean value property implies the smoothness as an immediate consequence of Proposition 2.4.

All results presented so far in this section have an important advantage, they can be adapted to other situations for instance by looking backward in time, or equivalently time reverting. This observation permits to study properties of the heat equation.

Definition 2.5. A function

h : A → R is said to be a temperature in A if h is a C 1,2 -function satisfying the heat equation: ∂ t h(t, x) -∆ x h(t, x) = 0, ∀(t, x) ∈ A. ( 2

Solving the Initial-Boundary Value Problem

This section deals with existence and uniqueness for solutions of the Initial-Boundary Value Problem (1.1) in a bounded domain D. These results are deeply related to the existence of a particular time-discrete martingale: we define M n := (T n , X n ) a sequence of R + × D-valued random variables. In order to define this sequence we introduce δ(x, ∂D) the Euclidean distance between the point x and the boundary of the domain. We also introduce the function α given by: α(u, v) = min u, e 2d δ 2 (v, ∂D) .

(3.1)

Let us consider:

• (U n ) n≥1 a sequence of independent uniformly distributed random vectors on [0, 1] ⌊d/2⌋+1 . We denote by Π U n the product of all coordinates of U n .

• (G n ) n≥1 a sequence of independent standard Gaussian r.v.

• (V n ) n≥1 a sequence of independent uniformly distributed random vectors on the unit sphere of dimension d, centered on the origin.

Further, we assume that these three sequences are independent. We define by F n the natural filtration generated by the sequences (U n ), (G n ) and (V n ). Let F 0 note the trivial σ-algebra.

Let us introduce:

R n+1 := Π U n+1 2/d exp -1 - 2 d ⌊ d 2 ⌋ G 2 n+1 .
The initial value of the sequence (T n , X n ) is then (T 0 , X 0 ) = (t, x) and the sequence is defined by recurrence as follows: for n ≥ 0,

T n+1 = T n -α(T n , X n )R n+1 , X n+1 = X n + 2 α(T n , X n )ψ d (R n+1 )V n+1 . (3.2)
Let us first note that, due to the definition, the sequence (T n , X n ) belongs always to the closed set [0, t] × D: the sequence is therefore bounded. Moreover as soon as M n reaches the boundary of [0, t] × D its value is frozen.

Lemma 3.1. If h belongs to C 1,2 ([0, t] × D) and if it is a temperature in [0, t] × D, then M n := h(T n , X n ) is a bounded F -martingale.
Proof. Since h is a continuous function on a compact set, it is bounded. Therefore the stochastic process M n itself is bounded. We obtain

E[M n+1 |F n ] = E[h(T n+1 , X n+1 )|F n ] =: H(T n , X n ),
where

H(u, v) = E h u -α(u, v)R n+1 , v + 2 α(u, v)ψ d (R n+1 )V n+1 . Since the pdf of R n+1 is given by f R (s) = 1 Γ(d/2)
ψ d d (s)
s 1 [0,1] (s) and since V n+1 is uniformly distributed on the sphere, we obtain:

H(u, v) = 1 Γ(d/2)σ(S d ) (s,y)∈]0,1[×S d 1 s h(t -α(u, v)s, x + 2 α(u, v)ψ d (s)y) ψ d d (s) dσ(y) ds.
If h belongs to C 1,2 ([0, t] × D) and if it is a temperature in [0, t] × D, then Theorem 2.6 implies the mean value property. Hence H(u, v) = h(u, v). We deduce easily that

E[M n+1 |F n ] = h(T n , X n ) = M n a.s. Lemma 3.2. The process M n = (T n , X n ) converges almost surely as n → ∞ to a limit (T ∞ , X ∞ ) that belongs to the set {0} × D ∪ ]0, t[×∂D.
Proof. Let us consider the function h(t, x) = x i the i-th coordinate of x ∈ R d . We observe that h is a temperature and belongs to C 1,2 (R + , R d ). By Lemma 3.1, we deduce that

M n := h(T n , X n ) = X n (i), the i-th coordinate of X n
, is a bounded martingale therefore it converges a.s. towards X ∞ (i). Since all coordinates converge we deduce that X n → X ∞ a.s. Moreover since T n is a non-increasing sequence of non negative random times, it converges a.s. towards a r.v. T ∞ which belongs to [0, t]. The sequence (T n , X n ) belongs to the closed set [0, t] × D, consequently its limit belongs to the same set.

Let us change the starting point of the Markov chain by replacing (t, x) by (T ∞ , X ∞ ) we obtain a constant Markov chain M n = (T ∞ , X ∞ ) for all n ≥ 0. Let us assume that this limit does not belong to {0} × D ∪ ]0, t[×∂D then α(T ∞ , X ∞ ) > 0 a.s. We deduce that M 1 = M 0 a.s. for this new Markov chain since their first coordinates are different a.s. This fact obviously cannot be satisfied by a constant Markov chain, therefore (T ∞ , X ∞ ) ∈ {0} × D ∪ ]0, t[×∂D.

u(t, x) = E (t,x) [f (T ∞ , X ∞ )1 {X∞∈∂D} ] + E (t,x) [f 0 (X ∞ )1 {X∞ / ∈∂D} ], for (t, x) ∈ [0, T ] × D. (3.3)
Proof. By Lemma 3.1, the process M n is a bounded martingale. Moreover Lemma 3.2 implies that (T n , X n ) converges to (T ∞ , X ∞ ). Since u is a continuous function, we deduce that M n converges a.s. and in L 2 towards u(T ∞ , X ∞ ). In particular, the martingale property leads to

u(t, x) = E[u(T ∞ , X ∞ )].
In order to conclude it suffices to use the initial and boundary conditions. Indeed Lemma 3.2 ensures that (T ∞ , X ∞ ) belongs to the set {0} × D ∪ ]0, t[×∂D.

We refer to Friedman [START_REF] Friedman | Partial differential equations of parabolic type[END_REF] for the existence of a solution to the Initial-Boundary Value Problem (1.1). More precisely, if the following particular conditions are fulfilled:

• f and f 0 are continuous functions such that f (0, x) = f 0 (x),

• the domain has an outside strong sphere property, then there exists a smooth solution u to (1.1): u ∈ C ∞ (R + × D, R). This statement results from a combination of Theorem 9 page 69 and Corollary 2 page 74 in [START_REF] Friedman | Partial differential equations of parabolic type[END_REF].

Approximation of the solution for an Initial-Boundary Value Problem

The aim of this section is to contruct an algorithm which approximates u(t, x), the solution of an Initial-Boundary Value Problem when (t, x) is given. For the Dirichlet problem such an algorithm was introduced by Müller [START_REF] Muller | Some continuous Monte Carlo methods for the Dirichlet problem[END_REF] and is called the Random Walk on spheres. We are concerned with the heat equation instead of the Laplace equation and therefore propose an adaptation of this algorithm in order to consider also the time variable. The algorithm is based on the sequence of random variables M n = (T n , X n ) defined by (3.2). We introduce a stopping rule: let ε > 0 be a small parameter, we define N ε the stopping time:

N ε := inf{n ≥ 0 : α(T n , X n ) ≤ ε}, (4.1)
where α is given by (3.1).

1. If δ 2 (X Nε , ∂D) ≤ 2εd e then we choose X ε ∈ ∂D such that δ(X Nε , X ε ) = δ(X Nε , ∂D)

and we denote by T ε := T Nε .

2. If δ 2 (X Nε , ∂D) > 2εd e then we set T ε = 0 and X ε := X Nε . We are now able to give an approximation of the solution to (1.1) by using the definition: 3). We introduce the approximation u ε given by (4.2). Then u ε converges towards u, as ε → 0, uniformly with respect to (t, x). Moreover there exist κ T,D (u) > 0 and ε > 0 such that

u ε (t, x) = E (t,x) [f (T ε , X ε )1 {Xε∈∂D} ] + E (t,x) [f 0 (X ε )1 {Xε / ∈∂D} ], for (t, x) ∈ [0, T ] × D.
|u(t, x) -u ε (t, x)| ≤ κ T,D (u) √ ε, ∀ε ≤ ε 0 , (t, x) ∈ [0, T ] × D.
Proof. Using the definition of u (resp. u ε ) in (3.3) (resp. (4.2)), we obtain:

|u(t, x) -u ε (t, x)| = |E[u(T ∞ , X ∞ )] -E[u(T ε , X ε )]|.
Since n → u(T n , X n ) is a bounded martingale and since N ε is a finite stopping time, we can apply the optimal stopping theorem leading to

|u(t, x) -u ε (t, x)| = |E[u(T Nε , X Nε )] -E[u(T ε , X ε )]| ≤ κT,D (u)E max(|T Nε -T ε |, |X Nε -X ε |)
where κT,D (u) := sup

(t,x)∈[0,T ]×D max ∂u ∂t (t, x) , ∂u ∂x (t, x) .
Taking into account the two different situations δ 2 (X Nε , ∂D) > 2εd e or δ 2 (X Nε , ∂D) ≤ 2εd e , we deduce that

max(|T Nε -T ε |, |X Nε -X ε |) ≤ max ε, 2εd e .
The statement follows with the particular choice κ T,D (u) = κT,D (u) 2d e .

Let us now focus our attention on the number of steps needed by the algorithm (3.2) before stopping. In order to present the main result, we need some particular properties on the domain D.

In the sequel, we shall assume that D is a 0-thick domain, that is: there exists a constant C > 0 (so-called the thickness of the domain) such that

H d (B(x, r) \ D) ≥ Cr d , ∀r < 1, ∀x ∈ ∂D. (4.3)
Here H d (S) denotes the d-dimensional Hausdorff content of the set S. This property is namely satisfied by

• convex domains;

• domains satisfying a cone condition;

• bounded domains with a smooth boundary ∂D.

We observe that the assumption is quite weak. For such domains, we can prove the following rate of convergence. Theorem 4.2. Let D ⊂ B(0, 1) be a 0-thick domain. The number of steps N ε , of the approximation algorithm, is almost surely finite. Moreover there exist constants C > 0 and

ε 0 > 0 such that E[N ε ] ≤ C| log ε|, for all ε ≤ ε 0 . (4.4)
The proof of this result is an adaptation of the classical random walk on spheres [START_REF] Binder | The rate of convergence of the walk on spheres algorithm[END_REF]. Nevertheless the dynamics of both coordinates of the random walk on spheres (T n , X n ) n≥0 being definitively different, this adaptation requires a quite tedious effort. In particular, we need to introduce a particular submartingale, based on the random walk, whose properties permit to prove the rate of convergence.

Submartingale related to the Riesz potential

We consider in this section the 0-thick domain D which is included in the unit ball of R d (assumption of Theorem 4.2). We introduce the set M of all Borel measures µ supported inside B(0, 2) and outside of D, satisfying the following condition:

µ(B(x, r)) ≤ r d , ∀x ∈ R d , ∀r > 0. (4.5)
Let us define the so-called energy function U:

U(x) = sup µ∈M U µ (x),
where U µ stands for the Riesz potential of the measure µ, that is,

U µ (x) = ∞ 0 µ(B(x, r)) r d+1 dr. (4.6)
The definition of U obviously implies that U(x) ≥ 0 for any x ∈ R d .

Remark 4.3. Binder and Bravermann [START_REF] Binder | The rate of convergence of the walk on spheres algorithm[END_REF] gave several properties of this energy function.

We just recall some of them:

1. Since the set of measures M is weakly * -compact, there exists a family of measures µ x , belonging to M, such that U(x) = U µx (x). This property will play a crucial role in the proof of Proposition B.2.

2. The energy function U is subharmonic in D. Consequently, due to the construction of the random walk (T n , X n ) which is based on uniform random variables on moving spheres, the process (U(X n )) n≥0 is a submartingale with respect to the filtration generated by (T n , X n ) n≥0 :

F n := σ{T 1 , . . . , T n , X 1 . . . , X n }. Hence E[U(X n+1 )|F n ] ≥ U(X n ) a.s. (4.7)
3. Easy computations on the Riesz potential permits to prove that

U(x) ≤ log 3 δ(x, ∂D) + 1 d , ∀x ∈ D. (4.8)
In particular, we obtain an important property of the submartingale (U(X n )) n≥0 :

If U(X n ) ≥ log 3 ε + 1 d than δ(X n , D) ≤ ε. (4.9)
These properties, of the energy function U, permit to sketch the proof of the convergence in the classical random walk on spheres case. Indeed we know that U(X n ) is a submartingale and the algorithm stops before U(X n ) becomes too large. So, it suffices to focus the attention on the time needed by the submartingale to exceed some given large threshold.

In the algorithm described by (3.2), a large value of U(X n ) is not sufficient to ensure that the stopping rule has been reached. Indeed the stopping procedure depends on both space and time variables, through the condition: α(T n , X n ) ≤ ε. That's why we need to adapt the classical study by considering a martingale based on the Riesz potential but taking also into account the decreasing time sequence (T n ) n≥0 .

Let us define the modified energy function

on R + × R d by U(t, x) = max 1 2 log 3 t , U(x) . (4.10)
This function will play a similar role as the energy function (in the classical case). In particular, if we apply U to the sequence (T n , X n ), we obtain a submartingale with nice properties.

Lemma 4.4. We define

U n := U(T n , X n ). Then the process (U n ) n≥0 is a F -submartingale.
Proof. First we can notice that

T n+k ≤ T n ⇒ log 3 T n+k ≥ log 3 T n .
By Jensen's inequality (see Lemma A.1) and using the submartingale property (4.7) of U(X n ), we obtain

E U n+1 σ(T n , X n ) = E max 1 2 log 3 T n+1 , U(X n+1 ) σ(T n , X n ) ≥ E max 1 2 log 3 T n , U(X n+1 ) σ(T n , X n ) ≥ max 1 2 log 3 T n , E[U(X n+1 )|σ(T n , X n )] ≥ max 1 2 log 3 T n , U(X n ) = U n .
We deduce that (U n ) is a submartingale as announced.

In order to describe an upper-bound for the sequence E[U 2 n ], we first point out an inequality relating U(t, x) to the function α(t, x), which plays an essential role in the algorithm (3.2). Lemma 4.5. There exists a constant κ > 0 (depending only on the space dimension d) such that

U(t, x) ≤ κ - 1 2 log(α(t, x)), ∀(t, x) ∈ R + × D. (4.11)
Proof. On one hand, the definition of α(t, x) in (3.1) implies that e 2d δ 2 (x, ∂D) ≥ α(t, x) and consequently

log 3 δ(x, ∂D) ≤ 1 2 log 9e 2d - 1 2 log α(t, x).
Using the property (4.8), we obtain

U(x) ≤ - 1 2 log α(t, x) + 1 2 log 9e 2d + 1 d .
On the other hand, the definition of α(t, x) also implies

1 2 log 3 t ≤ 1 2 log 3 α(t, x) = log 3 2 - 1 2 log α(t, x).
Combining both inequalities, we deduce that U(t, x) = max{ 1 2 log(3/t), U(x)} satisfies (4.11) with κ := max{ 1 2 log 3,

1 2 log 9e 2d + 1 d }. An immediate consequence of Lemma 4.5 is an L 2 -bound of U n , n fixed.
Proposition 4.6. Let (T 0 , X 0 ) = (t, x). There exist two constants C 1 and C 2 such that,

E[U 2 n ] ≤ (C 1 + C 2 n) 2 , for n ≥ κ - 1 2 log(α(t, x)).
Here κ stands for the constant defined in Lemma 4.5.

Proof. Let us first recall that U n := U(T n , X n ). Due to the definition of the function U, we observe that U(x) ≥ 0 for any x ∈ R d and consequently U(t, x) ≥ 0 and U n ≥ 0. Due to Lemma 4.5, we shall focus our attention on log α(T n , X n ). First we notice that (3.2) leads to

T n = T n-1 -α(T n-1 , X n-1 )R n ≥ T n-1 (1 -R n ). Hence -log(T n ) ≤ -log(T n-1 ) -log(1 -R n ) ≤ -log(α(T n-1 , X n-1 )) -log(1 -R n ). (4.12)
Moreover by (3.2),

δ(X n , ∂D) ≥ δ(X n-1 , ∂D) -2 α(T n-1 , X n-1 )ψ(R n ).
By its definition, α(t, x) ≤ e 2d δ 2 (x, ∂D), and we obtain

δ(X n , ∂D) ≥ δ(X n-1 , ∂D) 1 - 2e d ψ(R n ) ,
and therefore

-log e 2d δ 2 (X n , ∂D) ≤ -log(α(T n-1 , X n-1 )) -2 log 1 - 2e d ψ(R n ) . (4.13)
Let us define

W n := -2 log 1 -2e d ψ(R n ) -log(1 -R n ).
Combining (4.12) and (4.13) we finally obtain

-log(α(T n , X n )) ≤ -log(α(T n-1 , X n-1 )) + W n ≤ -log(α(t, x)) + n k=1 W k . (4.14)
Let us just note that (W n ) is a family of independent and identically distributed random variables and (t, x) is the starting position of the algorithm. Let us recall that U n ≥ 0. We obtain

E[U 2 n ] ≤ E κ - 1 2 log(α(T n , X n )) 2 ≤ E κ - 1 2 log(α(T 0 , X 0 )) + n k=1 W k 2 ≤ 2 κ - 1 2 log(α(t, x)) 2 + 2E n k=1 W k 2 ≤ 2 κ - 1 2 log(α(t, x)) 2 + 2Var n k=1 W k + 2E n k=1 W k 2 ≤ 2 κ - 1 2 log(α(t, x)) 2 + 2nVar(W 1 ) + 2n 2 (E[W 1 ]) 2 ≤ 2nVar(W 1 ) + 2n 2 (E[W 1 ] 2 + 1),
due to the hypothesis n ≥ κ -1 2 log(α(t, x)). So Lemma A.2 implies the statement of the Proposition 4.6: the upper-bound is quadratic with respect to n.

Let us now point out a lower-bound for the expected value of the submartingale: (E[U n ]) n≥0 . Proposition 4.7. There exist two constants C 3 ∈ R and C 4 > 0, such that

E[U n ] ≥ C 3 + C 4 n, n ≥ 1. (4.15)
Proof. Since (U n ) n≥0 is a submartingale we know that (E[U n ]) n≥0 is a non-decreasing sequence, but we need even more. In fact, due to the following lower-bound:

U n ≥ 1 4 log 3 T n + 1 2 U(X n ) =: V n ,
it suffices to point out the existence of a constant L 0 > 0 such that

E[V n+1 ] -E[V n ] ≥ L 0 , ∀n ≥ 0. (4.16)
In order to compute such a lower-bound, we consider two cases: either α

(T n , X n ) = T n (event denoted by T n ) or α(T n , X n ) = T n (event denoted by T n ).
Step 1. First case: α(T n , X n ) = T n . Then the definition of the random walk (3.2) implies that

T n+1 = T n -α(T n , X n )R n+1 = T n (1 -R n+1 ) on T n . Hence log 3 T n+1 = log 3 T n -log(1 -R n+1 ) on T n .
Let us denote by

L 1 = -1 4 E[log(1 -R 1 )] > 0. Since U(X n ) is a submartingale (U being subharmonic in D), we get E V n+1 1 Tn = E 1 4 log 3 T n+1 1 Tn + 1 2 E E[U(X n+1 )|F n ]1 Tn ≥ E 1 4 log 3 T n 1 Tn + L 1 P(T n ) + 1 2 E U(X n )1 Tn ≥ E V n 1 Tn + L 1 P(T n ). (4.17)
Step 2. Second case: α(T n , X n ) = T n . Let us recall that the random walk satisfies:

X n+1 = X n + 2 α(T n , X n )ψ d (R n+1 )V n+1
where R n+1 is a continuous random variable whose support is the whole interval [0, 1] and whose distribution does not depend on n. Observe also that

ρ(R n+1 ) := 2e d ψ d (R n+1 )
is also a continuous random variable with support [0, 1]. In other words, on the event α(T n , X n ) = T n and given R n+1 = r, the (n + 1)-th step of the random walk is exactly the same as the (n + 1)-th step of the classical random walk on spheres (see the Appendix B)

with radius β = ρ(r), for which we can obtain some lower-bound. So using Proposition B.2, we obtain

E V n+1 1 T n = E 1 4 log 3 T n+1 1 T n + 1 2 E[U(X n+1 )1 T n ] ≥ E 1 4 log 3 T n 1 T n + 1 2 E E[U(X n+1 )|σ(R n+1 , T n , X n )]1 T n ≥ E 1 4 log 3 T n 1 T n + 1 2 E[(U(X n ) + L1 {δ/4<1-ρ(R n+1 )<δ/2} )1 T n ] ≥ E V n 1 T n + L 2 P({δ/4 < 1 -ρ(R n+1 ) < δ/2} ∩ T n ) (4.18) = E V n 1 T n + L 2 P(T n ), (4.19)
where

L 2 = L 2 P(δ/4 < 1 -ρ(R 1 ) < δ/2
), R n+1 and (T n , X n ) being independent. Finally taking the sum of (4.17) and (4.18), we obtain (4.16) with L 0 = min(L 1 , L 2 ) > 0.

We end here the preliminary results concerning the submartingale (U n ) n≥0 . We are now ready to deal with the rate of convergence of the random walk on moving spheres.

Rate of convergence of the algorithm

Let us consider the algorithm (T n , X n ) n≥0 given by (3.2) and stopped as soon as α(T n , X n ) ≤ ε. We assume that the starting position satisfies (T 0 , X 0 ) = (t, x) ∈ R + × D. Then the mean number of steps is bounded and the bound depends on | log ε|.

Proof of Theorem 4.2. If the starting position (t, x) satisfies α(t, x) ≤ ε then the algorithm stops immediately (N ε = 0 a.s.) and the statement is satisfied. From now on, we assume that α(t, x) > ε.

Step 1. A remark on the stopping rule. The statement of Theorem 4.2 concerns N ε , see (4.1), the first time the random walk (T n , X n ) n≥0 hits a ε-neighborhood of the boundary. Let us introduce another stopping rule concerning U n := U(T n , X n ), U being defined by (4.10):

N ′ ε := inf n ≥ 0 : U n ≥ log 3 ε + 1 d .
Let us now point out that N ε ≤ N ′ ε a.s. for ε small enough (more precisely, we need ε ≤ 2d e ). Indeed, let us consider the first case: U(X n ) ≥ log 3 ε + 1 d , then (4.9) implies that δ(X n , ∂D) ≤ ε. Moreover, due to the condition ε ≤ 2d e , we get e 2d δ 2 (X n , ∂D) ≤ ε and therefore α(T n , X n ) ≤ ε. On the other side, if 1 2 log 3

Tn ≥ log 3 ε + 1 d ≥ 1 2 log 3 ε then T n ≤ ε and finally α(T n , X n ) ≤ ε. So we deduce that U n ≥ log 3 ε + 1 d implies that α(T n , X n ) ≤ ε.
In the sequel, we will find an upper-bound for the mean value of N ′ ε .

Step 2. The aim of the second step is to prove the existence of an integer η ∈ N and a constant p < 1 both independent with respect to the starting position of the random walk (T 0 , X 0 ) = (t, x) and independent of the parameter ε such that

P(N ′ ε > η⌊-log ε⌋) ≤ p, (4.20)
for ε small enough. Let us note that the complementary event satisfies, by definition,

P(N ′ ε ≤ η⌊-log ε⌋) ≥ P(U k ≥ β k ),
where k = η⌊-log ε⌋ and β k = log 3 + 1 + 1 d + k/η. We deduce that there exists a particular choice of the integer η such that, for ε small enough, β k < α k := (C 3 + C 4 k)/2 where C 3 and C 4 are defined in Proposition 4.7. So it is sufficient to find a lower-bound of P(U k > α k ) which should be positive when k is large. By Proposition 4.6, there exist two constants C 1 and C 2 such that

E[U 2 n ] ≤ (C 1 + C 2 n) 2 , for any n ≥ κ - 1 2 log(α(t, x)).
Due to the condition on the initial position α(t, x) > ε, the previous inequality is satisfied for n ≥ ⌊log ε⌋ when ε is small enough. In particular, it is satisfied for n = k = η⌊-log ε⌋.

We obtain

E[U k ] = E U k 1 {U k ≤α k } + E U k 1 {U k >α k } > C 3 + C 4 k. (4.21)
Then by an application of (4.15) and the Cauchy-Schwarz inequality, we get

α k + E[U 2 k ] P(U k > α k ) > C 3 + C 4 k.
Therefore, due to the upper-bound of the second moment,

α k + (C 1 + C 2 k) P(U k > α k ) > C 3 + C 4 k.
We deduce

P(U k > α k ) ≥ 1 4 C 3 + C 4 k C 1 + C 2 k 2 > 1 5 C 4 C 2 2 ,
for k large enough that is ε small enough. This implies the existence of the constant p > 0 in (4.20).

Step 3. Upper-bound of E[N ′ ε ]. Due to the first step it is sufficient to obtain an upper-bound of E[N ′ ε ] in order to prove the statement of the theorem. Such a result is essentially based on the Markov property of the sequence (T n , X n ) n≥0 : the second step implies in particular that P(N ′ ε > kη⌊-log ε⌋) ≤ p k , ∀k ≥ 1.

Hence

E[N ′ ε ] ≤ k≥1 kη⌊-log ε⌋P N ′ ε ≤ kη⌊-log ε⌋ N ′ ε > (k -1)η⌊-log ε⌋ ×P(N ′ ε > (k -1)η⌊-log ε⌋) ≤ η⌊-log ε⌋ k≥1 kp k-1 = η| log ε| (1 -p) 2 .

Examples and numerics

The aim of this section is to illustrate the random walk on spheres algorithm introduced in Section 4. Let us focus our attention on the numerical approximation of the solution to the value problem:

   ∂ t u(t, x) -∆ x u(t, x) = 0, ∀(t, x) ∈ R + × D, u(t, x) = f (t, x), ∀(t, x) ∈ R + × ∂D, u(0, x) = f 0 (x), ∀ x ∈ D, (5.1)
for particular domains D. First we shall present results obtained for the hypercube D =]0, 1[ d and secondly the half of a sphere

D = {x ∈ R d : x < 1, x 1 > 0}.
Of course these toy examples are not directly related to concrete situations in physics but they permit to emphasize the efficiency of the algorithm. Their advantage relies in the easy computation of the distance to the boundary. For more general situations, only this part of the procedure has to be modified and can sometimes become quite painful.

Hypercube

Let us first introduce the functions which take part to the boundary conditions. We choose a function with the following simple expression

f (t, x) = e t d i=1 x i (1 -x i )1 ]0,1[ (x i ), x ∈ D, ∀t ≥ 0. ( 5.2)
Setting f 0 (x) = f (0, x), we observe that both the compatibility and the continuity conditions are obviously satisfied. In this particular case, we have already pointed out, in the previous sections, that there exists a unique (smooth) solution to the Initial-Boundary Value Problem which can be approximated using the algorithm of moving spheres.

The solution can be approximated by u ǫ defined by (4.2), the error being directly related to the parameter ǫ. Since u ǫ (t, x) is the expectation of a random variable, we shall use a Monte-Carlo method in order to obtain an estimated value. Hence

u ε N (t, x) = 1 N N k=1 f (T ε,k , X ε,k )1 {X ε,k ∈∂D} + 1 N N k=1 f 0 (X ε,k )1 {X ε,k / ∈∂D} , (5.3)
where (T ε,k , X ε,k ) k≥0 is a sequence of independent and identically distributed couples of random variables, the distribution being defined at the begining of Section 4. The difference between u(t, x) and u ε N (t, x) actually relies on both the error described in Proposition 4.1 of order √ ε on one hand and the classical error of Monte Carlo methods of order N -1/2 on the other hand (the confidence interval depends as usual on the standard deviation of the underlying random variable). First let us present u ε N (t, x) for a particular point: the center of the hypercube (x = (0.5, . . . , 0.5) is the default setting in all this subsection) letting the time cross the whole interval [0, 2]. We present at the same time the associated Monte-Carlo 95%-confidence interval (Figure 3). Let us just notice that the choice N = 1 000 is not motivated by some computational facilities but rather to produce a clear picture, the confidence interval becoming very small for larger values of N. Of course the numerical method permits to observe directly the distribution of the random variable

Z ε = f (T ε , X ε )1 {Xε∈∂D} + f 0 (X ε )1 {X ε,k / ∈∂D} ,
which drastically changes as time elapses (Figure 4 and 5).

In our example, small values of Z ε are more frequently observed for small time values than for large ones. Such behaviour of the random variable is not linked to the particular boundary conditions we introduced, but relies on the following general argument. The random variable Z ε is obtained due to a stopping procedure on M n = (T n , X n ) defined by (3.2). The sequence is stopped as soon as either X n is ε-close to the boundary ∂D (we call this event stop due to space constraint) or T n is ε-close to 0 (stop due to time constraint). Then it seems quite obvious that stops due to time constraint are more likely to occur when t becomes small (see the proportion in Figure 6 left).

Let us now comment the algorithm efficiency by focusing our attention on the number of steps. The distribution of this random number depends on several parameters: the dimension d, the parameter ε and finally the choice of (t, x) (see histogram Figure 6 -right -for a particular choice of parameters). We have pointed out an upper bound for the average number of steps in Theorem 4.2. The numerics permit to present different curves illustrating all the dependences: the logarithm growth with respect to the parameter ε, the surprising behavior when the space position x varies and the influence of the dimension (Figures 7 and8). Let us notice that this algorithm is especially efficient (see the small number of steps) even in high dimensions.

Half sphere

All the studies developed in the hypercube case can also be considered for the half sphere. We introduce particular boundary conditions:

f (t, x) = (1 + cos(2πt)) x , ∀x ∈ D, ∀t ≥ 0, (5.4)
with f 0 (x) = f (0, x). Similarly as above, we present:

• the approximated solution u ε N (t, x) for the default value x = (0.5, 0 . . . , 0) and for t varying in the interval [0, 4] (Figure 9),

• the distribution of the Monte Carlo underlying variable Z ε for different values of t and different dimension values d (Figures 10 and11),

• different curves illustrating the influence of the parameter ε, the starting position x and the dimension d on the averaged number of steps (Figures 12 and13)

B Improvements for the classical random walk on spheres

In this section, we focus our attention to the classical random walk on spheres. We consider an 0-thick domain D, see the definition developed in (4.3), and the Euclidean distance to the boundary d(x) = δ(x, ∂D). The random walk is then defined as follows:

• we start with an initial condition X 0 and fix two parameters ε > 0 and β ∈]0, 1[.

• While d(X n ) > ε, we construct

X n+1 = X n + βd(X n )γ n (B.1)
where (γ n ) stands for a sequence of independent random variables uniformly distributed on the unit sphere in R d .

We adapt here several results of [START_REF] Binder | The rate of convergence of the walk on spheres algorithm[END_REF] to our particular situation. Let us recall that U is the energy function defined by (4.10) which is based on the set of measures M, defined by (4.5), and on the Riesz potential. Since D is a 0-thick domain, the following Lemma holds.

Lemma B.1. There exist two constants δ > 0 and η > 0, such that: for any y ∈ D (we define x the closest point of y belonging to the boundary) and any measure µ ∈ M, we have: This lemma, which is quite general and is not directly linked to the random walk, has an important consequence on it (for the proof of Lemma B.1, see [START_REF] Binder | The rate of convergence of the walk on spheres algorithm[END_REF]). In [START_REF] Binder | The rate of convergence of the walk on spheres algorithm[END_REF], the authors consider a general random walk defined by (B.1). They prove that there exist an interger k and a constant L > 0 such that E[U n+k -U n |U n ] > L. Here we adapt the proof by introducing a particular condition on the parameter β which permits in fact to set k = 1.

Proof. Let us consider X n . Due to the weak compactness of the set of measures M (see Remark 4.3), there exists a measure µ ∈ M such that

U n = U(X n ) = U µ (X n ).
For this particular measure, either the first or the second point of the previous lemma are satisfied.

Step 1. Let us assume that the first point is satisfied that is, U(z) > U µ (z) + 1 when zx < δd(y) and d(z) > δ/4 d(y). Since U µ (X n ) is a submartingale, we get

E[U n+1 -U n |X n ] = E[U n+1 -U µ n+1 |X n ] + E[U µ n+1 -U µ n |X n ] ≥ E[U n+1 -U µ n+1 |X n ] ≥ P X n+1 -x n < δd(X n ), d(X n+1 ) > (δ/4)d(X n ) X n ,
where x n is the closest point of X n on the boundary ∂D. We denote by u n = xn-Xn d(Xn) which belongs to the unit sphere. Using the definition of the random walk and the particular choice of the parameter β, we get immediately d(X n+1 ) > (1β)d(X n ) > δ/4d(X n ), and

X n+1 -x n = d(X n ) u n -βγ n ≤ d(X n ) (1 -β) γ n + γ n -u n = d(X n )(1 -β + γ n -u n ) < d(X n )(δ/2 + γ n -u n ). (B.2)
Let us recall that u n is a unit vector. Then we define the set Γ un of points u belonging to the unit sphere of dimension d such that uu n < δ/2. Let us just note that Γ un is a non empty open set. We observe that P(γ n ∈ Γ un ) =: p > 0 for any n ≥ 0 and does not depend on u n due to rotational invariant of the distribution. Furthermore, for any γ n ∈ Γ un , (B.2) implies that X n+1x n < δd(X n ). Therefore

E[U n+1 -U n |X n ] ≥ P X n+1 -x n < δd(X n ), d(X n+1 ) > (δ/4)d(X n ) X n ≥ P(γ n ∈ Γ un ) = p > 0.
Step dr outside the support of the measure µ (consequently U µ is a C 2 -function in the domain D), then, for any y satisfying y -X n ≤ βd(X n ), we get ∆U µ (y) ≥ 2(d + 2)µ(B(X n , 2d(X n )))

∞ (2+β)d(Xn ) dr r d+3 = 2µ(B(X n , 2d(X n )) ((2 + β)d(X n )) d+2 ≥ 2η (2 + β) d+2 d(X n ) -2 .

2 0 4 DFigure 1 :

 241 Figure 1: WOS algorithm

Proposition 2 . 4 .

 24 Let d > 4 and let h be a bounded function, defined on an open set A and satisfying the mean value property

Proposition 3 . 3 .

 33 (uniqueness) Set T > 0. Let u be a C 1,2 ([0, T ] × D)-function satisfying the Initial-Boundary Value Problem (1.1) and continuous with respect to both variables on [0, T ] × D. Then u is unique and given by the expression

(4. 2 ) 4 . 1 .

 241 Proposition Let us assume that the Initial-boundary Value Problem (1.1) admits an unique C 1,2 ([0, T ] × D)-solution u, defined by (3.

Figure 3 :

 3 Figure 3: u ε N (t, x) versus t for N = 1 000, ε = 0.001, d = 3.

Figure 4 :

 4 Figure 4: Histogram of the distribution of 10 000 random variables Z ε for various values of t: t = 0.05 (left), t = 0.1 (right) ε = 0.001, d = 3.

Figure 5 :

 5 Figure 5: Histogram of Z ε for d = 3 (left) and d = 10 (right), t = 4, ε = 0.001.

Figure 6 :

 6 Figure 6: Proportion of stops due to space constraint versus t for 10 000 trials, ε = 0.001, d = 3 (left); histogram of the number of steps t = 1, ε = 0.001, d = 3.

Figure 7 :

 7 Figure 7: Averaged number of steps versus n for ε = 0.5 n and x the center of the hypercube (left), averaged number of steps versus u for x = (u, u, u) and ε = 0.001. In both situations: 10 000 trials, d = 3, t = 1.

Figure 8 :

 8 Figure 8: Averaged number of steps versus the dimension d, ε = 0.001, 10 000 trials, t = 1.

Figure 9 :

 9 Figure 9: u ε N (t, x) versus t for N = 1 000, ε = 0.001, d = 3.

Figure 10 :

 10 Figure 10: Histogram of the distribution of 10 000 random variables Z ε for various values of t: t = 0.01 (left), t = 0.05 (right) ε = 0.001, d = 3.

Figure 11 :

 11 Figure 11: Histogram of Z ε for d = 3 and t = 1 (left), for d = 10 and t = 0.01 (right), both with ε = 0.001.

Figure 12 :

 12 Figure12: Averaged number of steps versus n for ε = 0.5 n and x = (0.5, 0, 0) (left), averaged number of steps versus u for x = (u, 0, 0) and ε = 0.001. In both situations: 10 000 trials, d = 3, t = 1.

Figure 13 :

 13 Figure 13: Averaged number of steps versus the dimension d, ε = 0.001, 10 000 trials, t = 1.

1 .

 1 either U(z) > U µ (z) + 1 whenever zx < δd(y) and d(z) > δ/4 d(y) 2. or µ(B(y, 2d(y))) ≥ ηd(y) d .

Proposition B. 2 .

 2 Let us consider δ the constant of Lemma B.1 and the random walk (X n ) n≥0 defined by (B.1) with β ∈]1δ/2, 1δ/4[. Then there exists a constant L > 0 such that the sequence U n := U(X n ) satisfies E[U n+1 -U n |U n ] > L, ∀n ≥ 0.

2 .

 2 The second case concerns the condition µ(B(y, 2d(y))) ≥ ηd(y) d .By the Green formula, for a C 2 -smooth function h,E[h(X n+1 )|X n ]h(X n ) = S(Xn,βd(Xn))h(y)dσ(y)h(X n )

  If h satisfies the mean value property (2.2) and is a C 1,2 -function for any (t, x) ∈ A and α > 0 such that [t, t + α] × B(x, 2 αd/(2e)) ⊂ A, then h is a reverse temperature in A.

		The reverse statement of the preceding result can also be proved. The first step consists
	in the following
	Proposition 2.3. Proof. Let us consider the function H : [0,	√ α] → R defined by
						H(r) = h(t + r 2 s, x + 2rψ d (s)y),
	for any (s, y) ∈ [0, 1] × S d . Using the Taylor expansion, we get
	√	2αψ d (s) where ψ d is defined by (2.3). We get		1[, and observe that ψ a,d (αs) =
		h(t, x) =	1 σ(S d )	0	1	S d	h(t + αs, x +	√	2ψ a,d (αs)y)	1 2as	ψ d a,d (αs) dσ(y) ds
		=	1 σ(S d )	0	1	S d	h(t + αs, x + 2	√	αψ d (s)y)	2 d/2 α d/2 2as	ψ d d (s) dσ(y) ds.
	Using both the explicit expression of α d/2 and the classical formula σ(S d ) = 2π d/2 /Γ(d/2)
	leads to (2.2).

  .15) By Proposition 2.2, Proposition 2.3 and Proposition 2.4, we obtain the following: Theorem 2.6. 1. Let A ⊂ R + × R d be a non empty open set. If a function h is a temperature in A then it has the following mean value property: where S d is the d-dimensional sphere of radius 1, σ is the Lebesgue measure on S d and ψ d is defined in (2.3). Equation (2.16) is satisfied for any α > 0 such that [tα, t] × B(x, 2 αd/(2e)) ⊂ A. 2. If h satisfies the mean value property (2.16) and is a C 1,2 -function for any (t, x) ∈ A and α > 0 such that [tα, t] × B(x, 2 αd/(2e)) ⊂ A then h is a temperature in A. 3. For d > 4, any bounded function h is a temperature iff it satisfies the mean value property (2.16).

	h(t, x) =	1 2π d/2	(s,y)∈]0,1[×S d	1 s	h(t -αs, x + 2	√ αψ d (s)y) ψ d d (s) dσ(y) ds,	(2.16)

A Technical results

We first start with Jensen's inequality: Lemma A.1. Let X and Y be two random variables and A a σ-algebra, then

Proof. We shall just prove the first inequality. The proof of the second one is similar. We get

Combining both inequalities leads to the result.

Let us now present properties concerning a particular probability distribution arising in the random walk on moving spheres.

where the function ψ d is defined by (2.3) and R is a random variable with the following probability density function:

Then W has its two first moments (denoted by M 1 and M 2 ) bounded.

Proof. Let us first note that W is a non-negative random variable, since R and 2e

.

Let us observe that f b (s) tends to 0 as s → 0 and in a neighborhood of 1,

.

In a neigborhood of 0, we have f a (s) ∼ C 2 s d/2 (log s) d/2+1 , in a neighborhood of 1, we observe

and finally in a neigborhood of 1/e, f a (s) ∼ C 4 log 2 |s -1 e |. We deduce that f a is integrable on the whole interval [0, 1] and E[W 2 a ] < ∞.

Applying the previous results to the particular regular function h = U µ , we deduce:


Maria Laura

Delle Monache

Paola Goatin

STABILITY ESTIMATES FOR SCALAR CONSERVATION LAWS WITH MOVING FLUX CONSTRAINTS

Keywords: 2010 Mathematics Subject Classification. Primary: 35L65; Secondary: 35B35 Scalar conservation laws, moving flux constraints, PDE-ODE system

We study well-posedness of scalar conservation laws with moving flux constraints. In particular, we show the Lipschitz continuous dependence of BV solutions with respect to the initial data and the constraint trajectory. Applications to traffic flow theory are detailed.

1. Introduction. Motivated by the modeling of moving bottlenecks in traffic flow, we consider the Cauchy problem for a scalar conservation law with moving flux constraint

∂ t ρ + ∂ x f (ρ) = 0, (t, x) ∈ R + × R, (1a)
ρ(0, x) = ρ 0 (x), x ∈ R, (1b)
f (ρ(t, y(t))) -ẏ(t)ρ(t, y(t)) ≤ F α ( ẏ(t)),

t ∈ R + , (1c)
where t → y(t) is a given trajectory, starting from y(0) = y 0 . Systems of the form [START_REF] Andreianov | Riemann problems with nonlocal point constraints and capacity drop[END_REF] arise in the modeling of moving bottlenecks in vehicular traffic [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF][START_REF] Villa | Moving bottlenecks for the Aw-Rascle-Zhang traffic flow model[END_REF]: ρ = ρ(t, x) ∈ [0, R] is the scalar conserved quantity and represents the traffic density, whose maximum attainable value is R. The flux function f : [0, R] → R + is assumed to be strictly concave and such that f (0) = f (R) = 0. The time-dependent variable y denotes the constraint position. In the present paper we consider a weakly coupled PDE-ODE system, in the sense that we assume that the constraint trajectory is given, and it does not depend on the solution of (1a).

Let us detail the meaning of inequality (1c). A moving flux constraint located at x = y(t) acts as an obstacle, thus hindering the flow as expressed by the unilateral constraint (1c). There, α ∈ ]0, 1[ is the dimensionless reduction rate of the road capacity (the maximal allowed density) at the bottleneck position. The inequality (1c) is derived by studying the problem in the constraint reference frame, i.e., setting ρ(t, x) := ρ(t, x + y(t)) and rewriting the conservation law (1a) as

∂ t ρ + ∂ x F (ρ) = 0, F (ρ) = f (ρ) -ẏ ρ. (2)
In fact, let f α : [0, αR] → R + be the rescaled flux function describing the reduced flow at x = y(t), i.e. f α (ρ) = αf (ρ/α), and ρ α ∈ ]0, αR[ such that

F α (ρ α ) = f α (ρ α ) -ẏρ α = 0 ⇔ f α (ρ α ) = ẏ with F α (ρ) = f α (ρ) -ẏρ, see Figure 1. Then setting F α ( ẏ) := F α (ρ α ) = f α (ρ α ) -ẏρ α
and imposing that in the obstacle reference frame the flux F is less than the maximum value of the flux of the reduced flow, one gets (1c). Notice that the inequality (1c) is always satisfied if

ẏ(t) = f (ρ(t, y(t)±)) ρ(t, y(t)±) ,
since the left-hand side is 0. Moreover, it is well defined even if ρ has a jump at y(t), because of the Rankine-Hugoniot conditions. Problem (2), (1b), (1c), can therefore be recast in the framework of conservation laws with fix local constraint, first introduced in [START_REF] Colombo | A well posed conservation law with a variable unilateral constraint[END_REF], then developed in [START_REF] Andreianov | Finite volume schemes for locally constrained conservation laws[END_REF][START_REF] Andreianov | Riemann problems with nonlocal point constraints and capacity drop[END_REF] for scalar equations and extended in [START_REF] Garavello | The Aw-Rascle traffic model with locally constrained flow[END_REF][START_REF] Andreianov | A second-order model for vehicular traffics with local point constraints on the flow[END_REF][START_REF] Garavello | The Cauchy problem for the Aw-Rascle-Zhang traffic model with locally constrained flow[END_REF] to systems. Following [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF]Definition 4.1] and [5, Definition 1 and 2], solutions of (1) are defined as follows.

ρ αR R f (ρ) 0 ρ( ẏ) ρ α ρ( ẏ) ẏ (a) Fixed reference frame ẏρ R F α f (ρ) -ẏρ 0 ρα ρ (b) Obstacle reference frame
Definition 1.1. Let y ∈ W 1,1 (R + ; R) and ρ 0 ∈ L 1 ∩ L ∞ (R; [0, R]) be given. A function ρ ∈ C 0 R + ; L 1 ∩ L ∞ (R; [0, R]) is a solution to (1) if 1. ρ satisfies Kružhkov entropy conditions [14] on (R + × R) \ {(t, y(t)) : t ∈ R + },
i.e. for every k ∈ [0, 1] and for all ϕ ∈ C 1 c (R 2 ; R + ) and ϕ(t, y(t)) = 0, t > 0,

R + R (|ρ -k|∂ t ϕ + sgn(ρ -k) (f (ρ) -f (k)) ∂ x ϕ) dx dt + R |ρ o -k|ϕ(0, x) dx ≥ 0 ; (3a)
2. for a. e. t > 0 the left and right traces of ρ at x = y(t) satisfy

(ρ(t, y(t)-), ρ(t, y(t)+)) ∈ G α ( ẏ(t)). ( 3b
)
The set G α ( ẏ) in (13b) is defined as follows, see [START_REF] Andreianov | A theory of L 1 -dissipative solvers for scalar conservation laws with discontinuous flux[END_REF][START_REF] Andreianov | Finite volume schemes for locally constrained conservation laws[END_REF][START_REF] Andreianov | Well-posedness for a onedimensional fluid-particle interaction model[END_REF].

Definition 1.2. The admissibility germ G α ( ẏ) ⊂ [0, R] 2 for (1a), (1c) is the union G α ( ẏ) := G 1 ( ẏ) ∪ G 2 ( ẏ) ∪ G 3 ( ẏ),
where

G 1 ( ẏ) := (c L , c R ) ∈ [0, R] 2 : c L > c R , f (c L ) -ẏc L = f (c R ) -ẏc R = F α ( ẏ) , G 2 ( ẏ) := (c, c) ∈ [0, R] 2 : f (c) -ẏc ≤ F α ( ẏ) , G 3 ( ẏ) := (c L , c R ) ∈ [0, R] 2 : c L < c R , f (c L ) -ẏc L = f (c R ) -ẏc R ≤ F α ( ẏ) .
The equivalence between Definition 1.1 and [11, Definition 4.1] can be proved as in [START_REF] Andreianov | Finite volume schemes for locally constrained conservation laws[END_REF]Proposition 2.6].

Systems of the type (1) arise in the modeling of moving bottlenecks in traffic flows, see [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF][START_REF] Villa | Moving bottlenecks for the Aw-Rascle-Zhang traffic flow model[END_REF] and Section 3 below, where they are coupled with an ODE depending on the downstream traffic velocity and describing the trajectory of a slow moving vehicle (a bus or a truck) acting as a bottleneck. This paper is a first step towards establishing well-posedness for the strongly coupled models [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF][START_REF] Villa | Moving bottlenecks for the Aw-Rascle-Zhang traffic flow model[END_REF]. Section 2 presents the main result, stating the L 1 Lipschitz continuous dependence of solutions of (1) from the initial data and the constraint trajectory. Section 3 describes in details the related traffic flow model with moving bottleneck. Technical proof details are deferred to Appendix A.

2. Lipschitz continuous dependence of BV solutions with respect to ẏ. Let us fix y ∈ W 1,1 (R + ) and ρ 0 ∈ BV(R; [0, R]), and let ρ ∈ L ∞ (R + ×R; [0, R]) be a solution of (1) in the sense of Definition 1.1. Applying the coordinate transformation ρ(t, x) := ρ(t, x + y(t)), ρ is a weak entropy solution of the problem

∂ t ρ + ∂ x F (ρ) = 0, (t, x) ∈ R + × R, (4a)
ρ(0, x) = ρ 0 (x + y 0 ), x ∈ R, (4b)
F (ρ(t, 0)) ≤ F α ( ẏ(t)), t ∈ R + , (4c)
in the sense of [3, Proposition 2.6(A)] where we have set F (ρ) := f (ρ) -ẏ ρ. Existence and uniqueness for (4) are proved in [START_REF] Andreianov | Finite volume schemes for locally constrained conservation laws[END_REF][START_REF] Colombo | A well posed conservation law with a variable unilateral constraint[END_REF]. In particular, we have that for every k ∈ [0, 1] and for all ϕ ∈ C 1 c (R 2 ; R + ) such that ϕ(t, 0) = 0, t > 0,

R + R (|ρ -k|∂ t ϕ + Φ( ẏ; ρ, k)∂ x ϕ) dx dt + R |ρ o (x + y 0 ) -k|ϕ(0, x) dx ≥ 0 ,
where we have set

Φ( ẏ; a, b) := sgn(a -b) (f (a) -f (b)) -ẏ|a -b|, for a, b ∈ R,

and

(ρ(t, 0-), ρ(t, 0+)) ∈ G α ( ẏ(t)) for a. e. t > 0, see [START_REF] Andreianov | Finite volume schemes for locally constrained conservation laws[END_REF]Proposition 2.6]. Moreover, we note that, since α < 1 and f is strictly concave, we have

F α ( ẏ) = f α (ρ α ) -ẏρ α < f (ρ α ) -ẏρ α = max ρ∈[0,R] F (ρ). ( 5
)
Remark 1. Solution to (4) are not in BV in general, since y ∈ W 1,1 (R + ) and thus [START_REF] Andreianov | Well-posedness for a onedimensional fluid-particle interaction model[END_REF].

ẏ ∈ BV(R + ). Yet, if ẏ ∈ BV([0, t]; R), ρ(t, •) ∈ BV(R; [0, R]) due to
In any case, the left and right traces at x = 0 do exist, see [START_REF] Andreianov | Finite volume schemes for locally constrained conservation laws[END_REF]Section 2].

We compare solutions of (4) corresponding to different constraint trajectories y and z.

Theorem 2.1. Assume y, z ∈ W 1,1 (R + ) with 0 ≤ ẏ, ż < f (0), ρ0 , σ0 ∈ L ∞ (R, [0, R]) and ρ0 -σ0 ∈ L 1 (R). Let ρ, σ ∈ L ∞ (R + × R; [0, R]
) be the solutions of (4) corresponding respectively to y, ρ0 and z, σ0 . Moreover we assume that TV(ρ(t, •)) ≤ C for all t > 0. Then we have

ρ(t, •) -σ(t, •) L 1 (R) ≤ ρ0 -σ0 L 1 (R) + (C + 2R) ẏ -ż L 1 ([0,t]) . (6)
Proof. The two solutions ρ, σ satisfy

∂ t |ρ -k| + ∂ x Φ( ẏ; ρ, k) ≤ 0, (7)
∂ t |σ -k| + ∂ x Φ( ż; σ, k) ≤ 0, (8)
in D (R + × R * ) (where we have noted R * = R \ {0}). Following the proofs of [5, Lemma 15] and [6, Theorem 3.1] we observe that

∂ t |ρ -k| + ∂ x Φ( ẏ; ρ, k) = ∂ t |ρ -k| + ∂ x Φ( ż; ρ, k) + ∂ x Φ( ẏ; ρ, k) -∂ x Φ( ż; ρ, k) therefore ∂ t |ρ -k| + ∂ x Φ( ż; ρ, k) ≤ ∂ x Φ( ż; ρ, k) -∂ x Φ( ẏ; ρ, k) ≤ | ẏ -ż||∂ x ρ|. (9)
Applying the classical Kružkov doubling of variables technique [START_REF] Kružkov | First order quasilinear equations with several independent variables[END_REF], with a test function ψ ∈ C 1 c (R 2 ; R + ) such that ψ(t, 0) = 0, we deduce the following Kato inequality where v(ρ) = V (1 -ρ/R) is the mean traffic speed, V being the maximal velocity allowed on the road, problem (1) can be used to describe the situation of a moving bottleneck along a road, see [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF]. In this case, we get f α (ρ) = V ρ 1 -ρ αR and

R + R (|ρ -σ|∂ t ψ + Φ( ż; ρ, σ)∂ x ψ) dx dt + R |ρ 0 -σ0 |ψ(0, x) dx + C ψ ∞ R + | ẏ -ż|dt ≥ 0 , We now choose the test function ψ(t, x) = θ ε (x)ξ(t, x), where ξ ∈ C 1 c (R 2 ; R + ) is an approximation of the characteristic function of the trapezoid {(s, x) : |x| ≤ M + L(t -s), 0 ≤ s ≤ t}
ρ α = αR 2 1 -ẏ V , so that F α ( ẏ) = αR 4V (V -ẏ) 2 .
Let us suppose that a slow and large vehicle, like for example a bus or a truck moves on the road. The slow vehicle, that in the following we will refer as bus, reduces the road capacity and moves with a trajectory given by the following ODE:

ẏ(t) = ω(ρ(t, y(t)+)), t ∈ R + , y(0) = y 0 , (11)
where the velocity of the bus is given by the following traffic density dependent function (see Figure 2)

ω(ρ) = V b if ρ ≤ ρ * . = R(1 -V b /V ), v(ρ) otherwise. ( 12
)
This means that if the traffic is not too congested, the bus moves at its own maximal speed V b < V . When the surrounding traffic density becomes too high, the bus adapts its velocity accordingly. In particular, it is not possible for the bus to overtake the cars. Solutions of the coupled system (1), [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF] for general initial data are defined as follows.

Definition 3.1. A couple (ρ, y) ∈ C 0 R + ; L 1 ∩ BV(R; [0, 1]) × W 1,1 (R + ; R) is a solution to (1) if 1. ρ satisfies Kružhkov entropy conditions [14] on (R + × R) \ {(t, y(t)) : t ∈ R + },
i.e. for every k ∈ [0, 1] and for all ϕ ∈ C 1 c (R 2 ; R + ) and ϕ(t, y(t)) = 0, t > 0,

R + R (|ρ -k|∂ t ϕ + sgn(ρ -k) (f (ρ) -f (k)) ∂ x ϕ) dx dt + R |ρ o -k|ϕ(0, x) dx ≥ 0 ; (13a
The proof of existence of solutions for the general Cauchy problem (1) strongly coupled with the bus trajectory [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF] with BV initial data can be found in [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF]. For completeness, we recall the definition the solution of a Riemann problem, as given in [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF]. Let us consider a Riemann type initial datum

ρ 0 (x) = ρ L if x < 0, ρ R if x > 0, y 0 = 0. ( 14
)
Denote by R the standard (i.e., without the constraint (1c)) Riemann solver for (1a)-(1b)-( 14), i.e., the (right continuous) map (t, x) → R(ρ L , ρ R )(x/t) given by the standard weak entropy solution, see for instance [START_REF] Bressan | Hyperbolic systems of conservation laws[END_REF]Chapter 5]. Moreover, assume that ẏ is constant and let ρ = ρ( ẏ) and ρ = ρ( ẏ), with ρ ≤ ρ, be the intersection points of the flux function f (ρ) with the line f α (ρ α ) + ẏ(ρ -ρ α ) (see Figure 1(a)): 1), ( 11), ( 14) is defined as follows.

ρ( ẏ) = R 2V 1 - √ 1 -α (V -ẏ), ρ( ẏ) = R 2V 1 + √ 1 -α (V -ẏ). ( 15
) Definition 3.2. The constrained Riemann solver R α : [0, R] 2 → L 1 loc (R; [0, R]) corresponding to (
1. If f (R(ρ L , ρ R )(V b )) > F α + V b R(ρ L , ρ R )(V b ), then R α (ρ L , ρ R )(x/t) = R(ρ L , ρ(V b ))(x/t) if x < V b t, R(ρ(V b ), ρ R )(x/t) if x ≥ V b t, and
y(t) = V b t. 2. If V b R(ρ L , ρ R )(V b ) ≤ f (R(ρ L , ρ R )(V b )) ≤ F α + V b R(ρ L , ρ R )(V b ), then R α (ρ L , ρ R ) = R(ρ L , ρ R ) and y(t) = V b t. 3. If f (R(ρ L , ρ R )(V b )) < V b R(ρ L , ρ R )(V b ), i.e., v(R(ρ L , ρ R )(V b )) < V b then R α (ρ L , ρ R ) = R(ρ L , ρ R ) and y(t) = v(ρ R )t.
Note that, when the constraint is enforced (point 1. in the above definition), a non-classical shock arises between ρ(V b ) and ρ(V b ), which satisfies the Rankine-Hugoniot condition but violates the Lax entropy condition, see Figure 3 for an example.

Remark 2. Unfortunately, no result about the Lipschitz continuous dependence of the solution y = y(t) of ( 11) from the solution ρ = ρ(t, x) of ( 1) is known at present. Related results [START_REF] Bressan | Uniqueness for discontinuous ODE and conservation laws[END_REF][START_REF] Colombo | A Hölder continuous ODE related to traffic flow[END_REF] concerning uniqueness and continuous dependence for ODEs of the form [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF] hold under hypothesis on the speed function ω that are not satisfied by [START_REF] Garavello | The Aw-Rascle traffic model with locally constrained flow[END_REF].

Appendix A. Proof of [START_REF] Colombo | A Hölder continuous ODE related to traffic flow[END_REF].

Lemma A.1. For any (ρ -, ρ + ) ∈ G α ( ẏ) and (σ -, σ + ) ∈ G α ( ż), it holds Φ( ż; ρ + , σ + ) -Φ( ż; ρ -, σ -) ≤ 2R| ẏ -ż|.
Proof. For simplicity, we consider the case

f (ρ) = V ρ(1 -ρ/R). First of all, let us remark that Φ( ż; ρ + , σ + ) -Φ( ż; ρ -, σ -) (16)
= sgn(ρ

+ -σ + ) (f (ρ + ) -żρ + -f (σ + ) + żσ + ) -sgn(ρ --σ -) (f (ρ -) -żρ --f (σ -) + żσ -) =(λ(ρ + , σ + ) -ż)|ρ + -σ + | -(λ(ρ -, σ -) -ż)|ρ --σ -|, (17)
where

λ(ρ, σ) = f (ρ) -f (σ) ρ -σ = V 1 - ρ + σ R .
Without loss of generality, we can assume ẏ < ż. Therefore we get

F α ( ẏ) -F α ( ż) = αR 4V (2V -ż -ẏ)( ż -ẏ) > 0. ( 18
)
We distinguish the following cases: 1. (ρ -, ρ + ) ∈ G 1 ( ẏ): we observe that

ρ -= ρ, ρ + = ρ and f (ρ -) -ẏρ -= f (ρ + ) -ẏρ + = F α ( ẏ).
Depending on the values of σ -, σ + , we different situations can occur as shown in Figure 4: 1.1 (σ -, σ + ) ∈ G 1 ( ż): in this case σ -= σ, σ + = σ, as shown in Figure 4a and

f (σ -) -żσ -= f (σ + ) - żσ + = F α ( ż), therefore (16) = (F α ( ẏ) + ẏρ + -żρ + -F α ( ż)) -(F α ( ẏ) + ẏρ --żρ --F α ( ż)) = (ρ --ρ + )( ż -ẏ) ≤ R| ẏ -ż|. 1.2 (σ -, σ + ) ∈ G 2 ( ż): we set σ := σ -= σ + and f (σ) -żσ =: F (σ) ≤ F α ( ż).
The following cases can occur:

• 0 ≤ σ ≤ σ: (16) = (F α ( ẏ) + ẏρ + -żρ + -F (σ)) -(F α ( ẏ) + ẏρ --żρ --F (σ)) = (ρ --ρ + )( ż -ẏ) ≤ R| ẏ -ż|. ρ α ρ * 1 f (ρ) 0 ρα ρ α ρα V b ρ V b x t ρ ρα ρα ρ (
ρ α ρ * 1 f (ρ) 0 ρα ρ α ρα V b ρ L ρ R V b x t ρ L ρ R (b)
Example of solution corresponding to case 2. of Definition 3.2 for 0 < ρL < ρα and ρα < ρR < ρ * .

ρ α ρ * 1 f (ρ) 0 ρα ρ α ρα V b ρ L ρ R v(ρ R ) x t ρ R ρ L (c)
σ -ρ + σ + ρ - ẏ ż (a) Case 1.1 ρ R f (ρ) σ ρ + σ ρ -σ ẏ ż (b) Case 1.2 ρ R f (ρ) σ -σ ρ σ ρ σ + ẏ ż(c) Case 1.3 Figure 4. Case 1 • σ ≤ σ ≤ ρ, see Figure 4b: (16) = -(F α ( ẏ) + ẏρ + -żρ + -F (σ)) -(F α ( ẏ) + ẏρ --żρ --F (σ)) = 2F (σ) -2F α ( ẏ) + (ρ + + ρ -)( ż -ẏ) ≤ 2F α ( ż) -2F α ( ẏ) + 2R( ż -ẏ) = - αR 2V (2V -ż -ẏ)( ż -ẏ) + 2R( ż -ẏ) ≤ 2R| ẏ -ż|.
• ρ ≤ σ ≤ R:

(16) = -(F α ( ẏ) + ẏρ + -żρ + -F (σ)) + (F α ( ẏ) + ẏρ --żρ --F (σ)) = (ρ + -ρ -)( ż -ẏ) ≤ 0. 1.3 (σ -, σ + ) ∈ G 3 ( ż): we set f (σ -) -żσ -= f (σ + ) -żσ + =: F (σ) ≤ F α ( ż).
We observe that sgn(ρ + -σ + ) < 0 and sgn(ρ --σ -) > 0, see Figure 4c. Therefore

(16) = -(F α ( ẏ) + ẏρ + -żρ + -F (σ)) -(F α ( ẏ) + ẏρ --żρ --F (σ)) = 2F (σ) -2F α ( ẏ) + (ρ + + ρ -)( ż -ẏ) ≤ 2F α ( ż) -2F α ( ẏ) + 2R( ż -ẏ) = - αR 2V (2V -ż -ẏ)( ż -ẏ) + 2R( ż -ẏ) ≤ 2R| ẏ -ż|. 2. (ρ -, ρ + ) ∈ G 2 ( ẏ): we set ρ := ρ -= ρ + and f (ρ) -ẏρ =: F (ρ) ≤ F α ( ẏ).
as illustrated in Figure 5.

ρ R f (ρ) 0 σ ρ σ ρ ρ ẏ ż (a) Case 2.1 ρ R f (ρ) 0 σ ρ σ ρ ρ σ ẏ ż (b) Case 2.2 ρ R f (ρ) σ -σ ρ σ ρ σ + ρ ẏ ż(c) Case 2.3 Figure 5. Case 2 2.1 (σ -, σ + ) ∈ G 1 ( ż): • 0 ≤ ρ ≤ σ: (16) = -(F (ρ) + ẏρ -żρ -F α ( ż)) + (F (ρ) + ẏρ -żρ -F α ( ż)) = 0.

STABILITY OF SCALAR CONSERVATION LAWS WITH MOVING FLUX CONSTRAINTS 11

• σ ≤ ρ ≤ ρ, this case is displayed in Figure 5a:

(16) = (F (ρ) + ẏρ -żρ -F α ( ż)) + (F (ρ) + ẏρ -żρ -F α ( ż)) = 2F (ρ) -2F α ( ż) + 2ρ( ẏ -ż) ≤ 2 (F α ( ẏ) -F α ( ż))
≤ αR| ẏ -ż|.

• ρ ≤ ρ ≤ R:

(16) = (F (ρ) + ẏρ -żρ -F α ( ż)) -(F (ρ) + ẏρ -żρ -F α ( ż)) = 0. 2.2 (σ -, σ + ) ∈ G 2 ( ż), illustrated in Figure 5b: We observe that sgn(ρ + -σ + ) = sgn(ρ --σ -) = sgn(ρ -σ), therefore (16) = sgn(ρ -σ) [(F (ρ) + ẏρ -żρ -F (σ)) -(F (ρ) + ẏρ -żρ -F (σ))] = 0. 2.3 (σ -, σ + ) ∈ G 3 ( ż) shown in Figure 5c: If sgn(ρ -σ + ) = sgn(ρ -σ -), we get (16) = sgn(ρ -σ + ) [(F (ρ) + ẏρ -żρ -F (σ)) -(F (ρ) + ẏρ -żρ -F (σ))] = 0.
Otherwise, we have that σ -≤ ρ ≤ ρ or ρ ≤ ρ ≤ σ + . In this case

λ(ρ + , σ + ) -ż = V 1 - ρ + σ + R -V 1 - σ + + σ - R = V R (σ --ρ) ≤ 0, λ(ρ -, σ -) -ż = V 1 - ρ + σ - R -V 1 - σ + + σ - R = V R (σ + -ρ) ≥ 0.
Therefore, (16) ≤ 0 by (17). 3. (ρ -, ρ + ) ∈ G 3 ( ẏ): we set f (ρ -) -ẏρ -= f (ρ + ) -ẏρ + =: F (ρ) ≤ F α ( ẏ).

See Figure 6, for a graphical representation.

3.1 (σ -, σ + ) ∈ G 1 ( ż), see Figure 6a.

≤ (ρ + -ρ -)( ż -ẏ) ≤ R| ẏ -ż|.
Otherwise as shown in Figure 6b, we have that ρ -≤ σ ≤ σ or σ ≤ σ ≤ ρ + . In this case

λ(ρ + , σ + ) -ż ≤ λ(ρ + , σ + ) - ẏ = V 1 - ρ + + σ R -V 1 - ρ + + ρ - R = V R (ρ --σ) ≤ 0.
Moreover, we observe that λ(ρ -, σ -) > ẏ. Indeed We observe that one of the following relations must hold

λ(ρ -, σ -) -ẏ = V 1 - ρ -+ σ R -V 1 - ρ + + ρ - R = V R (ρ + -σ) ≥ 0. ρ R f (ρ) σ + ρ - σ - ρ + ẏ ż (
ρ -≤ σ -< σ + ≤ ρ + , σ -≤ ρ -< σ + ≤ ρ + , σ -≤ ρ -< ρ + ≤ σ + .
For an example see Figure 6c. Therefore

λ(ρ + , σ + ) -ż = V 1 - ρ + + σ + R -V 1 - σ + + σ - R = V R (σ --ρ + ) ≤ 0, λ(ρ -, σ -) -ż = V 1 - ρ -+ σ - R -V 1 - σ + + σ - R = V R (σ + -ρ -) ≥ 0,
hence (16) ≤ 0 by (17).

Figure 1 .

 1 Figure 1. Graphical representation of the constraint action in the fixed (left) and moving (right) reference frames.

( 0 |

 0 where L ≥ |V -ż|)and θ ε a smooth approximation of x → min{|x|/ε, 1}. Following the proof of[START_REF] Colombo | A well posed conservation law with a variable unilateral constraint[END_REF] Proposition 4.4] and letting ε 0 we getM -M |ρ(t, x) -σ(t, x)|dx ≤ M +Lt -M -Lt |ρ 0 (x) -σ0 (x)|dx + C t 0 | ẏ(s) -ż(s)|ds + t 0 (Φ( ż; ρ(t, 0+), σ(t, 0+)) -Φ( ż; ρ(t, 0-), σ(t, 0-))) ds.STABILITY OF SCALAR CONSERVATION LAWS WITH MOVING FLUX CONSTRAINTS 5By Lemma A.1, the last integrand can be bounded by t 0 (Φ( ż; ρ(t, 0+), σ(t, 0+)) -Φ( ż; ρ(t, 0-), σ(t, 0-))) ds ≤ 2R t ẏ(s) -ż(s)|ds.(10) Letting M → ∞, we recover[START_REF] Bouchut | Kružkov's estimates for scalar conservation laws revisited[END_REF].

3 .

 3 Application to traffic modeling. Setting f (ρ) := ρv(ρ),

Figure 2 .

 2 Figure 2. Bus and cars speed.

  a) Example of solution corresponding to case 1. of Definition 3.2 for ρL = ρR = ρ, ρα < ρ < ρα.

  Example of solution corresponding to case 3. of Definition 3.2 for ρL, ρR > ρ * .

Figure 3 .

 3 Figure 3. Different solutions of the Riemann problem (14). Each subfigure illustrates a point of the Definition 3.2: fundamental diagram representation (left) and space-time diagram (right).

(

  16) = (F (ρ) + ẏρ + -żρ + -F α ( ż)) + (F (ρ) + ẏρ --żρ --F α ( ż)) = 2F (ρ) -2F α ( ż) + (ρ + + ρ -)( ẏ -ż) ≤ 2 (F α ( ẏ) -F α ( ż)) ≤ αR| ẏ -ż|. 3.2 (σ -, σ + ) ∈ G 2 ( ż): If sgn(ρ + -σ) = sgn(ρ --σ), we get (16) = sgn(ρ + -σ) [(F (ρ) + ẏρ + -żρ + -F (σ)) -(F (ρ) + ẏρ --żρ --F (σ))]

3 Figure 6 .

 36 Figure 6. Case 3

Acknowledgments. The authors are grateful to Boris Andreianov for suggesting links with the results in [5].

This research was supported by the INRIA associated team Optimal REroute Strategies for Traffic managEment (ORESTE).