Olivier Bokanowski

AND Yingda Cheng
email: ycheng@math.msu.edu

Chi-Wang Shu
email: shu@dam.brown.edu

CONVERGENCE OF DISCONTINUOUS GALERKIN SCHEMES FOR FRONT PROPAGATION WITH OBSTACLES

Keywords: February 21, 2015 Hamilton-Jacobi-Bellman equations, discontinuous Galerkin methods, level sets, front propagation, obstacle problems, dynamic programming principle, convergence

We study semi-Lagrangian discontinuous Galerkin (SLDG) and Runge-Kutta discontinuous Galerkin (RKDG) schemes for some front propagation problems in the presence of an obstacle term, modeled by a nonlinear Hamilton-Jacobi equation of the form min(u t +cu x , u-g(x)) = 0, in one space dimension. New convergence results and error bounds are obtained for Lipschitz regular data. These "low regularity" assumptions are the natural ones for the solutions of the studied equations. Numerical tests are given to illustrate the behavior of our schemes.

Introduction

In this paper, we establish convergence of a class of discontinuous Galerkin (DG) methods for the one-dimensional Hamilton-Jacobi (HJ) equation below, hereafter also called the "obstacle" equation, min(u t + c u x , ug(x)) = 0, x ∈ I = (0, 1), t > 0, (1a) u(0, x) = u 0 (x), x ∈ (0, 1), (1b) with periodic boundary conditions on I and a constant c ∈ R. In [START_REF] Abgrall | Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes[END_REF], the function g is called the "obstacle" function.

It is well known that taking constraints in optimal control problems is not an obvious task. Within the viscosity theory, it is possible to devise schemes for obstacle equations such as [START_REF] Abgrall | Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes[END_REF]. However a monotonicity condition is needed in general for proving the convergence of the scheme. The monotonicity condition can yield a convergence proof of one-half order in the mesh size [START_REF] Crandall | Two approximations of solutions of Hamilton Jacobi equations[END_REF] (see also [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF] for more specific partial differential equations (PDEs) with an obstacle term and related error estimates). However, a serious restriction of such monotonicity condition is that the schemes become at most first order accurate for smooth solutions or in smooth regions, thus making the schemes highly inefficient for practical computation. On the other hand, it is very difficult to show convergence of formally higher order schemes, such as the ones studied in this paper, when the solution is not regular enough.

In our previous work [START_REF] Bokanowski | A discontinuous galerkin scheme for front propagation with obstacles[END_REF], we proposed a class of Runge-Kutta DG (RKDG) methods adapted to front propagation problems with obstacles. The DG methods under consideration were originally devised to solve conservation laws, see for example the review paper [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convectiondominated problems[END_REF]. As for the DG-HJ solvers, in [START_REF] Hu | A discontinuous Galerkin finite element method for Hamilton-Jacobi equations[END_REF][START_REF] Li | Reinterpretation and simplified implementation of a discontinuous Galerkin method for Hamilton-Jacobi equations[END_REF], the first efforts relied on solving the conservation law system satisfied by the derivative of the solution. See also [START_REF] Chen | An adaptive high-order discontinuous Galerkin method with error control for the Hamilton-Jacobi equations. part I: The one-dimensional steady state case[END_REF] for an adaptive version of this scheme. In [START_REF] Cheng | A discontinuous Galerkin finite element method for directly solving the Hamilton-Jacobi equations[END_REF], a DG method for directly solving the Hamilton-Jacobi equation was developed and was later generalized to solve front propagation problems [START_REF] Bokanowski | A discontinuous Galerkin solver for front propagation[END_REF] and obstacle problems [START_REF] Bokanowski | A discontinuous galerkin scheme for front propagation with obstacles[END_REF]. Other direct DG solvers include the central DG scheme [START_REF] Li | A central discontinuous Galerkin method for Hamilton-Jacobi equations[END_REF] and the local DG scheme [START_REF] Yan | A local discontinuous Galerkin method for directly solving Hamilton-Jacobi equations[END_REF]. The schemes proposed in [START_REF] Bokanowski | A discontinuous galerkin scheme for front propagation with obstacles[END_REF] feature a simple treatment of the obstacle functions. Stability analysis is performed with forward Euler, a Heun scheme and a TVD Runge-Kutta third order (TVD-RK3) time discretization using the techniques developed in Zhang and Shu [START_REF] Zhang | Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF].

On the other hand, the semi-Lagrangian DG (SLDG) methods were proposed in [START_REF] Restelli | A semi-Lagrangian discontinuous Galerkin method for scalar advection by incompressible flows[END_REF][START_REF] Rossmanith | A positivity-preserving high-order semi-Lagrangian discontinuous Galerkin scheme for the Vlasov-Poisson equations[END_REF][START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: theoretical analysis and application to the Vlasov-Poisson system[END_REF] to compute incompressible flow and Vlasov equations, as well as in [START_REF] Bokanowski | Semi-Lagrangian discontinuous Galerkin schemes for some first and second order partial differential equations[END_REF] for some general linear first and second order PDEs. The advantage of SLDG is its ability to take large time steps without a CFL restriction. However, it is difficult to design SLDG methods for nonlinear problems (some SL schemes for Hamilton-Jacobi-Bellman equations were proposed in [START_REF] Bokanowski | An adaptive sparse grid semi-Lagrangian scheme for first order Hamilton-Jacobi Bellman equations[END_REF], but without convergence proof). For general SL methods, we refer to the works of Falcone and Ferretti [START_REF] Falcone | Discrete time high-order schemes for viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Falcone | Convergence analysis for a class of high-order semi-Lagrangian advection schemes[END_REF][START_REF] Falcone | Semi-Lagrangian schemes for Hamilton-Jacobi equations, discrete representation formulae and Godunov methods[END_REF], see also the textbook [START_REF] Falcone | Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations[END_REF]. There is also a vast literature on high order finite difference schemes for solving HJ equations, see, e.g. [START_REF] Osher | High order essentially non-oscillatory schemes for Hamilton-Jacobi equations[END_REF][START_REF] Abgrall | Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes[END_REF][START_REF] Jiang | Weighted ENO schemes for Hamilton-Jacobi equations[END_REF][START_REF] Kurganov | New high-resolution semi-discrete central schemes for Hamilton-Jacobi equations[END_REF].

Beyond the scope of the present paper, yet of great interest, we also mention the second order PDE with obstacle terms, such as min(u tcu xx , ug(x)) = 0, with c > 0. This is the case of the so-called "American options" in mathematical finance [START_REF] Achdou | Computational Methods for Option Pricing[END_REF]. Explicit schemes were proposed and proved to converge within the viscosity theory (see for instance [START_REF] Jakobsen | On the rate of convergence of approximation schemes for Bellman equations associated with optimal stopping time problems[END_REF]) yet with a reduced rate of convergence. Variational methods for nonlinear obstacle equations can also be devised [START_REF] Hintermüller | An inverse problem in American options as a mathematical program with equilibrium constraints: C-stationarity and an active-set-Newton solver[END_REF] but will lead in general to nonlinear implicit schemes which are computationally more demanding.

The scope of the present paper is to study convergence of the SLDG and RKDG schemes for the obstacle problem [START_REF] Abgrall | Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes[END_REF]. The main challenges include the low regularity of the solution and the nonlinear treatment needed to obtain the obstacle solution. Due to the fully discrete nature of the method, traditional techniques for obtaining semi-discrete error estimates of DG methods for hyperbolic problems cannot directly apply here. Therefore, fully discrete analysis is necessary. Fully discrete analysis of RKDG methods for conservation laws has been performed in the literature. In [START_REF] Zhang | Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF], error estimates for RKDG methods for scalar conservation laws were provided for smooth solutions. In [START_REF] Cockburn | Error estimates for the Runge-Kutta discontinuous Galerkin method for the transport equation with discontinuous initial data[END_REF][START_REF] Zhang | Error estimates for the third order explicit Runge-Kutta discontinuous Galerkin method for linear hyperbolic equation in one-dimension with discontinuous initial data[END_REF], discontinuous solutions with RK2 and linear polynomials and RK3 time discretizations with general polynomials were studied for linear conservation laws. However, to our best knowledge, convergence results for second and higher order DG schemes solving nonlinear hyperbolic equations with irregular solutions are not available.

As a consequence of our results, assuming that h is a space step and ∆t a time step, we shall show error bounds of the order of O(h 9/10 ) for the SLDG schemes (under time stepping ∆t ≡ Ch 3/5 , larger time step can be taken with a lower convergence rate), and of order O(h 1/2 ) for the RKDG schemes (under time stepping ∆t ≡ Ch). We will need a natural "no shattering" regularity assumption on the exact solution that will be made precise in the sequel, otherwise we will typically assume the exact solution to be Lipschitz regular and piecewise C q regular for some q ≥ 1 for SLDG (resp. q ≥ 2 for RKDG).

The main idea of our proof is based on the dynamic programming principles as illustrated below. The viscosity solution of (1) also corresponds to the following optimal control problem:

u(t, x) = max u 0 (x -ct), max θ∈[0,t] g(x -cθ) . (2)
The function u is also the solution of the Bellman's dynamic programming principle (DPP): for any ∆t > 0,

u(t + ∆t, x) = max u(t, x -c∆t), max θ∈[0,∆t] g(x -cθ) , ∀ t ≥ 0. (3)
Notice conversely that the DPP (3), together with u(0, x) = u 0 (x), implies [START_REF] Achdou | Computational Methods for Option Pricing[END_REF].

If we denote u n (x) := u(t n , x) and

g t (x) := max θ∈[0,t] g(x -cθ),
then the DPP implies in particular for any x and n ≥ 0:

u n+1 (x) = max(u n (x -c∆t), g ∆t (x)). ( 4
)
Using formula (2), we can see that when u 0 and g are Lipschitz regular, then u is also Lipschitz regular in space and time, and in general no more regularity can be assumed (the maximum of two regular functions is in general no more than Lipschitz regular). When u 0 is a discontinuous function (otherwise piecewise regular), formula (2) implies also some regularity on the solution u.

The rest of the paper is organized as follows. In Section 2, we describe the DG methods for the obstacle equations. In Section 3, we collect some lemmas that will be used in our convergence proofs. Section 4 and Section 5 are devoted to the convergence analysis of SLDG and RKDG methods, respectively. In both sections, we will first establish error estimates for SLDG and RKDG schemes for linear transport equations without obstacles. Then we will use DPP illustrated above to prove convergence of the numerical solution in the presence of obstacles. Numerical examples are given in Section 6. We conclude with a few remarks for the multi-dimensional case in Section 7.

DG schemes for the obstacle equation

In this section, we will introduce the SLDG and RKDG methods for the obstacle equation [START_REF] Abgrall | Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes[END_REF]. For simplicity of discussion, in the rest of the paper, we will assume c to be a positive constant.

Let

I j := (x j-1 2 , x j+ 1 2
), j = 1, . . . , N be a set of intervals forming a partition of I = (0, 1). We denote h := max j h j where h j = |I j | is the length of the interval I j . For a given integer k ≥ 0, let V h be the DG space of piecewise polynomial of degree at most k on each interval I j :

V h := {v : I → R, v |I j ∈ P k , ∀j}. (5)
To introduce the methods for the obstacle problems, we follow two steps. Firstly, we describe the DG solvers for the linear advection equation v t + cv x = 0.

To advance the numerical solution in one step from v n h ∈ V h to v n+1 h ∈ V h , we consider one of the two DG methods described below.

SLDG Scheme:

v n+1 h := Π h v n h (• -c∆t) , (6)
where Π h is the L 2 projection onto the space V h . We shall denote this SLDG solver by

v n+1 h = G SL ∆t (v n h ).
RKDG Scheme (by TVD-RK3 time stepping):

find v n,1 h , v n,2 h , v n+1 h ∈ V h , such that (v n,1 h -v n h , ϕ h ) = ∆tH(v n h , ϕ h ), ∀ϕ h ∈ V h (7a) (v n,2 h - 3 4 v n h - 1 4 v n,1 h , ϕ h ) = ∆t 4 H(v n,1 h , ϕ h ), ∀ϕ h ∈ V h (7b) (v n+1 h - 1 3 v n h - 2 3 v n,2 h , ϕ h ) = 2∆t 3 H(v n,2 h , ϕ h ), ∀ϕ h ∈ V h (7c) where (φ, ϕ) = I φ ϕ dx and H j (φ h , ϕ h ) = I j cφ h (ϕ h ) x dx -c((φ h ) - j+ 1 2 (ϕ h ) - j+ 1 2 -(φ h ) - j-1 2 (ϕ h ) + j-1 2 ), H(φ h , ϕ h ) = j H j (φ h , ϕ h ).
We shall denote this RKDG solver by

v n+1 h = G RK ∆t (v n h ).
After introducing DG schemes for the linear transport problem, for the obstacle equation (1), we shall consider two approaches: one by using the L 2 projection

u n+1 h := Π h max(G ∆t (u n h ), g) , (8)
where G ∆t = G SL ∆t or G RK ∆t and g ≡ g ∆t or g ≡ max(g(x), g(x -c∆t)).

The idea of the formulation above is to try to follow the relation (4), and the projection step is to project the function into the piecewise polynomial space V h .

Unfortunately, the scheme ( 8) is difficult to implement, because we need to compute the maximum of two functions, which requires locating the roots of G ∆t (u n h )-g. Another more practical approach is to define u n+1 h as the unique polynomial in V h such that

u n+1 h (x j α ) := max G ∆t (u n h )(x j α ), g(x j α ) , ∀ j = 1, . . . , N, α = 0, . . . , k, (10)
where (x j α ) α=0,...,k are the k + 1 Gauss-Legendre quadrature points on the interval I j , and w j α are the corresponding quadrature weights. Those schemes were studied in details and stability was established in [START_REF] Bokanowski | A discontinuous galerkin scheme for front propagation with obstacles[END_REF]. We shall see in later sections that definitions [START_REF] Chen | An adaptive high-order discontinuous Galerkin method with error control for the Hamilton-Jacobi equations. part I: The one-dimensional steady state case[END_REF] or [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] lead to similar error estimates, although the second approach is much easier to implement.

Finally, we remark that in [START_REF] Bokanowski | A discontinuous galerkin scheme for front propagation with obstacles[END_REF], forward Euler and TVD-RK2 temporal discretizations are also considered. However, the stability restriction for the time step for the forward Euler method is rather severe as ∆t ≤ Ch 2 , and the stability proof of TVD-RK2 only works for piecewise linear polynomials. Therefore, in this paper, we will only consider TVD-RK3 time discretizations.

Preliminaries

In this section, we collect some lemmas which will be used in our convergence proof, and discuss properties of the obstacle solutions. Here and below, we use C (possibly with subscripts) to denote a positive constant depending solely on the exact solution, which may have a different value in each occurrence.

Let us introduce, for ℓ ≥ 0, the following function sets:

C ℓ+1 p,L,c 0 (0, 1) := v : (0, 1) → R, v Lipschitz continuous with v ′ L ∞ ≤ L, v piecewise C ℓ+1 with v (ℓ+1) L ∞ ≤ c 0 ,
and v admits at most p ≥ 0 non regular points.

where v (ℓ+1) denotes the (ℓ + 1)-th derivative almost everywhere, and ∆ 2 q (0, 1) := g : (0, 1) → R, g has at most q local maxima points and g is twice differentiable at each local maxima.

The following ℓ 2 pseudo-norm definition will also be used:

f ℓ 2 := i,α w i α |f (x i α )| 2 h i 1/2 . ( 11
)
In particular, using the Gauss-Legendre quadrature rule, for any f ∈ V h we have

f ℓ 2 = f L 2 .
From this point on, we will use f to denote f L 2 , and f D to denote f L 2 (D) for a given domain D.

3.1. Properties of projections and the obstacle function. For the RKDG method, it is necessary to consider the following Legendre-Gauss-Radau projection P h . For any function ϕ, P h ϕ ∈ V h , and for any element I j , it holds that

(P h ϕ) - j+1/2 = ϕ - j+1/2 , I j (P h ϕ -ϕ) ψ h dx = 0, ∀ ψ h ∈ P k-1 (I j ).
In the lemma below, we will first establish the projection properties for functions in the space C ℓ+1 p,L,c 0 .

Lemma 3.1. Let ℓ ≥ 0 and let ϕ be in C ℓ+1 p,L,c 0 : ϕ is Lipschitz continuous, piecewise C ℓ+1 , with at most p ≥ 0 non regular points. Then there exists a constant C ≥ 0, depending only on ℓ, p, L, c 0 such that

ϕ -P h ϕ ≤ Ch q k,ℓ ,
where the projection P h can be either P h or Π h , and q k,ℓ is defined as

q k,ℓ := min min(k, ℓ) + 1, 3 2 ≡ 3/2 if k, ℓ ≥ 1 1 if k = 0 or ℓ = 0 (12)
Proof. Using the property of the projections [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF], on each regular cell I j , we have

ϕ -P h ϕ I j ≤ Ch min(ℓ,k)+1 ϕ H ℓ+1 (I j ) ≤ Ch min(ℓ,k)+1
since ϕ ∈ C l+1 on I j . Therefore,

ϕ -P h ϕ L 2 (∪I j , ϕ |I j regular) = ∪I j , ϕ |I j regular |ϕ -P h ϕ| 2 dx 1/2 ≤ Ch min(ℓ,k)+1 .
On the other hand, on the intervals I j where ϕ is not regular, we have

ϕ -P h ϕ I j ≤ Ch min(0,k)+1 ϕ H 1 (I j ) ≤ Ch ϕ H 1 (I j ) ≤ Ch 3/2
because ϕ is Lipschitz continuous. Hence summing up on the "bad" intervals I j (with at most p such intervals),

ϕ -P h ϕ L 2 (∪ bad I) ≤ Cp 1/2 h 3/2 .
Summing up the bounds with bad intervals and good ones, we prove the desired result.

We now state a similar estimate for the ℓ 2 norm: Lemma 3.2. Assume that ϕ belongs to C ℓ+1 p,L,c 0 , ℓ ≥ 0, then we have

ϕ -P h ϕ ℓ 2 = j,α w j α v(x j α ) -(P h ϕ)(x j α ) 2 h j 1/2 ≤ Ch q k,ℓ (13)
where P h = P h or Π h , q k,ℓ is defined as in [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convectiondominated problems[END_REF], and the constant C depends only on p, L and c 0 .

Proof. On each regular cell I j , we have [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF],

ϕ -P h ϕ L ∞ (I j ) ≤ Ch min(ℓ,k)+1/2 ϕ H ℓ+1 (I j ) ≤ Ch min(ℓ,k)+1 .
Then using the fact that w j α ≥ 0, α w j α = 1 we obtain that

   I j , ϕ |I j regular α w j α ϕ n (x j α ) -(P h ϕ n )(x j α ) 2 h j    1/2 ≤ Ch min(ℓ,k)+1 .
On the other hand, when the interval I j is such that ϕ contains a non-regular point, we can write

ϕ -P h ϕ L ∞ (I j ) ≤ Ch min(0,k)+1/2 ϕ H 1 (I j ) ≤ Ch. Therefore    I j , ϕ |I j not regular α w j α ϕ n (x j α ) -(P h ϕ n )(x j α ) 2 h j    1/2 ≤    I j , ϕ |I j not regular h(Ch) 2    1/2 ≤ Cp 1/2 h 3/2 .
This concludes our proof.

We now turn to some estimates related to the obstacle function g ∆t .

Lemma 3.3. Assume that g ∈ ∆ 2 q (0, 1) for some integer q ≥ 1. Let g(x) := max(g(x), g(x -c∆t)). There exists a constant C ≥ 0 such that

g -g ∆t ≤ C √ q∆t 5/2 . ( 14
)
and

g -g ∆t ℓ 2 ≤ C √ q √ ∆t + h ∆t 2 . (15) Proof. Denote M g the set of local maximum points of g, if [x -c∆t, x] ∩ M g = ∅, then we see that g(x) -g ∆t (x) = 0. Furthermore {x, g ∆t (x) =g(x)} dx ≤ {x, [x-c∆t,x]∩Mg =∅} dx ≤ q c∆t (16)
since there are at most q local maxima.

Consider the case when [x-c∆t, x] contains at least a local maxima of g. Assuming that ∆t is small enough we can assume that x * is the only local maximum of g on the interval [x -c∆t, x], so that g ∆t (x) = g(x * ). Then,

|g(x) -g ∆t (x)| = |g(x) -g(x * )| ≤ C|x -x * | 2 ≤ C∆t 2 , ( 17
) since g is twice differentiable at x * .
Combining ( 17) and ( 16) we obtain the bound [START_REF] Falcone | Discrete time high-order schemes for viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF].

For the second estimate [START_REF] Falcone | Convergence analysis for a class of high-order semi-Lagrangian advection schemes[END_REF], we make use of a minimal covering ∪I of the set

x, [x -c∆t, x] ∩ M g = ∅ ,
using mesh intervals. The length of this covering is bounded by c∆t + 2h for each maximum point, since in order to cover any interval [a, b] we may need two more mesh intervals I, of length ≤ 2h than the minimum required length ba. Overall the length of the total covering is bounded by q(c∆t + 2h), hence we obtain the desired result.

3.2.

Properties of the obstacle solutions. We shall impose some restrictions on the regularity of the obstacle solutions as described below.

Definition 3.1. ["no shattering" property] For a given T ≥ 0, we will say that the exact solution u of the problem (1) is "not shattering" if there exists some ℓ ≥ 0 and constants p, L, c 0 such that the exact solution satisfies, for any t

∈ [0, T ], u(t, .) ∈ C ℓ+1 p,L,c 0 . Recall that the exact solution satisfies u(t, x) = max(u 0 (x -ct), g t (x)) = u 0 (x - ct) + max(0, g t (x) -u 0 (x -ct))
. Therefore if u is not shattering, it implies that u 0 (xct)g t (x) has bounded number of zeros (since otherwise x → u(t, x) would have an unbounded number of singularities).

A typical example where shattering occurs (therefore not satisfying definition 3.1), can be constructed as follows. Suppose the domain Ω contains the interval (-2, 2), let u 0 (x) ≡ (x -1) + (x -1) 3 sin(1/(x -1)) and g(x) ≡ x on [-2, 2], together with velocity constant c = 1. The function u 0 is of class C 2 on the interval (-2, 2). Notice then, g t (x) = max θ∈[0,t] g(xθ) = g(x) for t ∈ [0, 1] and x ∈ [-1, 1], and the exact solution is therefore u

(t, x) = max(u 0 (x-t), g(x)) for t ∈ [0, 1] and x ∈ [-1, 1]. So at time t = 1, u(1, x) = max(u 0 (x-1), g(x)) ≡ x+max(x 3 sin(1/x), 0) (for x ∈ [-1, 1]).
This function has an infinite number of non regular points in the interval [-1, 1], and therefore does not satisfy the "no shattering" property at time t = 1.

It is not easy to state precise conditions on the initial data u 0 and g to ensure that the no shattering property will be satisfied. Mainly, ∀t, x → g t (x)u 0 (xct) should have only a finitely bounded number of zeros, as is detailed below. However, it is clear that shattering will not occur for generic data u 0 and g. This definition still allows for a finite (bounded) number of singularities in u(t n , .), as is generally the case when taking the maximum of two regular functions. Finally we also give an example (see Lemma 3.5 below) where we can prove that the "no shattering" condition is satisfied. Lemma 3.4. Assume that, for a given T > 0,

g ∈ C ℓ+1
pg,Lg,cg (0, 1), (18a)

u 0 ∈ C ℓ+1 pu 0 ,Lu 0 ,cu 0 (0, 1), (18b) and that, ∀ t ∈ [0, T ],
x → g t (x)u 0 (xct) has a finitely bounded number of zeros in (0, 1), (18c) the bound being independent of t.

Then there exists constants p, L, c 0 such that, ∀t ∈ [0, T ], u(t, x) belongs to C ℓ+1 p,L,c 0 (with L = max(L g , L u 0 ), c 0 = max(c g , c u 0 ), and p that are independent of t.)

Proof. First, for g ∈ C ℓ+1 pg,Lg,cg (0, 1), if we denote M (resp. m) to be the number of local maxima (resp. minima) of g, then we have g t ∈ C ℓ+1 pg+2M +m,Lg,cg (0, 1). This is because the Lipschitz constant of g t is bounded by L g by an elementary verification. Each local maxima of g may develop into two singularities in g t , each local minima may develop into one singularity in g t (•), and each singular point of g may continue to be a singularity in g t (•). Hence the total number of non-regular points in g t (•) will be bounded by 2M + m + p g . The bound of the (ℓ + 1) derivative can also be obtained easily.

Because the exact solution is given by u

(t, x) = max(u 0 (x -ct), g t (x)) = u 0 (x - ct) + max(0, g t (x) -u 0 (x -ct)), the Lipschitz constant of u(t, .) is therefore bounded by max(L u 0 (•-ct) , L g ) = max(L u 0 , L g ).
On the other hand, the number of singular points of max(0, g t (x)u 0 (xct)) is bounded by the sum of the number of singular points of the function x → g t (x)u 0 (xct)), plus the number of zeros of the same function, which is assumed to be bounded independently of t ≥ 0. Hence the the number of singular points of x → u(t, x) is bounded independently of t ≥ 0.

Finally the bound on the x-partial derivative u (ℓ+1) (t, •) L ∞ , in the regular region of u, is easily obtained, because u is then locally one of the two functions g t or u 0 (•ct). Lemma 3.5. Assume that g and u 0 satisfy the regularity assumptions (18a) and (18b) of Lemma 3.4. Assume furthermore that there exists a constant L 1 ≥ 0 such that |u ′ 0 (x)| ≥ L 1 > L g for a.e. x ∈ (0, 1). Then (18c) holds and u satisfies the "no shattering" assumption.

Proof. As we can see in Lemma 3.4, g t also has a Lipschitz constant ≤ L g . On the other hand on each regular part of u 0 there is a slope ≥ L 1 which is strictly greater that L g . By an elementary verification, one can show that assumption (18c) is satisfied and that the result of Lemma 3.4 can be applied.

Convergence of the SLDG schemes

In this section, we will provide the convergence proof of the SLDG scheme. In particular, we will proceed in three steps. First, we will establish error estimates of the SLDG methods for the linear transport equation

v t + c v x = 0, v(0, x) = v 0 (x).
We will then generalize the results to scheme [START_REF] Chen | An adaptive high-order discontinuous Galerkin method with error control for the Hamilton-Jacobi equations. part I: The one-dimensional steady state case[END_REF], and finally to scheme [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] for the obstacle problem.

4.1. Convergence of the SLDG scheme for the linear advection equation.

We first consider the linear equation v t + cv x = 0, for which

v(t + ∆t) = v(t, x -c∆t).
We denote v n (•) = v(t n , •), and we define the numerical solution of the SLDG method at t n to be v n h . In particular, the scheme writes: initialize with v 0 h := Π h v 0 , and

v n+1 h = G SL ∆t (v n h ) = Π h (v n h (• -c∆t)) for n ≥ 0. Theorem 4.1. We consider v t + cv x = 0, v(0, x) = v 0 (x). If v 0 ∈ C ℓ+1 p,L,c 0 (0, 1), then we have for all n such that n∆t ≤ T , v n h -v n ≤ CT h q k,ℓ ∆t
for some constant C ≥ 0 independent of n, and q k,ℓ is defined in [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convectiondominated problems[END_REF].

Proof. If v 0 ∈ C ℓ+1 p,L,c 0 (0, 1), then v n ∈ C ℓ+1 p,L,c 0 (0, 1), and due to Lemma 3.1, we have:

v n (• -c∆t) -Π h (v n (• -c∆t)) ≤ Ch q k,ℓ ,
for some constant C independent of n. Hence

v n+1 h -v n+1 = Π h (v n h (• -c∆t) -v n (• -c∆t) ≤ Π h (v n h (• -c∆t) -Π h (v n (• -c∆t) + Π h (v n (• -c∆t) -v n (• -c∆t) ≤ Π h (v n h (• -c∆t) -Π h (v n (• -c∆t) + Ch q k,ℓ ≤ v n h (• -c∆t) -v n (• -c∆t) + Ch q k,ℓ ≤ v n h -v n + Ch q k,ℓ
, where we have used the fact Π h u ≤ u in the fourth row, and the periodic boundary conditions in the last row. As for the initial condition, we have

v 0 h -v 0 = Π h v 0 -v 0 ≤ Ch q k,ℓ .
Finally we obtain for any given T ≥ 0 the existence of a constant C ≥ 0 (independent of T and of n), such that

v n h -v n ≤ C(n + 1) h q k,ℓ ≤ CT h q k,ℓ ∆t , for all n such that n∆t ≤ T . Therefore, if min(k, ℓ) = 0, then v n h -v n ≤ CT h ∆t , and we need h = o(∆t) for the convergence. If otherwise min(k, ℓ) ≥ 1, it holds v n h -v n ≤ CT h 3/2
∆t and we would only need h = o(∆t 2/3 ) for the convergence.

4.2.

Convergence of the first SLDG scheme in the obstacle case. Now we turn to scheme [START_REF] Chen | An adaptive high-order discontinuous Galerkin method with error control for the Hamilton-Jacobi equations. part I: The one-dimensional steady state case[END_REF] for the nonlinear equation [START_REF] Abgrall | Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes[END_REF]. In particular, the scheme writes: initialize with u 0 h := Π h u 0 , and u n+1 h = Π h (max(G SL ∆t (u n h ), g)) for n ≥ 0. Theorem 4.2. Assume that the exact solution u is not shattering in the sense of Definition 3.1, for some integer ℓ ≥ 0. Then the following error bound holds:

u n h -u n ≤ CT h q k,ℓ ∆t + CT g -g ∆t ∆t for some constant C ≥ 0 independent of n. In particular, (i) If g = g ∆t , then ∀ t n ≤ T : u n h -u n ≤ CT h q k,ℓ ∆t (ii) If g(x) := max g(x)
, g(x -c∆t) , and if g ∈ ∆ 2 q (0, 1), then ∀ t n ≤ T :

u n h -u n ≤ CT h q k,ℓ ∆t + CT ∆t 3/2 .
Proof. By Lemma 3.1 and the "no shattering" assumption,

Π h u n -u n ≤ Ch q k,ℓ . (19)
Hence, from the definitions,

u n+1 h -u n+1 ≤ u n+1 h -Π h u n+1 + Π h u n+1 -u n+1 ≤ Π h (max(G SL ∆t (u n h ), g)) -Π h (max(u n (• -c∆t), g ∆t )) + Ch q k,ℓ ≤ max(G SL ∆t (u n h ), g ∆t ) -max(u n (• -c∆t), g) + Ch q k,ℓ . Using the fact that | max(a 1 , b 1 ) -max(a 2 , b 2 )| ≤ |a 1 -a 2 | + |b 1 -b 2 |, we obtain u n+1 h -u n+1 ≤ G SL ∆t (u n h ) -u n (• -c∆t) + g -g ∆t + Ch q k,ℓ ≤ Π h (u n h (• -c∆t)) -Π h (u n (• -c∆t)) + g -g ∆t + Ch q k,ℓ ≤ u n h -u n + g -g ∆t + Ch q k,ℓ (20)
where we have used again [START_REF] Hu | A discontinuous Galerkin finite element method for Hamilton-Jacobi equations[END_REF] and the periodic boundary condition. Using Lemma 3.3 and by induction on n, we are done. ∆t )+O(∆t 3/2 ). Therefore the optimal estimate is obtained when h 3/2 ≡ ∆t 5/2 , or ∆t ≡ h 3/5 , and the error is of order O(h 9/10 ). 4.3. Convergence of the SLDG scheme defined with Gauss-Legendre quadrature points. Now we turn to scheme [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] for the nonlinear equation [START_REF] Abgrall | Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes[END_REF]. In particular, the scheme writes: initialize with u 0 h := Π h u 0 , and u n+1 h is defined as the unique polynomial in V h such that:

u n+1 h (x j α ) := max(G SL ∆t (u n h )(x j α ), g(x j α )), ∀j, α. (21)
Theorem 4.3. Let ℓ ≥ 0 and assume that the exact solution is not shattering in the sense of Definition 3.1. The following error bound holds:

u n h -u n ≤ CT h q k,ℓ ∆t + CT g -g ∆t ∆t ,
for some constant C ≥ 0 independent of n. In particular, (i) If g = g ∆t , then ∀t n ≤ T :

u n h -u n ≤ CT h q k,ℓ ∆t .
(ii) If g(x) := max g(x), g(x -c∆t) and g ∈ ∆ 2 q (0, 1), then ∀t n ≤ T :

u n h -u n ≤ CT h q k,ℓ ∆t + CT ∆t √ ∆t + h.
Proof. In view of Lemma 3.1 and the "no shattering" property, we have

u n+1 -Π h u n+1 ≤ Ch q k,ℓ .
Now we turn back to the error estimate, and consider the case of g = g ∆t :

u n+1 h -u n+1 ≤ u n+1 h -Π h u n+1 + Ch q k,ℓ = u n+1 h -Π h u n+1 ℓ 2 + Ch q k,ℓ ≤ u n+1 h -u n+1 ℓ 2 + Ch q k,ℓ ≤ j,α w j α u n+1 h (x j α ) -u n+1 (x j a ) 2 h j 1/2 + Ch q k,ℓ
where in the third line we have used ( 13) for u n+1 . Because of the DPP, for all points x, the exact solution satisfies:

u n+1 (x) = max(u n (x -c∆t), g ∆t (x)).
Therefore, for ∀j, α:

(u n+1 h -u n+1 )(x j α ) ≤ G SL ∆t (u n h )(x j α ) -u n (x j α -c∆t) + g(x j α ) -g ∆t (x j α )
and

u n+1 h -u n+1 ≤ j,α w j α G SL ∆t (u n h )(x j α ) -u n (x j α -c∆t) 2 h j 1/2 + j,α w j α g(x j α ) -g ∆t (x j α ) 2 h j 1/2 + Ch q k,ℓ ≤ G SL ∆t (u n h ) -u n (• -c∆t) ℓ 2 + g -g ∆t ℓ 2 + Ch q k,ℓ ≤ G SL ∆t (u n h ) -Π h u n (• -c∆t) ℓ 2 + g -g ∆t ℓ 2 + Ch q k,ℓ = Π h u n h (• -c∆t) -Π h u n (• -c∆t) ℓ 2 + g -g ∆t ℓ 2 + Ch q k,ℓ = Π h u n h (• -c∆t) -Π h u n (• -c∆t) + g -g ∆t ℓ 2 + Ch q k,ℓ ≤ u n h (• -c∆t) -u n (• -c∆t) + g -g ∆t l 2 + Ch q k,ℓ ≤ u n h -u n + g -g ∆t ℓ 2 + Ch q k,ℓ
where in the fourth line we have used Lemma 3.2 for the function u n (• -c∆t). Using Lemma 3.3 and by induction on n, we are done.

Remark 4.2. Defining g as in (ii), and assuming min(ℓ, k) ≥ 1, the error is bounded by O( h 3/2 ∆t )+O(∆t √ ∆t + h). Therefore the optimal estimate is obtained when h 3/2 ≡ ∆t 2 √ ∆t + h. So h ∆t → 0, h 3/2 ≡ ∆t 5/2 , or ∆t ≡ h 3/5 (as in Remark 4.1), and the error of the scheme defined with Gauss-Legendre quadrature points is again of order O(h 9/10 ) for this particular time stepping.

Convergence of RKDG schemes

In this section, we will prove convergence for the RKDG schemes. We will proceed in three steps similar to the previous section.

Firstly, let us recall the following properties for the bilinear operator H. Lemma 5.1. [START_REF] Zhang | Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF] For any φ h , ϕ h ∈ V h , we have

H(φ h , ϕ h ) + H(ϕ h , φ h ) = - j c[φ h ] j+1/2 • [ϕ h ] j+1/2 H(φ h , φ h ) = - 1 2 j c[φ h ] 2 j+1/2 .
We also recall inverse inequalities [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] for the finite element space V h . In particular, there exists a constant C (independent of h), such that, for any

ϕ h ∈ V h , (ϕ h ) x ≤ Ch -1 ϕ h , ϕ h L ∞ ≤ Ch -1/2 ϕ h .
5.1. Convergence of the RKDG scheme for the linear advection equation.

We first consider the linear equation v t + cv x = 0, for which

v(t + ∆t) = v(t, x -c∆t).
We still denote v n (•) = v(t n , •). In particular, the scheme writes: initialize with v 0 h := Π h v 0 , and

v n+1 h = G RK ∆t (v n h )
for n ≥ 0. This convergence proof closely follows the work in [START_REF] Zhang | Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF] for smooth solution, but additional difficulties are encountered because we consider solutions with less regularity. The main technique is to introduce piecewisely defined intermediate stage functions and the careful treatment of intervals containing irregular points.

Theorem 5.1. We consider v t + cv x = 0, v(0, x) = v 0 (x). Let v 0 be in C ℓ+1 p,L,c 0 (0, 1), ℓ ≥ 2, k ≥ 1, and assume the CFL condition ∆t ≤ C 0 h for C 0 small enough (the usual CFL condition for stability of the RKDG scheme). The following bound holds:

v n h -v n ≤ C 1 (h + h 3 ∆t 2 ) 1/2 , for some constant C 1 ≥ 0 independent of h, ∆t, v h .
In particular, if ∆t/h is bounded from below ( ∆t h ≥ C0 for some constant C0 > 0), then

v n h -v n ≤ C 1 h 1/2 . Proof.
We need to introduce some intermediate stages of the exact solution. Firstly we define v (1)

:= v n -c∆t (v n ) x
Now we are ready to define the errors e (1)

:= v (1) -v n,1 h , ξ (1) := P h v (1) -v n,1 h , η (1)
:= P h v (1)v (1) , e (2)

:= v (2) -v n,2 h , ξ (2) := P h v (2) -v n,2 h , η (2) := P h v (2) -v (2) , (25)
e

n := v n -v n h , ξ n := P h v n -v n h , η n := P h v n -v n .
Clearly, e (1) = ξ (1)η (1) , e (2) = ξ (2)η (2) , e n = ξ nη n .

Our next step is to establish the error equations. First, let us recall that the numerical solution satisfies:

I j v n,1 h ϕ h dx = I j v n h ϕ h dx + ∆tH j (v n h , ϕ h ), ∀ϕ h ∈ V h I j v n,2 h ϕ h dx = 3 4 I j v n h ϕ h dx + 1 4 I j v n,1 h ϕ h dx + ∆t 4 H j (v n,1 h , ϕ h ), ∀ϕ h ∈ V h I j v n+1 h ϕ h dx = 1 3 I j v n h ϕ h dx + 2 3 I j v n,2 h ϕ h dx + 2∆t 3 H j (v n,2 h , ϕ h ). ∀ϕ h ∈ V h .
From the definitions of v (1) , v (2) , v (3) , we can verify

I j v (1) ϕ h dx = I j v n ϕ h dx + ∆tH j (v n , ϕ h ), ∀ϕ h ∈ V h I j v (2) ϕ h dx = 3 4 I j v n ϕ h dx + 1 4 I j v (1) ϕ h dx + ∆t 4 H j (v (⋆1) , ϕ h ), ∀ϕ h ∈ V h + 0, j / ∈ I n 1 4 (L a , ϕ h ), j ∈ I n I j v (3) ϕ h dx = 1 3 I j v n ϕ h dx + 2 3 I j v (2) ϕ h dx + 2∆t 3 H j (v (⋆2) , ϕ h ), ∀ϕ h ∈ V h
where for j ∈ I n , (⋆1) = n, (⋆2) = n; otherwise, (⋆1) = (1), (⋆2) = (2). Notice that the formulations above are correct because we have enforced continuity of the first function appearing in operator H j in all cases. In particular, the procedure to enforce continuity of v (2) at ∂B n turns out to be necessary here. Combining the previous two relations, we derive the error equations

I j
e (1) 1) ϕ h dx + ∆t 4 H j (e (1) , ϕ h ), 2) , ϕ h ), j ∈ I n where Υ = v n+1v (3) . Using the decomposition of errors ( 25), we get

ϕ h dx = I j e n ϕ h dx + ∆tH j (e n , ϕ h ), ∀ϕ h ∈ V h I j e (2) ϕ h dx = 3 4 I j e n ϕ h dx + 1 4 I j e (
∀ϕ h ∈ V h + 0, j / ∈ I n ∆t 4 H j (v n -v (1) , ϕ h ) + 1 4 (L a , ϕ h ), j ∈ I n I j e n+1 ϕ h dx = I j Υϕ h dx + 1 3 I j e n ϕ h dx + 2 3 I j e (2) ϕ h dx, ∀ϕ h ∈ V h + 2∆t 3 H j (e (2) , ϕ h ) + 0, j / ∈ I n 2∆t 3 H j (v n -v (
I j ξ (1) ϕ h dx = I j ξ n ϕ h dx + ∆tJ j (ϕ h ), ∀ϕ h ∈ V h (26a) I j ξ (2) ϕ h dx = 3 4 I j ξ n ϕ h dx + 1 4 I j ξ (1) ϕ h dx + ∆t 4 K j (ϕ h ), ∀ϕ h ∈ V h (26b) I j ξ n+1 ϕ h dx = 1 3 I j ξ n ϕ h dx + 2 3 I j ξ (2) ϕ h dx + 2∆t 3 L j (ϕ h ), ∀ϕ h ∈ V h (26c)
where

J j (ϕ h ) = I j 1 ∆t (η (1) -η n )ϕ h dx + H j (e n , ϕ h ), ∀ϕ h ∈ V h (27a) K j (ϕ h ) = I j 1 ∆t (4η (2) -3η n -η (1) )ϕ h dx + H j (e (1) , ϕ h ), ∀ϕ h ∈ V h (27b) + 0, j / ∈ I n H j (v n -v (1) , ϕ h ) + 1 ∆t (L a , ϕ h ), j ∈ I n L j (ϕ h ) = I j 1 2∆t (3η n+1 -η n -2η (2) + 3Υ)ϕ h dx + H j (e (2) , ϕ h ), ∀ϕ h ∈ V h (27c) + 0, j / ∈ I n H j (v n -v (2) , ϕ h ), j ∈ I n We further denote J (ϕ h ) = j J j (ϕ h ), K(ϕ h ) = j K j (ϕ h ), L(ϕ h ) = j L j (ϕ h ).
By letting ϕ h = ξ n , 4ξ (1) , 6ξ (2) in (27a), (27b), (27c), respectively, we get the following energy equation for ξ n , [START_REF] Zhang | Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF] 3 ξ n+1 2 -3 ξ n 2 = ∆t[J (ξ n ) + K(ξ (1) ) + L(ξ (2) ] (28) + 2ξ (2) ξ (1) ξ n 2 + 3(ξ n+1ξ n , ξ n+1 -2ξ (2) + ξ n ) Now we define Π 1 := ∆t[J (ξ n ) + K(ξ (1) ) + L(ξ (2) ], Π 2 := 2ξ (2) ξ (1) ξ n 2 + 3(ξ n+1ξ n , ξ n+1 -2ξ (2) + ξ n ). We will estimate those two terms separately.

Estimate of Π 1

Firstly, we notice that

∆tJ (ξ n ) = (η (1) -η n , ξ n ) + ∆tH(e n , ξ n ) = (η (1) -η n , ξ n ) + ∆tH(ξ n , ξ n ) = (η (1) -η n , ξ n ) - ∆t 2 j c[ξ n ] 2 j+1/2 ≤ η (1) -η n • ξ n - ∆t 2 j c[ξ n ] 2 j+1/2 ≤ 1 4∆tǫ η (1) -η n 2 + ǫ∆t ξ n 2 - ∆t 2 j c[ξ n ] 2 j+1/2
where in the second line we have used the property of the Legendre-Gauss-Radau projection to get H(η n , ϕ h ) = 0. In the formulas above, ǫ is a positive constant of order 1. Since

η (1) -η n = P h (v (1) -v n ) -(v (1) -v n ) = -∆t c (P h (v n ) x -(v n ) x ) ,
similar to Lemma 3.1, we get

P h (v n ) x -(v n ) x ≤ P h (v n ) x -(v n ) x B n + P h (v n ) x -(v n ) x I\B n ≤ Ch 1/2 + Ch min(ℓ,k+1) ≤ Ch 1/2 , Therefore η (1) -η n ≤ C∆th 1/2 and ∆tJ (ξ n ) ≤ C∆th + ǫ∆t ξ n 2 - ∆t 2 j c[ξ n ] 2 j+1/2 .
Similarly, ∆tK(ξ (1) ) = (4η (2) η (1) -3η n + L a , ξ (1) ) + ∆tH(e (1) , ξ (1) ) + ∆t 1) , ξ (1) ) 1) , ξ (1) )

j∈I n H j (v n -v (1) , ξ (1) ) = (4η (2) -η (1) -3η n + L a , ξ (1) ) - ∆t 2 j c[ξ (1) ] 2 j+1/2 + ∆t j∈I n H j (v n -v (
≤ 1 4∆tǫ 4η (2) -η (1) -3η n 2 + 1 4∆tǫ L a 2 + ǫ 2 ∆t ξ (1) 2 - ∆t 2 j c[ξ (1) ] 2 j+1/2 +∆t j∈I n H j (v n -v (
Since ||L a || ∞ ≤ C∆t 2 , and L a = 0 only on B n , therefore ||L a || ≤ C∆t 2 h 1/2 . Next, we will estimate the term 4η (2) η (1) -3η n and ∆t j∈I H j (v nv (1) , ξ (1) ). We can derive that 4η (2) η (1) -3η n = -∆tc(P h v

(1)

xv

(1)

x ), x / ∈ B n -∆tc(P h v n x -v n x ) + (P h L a -L a ), x ∈ B n
Because L a is a linear polynomial and k ≥ 1, P h L a -L a = 0, and similar to the previous argument, we get 4η (2) η (1) -

3η n ≤ C∆th 1/2 .
As for ∆t j∈I H j (v nv (1) , ξ (1) ), we have

v n -v (1) = c∆t(v n ) x ,
therefore for any j

|(v n -v (1) ) j±1/2 | ≤ C∆t v W 1,∞ and v n -v (1) I j ≤ C∆th 1/2 v W 1,∞ .
Hence

H j (v n -v (1) , ξ (1) ) = I j c(v n -v (1) )ξ (1) x dx -c(v n -v (1) ) -(ξ (1) ) - j+1/2 + c(v n -v (1) ) -(ξ (1) ) + j-1/2 ≤ C∆th 1/2 ξ (1)
x

I j + C∆t|(ξ (1) ) - j+1/2 | + C∆t|(ξ (1) ) + j-1/2 | ≤ C∆th -1/2 ξ (1)
I j by inverse inequalities, and ∆t j∈I n H j (v nv (1) , ξ (1)

) ≤ C∆t 2 h -1/2 ξ (1) B n ≤ C∆t 2 + ǫ 2 ∆t ξ (1) 2 B n ≤ C∆t 2 + ǫ 2 ∆t ξ (1) 2 ,
where in the second line we have used the CFL condition ∆t ≤ C cf l h. Putting everything together, and using the CFL condition again, we have ∆tK(ξ (1) )

≤ C∆th + ǫ∆t ξ (1) 2 - ∆t 2 j c[ξ (1) ] 2 j+1/2 Finally, ∆tL(ξ (2) ) = 1 2 (3η n+1 -2η (2) -η n + 3Υ, ξ (2)
) + ∆tH(e (2) , ξ (2) ) + ∆t

j∈I n H j (v n -v (2) , ξ (2) ) = 1 2 (3η n+1 -2η (2) -η n + 3Υ, ξ (2) ) - ∆t 2 j c[ξ (2) ] 2 j+1/2 + ∆t j∈I n H j (v n -v (2) , ξ (2)
) 2) , ξ (2) ).

≤ 1 8∆tǫ 3η n+1 -2η (2) -η n 2 + 9 8∆tǫ Υ 2 + ǫ∆t ξ (2) 2 - ∆t 2 j c[ξ (2) ] 2 j+1/2 +∆t j∈I n H j (v n -v (
Using the same argument as the previous terms, we have 3η n+1 -2η (2)η n ≤ C∆th 1/2 and ∆t j∈I n H j (v nv (2) , ξ (2) ) ≤ C∆t 2 + ǫ∆t ξ (2) 2 . As for Υ, we have

Υ 2 = B n Υ 2 dx + I\B n Υ 2 dx = B n (v n+1 -v (3) ) 2 dx + I\B n Υ 2 dx ≤ B n (v n+1 -v n + c∆t(v n ) x - 1 6 L a ) 2 dx + (C∆t 4 ) 2 = B n (v n+1 -v n + c∆t(v n ) x ) 2 dx + C∆t 4 h + (C∆t 4 ) 2 ≤ C∆t 2 h + C∆t 4 h + (C∆t 4 ) 2 ≤ C∆t 2 h.
Finally we obtain ∆tL(ξ (2)

) ≤ C∆th + ǫ∆t ξ (2) 2 - ∆t 2 j c[ξ (2) ] 2 j+1/2 .
Putting everything together, we have

Π 1 ≤ C∆th + ǫ∆t ξ n 2 + ǫ∆t ξ (1) 2 + ǫ∆t ξ (2) 2 - ∆t 2 j c[ξ n ] 2 j+1/2 - ∆t 2 j c[ξ (1) ] 2 j+1/2 - ∆t 2 j c[ξ (2) ] 2 j+1/2
Estimate of Π 2 To estimate Π 2 , we first introduce

G 1 := ξ (1) -ξ n G 2 := 2ξ (2) -ξ (1) -ξ n G 3 := ξ n+1 -2ξ (2) + ξ n .
From the error equation ( 26), we can deduce

I j G 1 ϕ h dx =∆tJ j (ϕ h ), ∀ϕ h ∈ V h (29a) I j G 2 ϕ h dx = ∆t 2 (K j (ϕ h ) -J j (ϕ h )), ∀ϕ h ∈ V h (29b) I j G 3 ϕ h dx = ∆t 3 (2L j (ϕ h ) -K j (ϕ h ) -J j (ϕ h )), ∀ϕ h ∈ V h (29c) Now, Π 2 = (G 2 , G 2 ) + 3(G 1 , G 3 ) + 3(G 2 , G 3 ) + 3(G 3 , G 3 ). First, let us estimate (G 2 , G 2 ) + 3(G 1 , G 3 ). (G 2 , G 2 ) + 3(G 1 , G 3 ) = -G 2 2 + 2(G 2 , G 2 ) + 3(G 1 , G 3 ) = -G 2 2 + ∆t [K(G 2 ) -J (G 2 ) + 2L(G 1 ) -K(G 1 ) -J (G 1 )] We have ∆t(K(G 2 ) -J (G 2 )) = (4η (2) -3η n -η (1) -(η (1) -η n ) + L a , G 2 ) + ∆tH(e (1) -e n , G 2 ) +∆t j∈I n H j (v n -v (1) , G 2 ) = (4η (2) -3η n -η (1) -(η (1) -η n ) + L a , G 2 ) + ∆tH(G 1 , G 2 ) +∆t j∈I n H j (v n -v (1) , G 2 ) ≤ 4η (2) -3η n -η (1) -(η (1) -η n ) + L a 2 + 1 4 G 2 2 +∆tH(G 1 , G 2 ) + ∆t j∈I n H j (v n -v (1) , G 2 ) ≤ C∆t 2 h + 1 4 G 2 2 + ∆tH(G 1 , G 2 ) + ∆t j∈I n H j (v n -v (1) , G 2 ) ≤ C∆t 2 h + 1 4 G 2 2 + ∆tH(G 1 , G 2 ) + C∆t 2 h -1/2 G 2 ≤ C∆t 2 h + 1 4 G 2 2 + ∆tH(G 1 , G 2 ) + C∆t 4 h -1 + 1 4 G 2 2 ≤ C∆t 2 h + 1 2 G 2 2 + ∆tH(G 1 , G 2 ).
For the other term, similarly we get

∆t(2L(G 1 ) -K(G 1 ) -J (G 1 )) = (3η n+1 -2η (2) -η n + 3Υ -(4η (2) -3η n -η (1) ) -(η (1) -η n ) -L a , G 1 )
+∆tH(2e (2) e (1) e n , G 1 ) + ∆t

j∈I n H j (2(v n -v (2) ) -(v n -v (1) ), G 1 ) ≤ C∆th + ∆tH(G 2 , G 1 ) + ǫ 2 ∆t G 1 2 ≤ C∆th + ∆tH(G 2 , G 1 ) + ǫ∆t ξ n 2 + ǫ∆t ξ (1) 2 .
Therefore,

(G 2 , G 2 ) + 3(G 1 , G 3 ) ≤ - 1 2 G 2 2 + C∆th + ∆tH(G 2 , G 1 ) + ∆tH(G 1 , G 2 ) + ǫ∆t ξ n 2 + ǫ∆t ξ (1) 2 ≤ - 1 2 G 2 2 + C∆th -∆t j c[G 1 ][G 2 ] + ǫ∆t ξ n 2 + ǫ∆t ξ (1) 2 ≤ - 1 2 G 2 2 + C∆th + ∆t 4 j c|[G 1 ]| 2 + ∆t j c|[G 2 ]| 2 + ǫ∆t ξ n 2 + ǫ∆t ξ (1) 2 . Also 3(G 3 , G 2 ) = ∆t(2L(G 2 ) -K(G 2 ) -J (G 2 )) = (3η n+1 -2η (2) -η n + 3Υ -(4η (2) -3η n -η (1) ) -(η (1) -η n ) -L a , G 2 )
+∆tH(2e (2) e (1) e n , G 2 ) + ∆t

j∈I n H j (2(v n -v (2) ) -(v n -v (1) ), G 2 ) ≤ C∆th + ∆tH(G 2 , G 2 ) + ǫ∆t G 2 2 ≤ C∆th - ∆t 2 j c[G 2 ] 2 j+1/2 + ǫ∆t G 2 2 . and 3 G 3 2 = 3(G 3 , G 3 ) = ∆t(2L(G 3 ) -K(G 3 ) -J (G 3 )) = (3η n+1 -2η (2) -η n + 3Υ -(4η (2) -3η n -η (1) ) -(η (1) -η n ) -L a , G 3 ) +∆tH(2e (2) -e (1) -e n , G 3 ) + ∆t j∈I n H j (2(v n -v (2) ) -(v n -v (1) ), G 3 ) ≤ C∆th 1/2 G 3 + C ∆t h G 2 • G 3 + C∆t 2 h -1/2 G 3 Therefore G 3 ≤ C∆th 1/2 + C G 2
due to the CFL condition and

3 G 3 2 ≤ C∆t 2 h + 1 4 G 2 2 .
Putting everything together, we have

Π 2 ≤ (- 1 4 + ǫ∆t) G 2 2 + C∆th + ∆t 4 j c[G 1 ] 2 j+1/2 + ∆t 2 j c[G 2 ] 2 j+1/2 + ǫ∆t ξ n 2 + ǫ∆t ξ (1) 2 ≤ (- 1 4 + ǫ∆t) G 2 2 + C∆th + ∆t 2 j c([ξ n ] 2 j+1/2 + [ξ (1) ] 2 j+1/2 ) +C ∆t h G 2 2 + ǫ∆t ξ n 2 + ǫ∆t ξ (1) 2 ≤ (- 1 4 + ǫ∆t + CC cf l ) G 2 2 + C∆th + ∆t 2 j c([ξ n ] 2 j+1/2 + [ξ (1) ] 2 j+1/2 ) +ǫ∆t ξ n 2 + ǫ∆t ξ (1) 2
When ǫ and C cf l are small enough, -1 4 + ǫ∆t + CC cf l ≤ 0, and

Π 2 ≤ C∆th + ∆t 2 j c([ξ n ] 2 j+1/2 + [ξ (1) ] 2 j+1/2 ) + ǫ∆t ξ n 2 + ǫ∆t ξ (1) 2 .
Finally

Π 1 + Π 2 ≤ C∆th + 2ǫ∆t ξ n 2 + 2ǫ∆t ξ (1) 2 + ǫ∆t ξ (2) 2 .
At this point, we need to provide an estimate of ξ (1) , ξ (2) to finish the proof. Plug in the error equation [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: theoretical analysis and application to the Vlasov-Poisson system[END_REF], ξ (1) 2 = (ξ (1) , ξ (1) ) = (ξ n , ξ (1) ) + ∆tJ (ξ (1) ) ≤ ξ n • ξ (1) + (η (1) η n , ξ (1) ) + ∆tH(ξ nη n , ξ (1) )

≤ ξ n • ξ (1) + C∆th 1/2 ξ (1) + C ∆t h ( ξ n + η n ) ξ (1) ≤ C ξ n • ξ (1) + C∆th 1/2 ξ (1)
Therefore,

ξ (1) ≤ C ξ n + C∆th 1/2 .
Similarly, i.e.

ξ (2) ≤ C ξ n + C ξ (1) + C∆th 1/2 ≤ C ξ n + C∆th 1/2 . Overall, Π 1 + Π 2 ≤ C∆th + C∆t ξ n 2 ,
ξ n+1 2 ≤ (1 + C∆t) ξ n 2 + C∆th (30)
and by induction with the initial condition satisfying ||ξ 0 || ≤ Ch q k,ℓ ,

ξ n ≤ Ch 1/2
and we are done using the projection property η n ≤ Ch q k,ℓ , since q k,ℓ = 3 2 in this case.

The final bound is

u n h -u n ≤ C (h + h 3 ∆t 2 ) 1/2 .
In the case h/∆t is bounded from below, we obtain a bound of order h 1/2 . 5.2. Convergence of RKDG scheme in the obstacle case. Now we turn to scheme [START_REF] Chen | An adaptive high-order discontinuous Galerkin method with error control for the Hamilton-Jacobi equations. part I: The one-dimensional steady state case[END_REF] for the obstacle equation [START_REF] Abgrall | Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes[END_REF]. In particular, the scheme writes: initialize with u 0 h := Π h u 0 , and

u n+1 h = Π h (max(G RK ∆t (u n h ), g))
for n ≥ 0. The main idea follows closely the proof of Theorem 4.2, but utilizes the estimates in Theorem 5.1.

Theorem 5.2. Assume that the exact solution is not shattering in the sense of Definition 3.1. Under the same assumption as in Theorem 5.1 (in particular the assumption on the CFL condition) with both g = g ∆t and g = max g(•), g(• -c∆t) , scheme [START_REF] Chen | An adaptive high-order discontinuous Galerkin method with error control for the Hamilton-Jacobi equations. part I: The one-dimensional steady state case[END_REF] with RKDG solver G RK ∆t satisfies

u n h -u n ≤ C (h + h 3 ∆t 2 ) 1/2 for some constant C ≥ 0 independent of h, ∆t, u h .
In particular, if ∆t/h is bounded from below ( ∆t h ≥ C0 for some constant C0 > 0), we have

v n h -v n ≤ C 1 h 1/2 . Proof.
Similar to the proof of Theorem 4.2, we can obtain

e n+1 = u n+1 h -u n+1 ≤ G RK ∆t (u n h ) -u n (• -c∆t) + g -g ∆t + Ch q k,ℓ . We decompose the error u n h -u n = η n -ξ n , where ξ n = P h u n -u n h , η n = P h u n -u n , and G RK ∆t (u n h ) -u n (• -c∆t) = ξ ′ -η ′ where ξ ′ = P h u n (• -c∆t) -G RK ∆t (u n h ), η ′ = P h u n (• -c∆t) -u n (• -c∆t).
Therefore

e n+1 ≤ ξ ′ + η ′ + g -g ∆t + Ch q k,ℓ ≤ ξ ′ + g -g ∆t + Ch q k,ℓ (31)
by using the projection property again. Hence

e n+1 2 ≤ (1 + ǫ) ξ ′ 2 + (1 + 1 ǫ ) ( g -g ∆t + Ch q k,ℓ ) 2 (32)
Using [START_REF] Zhang | Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF] in the proof of Theorem 5.1, we obtain

ξ ′ 2 ≤ (1 + C∆t) ξ n 2 + C∆th Hence e n+1 2 ≤ (1 + ǫ) (1 + C∆t) ξ n 2 + C∆th + (1 + 1 ǫ ) ( g -g ∆t + Ch q k,ℓ ) 2
Now we take ǫ = C∆t,

e n+1 2 ≤ (1 + C∆t) ξ n 2 + C∆th + C∆t -1 ( g -g ∆t + Ch q k,ℓ ) 2 ≤ (1 + C∆t) ξ n 2 + C∆th + C∆t -1 h 3 ≤ (1 + ǫ)(1 + C∆t) e n 2 + (1 + 1 ǫ )(1 + C∆t) η n 2 + C∆th + C∆t -1 h 3 ≤ (1 + C∆t) e n 2 + C∆th + C∆t -1 h 3
where in the last line we have taken ǫ = C∆t again. By induction on n, we are done.

5.3. Convergence of the RKDG scheme defined with Gauss-Legendre quadrature points. Now we turn to scheme [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] for the nonlinear equation [START_REF] Abgrall | Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes[END_REF]. In particular, the scheme writes: initialize with u 0 h := Π h u 0 , and u n+1 h is defined as the unique polynomial in V h such that :

u n+1 h (x j α ) := max(G RK ∆t (u n h )(x j α ), g(x j α )), ∀j, α. (33)
Theorem 5.3. Assume that the exact solution is not shattering in the sense of Definition 3.1. Under the same assumption as in Theorem 5.1 (in particular assuming the CFL condition) with both g = g ∆t and g(x) := max g(x), g(x -c∆t) , scheme [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] with RKDG solver G RK ∆t satisfies

u n h -u n ≤ C (h + h 3 ∆t 2 ) 1/2 for some constant C ≥ 0 independent of h, ∆t, u h .
In particular, if ∆t/h is bounded from below ( ∆t h ≥ C0 for some constant C0 > 0), we have

v n h -v n ≤ C 1 h 1/2 .
Proof. Following the same lines as the proof for Theorem 4.3, we obtain

e n+1 = u n+1 h -u n+1 ≤ G RK ∆t (u n h ) -Π h u n (• -c∆t) ℓ 2 + g -g ∆t ℓ 2 + Ch q k,ℓ = G RK ∆t (u n h ) -Π h u n (• -c∆t) + g -g ∆t ℓ 2 + Ch q k,ℓ ≤ G RK ∆t (u n h ) -u n (• -c∆t) + g -g ∆t ℓ 2 + Ch q k,ℓ
. By the same argument as in Theorem 5.2

e n+1 2 ≤ (1 + C∆t) ξ n 2 + C∆th + C∆t -1 ( g -g ∆t ℓ 2 + Ch q k,ℓ ) 2 ≤ (1 + C∆t) ξ n 2 + C∆th + C∆t -1 h 3 ≤ (1 + C∆t) e n 2 + C∆th + C∆t -1 h 3 ,
and we are done.

Numerical results

In this section we consider one-and a two-dimensional examples to validate our results. In the one dimensional setting the SLDG scheme and the RKDG scheme are tested and in the two dimensional setting only the RKDG scheme is tested. Example 1 (1-d). This is a one-dimensional test (same as [4, Example 1]):

min(u t + u x , u -g(x)) = 0, t > 0, x ∈ [-1, 1], (34) u(0, x) = u 0 (x), x ∈ [-1, 1], (35)
with periodic boundary conditions and g(x) := sin(πx), u 0 (x) := 0.5 + sin(πx). In that case, for times 0 ≤ t ≤ 1, the exact solution is given by :

u (1) (t, x) =    max(u 0 (x -t), g(x)) if t < 1 3 , max(u 0 (x -t), g(x), 1 x∈[0.5,1] ) if t ∈ [ 1 3 , 1 3 + 1 2 ], max(u 0 (x -t), g(x), 1 x∈[-1,t-1 3 -1 2 ]∪[0.5,1] ) if t ∈ [ 1 3 + 1 2 , 1]
. We first consider the RKDG scheme [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. The choice ∆t = 0.2 h is made considering the stability of the RKDG scheme with P 2 elements.

In Figure 1 the numerical solution is shown, which agrees well with the exact solution everywhere. Table 1 contains the numerical errors at time t = 0.5. The errors are computed globally (a uniform grid mesh of 10 4 points is used in each mesh cell to estimate the errors). In this example there are three singular points (s 1 ≃ -0.1349, s 2 := 0.5 and s 3 = 2/3) where the solution is continuous but with different left and right derivatives. We use a least square procedure to calculate the approximate order of the scheme, see Figure 2. From the calculation, the orders for L 1 , L 2 , L ∞ errors are 1.75, 1.35 and 0.97, respectively.

Next, in Table 2 and in Table 3, numerical errors for the SLDG scheme are given at time t = 0.5. Table 2 shows the results for the choice of time step ∆t = h/2. Then, in Table 3 the choice of ∆t of the order of h 3/5 is made, as suggested in the theoretical study (number of time steps is 10(N x /10) 3/5 ) and since there is no CFL restriction on this scheme. The L 1 , L 2 , L ∞ order of the methods based on the least square plots Figure 3 are 1.66, 1.21, 0.72 (when ∆t = h/2) and 1.48, 1.37, 0.95 (when ∆t = C h 3/5 ).

We can clearly observe from this example that the numerical orders of convergence in the L ∞ norm are at least as good as and often better than those obtained by our analysis.

-1 0 1 -2 -1 0 1 2 t=0 Exact Obstacle DG -1 0 1 -2 -1 0 1 2 t=0.5 Exact Obstacle DG -1 0 1 -2 -1 0 1 2 t=1 Exact Obstacle DG Figure 1.
Example 1, RKDG scheme, times t = 0, t = 0.5 and t = 1, using P 2 elements with N x = 20 mesh cells (obstacle : green dotted line).

(u t + 1 2 u x + 1 2 u y , u -g(x, y)) = 0, t > 0, (x, y) ∈ Ω, (36)
u(0, x, y) = u 0 (x, y), (x, y) ∈ Ω, (37) where g(x, y) := sin(π(x+y)), u 0 (x, y) = 0.5+g(x, y), and Ω = [-1, 1] 2 with periodic boundary conditions. The exact solution is known and is obtained as in Example 1: u(t, x, y) = u (1) (t, x + y). We now consider the two-dimensional version of the RKDG scheme (10), using Q 2 elements (tensor product P 2 ⊗ P 2 ).

Accuracy results are shown in Table 4 and Figure 4 for time t = 0.5. The errors are computed globally (a uniform grid mesh of 50 2 points is used in each mesh cell to estimate the errors). The results in this example demonstrate that our scheme is also convergent in twodimensions. We will comment upon the two-dimensional case in the next section. The numerical solution and the exact solution are also plotted in Figure 5, showing good agreements. In this paper, we prove convergence of the SLDG and RKDG methods for the obstacle problem under the "no shattering" assumption of the exact solution. We utilize the DPP of the obstacle solutions. The proof of the SLDG methods relies on the property of the L 2 projection, while new techniques of devising piecewise intermediate stage functions are developed for the convergence of the RKDG methods.

We remark that the proposed methods can be easily extended to treat multidimensional obstacle equations. The simplest way is to define a multi-dimensional DG basis obtained as a tensor product of the one-dimensional DG basis. The Gaussian points on each cell can then be defined accordingly. The operation of taking the maximum at the Gaussian points is therefore straightforward as in the onedimensional case. The definition of the RKDG scheme in the multi-dimensional case is well known and is an extension of the one-dimensional case. On the other hand, the definition of the SLDG scheme in multi-dimensions is not straightforward. However, high-order stable splittings methods for SLDG for the linear advection equation can be devised (see for example Bokanowski and Simarmata [START_REF] Bokanowski | Semi-Lagrangian discontinuous Galerkin schemes for some first and second order partial differential equations[END_REF]). Finally the "no-shattering" property can be easily extended to multi-dimensions by demanding that the exact solution be piecewise regular except on a finite union of compact submanifolds, and error estimates of the same order as in the one-dimensional case will then hold. We refer to Example 2 in the previous section for the numerical performance of our RKDG methods in two-dimensions.

Remark 4 . 1 .

 41 Defining g as in (ii), and assuming min(ℓ, k) ≥ 1, the error is bounded by O( h 3/2

and 3 ξ n+1 2 -3 ξ n 2 ≤

 22 C∆th + C∆t ξ n 2

Figure 2 .

 2 Figure 2. Example 1, RKDG scheme with P 2 elements.

Figure 3 .

 3 Figure 3. Example 1, SLDG scheme, with P 2 elements. Left: ∆t = h/2, Right: ∆t = C h 3/5 .

Figure 4 .

 4 Figure 4. Example 2, Q 2 elements.

Figure 5 .

 5 Figure 5. Example 2 at time t = 0.5, numerical (left) and exact (right) data, using Q 2 elements with N x = 20 mesh cells

Table 1 .

 1 Example 1, RKDG scheme with P 2 elements.

	N x L 1 -error order L 2 -error order L ∞ -error order
	80 2.68E-04 1.52 1.01E-03 1.07 1.25E-02 0.27
	160 6.47E-05 2.05 3.46E-04 1.55 6.80E-03 0.87
	320 1.96E-05 1.72 1.30E-04 1.41 2.52E-03 1.43
	640 6.40E-06 1.62 5.65E-05 1.21 1.43E-03 0.82
	1280 2.10E-06 1.61 2.54E-05 1.15 8.26E-04 0.79
	2560 6.14E-07 1.77 1.02E-05 1.32 4.76E-04 0.79
	5120 1.98E-07 1.63 4.35E-06 1.22 2.75E-04 0.79
	10240 6.19E-08 1.68 1.88E-06 1.21 1.61E-04 0.78
	Example 2 (2-d, RKDG). This is a two-dimensional test, same as [4, Example
	3]). The equation solved is
	min

Table 2 .

 2 Example 1, SLDG scheme, with P 2 elements and ∆t = h/2. 38E-05 2.86 1.21E-04 2.35 8.58E-04 1.86 320 1.56E-05 0.61 9.77E-05 0.31 8.26E-04 0.05 640 4.73E-06 1.72 3.54E-05 1.47 4.30E-04 0.94 1280 1.73E-06 1.45 1.94E-05 0.87 3.25E-04 0.41 2560 4.11E-07 2.07 6.79E-06 1.52 1.76E-04 0.88 5120 1.31E-07 1.65 3.03E-06 1.16 1.19E-04 0.57 10240 4.03E-08 1.70 1.22E-06 1.32 6.67E-05 0.83

	N x L 1 -error order L 2 -error order L ∞ -error order
	80 1.73E-04	-	6.21E-04	-	3.11E-03	-
	160 2.

Table 3 .

 3 Example 1, SLDG scheme, with P 2 elements and ∆t = C h 3/5 .N x N L 1 -error order L 2 -error order L ∞ -error order 80 34 6.04E-05 1.36 1.79E-04 1.20 8.26E-04 0.94 160 52 1.99E-05 1.60 9.34E-05 0.94 7.58E-04 0.12 320 79 7.95E-06 1.33 2.66E-05 1.81 2.70E-04 1.49 640 121 2.50E-06 1.67 1.43E-05 0.90 1.86E-04 0.54 1280 183 9.80E-07 1.35 7.45E-06 0.94 1.46E-04 0.34 2560 278 3.66E-07 1.42 2.93E-06 1.35 8.03E-05 0.86 5120 422 1.22E-07 1.58 7.41E-07 1.98 2.74E-05 1.55 10240 639 4.48E-08 1.45 1.83E-07 2.02 5.82E-06 2.23

Table 4 .

 4 Example 2, Q 2 elements. N x = N y h x = h y L 1 -error order L 2 -error order L ∞ -error order

	10	2.00e-1 2.25E-02	-	2.27E-02	-	1.68E-01	-
	20	1.00e-1 6.70E-03 1.75 9.85E-03 1.20 1.24E-01 0.44
	40	5.00e-2 1.70E-03 1.98 3.29E-03 1.58 4.35E-02 1.51
	80	2.50e-2 5.11E-04 1.73 1.31E-03 1.33 2.67E-02 0.70
	160	1.25e-2 1.46E-04 1.80 5.19E-04 1.33 9.28E-03 1.52

Research of the first author is supported by the EU under the 7th Framework Programme Marie Curie Initial Training Network "FP7-PEOPLE-2010-ITN", SADCO project, GA number 264735-SADCO.. Research of the second author is supported by NSF grant DMS-1217563 and the start-up grant from Michigan State University. Research of the third author is supported by ARO grant W911NF-11-1-0091 and NSF grants DMS-1112700 and DMS-1418750.

where the spatial derivative (v n ) x should be understood in the weak sense. We notice that v (1) may become discontinuous at the irregular points of v n .

To define the second intermediate stage v (2) , we need to distinguish the "good" and "bad" intervals. Since v n ∈ C ℓ+1 p,L,c 0 (0, 1), when the mesh is fine enough, there are at most p irregular intervals. Because of the CFL condition (which we assume implies in particular that c△t ≤ h), each irregular point at t n may influence at most three intervals at time t n+1 . Now we introduce sets B n and I n such that B n := j I j , s.t. I j or its immediate neighbors contain an irregular point of v n and the corresponding set of indices:

Therefore meas(B n ) ≤ 3ph and Card(I n ) ≤ 3p. (In the case of the irregular points located exactly at the cell interface x j+1/2 , we include the point's neighboring cells I j , I j+1 in B n , and j, j + 1 in I n .)

Now, we define ṽ(2) :=

For points not located in B n , the definition coincides with [START_REF] Zhang | Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF]. For points in B n , v n is used instead of v (1) to avoid discontinuity at the irregular points. Notice that this causes ṽ(2) to be discontinuous at ∂B n . For example, if x a ∈ ∂B n , then the jump of ṽ(2) at x a is of magnitude c ∆t 4 (v (1) v n ) x (x a ), and this is bounded by CL∆t 2 . By these arguments, we could add a linear interpolating function defined by 1 4 L a (x) which is nonzero only on B n to enforce continuity at ∂B n , and ||L a || ∞ < C∆t 2 , i.e. we introduce

and L a is chosen to be a linear polynomial so that v (2) is continuous at ∂B n . We can now define

Notice that for x / ∈ B n , the definition is still consistent with [START_REF] Zhang | Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF] for smooth solutions, and it is well defined because ℓ ≥ 2. However, for irregular intervals, the definition is modified due to the lower regularity of the solution.