Alternating evolution discontinuous Galerkin methods for Hamilton ...

Journal of Computational Physics 258 (2014) 31–46

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Alternating evolution discontinuous Galerkin methods for Hamilton–Jacobi equations Hailiang Liu, Michael Pollack Iowa State University, Mathematics Department, Ames, IA 50011, United States

a r t i c l e

i n f o

Article history: Received 27 June 2013 Accepted 21 September 2013 Available online 8 October 2013 Keywords: Alternating evolution Hamilton–Jacobi equations Viscosity solution

a b s t r a c t In this work, we propose a high resolution Alternating Evolution Discontinuous Galerkin (AEDG) method to solve Hamilton–Jacobi equations. The construction of the AEDG method is based on an alternating evolution system of the Hamilton–Jacobi equation, following the previous work Liu et al. (2013) [31] on AE schemes for Hamilton–Jacobi equations. A semidiscrete AEDG scheme derives directly from a sampling of this system on alternating grids. Higher order accuracy is achieved by a combination of high-order polynomial approximation near each grid and a time discretization with matching accuracy. The AEDG methods have the advantage of easy formulation and implementation, and efficient computation of the solution. For the linear equation, we prove the L 2 stability of the method. Numerical experiments for a set of Hamilton–Jacobi equations are presented to demonstrate both accuracy and capacity of these AEDG schemes. © 2013 Elsevier Inc. All rights reserved.

1. Introduction In this paper, we develop a new alternating evolution discontinuous Galerkin (AEDG) method to solve time-dependent Hamilton–Jacobi (H–J) equations. We describe the designing principle of our AEDG scheme through the following form:

φt + H (x, φ, ∇x φ) = 0,

φ(x, 0) = φ0 (x),

x ∈ Rd , t > 0.

(1.1)

Here, d is the space dimension, the unknown φ is scalar, and H : R → R is a nonlinear Hamiltonian. The Hamilton–Jacobi equation arises in many applications ranging from geometrical optics to differential games. These nonlinear equations typically develop discontinuous derivatives even with smooth initial conditions, the solutions of which are non-unique. In this paper, we are only interested in the viscosity solution [11,12,34], which is the physically relevant solution in some important applications. The difficulty encountered for the satisfactory approximation of the exact solutions of these equations lies in the presence of discontinuities in the solution derivatives. An important class of finite difference methods for computing the viscosity solution is the class of monotone schemes introduced by Crandall–Lions [13]. Unfortunately, monotone schemes are at most first-order accurate. The need for devising more accurate numerical methods for Hamilton–Jacobi equations has prompted the abundant research in this area in the last two decades, including essentially non-oscillatory (ENO) or weighted ENO (WENO) finite difference schemes, see e.g., [15,17,21,27,28,33,36], and Central or Central-Upwind finite difference schemes, see e.g., [2–4,19,20,22,23], as well as discontinuous Galerkin methods [5,16,25,35]. The preliminary goal of this paper is to design high-order discontinuous Galerkin methods for Hamilton–Jacobi equations by following the alternating evolution framework introduced in [24] for hyperbolic conservation laws; see local AE schemes developed by Saran and Liu [30], see also [1]. The high resolution finite difference scheme under this AE framework was d

E-mail addresses: [email protected] (H. Liu), [email protected] (M. Pollack). 0021-9991/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.09.038

1

32

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

further developed for Hamilton–Jacobi equations in [31]. In such a framework, the general setting is to refine the original PDE by an alternating evolution (AE) system which involves two representatives: {u , v }. To update u, terms involving spatial derivatives are replaced by v’s derivatives and augmented by an additional relaxation term ( v − u )/ , which serves to communicate the two representatives. The AE system for scalar hyperbolic conservation laws was shown in [24] to be capable of capturing the exact solution when initially both representatives are chosen as the given initial data. Such a feature allows for a sampling of two representatives over alternating grids. Using this alternating system as a ‘building base’, we apply standard approximation techniques to the AE system. For Hamilton–Jacobi equations, we consider the following alternating evolution system

ut + H (x, u , ∇x v ) = v t + H (x, v , ∇x u ) = Here,

1

 1



( v − u ), (u − v ).

 > 0 is a scale parameter of user’s choice. At t = 0, we take the initial data u 0 (x) = v 0 (x) = φ0 (x).

We sample the above system over alternating grids when performing the spatial discretization, leading to a class of high order semi-discrete AEDG schemes. The semi-discrete scheme thus obtained has the form

d





dt Iα



H x, Φ, ∇x Φ S N

Φ η dx +



η dx + λ[Φ]η|xα +

Iα

1



Φ η dx =

 Iα

1



Φ S N η dx,

 Iα

where λ = H p (x, Φ(x), ∇x Φ), and Φ are polynomial elements sampled from neighboring cells I α ±e j . The scale parameter  and the time step size t are chosen to stabilize the time discretization. In the one-dimensional case with H = β p, the semi-discrete AEDG scheme of first-order becomes SN

d dt



Φα + β

Φ α + 1 − Φα − 1 2 x



=

1 2

(Φα +1 − 2Φα + Φα −1 ),

 = t, reduces to the celebrated Lax–Friedrichs scheme   n Φα +1 − Φαn −1  1 . Φαn+1 = Φαn +1 + Φαn −1 −  β 2 2 x

which, using Euler-forward with

We recall that the success of finite difference schemes for Hamilton–Jacobi equations is due to two factors: the local enforcement of the equation in question and the non-oscillatory piecewise polynomial interpolation from evolved grid values. The difference in the AE schemes proposed in [31] and the existing finite difference schemes lies in the local enforcement. The ENO/WENO type schemes [17,21,27,28,32,33,36] are based mainly on some local refinement of H–J equations by

  φt + Hˆ x, φx+ , φx− = 0,

(1.2)

ˆ is the numerical Hamiltonian which needs to be carefully chosen to ensure that the viscosity solution is capwhere H tured when φx becomes discontinuous. In comparison, the central type schemes [2–4,19,20,22,23] choose to evolve the constructed polynomials in smooth regions such that the Taylor expansion may be used in the scheme derivation. Finite difference methods are quite efficient for Cartesian meshes. However, on unstructured meshes, the scheme is quite complicated. Alternatively, the discontinuous Galerkin (DG) method, originally devised to solve first-order partial differential equations such as conservation laws [6–10,29], has the advantage of flexibility for arbitrarily unstructured meshes, and with the ability to easily achieve arbitrary order of accuracy. However, some new difficulties occur when the existing ideas with finite difference methods are applied toward some discontinuous Galerkin discretization. One main difficulty comes from the non-conservative form of the equation, which precludes the use of integration by parts to establish the cell to cell communication via numerical fluxes as usually done with the DG methods for hyperbolic conservation laws. In spite of this difficulty, some progress has been made. In [16], Hu and Shu proposed a discontinuous Galerkin method to solve the Hamilton–Jacobi equation φt + H (∇x φ) = 0. They use the fact that the derivative of the solution φ satisfies a conservation law system, and apply the usual discontinuous Galerkin method on this system to advance the derivatives of φ . The solution φ itself is then recovered from these derivatives by a least square procedure for multi-dimensional cases and with an independent evolution of the cell averages of φ . However, in the multi-dimensional case, a single equation of φ is converted to a system of w = ∇x φ , which is only weakly hyperbolic at some points. Also, the resulting algorithm appears indirect and complicated. This motivates the development of another DG method in [5] for directly solving H–J equations; for nonlinear convex Hamiltonians, their numerical experiments demonstrate that the method is stable and accurate, however, entropy correction is needed near the kink formation in some cases. In [35], Yan and Osher apply the LDG discretization to (1.2)

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

33

by introducing two additional variables p 1 = φx+ and p 2 = φx− and then apply the usual discontinuous Galerkin method to advance (φ, p 1 , p 2 ) all together in the one-dimensional case, with an extension to the two-dimensional case. In [25], Li and Yakovlev proposed a central discontinuous Galerkin method that directly solves the Hamilton–Jacobi equations using two polynomial representatives solved on two sets of overlapping meshes [26]. The AEDG scheme proposed in this paper only has one polynomial representative associated with each grid point, even in multi-dimensional case. The scheme that we design here is also different in that it is based on alternatively sampling the AE system of the Hamilton– Jacobi equation with spatial accuracy enhanced by interlaced polynomials. However, the AEDG is also similar to the central DG scheme in the sense that neighboring polynomials in the AEDG (the other representative polynomials in the central DG) are used whenever a calculation of the Hamiltonian is needed. For the one-dimensional Hamiltonian H (x, φ, φx ), the two schemes are almost equivalent when rescaling the mesh size x to 2x, except that in the AEDG scheme the neighboring polynomials are used to compute only the solution gradient φx , while the central DG method uses the other representative polynomials to compute the entire H (x, φ, φx ). The main advantage of the AEDG scheme is that in the multi-dimensional case, the alternating evolution framework results in easy formulation of the computation of the Hamiltonian. For example, in the two-dimensional case, AEDG evaluates the gradient using the closest polynomials. It should be noted that even though our AEDG schemes are derived based on sampling the alternating evolution system, we do not solve the system directly. The AE system simply provides a systematic way for developing numerical schemes of both semi-discrete and fully discrete form for the original problem, instead of as an approximation system at the continuous level. The article is organized as follows: in Section 2, we formulate the AEDG method for one-dimensional Hamilton–Jacobi equations. We also discuss the use of kink detector, and reconstruction of polynomial elements to make the approximation less oscillatory, so that the computed solution is the viscosity solution as desired. We then give a rigorous proof of L 2 stability by the AEDG method for H = α p in Section 3. Extensions to the multi-D case is given in Section 4. In Section 5, we show numerical results, which includes both one and two-dimensional problems. The results illustrate accuracy, efficiency, and high resolution near kinks. The last Section 6 ends this paper with our concluding remarks. 2. Alternating evolution DG methods Our numerical schemes for H–J equations consist of a semi-discrete formulation based on sampling of the AE system on alternating grids and a fully discrete version by using an appropriate Runge–Kutta solver. To illustrate, we start with the one-dimensional Hamilton–Jacobi equation of the form

φt + H (x, φ, φx ) = 0. The ‘building base’ is the following AE system

ut + H (x, u , v x ) = v t + H (x, v , u x ) =

1

 1



( v − u ),

(2.1)

(u − v ).

(2.2)

2.1. Semi-discrete AEDG formulation We develop an AE discontinuous Galerkin (AEDG) method for the H–J equation subject to initial data φ0 (x). We divide the spatial domain to form a uniform grid, and let {x j } be the grid point, with x = x j +1 − x j . We denote I j = (x j −1 , x j +1 ) for j = 2, . . . , N − 1 with I 1 = (x1 , x2 ) and I N = (x N −1 , x N ). Centered at each grid {x j }, the numerical approximation is a polynomial Φ| I j = Φ j (x) ∈ P k , where P k denotes a linear space of all polynomials of degree at most k:

  

P k := p  p (x)| I j =



a i ( x − x j )i , a i ∈ R .

0i k

Note dim( P k ) = k + 1. We denote v (x± ) = lim →0± v (x +  ), and v ± = v (x±j ). The jump at x j is [ v ] j = v (x+j ) − v (x−j ), and j

− the average { v } j = 12 ( v + j + v j ). Integrating the AE system (2.1) over I j against the test function semi-discrete AEDG scheme





    1 ∂t Φ j + H x, Φ j , ∂x Φ Sj N η dx + λ Φ Sj N x η(x j ) = j

Ij



 Φ Sj N η dx −



Ij

η ∈ P k ( I j ), we obtain the

 Φ j η dx ,

Ij

where λ = H p (x j , Φ j (x j ), ∂x Φ j (x j )) and Φ Sj N are sampled from neighboring polynomials Φ j ±1 in the following way:

(2.3)

34

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

 H



x, Φ j , ∂x Φ Sj N



x j

η dx =

x j +1 

H (x, Φ j , ∂x Φ j −1 )η dx +

x j −1

Ij





x j +1

Φ j −1 η dx +

x j −1

Ij

xj



xj

Φ Sj N η dx =

H (x, Φ j , ∂x Φ j +1 )η dx,

Φ j +1 η dx, xj

     S N  Φ j x = Φ j +1 x+j − Φ j −1 x−j . j

To update each grid-centered polynomial element Φ , we write the compact form of the semi-discrete scheme

d



 Φ j η dx = L Φ j ; Φ Sj N , η (t ),

dt

(2.4)

Ij

where

L [Φ j ; Φ j ±1 , η](t ) =

1







 Φ Sj N − Φ j η dx −

Ij





H x, Φ j , ∂x Φ Sj N







η dx − λ Φ Sj N x j η(x j ).

(2.5)

Ij

For cells near the boundary we treat with care. For a computational domain [a, b] with x1 = a, x N = b and x = (b − a)/( N − 1), the two equations near boundary may be given as

x2 x2   λ 1 ∂t Φ1 + H (x, Φ1 , ∂x Φ2 ) η dx + [Φ]|x1 η(x1 ) = (Φ2 − Φ1 )η dx,



2

x1

xN

x1

xN



 λ 1 ∂t Φ N + H (x, Φ N , ∂x Φ N −1 ) η dx + [Φ]|xN η(xN ) = 2

x N −1

(2.6)



(Φ N −1 − Φ N )η dx.

(2.7)

x N −1

At the left boundary, η ∈ P (x1 , x2 ), and on the right boundary, η ∈ P k (x N −1 , x N ), so that the length of the interval is x instead of 2x. If the flow is incoming at x = a, one has to impose a boundary condition φ(a, t ) = g 1 (t ). As a result, one is required to modify (2.6) by changing [Φ] to Φ2 (x+ 1 ) − g 1 (t ); for the outflow case, one may take [Φ] = 0 in (2.6). Similarly, at x = b, the inflow boundary condition φ(b, t ) = g 2 (t ) can be incorporated in (2.7) by replacing [Φ] by g 2 (t ) − Φ N −1 (x− N , t ); and for outgoing flow, replacing [Φ] by 0. The determination of the inflow or outflow of the boundary condition may be obtained by checking the sign of the characteristic speed ∂x H (x, φ, p ). In the case of periodic boundary conditions, [Φ] can be computed − as Φ2 (x+ 1 ) − Φ N −1 (x N ) at x = x1 , x N , and Φ N (x) is regarded to be identical to Φ1 (x). The fully discrete scheme follows from applying an appropriate Runge–Kutta solver to (2.4). We summarize the algorithm as follows. k

Algorithm. 1. Initialization: in any cell I j , compute the initial data by the local L 2 -projection





 Φ 0 − φ0 η dx = 0,

η ∈ P k ( I j ).

(2.8)

Ij

2. Alternating evaluation: take polynomials Φ j ±1 (x) = Φ| I j±1 , and then sample in I j to get L [Φ j ; Φ j ±1 , η] as defined in (2.5). 3. Evolution: obtain Φ n+1 from Φ n by some Runge–Kutta type procedure to solve the ODE system (2.4). In the AEDG schemes,

t    Q

 is chosen such that the stability condition,

x , max | H p (x, p )|

is satisfied. The choice of Q depends on the order of the scheme, similar to that used in the AE finite difference scheme introduced in [31].

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

35

2.2. Discrete ODEs We now present AE schemes with polynomial elements of degree k for H–J equations. Let {φl (ξ )}lk=+11 be the basis in the master domain [−1, 1], then in each cell I j we can express k +1

Φ j (x) =

alj φl (ξ ) =: φ  (ξ )a j ,

ξ=

l =1

x − xj

x

,

then

Φ j ±1 (x) =

k +1 l =1

alj ±1 φl (ξ ∓ 1) = φ  (ξ ∓ 1)a j ±1 .

A simple calculation of (2.3) gives

1 1 M a˙ j = − Ma j + L j ,



(2.9)



where

1 M =h

φ(ξ )φ  (ξ ) dξ,

L j = C − a j −1 + C + a j +1 −  hH j

−1

with

0 C− = h

φ(ξ + 1)φ  (ξ ) dξ +  λφ(0)φ  (1),

−1

1 C+ = h

φ(ξ − 1)φ  (ξ ) dξ −  λφ(0)φ  (−1),

0

0 Hj =

  φ(ξ ) H x j + hξ, φ  (ξ )a j , h−1 ∂ξ φ  (ξ + 1)a j −1 dξ

−1

1 +

  φ(ξ ) H x j + hξ, φ  (ξ )a j , h−1 ∂ξ φ  (ξ − 1)a j +1 dξ.

0

The Hamiltonian H is generally nonlinear, some quadrature rule must be applied to evaluate H j . The parameter λ is used to ensure the stability. The following two options are acceptable:

  1 λ = λ j = H p x j , Φ j (x j ), (∂x Φ j +1 + ∂x Φ j −1 )(x j ) 2

or

  λ = H p x j , Φ j (x j ), ∂x Φ j (x j ) . 2.3. Troubled cell detection and limiters It is known that near singularities one may need to impose some limiters so that the proposed scheme can capture the viscosity solution. To determine which cells are troubled cells and require an application of a limiter, we use a strategy that is based on a strong superconvergence as in [18]. We construct the troubled cell detector in one-dimensional case:

Ij =

|∂x φ j −1 (x−j ) − ∂x φ j +1 (x+j )| h

k −1 2

φL 2 ( I j )

.

The order of the numerator can be determined by the continuity of φx at the interface x = x j ,

 k      ∂x φ j −1 x− − ∂x φ j +1 x+  = O(h ), if φx is smooth, j j O(1), if φx is discontinuous,

36

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

so that

Ij =

k

O(h 2 ), O(h

− 2k

if φx is smooth,

), if φx is discontinuous,

for k  1. Then, I j → 0 as h → 0 in smooth regions and I j → ∞ as h → 0 in nonsmooth regions. So, if I j > 1, then φx is . If I j < 1, discontinuous and the cell I j is identified as a troubled cell. A limiter is then applied to that cell to obtain φ new j

then φx is continuous so that φ new = φj. j In order to capture the viscosity solution, nonlinear limiters are applied (see [25]) to the detected troubled cells. For k = 1, φ j is written as

 φ j = a1j + a2j

x − xj

 .

x

In order to update the a2j term, we consider two new candidates for the slope of φ j ,

s Lj = a1j − a1j −1 ,

s Rj = a1j +1 − a1j .

We then compute

   a2j = MM a2j , s Lj , s Rj . Here M M denotes the Min–Mod nonlinear limiter

⎧ ⎨ min j {b j } if b j > 0, ∀ j , MM{b1 , b2 , . . .} = max j {b j } if b j < 0, ∀ j , ⎩ 0 otherwise.

If |a2j − a2j |   for small For k = 2,

 φ j = a1j + a2j

 , then no change is needed. Otherwise, a2j = a2j .

x − xj



 + a3j

x

x − xj

x

2 .

The application of the limiter to a2j is the same as in the k = 1 case. Once the a2j terms have been updated, the a3j terms x−x

are updated in a similar fashion, but for the polynomial x∂x φ j = a2j + 2a3j ( x j ). The second derivative approximations can be computed as

 L sj

 R

= a2j − a2j −1 ,

sj

= a2j +1 − a2j .

We then compute

  L   R   . a3j = MM 2a3j , s j , s j If |a3j − a3j |   for small

 , then no change is needed. Otherwise, a3j = a3j /2.

2.4. Time discretization We now turn to time discretization of (2.9), following [28]. Let {t n }, n = 0, 1, . . . , K be a uniform partition of the time interval [0, T ]. Let Φ 0 = P φ0 be the piecewise L 2 projection of φ0 (x) defined in (2.8), and A = [a1 , . . . , a N ] , then

A t = G ( A ). We use the third order explicit SSP Runge–Kutta method [14] for time discretization. In details, let A n be the solution at time level n,





A (1) = A n + tG A n ,

  3 ( 0) 1  ( 1 ) A + tG A (1) , A + 4 4   1 ( 0) 2  ( 2 ) n +1 A = A + A + tG A (2) . 3 3 A (2 ) =

(2.10)

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

37

3. Stability analysis In this section we show that the semi-discrete AEDG scheme is L 2 stable for linear Hamiltonian. In particular, when H = α p, α = const., we have the following. Theorem 3.1. Let Φ be computed from the AEDG scheme (2.4) for the Hamilton–Jacobi equation φt + H (φx ) = 0 with linear Hamiltonian H ( p ) = α p and periodic boundary conditions. Then

d

x j +1  N −1  Φ 2j +1 + Φ 2j

dt

2

j =1 x j

 dx = −

x N −1 j +1 1



(Φ j +1 − Φ j )2 dx.

j =1 x j

η = Φ2i =: u in (2.4) for j = 2i, and η = Φ2i+1 =: v in (2.4) for j = 2i + 1, respectively, it gives     u2 1 H (x, u , v x )u dx + λ[ v ]u |x2i = uv dx − u 2 dx , dx +

Proof. Taking



d dt



2

I 2i

I 2i



d

v

dt

2

I 2i



dx +

2 I 2i +1

H (x, v , u x ) v dx + λ[u ] v |x2i+1 =

1

I 2i





I 2i +1



I 2i +1

 v 2 dx .

uv dx − I 2i +1

Note that N /2−1



[a, b] =

N /2−1

I 2i + [x N −1 , x N ] = [x1 , x2 ] +

i =1



I 2i +1 ,

i =1

we sum over all index and use (2.6)–(2.7) to obtain

d

b

u2

dt

2

dx +

 N /2−1 

d

v2

dt

2

+

i =1 I 2i

a

b

xN 

 x2 dx +

+

a

N /2−1

H (x, u , v x )u dx +

H (x, v , u x ) v dx +

i =1 I 2i +1

1 2

1

1

2



λ[ v ]u |x2i + λ[ v ]u |xN =

i =1

x N −1

 N /2−1 

x1



N /2−1

λ[u ] v |x1 +

i =1

λ[u ] v |x2i+1 =

b



a

1



b a

Adding these two relations up, we obtain

d

b

u2 + v 2

dt

2

dx + J 1 + J 2 = −

a

1



b (u − v )2 dx, a

where for H ( p ) = α p all other terms in J 1 and J 2 are expressed as follows:

J1 =

 N /2 x2i i =1x2i −1

x2i +1  

N /2−1

+

i =1



H (x, u , v x )u + H (x, v , u x ) v dx

x2i

 N /2−1  N /2    − = −α [u ] v |x2i+1 + [ v ]u |x2i − α Φ2 x+ 1 Φ1 (x1 ) + α Φ N (x N )Φ N −1 x N , i =1

i =1

and

 N /2−1  N /2  −     λ J2 = λ [u ] v |x2i+1 + [ v ]u |x2i + Φ1 (x1 ) + Φ N (xN ) Φ2 x+ 1 − Φ N −1 x N . i =1

i =1

2

Using the periodic boundary condition Φ1 (x1 ) = Φ N (x N ), we have

J1 + J2 = 0 when λ = α .

2



uv − u 2 dx,





uv − v 2 dx.

38

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

4. The multi-dimensional case By similar procedures we can construct AEDG schemes for multi-dimensional H–J equations:

φt + H (x, φ, ∇x φ) = 0,





x := x1 , . . . , xd ∈ Rd .

We start with the AE formulation

ut +

1



u=

1



v − H (x, u , ∇x v ).

(4.1)

Let {xα } be distributed grids in Rd . Consider I α be a hypercube centered at xα with vertices at xα ±1 where the number of vertices is 2d . Centered at each grid {xα }, the numerical approximation is a polynomial Φ| I α = Φα (x) ∈ P r , where P r denotes a linear space of all polynomials of degree at most r in all xi :

  

P r := p  p (x)| I α =





aβ (x − xα )β , 1  i  d, aβ ∈ R .

0βi r

Note dim( P r ) = (r + 1) . Sampling the AE system (4.1) in I α , we obtain the semi-discrete AEDG scheme d



   1 ∂t Φα + H x, Φα , ∇x ΦαS N η dx +



Φα η dx =



Iα

1

d  j =1 I

SN

where Φα shall take



ΦαS N η dx − B ,



Iα

B=



(4.2)

Iα j



H j (x, Φα , ∇Φα )[Φ]η(x)δ x j − xα j dx,

(4.3)

α

are sampled from neighboring polynomials, and H j = ∂ p j H (x, z, p ). The choice of ΦαS N is not unique, and we

ΦαS N ∈ span{Φα ±e j }dj=1 . In the two-dimensional case with y j +1 xi+1 



H x, y , Φ, ∂x Φ

SN

α = (i , j ), the terms involving neighboring polynomials are as follows:

, ∂yΦ

x i −1 y j −1

SN



y j +1 xi+1 

η dy dx =

H (x, y , Φi , j , ∂x Φi +1, j , ∂ y Φi , j +1 )η(x, y ) dy dx xi

yj

xi+1 y j +

H (x, y , Φi , j , ∂x Φi +1, j , ∂ y Φi , j −1 )η(x, y ) dy dx x i y j −1 y j +1 xi 

+

H (x, y , Φi , j , ∂x Φi −1, j , ∂ y Φi , j +1 )η(x, y ) dy dx

x i −1 y j

xi y j +

H (x, y , Φi , j , ∂x Φi −1, j , ∂ y Φi , j −1 )η(x, y ) dy dx.

x i −1 y j −1

An average of two neighboring polynomials Φi ±1, j and Φi , j ±1 will be used to evaluate

 Φ

SN

η dx =

1 2

y j +1 xi+1  (Φi +1, j + Φi , j +1 )η(x, y ) dy dx xi

Iα

+

1

yj

xi+1 y j (Φi +1, j + Φi , j −1 )η(x, y ) dy dx

2 x i y j −1

+

1 2

y j +1 xi  (Φi −1, j + Φi , j +1 )η(x, y ) dy dx x i −1 y j

 Iα

Φ S N η dx, that is

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

+

39

xi y j

1

(Φi −1, j + Φi , j −1 )η(x, y ) dy dx.

2 x i −1 y j −1

The B term in (4.3) can then be computed as y j +1 

B=









H 1 (x, y , Φi , j , ∂x Φi , j , ∂ y Φi , j ) Φi +1, j x+ , y − Φi −1, j x− ,y i i



η(xi , y ) dy

y j −1

xi+1 +









− Φi , j −1 x, y −j H 2 (x, y , Φi , j , ∂x Φi , j , ∂ y Φi , j ) Φi , j +1 x, y + j



η(x, y j ) dx.

x i −1

For the boundaries, we first consider the side boundary along x = x1 . Then y j +1 x2 



H x, y , Φ, ∂x Φ S N , ∂ y Φ

 SN

y j +1 x2 

η dy dx =

x1 y j −1

H (x, y , Φ1, j , ∂x Φ2, j , ∂ y Φ1, j +1 )η(x, y ) dy dx x1

yj

x2 y j +

H (x, y , Φ1, j , ∂x Φ2, j , ∂ y Φ1, j −1 )η(x, y ) dy dx, x1 y j −1

and y j +1 y j +1 x2  x2  1 SN Φ η dy dx = (Φ2, j + Φ1, j +1 )η(x, y ) dy dx

2

x1 y j −1

x1

+

yj

x2 y j

1

(Φ2, j + Φ1, j −1 )η(x, y ) dy dx,

2 x1 y j −1

and

B=

y j +1 

1 2









− H 1 (x, y , Φ1, j , ∂x Φ1, j , ∂ y Φ1, j ) Φ2, j x+ 1 , y − Φ0, j x1 , y



η(x1 , y ) dy

y j −1

x2 +









− H 2 (x, y , Φ1, j , ∂x Φ1, j , ∂ y Φ1, j ) Φ1, j +1 x, y + j − Φ1, j −1 x, y j



η(x, y j ) dx.

x1

In the above, Φ0, j (x− 1 , y ) may be taken different ways: the given boundary data at x1 for inflow boundary; Φ2, j (x1 , y ) for outgoing flow; and Φ N x −1, j (x− N x , y ) for periodic boundary conditions. Similar computations are made along the other side boundaries x = x Nx , y = y 1 , y = y N y . At the corners, we first consider the southwest corner (x1 , y 1 ). Then

x2 y 2



H x, y , Φ, ∂x Φ

SN

, ∂yΦ

x1 y 1

SN



x2 y 2

η dy dx =

H (x, y , Φ1,1 , ∂x Φ2,1 , ∂ y Φ1,2 )η(x, y ) dy dx, x1 y 1

and

x2 y 2 Φ x1 y 1

and

SN

η dy dx =

1 2

x2 y 2 (Φ2,1 + Φ1,2 )η(x, y ) dy dx, x1 y 1

40

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

B=

1

y 2









− H 1 (x, y , Φ1,1 , ∂x Φ1,1 , ∂ y Φ1,1 ) Φ2,1 x+ 1 , y − Φ0,1 x1 , y

2



η(x1 , y ) dy

y1

+

x2

1 2









− H 2 (x, y , Φ1,1 , ∂x Φ1,1 , ∂ y Φ1,1 ) Φ1,2 x, y + 1 − Φ1,0 x, y 1



η(x, y 1 ) dx.

x1

Similar computations can be made for the other corners (x1 , y N y ), (x Nx , y 1 ), (x Nx , y N y ). Boundary conditions are incorporated in the following ways:

⎧ ⎪ boundary data inflow,  −  ⎨ + outflow, Φ0,1 x1 , y = Φ2,1 (x1 , y ) ⎪ − ⎩Φ N x −1,1 (x N x , y ) periodic,

⎧ inflow, ⎪ ⎨ boundary+data  − ( x , y ) outflow , Φ 1 , 2 Φ1,0 x, y 1 = 1 ⎪ ⎩ Φ1, N −1 (x1 , y − ) periodic. y Ny

The SSP Runge–Kutta method [14] can be used to achieve a time discretization with matching accuracy. 5. Numerical tests In this section, we use some model problems to numerically test the performance of the proposed AEDG scheme in both one dimension and two dimensions. In what follows, we use N to denote the number of grid points the domain is divided into, i.e., x = (b − a)/( N − 1) for domain [a, b]. The same notation will be also used in two-dimensional case as long as the same partition number is used in both x and y direction. Letting the time step at tn be tn and Q to be the CFL number, we then take  as

tn <  =

Q x max | H p |

.

In the two-dimensional case,

 tn <  = Q

max | H 1 |

x

+

max | H 2 |

−1 .

y

The practical choice for Q is about 0.5 or smaller as k increases. We compute the errors as

e L 1 =

     φi (x) − φ ex (x) dx ,

e L ∞

(5.1)

j

i

j

j

Ki

   = max maxφi (x j ) − φ ex j (x j ) , i

(5.2)

j

where the indice i corresponds to the computed solution for x ∈ K i = [xi −1/2 , xi +1/2 ] and j corresponds to the exact solution, j

or the reference solution with a finer mesh K i such that

N

j =1

j

Ki = Ki .

Example 5.1. We first test the numerical accuracy of the designed schemes by using the Hamilton–Jacobi equation with convex Hamiltonian and with smooth initial data. The equation is

Φt +

(Φx + 1)2 2

= 0,

−1  x  1,

with initial data

Φ(x, 0) = − cos(π x). In Tables 1–2, the numerical accuracy results using P 1 and P 2 polynomials are presented when the solution is smooth (time T = 0π.52 ). All the proposed methods give the desired order of accuracy.

The solution becomes discontinuous at time T = π12 , and in Fig. 1, we plot the solutions at time T = 1π.52 . For the numerical experiments the CFL number is taken to be 0.5 and t = 0.5 .

Example 5.2. The following example is with a non-convex Hamiltonian,

1





Φ 2 − 1 Φx2 − 4 = 0, 4 x Φ(x, 0) = −2|x|. Φt +

−1  x  1,

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

41

Table 1 The L 1 and L ∞ error for convex Hamiltonian, Example 5.1 using N equally spaced grid points for AEDG scheme using P 1 polynomials at T = 0π.25 , when the solution is continuous. N −1

Scheme

L 1 error

L 1 order

L ∞ error

L ∞ order

20 40 80 160 320 640 1280

AEDG AEDG AEDG AEDG AEDG AEDG AEDG

0.019367230146 0.005086873182 0.001192497567 0.000282925725 0.000068076695 0.000016043798 0.000003700941

1.92876661 2.09279280 2.07549106 2.05519046 2.08514528 2.11605179

0.029367856297 0.012045086833 0.003473183828 0.000859545520 0.000214389218 0.000053242772 0.000013150665

1.28579317 1.79411412 2.01461282 2.00334169 2.00957476 2.01744992

Table 2 The L 1 and L ∞ error for convex Hamiltonian, Example 5.1 using N equally spaced grid points for AEDG scheme using P 2 polynomials at T = 0π.25 , when the solution is continuous. N −1 20 40 80 160 320 640

Scheme

L 1 error

L 1 order

L ∞ error

L ∞ order

AEDG AEDG AEDG AEDG AEDG AEDG

0.001261711385 0.000143734248 0.000016364650 0.000002050953 0.000000234804 0.000000027526

3.13390617 3.13474920 2.99621664 3.12676552 3.09258614

0.008085223415 0.000930253332 0.000111661585 0.000022752723 0.000003160315 0.000000362958

3.11959208 3.05849071 2.29502185 2.84789905 3.12219537

Fig. 1. Comparison of plots for H–J equation with convex Hamiltonian, Example 5.1, at discontinuity on [−1, 1], T = 1.5/π 2 , N = 81, using P 1 and P 2 polynomials.

This test problem is used to show resolution of discontinuities when the initial data has discontinuous derivatives. In testing our AEDG methods, the CFL number used for all order schemes was 0.5. Linear extension boundary conditions are used for the numerical tests. In Fig. 2, we can see that the AEDG scheme with the application of a limiter provides the desired viscosity solution. In Fig. 3, a comparison of the scheme with the applied limiter is shown for P 1 and P 2 polynomials. Example 5.3. We next solve a linear Hamiltonian with kink solutions,

Φt + Φx = 0, Φ(x, 0) = Φ0 (x), with Φ0 (x) = f (x − 0.5) where

⎧ √ 2 ⎪ 2 cos( 3π2x ) − 3, −1  x < − 13 , ⎪ ⎪ ⎪  √ ⎨ 3 + 3 cos(2π x), − 13  x < 0, 3 9 2π (x + 1) + 215 f (x) = − + + ⎪ 2 2 3 − 3 cos(2π x), 0  x < 13 , ⎪ 2 ⎪ ⎪ ⎩ 28+4π +cos(3π x) + 6π x(x − 1), 13  x < 1. 3

42

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

Fig. 2. Comparison of plots without and with an applied limiter for H–J equation with Riemann initial data, Example 5.2, at discontinuity on [−1, 1], T = 1, N = 81, t = 0.8 , using P 1 polynomials.

Fig. 3. Comparison of plots for H–J equation with Riemann initial data, Example 5.2, at discontinuity on [−1, 1], T = 1, N = 81, using P 1 and P 2 polynomials.

We incorporate periodic boundary conditions on the domain [−1, 1]. This example was tested in [35] and [17]. We run the scheme to times t = 2 and t = 32 for P 1 and P 2 polynomials, see Figs. 4 and 5. Example 5.4. We solve the two-dimensional Burgers’ equation using P 1 and P 2 polynomials

Φt +

(Φx + Φ y + 1)2 2

= 0,

with initial data

Φ(x, y , 0) = − cos(x + y ). The computational domain is [0, 2π ]2 and we impose periodic boundary conditions. The solution is still smooth at time t = 0π.52 and we test the order of accuracy at this time in Tables 3 and 4. For P 1 and P 2 polynomials, we get the expected order in all norms. Example 5.5.

Φt +

 Φx2 + Φ 2y + 1 = 0,

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

43

Fig. 4. Plots for Example 5.3 at times t = 2 and t = 32 with N = 201 using P 1 polynomials.

Fig. 5. Plots for Example 5.3 at times t = 2 and t = 32 with N = 201 using P 2 polynomials.

Table 3 L 1 and L ∞ comparison for Example 5.4 at time t = 0.5/π 2 using P 1 polynomials. Nx − 1 = N y − 1 20 40 80 160 320 640

Scheme

L 1 error

L 1 order

L ∞ error

L ∞ order

AEDG AEDG AEDG AEDG AEDG AEDG

0.895488976035 0.213937258115 0.047592740748 0.010144117991 0.002444085288 0.000523995236

2.06548792 2.16837431 2.23009811 2.05327690 2.22166903

0.217832830215 0.067162140900 0.016019009188 0.003553939842 0.000778890123 0.000152697730

1.697501279 2.067863303 2.172293751 2.189927533 2.350741203

Table 4 L 1 and L ∞ comparison for Example 5.4 at time t = 0.5/π 2 using P 2 polynomials Nx − 1 = N y − 1 20 40 80 160 320

Scheme

L 1 error

L 1 order

L ∞ error

L ∞ order

AEDG AEDG AEDG AEDG AEDG

0.115413532776 0.014623191059 0.002217643937 0.000459523717 0.000062313598

2.980482322 2.721158518 2.270816516 2.882520402

0.02908625138 0.00295732018 0.00030991542 0.00005611974 0.00000824372

3.297975021 3.254344004 2.465294309 2.767140916

44

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

Fig. 6. Surface and Contour plots for Example 5.5 at time t = 0.5 with N x = N y = 81 using P 1 polynomials. Table 5 L 1 and L ∞ comparison for Example 5.6 at time t = 0π.25 using P 1 polynomials. Nx − 1 = N y − 1 20 40 80 160 320 640

Scheme

L 1 error

L 1 order

L ∞ error

L ∞ order

AEDG AEDG AEDG AEDG AEDG AEDG

0.877920684785 0.206884240780 0.046040667757 0.010082112689 0.002482044498 0.000537769968

2.085266853 2.167843090 2.191110773 2.022197099 2.206467886

0.273052399152 0.091023173978 0.026474956541 0.005777262280 0.001550191141 0.000355290164

1.584872034 1.781605578 2.196170418 1.897939880 2.125376458

Table 6 L 1 and L ∞ comparison for Example 5.6 at time t = 0π.25 using P 2 polynomials. Nx − 1 = N y − 1 20 40 80 160 320

Scheme

L 1 error

L 1 order

L ∞ error

L ∞ order

AEDG AEDG AEDG AEDG AEDG

0.047529173052 0.006912332179 0.000923764651 0.000156130943 0.000022947196

2.781568846 2.903575303 2.564768851 2.766366722

0.017053409963 0.001645954205 0.000192635173 0.000034837194 0.000005077414

3.373064144 3.094981149 2.467170832 2.778462664

with smooth initial data

Φ(x, y , 0) =

1 4





cos(2π x) − 1 cos(2π y ) − 1 − 1,

and computational domain [0, 1]2 and periodic boundary conditions. The results at time t = 0.5 are shown in Fig. 6. The AEDG scheme provides high resolution in the formation of the singularity. Example 5.6. We now consider an example with a non-convex Hamiltonian:

Φt − cos(Φx + Φ y + 1) = 0, with initial data

Φ(x, y , 0) = − cos



π (x + y ) 2

 ,

on the computation domain [−2, 2]2 with periodic boundary conditions. At time t = 0π.52 , the solution is still smooth and we

test the order of accuracy using P 1 and P 2 polynomials. The results in Tables 5 and 6 show that we get the desired order in all norms. At time t = π12 , the solution develops discontinuous derivative and results at time t = 1π.52 are shown in Fig. 7.

Example 5.7. We now consider an example relating to controlling optimal cost determination from [25,16] with a nonsmooth Hamiltonian:

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

Fig. 7. Surface plot for Example 5.6 at time t = 1π.25 using P 2 polynomials, with uniform mesh x =  y =

45

1 . 20

Fig. 8. Numerical solution and optimal control plots for Example 5.7 at time t = 1 with N x = N y = 81 using P 2 polynomials.

  1 Φt + sin( y )Φx + sin(x) + sign(Φ y ) Φ y − sin2 ( y ) + cos(x) − 1 = 0, 2

with initial data

Φ(x, y , 0) = 0, on the computation domain [−π , π ]2 with periodic boundary conditions. The results for the numerical solution and optical control sign(φ y ) are shown in Fig. 8 using P 2 polynomials and are in agreement with those found in [25,16]. 6. Concluding remarks In this work, we have developed a high-order AEDG scheme to solve the Hamilton–Jacobi equations by following the previous work on high resolution AE finite difference schemes for the Hamilton–Jacobi equation. The AEDG methods have the advantage of easy formulation and implementation, and efficient computation of the solution. High-order accuracy is achieved through the high-order non-oscillatory polynomial approximation and a matching time discretization, numerical experiments in one and two dimensions for P 1 and P 2 polynomials are presented. The L 2 stability of the method was proven for the linear equation.

46

H. Liu, M. Pollack / Journal of Computational Physics 258 (2014) 31–46

Acknowledgements This research was supported by the National Science Foundation under grant DMS13-12636 and the KI-Net research network. References [1] Ahmed Haseena, High-resolution alternating evolution schemes for hyperbolic conservation laws and Hamilton–Jacobi equations, PhD thesis, Department of Mathematics, Iowa State University, 2008. [2] S. Bryson, A. Kurganov, D. Levy, G. Petrova, Semi-discrete central-upwind schemes with reduced dissipation for Hamilton–Jacobi equations, IMA J. Numer. Anal. 25 (1) (2005) 113–138. [3] S. Bryson, D. Levy, High-order central WENO schemes for multidimensional Hamilton–Jacobi equations, SIAM J. Numer. Anal. 41 (4) (2003) 1339–1369. [4] S. Bryson, D. Levy, High-order semi-discrete central-upwind schemes for multi-dimensional Hamilton–Jacobi equations, J. Comput. Phys. 189 (1) (2003) 63–87. [5] Y. Cheng, C.-W. Shu, A discontinuous Galerkin finite element method for directly solving the Hamilton–Jacobi equations, J. Comput. Phys. 223 (1) (2007) 398–415. [6] B. Cockburn, S. Hou, C.-W. Shu, The Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV. The multidimensional case, Math. Comput. 54 (190) (1990) 545–581. [7] B. Cockburn, S.Y. Lin, C.-W. Shu, TVB Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws. III. Onedimensional systems, J. Comput. Phys. 84 (1) (1989) 90–113. [8] B. Cockburn, C.-W. Shu, TVB Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework, Math. Comput. 52 (186) (1989) 411–435. [9] B. Cockburn, C.-W. Shu, The Runge–Kutta local projection P 1 discontinuous Galerkin finite element method for scalar conservation laws, RAIRO Model. Math. Anal. Numer. 25 (3) (1991) 337–361. [10] B. Cockburn, C.-W. Shu, The Runge–Kutta discontinuous Galerkin method for conservation laws. V. Multidimensional systems, J. Comput. Phys. 141 (2) (1998) 199–224. [11] M.G. Crandall, L.C. Evans, P.L. Lions, Some properties of viscosity solutions of Hamilton–Jacobi equations, Trans. Am. Math. Soc. 282 (1984) 487–502. [12] M.G. Crandall, P.-L. Lions, Viscosity solutions of Hamilton–Jacobi equations, Trans. Am. Math. Soc. 277 (1) (1983) 1–42. [13] M.G. Crandall, P.-L. Lions, Two approximations of solutions of Hamilton–Jacobi equations, Math. Comput. 43 (167) (1984) 1–19. [14] S. Gottlieb, C.-W. Shu, E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev. 43 (1) (2001) 89–112. [15] A. Harten, B. Engquist, S. Osher, S.R. Chakravarthy, Uniformly high-order accurate essentially nonoscillatory schemes. III, J. Comput. Phys. 71 (2) (1987) 231–303. [16] C. Hu, C.-W. Shu, A discontinuous Galerkin finite element method for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (2) (1999) 666–690. [17] G.-S. Jiang, D. Peng, Weighted ENO schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (6) (2000) 2126–2143. [18] L. Krivodonova, J. Xin, J.F. Remacle, N. Chevaugeon, J.E. Flaherty, Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws, Appl. Numer. Math. 48 (2004) 323–338. [19] A. Kurganov, S. Noelle, G. Petrova, Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton–Jacobi equations, SIAM J. Sci. Comput. 23 (3) (2001) 707–740. [20] A. Kurganov, E. Tadmor, New high-resolution semi-discrete central schemes for Hamilton–Jacobi equations, J. Comput. Phys. 160 (2) (2000) 720–742. [21] D. Levy, S. Nayak, C.-W. Shu, Y.-T. Zhang, Central WENO schemes for Hamilton–Jacobi equations on triangular meshes, SIAM J. Sci. Comput. 28 (6) (2006) 2229–2247. [22] C.-T. Lin, E. Tadmor, High-resolution nonoscillatory central schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (6) (2000) 2163–2186. [23] C.-T. Lin, E. Tadmor, L 1 -stability and error estimates for approximate Hamilton–Jacobi solutions, Numer. Math. 87 (2001) 701–735. [24] H. Liu, An alternating evolution approximation to systems of hyperbolic conservation laws, J. Hyperbolic Differ. Equ. 5 (2) (2008) 1–27. [25] F.Y. Li, S. Yakovlev, A central discontinuous Galerkin method for Hamilton–Jacobi equations, J. Sci. Comput. 45 (1–3) (2010) 404–428. [26] Y. Liu, Central schemes on overlapping cells, J. Comput. Phys. 209 (1) (2005) 82–104. [27] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1) (1988) 12–49. [28] S. Osher, C.-W. Shu, High-order essentially nonoscillatory schemes for Hamilton–Jacobi equations, SIAM J. Numer. Anal. 28 (4) (1991) 907–922. [29] W.H. Reed, T.R. Hill, Triangular mesh methods for the neutron transport equation, Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973. [30] H. Saran, H. Liu, Formulation and analysis of alternating evolution schemes for conservation laws, SIAM J. Sci. Comput. 33 (2011) 3210–3240. [31] H. Liu, M. Pollack, H. Saran, Alternating evolution schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 35 (1) (2013) 122–149. [32] C.-W. Shu, High order ENO and WENO schemes for computational fluid dynamics, in: High-Order Methods for Computational Physics, in: Lect. Notes Comput. Sci. Eng., vol. 9, Springer, Berlin, 1999, pp. 439–582. [33] J. Qiu, C.W. Shu, Hermite WENO schemes for Hamilton–Jacobi equations, J. Comput. Phys. 204 (2005) 82–99. [34] P.E. Souganidis, Approximation schemes for viscosity solutions of Hamilton–Jacobi equations, J. Differ. Equ. 59 (1985) 1–43. [35] J. Yan, S. Osher, A new discontinuous Galerkin method for Hamilton–Jacobi equations, J. Comput. Phys. 230 (1) (2011) 232–244. [36] Y.-T. Zhang, C.-W. Shu, High order WENO schemes for Hamilton–Jacobi equations on triangular meshes, SIAM J. Sci. Comput. 24 (2003) 1005–1030.