On numerical approximation of the Hamilton-Jacobi-transport system arising in high frequency approximations

In the present article, we study the numerical approximation of a system of Hamilton-Jacobi and transport equations arising in geometrical optics. We consider a semi-Lagrangian scheme. We prove the well posedness of the discrete problem and the convergence of the approximated solution toward the viscosity-measure valued solution of the exact problem.


1.
Introduction. In this article, we consider the following system of an Hamilton-Jacobi type equation (HJ) and a continuity equation (CE) describing the transport of the conserved measure m 0 . Even if the vector field a can be smooth (e.g. in the simplest and reference case, a(x, p) = p), the scalar function u as a solution of (1.1) is intended in the viscosity sense. Therefore, ∇u is at most L ∞ , even for smooth initial data u 0 , and the regularity that we can expect for the vector field a(·, ∇u) is the very weak regularity a(· , ∇u(·, ·)) ∈ (L ∞ (R d × [0, T ])) d . (1.4) As a consequence, despite the fact that the (HJ) equation can be solved independently from (CE), system (1.1)-(1.2) leads to some interesting mathematical issues as well as numerical challenges.
In order to obtain global solutions in the general multi-dimensional case, a well adapted notion of solution for (1.2) under hypothesis (1.4) is that of measure solutions introduced by Poupaud-Rascle [26] (see also [25] for the non-conservative transport equation). Indeed, the latter perfectly makes sense if the vector field a(·, ∇u) satisfies a one-sided Lipschitz condition (OSL), (see (2.3) below), and this condition can be obtained by the semiconcavity of the viscosity solution of (1.1), at least in the reference case a(x, p) = p, or in the one dimensional case for a class of functions a.
Existence, uniqueness and stability results for problem (1.1)-(1.2)-(1.3) in the framework of viscosity-measure solutions have been given for example in [2,29]. Therefore, our goal here is not to refine these previous results but to construct consistent numerical approximations.
Finite difference schemes for systems of type (1.1)-(1.2)-(1.3) have been given in [18]. There the authors make use of the notion of duality solution for the (CE) equation, as introduced in [4]. Consequently, they have to restrict the analysis to the one dimensional case. Moreover, in order to obtain discrete semiconcave estimates for the numerical solution of the (HJ) equation, they consider Lax-Friedrichs type schemes for the latter and the strict convexity in p for the hamiltonian H is required. A general class of upwind schemes previously analyzed in [17] has been considered for the approximation of the (CE) equation.
In this paper, the considered scheme is based on a semi-lagrangian discretization of (1.1)-(1.2) for the time approximation coupled with a finite element discretization for the space variable. It requires a convex hamiltonian H and no restriction on the spatial dimension is necessary. Since the semiconcavity of the initial data u 0 is conserved by the scheme, it is sufficient to require only a weak (OSL) condition at the discrete level (automatically satisfied in the reference case a(x, p) = p, see Remark 6) without semiconcavity requirement for the viscosity solution u.
The stability of the Filippov characteristics and of the corresponding measure solution of the transport equation (CE) will be the key tool to prove the convergence of the numerical schemes. As a by-product, we obtain of course a new existence and uniqueness result for the viscosity-measure solution of (1.1)-(1.2)-(1.3).
Regarding the applications of our numerical approximations of (1.1)-(1.2)-(1.3), it is worth to recall that these types of systems arise for example in the semi-classical limit for the Schrödinger equation [8] and of the spinless Bethe-Salpeter equation (the relativistic Schrödinger equation, [3]), or in the high frequency approximation of the Helmholtz equation. In these cases, the hamiltonian H and the vector field a take the forms ; a(x, p) = p , (Schrödinger and Helmholtz), where V (x, t) is the potential (see [2,18,29] for a rapid derivation of (1.5) and (1.6)). It is worth noticing that in the examples above a = c ∇ p H. However, this condition is not required in the rest of the paper. Coupling between first or second order Hamilton-Jacobi equations and transport equations has been also recently considered in the framework of Mean Field Games theory [7,19]. Nonetheless, in this case, the system is given by a backward APPROXIMATION THE HAMILTON-JACOBI-TRANSPORT   3 Hamilton-Jacobi equation and a forward transport equation, with initial condition for the transport equation and terminal condition for the Hamilton-Jacobi equation, plus coupling terms in both equations. Moreover, while MFG preserves the regularity of the initial data, for system (1.1)-(1.2)-(1.3) the use of the measure valued solution cannot be avoided.
The paper is organized as follows. Section 2 is devoted to the preliminary definitions and known results concerning system (1.1)-(1.2)-(1.3). In Section 3 we construct a semi-lagrangian scheme for the time discretization of system (1.1)-(1.2)-(1.3) and we prove its convergence. The corresponding fully discrete scheme is given and analyzed in Section 4. Finally, the appendix is devoted to the proof of some technical lemmas for the reader convenience.
Throughout the paper, we will denote by C 0 0 (R d ) (resp. C 0 c (R d )) the space of continuous functions which tend to 0 at infinity (resp. with compact support); by ρ ε a standard mollifier, i.e. ρ ε ( ; by * the convolution with respect the space variable and by C any numerical constant that can vary from line to line in the computations.

2.
Preliminaries: The viscosity-measure solutions. As mentioned in the introduction, a solution of (1.1)-(1.2) is intended in the viscosity sense for (1.1), while in the sense of Poupaud-Rascle [26] for (1.2). Concerning the definition of viscosity solutions, we refer to the pioneering articles [12,11]. Here, for the reader's convenience, we recall the usual assumptions on H and the consequent general existence and uniqueness result for (1.1) that we shall use in the sequel.
Let us define Q T := R d × [0, T ], T > 0, d ≥ 1, and let the hamiltonian H satisfy the following: As stated in Theorem 2.1, the expected regularity for ∇u is L ∞ (Q T ). Therefore, the characteristics X(t; x) associated to the conservative transport equation (1.2) cannot be defined as classical or distributional solutions of the system below ∂ t X(t; x) = a(X(t; x), ∇u(X(t; x), t)) but have to be understood in a generalized sense. Once the flow X(t; x) is uniquely defined and continuous on [0, T ] × R d , the natural definition of a solution of the conservation law (1.2) with a given initial data m 0 belonging to the space of bounded measures M 1 (R d ) is, see [26], Equation (2.2) means that the measure solution m(t) is the image (or the pushforward) of m 0 by the flow X(t ; ·), i.e. for any Borel set B ⊂ R d , or equivalently m(t), φ = m 0 , φ(X(t; ·)) , for any φ ∈ C 0 0 (R d ) . It turns out that the definition of Filippov characteristics is well suited and that with it, (2.2) perfectly make sense as soon as the characteristics are unique. A definition of Filippov characteristics solution of (2.1) equivalent to the original one and easier to handle is given in [26], as follows: Then, assuming , the vector field a(·, ∇u) satisfies (1.4) and the Filippov characteristics of system (2.1) exist. Moreover, if the following one-sided Lipschitz condition (OSL) holds true (a(x, ∇u(x, t)) − a(y, ∇u(y, t))) · (x − y) ≤ γ(t)|x − y| 2 a.e. t ∈ [0, T ], x, y ∈ R d (2.3) for a function γ ∈ L 1 ([0, T ]), these characteristics are unique, the flow (t, x) → X(t; x) is continuous on [0, T ] × R d , the map from R d to R d , x → X(t; x) is onto; the following existence and uniqueness result of a viscosity-measure solution of (1.1)-(1.2)-(1.3) can be stated, following [2,26]. In the reference case a(x, p) = p, as in (1.5), requiring the (OSL) condition (2.3) is equivalent to require that the viscosity solution u be semiconcave with respect to x, i.e.
for all x, y ∈ R d , t ∈ [0, T ] and for some β ∈ L 1 ([0, T ]). The semiconcavity property (2.4) often characterizes the viscosity solution of Hamilton-Jacobi equations arising in control problems. However, for general vector fields a, it does not imply the necessary condition (2.3). Indeed, for u ε = u * ρ ε , where ρ ε is the mollifier introduced above, the semiconcavity (2.4) implies: ν T D 2 u ε (x, t) ν ≤ β(t), for all ν ∈ R d s.t. |ν| = 1 and for all (x, t) ∈ Q T , (see [6]). Moreover, ∇u ε (x, t) → ∇u(x, t) as ε → 0, for a.a. (x, t) ∈ Q T . Owing to this convergence property, it is sufficient to obtain (2.3) for u ε , whenever a is at least continuous w.r.t. p. Assuming in addition that a is differentiable w.r.t. p and satisfies a one sided Lipschitz condition w.r.t. x locally uniformly in p, we have ). Therefore, the (OSL) condition follows easily from (2.4), if either D p a(x, p) = I, as expected, or d = 1 and ∂ p a(x, p) is non negative and upper bounded, as for (1.6). Apart from these two cases (already considered in [2]), the (OSL) condition is an assumption and has to be verified for the specific model at hand.
We conclude this section by listing the properties that the measure solution (2.2) satisfies, see [26] : It is worth to underline that the uniform compactness at infinity property (v) is fundamental to prove the convergence of approximate measure solutions toward the exact one in 3. The semi-lagrangian scheme. This section is devoted to the construction and the convergence analysis of a semi-lagrangian scheme for the time discretization of system (1.1)-(1.2) over the time interval [0, T ]. For an introduction to this class of schemes we refer to [1,Appendix B] and [10,14,15,16]. To proceed, we need to refine the previous hypothesis on the hamiltonian H and to add assumptions on the growth of H w.r.t. p. Therefore, let us assume in this section (H 1 ), (H 3 ) and (H 2 ) for any R > 0, sup Q T ×B(0,R) |H(x, t, p)| ≡ M (R) < +∞; (H 4 ) H is convex in p and either linear at infinity (i) or superlinear (ii): Remark 1. It is worth noticing that the hamiltonians (1.5) and (1.6) satisfy all the required assumptions as soon as the potential V is uniformly continuous and bounded over Q T and Lipschitz continuous in x uniformly in t. The growth assumption (H 4 )-(i) can be relaxed to include also hamiltonians that are not uniformly positive at infinity. But we leave this generalization to the reader.
3.1. The semi-lagrangian scheme for the (HJ) equation. Let N ∈ N be fixed and let h = T N −1 be the associated time step used in the semi-discrete scheme to be defined. Since the hamiltonian H is convex in p and continuous (in fact lower semi-continuity would be sufficient [27]), then H(x, t, p) = H * * (x, t, p), where H * is the Legendre transform of H with respect to p, i.e. H * (x, t, ξ) = sup p∈R d {ξ · p − H(x, t, p)}, and the (HJ) equation can be written as Next, plugging the first order forward finite difference for the approximation of ∂ t u into (3.1), one easily obtains the following first order approximation for u(x, t + h) For general hamiltonian H whose Legendre transform H * cannot be computed explicitly see [5,9,22,23].
Let us observe that whenever the exact solution u of (HJ) is semiconcave in x, then it possesses one-sided directional derivatives at any x and in any direction, i.e. the limit as h 0 of the right hand side of (3.2) always exists. Moreover, if u is differentiable w.r.t. x, the previous limit coincides with the left hand side of (3.2). In other words, the approximation (3.2) is consistent.
For the semi-discrete scheme, it is sufficient to inductively define the approximation u n = u n (x) of the exact solution u of (HJ) at time t n = n h for n = 0, . . . , N , by with u 0 = u 0 , the initial data for (HJ) in (1.3). Next, we need to show that u n as given by (3.3) shares the properties of u 0 , in the same way as the exact viscosity solution u does. These are well established facts in the context of the approximation of functions by the inf-convolution operator (see for example [6]). However, the difficulties here are on the one hand to show that the properties of the initial data u 0 are propagated to u n uniformly with respect to n and h, and on the other hand to handle the dependency of H on (x, t).
Let us define Q h := R d × {0, . . . , N } and the set of the arguments associated to the infimum in (3.3) The proof of the Lemma above is given in the Appendix. Let us just point out here that we have chosen to take the hamiltonian H to be Lipschitz continuous in x uniformly in p and t when H is linear in p for |p| → +∞, in order to obtain the Lipschitz continuity of H * w.r.t. x. Indeed, in this case the supremum defining H * may not be reached. Therefore we cannot conclude as in case (H 4 )-(ii).  by (3.3), i.e. A n (x) is a non empty bounded set of R d uniformly in (x, n) ∈ Q h and the infimum is a minimum. Moreover, u n ∈ W 1,∞ (R d ) and u n W 1,∞ (R d ) is bounded uniformly in n and, assuming H * ∈ L ∞ (Q T × B(0, K)) when the growth condition (H 4 )-(i) is satisfied, u n is Lipschitz continuous with respect to n uniformly in x, i.e. there exist three positive constants C 0 , C 1 and C 2 independent of h and n, s.t.
where C n 0 and C n 1 are independent of h. Then, the proof will follow by induction on n.
It follows immediately by Lemma 3.1 that u n+1 is upper bounded since Moreover, u n+1 is obviously lower bounded since u n is bounded and H * (x, t, ξ) ≥ −M , for any ( Next, if H satisfies (H 4 )-(i), A n (x) ⊂ B(0, K) for any x ∈ R d and the infimum in (3.3) is attained due to the continuity of u n (x − ξh) + h H * (x, t n , ξ) w.r.t. ξ.
On the other hand, if H satisfies (H 4 )-(ii), since u n and H * (x, t, 0) are bounded and H * is superlinear, there exists R n > 0 (increasing w.r.t. n but upper bounded uniformly in n and h) s.t. A n (x) ⊂ B(0, R n ) for any x ∈ R d and again the infimum in (3.3) is attained.
Let us prove now that u n+1 is Lipschitz continuous. We have for any x, y ∈ R d and for α n (x) ∈ A n (x) and u n+1 (y) ≤ u n (y − h α n (x)) + h H * (y, t n , α n (x)) . Therefore, where L H * is the Lipschitz constant of H * w.r.t. x, and the statement follows exchanging the role of x and y. Consequently, u n+1 ∈ W 1,∞ (R d ).

YVES ACHDOU, FABIO CAMILLI AND LUCILLA CORRIAS
It remains to prove that u n is Lipschitz continuous with respect to n. As before, we have for any x ∈ R d and the corresponding argument α n (x) ∈ A n (x) , and the proof is completed, thanks to the uniform boundedness of A n and H * .
Let us observe that when the hamiltonian H grows linearly at infinity, i.e. when the growth condition (H 4 )-(i) is satisfied, then H * is not necessarily upper bounded. Therefore, the additional hypothesis H * ∈ L ∞ (Q T × B(0, K)) in Lemma 3.2 is necessary. It is easily seen that this additional hypothesis is satisfied by the hamiltonian (1.6) for K = 1 √ 2 (see also Remark 2). Finally, in order to obtain semiconcavity, we need the following semiconcavity hypothesis on H * . (H 5 ) H * = H * (x, t, ξ) is semiconcave in x uniformly in t and ξ, with semiconcavity constant C H * conc . Proof. We proceed as in Lemma 3.2 supposing that u n is semiconcave with semiconcavity constant C n conc . Then, for any x ∈ R d , any corresponding argument α n (x) ∈ A n (x) and any y ∈ R d , we have that ≤ u n (x + y − h α n (x)) − 2 u n (x − h α n (x)) + u n (x − y − h α n (x)) + h [H * (x + y, t n , α n (x)) − 2 H * (x, t n , α n (x)) + H * (x − y, t n , α n (x))] Hence, u n+1 is semiconcave with C n+1 conc ≤ C n conc + h C H * conc , and by induction u n is semiconcave for all n = 1, . . . , N with : C n conc ≤ C u0 conc + T C H * conc . Remark 3. Assumption (H 5 ) is necessary since the hamiltonian H depends on x and t. It can be satisfied under regularity hypothesis on H. In the case of the Hamiltonians (1.5) and (1.6), (H 5 ) is satisfied if the potential V is convex in x uniformly in t.

3.2.
The semi-lagrangian scheme for the (CE) equation. Now we can proceed to the construction of a semi-discrete scheme for the transport equation (1.2), replacing the time continuous flow X(t; x) with a time discrete one, (see [13,24] for similar schemes in the framework of crowd dynamics). Set u n ε = u n * ρ ε and let us consider the following implicit backward Euler scheme X n+1 = X n + h a(X n+1 , ∇u n+1 ε (X n+1 )) , n = 0, · · · , N − 1, Next, given the initial measure m 0 and replacing (2.1) with (3.7), we define the semi-discrete approximation of the measure solution (2.2) as m n = X n (·) # m 0 , for n = 0, . . . , N , i.e. for any Borel set or equivalently m n , φ = m 0 , φ(X n (·)) for any φ ∈ C 0 0 (R d ). Above, the dependence of X n on x and ε and consequently of m n on ε has been skipped to simplify the notations, but it will be made explicit when necessary below.
It is worth noticing that it is not possible to define the semi-discrete flow (3.7) with ∇u n+1 instead of the gradient of the regularized function u n+1 ε , since we need X n to be defined for all (x, n) ∈ Q h . Moreover, it is also not possible to use an explicit forward Euler scheme to define the trajectories X n , because the discrete (OSL) condition (3.9) does not provide the necessary equicontinuity of these trajectories, in contrast with the implicit one (3.7), (see Lemma 3.4 and Remark 4 below).
The following Lemma gives us the existence, uniqueness and the regularity properties of X n necessary to the well posedness for m n and to pass to the limit as h and ε go to 0. Owing to the implicit nature of the Euler scheme (3.7), we have to assume an upper bound on the space-time step h.
Proof. The existence of X n+1 given X n , is insured by the Brouwer fixed point theorem, the map y → X n + h a(y, ∇u n+1 ε (y)) being continuous and bounded. The uniqueness of X n+1 and the Lipschitz continuity w.r.t. x follow both from (3.9) and the upper bound on h. Indeed, for X n+1 and Y n+1 defined by (3.7) we have Taking X n = Y n we get the uniqueness, while iterating over n we get for two starting points x, y ∈ R d Therefore, for δ ∈ (0, 1) s.t. C h ≤ 1 − δ, we have that 1 + Ch/(1 − Ch) ≤ e Ch/(1−Ch) ≤ e Chδ −1 and |X n (x) − X n (y)| ≤ exp(n C h δ −1 )|x − y| ≤ exp(CT δ −1 )|x − y| , n = 0, · · · , N .
Finally, for any x ∈ R d it holds true that where C 1 is given in (3.6), and we obtain the Lipschitz continuity w.r.t. n. The above estimate also gives us that the image by X n of B(0, R) is contained in the ball of radius (R + T sup |a(x, p)|). Therefore, X n is locally bounded uniformly in ε and n and the Lemma is proved. Proof. By Lemma 3.4, the map x → X n (x) is uniquely defined and continuous over R d . Moreover, it is onto from R d to R d . Therefore, the approximated solution m n is well defined in M 1 (R d ). Furthermore, it is easy to see that m n satisfies all the properties (i)-(v) of the exact measure solution m, uniformly in n. Concerning (v), it follows from the local and uniform in n boundedness of X n proved above and the contraction property (iii) (see Lemma 3.1 in [26]).
3.3. The convergence. We can finally discuss the convergence of the semi-discrete scheme The key tools to prove the convergence result have been given in [26], where the authors analyzed the stability of the measure solution with respect to perturbations of the vector field a. The stability of the Filippov characteristics is of course the basic stone. Later, these tools have been also used in [2] to analyze the stability of the viscosity-measure solution of (1.1)-(1.2)-(1.3) with a(x, p) = p, with respect to the vanishing viscosity perturbation of the (HJ) equation.
To discuss the convergence, we define the piecewise constant w.r.t. time approximated solutions , and u ε h = u h * ρ ε . Moreover, we also need to define the following time continuous trajectories by linear interpolation of X n (x) and X n+1 (x), for any x ∈ R d , n = 0, · · · , N Note that (3.11) gives us time continuous trajectories X ε h with the same regularity as X n , i.e. X ε h is locally bounded and Lipschitz continuous w.r.t. (t, x) ∈ [0, T ]×R d , uniformly in ε and h. . Assume in addition that a ∈ (C 0 (R d × R d )) d , u n ε satisfies the (OSL h ) condition (3.9) and Ch < 1. Let ε = ε(h) be s.t. ε(h) → 0 as h → 0. Then, as h → 0, h (x, t) → ∇u(x, t) at any point (x, t) ∈ Q T of differentiability of u, where u is the unique viscosity solution of (1.1); (ii) X ε h converges locally uniformly in Q T toward the unique Filippov characteristic X associated to the vector field a(·, ∇u); , where m is the unique measure solution of (1.2).
Proof. Following standard results in the viscosity solution theory (see for instance [28]), it can be proved that there exists a constant C independent of h, s.t.
Next, let (x, t) ∈ Q T be a point of differentiability of u. For n = t h , by Taylor expansion and the semiconcavity of u n ε , we have that By the uniform bound on ∇u n , there exists p such that, up to a subsequence extraction, ∇u n → p. Therefore, the previous inequality and the convergence of u n ε (t) toward u(t) give us i.e. p belongs to the subdifferential of u at x. Being u differentiable at (x, t), p = ∇u(x, t) and we have obtained the proof of (i). We now prove the second part of the statement. From Lemma 3.4 it immediately follows that, up to a subsequence, X ε h converges uniformly on every compact set of Q T to an absolutely continuous function Y , as h → 0. Next, since by (i) and the (OSL h ) condition (3.9), the (OSL) condition (2.3) is also satisfied, the Filippov characteristic X associated to the vector field a(·, ∇u) is unique for any x ∈ R d . In order to prove that Y = X, let us define ∆ ε h (t; x) := X ε h (t; x) − X(t; x). We are going to prove that ∆ ε h (t; x) → 0 as h → 0. Indeed, we have By the Definition 2.2 of Filippov characteristics it holds true, for all x ∈ R d , a.e. t ∈ [0, T ] and all r > 0, that a(y, ∇u(y, t)) · ∆ ε h (t; x) .
(3.15) In order to estimate I 1 , we decompose ∆ ε h (t; x) as |a(x, p)|, where C 1 is given in (3.6), and we make use of the (OSL h ) condition (3.9) to get (3.16) From (3.16), (3.15) and the uniform boundedness of ∆ ε h (t; x) on every compact subset of Q T , we have obtained Using the convergence results obtained above and passing to the limit into the previous inequality as h → 0 and δ → 0, we get that Finally, from the continuity of Y and X we deduce that Y = X. By the uniqueness of X and the uniform boundedness of X ε h , we have that all the sequence X ε h converge. It remains to prove the convergence of m ε h . This can be obtained exactly as in Proposition 3.1 in [26] and we give the proof for the reader convenience. By the strong convergence of the trajectories X ε h and the Lebesgue dominated convergence theorem, the sequence m ε h (t), φ converges toward m(t), φ as h → 0, for any φ ∈ C 0 0 (R d ). By the uniform compactness at infinity of m ε h (property (v), Lemma 3.5), it is sufficient to consider test functions φ ∈ C 0 c (R d ). Then, (iii) follows by the Ascoli-Arzela theorem since the time equicontinuity of X ε h yields the time equicontinuity of m ε h (t), φ .
Proceeding as in (2.5), one can obtain the (OSL h ) condition (3.9) for special a. This is summarized in the following Corollary which is a straightforward consequence of Theorem 3.6. It is worth noticing that the Hamiltonians (1.5) and (1.6) enter in the framework of this Corollary. Corollary 1. Assume the same hypothesis of Theorem 3.6, except the (OSL h ) condition (3.9). If in addition, a is differentiable w.r.t. p, satisfies a one sided Lipschitz condition w.r.t. x locally uniformly in p and either D p a(x, p) = I or d = 1 and ∂ p a(x, p) is nonnegative and upper bounded, then the same conclusions as in Theorem 3.6 hold true.
Remark 4 (The explicit Euler scheme). The explicit Euler scheme X n+1 = X n + h a(X n , ∇u n ε (X n )) , (3.17) does not provide a family of equicontinuous characteristics under the (OSL h ) condition (3.9). This is a difficulty naturally intrinsic to the (OSL) condition that allows discontinuity of compressive type only, while the map x → X n (x) given by (3.17) can be expansive. Indeed, we have + h 2 |a(X n , ∇u n ε (X n )) − a(Y n , ∇u n ε (Y n ))| 2 , i.e., with the same constant A as in Theorem 3.6, and iterating over n : |X n (x) − X n (y)| 2 ≤ e 2C T |x − y| 2 + 8 A 2 T h .

4.
The fully discrete semi-lagrangian scheme. In this section we introduce a finite element discretization of (3.3) and an approximation of (3.8) by a bounded discrete measure, yielding a fully discrete scheme for (1.1)-(1.2)-(1.3).
For an arbitrarily fixed space step k > 0, we consider the regular uniform grid of R d given by X k := {x i = ik , i ∈ Z d } . Let T k = {S k j } j∈J k be an associated collection of non-degenerate, pairwise disjoint and uniform simplices having as vertices lattice points x i ∈ X k and covering R d , (S k j are triangles in dimension 2 and tetrahedra in dimension 3). We denote also by w is linear on S k j , j ∈ J k }, the space of continuous and piecewise linear functions on T k . Then, each w ∈ W k can be expressed as for basis functions β k i ∈ W k satisfying β k i (x j ) = δ ij for i, j ∈ Z d . It immediately follows that any β k i has compact support, 0 ≤ β k i (x) ≤ 1, i∈Z d β k i (x) = 1 and at any x ∈ R d at most (d+1) functions β k i are non-zero. In the sequel, the interpolation operator defined by (4.1) will be denoted P k .
The fully discrete approximation of (3.3) based on the above space discretization is naturally given by    u n k (x) = i∈Z d β k i (x) u n k,i , n = 0, · · · , N , where u 0 k,i = u 0 (x i ). It is obvious that, thanks to the continuous piecewise linear interpolation of the discrete approximation (u n k,i ) i∈Z d on the lattice at any time step, the continuous function u n k ∈ W k shares the properties of the semi-discrete approximation u n given in Lemma 3.2. Therefore, the equivalent of Lemma 3.2 for the fully discrete approximation u n k will be skipped here. On the other hand, the semiconcavity of u n k , given a semiconcave initial data u 0 , is not straightforward. Therefore, we shall give the equivalent of Lemma 3.3 here and the proof in the Appendix. and (H 5 ), u n k is discretely semiconcave for all n = 0, · · · N , i.e. it verifies both u n k (x+x j )−2 u n k (x)+u n k (x−x j ) ≤ (C u0 conc +T C H * conc )|x j | 2 , x ∈ R d , j ∈ Z d , (4.3) Furthermore, property (4.13) is also satisfied by m n k uniformly in n = 0, . . . , N , by the definition of the coefficients m n k,i and (4.12). With the above definitions, set U n k = (u n k,i ) i∈Z d and M n k = (m n k,i ) i∈Z d . The fully discrete scheme (4.2), (4.16), reads U n+1 k = inf ξ∈R d {B k (ξ) U n k + h L n (ξ)} , n = 0, · · · , N − 1 , M n k = Λ n k M 0 k , n = 1, · · · , N , (4.17) with B k (ξ) = (β k i (x j − hξ)) i,j∈Z d , L n (ξ) = (H * (x i , t n , ξ)) i∈Z d Λ n k = (β k i (X n k (x j ))) i,j∈Z d .
It is worth noticing that B k is a stochastic matrix, while Λ n k is the transpose of a stochastic matrix.
Defining again the following time piecewise constant approximations (4.18) as well as u ε h,k = u h,k * ρ ε and the time linear interpolated trajectories X ε h,k exactly as in (3.11), we can prove the convergence of the fully discrete scheme (4.17). . Assume in addition that the (OSL k h ) condition (4.8) is satisfied with h sufficiently small such that C h < 1 and that (4.11) is satisfied with ε = λh, (λ > C ). Then, as h → 0 and k h 2 → 0, (i) u h,k → u in L ∞ (Q T ) and ∇u ε h,k (x, t) → ∇u(x, t) in any point (x, t) ∈ Q T of differentiability of u, where u is the unique viscosity solution of (1.1); (ii) X ε h,k converges locally uniformly in Q T toward the unique Filippov characteristic X associated to the vector field a(·, ∇u); Proof. We shall prove that there exists a constant C independent of h and k, such that u n k − u n L ∞ (R d ) ≤ C(n + 1)k , n = 0, . . . , N . As a consequence, by the time regularity of u and by estimate (3.12), it follows that Estimate (4.19) is obvious for n = 0, being u 0 k = P k u 0 and u 0 = u 0 . Let us suppose that it holds true for a given n. Then, for any argument α n (x i ) ∈ A n (x i ), we have by (3.3), (3.4) and (4.2) u n+1 k,i − u n+1 (x i ) ≤ u n k (x i − h α n (x i )) − u n (x i − h α n (x i )) ≤ C(n + 1)k .
Exchanging the roles of u n+1 k,i and u n+1 (x i ) we obtain so that for any x ∈ R d , Since the convergence of ∇u ε h,k can be proved exactly as in Theorem 3.6, from the weak-semiconcavity of u n k,ε , statement (i) is proved. Next, from (4.19), we have the same estimate for the regularized approximation u n k,ε and u n ε , yielding ∇u n k,ε − ∇u n ε L ∞ (R d ) ≤ C k ε h = C λ k h 2 , n = 0, . . . , N . (4.20)