- Email: [email protected]
A local piecewise parabolic method for Hamilton–Jacobi equations
Applied Numerical Mathematics 45 (2003) 499–511 www.elsevier.com/locate/apnum
A local piecewise parabolic method for Hamilton–Jacobi equations Youssef Stiriba 1 CERFACS/CFD Team, 42 avenue Gaspard Coriolis, 31057 Toulouse cedex 1, France
Abstract In this paper we present, analyze and implement a new local piecewise parabolic method for nonlinear Hamilton– Jacobi equations. The scheme is third-order accurate in smooth regions, uses a concept of local smoothing to prevent the excessive increase of the total variation at discontinuity, and has a local stencil in the sense that it does not extrapolate from data of the smoothest neighboring cells. One and two-dimensional numerical experiments, accuracy tests, and the behavior of the total variation of the approximate solution are presented to prove the accuracy and good resolution of our method. 2002 IMACS. Published by Elsevier Science B.V. All rights reserved.
1. Introduction We consider numerical solutions of Hamilton–Jacobi (H–J) equations φt + H (x1 , . . . , xd , t, φx1 , . . . , φxd ) = 0, φ(x, 0) = φ 0 (x),
(1)
by a third-order local parabolic method. Solutions of H–J equations are continuous, but may have discontinuous derivatives even with smooth initial condition φ 0 (x), and these solutions are in general not unique. Therefore, it is necessary to introduce the notion of entropy condition to select a unique relevant solution [1]. Recently a wide variety of finite difference schemes have been developed. For instance, in [2] the authors introduced a class of monotone schemes and proved that these schemes converge to the viscosity solution of (1). However, such methods are at most first-order accurate and too dissipative for problems of practical interest. In [9,10] the authors have used the fact that the 1D Hamilton–Jacobi equation (1) is closely related to the 1D conservation law ut + f (u)x = 0, to adapt the high-order essentially E-mail address: [email protected] (Y. Stiriba). 1 Current address: Abteilung für Numerische Mathematik; IGPM, RWTH Aachen, Templergraben 55, 52062 Aachen,
Germany. 0168-9274/02/$30.00 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. doi:10.1016/S0168-9274(02)00211-8
500
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
nonoscillatory (ENO) schemes [3,11] to solve Hamilton–Jacobi equations in conjunction with several monotone building blocks. Their numerical experiments indicate good accuracy and better resolution of the corner-like discontinuity. The same idea was reported in [4] to approximate solution of (1) by Weighted ENO (WENO) schemes. In [7,5], the authors extended their high-resolution central schemes for hyperbolic conservation laws to the Hamilton–Jacobi equations. In this paper we construct a piecewise parabolic method for approximating Hamilton–Jacobi equations (1), by applying a preprocessing of derivatives in each cell to prevent the excessive increase of the total variation of the solution. We adapt this method in conjunction with Osher and Shu local Lax– Friedrichs flux function [10] to solve H–J equations. The resulting scheme is more local in the sense that it uses only four points to reconstruct piecewise smooth functions in contrast to ENO-3 that depends on six point values, gives sharp resolution of jump derivatives, and is total variation bounded even when approximating discontinuous solutions. The outline of the paper is the following: In Section 2, we present the reconstruction procedure. In Section 3, we describe the numerical scheme for H–J equations in two space dimension. Numerical results are presented in Section 4. Finally, we close the paper with some concluding remarks.
2. The local piecewise parabolic reconstruction Before we construct a local piecewise parabolic (LPP) scheme for approximating H–J equations, we first present the LPP reconstruction procedure. Let f (x) be a piecewise smooth function. We consider a uniform grid in the x − t plane, with mesh spacing x and time step t. We denote the grid points by xi = i x and tn = t, and the computational cells by Ci = [xi−1/2 , xi+1/2 ] where xi±1/2 = (i ± 1/2) x, ∀i ∈ Z. We will also use: (1) The cell averaged of f : 1 f (ξ ) dξ, vi = x
i ∈ Z.
(2)
Ci
(2) The forward/backward differences ±(vi±1 − vi ) , i ∈ Z. di±1/2 = x The reconstruction problem is to find a parabola ci Pi (x) = ai + bi · (x − xi ) + · (x − xi )2 , 2 defined on Ci , such that it is a third-order approximation to f (x) on Ci : f (x) = Pi (x) + O x 3 every time f (x) is smooth in Ci . The parabola Pi (x) is required to satisfy 1 (1) vi = x Ci Pi (ξ ) dξ , (2) di+1/2 = Pi (xi+1/2 ),
(3)
(4)
(5)
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
501
for every i. These imply that the reconstruction is third-order accurate. For each cell, there exists a unique parabola of the form (4) determined from vi , di−1/2 , and di+1/2 by di+1/2 − di−1/2 di+1/2 + di−1/2 , bi = , x 2 and satisfies di+1/2 + di−1/2 . P (xi ) = 2 The lateral derivatives of Pi are ci =
P (xi+1/2 ) = bi + ci · x
and
and
ai = vi − ci
x 2 , 24
P (xi−1/2 ) = bi − ci · x.
(6)
(7)
(8)
Since the reconstruction is confined in each computational cell and repeated every time step, the local total variation (LTV), defined by LTV i = TV(Pi ) (the total variation of the function P in the cell Ci , see [6]), must be controlled. A method is called local total variation bounded (LTVB) if LTV(Pi ) M · x
∀i,
for some fixed positive M that depends only on the initial conditions. The LTV for our reconstruction is x · |d | + |d | if di+1/2 · di−1/2 0, i+1/2 i−1/2 2 LTV(Pi ) = x 2·|di+1/2 |·|di−1/2 | · |di+1/2 | + |di−1/2 | − |di+1/2 |+|di−1/2 | elsewhere, 2
(9)
(10)
for every i. It is clear that the present parabola is not LTVB. So, the LTV may blow up near discontinuities (i.e., for some i such that di+1/2 = O( x −1 )). This behavior has been observed when approximating Burgers’ equation (when the waves break into discontinuities). The total variation of the function f (x) on a cell Ci where it is smooth depends essentially, up to O( x 3 ), on the size of di±1/2 (this follows from trapezoidal rule). Thus to prevent the increase of LTV, we shall correct the lateral derivatives di+1/2 and di−1/2 by a pair of local consistent numerical derivatives {dli , dri } such that (11) dli = di−1/2 + O x 2 and dri = di+1/2 + O x 2 , which preserve third order accuracy. Below we outline the LPP reconstruction. To construct the parabola in a nontransition cell Ci (i.e., when di+1/2 · di−1/2 > 0), we can use (6) and assign to b the average di+1/2 + di−1/2 . (12) 2 The resulting parabola has as derivatives at the cell interfaces di+1/2 and di−1/2 . However, with this choice the LTV may increase. So, for transition cells Ci with di+1/2 = O( x −1 ), the LTV is (13) LTV(Pi ) = O x −1 . b=
Thus, the method becomes unstable, and more smoothing is necessary to prevent the excessive increase of the local total variation of the solution, as observed in our numerical experiments. To provide more smoothing, we can use the harmonic mean of the lateral derivatives 2 · di+1/2 · di−1/2 , (14) b= di+1/2 + di−1/2
502
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
as in [8], used to construct local hyperbola, such that the derivatives of the parabola interpolates b at the cell center xi and the lateral grid derivative with smallest absolute value. This gives a parabola with lateral derivative smaller than the original one in nonsmooth side. The harmonic mean is smaller in absolute value than (12), is a second-order approximation of the derivative at the cell center xi , and gives more smoothing and good results according to our numerical experiments (see also Section 4 and [12]). It is possible to choose other means such that the method of reconstruction is LTVB, but we have found, through our numerical experiments, that our choice gives good results. In transition cells (i.e., when di+1/2 ·di−1/2 < 0), we follow the ENO idea [11] and change the derivative with largest absolute value by the other one multiplied by x 2 . The algorithm in the flavor of FORTRAN 77 is the following Algorithm1 (LPP reconstruction). *tol = O x 2 if (|di−1/2 | tol) and (|di+1/2 | tol) then bi = 0 and ci = 0 else if (|di−1/2 | tol) or (Ci is a transition cell with |di−1/2 | |di+1/2 |) bi = 2 · di+1/2 · else
x 2 1+ x 2
and ci =
di+1/2 · x 2 −bi x
if (|di+1/2 | tol) or (Ci is a transition cell with |di+1/2 | |di−1/2 |) bi = 2 · di−1/2 · else
2·d
x 2 1+ x 2
and ci =
bi −di−1/2 · x 2 x
·d
i−1/2 i+1/2 bi = di−1/2 +di+1/2 if (|di−1/2 | |di+1/2 |) then b −di−1/2 ci = i 2· x else d −bi ci = i+1/2 2· x
We compute ai from vi and ci using (6). We note that we can also implement our algorithm by using the harmod slope limiter, which consist of the harmonic mean on the lateral slopes when they have the same sign and zero when they have different signs. The algorithm is the following: Algorithm 2 (LPP-HARMOD reconstruction). if (di−1/2 · di+1/2 tol) then bi = 0 and ci = 0 else if (|di−1/2 | tol) and (|di+1/2 | tol) then d +d d −d bi = i−1/2 2 i+1/2 and ci = i+1/2 x i−1/2 else 2·di−1/2 ·di+1/2 bi = di−1/2 +di+1/2 if (|di−1/2 | |di+1/2 |) then b −di−1/2 ci = i 2· x
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
else ci =
503
di+1/2 −bi 2· x
The LPP-HARMOD has a much simpler form, and gives similar results to those obtained with LPP approximation, according to our numerical experiments. Remark 1. The LPP reconstructions is more local in the sense that it uses only four points to reconstruct piecewise smooth functions in contrast with other reconstruction like ENO-3 or WENO-5 that depend on six point values. This has good consequence in reducing numerical noise generated by singularities during numerical process (see Section 4). Remark 2. We note also that our methods does not use second differences to compute the corner like discontinuity (a jump in the first derivative). The following proposition shows that the LPP reconstruction is LTVB. Proposition 1. The LPP method of reconstruction of the function f (x) is LTVB. Proof. Since f (x) is a piecewise smooth function, we can show, as in Theorem 2.1 of [8], that there exists a constant M > 0 that depends on f (x) in smooth regions except for finite number of i for which di+1/2 = O( x −1 ), such that |di+1/2 | M.
(15)
If Ci a nontransition cell and |di−1/2 | |di+1/2 |, then the preprocessed derivatives at cell interfaces are: dli = di−1/2
and
dri = 2bi − di−1/2 .
(16)
Hence from (10) x |di−1/2 | + |2bi − di−1/2 | 3 · M · x. 2 We use the same argument for transition cell. LTV(Pi )
(17)
3. The LPP method for Hamilton–Jacobi equations In this section, we consider the 2D Hamilton–Jacobi equation φt + H (φx , φy ) = 0, (18) φ(x, y, 0) = φ 0 (x, y) where (x, y) ∈ [a, b] × [c, d], t 0, H is a locally Lipschitz continuous Hamiltonian, and the initial data φ 0 (x, y) is a locally Lipschitz continuous. n denote the numerical Let x and y be the grid spacing in x and y directions, respectively. Let φi,j approximation to the viscosity solution φ(i x, j y, n t) of (18). We use the method of line, by discretizing only in space and leaving the problem continuous in time, in designing our algorithm. In such approach the semi-discrete form of the LPR method is + dφi,j − + − φx,i,j = −H , φx,i,j , φy,i,j , φy,i,j . (19) dt
504
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
± Here φx,i,j are LPP approximations to the left and right x-derivative of φ(x, y, t) at (xi , yj , tn ). Similarly ± for φy,i,j . The numerical Hamiltonian is a Lipschitz continuous monotone flux (i.e., nonincreasing in its fifth and seventh arguments and nondecreasing in its sixth and eighth arguments), and consistent with H
(u, u, v, v) = H (u, v). H
(20)
There are many Hamiltonian that may be used as a building block (see [10]). In our numerical experiments we use the local Lax–Friedrichs (LLF) flux + − + − 1 + − + − v 1 + u + v u u ,u ,v ,v = H , − α u+ − u− − β v + − v − , (21) H 2 2 2 2 where
H1 (u, v) , max u∈I (u− ,u+ ),v∈[C,D] I (u, v) = min(u, v), max(u, v) ,
α=
H1 (u, v) =
∂H (u, v) , ∂u
β=
H2 (u, v) =
max
H2 (u, v) ,
(22)
v∈I (v − ,v + ),u∈[A,B]
(23) ∂H (u, v) . ∂v
(24)
Here [A, B] is the value range for u± and [C, D] is the value range for v ± . If we replace I (u− , u+ ) and I (v − , v + ) by [A, B] and [C, D], respectively, we recover the Lax–Friedrichs flux. Time discretization is done by the third-order TVD (total variation diminishing) Runge–Kutta of Shu and Osher [11]. Given the numerical solution at time n · t, {φin }, then the numerical solution at the new time (n + 1) · t, {φin+1 }, is computed by the following recursive formula φi(0) , φi(1) = φi(0) + t H
φi(0) = φin ,
1 2 2 (2) φi , φi(3) = φi(0) + φi(2) + t H 3 3 3
3 1 1 (1) φi , φi(2) = φi(0) + φi(1) + t H 4 4 4 φi(n+1) = φi(3) ,
i denotes the numerical Hamiltonians computed from the point values φi for all i. Here H sponding to substep j , where j = 1, 2, 3. (j )
(j )
corre-
4. Numerical examples In this section we study the performance of our numerical scheme by applying it to some test problems. 4.1. Accuracy tests Example 1. We solve the one-dimensional H–J equations: φt + H (φx ) = 0, −1 x 1, φ(x, 0) = − cos(π x),
(25)
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
505
Table 1 L1 and L∞ errors and orders for 1D Burgers’ equation N
L1 error
L1 order
L∞ error
L∞ order
0.331E–01 0.133E–01 0.261E–02 0.441E–03
– 1.31 2.34 2.56
0.449E–01 0.152E–01 0.281E–02 0.691E–03
– 1.70 2.43 2.18
L∞ error
L∞ order
0.824E–01 0.225E–02 0.514E–03 0.118E–03
– 1.87 2.13 2.12
0.224E–01 0.510E–02 0.127E–02 0.288E–03
– 2.13 2.00 2.14
T = 0.5/π 2 10 20 40 80
0.154E–01 0.265E–02 0.432E–03 0.635E–04
– 2.53 2.61 2.76 T = 1.5/π 2
10 20 40 80
0.216E–01 0.338E–02 0.454E–03 0.441E–04
– 2.67 2.89 3.36
Table 2 L1 and L∞ errors and orders for 1D nonconvex problem N
L1 error
L1 order T = 0.5/π 2
10 20 40 80
0.305E–02 0.878E–03 0.117E–03 0.177E–04
– 1.79 2.90 2.72 T = 1.5/π 2
10 20 40 80
0.868E–02 0.176E–02 0.292E–03 0.441E–04
– 2.30 2.59 2.72
with a convex Hamiltonian H (Burgers’ equation): (u + α)2 , 2 and a nonconvex Hamiltonian H : H (u) =
(26)
H (u) = − cos(u + α).
(27)
For simplicity, we take α = 1. For both cases the solution is smooth at time T = 0.5/π 2 , and develop a discontinuous derivative at time T = 1.5/π 2 . Let φin be the reconstructed solution at (xi , t n ). The discrete L1 and L∞ norms are given by N
n
φ(xi , tn ) − φ n ,
e ·, t 1 = x i L
(28)
i=1
n
e ·, t ∞ = max φ(xi , tn ) − φ n . i L 1iN
(29)
506
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
Fig. 1. One-dimensional Burgers’ equation; (a) before singularity, T = 0.5/π 2 ; (b) after singularity, T = 1.5/π 2 .
Fig. 2. One-dimensional nonconvex problem; (a) before singularity, T = 0.5/π 2 ; (b) after singularity, T = 1.5/π 2 .
Here N is the number of grid points used. In Tables 1 and 2, we list L1 - and L∞ -errors and convergence rates for both convex and nonconvex Hamiltonians. The exact solution can be obtained using the method of characteristic (see [9]). We observe L1 -third-order-accuracy and more than second-order accuracy in L∞ -norm. Numerical results together with exact solutions are displayed in Figs. 1 and 2. We see that the solution is free from spurious oscillations, and the scheme has good behavior at the corner-likediscontinuity. 4.2. The behavior of the total variation Example 2. We numerically examine the behavior of the total variation N
n
φ − φ n
TV φ n := i+1 i i=1
(30)
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
507
Fig. 3. Total variation of the solution; N = 40 (‘·’), N = 80 (‘- · -’), N = 160, (‘-’), exact TV (‘- -’); (a) Burgers’ equation; (b) Nonconvex problem.
of the discrete solution. Figs. 3 and 4 display the TV of the solution with N = 40, 80, 160. We observe that approximate TV have the same behavior as the TV of the exact solution, so it begin to decrease after some time. The amplitude of the oscillations also decreases as the grid points are refined. We note that all the TV computed by our scheme are free from undesirable behavior even when the solution develops a discontinuity in the derivative. The TV of the approximate solution, for both convex and nonconvex Riemann problem never increase above the exact TV, and approaches it as the mesh is refined. These prove the TVB character of our scheme even when approximating discontinuous solution.
508
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
4.3. Two-dimensional experiments Example 3. Two-dimensional Burgers’ equation: (φ +φ +1)2 = 0, −2 x 2, −2 y 2, φt + x 2y , φ(x, y, 0) = − cos π(x+y) 2
(31)
with periodic boundary condition. The results at T = 0.5/π 2 , 1.5/π 2 are displayed in Fig. 4, with mesh 40 × 40. The solution is smooth at T = 0.5/π 2 . At time T = 1.5/π 2 , the solution has developed discontinuity in the derivative. We see that the solution converges to the viscosity solution with sharp for the corner-like discontinuity.
Example 4. Two-dimensional equation: φt − cos(φx + φy + 1) = 0, −2 x 2, −2 y 2, φ(x, y, 0) = − cos π(x+y) 2
Fig. 4. Two-dimensional Burgers’ equation; Mesh 40 × 40; (a) T = 0.5/π 2 , (b) T = 1.5/π 2 .
Fig. 5. Two-dimensional nonconvex Riemann problem; Mesh 40 × 40; (a) T = 0.5/π 2 , (b) T = 1.5/π 2 .
(32)
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
509
with periodic boundary condition and mesh 40 × 40. The result at T = 0.5/π 2 , when the solution is still smooth is displayed in Fig. 5(a). The solution has developed a discontinuous derivative at T = 1.5/π 2 , and is shown in Fig. 5(b).
Example 5 (Propagating surface). We solve the following problems of propagating surfaces [10]: φx2 + φy2 + 1 = 0, 0 x 1, 0 y 1, φ(x, y, 0) = 1 − 14 cos(2π x) − 1 cos(2πy) − 1 , φt − φx2 + φy2 + 1 = 0, 0 x 1, 0 y 1, φ(x, y, 0) = cos(2πy) − cos(2π x),
φt −
(33)
(34)
with periodic boundary condition. Figs. 6 and 7 display the evolution of the surfaces for the first and second problem, respectively, with mesh 40 × 40 at times T = 0, 0.6. The singularity in derivative has been solved with accuracy and sharpness.
Fig. 6. Propagating surface problem; First example; Mesh 40 × 40; (a) T = 0, (b) T = 0.6.
Fig. 7. Propagating surface problem; Second example; Mesh 40 × 40; (a) T = 0, (b) T = 0.6.
510
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
Fig. 8. Optimal control problem; Mesh 40 × 40; (a) T = 0.5, (b) T = 1.
Example 6 (Optimal control). We solve the following problem from optimal control [9] φt − sin(y) φx + sin(x) + sign(φy ) φy − 12 sin2 (y) − 1 − cos(x) = 0, φ(x, y, 0) = 0
(35)
in (x, y) ∈ [−π, π ] × [−π, π ] with periodic boundary conditions. Numerical solutions at times T = 0.5, 1 are presented in Fig. 8.
5. Conclusion We have introduced and tested a parabolic reconstruction for solving Hamilton–Jacobi equations. The scheme is formally third-order accurate in smooth regions, and more local in the sense that it uses only four points to reconstruct a piecewise smooth function. Numerical experiments in 1D and 2D show that our approximation is less viscous, and has good behavior in the presence of discontinuity.
Acknowledgements Professor Antonio Marquina (Universidad de Valencia) is gratefully acknowledged for support and for the critical reading of the manuscript.
References [1] M. Crandall, P. Lion, Viscosity solutions of Hamilton–Jacobi equations, Trans. Amer. Math. Soc. 277 (1983) 1–42. [2] M. Crandall, P. Lion, Two approximations of solutions of Hamilton–Jacobi equations, Math. Comp. 43 (1984) 1–19. [3] A. Harten, B. Engquist, S. Osher, S. Chakravarty, Uniformly high-order accurate essentially non-oscillatory scheme II, J. Comput. Phys. 71 (1987) 231–303. [4] G. Jiang, D. Peng, Weighted ENO schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (2000) 2126–2143.
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
511
[5] A. Kurganov, E. Tadmor, New high-resolution semi-discrete central schemes for Hamilton–Jacobi equations, J. Comput. Phys. 162 (2000) 720–742. [6] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser, Basel, 1990. [7] C.T. Lin, E. Tadmor, High-resolution non-oscillatory central schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (2000) 2163–2186. [8] A. Marquina, Local piecewise hyperbolic reconstruction of numerical fluxes for nonlinear scalar conservation laws, SIAM J. Sci. Comput. 15 (1994) 892. [9] S. Osher, J. Sethian, Front propagating with curvature dependent speed: Algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [10] S. Osher, C.W. Shu, High-order essentially non-oscillatory schemes for Hamilton–Jacobi equations, SIAM J. Numer. Anal. 28 (1991) 907–922. [11] C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes II, J. Comput. Phys. 83 (1989) 32–78. [12] Y. Stiriba, A. Marquina, R. Donat, Equilibrium real gas computations using Marquina’s scheme, Internat. J. Numer. Methods Fluids 41 (2003) 275–301.
A local piecewise parabolic method for Hamilton–Jacobi equations Youssef Stiriba 1 CERFACS/CFD Team, 42 avenue Gaspard Coriolis, 31057 Toulouse cedex 1, France
Abstract In this paper we present, analyze and implement a new local piecewise parabolic method for nonlinear Hamilton– Jacobi equations. The scheme is third-order accurate in smooth regions, uses a concept of local smoothing to prevent the excessive increase of the total variation at discontinuity, and has a local stencil in the sense that it does not extrapolate from data of the smoothest neighboring cells. One and two-dimensional numerical experiments, accuracy tests, and the behavior of the total variation of the approximate solution are presented to prove the accuracy and good resolution of our method. 2002 IMACS. Published by Elsevier Science B.V. All rights reserved.
1. Introduction We consider numerical solutions of Hamilton–Jacobi (H–J) equations φt + H (x1 , . . . , xd , t, φx1 , . . . , φxd ) = 0, φ(x, 0) = φ 0 (x),
(1)
by a third-order local parabolic method. Solutions of H–J equations are continuous, but may have discontinuous derivatives even with smooth initial condition φ 0 (x), and these solutions are in general not unique. Therefore, it is necessary to introduce the notion of entropy condition to select a unique relevant solution [1]. Recently a wide variety of finite difference schemes have been developed. For instance, in [2] the authors introduced a class of monotone schemes and proved that these schemes converge to the viscosity solution of (1). However, such methods are at most first-order accurate and too dissipative for problems of practical interest. In [9,10] the authors have used the fact that the 1D Hamilton–Jacobi equation (1) is closely related to the 1D conservation law ut + f (u)x = 0, to adapt the high-order essentially E-mail address: [email protected] (Y. Stiriba). 1 Current address: Abteilung für Numerische Mathematik; IGPM, RWTH Aachen, Templergraben 55, 52062 Aachen,
Germany. 0168-9274/02/$30.00 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. doi:10.1016/S0168-9274(02)00211-8
500
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
nonoscillatory (ENO) schemes [3,11] to solve Hamilton–Jacobi equations in conjunction with several monotone building blocks. Their numerical experiments indicate good accuracy and better resolution of the corner-like discontinuity. The same idea was reported in [4] to approximate solution of (1) by Weighted ENO (WENO) schemes. In [7,5], the authors extended their high-resolution central schemes for hyperbolic conservation laws to the Hamilton–Jacobi equations. In this paper we construct a piecewise parabolic method for approximating Hamilton–Jacobi equations (1), by applying a preprocessing of derivatives in each cell to prevent the excessive increase of the total variation of the solution. We adapt this method in conjunction with Osher and Shu local Lax– Friedrichs flux function [10] to solve H–J equations. The resulting scheme is more local in the sense that it uses only four points to reconstruct piecewise smooth functions in contrast to ENO-3 that depends on six point values, gives sharp resolution of jump derivatives, and is total variation bounded even when approximating discontinuous solutions. The outline of the paper is the following: In Section 2, we present the reconstruction procedure. In Section 3, we describe the numerical scheme for H–J equations in two space dimension. Numerical results are presented in Section 4. Finally, we close the paper with some concluding remarks.
2. The local piecewise parabolic reconstruction Before we construct a local piecewise parabolic (LPP) scheme for approximating H–J equations, we first present the LPP reconstruction procedure. Let f (x) be a piecewise smooth function. We consider a uniform grid in the x − t plane, with mesh spacing x and time step t. We denote the grid points by xi = i x and tn = t, and the computational cells by Ci = [xi−1/2 , xi+1/2 ] where xi±1/2 = (i ± 1/2) x, ∀i ∈ Z. We will also use: (1) The cell averaged of f : 1 f (ξ ) dξ, vi = x
i ∈ Z.
(2)
Ci
(2) The forward/backward differences ±(vi±1 − vi ) , i ∈ Z. di±1/2 = x The reconstruction problem is to find a parabola ci Pi (x) = ai + bi · (x − xi ) + · (x − xi )2 , 2 defined on Ci , such that it is a third-order approximation to f (x) on Ci : f (x) = Pi (x) + O x 3 every time f (x) is smooth in Ci . The parabola Pi (x) is required to satisfy 1 (1) vi = x Ci Pi (ξ ) dξ , (2) di+1/2 = Pi (xi+1/2 ),
(3)
(4)
(5)
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
501
for every i. These imply that the reconstruction is third-order accurate. For each cell, there exists a unique parabola of the form (4) determined from vi , di−1/2 , and di+1/2 by di+1/2 − di−1/2 di+1/2 + di−1/2 , bi = , x 2 and satisfies di+1/2 + di−1/2 . P (xi ) = 2 The lateral derivatives of Pi are ci =
P (xi+1/2 ) = bi + ci · x
and
and
ai = vi − ci
x 2 , 24
P (xi−1/2 ) = bi − ci · x.
(6)
(7)
(8)
Since the reconstruction is confined in each computational cell and repeated every time step, the local total variation (LTV), defined by LTV i = TV(Pi ) (the total variation of the function P in the cell Ci , see [6]), must be controlled. A method is called local total variation bounded (LTVB) if LTV(Pi ) M · x
∀i,
for some fixed positive M that depends only on the initial conditions. The LTV for our reconstruction is x · |d | + |d | if di+1/2 · di−1/2 0, i+1/2 i−1/2 2 LTV(Pi ) = x 2·|di+1/2 |·|di−1/2 | · |di+1/2 | + |di−1/2 | − |di+1/2 |+|di−1/2 | elsewhere, 2
(9)
(10)
for every i. It is clear that the present parabola is not LTVB. So, the LTV may blow up near discontinuities (i.e., for some i such that di+1/2 = O( x −1 )). This behavior has been observed when approximating Burgers’ equation (when the waves break into discontinuities). The total variation of the function f (x) on a cell Ci where it is smooth depends essentially, up to O( x 3 ), on the size of di±1/2 (this follows from trapezoidal rule). Thus to prevent the increase of LTV, we shall correct the lateral derivatives di+1/2 and di−1/2 by a pair of local consistent numerical derivatives {dli , dri } such that (11) dli = di−1/2 + O x 2 and dri = di+1/2 + O x 2 , which preserve third order accuracy. Below we outline the LPP reconstruction. To construct the parabola in a nontransition cell Ci (i.e., when di+1/2 · di−1/2 > 0), we can use (6) and assign to b the average di+1/2 + di−1/2 . (12) 2 The resulting parabola has as derivatives at the cell interfaces di+1/2 and di−1/2 . However, with this choice the LTV may increase. So, for transition cells Ci with di+1/2 = O( x −1 ), the LTV is (13) LTV(Pi ) = O x −1 . b=
Thus, the method becomes unstable, and more smoothing is necessary to prevent the excessive increase of the local total variation of the solution, as observed in our numerical experiments. To provide more smoothing, we can use the harmonic mean of the lateral derivatives 2 · di+1/2 · di−1/2 , (14) b= di+1/2 + di−1/2
502
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
as in [8], used to construct local hyperbola, such that the derivatives of the parabola interpolates b at the cell center xi and the lateral grid derivative with smallest absolute value. This gives a parabola with lateral derivative smaller than the original one in nonsmooth side. The harmonic mean is smaller in absolute value than (12), is a second-order approximation of the derivative at the cell center xi , and gives more smoothing and good results according to our numerical experiments (see also Section 4 and [12]). It is possible to choose other means such that the method of reconstruction is LTVB, but we have found, through our numerical experiments, that our choice gives good results. In transition cells (i.e., when di+1/2 ·di−1/2 < 0), we follow the ENO idea [11] and change the derivative with largest absolute value by the other one multiplied by x 2 . The algorithm in the flavor of FORTRAN 77 is the following Algorithm1 (LPP reconstruction). *tol = O x 2 if (|di−1/2 | tol) and (|di+1/2 | tol) then bi = 0 and ci = 0 else if (|di−1/2 | tol) or (Ci is a transition cell with |di−1/2 | |di+1/2 |) bi = 2 · di+1/2 · else
x 2 1+ x 2
and ci =
di+1/2 · x 2 −bi x
if (|di+1/2 | tol) or (Ci is a transition cell with |di+1/2 | |di−1/2 |) bi = 2 · di−1/2 · else
2·d
x 2 1+ x 2
and ci =
bi −di−1/2 · x 2 x
·d
i−1/2 i+1/2 bi = di−1/2 +di+1/2 if (|di−1/2 | |di+1/2 |) then b −di−1/2 ci = i 2· x else d −bi ci = i+1/2 2· x
We compute ai from vi and ci using (6). We note that we can also implement our algorithm by using the harmod slope limiter, which consist of the harmonic mean on the lateral slopes when they have the same sign and zero when they have different signs. The algorithm is the following: Algorithm 2 (LPP-HARMOD reconstruction). if (di−1/2 · di+1/2 tol) then bi = 0 and ci = 0 else if (|di−1/2 | tol) and (|di+1/2 | tol) then d +d d −d bi = i−1/2 2 i+1/2 and ci = i+1/2 x i−1/2 else 2·di−1/2 ·di+1/2 bi = di−1/2 +di+1/2 if (|di−1/2 | |di+1/2 |) then b −di−1/2 ci = i 2· x
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
else ci =
503
di+1/2 −bi 2· x
The LPP-HARMOD has a much simpler form, and gives similar results to those obtained with LPP approximation, according to our numerical experiments. Remark 1. The LPP reconstructions is more local in the sense that it uses only four points to reconstruct piecewise smooth functions in contrast with other reconstruction like ENO-3 or WENO-5 that depend on six point values. This has good consequence in reducing numerical noise generated by singularities during numerical process (see Section 4). Remark 2. We note also that our methods does not use second differences to compute the corner like discontinuity (a jump in the first derivative). The following proposition shows that the LPP reconstruction is LTVB. Proposition 1. The LPP method of reconstruction of the function f (x) is LTVB. Proof. Since f (x) is a piecewise smooth function, we can show, as in Theorem 2.1 of [8], that there exists a constant M > 0 that depends on f (x) in smooth regions except for finite number of i for which di+1/2 = O( x −1 ), such that |di+1/2 | M.
(15)
If Ci a nontransition cell and |di−1/2 | |di+1/2 |, then the preprocessed derivatives at cell interfaces are: dli = di−1/2
and
dri = 2bi − di−1/2 .
(16)
Hence from (10) x |di−1/2 | + |2bi − di−1/2 | 3 · M · x. 2 We use the same argument for transition cell. LTV(Pi )
(17)
3. The LPP method for Hamilton–Jacobi equations In this section, we consider the 2D Hamilton–Jacobi equation φt + H (φx , φy ) = 0, (18) φ(x, y, 0) = φ 0 (x, y) where (x, y) ∈ [a, b] × [c, d], t 0, H is a locally Lipschitz continuous Hamiltonian, and the initial data φ 0 (x, y) is a locally Lipschitz continuous. n denote the numerical Let x and y be the grid spacing in x and y directions, respectively. Let φi,j approximation to the viscosity solution φ(i x, j y, n t) of (18). We use the method of line, by discretizing only in space and leaving the problem continuous in time, in designing our algorithm. In such approach the semi-discrete form of the LPR method is + dφi,j − + − φx,i,j = −H , φx,i,j , φy,i,j , φy,i,j . (19) dt
504
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
± Here φx,i,j are LPP approximations to the left and right x-derivative of φ(x, y, t) at (xi , yj , tn ). Similarly ± for φy,i,j . The numerical Hamiltonian is a Lipschitz continuous monotone flux (i.e., nonincreasing in its fifth and seventh arguments and nondecreasing in its sixth and eighth arguments), and consistent with H
(u, u, v, v) = H (u, v). H
(20)
There are many Hamiltonian that may be used as a building block (see [10]). In our numerical experiments we use the local Lax–Friedrichs (LLF) flux + − + − 1 + − + − v 1 + u + v u u ,u ,v ,v = H , − α u+ − u− − β v + − v − , (21) H 2 2 2 2 where
H1 (u, v) , max u∈I (u− ,u+ ),v∈[C,D] I (u, v) = min(u, v), max(u, v) ,
α=
H1 (u, v) =
∂H (u, v) , ∂u
β=
H2 (u, v) =
max
H2 (u, v) ,
(22)
v∈I (v − ,v + ),u∈[A,B]
(23) ∂H (u, v) . ∂v
(24)
Here [A, B] is the value range for u± and [C, D] is the value range for v ± . If we replace I (u− , u+ ) and I (v − , v + ) by [A, B] and [C, D], respectively, we recover the Lax–Friedrichs flux. Time discretization is done by the third-order TVD (total variation diminishing) Runge–Kutta of Shu and Osher [11]. Given the numerical solution at time n · t, {φin }, then the numerical solution at the new time (n + 1) · t, {φin+1 }, is computed by the following recursive formula φi(0) , φi(1) = φi(0) + t H
φi(0) = φin ,
1 2 2 (2) φi , φi(3) = φi(0) + φi(2) + t H 3 3 3
3 1 1 (1) φi , φi(2) = φi(0) + φi(1) + t H 4 4 4 φi(n+1) = φi(3) ,
i denotes the numerical Hamiltonians computed from the point values φi for all i. Here H sponding to substep j , where j = 1, 2, 3. (j )
(j )
corre-
4. Numerical examples In this section we study the performance of our numerical scheme by applying it to some test problems. 4.1. Accuracy tests Example 1. We solve the one-dimensional H–J equations: φt + H (φx ) = 0, −1 x 1, φ(x, 0) = − cos(π x),
(25)
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
505
Table 1 L1 and L∞ errors and orders for 1D Burgers’ equation N
L1 error
L1 order
L∞ error
L∞ order
0.331E–01 0.133E–01 0.261E–02 0.441E–03
– 1.31 2.34 2.56
0.449E–01 0.152E–01 0.281E–02 0.691E–03
– 1.70 2.43 2.18
L∞ error
L∞ order
0.824E–01 0.225E–02 0.514E–03 0.118E–03
– 1.87 2.13 2.12
0.224E–01 0.510E–02 0.127E–02 0.288E–03
– 2.13 2.00 2.14
T = 0.5/π 2 10 20 40 80
0.154E–01 0.265E–02 0.432E–03 0.635E–04
– 2.53 2.61 2.76 T = 1.5/π 2
10 20 40 80
0.216E–01 0.338E–02 0.454E–03 0.441E–04
– 2.67 2.89 3.36
Table 2 L1 and L∞ errors and orders for 1D nonconvex problem N
L1 error
L1 order T = 0.5/π 2
10 20 40 80
0.305E–02 0.878E–03 0.117E–03 0.177E–04
– 1.79 2.90 2.72 T = 1.5/π 2
10 20 40 80
0.868E–02 0.176E–02 0.292E–03 0.441E–04
– 2.30 2.59 2.72
with a convex Hamiltonian H (Burgers’ equation): (u + α)2 , 2 and a nonconvex Hamiltonian H : H (u) =
(26)
H (u) = − cos(u + α).
(27)
For simplicity, we take α = 1. For both cases the solution is smooth at time T = 0.5/π 2 , and develop a discontinuous derivative at time T = 1.5/π 2 . Let φin be the reconstructed solution at (xi , t n ). The discrete L1 and L∞ norms are given by N
n
φ(xi , tn ) − φ n ,
e ·, t 1 = x i L
(28)
i=1
n
e ·, t ∞ = max φ(xi , tn ) − φ n . i L 1iN
(29)
506
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
Fig. 1. One-dimensional Burgers’ equation; (a) before singularity, T = 0.5/π 2 ; (b) after singularity, T = 1.5/π 2 .
Fig. 2. One-dimensional nonconvex problem; (a) before singularity, T = 0.5/π 2 ; (b) after singularity, T = 1.5/π 2 .
Here N is the number of grid points used. In Tables 1 and 2, we list L1 - and L∞ -errors and convergence rates for both convex and nonconvex Hamiltonians. The exact solution can be obtained using the method of characteristic (see [9]). We observe L1 -third-order-accuracy and more than second-order accuracy in L∞ -norm. Numerical results together with exact solutions are displayed in Figs. 1 and 2. We see that the solution is free from spurious oscillations, and the scheme has good behavior at the corner-likediscontinuity. 4.2. The behavior of the total variation Example 2. We numerically examine the behavior of the total variation N
n
φ − φ n
TV φ n := i+1 i i=1
(30)
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
507
Fig. 3. Total variation of the solution; N = 40 (‘·’), N = 80 (‘- · -’), N = 160, (‘-’), exact TV (‘- -’); (a) Burgers’ equation; (b) Nonconvex problem.
of the discrete solution. Figs. 3 and 4 display the TV of the solution with N = 40, 80, 160. We observe that approximate TV have the same behavior as the TV of the exact solution, so it begin to decrease after some time. The amplitude of the oscillations also decreases as the grid points are refined. We note that all the TV computed by our scheme are free from undesirable behavior even when the solution develops a discontinuity in the derivative. The TV of the approximate solution, for both convex and nonconvex Riemann problem never increase above the exact TV, and approaches it as the mesh is refined. These prove the TVB character of our scheme even when approximating discontinuous solution.
508
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
4.3. Two-dimensional experiments Example 3. Two-dimensional Burgers’ equation: (φ +φ +1)2 = 0, −2 x 2, −2 y 2, φt + x 2y , φ(x, y, 0) = − cos π(x+y) 2
(31)
with periodic boundary condition. The results at T = 0.5/π 2 , 1.5/π 2 are displayed in Fig. 4, with mesh 40 × 40. The solution is smooth at T = 0.5/π 2 . At time T = 1.5/π 2 , the solution has developed discontinuity in the derivative. We see that the solution converges to the viscosity solution with sharp for the corner-like discontinuity.
Example 4. Two-dimensional equation: φt − cos(φx + φy + 1) = 0, −2 x 2, −2 y 2, φ(x, y, 0) = − cos π(x+y) 2
Fig. 4. Two-dimensional Burgers’ equation; Mesh 40 × 40; (a) T = 0.5/π 2 , (b) T = 1.5/π 2 .
Fig. 5. Two-dimensional nonconvex Riemann problem; Mesh 40 × 40; (a) T = 0.5/π 2 , (b) T = 1.5/π 2 .
(32)
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
509
with periodic boundary condition and mesh 40 × 40. The result at T = 0.5/π 2 , when the solution is still smooth is displayed in Fig. 5(a). The solution has developed a discontinuous derivative at T = 1.5/π 2 , and is shown in Fig. 5(b).
Example 5 (Propagating surface). We solve the following problems of propagating surfaces [10]: φx2 + φy2 + 1 = 0, 0 x 1, 0 y 1, φ(x, y, 0) = 1 − 14 cos(2π x) − 1 cos(2πy) − 1 , φt − φx2 + φy2 + 1 = 0, 0 x 1, 0 y 1, φ(x, y, 0) = cos(2πy) − cos(2π x),
φt −
(33)
(34)
with periodic boundary condition. Figs. 6 and 7 display the evolution of the surfaces for the first and second problem, respectively, with mesh 40 × 40 at times T = 0, 0.6. The singularity in derivative has been solved with accuracy and sharpness.
Fig. 6. Propagating surface problem; First example; Mesh 40 × 40; (a) T = 0, (b) T = 0.6.
Fig. 7. Propagating surface problem; Second example; Mesh 40 × 40; (a) T = 0, (b) T = 0.6.
510
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
Fig. 8. Optimal control problem; Mesh 40 × 40; (a) T = 0.5, (b) T = 1.
Example 6 (Optimal control). We solve the following problem from optimal control [9] φt − sin(y) φx + sin(x) + sign(φy ) φy − 12 sin2 (y) − 1 − cos(x) = 0, φ(x, y, 0) = 0
(35)
in (x, y) ∈ [−π, π ] × [−π, π ] with periodic boundary conditions. Numerical solutions at times T = 0.5, 1 are presented in Fig. 8.
5. Conclusion We have introduced and tested a parabolic reconstruction for solving Hamilton–Jacobi equations. The scheme is formally third-order accurate in smooth regions, and more local in the sense that it uses only four points to reconstruct a piecewise smooth function. Numerical experiments in 1D and 2D show that our approximation is less viscous, and has good behavior in the presence of discontinuity.
Acknowledgements Professor Antonio Marquina (Universidad de Valencia) is gratefully acknowledged for support and for the critical reading of the manuscript.
References [1] M. Crandall, P. Lion, Viscosity solutions of Hamilton–Jacobi equations, Trans. Amer. Math. Soc. 277 (1983) 1–42. [2] M. Crandall, P. Lion, Two approximations of solutions of Hamilton–Jacobi equations, Math. Comp. 43 (1984) 1–19. [3] A. Harten, B. Engquist, S. Osher, S. Chakravarty, Uniformly high-order accurate essentially non-oscillatory scheme II, J. Comput. Phys. 71 (1987) 231–303. [4] G. Jiang, D. Peng, Weighted ENO schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (2000) 2126–2143.
Y. Stiriba / Applied Numerical Mathematics 45 (2003) 499–511
511
[5] A. Kurganov, E. Tadmor, New high-resolution semi-discrete central schemes for Hamilton–Jacobi equations, J. Comput. Phys. 162 (2000) 720–742. [6] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser, Basel, 1990. [7] C.T. Lin, E. Tadmor, High-resolution non-oscillatory central schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (2000) 2163–2186. [8] A. Marquina, Local piecewise hyperbolic reconstruction of numerical fluxes for nonlinear scalar conservation laws, SIAM J. Sci. Comput. 15 (1994) 892. [9] S. Osher, J. Sethian, Front propagating with curvature dependent speed: Algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [10] S. Osher, C.W. Shu, High-order essentially non-oscillatory schemes for Hamilton–Jacobi equations, SIAM J. Numer. Anal. 28 (1991) 907–922. [11] C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes II, J. Comput. Phys. 83 (1989) 32–78. [12] Y. Stiriba, A. Marquina, R. Donat, Equilibrium real gas computations using Marquina’s scheme, Internat. J. Numer. Methods Fluids 41 (2003) 275–301.