- Email: [email protected]
A TVD-type method for 2D scalar Hamilton–Jacobi equations on unstructured meshes
Journal of Computational and Applied Mathematics 195 (2006) 182 – 191 www.elsevier.com/locate/cam
A TVD-type method for 2D scalar Hamilton–Jacobi equations on unstructured meshes夡 Lingyan Tang∗ , Songhe Song Department of Mathematics and System Science, Science School, National University of Defense Technology, Changsha, Hunan 410073, China Received 15 August 2004; received in revised form 17 March 2005
Abstract In this paper, a TVD-type difference scheme which satisfies maximum principle is developed for 2D scalar Hamilton–Jacobi equations on unstructured triangular meshes. The main ideas are node-based approximations and derivative-limited reconstruction with quadratic interpolation polynomial. The solution’s slope satisfies maximum principle. Numerical experiments are performed to demonstrate high-order accuracy in smooth fields and good resolution of derivative singularities. The new method is simpler than WENO. © 2005 Elsevier B.V. All rights reserved. Keywords: Hamilton–Jacobi equation; Unstructured mesh; Maximum principle
1. Introduction Hamilton–Jacobi (H–J) equations are of special interest in a variety of applications, such as optimal control, differential games, image processing and geometric optics. With the popularity of level set methods, the necessity to find good algorithms to solve H–J equations becomes even more obvious.
夡
Supported in part by the NFS of China (10001038).
∗ Corresponding author.
E-mail address: [email protected] (L. Tang). 0377-0427/$ - see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2005.03.082
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
183
As we all know, solutions of H–J equations are very complex. They are typically continuous but with discontinuous derivatives even when the initial condition is smooth, and unfortunately, such kinds of discontinuities are in general not unique. Therefore, a mechanism is required for singling out a “physically relevant solution”, the viscosity solution. For convex Hamiltonians the viscosity solution coincides with the limit solution obtained by the vanishing viscosity method [5]. Extensions to general Hamiltonians were introduced by Crandall and Lions in [3]. H–J equations are closely related to hyperbolic conservation laws. In certain sense H–J equations are easier to solve than their conservation laws counterparts, because the solutions here are smoother. Up to now, most numerical experiments have used Cartesian-based algorithms and regular meshes restricting the shape of calculational domain to simple configurations like squares or rectangles, for example, the earlier second-order ENO schemes of Osher and Sethian [7], higher-order essentially nonoscillatory (ENO) schemes of Osher and Shu [8], higher-order weighted ENO (WENO) schemes of Jiang and Peng [4], and central high-resolution schemes of Lin and Tadmor [6], among many others. However, for unstructured meshes, algorithms, especially high-order algorithms are relatively few. The problem is that the H–J equations cannot be written in a “conservation form”, suitable for integration by parts which is the backbone of finite volume method and finite element method. We have the very good first- and second-order finite volume type schemes of Abgrall [1] and ENO or limiter-type high-order version of Augoula and Abgrall [2]. The monotone Hamiltonians developed in [1] are used in this paper as building blocks. In [11], Zhang and Shu constructed WENO schemes for solving the nonlinear H–J equations on two-dimensional unstructured meshes. By using a strategy to choose diversified smaller stencils to make up the bigger stencil in the WENO procedure, their schemes which combine secondorder approximations with nonlinear weights can get third/fourth-order accuracy. Numerical experiments show that the WENO scheme is a high-order, high-resolution scheme, but its algorithm is very complex and it is not convenient for application. In this paper, we develop a TVD-type difference scheme on two-dimensional unstructured triangular meshes for scalar H–J equations. The scheme is based on second-order central scheme and locally derivative-limited reconstruction for quadratic interpolation polynomial. As the solution’s slope is locally limited, the new scheme has high resolution of derivative singularities. The derivative limited methodology adopted in this paper can be traced back to the earlier work for 2D conservation laws started in [10] for a second-order finite-volume version which satisfies maximum principle. Because H–J equations are in some sense the conservation laws “integrated once”, quadratic reconstruction here uses derivatives one order higher than those for the latter. Compared with the WENO scheme [11], our scheme avoids complex stencil-searching procedure and so saves vast CPU times. The new scheme is uniformly second-order accuracy which is sufficient for practice in smooth region, even near extremum. Numerical experiments demonstrate its convergence to viscosity solutions of H–J equation. The rest of the paper is organized as follows. In Section 2, a detailed description of algorithm is presented, we use central scheme to discretize spatial term in Section 2.1. In Section 2.2, a derivative limiter is constructed to modify the quadratic interpolation polynomial obtained in Section 2.1 and time discretization scheme is given in Section 2.3. In Section 3, some numerical examples are performed. The conclusion is given in Section 4.
184
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
Fig. 1. Node i and its angular sectors.
2. Numerical method 2.1. Spatial discretization In this paper, we consider the numerical solution of H–J equations ju j u ut + H , = 0. jx jy
(1)
Eq. (1) is solved in the domain , which has a triangulation h consisting of triangles. The nodes will be named by their indices 0 i N, with a total number N +1. As shown in Fig. 1, for every node i, we define the ki + 1 angular sectors meeting at the point i as T0 , . . . , Tki ; They are the inner angles of the triangles which have i as a vertex. At node i, the indexing of the angular sectors is ordered counterclockwise. − → n l+1/2 denotes the unit vector of the half-line Dl+1/2 = Tl ∩ Tl+1 , and l the inner angle of sector Tl , 0 l ki . We will denote by ui the numerical approximation to the viscosity solution of (1) at node i. (∇u)0 , . . . , (∇u)ki will, respectively, represent the numerical approximation of ∇u at node i in each angular sector T0 , . . . , Tki . An important building block for the algorithms in this paper is the Lax–Friedrichs-type monotone Hamiltonian for arbitrary triangulations developed by Abgrall in [1]. It is given by ((∇u)0 , . . . , (∇u)k ) = H H i
k
l=0 l (∇u)l i
2
−
ki
l=0
l+1/2
(∇u)l + (∇u)l+1 2
→ ·− n l+1/2 , (2)
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
185
Fig. 2. Triangle l and its neighbor.
where l+1/2 = tan(l /2) + tan(l+1 /2), = max{max |H1 (u, v)|, max |H2 (u, v)|}. Here H1 and H2 are the partial derivatives of H with respect to ux and uy , respectively, or the Lipschitz constants of H with respect to ux and uy , if H is not differentiable. in (2) defines a monotone Hamiltonian. It is Lipschitz continuous in all arguments and is The H (∇u, . . . , ∇u) = H (∇u). Hence if we have high-order approximations to ∇u consistent with H, i.e., H will be a high-order approximation to H. at node i in every angular sector, the numerical Hamiltonian H Using the Lax–Friedrichs Hamiltonian (2), we can write (1) into the following semi-discrete form: d ((∇u)0 , . . . , (∇u)k ) = 0. ui (t) + H i dt
(3)
For time evolution, TVD Runge–Kutta scheme of Shu and Osher [9] is adopted. This article focuses for (3). primarily on designing an appropriate approximation of ∇u which can fetch out a gratifying H First, we use the central stencil to construct a second-order interpolation polynomial on each triangle. Namely, for triangle l , see Fig. 2, we interpolate on the six nodes: Ai , Bi (i = 1, 2, 3). By solving a linear system, we obtain v l (x, y) = a1 (x − xcl )2 + a2 (x − xcl )(y − ycl ) + a3 (y − ycl )2 + a4 (x − xcl ) + a5 (y − ycl ) + a6 ,
(4)
where Cl = (xcl , ycl ) is the barycenter of l . Because central stencil lacks the capability of treating with singularity, spurious oscillations will occur if we use v l (x, y) directly in (3). Next, we will introduce our novel reconstruction procedure which is the key component of our scheme.
186
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
2.2. The quadratic reconstruction In this section, we will construct a new polynomial pl (x, y) in each triangular which is a modification of v l (x, y) obtained in Section 2.1. In order to exclude oscillation, the partial derivatives of p l (x, y) are subjected to a so-called “limiter”. Due to the operation of this limiter, the slope of pl (x, y) satisfies the maximum principle (as proven in [10]). In smooth region, pl (x, y) is almost the same as v l (x, y). While at the discontinuous points of ∇u, pl (x, y) is modified to avoid new extremum. In the following, we will describe the reconstruction of pl (x, y) in detail. For every triangle l , let N(l ) denote the set of barycenters of triangles which have at least one point meeting with l . Define Mx l =
max
Cj ∈N (j )
(vxl (Cj )),
mx l =
min
Cj ∈N (j )
(vxl (Cj )).
From (4), the x-directional derivative of v l (x, y) is vxl = 2a1 (x − xcl ) + a2 (y − ycl ) + a4
on l .
(5)
To guarantee that the x-directional derivative of pl (x, y) does not suffer from new extremum, we construct a x-directional derivative limiter x for pl (x, y). Thus we obtain pxl (x, y) = x · (2a1 (x − xcl ) + a2 (y − ycl )) + a4 where x = min(x (Ai ))(i = 1, 2, 3). ⎧ Mx l − vxl (Cl ) ⎪ ⎪ min 1, l , ⎪ ⎪ vx (Ai ) − vxl (Cl ) ⎨ mx l − vxl (Cl ) x (Ai ) = ⎪ , min 1, l ⎪ ⎪ vx (Ai ) − vxl (Cl ) ⎪ ⎩ 1,
on l ,
(6)
vxl (Ai ) > vxl (Cl ), vxl (Ai ) < vxl (Cl ),
(7)
vxl (Ai ) = vxl (Cl ).
It is obvious that pxl (x, y) satisfies the maximum principle pxl (x, y) ∈ [mx l , Mx l ]
on l .
(8)
In the same way, We can define My l and my l , and construct a y-directional derivative limiter y to modify vyl (x, y). Hence, the y-directional derivative of pl (x, y) is given by pyl (x, y) = y · (a2 (x − xcl ) + 2a3 (y − ycl )) + a5
on l
(9)
and pyl (x, y) satisfies the local maximum principle too, i.e., pyl (x, y) ∈ [my l , My l ]
on l .
(10)
Gathering from (4), (6) and (9), we obtain pl (x, y) = x · a1 (x − xcl )2 + m · a2 (x − xcl )(y − ycl ) + y · a3 (y − ycl )2 + a4 (x − xcl ) + a5 (y − ycl ) + a6 , where m = min(x , y ).
(11)
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
187
Remark. From (6) and (9), we get two limiters (x and y ) of coefficient a2 in (4). Here we choose the smaller one for (11) to prevent new extremum of p l (x, y)’s two partial derivatives. So far, we get our final approximation p l (x, y) for u and (3) can be written as d ((∇p 0 ), . . . , (∇p k )). ui (t) = L(ui ) = −H dt
(12)
From (8) and (10), we can see that p l (x, y) satisfies the maximum principle; on the other hand, p l (x, y) has no numerical accuracy loss compared with v l (x, y), so our scheme is a second-order TVD scheme. Besides a second-order central stencil interpolation, our reconstruction procedure only involves simple modification, so it is very fast. 2.3. Temporal discretization To keep the same numerical accuracy of spatial discretization, we use the second-order TVD Runge–Kutta scheme of Shu and Osher [9] for the temporal discretization. The algorithm is given by (1)
(0)
(0)
ui = ui + t L(ui ), (2)
(0)
(13)
(1)
(1)
ui = 21 ui + 21 ui + 21 t L(ui )
(14)
with u(0) = un and un+1 = u(2) . 3. Numerical examples In this section, we apply the new scheme developed in Section 2 to some two-dimensional model problems. The triangular mesh is shown in Fig. 3. Let N be the total number of calculational grids. Example 1. Two-dimensional Burgers equation: ut +
(u2x + u2y + 1)2 2
u(x, y, 0) = − cos
= 0,
−2 x 2, −2 y 2,
(x + y)
2
(15)
with periodic boundary conditions. We take = 6.0 and N = 512 in (2) and show the results in Fig. 4. We can see the solution develops a discontinuous derivative when t = 1.5/2 . Example 2. Two-dimensional burning problem: ut − u2x + u2y + 1 = 0, 0 x 1, 0 y 1, u(x, y, 0) = cos(2x) − cos(2y).
(16)
188
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
U
Fig. 3. The calculational triangular meshes.
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 2
-2 -1
1
0 X
0 -1
-2
1 2
Y
Fig. 4. The two-dimensional Burgers equation, H (u, v) = (u2 + u2 + 1)2 /2. Top: t = 0; bottom: t = 1.5/2 .
We take =
√
2/2 and N = 512 in (2). The results are shown in Figs. 5 and 6 .
Example 3. Two-dimensional eikonal equation with a nonconvex Hamiltonian, which arises in geometric optics: ut + u2x + u2y + 1 = 0, 0 x 1, 0 y 1, u(x, y, 0) = 0.25(cos(2x) − 1)(cos(2y) − 1) − 1. √ We take = 2/2 and N = 2048 in (2). The results are shown in Fig. 7.
(17)
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
189
2 1.5 1 0.5 U
0 -0.5 -1 -1.5 1 0.5 Y
0
0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 X
0
2 1.5 1 U
0.5 0 -0.5 -1
Y
-1.5 1 0.8 0.6 0.4 0.2 0
0
0.9 1 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 X
Fig. 5. The burning problem, H (u, v) = − u2 + v 2 + 1. Top: t = 0.5/2 ; bottom: t = 1.5/2 . 16
16
14
14
12
12
10
10
8
8
6
6
4
4
2
2 2
4
6
8
10
12
14
16
2
4
6
8
Fig. 6. Isoline of the burning problem. Left: t = 0.5/2 ; right: t = 1.5/2 .
10
12
14
16
190
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191 0 -0.2 T he
U
-0.4 -0.6 -0.8 -1 1 0.8 0.6 0.4 0.2
Y
0
0
0.2
0.4
0.6
0.8
1
X
-1.4
U
-1.45 -1.5 -1.55 -1.6 1 0.8 0.6 0.4 Y
0.2 0
0
0.2
0.4
0.6
0.8
1
X
Fig. 7. Two-dimensional eikonal equation, H (u, v) = u2 + v 2 + 1. Top: t = 0; bottom: t = 0.6.
4. Concluding remark We developed a TVD-type difference scheme on unstructured triangular meshes for 2D scalar H–J equations. With its solution’s slope satisfying maximum principle, the scheme has two-order numerical accuracy in smooth region. Besides, no parameters have to be turned in our method. Numerical experiments show that the new method has high resolution of derivative singularities and is significantly simpler than WENO. The method should be useful if geometry or adaptivity requires unstructured meshes for the solution of H–J equations. References [1] R. Abgrall, Numerical discretization of the first-order Hamilton–Jacobi equation on triangular meshes, Comm. Pure Appl. Math. 49 (1996) 1339–1373. [2] S. Augoula, R. Abgrall, High order numerical discretization for Hamilton–Jacobi equations on triangular meshes, J. Sci. Comput. 15 (2000) 197–229. [3] M. Crandall, P. Lions, Viscosity solutions of Hamilton–Jacobi equations, Trans. Amer. Math. Section 277 (2) (1983) 1–42. [4] G. Jiang, D.P. Peng, Weighted ENO schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (2000) 2126–2143.
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
191
[5] S.N. Kruzkov, The Cauchy problem in the large for nonlinear equations and for certain quasilinear systems of the first order with several variables, Soviet Math. Dokl. 5 (1964) 493–496. [6] C.T. Lin, E. Tadmor, High-resolution non-oscillatory central schemes for approximate Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (2000) 2163–2186. [7] S. Osher, J. Sethian, Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [8] S. Osher, C.W. Shu, High-order essentially schemes for Hamilton–Jacobi equations, SIAM J. Numer. Anal. 28 (1991) 907–922. [9] C. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys. 77 (1988) 439–471. [10] S. Song, Y. Li, A class of finite volume scheme satisfying maximum principle for 2D scalar hyperbolic conservation laws on unstructured meshes, Chinese J. Numer. Math. Appl. 19 (4) (1997). [11] Y. Zhang, C. Shu, High order WENO schemes for Hamilton–Jacobi equations on triangular meshes, SIAM J. Sci. Comput. 24 (2003) 1005–1030.
A TVD-type method for 2D scalar Hamilton–Jacobi equations on unstructured meshes夡 Lingyan Tang∗ , Songhe Song Department of Mathematics and System Science, Science School, National University of Defense Technology, Changsha, Hunan 410073, China Received 15 August 2004; received in revised form 17 March 2005
Abstract In this paper, a TVD-type difference scheme which satisfies maximum principle is developed for 2D scalar Hamilton–Jacobi equations on unstructured triangular meshes. The main ideas are node-based approximations and derivative-limited reconstruction with quadratic interpolation polynomial. The solution’s slope satisfies maximum principle. Numerical experiments are performed to demonstrate high-order accuracy in smooth fields and good resolution of derivative singularities. The new method is simpler than WENO. © 2005 Elsevier B.V. All rights reserved. Keywords: Hamilton–Jacobi equation; Unstructured mesh; Maximum principle
1. Introduction Hamilton–Jacobi (H–J) equations are of special interest in a variety of applications, such as optimal control, differential games, image processing and geometric optics. With the popularity of level set methods, the necessity to find good algorithms to solve H–J equations becomes even more obvious.
夡
Supported in part by the NFS of China (10001038).
∗ Corresponding author.
E-mail address: [email protected] (L. Tang). 0377-0427/$ - see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2005.03.082
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
183
As we all know, solutions of H–J equations are very complex. They are typically continuous but with discontinuous derivatives even when the initial condition is smooth, and unfortunately, such kinds of discontinuities are in general not unique. Therefore, a mechanism is required for singling out a “physically relevant solution”, the viscosity solution. For convex Hamiltonians the viscosity solution coincides with the limit solution obtained by the vanishing viscosity method [5]. Extensions to general Hamiltonians were introduced by Crandall and Lions in [3]. H–J equations are closely related to hyperbolic conservation laws. In certain sense H–J equations are easier to solve than their conservation laws counterparts, because the solutions here are smoother. Up to now, most numerical experiments have used Cartesian-based algorithms and regular meshes restricting the shape of calculational domain to simple configurations like squares or rectangles, for example, the earlier second-order ENO schemes of Osher and Sethian [7], higher-order essentially nonoscillatory (ENO) schemes of Osher and Shu [8], higher-order weighted ENO (WENO) schemes of Jiang and Peng [4], and central high-resolution schemes of Lin and Tadmor [6], among many others. However, for unstructured meshes, algorithms, especially high-order algorithms are relatively few. The problem is that the H–J equations cannot be written in a “conservation form”, suitable for integration by parts which is the backbone of finite volume method and finite element method. We have the very good first- and second-order finite volume type schemes of Abgrall [1] and ENO or limiter-type high-order version of Augoula and Abgrall [2]. The monotone Hamiltonians developed in [1] are used in this paper as building blocks. In [11], Zhang and Shu constructed WENO schemes for solving the nonlinear H–J equations on two-dimensional unstructured meshes. By using a strategy to choose diversified smaller stencils to make up the bigger stencil in the WENO procedure, their schemes which combine secondorder approximations with nonlinear weights can get third/fourth-order accuracy. Numerical experiments show that the WENO scheme is a high-order, high-resolution scheme, but its algorithm is very complex and it is not convenient for application. In this paper, we develop a TVD-type difference scheme on two-dimensional unstructured triangular meshes for scalar H–J equations. The scheme is based on second-order central scheme and locally derivative-limited reconstruction for quadratic interpolation polynomial. As the solution’s slope is locally limited, the new scheme has high resolution of derivative singularities. The derivative limited methodology adopted in this paper can be traced back to the earlier work for 2D conservation laws started in [10] for a second-order finite-volume version which satisfies maximum principle. Because H–J equations are in some sense the conservation laws “integrated once”, quadratic reconstruction here uses derivatives one order higher than those for the latter. Compared with the WENO scheme [11], our scheme avoids complex stencil-searching procedure and so saves vast CPU times. The new scheme is uniformly second-order accuracy which is sufficient for practice in smooth region, even near extremum. Numerical experiments demonstrate its convergence to viscosity solutions of H–J equation. The rest of the paper is organized as follows. In Section 2, a detailed description of algorithm is presented, we use central scheme to discretize spatial term in Section 2.1. In Section 2.2, a derivative limiter is constructed to modify the quadratic interpolation polynomial obtained in Section 2.1 and time discretization scheme is given in Section 2.3. In Section 3, some numerical examples are performed. The conclusion is given in Section 4.
184
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
Fig. 1. Node i and its angular sectors.
2. Numerical method 2.1. Spatial discretization In this paper, we consider the numerical solution of H–J equations ju j u ut + H , = 0. jx jy
(1)
Eq. (1) is solved in the domain , which has a triangulation h consisting of triangles. The nodes will be named by their indices 0 i N, with a total number N +1. As shown in Fig. 1, for every node i, we define the ki + 1 angular sectors meeting at the point i as T0 , . . . , Tki ; They are the inner angles of the triangles which have i as a vertex. At node i, the indexing of the angular sectors is ordered counterclockwise. − → n l+1/2 denotes the unit vector of the half-line Dl+1/2 = Tl ∩ Tl+1 , and l the inner angle of sector Tl , 0 l ki . We will denote by ui the numerical approximation to the viscosity solution of (1) at node i. (∇u)0 , . . . , (∇u)ki will, respectively, represent the numerical approximation of ∇u at node i in each angular sector T0 , . . . , Tki . An important building block for the algorithms in this paper is the Lax–Friedrichs-type monotone Hamiltonian for arbitrary triangulations developed by Abgrall in [1]. It is given by ((∇u)0 , . . . , (∇u)k ) = H H i
k
l=0 l (∇u)l i
2
−
ki
l=0
l+1/2
(∇u)l + (∇u)l+1 2
→ ·− n l+1/2 , (2)
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
185
Fig. 2. Triangle l and its neighbor.
where l+1/2 = tan(l /2) + tan(l+1 /2), = max{max |H1 (u, v)|, max |H2 (u, v)|}. Here H1 and H2 are the partial derivatives of H with respect to ux and uy , respectively, or the Lipschitz constants of H with respect to ux and uy , if H is not differentiable. in (2) defines a monotone Hamiltonian. It is Lipschitz continuous in all arguments and is The H (∇u, . . . , ∇u) = H (∇u). Hence if we have high-order approximations to ∇u consistent with H, i.e., H will be a high-order approximation to H. at node i in every angular sector, the numerical Hamiltonian H Using the Lax–Friedrichs Hamiltonian (2), we can write (1) into the following semi-discrete form: d ((∇u)0 , . . . , (∇u)k ) = 0. ui (t) + H i dt
(3)
For time evolution, TVD Runge–Kutta scheme of Shu and Osher [9] is adopted. This article focuses for (3). primarily on designing an appropriate approximation of ∇u which can fetch out a gratifying H First, we use the central stencil to construct a second-order interpolation polynomial on each triangle. Namely, for triangle l , see Fig. 2, we interpolate on the six nodes: Ai , Bi (i = 1, 2, 3). By solving a linear system, we obtain v l (x, y) = a1 (x − xcl )2 + a2 (x − xcl )(y − ycl ) + a3 (y − ycl )2 + a4 (x − xcl ) + a5 (y − ycl ) + a6 ,
(4)
where Cl = (xcl , ycl ) is the barycenter of l . Because central stencil lacks the capability of treating with singularity, spurious oscillations will occur if we use v l (x, y) directly in (3). Next, we will introduce our novel reconstruction procedure which is the key component of our scheme.
186
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
2.2. The quadratic reconstruction In this section, we will construct a new polynomial pl (x, y) in each triangular which is a modification of v l (x, y) obtained in Section 2.1. In order to exclude oscillation, the partial derivatives of p l (x, y) are subjected to a so-called “limiter”. Due to the operation of this limiter, the slope of pl (x, y) satisfies the maximum principle (as proven in [10]). In smooth region, pl (x, y) is almost the same as v l (x, y). While at the discontinuous points of ∇u, pl (x, y) is modified to avoid new extremum. In the following, we will describe the reconstruction of pl (x, y) in detail. For every triangle l , let N(l ) denote the set of barycenters of triangles which have at least one point meeting with l . Define Mx l =
max
Cj ∈N (j )
(vxl (Cj )),
mx l =
min
Cj ∈N (j )
(vxl (Cj )).
From (4), the x-directional derivative of v l (x, y) is vxl = 2a1 (x − xcl ) + a2 (y − ycl ) + a4
on l .
(5)
To guarantee that the x-directional derivative of pl (x, y) does not suffer from new extremum, we construct a x-directional derivative limiter x for pl (x, y). Thus we obtain pxl (x, y) = x · (2a1 (x − xcl ) + a2 (y − ycl )) + a4 where x = min(x (Ai ))(i = 1, 2, 3). ⎧ Mx l − vxl (Cl ) ⎪ ⎪ min 1, l , ⎪ ⎪ vx (Ai ) − vxl (Cl ) ⎨ mx l − vxl (Cl ) x (Ai ) = ⎪ , min 1, l ⎪ ⎪ vx (Ai ) − vxl (Cl ) ⎪ ⎩ 1,
on l ,
(6)
vxl (Ai ) > vxl (Cl ), vxl (Ai ) < vxl (Cl ),
(7)
vxl (Ai ) = vxl (Cl ).
It is obvious that pxl (x, y) satisfies the maximum principle pxl (x, y) ∈ [mx l , Mx l ]
on l .
(8)
In the same way, We can define My l and my l , and construct a y-directional derivative limiter y to modify vyl (x, y). Hence, the y-directional derivative of pl (x, y) is given by pyl (x, y) = y · (a2 (x − xcl ) + 2a3 (y − ycl )) + a5
on l
(9)
and pyl (x, y) satisfies the local maximum principle too, i.e., pyl (x, y) ∈ [my l , My l ]
on l .
(10)
Gathering from (4), (6) and (9), we obtain pl (x, y) = x · a1 (x − xcl )2 + m · a2 (x − xcl )(y − ycl ) + y · a3 (y − ycl )2 + a4 (x − xcl ) + a5 (y − ycl ) + a6 , where m = min(x , y ).
(11)
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
187
Remark. From (6) and (9), we get two limiters (x and y ) of coefficient a2 in (4). Here we choose the smaller one for (11) to prevent new extremum of p l (x, y)’s two partial derivatives. So far, we get our final approximation p l (x, y) for u and (3) can be written as d ((∇p 0 ), . . . , (∇p k )). ui (t) = L(ui ) = −H dt
(12)
From (8) and (10), we can see that p l (x, y) satisfies the maximum principle; on the other hand, p l (x, y) has no numerical accuracy loss compared with v l (x, y), so our scheme is a second-order TVD scheme. Besides a second-order central stencil interpolation, our reconstruction procedure only involves simple modification, so it is very fast. 2.3. Temporal discretization To keep the same numerical accuracy of spatial discretization, we use the second-order TVD Runge–Kutta scheme of Shu and Osher [9] for the temporal discretization. The algorithm is given by (1)
(0)
(0)
ui = ui + t L(ui ), (2)
(0)
(13)
(1)
(1)
ui = 21 ui + 21 ui + 21 t L(ui )
(14)
with u(0) = un and un+1 = u(2) . 3. Numerical examples In this section, we apply the new scheme developed in Section 2 to some two-dimensional model problems. The triangular mesh is shown in Fig. 3. Let N be the total number of calculational grids. Example 1. Two-dimensional Burgers equation: ut +
(u2x + u2y + 1)2 2
u(x, y, 0) = − cos
= 0,
−2 x 2, −2 y 2,
(x + y)
2
(15)
with periodic boundary conditions. We take = 6.0 and N = 512 in (2) and show the results in Fig. 4. We can see the solution develops a discontinuous derivative when t = 1.5/2 . Example 2. Two-dimensional burning problem: ut − u2x + u2y + 1 = 0, 0 x 1, 0 y 1, u(x, y, 0) = cos(2x) − cos(2y).
(16)
188
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
U
Fig. 3. The calculational triangular meshes.
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 2
-2 -1
1
0 X
0 -1
-2
1 2
Y
Fig. 4. The two-dimensional Burgers equation, H (u, v) = (u2 + u2 + 1)2 /2. Top: t = 0; bottom: t = 1.5/2 .
We take =
√
2/2 and N = 512 in (2). The results are shown in Figs. 5 and 6 .
Example 3. Two-dimensional eikonal equation with a nonconvex Hamiltonian, which arises in geometric optics: ut + u2x + u2y + 1 = 0, 0 x 1, 0 y 1, u(x, y, 0) = 0.25(cos(2x) − 1)(cos(2y) − 1) − 1. √ We take = 2/2 and N = 2048 in (2). The results are shown in Fig. 7.
(17)
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
189
2 1.5 1 0.5 U
0 -0.5 -1 -1.5 1 0.5 Y
0
0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 X
0
2 1.5 1 U
0.5 0 -0.5 -1
Y
-1.5 1 0.8 0.6 0.4 0.2 0
0
0.9 1 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 X
Fig. 5. The burning problem, H (u, v) = − u2 + v 2 + 1. Top: t = 0.5/2 ; bottom: t = 1.5/2 . 16
16
14
14
12
12
10
10
8
8
6
6
4
4
2
2 2
4
6
8
10
12
14
16
2
4
6
8
Fig. 6. Isoline of the burning problem. Left: t = 0.5/2 ; right: t = 1.5/2 .
10
12
14
16
190
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191 0 -0.2 T he
U
-0.4 -0.6 -0.8 -1 1 0.8 0.6 0.4 0.2
Y
0
0
0.2
0.4
0.6
0.8
1
X
-1.4
U
-1.45 -1.5 -1.55 -1.6 1 0.8 0.6 0.4 Y
0.2 0
0
0.2
0.4
0.6
0.8
1
X
Fig. 7. Two-dimensional eikonal equation, H (u, v) = u2 + v 2 + 1. Top: t = 0; bottom: t = 0.6.
4. Concluding remark We developed a TVD-type difference scheme on unstructured triangular meshes for 2D scalar H–J equations. With its solution’s slope satisfying maximum principle, the scheme has two-order numerical accuracy in smooth region. Besides, no parameters have to be turned in our method. Numerical experiments show that the new method has high resolution of derivative singularities and is significantly simpler than WENO. The method should be useful if geometry or adaptivity requires unstructured meshes for the solution of H–J equations. References [1] R. Abgrall, Numerical discretization of the first-order Hamilton–Jacobi equation on triangular meshes, Comm. Pure Appl. Math. 49 (1996) 1339–1373. [2] S. Augoula, R. Abgrall, High order numerical discretization for Hamilton–Jacobi equations on triangular meshes, J. Sci. Comput. 15 (2000) 197–229. [3] M. Crandall, P. Lions, Viscosity solutions of Hamilton–Jacobi equations, Trans. Amer. Math. Section 277 (2) (1983) 1–42. [4] G. Jiang, D.P. Peng, Weighted ENO schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (2000) 2126–2143.
L. Tang, S. Song / Journal of Computational and Applied Mathematics 195 (2006) 182 – 191
191
[5] S.N. Kruzkov, The Cauchy problem in the large for nonlinear equations and for certain quasilinear systems of the first order with several variables, Soviet Math. Dokl. 5 (1964) 493–496. [6] C.T. Lin, E. Tadmor, High-resolution non-oscillatory central schemes for approximate Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (2000) 2163–2186. [7] S. Osher, J. Sethian, Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [8] S. Osher, C.W. Shu, High-order essentially schemes for Hamilton–Jacobi equations, SIAM J. Numer. Anal. 28 (1991) 907–922. [9] C. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys. 77 (1988) 439–471. [10] S. Song, Y. Li, A class of finite volume scheme satisfying maximum principle for 2D scalar hyperbolic conservation laws on unstructured meshes, Chinese J. Numer. Math. Appl. 19 (4) (1997). [11] Y. Zhang, C. Shu, High order WENO schemes for Hamilton–Jacobi equations on triangular meshes, SIAM J. Sci. Comput. 24 (2003) 1005–1030.