Compressed Semi-Discrete Central-Upwind Schemes for Hamilton- Jacobi Equations Steve Bryson*, Alexander Kurganovt, Doron Levy: and Guergana Petrovas Abstract We introduce a new family of Godunov-type semi-discrete central schemes for multidimensional Hamilton-Jacobi equations. These schemes are a less dissipative generalization of the central-upwind schemes that have been recently proposed in series of works. We provide the details of the new family of methods in one, two, and three space dimensions, and then verify their expected low-dissipative property in a variety of examples. AMS subject classification: Primary 65M06; Secondary 35L65 Key Words: Hamilton-Jacobi equations, central-upwind schemes, semi-discrete methods. 1 Introduction We consider the multidimensional Hamilton-Jacobi equation, + pt H(V,p) = 0, x E Rd, with Hamiltonian H. First-order numerical schemes that converge to the viscosity solution of (1.1) were first introduced by Crandall and Lions in [5] and by Souganidis in [17]. Recent attempts to obtain higher-order approximate solutions of (1.1) include upwind methods, discontinuous Galerkin methods, and others. Here, we study a class of projection-evolution methods, called Godunov-type schemes. The main structure of these schemes is as follows: one starts with the point values of the solution, constructs an (essentially) non-oscillatory continuous piecewise polynomial interpolant, and then evolves it to the next time level while projecting the solution back onto the computational grid. The key idea in Godunov-type central schemes is to avoid solving (generalized) Riemann problems, by evolving (locally) smooth parts-of the solution. Second-order staggered Godunov-type central schemes were introduced by Lin and Tadmor in [14, 151. L'-convergence results for these schemes were obtained in [14]. More efficient non-staggered central schemes as well as genuinely multidimensional generalizations of the schemes in [15] were presented in [l], with high-order extensions (up to fifth-order) proposed in [2, 31. 'Program in Scientific Cornputing/Computazonal Mathematics, Stanford University and the NASA Advanced Super- computing Division, NASA Ames Research Center, Moffett Field, CA 94035-1000; [email protected] +Department of Mathematics, Tulane University, 6823 St. Charles Avenue, New Orleans, LA 70115; [email protected], tel: +1-504-862-3443, fax: +1-504-865-5063 fDepartment of Mathematics, Stanford University, Stanford, CA 94305-2125; dlevyamath. stanf ord. edu, tel: +1-650-723-4157, fax: +1-650-725-4066 §Department of Mathematics, Texas A&M University, College Station, TX 77843-3368; gpetrovaQmath.t am. edu, tel: +1-979-845-5298, fax: +1-979-845-6028 2 S. BRYSONA, . KURGANODV,. LEVY & G. PETROVA Second-order semi-discrete Godunov-type central schemes were introduced in [12], where local speeds of propagation were employed to reduce the numerical dissipation. The numerical viscosity was fur- ther reduced in the central-upwind schemes [lo] by utilizing one-sided estimates of the local speeds of propagation. Higher-order extensions of these schemes were introduced in [4], where weighted es- sentially non-oscillatory (WENO) interpolants were used to increase accuracy. WENO interpolants were originally developed for numerical methods for hyperbolic conservation laws [16, 81, and were first implemented in the context of upwind schemes for Hamilton-Jacobi equations in [7]. Godunov-type central-upwind schemes are constructed in two steps. First, the solution is evolved to the next time level on a nonuniform grid (the location of the grid points depends on the local speeds, and thus can vary at every time step). The solution is then projected back onto the original grid. The projection step requires an additional piecewise polynomial reconstruction over the nonuniform grid. In this paper we show that in the semi-discrete setting different choices of such a reconstruction lead to different numerical Hamiltonians, and thus to different schemes. In particular, we can recover the scheme from [lo]. A more careful selection of the reconstruction results in a new central-upwind scheme with smaller numerical dissipation. This approach was originally proposed in [ll],w here it was applied to one-dimensional (1-D) systems of hyperbolic conservation laws. It has been recently generalized and implemented for multidimensional systems of hyperbolic conservation laws in [9]. The paper is organized as follows. In $2, we develop new semi-discrete centraI-upwind schemes for 1-D Hamilton-Jacobi equations. We also review the interpolants that are required to complete the construction of the second- and fifth-order schemes. Generalizations to more than one space dimensions (with special emphasis on the two-dimensional setup) are then presented in $3, where the corresponding multidimensional interpolants are also discussed. In $4, we evaluate the performance of the new schemes with a series of numerical tests. Finally, in the Appendix, we prove the monotonicity of the new numerical Hamiltonian. 2 One-Dimensional Schemes 2.1 Semi-Discrete Central-Upwind Schemes for Hamilton-Jacobi Equations In this section, we describe the derivation of a new family of semi-discrete central-upwind schemes for the 1-D Hamilton-Jacobi equation, + Pt H (%> = 0, 2 E R, (2.1) subject to the initial data cp(x,t = 0) = (po(x). We follow the approach in [lo] (see also [12]). For simplicity we assume a uniform grid in space and time with grid spacing Ax and At, respectively. The grid points are denoted by xJ := jAx, t" := nAt, and the approximate value of cp (xJt," ) is denoted by 97. Assume that the approximate solution at time t", (p3n, is given, and that a continuous piecewise- polynomial interpolant @(xt, ") is reconstructed from $'. At every grid point, the maximal right and left speeds of propagation, a; and aJ-, are then estimated by I a: = max {"(4,0> > a; min min{vF ,v$}<uSmax{vp; ,vf} = min{vp;,vf)lu<max{vp;,vf) where cpz are the one-sided derivatives at x = xJ,t hat is * cp : := &(XJ 0,t"). COMPRESSEDC ENTRAL-UPWINSDC EEMES 3 y+n+l Figure 2.1: Central-upwind differencing: 1-D If the Hamiltonian is convex, (2.2) reduces to a: =max{H'(p,),H'(cp,+),Oa}; ,= Imin{H'(p,),H'(p,+),O}). (2.3) We then proceed by evolving the reconstruction (p at the evolution points xyf := xj fa FAt, to the next time level according to (2.1). The time step At is chosen so that xy+ < x;+~)- for all j. Therefore the solution remains smooth at xy* for t E [t",t" +'] (see Figure 2.1) and we can compute the values of the evolved solution {pyzl} by the Taylor expansion: + p?zl = @(xCj",t," ) - AtH (CpZ(xy*, t")) 0 (At') . (2.4) Using the values {p?:'} on the nonuniform grid {xy*}, we construct a new quadratic interpolant .?+I, &(z, tn+') on the interval [xy-, (&,)yc1 pzz(q7 q + where is yet to be determined and is an approximation to tnfl), = (x?+ zjn-)/2. The projection back onto the original grid is then carried out by evaluating &(x, t"+l) at xj, q Note that if the Riemann fan is symmetric, that is, if a: = a;, then = xj. Substituting (2.4) in (2.6) yields Using the Taylor expansion, + @(xykt7" ) = p? f Ata,f pZ* o(At)2, 4 S. BRYSONA, . KURGANOVD,. LEVY& G. PETROVA we arrive at 1 &5)y+1 + - - ( uj' a; (At)' 0 (A t)2. (2.9) 2 We then let At -+ 0, and end up with the (family of) semi-discrete central-upwind schemes: (2.10) & Here, the one-sided speeds of propagation, af, are given by (2.2), and are the left and right derivatives at the point 2 = zj of the reconstruction p(.,t)a t time t. Finally, in order to complete the construction of the scheme, we must determine (@,,);+. ' For (@,z)y+l example, selecting to be independent of At gives - and then (2.10) recovers the centraI-upwind scheme in [lo]. However, since the interpolant $(., tn+') (sZz)7'l is defined on the intervals [xjn_,2 7+],w hose size is proportional to At, it is natural to choose to be proportional to l/At. In this case, the approximation of the second derivative in (2.10) will add a non-zero contribution to the limit as At -+ 0. At the same time, to guarantee a non-oscillatory reconstruction, we should use a nonlinear limiter. For example, one can use the minmod limiter: & , where is the derivative of (2.5), &(z?*, tn+')a re the values of the derivative of the evolved recon- struction ?(., t") at t = tnfl,a nd the multivariate minmod function is defined by min{zj}, if xj > 0 'dj, minmod(zl,x2 ,...) := if xj < 0 Vj, (2.12) otherwise. A different choice of limiter will result in a different scheme from the same family of central-upwind schemes. All that remains is to determine the quantities used in (2.11). Since all data are smooth along the line segments (xY*, t),t n I t < tn+l,w e can use a Taylor expansion to obtain @z(~y*, + @,(z;*, tn+l)= t") O(At). (2.13) - &(q, According to (2.5), the derivative tn+l)o f the new reconstruction $ at time level t*+l is (2.14) and after substituting (2.4) and (2.8) into (2.14), we obtain (2.15) COMPRESSECDE NTRAL-UPWINSDC HEMES 5 Passing to the semi-discrete limit (At + 0) in (2.11), (2.13), and (2.15) gives $Ft,$ Ft At(FXx);+']= +2 minmod (cp: - - c)p; , (2.16) At-0 (a+i a3: ) where (2.17) Finally, substituting (2.16) into (2.10), we obtain the 1-D low-dissipative semi-discrete central-upwind scheme: $int where is given by (2.17). For future reference, we denote the RHS of (2.18) by --ETBKLP. G(., Notice, that in the fully-discrete setting the use of the intermediate quadratic reconstruction tn") at level tn+l (as opposed to the intermediate piecewise linear reconstruction in [lo]) increases the + + accuracy of the resulting fully-discrete scheme: O((At)2 (Ax)') versus O(At (Az)'),w here T is the (formal) order of accuracy of the continuous piecewise polynomial reconstruction @(.,t "). When we pass to the semi-discrete limit (At 0), both quadratic and linear interpolation errors go to 0, and therefore --f the (formal) order of accuracy of both (2.17)-(2.18) and the semi-discrete scheme in [lo] is AS)'), and the temporal error is determined solely by the (formal) order of accuracy of the ODE solver used to integrate (2.18). However, the minmod limiter introduces a new term that leads to a reduction of the numerical dissipation without affecting the accuracy of the scheme. To demonstrate this, we show $pt that is always in the interval [min{(pz, pi}, max{cp$, cp;}], and therefore the absolute value of the term (cp: - c)p; in the numerical dissipation in the scheme from [lo] is always greater than the absolute value of the new term, that is 19 : - cp,l 2 /cp$ - cp ; - minmod (p: - - cp;)l. Indeed, we have where 5 E (min{cp:, cp },; max{cp,f, cp;}). It follows from the definition of the local speeds (2.2) that + aj' - ~ ' (20 0, H'(J) a; 2 0. $int Thus, (2.19) is a convex combination of cp ; and (pi,a nd therefore E [min{cp$, cp ,}, max{cp$, p;}] It was shown in [4]t hat the numerical Hamiltonian HKNPf rom [lo] is monotone, provided that the Hamiltonian H is convex. Here, we state a theorem about the monotonicity of HBKLP- the new, less dissipative Hamiltonian in (2.18). The proof is left to the Appendix. We will consider only Hamiltonians for which H' changes sign, because otherwise either a- -= 0 or a+ E 0 and the Hamiltonian in (2.18) reduces to the upwind one for which such a theorem is known. Theorem 2.1 Let the Hamiltonian H E C2 be convex and satisfy the following two assumptions: (Al) The function + G(u,v ) := 2HN(u)[ (u- v)H'(v)- (H(u)- H(v))] [H'(u)- H'(v)125 0 (2.20) 6 S. BRYSONA, . KURGANOVD,. LEVY& G. PETROVA for all u and w in the set S-(u,w ) u S+(u,w ), where + H(u) - H(w) H’(u) H’(w) I u-v 2 (2.21) and u* is the only point such that H‘(u*) = 0; (A21 For any w and for an arbitrary interval [u,b ], the sets S-(u, w) n [a,b ] and S+(u,v ) n [a,b] are either the empty set or finite unions of closed intervals and/or points. Then the numerical Hamiltonian in (2.18): HBKLP + a-H(u+) + a+H(u-) (u ,u-) := + a+ a- [ u+ -u- u+ - uint Uint - - a+a- + - minmod (2.22) ’ a+ a- a++a- where + ,@t ..-_ a+u+ + a-u- - H(u+)+- H(u-) (a+ a-) (a+ a-) ’ and a+ := a+(u+, u-)= max{H’(u+), H’(u-), 0}, a- := a-(u+, u-) = 1 min{H’(u+), H’(u-), O}( is monotone, that is, HBKLPis a non-increasing function of u+ and a non-decreasing function of u-. Remarks. 1. The classification of all Hamiltonians that satisfy conditions (2.20)-(2.21) is an open problem. However, we were able to verify these conditions for certain Hamiltonians. For example, a + + straightforward computation shows that for any convex quadratic Hamiltonians H(u)= au2 bu c, = the function G(u,v ) 0, the sets in (A2) are either 0 or one closed interval, or one point, and therefore the theorem holds. Another example, for which Theorem 2.1 is valid, is H(u) = u4. In this case, the sets (2.21) are + + S-(u, v) = {(u,v) : u w 2 0 2 v} , S+(U, v) = {(u,v ) : u v 5 0 I v) , and, as one can easily verify, the function + + + G(u,w ) = -8(u - I J ) ~ ( u3u~2w 6uw2 2w3) 5 0, in S-(u, w) U S+(u,v ). As for the sets in assumption (A2), they are either 0 or one closed interval, or one point. 2. Notice that assumption (A2)i n Theorem 2.1 is needed onIy for technical purposes and in fact it is satisfied by (almost) every Hamiltonian H that arises in applications. 2.2 A Second-Order Scheme A non-oscillatory second-order scheme can be obtained if one uses a non-oscillatory continuous piecewise quadratic interpolant @. The values of the one-sided derivatives of @ at (xj,t”) in (2.17) and (2.18), are given by (2.23) COMPRESSECDE NTRAL-UPWINSDC HEMES 7 where the second derivative is computed with a nonlinear limiter. For example, Here, 8 E [l,21 and the minmod function is given by (2.12). The scheme requires a second-order ODE solver. 2.3 Higher-Order Schemes In this section, we briefly describe the third- and fifth-order weighted essentially non-oscillatory (WENO) reconstructions. They were derived in [4] in the context of central-upwind schemes, and are similar to those used in high-order upwind schemes [7]. In smooth regions, the WENO reconstructions use a convex combination of multiple overlapping reconstructions to attain high order accuracy. In nonsmooth regions, a smoothness measure is em- ployed to increase the weight of the least oscillatory reconstruction. Here, we reconstruct the one-sided derivatives (~p:)~,~a t x = xj for k = 1,.. . , d stencils, and write the convex combination: d d * 9: = w&(’P~) k,j> w& = ‘7 Wfk, j 2 O, (2.25) k=l k=l where the values cp$ are to be used in the scheme (2.17)-(2.18). The weights w $ ~ar e defined as (2.26) The constants ckf are set so that the convex combination in (2.25) is of the maximal possible order of accuracy in smooth regions. We take p = 2 and choose E = to prevent the denominator in (2.26) from vanishing. A third-order WENO reconstruction is obtained in the case d = 2 with The constants cf are given by c+- -- 2 +- - - _1 1 -c2 --> 3 c, -c1 -3’ and the smoothness measures are Here, 8 S. BRYSONA, . KURGANOVD,. LEVY& G. PETROVA A fifth-order WEN0 reconstruction is obtained when d = 3. In this case, The constants ct are given by c+---- 3 C + - - = - 1 c*-- 3 1 -c3 - 10' 3 -c1 10' 2 - 5' and the smoothness measures are S1-:j t-- sj L-2, 01 , S2j = sj [-1,1J , sZ~= Sj IO, 21 , SU = Sj [-3, -11 S<j = Sj [-2,O] , S<j = Sj [-II 11 . ) The time evolution of (2.18) should be performed with an ODE solver whose order of accuracy is compatible with the spatial order of the scheme. In our numerical exampIes, we use the strong stability preserving (SSP) Runge-Kutta methods from [6]. 3 Multidimensional Schemes In this section, we derive the two-dimensional (2-D) generalization of the compressed semi-discrete central-upwind scheme (2.17)-(2.18), and then extend it to three space dimensions. We also comment on the multidimensional interpolants that these extensions require. 3.1 A Two-Dimensional Scheme We 'consider the 2-D Hamilton- Jacobi equation, and proceed as in [lo]. we assume that at time t = t" the approximate point Values Pj"k M 'p(xc3y,k , t") are given, and constreuct a 2 -D continuous piecewise-quadratic interpolant, p(z,y , t"), defined on the + cells sjk := {(x, y) : 5 1). On each cell sjk there will be four such interpolants (labeled NW, NE, SE, and SW), one for each triangle that constitutes sjk (see Figure 3.1). Specific examples of @(x,y , t") are discussed in $3.3. Similarly to the 1-D case, we use the maximal values of the one-sided local speeds of propagation in the 2- and y-directions to estimate the widths of the local Riemann fans. These values at any grid point (qy,k) are given by where cjk := [zi-;, xi+$] x [yk-;, yk+;], ( . )+ := ma(.,0 ), ( .)- := mi.(., 0)) and (&,flu)* is the gradient of H. COMPRESSECDE NTRAL-UPWINSDC HEMES 9 Li, k+l) Figure 3.1: Central-upwind differencing: 2-D The reconstruction Cp(z, y,t ") is then evolved according to the Hamilton-Jacobi equation (3.1). Due to the finite speed of propagation, for sufficiently small At, the solution of (3.1) with initial data Cp is smooth around (zyi, y&) where xy* := Xj f a:,At, y;* := Yk f bTkAt, see Figure 3.1. We denote by e*,k*:= Cp(zy*, y;?,t "), and use the Taylor expansion to calculate the intermediate values at the next time level t = tn+l $$,;* = q*,kf - At. W&kq*, y;*, tn>,C p&y*,Y;*, t"))+ Ww2. (3.3) We now project the intermediate values qYz,;* onto the original grid points (zj,Y k). First, similarly to (2.5), we use new 1-D quadratic interpolants in the y-direction, $(xy*, .,tn+'),t o obtain (3.4) - where (Gyy);$ M 'pyy(xy*$, , tn+') and $? := (y:+ + y;-)/2. Next, we use the values '$(zy*, yk, tn+l) to construct another 1-D quadratic interpolant &(., Yk,t nfl)t,h is time in the z-direction, whose values at the original grid points are pjnk+ l := $(zj,y k, tn+') (3.5) (@zz):il + Here, (pzzy(kq, tn,+ ') and := (xy+ zy-)/2. We choose '&(:;). to be the weighted average 10 S. BRYSONA, . KURGANOVD, . LEVY & G. PETROVA We then substitute (3.4) and (3.6) into (3.5), and obtain Substituting (3.3) into (3.7) yields The values q*,kf are computed by the Taylor expansions: * * * * + (3.9) qf,k& = (Pxk Atajk'P+ Atbjk'Py o(At)2> where := (PZ(xj~0Yk, ,t ") and := FY(zj,Y kzto,t ") are the corresponding right and left derivatives of the continuous piecewise quadratic reconstruction at (zj,Y k). Next, substituting (3.9) into (3.8) gives (3.10)