Convergence of DG schemes for obstacle problems

# Convergence of discontinuous Galerkin schemes for front propagation with obstacles

Olivier Bokanowski Laboratoire Jacques-Louis Lions, Université Pierre et Marie Curie 75252 Paris Cedex 05 France, and Université Paris-Diderot (Paris 7), 5 Rue Thomas Mann 75205, Paris CEDEX 13, France. Yingda Cheng Department of Mathematics, Michigan State University, East Lansing, MI, 48824 USA  and  Chi-Wang Shu Division of Applied Mathematics, Brown University, Providence, RI 02912, USA
February 21, 2015
###### Abstract.

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 , 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.

###### Key words and phrases:
Hamilton-Jacobi-Bellman equations; discontinuous Galerkin methods; level sets; front propagation; obstacle problems; dynamic programming principle; convergence
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.

## 1. 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,

 (1a) min(ut+cux, u−g(x))=0,x∈I=(0,1),t>0, (1b) u(0,x)=u0(x),x∈(0,1),

with periodic boundary conditions on and a constant . In (1), the function 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 (1). 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 [13] (see also [5] 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 [4], 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 [12]. As for the DG-HJ solvers, in [19, 23], the first efforts relied on solving the conservation law system satisfied by the derivative of the solution. See also [8] for an adaptive version of this scheme. In [9], a DG method for directly solving the Hamilton-Jacobi equation was developed and was later generalized to solve front propagation problems [3] and obstacle problems [4]. Other direct DG solvers include the central DG scheme [24] and the local DG scheme [29]. The schemes proposed in [4] 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 [30].

On the other hand, the semi-Lagrangian DG (SLDG) methods were proposed in [27, 28, 26] to compute incompressible flow and Vlasov equations, as well as in [7] 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 [6], but without convergence proof). For general SL methods, we refer to the works of Falcone and Ferretti [14, 15, 16], see also the textbook [17]. There is also a vast literature on high order finite difference schemes for solving HJ equations, see, e.g. [25, 1, 21, 22].

Beyond the scope of the present paper, yet of great interest, we also mention the second order PDE with obstacle terms, such as , with . This is the case of the so-called “American options” in mathematical finance [2]. Explicit schemes were proposed and proved to converge within the viscosity theory (see for instance [20]) yet with a reduced rate of convergence. Variational methods for nonlinear obstacle equations can also be devised [18] 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 (1). 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 [30], error estimates for RKDG methods for scalar conservation laws were provided for smooth solutions. In [11, 31], 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 is a space step and a time step, we shall show error bounds of the order of for the SLDG schemes (under time stepping , larger time step can be taken with a lower convergence rate), and of order for the RKDG schemes (under time stepping ). 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 regular for some for SLDG (resp. 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:

 (2) u(t,x)=max(u0(x−ct), maxθ∈[0,t]g(x−cθ)).

The function is also the solution of the Bellman’s dynamic programming principle (DPP): for any ,

 (3) u(t+Δt,x)=max(u(t,x−cΔt),maxθ∈[0,Δt]g(x−cθ)),∀t≥0.

Notice conversely that the DPP (3), together with , implies (2).

If we denote

 un(x):=u(tn,x)

and

 gt(x):=maxθ∈[0,t]g(x−cθ),

then the DPP implies in particular for any and :

 (4) un+1(x)=max(un(x−cΔt),gΔt(x)).

Using formula (2), we can see that when and are Lipschitz regular, then 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 is a discontinuous function (otherwise piecewise regular), formula (2) implies also some regularity on the solution .

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.

## 2. DG schemes for the obstacle equation

In this section, we will introduce the SLDG and RKDG methods for the obstacle equation (1). For simplicity of discussion, in the rest of the paper, we will assume to be a positive constant.

Let be a set of intervals forming a partition of . We denote where is the length of the interval . For a given integer , let be the DG space of piecewise polynomial of degree at most on each interval :

 (5) Vh:={v:I→R,  v|Ij∈Pk, ∀j}.

To introduce the methods for the obstacle problems, we follow two steps. Firstly, we describe the DG solvers for the linear advection equation .

To advance the numerical solution in one step from to , we consider one of the two DG methods described below.

SLDG Scheme:

 (6) vn+1h:=Πh(vnh(⋅−cΔt)),

where is the projection onto the space . We shall denote this SLDG solver by .

RKDG Scheme (by TVD-RK3 time stepping): find , , , such that

 (7a) (vn,1h−vnh,φh)=ΔtH(vnh,φh),∀φh∈Vh (7b) (vn,2h−34vnh−14vn,1h,φh)=Δt4H(vn,1h,φh),∀φh∈Vh (7c) (vn+1h−13vnh−23vn,2h,φh)=2Δt3H(vn,2h,φh),∀φh∈Vh

where

 (ϕ,φ)=∫Iϕφdx

and

 Hj(ϕh,φh)=∫Ijcϕh(φh)xdx−c((ϕh)−j+12(φh)−j+12−(ϕh)−j−12(φh)+j−12),
 H(ϕh,φh)=∑jHj(ϕh,φh).

We shall denote this RKDG solver by .

After introducing DG schemes for the linear transport problem, for the obstacle equation (1), we shall consider two approaches: one by using the projection

 (8) un+1h:=Πh(max(GΔt(unh), ~g)),

where or and

 (9) ~g≡gΔtor~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 .

Unfortunately, the scheme (8) is difficult to implement, because we need to compute the maximum of two functions, which requires locating the roots of . Another more practical approach is to define as the unique polynomial in such that

 (10) un+1h(xjα):=max(GΔt(unh)(xjα), ~g(xjα)),∀j=1,…,N,α=0,…,k,

where are the Gauss-Legendre quadrature points on the interval , and are the corresponding quadrature weights. Those schemes were studied in details and stability was established in [4]. We shall see in later sections that definitions (8) or (10) lead to similar error estimates, although the second approach is much easier to implement.

Finally, we remark that in [4], 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 , 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.

## 3. 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 (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 , the following function sets:

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

where denotes the -th derivative almost everywhere, and

 Δ2q(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 pseudo-norm definition will also be used:

 (11) ∥f∥ℓ2:=(∑i,αwiα|f(xiα)|2hi)1/2.

In particular, using the Gauss-Legendre quadrature rule, for any we have . From this point on, we will use to denote , and to denote 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 . For any function , , and for any element , it holds that

 (Phφ)−j+1/2=φ−j+1/2,∫Ij(Phφ−φ)ψhdx=0,∀ψh∈Pk−1(Ij).

In the lemma below, we will first establish the projection properties for functions in the space .

###### Lemma 3.1.

Let and let be in : is Lipschitz continuous, piecewise , with at most non regular points. Then there exists a constant , depending only on such that

 ∥φ−Phφ∥≤Chqk,ℓ,

where the projection can be either or , and is defined as

 (12) qk,ℓ:=min(min(k,ℓ)+1,32)≡{3/2if k,ℓ≥11if k=0 or ℓ=0
###### Proof.

Using the property of the projections [10], on each regular cell , we have

 ∥φ−Phφ∥Ij≤Chmin(ℓ,k)+1∥φ∥Hℓ+1(Ij)≤Chmin(ℓ,k)+1

since on . Therefore,

 ∥φ−Phφ∥L2(∪Ij, φ|Ijregular)=⎛⎝∫∪Ij, φ|Ijregular|φ−Phφ|2dx⎞⎠1/2≤Chmin(ℓ,k)+1.

On the other hand, on the intervals where is not regular, we have

 ∥φ−Phφ∥Ij≤Chmin(0,k)+1∥φ∥H1(Ij)≤Ch∥φ∥H1(Ij)≤Ch3/2

because is Lipschitz continuous. Hence summing up on the “bad” intervals (with at most such intervals),

Summing up the bounds with bad intervals and good ones, we prove the desired result. ∎

We now state a similar estimate for the norm:

###### Lemma 3.2.

Assume that belongs to , , then we have

 (13) ∥φ−Phφ∥ℓ2=(∑j,αwjα∣∣v(xjα)−(Phφ)(xjα)∣∣2hj)1/2≤Chqk,ℓ

where or , is defined as in (12), and the constant depends only on and .

###### Proof.

On each regular cell , we have [10],

 ∥φ−Phφ∥L∞(Ij)≤Chmin(ℓ,k)+1/2∥φ∥Hℓ+1(Ij)≤Chmin(ℓ,k)+1.

Then using the fact that , we obtain that

 ⎛⎜ ⎜⎝∑Ij, φ|Ij regular∑αwjα∣∣φn(xjα)−(Phφn)(xjα)∣∣2hj⎞⎟ ⎟⎠1/2≤Chmin(ℓ,k)+1.

On the other hand, when the interval is such that contains a non-regular point, we can write

 ∥φ−Phφ∥L∞(Ij)≤Chmin(0,k)+1/2∥φ∥H1(Ij)≤Ch.

Therefore

 ⎛⎜ ⎜⎝∑Ij, φ|Ij not regular∑αwjα∣∣φn(xjα)−(Phφn)(xjα)∣∣2hj⎞⎟ ⎟⎠1/2 ≤ ⎛⎜ ⎜⎝∑Ij, φ|Ij not regularh(Ch)2⎞⎟ ⎟⎠1/2 ≤ Cp1/2h3/2.

This concludes our proof. ∎

We now turn to some estimates related to the obstacle function .

###### Lemma 3.3.

Assume that for some integer . Let . There exists a constant such that

 (14) ∥~g−gΔt∥≤C√qΔt5/2.

and

 (15) ∥~g−gΔt∥ℓ2≤C√q √Δt+hΔt2.
###### Proof.

Denote the set of local maximum points of , if , then we see that . Furthermore

 (16) ∫{x, gΔt(x)≠~g(x)}dx ≤ ∫{x, [x−cΔt,x]∩Mg≠∅}dx ≤ qcΔt

since there are at most local maxima.

Consider the case when contains at least a local maxima of . Assuming that is small enough we can assume that is the only local maximum of on the interval , so that . Then,

 (17) |g(x)−gΔt(x)|=|g(x)−g(x∗)|≤C|x−x∗|2≤CΔt2,

since is twice differentiable at .

Combining (17) and (16) we obtain the bound (14).

For the second estimate (15), we make use of a minimal covering of the set

 {x, [x−cΔt,x]∩Mg≠∅},

using mesh intervals. The length of this covering is bounded by for each maximum point, since in order to cover any interval we may need two more mesh intervals , of length than the minimum required length . Overall the length of the total covering is bounded by , 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 , we will say that the exact solution of the problem (1) is “not shattering” if there exists some and constants such that the exact solution satisfies, for any , .

Recall that the exact solution satisfies . Therefore if is not shattering, it implies that has bounded number of zeros (since otherwise 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 , let and on , together with velocity constant . The function is of class on the interval . Notice then, for and , and the exact solution is therefore for and . So at time , (for ). This function has an infinite number of non regular points in the interval , and therefore does not satisfy the “no shattering” property at time .

It is not easy to state precise conditions on the initial data and to ensure that the no shattering property will be satisfied. Mainly, , 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 and . This definition still allows for a finite (bounded) number of singularities in , 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 ,

 (18a) g∈Cℓ+1pg,Lg,cg(0,1), (18b) u0∈Cℓ+1pu0,Lu0,cu0(0,1), and that, ∀t∈[0,T], (18c) x→gt(x)−u0(x−ct)has a finitely % bounded number of zeros in (0,1),

the bound being independent of .

Then there exists constants such that, , belongs to (with , , and that are independent of .)

###### Proof.

First, for , if we denote (resp. ) to be the number of local maxima (resp. minima) of , then we have . This is because the Lipschitz constant of is bounded by by an elementary verification. Each local maxima of may develop into two singularities in , each local minima may develop into one singularity in , and each singular point of may continue to be a singularity in . Hence the total number of non-regular points in will be bounded by . The bound of the derivative can also be obtained easily.

Because the exact solution is given by , the Lipschitz constant of is therefore bounded by .

On the other hand, the number of singular points of is bounded by the sum of the number of singular points of the function , plus the number of zeros of the same function, which is assumed to be bounded independently of . Hence the the number of singular points of is bounded independently of .

Finally the bound on the -partial derivative , in the regular region of , is easily obtained, because is then locally one of the two functions or . ∎

###### Lemma 3.5.

Assume that and satisfy the regularity assumptions (18a) and (18b) of Lemma 3.4. Assume furthermore that there exists a constant such that for a.e. . Then (18c) holds and satisfies the “no shattering” assumption.

###### Proof.

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

## 4. 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

 vt+cvx=0,v(0,x)=v0(x).

We will then generalize the results to scheme (8), and finally to scheme (10) for the obstacle problem.

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

We first consider the linear equation , for which

 v(t+Δt)=v(t,x−cΔt).

We denote , and we define the numerical solution of the SLDG method at to be . In particular, the scheme writes: initialize with , and for .

###### Theorem 4.1.

We consider , . If , then we have for all such that ,

 ∥vnh−vn∥≤CT hqk,ℓΔt

for some constant independent of , and is defined in (12).

###### Proof.

If , then , and due to Lemma 3.1, we have:

 ∥vn(⋅−cΔt)−Πh(vn(⋅−cΔt))∥≤Chqk,ℓ,

for some constant independent of . Hence

 ∥vn+1h−vn+1∥=∥Πh(vnh(⋅−cΔt)−vn(⋅−cΔt)∥ ≤∥Πh(vnh(⋅−cΔt)−Πh(vn(⋅−cΔt)∥+∥Πh(vn(⋅−cΔt)−vn(⋅−cΔt)∥ ≤∥Πh(vnh(⋅−cΔt)−Πh(vn(⋅−cΔt)∥+Chqk,ℓ ≤∥vnh(⋅−cΔt)−vn(⋅−cΔt)∥+Chqk,ℓ ≤∥vnh−vn∥+Chqk,ℓ,

where we have used the fact in the fourth row, and the periodic boundary conditions in the last row. As for the initial condition, we have

 ∥v0h−v0∥=∥Πhv0−v0∥≤Chqk,ℓ.

Finally we obtain for any given the existence of a constant (independent of and of ), such that

 ∥vnh−vn∥≤C(n+1)hqk,ℓ≤CT hqk,ℓΔt,

for all such that . ∎

Therefore, if , then , and we need for the convergence. If otherwise , it holds and we would only need for the convergence.

### 4.2. Convergence of the first SLDG scheme in the obstacle case

Now we turn to scheme (8) for the nonlinear equation (1). In particular, the scheme writes: initialize with , and for .

###### Theorem 4.2.

Assume that the exact solution is not shattering in the sense of Definition 3.1, for some integer . Then the following error bound holds:

 ∥unh−un∥≤CT hqk,ℓΔt+CT∥~g−gΔt∥Δt

for some constant independent of . In particular,

If , then :

 ∥unh−un∥≤CT hqk,ℓΔt

If , and if , then :

 ∥unh−un∥≤CT hqk,ℓΔt+CTΔt3/2.
###### Proof.

By Lemma 3.1 and the “no shattering” assumption,

 (19) ∥Πhun−un∥≤Chqk,ℓ.

Hence, from the definitions,

 ∥un+1h−un+1∥≤∥un+1h−Πhun+1∥+∥Πhun+1−un+1∥ ≤∥Πh(max(GSLΔt(unh),~g))−Πh(max(un(⋅−cΔt),gΔt))∥+Chqk,ℓ ≤∥max(GSLΔt(unh),gΔt)−max(un(⋅−cΔt),~g)∥+Chqk,ℓ.

Using the fact that

 |max(a1,b1)−max(a2,b2)|≤|a1−a2|+|b1−b2|,

we obtain

 ∥un+1h−un+1∥ ≤ ∥GSLΔt(unh)−un(⋅−cΔt)∥+∥~g−gΔt∥+Chqk,ℓ ≤∥Πh(unh(⋅−cΔt))−Πh(un(⋅−cΔt))∥+∥~g−gΔt∥+Chqk,ℓ (20) ≤∥unh−un∥+∥~g−gΔt∥+Chqk,ℓ

where we have used again (19) and the periodic boundary condition. Using Lemma 3.3 and by induction on , we are done. ∎

###### Remark 4.1.

Defining as in , and assuming , the error is bounded by . Therefore the optimal estimate is obtained when , or , and the error is of order .

### 4.3. Convergence of the SLDG scheme defined with Gauss-Legendre quadrature points

Now we turn to scheme (10) for the nonlinear equation (1). In particular, the scheme writes: initialize with , and is defined as the unique polynomial in such that:

 (21) un+1h(xjα):=max(GSLΔt(unh)(xjα), ~g(xjα)),∀j,α.
###### Theorem 4.3.

Let and assume that the exact solution is not shattering in the sense of Definition 3.1. The following error bound holds:

 ∥unh−un∥≤CT hqk,ℓΔt+CT∥~g−gΔt∥Δt,

for some constant independent of . In particular,

If , then :

 ∥unh−un∥≤CT hqk,ℓΔt.

If and , then :

 ∥unh−un∥≤CT hqk,ℓΔt+CTΔt√Δt+h.
###### Proof.

In view of Lemma 3.1 and the “no shattering” property, we have

 ∥un+1−Πhun+1∥≤Chqk,ℓ.

Now we turn back to the error estimate, and consider the case of :

 ∥un+1h−un+1∥≤∥un+1h−Πhun+1∥+Chqk,ℓ =∥un+1h−Πhun+1∥ℓ2+Chqk,ℓ ≤∥un+1h−un+1∥ℓ2+Chqk,ℓ ≤(∑j,αwjα∣∣un+1h(xjα)−un+1(xja)∣∣2hj)1/2+Chqk,ℓ

where in the third line we have used (13) for . Because of the DPP, for all points , the exact solution satisfies:

 un+1(x)=max(un(x−cΔt), gΔt(x)).

Therefore, for :

 ∣∣∣(un+1h−un+1)(xjα)∣∣∣≤