Local discontinuous Galerkin method for the Keller ... - Brown University
Local discontinuous Galerkin method for the Keller-Segel chemotaxis model Xingjie Helen Li1 , Chi-Wang Shu2 and Yang Yang3 Abstract In this paper, we apply the local discontinuous Galerkin (LDG) method to 2D Keller– Segel (KS) chemotaxis model. We improve the results upon (Y. Epshteyn and A. Kurganov, SIAM Journal on Numerical Analysis, 47 (2008), 368-408) and give optimal rate of convergence under special finite element spaces. Moreover, to construct physically relevant numerical approximations, we develop a positivity-preserving limiter to the scheme, extending the idea in (Y. Zhang, X. Zhang and C.-W. Shu, Journal of Computational Physics, 229 (2010), 8918-8934). With this limiter, we can prove the L1 -stability of the numerical scheme. Numerical experiments are performed to demonstrate the good performance of the positivity-preserving LDG scheme. Moreover, it is known that the chemotaxis model will yield blow-up solutions under certain initial conditions. We numerically demonstrate how to find the numerical blow-up time by using the L2 -norm of the L1 -stable numerical approximations.
Keywords: Local discontinuous Galerkin method, Keller-Segel chemotaxis model, positivity preserving, error estimate, Neumann boundary condition, blow-up, L1 stability 1
Department of Mathematics and Statistics, University of North Carolina at Charlotte, Charlotte, NC 28223. E-mail: [email protected] 2 Division of Applied Mathematics, Brown University, Providence, RI 02912. E-mail: [email protected]. Research supported by ARO grant W911NF-15-1-0226 and NSF grant DMS1418750. 3 Department of Mathematical Sciences, Michigan Technological University, Houghton, MI 49931. E-mail: [email protected]
1
1
Introduction
In this paper, we study the Keller-Segel (KS) chemotaxis model in two space dimensions [31, 27] and focus on the following common formulation [3], ut − div(∇u − χu∇v) = 0, vt − △v = u − v,
x ∈ Ω, t > 0, x ∈ Ω, t > 0,
(1.1)
where Ω is assumed to be a convex, bounded and open set in R2 . Chemotaxis is the highly nonlinear term which indicates movements by cells in reaction to a chemical substance, where cells approach chemically favorable environments and avoid unpleasant ones. In (1.1), u and v denote the densities of cells and the chemical concentration, respectively. The chemotactic sensitivity function χ is assumed to be a positive constant. For simplicity, we take χ ≡ 1. The boundary conditions are set to be ∇u · n = ∇v · n = 0, where n is the outer normal of the boundary ∂Ω. With this boundary condition, constant during the time evolution and the system is thus isolated.
(1.2) R
Ω
u is a
The exact solutions of the KS chemotaxis model are always positive. Moreover, the model exhibits blow-up patterns with certain initial conditions [18, 17]. When the blow-up patterns occur, the density u of cells will strengthen in the neighborhood of isolated points, and these regions become sharper and eventually result in finite time point-wise blow-up. Mathematical proofs for spherically symmetric solutions in a ball have been investigated in [23]. The blow-up point in this case is the center of the ball, and it is proved to be the only possible singularity. For nonsymmetric cases, blow-up results have also been studied in [24, 25], and the blow-up points will be located on the boundary. More theoretical works can be found in [28, 18]. It is difficult to construct numerical schemes for (1.1), and most of the previous works are for the following simplified system ut − div(∇u − χu∇v) = 0, −△v = u − v, 2
x ∈ Ω, t > 0, x ∈ Ω, t > 0,
(See, for example, [38, 29, 16, 22] and the references therein). Recently, there are some significant works designed to solve (1.1) directly [15, 37]. In [37], the authors constructed implicit second-order positivity preserving finite-volume schemes in three-dimensional space, and their technique requires solving a large linear system of equations coupling together all grid points at each stage of the two stage TR-BDF2 method when updating the diffusion terms at each time step. In [15], the authors applied the interior penalty discontinuous Galerkin (IPDG) method on rectangular meshes to obtain suboptimal rate of convergence, and the finite element space is assumed to be piecewise polynomials of degree k ≥ 2. Other related works in this direction include [14, 12, 13] and the positivity-preserving property was demonstrated by numerical experiments only. Besides the above, in [2], the authors constructed a second order positivity-preserving scheme to a revised system by differentiating (1.1) with respect to x and y. Hence, the scheme was not designed to solve (1.1) directly, similar work can also be found in [34]. Recently, Zhang and Shu introduced positivity-preserving limiter to hyperbolic equations in [45, 46, 47]. Subsequently, the idea was extended to parabolic problems in [48]. In this paper, we will apply the positivity-preserving limiter introduced in [48] to construct second-order positivity-preserving local discontinuous Galerkin (LDG) schemes to obtain physically relevant numerical approximations. The method we plan to use preserves the positivity of the numerical solutions, and can be applied to unstructured meshes. The DG method was first introduced in 1973 by Reed and Hill [33] in the framework of neutron linear transport. Subsequently, Cockburn et al. developed Runge-Kutta discontinuous Galerkin (RKDG) methods for hyperbolic conservation laws in a series of papers [9, 6, 8, 10]. In [11], Cockburn and Shu introduced the LDG method to solve the convectiondiffusion equations. Their idea was motivated by Bassi and Rebay [1], where the compressible Navier-Stokes equations were successfully solved. Recently, the DG methods were applied to linear hyperbolic equations with δ-singularities [41] to obtain high-order approximations under suitable negative-order norms. Subsequently, the methods have also been applied to
3
nonlinear hyperbolic equations with δ-singularities [42, 49]. Combined with special limiters, the schemes were proved to be L1 stable [42, 32]. Recently, the idea has been extended to parabolic equations with blow-up solutions by using the LDG method [21]. In this paper, we follow the same direction and employ the LDG method to capture the blow-up phenomenon. In the LDG method, we introduce auxiliary variables p = ux and q = uy . Numerical experiments will be given to demonstrate the optimal rate of convergence. However, in this problem the approximations of p and q are discontinuous across the cell interfaces and we can hardly obtain error estimates if we analyze the convection and diffusion terms separately. To explain this point, let us consider the following hyperbolic equation ut + (a(x)u)x = 0, where a(x) is discontinuous at x = x0 . In [19, 26], the authors studied such a problem and defined Q=
a(x0 + b) − a(x0 ) . b
If Q is bounded from below for all b, then the solution exists, but may not be unique. If Q is bounded from above for all b, we can guarantee the uniqueness, but the solution may not exist. To bridge this gap, most of the previous works for similar problems are based on mixed finite element methods (see e.g. [29]). In this paper, we consider a new technique for error analysis following the idea in [39, 40]. Moreover, we also apply the positivity-preserving limiter to guarantee positivity of the numerical approximation. With this special technique, the L1 norm of the numerical approximation is a constant during the time evolution due to the mass conservation. However, the L2 norm might still be unbounded when blow-up patterns occur. We thus introduce a new idea to capture the numerical blow-up time based on the L2 norm of the numerical approximations. To our best knowledge, this is the first paper studying the numerical blow-up time with a concrete theoretical background. The organization of this paper is as follows. In Section 2, we construct the LDG scheme. In Section 3, we give the error estimates based on two different finite element spaces. In 4
Section 4, we discuss the positivity-preserving technique, the foundation of limiters and high order time discretizations. Numerical experiments are given in Section 5. Finally, we will end in Section 6 with concluding remarks and remarks for future work.
2
LDG scheme
In this section, we define the finite element spaces and proceed to construct the LDG scheme. In this paper, we consider the computational domain to be Ω = [0, 1] × [0, 1] and use rectangular meshes for simplicity. The analyses for triangular meshes can be extended with some minor changes, so we will skip this part. Let Kij = (xi− 1 , xi+ 1 ) × (yj− 1 , yj+ 1 ), i = 2
2
2
2
1, . . . , Nx , j = 1, . . . , Ny , be a partition of the computational domain Ω, and the mesh sizes in the x and y directions are given as ∆xi = xi+ 1 −xi− 1 and ∆yj = yj+ 1 −yj− 1 , respectively. For 2
2
2
2
simplicity, in this paper, we assume uniform meshes and denote h = ∆xi = ∆yj . However, this assumption is not essential in the proof. Moreover, we also denote Ii = [xi− 1 , xi+ 1 ] and 2
2
Jj = [yj− 1 , yj+ 1 ] to be the cells in x and y directions. 2
2
We define the finite element spaces Vhk as n o Vhk = z : z Kij ∈ P k (Kij ) ,
where P k (Kij ) denotes the set of polynomials of degree up to k in cell Kij . To construct the LDG scheme, we introduce the axillary variables p = ux , q = uy , r = vx and s = vy , then (1.1) can be written as ut = −(ru)x − (su)y + px + qy , p = ux , q = uy , vt = rx + sy + u − v, r = vx , s = vy . 5
For simplicity, we use notations o to represent the exact solutions of the six variables u, p, q, v, r, s and oh for the corresponding numerical approximations. Then the LDG scheme is to find uh , ph , qh ∈ Vhk1 and vh , rh , sh ∈ Vhk2 , such that for any test functions w u , w p, w q ∈ Vhk1 and w v , w r , w s ∈ Vhk2 u
(uh t , w )ij =
(rh uh , wxu )ij
+
Z
Jj
Z
u + rd h uh (w )i− 1 ,j dy 2
−
Z
Z
Jj
u − rd h uh (w )i+ 1 ,j dy 2
u − u + sd + (sh uh , wyu )ij + sd h uh (w )i,j+ 1 dx h uh (w )i,j− 1 dx − 2 2 Ii Z Ii Z − (ph , wxu )ij − pbh (w u )+ dy + pbh (w u )− dy i− 12 ,j i+ 21 ,j Jj Jj Z Z u u + qbh (w u )− dx, − (qh , wy )ij − qbh (w )i,j− 1 dx + i,j+ 21 2 Ii Ii Z Z p p p + uh (w p )− c dy, (ph , w )ij = −(uh , wx )ij − uh (w )i− 1 ,j dy + c i+ 21 ,j 2 Jj Jj Z Z q q q + uh (w q )− c dx, (qh , w )ij = −(uh , wy )ij − uh (w )i,j− 1 dx + c i,j+ 21 2 Ii Ii Z Z v v v + rbh (w v )− dy rbh (w )i− 1 ,j dy + (vht , w )ij = −(rh , wx )ij − i+ 21 ,j 2 Jj Jj Z Z v v + − (sh , wy )ij − sbh (w )i,j− 1 dx + sbh (w v )− dx i,j+ 1 2
Ii
2
where (u, v)ij = for
− wi+ 1 , ,j 2
+ wi,j− 1 2
R
Kij
and
(2.2) (2.3)
2
Ii
+ (uh − vh , w v )ij , Z Z r + r r vbh (w r )− vbh (w )i− 1 ,j dy + dy, (rh , w )ij = −(vh , wx )ij − i+ 12 ,j 2 Jj Jj Z Z s s s + (sh , w )ij = −(vh , wy )ij − vbh (w )i,j− 1 dx + vbh (w s )− dx, i,j+ 1 Ii
(2.1)
(2.4) (2.5) (2.6)
2
Ii
+ is the trace of w at x = xi− 1 in cell Kij . Likewise uvdxdy, and wi− 1 ,j 2
2
− wi,j+ 1. 2
We also denote rd h uh , sd h uh and obh to be the numerical fluxes.
In this paper, we use the alternating fluxes for the diffusion terms and Lax-Friedrichs fluxes for the convection terms. For simplicity, if we consider the generic trace at the interior cell interfaces, the subscript will be omitted. Due to the Neumann boundary condition (1.2), we take
6
1. For j = 1, 2, · · · , Ny , + , i = 0, u 1 ,0 h i+ 2 ,j − uh + , i = 1, · · · , Nx − 1, u ch i+ 1 ,j , pbh i+ 1 ,j = 1 , ph 1 i+ ,j 2 2 2 i+ 2 ,j uh − 1 , 0 , i = Nx , i+ ,j
(2.7)
2
+ , 0 , i = 0, v h i+ 1 ,j 2 − vh + , i = 1, · · · , Nx − 1, vbh i+ 1 ,j , rbh i+ 1 ,j = 1 , rh 1 i+ ,j 2 2 2 i+ 2 ,j vh − 1 , 0 , i = Nx . i+ ,j
(2.8)
2
2. For i = 1, 2, · · · , Nx
+ , 0 , j = 0, u h 1 i+ 2 ,j − uh + , j = 1, · · · , Ny − 1, uh i,j+ 1 , qbh i,j+ 1 = c 1 , qh 1 i,j+ 2 2 2 i,j+ 2 uh − 1 , 0 , j = Ny , i,j+
(2.9)
2
+ v , 0 , j = 0, h i+ 1 ,j 2 , sh − , j = 1, · · · , Ny − 1, vh + vbh i,j+ 1 , sbh i,j+ 1 = 1 i,j+ 21 2 2 i,j+ 2 vh − 1 , 0 , j = Ny . i,j+
(2.10)
2
For the convection term, we consider Lax-Friedrichs flux rd h uh =
1 + + + − rh uh + rh− u− h − α(uh − uh ) , 2
(2.11)
where α is chosen by the positivity-preserving technique. Due to the mass conservation, we also take rd h uh Nx + 1 ,j = 0, h uh 1 ,j = rd 2
2
j = 1, · · · , Ny .
(2.12)
We also define sd h uh is a similar way.
Also let [w] = w + − w − denote the jump of w at the interface of each cell and (u, v) =
Ny Nx X X
(u, v)ij ,
i=1 j=1
then (2.1)-(2.6) can be written as (uht , w u ) = Lcx (rh , uh , w u ) + Lcy (sh , uh , w u ) − Ldx (ph , w u ) − Ldy (qh , w u ), 7
(2.13)
(ph , w p ) = −Qx (uh , w p ),
(2.14)
(qh , w q ) = −Qy (uh , w q ),
(2.15)
(vht , w v ) = −Ldx (rh , w v ) − Ldy (sh , w v ) + (uh − vh , w v ),
(2.16)
(rh , w r ) = −Qx (vh , w r ),
(2.17)
(sh , w s ) = −Qy (vh , w s ),
(2.18)
where Lcx (r, u, w)
Ny Z y Nx −1 X 1 j+ 2 1 X r + u+ + r − u− − α[u] [w]i+ 1 ,j dy, = (ru, wx ) + 2 2 i=1 j=1 y 1
(2.19)
j− 2
Lcy (s, u, w) = (su, wy ) + Ldx (p, w) = (p, wx ) +
1 2
y −1 Z x Nx N X X i+ 1
2
i=1 j=1
2
Q (u, w) = (u, wx ) +
yj− 1
2
xi− 1
q − [w]i,j+ 1 dx,
(2.22)
2
2
u [w]i+ 1 ,j dy + 2
Ny Z y N x −1 X X j+ 1 2
yj− 1
xi− 1
y −1 Z Nx N X X
i=1 j=1
yj+ 1 2
u− [w]Nx + 1 ,j dy (2.23) 2
yj− 1
2
− [u]wi+ 1 dy, ,j
(2.24)
2
2
2
i=1 j=0
Ny Z X j=1
2
y −1 Z x Nx N X X i+ 1
= −(uy , w) −
+
yj− 1
i=1 j=1
Q (u, w) = (u, wy ) +
(2.21)
2
2
i=0 j=1
y
p− [w]i+ 1 ,j dy,
Ny Z y N x −1 X X j+ 1
= −(ux , w) −
(2.20)
2
2
y −1 Z x Nx N X X i+ 1
i=1 j=1
x
2
Ny Z y 1 N x −1 X X j+ i=1 j=1
Ldy (q, w) = (q, wy ) +
xi− 1
s+ u+ + s− u− − α[u] [w]i,j+ 1 dx,
+
u [w]i,j+ 1 dx + 2
Nx Z X i=1
2
xi+ 1
2
xi+ 1
2
xi− 1
u− [w]i,Ny + 1 dx (2.25) 2
2
− [u]wi,j+ 1 dx.
(2.26)
2
xi− 1 2
+ + Here, for simplicity, we take w −1 ,j = wN = wi,−1 = wi,N 1 1 = 0. It is easy to check the x + ,j y+ 2
2
2
2
following identities by integration by parts on each cell: for any functions u and w, Ldx (u, w) + Qx (w, u) = 0,
(2.27)
Ldy (u, w) + Qy (w, u) = 0.
(2.28)
8
3
Error estimates
In this section, we analyze the error between the numerical and exact solutions. We first introduce the norms, construct the error equations and state the main theorem. Then we proceed to the proof. The whole procedure can be divided into several steps. We first construct the special projections that will be used in this section. Then we give an a priori error estimate and its verification. Finally, we will consider another finite element space and the corresponding error estimate. Let us define the norms first.
3.1
Norms
In this subsection ,we define several norms that will be used throughout the paper. Denote kuk0,Kij to be the standard L2 norm of u in cell Kij . For any natural number ℓ, we consider the norm of the Sobolev space H ℓ (Kij ), defined by ) 12 (
X
∂ α+β u 2
. kukℓ,Kij =
∂ α x∂ β y 0,Kij 0≤α+β≤ℓ
Moreover, we define the norms on the whole computational domain as ! 12 Ny Nx X X kukℓ = . kuk2ℓ,Kij i=1 j=1
For convenience, if we consider the standard L2 norm, then the corresponding subscript will be omitted. Let Γij be the edges of Kij , and we define Z Z 2 − 2 2 − 2 2 + (u+ dx. kukΓij = (ui− 1 ,j ) + (ui+ 1 ,j ) dy + 1 ) + (u 1) i,j− i,j+ Jj
2
2
Ii
2
2
Moreover, we also define
kuk2Γh
=
Ny Nx X X
kuk2Γij .
i=1 j=1
Finally, we define the standard L∞ norm of u in Kij as kuk∞,Kij , and define the L∞ norm on the whole computational domain as kuk∞ = max kuk∞,Kij . i,j
9
3.2
Error equations and the main theorem
In this subsection, we proceed to construct the error equations. Denote the error between the exact and numerical solutions to be eo = o − oh , then we have the equations satisfied by the errors as (eut , w u ) =Lcx (r, u, w u) − Lcx (rh , uh , w u ) + Lcy (s, u, w u) − Lyc (sh , uh , w u ) − Ldx (ep , w u ) − Ldy (eq , w u ),
(3.1)
(ep , w p ) = − Qx (eu , w p ),
(3.2)
(eq , w q ) = − Qy (eu , w q ),
(3.3)
(ev t , w v ) = − Ldx (er , w v ) − Ldy (es , w v ) + (eu − ev , w v ),
(3.4)
(er , w r ) = − Qx (ev , w r ),
(3.5)
(es , w s ) = − Qy (ev , w s).
(3.6)
Now, we can state our main theorem. Theorem 3.1. Suppose u and v are smooth functions, r and s are uniformly bounded for t ≤ T . The numerical approximations uh , ph , qh ∈ Vhk1 and vh , sh , rh ∈ Vhk2 . The initial discretization is given as the standard L2 -projection. If we take k1 ≥ 1 and k2 ≥ 2, then we have ku − uh k + kv − vh k ≤ Chmin{k1 +1,k2 } , where the positive constant C does not depend on h.
3.3
Preliminaries and projections
In this subsection, we study the basic properties of the finite element space. Let us start with the classical inverse properties. Lemma 3.1. Assuming u ∈ Vhk , there exists a positive constant C independent of h and u such that hkuk∞,Kij + h1/2 kukΓij ≤ CkukKij . 10
In this paper, we consider several special projections. We define the elliptic projections Pu, Pp, Pq ∈ Vhk1 such that for any w u , w p, w q ∈ Vhk1 Ldx (p, w u ) + Ldy (q, w u) =Ldx (Pp, w u ) + Ldy (Pq, w u ),
(3.7)
(Pp, w p ) + (Pq, w q ) = − Qx (Pu, w p ) − Qy (Pu, w q ).
(3.8)
Similarly, we also define Pv, Pr, Ps ∈ Vhk2 such that for any w v , w r , w s ∈ Vhk2 Ldx (r, w v ) + Ldy (s, w v ) =Ldx (Pr, w v ) + Ldy (Ps, w v ),
(3.9)
(Pr, w r ) + (Ps, w s ) = − Qx (Pv, w r ) − Qy (Pv, w s ).
(3.10)
The existence of the elliptic projections will be given in the following lemma [5], Lemma 3.2. Suppose
R
Pudxdy = Ω
R
udxdy and Ω
R
Pvdxdy = Ω
R
Ω
vdxdy, then Po are
uniquely determined by (3.7)-(3.10). Moreover, we have the following estimates kp − Ppk + kq − Pqk ≤Chk1 , ku − Puk ≤Chk1 +1 , kr − Prk + ks − Psk ≤Chk2 , kv − Pvk ≤Chk2 +1 ,
(3.11) (3.12) (3.13) (3.14)
where C is independent of h. In addition, we also define the standard L2 projection P by (P u, v)ij = (u, v)ij ,
∀v ∈ P k (Kij ).
The L2 projection satisfies the following property [4]. Lemma 3.3. Suppose u ∈ C k+1 (Ω). Then there exists a positive constant C independent of h and u such that ku − P uk + hku − P uk∞ + h1/2 ku − P ukΓh ≤ Chk+1 kukk+1. 11
Let us finish this section by proving the following lemma. Lemma 3.4. Let u ∈ C k+1(Ω) and Πu ∈ Vhk . Suppose ku − Πuk ≤ Chκ for some positive constant C and κ ≤ k + 1. Then hku − Πuk∞ + h1/2 ku − ΠukΓh ≤ Chκ , where the positive constant C does not depend on h. Proof. Actually, hku − Πuk∞ + h1/2 ku − ΠukΓh ≤ hkP u − Πuk∞ + hku − P uk∞ + h1/2 kP u − ΠukΓh + h1/2 ku − P ukΓh ≤ CkP u − Πuk + Chk+1 + CkP u − Πuk + Chk+1 ≤ Cku − Πuk + Cku − P uk + Chk+1 ≤ Chκ , where the first and third steps follow from triangle inequality, the second step requires Lemma 3.1 and Lemma 3.3, and the last step we use Lemma 3.3.
3.4
A priori error estimate
In this subsection, we would like to make an a priori error estimate that ku − uh k ≤ h,
(3.15)
which further implies ku − uh k∞ ≤ C by Lemma 3.4. Moreover, if u is bounded, then kuh k∞ ≤ C. This assumption, which will be verified in Section 3.6, is used for the estimate of the convection terms. This idea has been used to obtain the error estimates for nonlinear hyperbolic equations [43, 44, 30]. With this assumption, we can proceed to the main proof of the theorem.
12
3.5
Proof of Theorem 3.1
In this subsection, we give the proof of Theorem 3.1. As the general treatment of the finite element methods, we divide the error as eo = ηo − ξo , where ηo = o − Po and ξo = oh − Po . Clearly, ξo is chosen from the desired finite element space, and following [39, 40] we have Lemma 3.5. k(ξu )x k2 + k(ξu )y k2 + h−1 k[ξu ]k2Γh ≤ C(kξp k2 + kξq k2 ). Proof. By (2.14), (2.15) and (3.8), we have (ξp , w p ) + (ξq , w q ) = −Qx (ξu , w p ) − Qy (ξu , w q ) We take w p (x, y) =
xi+ 1 −x 2
∆xi
(3.16)
(ξu )x , (x, y) ∈ Kij and w q = 0, then by (2.24) and (2.14) it is
easy to show |||(ξu )x |||2Kij where
xi+ 1 − x 2 (ξu )x ≤ k(ξu )x kKij kξp kKij . = ξp , ∆xi ij |||u|||2Kij
=
Z
xi+ 1 − x 2
∆xi
Kij
u2 dxdy.
2
Clearly, ||| · |||Kij is a weighted L -norm in P k1 (Kij ), and is equivalent to the standard L2 -norm, i.e. for any u ∈ P k1 (Kij ), kukKij ≤ C|||u|||Kij , where the constant C does not depend on Kij . Therefore, k(ξu )x k2Kij ≤ C|||(ξu )x |||2Kij ≤ Ck(ξu )x kKij kξp kKij , which further yields k(ξu )x kKij ≤ Ckξp kKij . Similarly, we can also prove k(ξu )y k ≤ Ckξq k 13
following the same line. Next, we will take care of interfacial terms. We take w p (x, y) = ξu (x, y) − ξu+ (xi+ 1 , y), 2
(x, y) ∈ Kij and w q = 0 in (3.16) to obtain in cell Kij , Z
Jj
[ξu (xi+ 1 , y)]2 dy = (ξp − (ξu )x , w p )ij ≤ Ckw p kKij kξp kKij . 2
It is easy to see that p
w (x, y) = −[ξu (xi+ 1 , y)] − 2
Z
x+ 12
(ξu )x dx. x
Therefore, by Minkowski’s inequality, Z
[ξu (xi+ 1 , y)] dy ≤ Ckξp kKij [ξu (xi+ 1 , y)] + ∆xi k(ξu )x kKij 2 2 Kij Jj !1/2 Z 1/2 ≤ Ckξp kKij ∆xi [ξu (xi+ 1 , y)]2dy + ∆xi kξp kKij 2
1 ≤ 2
which further yields
Similarly, we have
Jj
Z
Z
Jj
Ii
[ξu (xi+ 1 , y)]2dy + C∆xi kξp k2Kij 2
[ξu (xi+ 1 , y)]2 dy ≤ C∆xi kξp k2Kij .
(3.17)
[ξu (x, yj+ 1 )]2 dx ≤ C∆yj kξq k2Kij .
(3.18)
Jj
Z
2
2
2
Hence, h−1 k[ξu ]k2Γh ≤ C(kξp k2 + kξq k2 ).
With the above lemma, we can proceed to the proof of the main theorem. We take w o = ξo in (3.1)-(3.6) to obtain 1 dkξu k2 + kξp k2 + kξq k2 = T1 − T2 − T3 , 2 dt 1 dkξv k2 + kξr k2 + kξs k2 = T4 , 2 dt
14
(3.19) (3.20)
where T1 = ((ηu )t , ξu ), T2 = (ru − rh uh , (ξu )x ) + hru − rd h uh , [ξu ]ix , T3 = (su − sh uh , (ξu )y ) + hsu − sd h uh , [ξu ]iy , T4 = ((ηv )t , ξv ) − (ηu − ξu − ηv + ξv , ξv ),
with hu, vix =
Ny Z N x −1 X X i=1 j=1
Jj
uvi+ 1 ,j dy 2
and hu, viy =
y −1 Z Nx N X X
i=1 j=1
Furthermore, we also define
kukΓxh = hu, uix
Ii
uvi,j+ 1 dx. 2
and kukΓyh = hu, uiy .
Then by (3.18) and (3.17), we have k[ξu ]kΓxh ≤ Ch1/2 kξp k and k[ξu ]kΓyh ≤ Ch1/2 kξq k.
(3.21)
Now we give the estimates of Ti′ s. By Cauchy-Schwarz inequality and Lemma 3.2 T1 ≤ k(ηu )t kkξu k ≤ Chk1 +1 kξu k.
(3.22)
The estimates of T2 and T3 are similar, and we consider T2 only. T2 = (r(u − uh ) + (r − rh )uh , (ξu )x ) + hru − rd h uh , [ξu ]ix , = T21 + T22 + T23 ,
where T21 = (r(u − uh ) + (r − rh )uh , (ξu )x ) T22 =
1 − − h2ru − rh+ u+ h − rh uh , [ξu ]ix 2
T23 = −hα[ηu − ξu ], [ξu ]ix , Using the fact that krk∞ ≤ C and the a priori error estimate kuh k∞ ≤ C, we have T21 ≤ C (ku − uh k + kr − rh k) k(ξu )x k 15
(3.23)
≤ C(kηu k + kξu k + kηr k + kξr k) kξp k ≤ C(hk1 +1 + kξu k + hk2 + kξr k) kξp k,
(3.24)
where the first step is based on Cauchy-Schwarz inequality, the second step follows from Lemma 3.5 and triangle inequality, and in the last one we use Lemma 3.2. The estimate of T22 also requires the boundedness of r and uh , T22 =
1 + + − − − hr(u − u+ h ) + (r − rh )uh + r(u − uh ) + (r − rh )uh , [ξu ]ix 2
≤ C (ku − uh kΓh + kr − rh kΓh ) k[ξu ]kΓxh ≤ Ch1/2 (kηu kΓh + kξu kΓh + kηr kΓh + kξr kΓh ) kξp k ≤ C(hk1 +1 + kξu k + hk2 + kξr k) kξp k.
(3.25)
Here in the second step we use Cauchy-Schwarz inequality, the third step follows from triangle inequality and (3.21), the last one requires Lemma 3.4 and Lemma 3.2. Now we proceed to the estimate of T23 , T23 ≤ C(k[ηu ]kΓh + k[ξu ]kΓxh ) k[ξu ]kΓxh ≤ Ch1/2 (k[ηu ]kΓh + k[ξu ]kΓh ) kξp k ≤ C(hk1 +1 + kξu k) kξp k,
(3.26)
where the first step follows from Cauchy-Schwarz inequality, the second step is based on (3.21), the third one requires Lemma 3.4 and Lemma 3.2. Plug (3.24), (3.25) and (3.26) into (3.23) to obtain T2 ≤ (Chk1 +1 + Chk2 + Ckξu k + Ckξr k) kξp k.
(3.27)
Following the same analysis of T2 , we obtain a similar estimate of T3 : T3 ≤ (Chk1 +1 + Chk2 + Ckξu k + Ckξs k) kξq k.
(3.28)
The estimate of T4 is simple, we only use Cauchy-Schwartz inequality and Lemma 3.2, T4 ≤ (k(ηv )t k + kηu k + kξu k + kηv k + kξv k) kξv k 16
≤ C hmin(k1 ,k2 )+1 + kξu k + kξv k kξv k
(3.29)
Plug (3.22), (3.27) and (3.28) into (3.19) to obtain 1 dkξu k2 + kξp k2 + kξq k2 ≤ Chk1 +1 kξu k + C(hmin(k1 +1,k2 ) + kξu k)(kξp k + kξq k) 2 dt +Ckξr kkξp k + Ckξs kkξq k ≤ C h2 min(k1 +1,k2 ) + kξu k2 + kξs k2 + kξr k2 + kξp k2 + kξq k2 , which further implies 1 dkξu k2 ≤ C h2 min(k1 +1,k2 ) + kξu k2 + kξs k2 + kξr k2 . 2 dt
(3.30)
Moreover, plug (3.29) into (3.20) to yield 1 dkξv k2 + kξr k2 + kξs k2 ≤ C hmin(k1 ,k2 )+1 + kξu k + kξv k kξv k. 2 dt
(3.31)
Combining (3.30) and (3.31), we can find a constant γ such that 1 dkξu k2 +γ 2 dt
1 dkξv k2 + kξr k2 + kξs k2 2 dt
≤ C h2 min(k1 +1,k2 ) + kξu k2 + kξs k2 + kξr k2 + γC hmin(k1 ,k2 )+1 + kξu k + kξv k kξv k Take γ to be sufficiently large, we have dkξu k2 dkξv k2 +γ ≤ Ch2 min(k1 +1,k2 ) + Ckξu k2 + Ckξv k2 . dt dt Therefore, by the Gronwall’s inequality and use L2 -projection as the initial discretization, kξu k2 + kξv k2 ≤ Ch2 min(k1 +1,k2 ) , which further yield ku − uh k + kv − vh k ≤ Chmin(k1 +1,k2 ) by Lemma 3.2.
17
3.6
Verification of the a priori error estimate
In this section, we proceed to verify the a priori error estimate (3.15). Actually, (3.15) is true at t = 0. Denote t⋆ = inf{t|k(u − uh )(t)k > h}. If t⋆ < T , then at t = t⋆ , we have h = ku − uh k ≤ Chmin(k1 +1,k2 ) , which is a contradiction if k1 ≥ 1, k2 ≥ 2 and h is small. Therefore, we have ku − uh k ≤ Chmin(k1 +1,k2 ) . Now we finish the whole proof.
3.7
Qk polynomials
In this section, we consider another finite element space and give the error estimates. If not otherwise stated, we follow the same notations used before. We first introduce the finite element space Whk defined as n o Whk = z : z Kij ∈ Qk (Kij ) ,
where Qk (Kij ) denotes the set of tensor product polynomials of degree up to k in cell Kij . In this section, we will take vh , rh , sh ∈ Whk and uh , ph , qh ∈ Vhk to obtain optimal error estimate. We first define several special projections. We define P+ into Whk which is, for each cell Kij , Z k−1 (P+ u − u, v)ij = 0, ∀v ∈ Q (Kij ), (P+ u − u)(xi− 1 , y)v(y)dy = 0, ∀v ∈ P k−1(Jj ), 2 Jj Z (P+ u − u)(xi− 1 , yj− 1 ) = 0, (P+ u − u)(x, yj− 1 )v(x)dx = 0, ∀v ∈ P k−1(Ii ). (3.32) 2
2
2
Ii
Moreover, we also define Πxk and Πyk into Whk which are, for each cell Kij , Z x k (Πk u − u, vx )ij = 0, ∀v ∈ Q (Kij ), (Πxk u − u)(xi+ 1 , y)v(y)dy = 0, ∀v ∈ P k (Jj ), 2 J Zj (3.33) y y k k (Πk u − u, vy )ij = 0, ∀v ∈ Q (Kij ), (Πk u − u)(x, yj+ 1 )v(x)dx = 0, ∀v ∈ P (Ii ). 2
Ii
Following the same analysis before, we define ηv = v − P+ v,
ηr = r − Πxk r,
ηs = s − Πyk s,
ξv = vh − P+ v,
ξr = rh − Πxk r,
ξs = sh − Πyk s.
and
Following the proof of Lemma 3.5, we obtain another similar one 18
Lemma 3.6. Suppose ξv , ξr and ξs are defined above, we have k(ξv )x k ≤ Ckξr k,
k(ξv )y k ≤ Ckξs k,
h−1 k[ξv ]k2Γxh ≤ Ckξr k2 ,
h−1 k[ξv ]k2Γy ≤ Ckξs k2 . h
Moreover, we will use the following lemma [4] Lemma 3.7. Suppose v is sufficiently smooth, then we have kηv k + hk∇ηv k + h1/2 kηv kΓh ≤Chk+1 ,
(3.34)
kηr k + hk∇ηr k + h1/2 kηr kΓh ≤Chk+1 ,
(3.35)
kηs k + hk∇ηs k + h1/2 kηs kΓh ≤Chk+1 .
(3.36)
Finally, we also need the following superconvergence property [7]. Lemma 3.8. For any v ∈ Whk , we have |Qx (ηu , v)| ≤ Chk+1 kukk+2kvk,
|Qy (ηu , v)| ≤ Chk+1 kukk+2kvk.
Now we can state the main theorem. Theorem 3.2. Suppose the exact solutions u, v are smooth and the derivatives r, s are uniformly bounded. The LDG scheme is defined as (2.1)-(2.6) with uh , qh , qh ∈ Vhk and vh , rh , sh ∈ Whk for k ≥ 1. Moreover, the initial discretization is also given as the L2 projection. Then we have ku − uh k + kv − vh k ≤ Chk+1 . Proof. In the proof, we skip the a posterior error estimate, its verification in Section 3.5 and most of the derivation steps, since they can be obtained from Section 3.5 with some minor changes. We take w o = ξo in (3.1)-(3.6) to obtain 1 dkξu k2 + kξp k2 + kξq k2 = T1 − T2 − T3 , 2 dt 1 dkξv k2 + kξr k2 + kξs k2 = T4 + T5 , 2 dt
19
(3.37) (3.38)
where T1 = ((ηu )t , ξu ), T2 = (ru − rh uh , (ξu )x ) + hru − rd h uh , [ξu ]ix , T3 = (su − sh uh , (ξu )y ) + hsu − sd h uh , [ξu ]iy , T4 = ((ηv )t , ξv ) − (ηu − ξu − ηv + ξv , ξv ), T5 = Qx (ηu , ξv ) + Qy (ηu , ξv ) Following the analyses in Section 3.5, we obtain the estimates of Ti′ s. T1 ≤ Chk1 +1 kξu k.
(3.39)
T2 ≤ (Chk+1 + Ckξu k + Ckξr k) kξp k
(3.40)
T3 ≤ (Chk+1 + Ckξu k + Ckξs k) kξq k T4 ≤ C hk+1 + kξu k + kξv k kξv k
(3.41) (3.42)
Finally, the estimate of T5 follows from Lemma 3.8 directly, T5 ≤ Chk+1 (kξr k + kξs k).
(3.43)
Plug (3.39), (3.40) and (3.41) into (3.37) to obtain 1 dkξu k2 + kξp k2 + kξq k2 ≤ Chk+1 kξu k + C(hk+1 + kξu k)(kξpk + kξq k) 2 dt +Ckξr kkξp k + Ckξs kkξq k
which further implies
≤ C hk+1 + kξu k2 + kξs k2 + kξr k2 + kξp k2 + kξq k2 , 1 dkξu k2 ≤ C h2k+2 + kξu k2 + kξs k2 + kξr k2 . 2 dt
Moreover, plug (3.42) and (3.43) into (3.38) to yield
1 dkξv k2 + kξr k2 + kξs k2 ≤ C hk+1 + kξu k + kξv k kξv k + Chk+1(kξr k + kξs k) 2 dt 20
(3.44)
which further implies
1 1 ≤ C hk+1 + kξu k + kξv k kξv k + Ch2k+2 + kξr k2 + kξs k2 , 2 2
dkξv k2 + kξr k2 + kξs k2 ≤ C h2k+2 + kξu k2 + kξv k2 dt
(3.45)
Combining (3.44) and (3.45), we can find a constant γ such that dkξu k2 dkξv k2 2 2 +γ + kξr k + kξs k dt dt ≤ C h2k+2 + kξu k2 + kξs k2 + kξr k2 + γC h2k+2 + kξu k2 + kξv k2 Take γ to be sufficiently large, then
dkξu k2 dkξv k2 +γ ≤ Ch2k+2 + Ckξu k2 + Ckξv k2 . dt dt Therefore, by the Gronwall’s inequality and L2 initial discretization kξu k2 + kξv k2 ≤ Ch2k+2 , which further yield ku − uh k + kv − vh k ≤ Chk+1 by Lemma 3.2 and Lemma 3.7.
4
Positivity preserving property
In this subsection, we apply the positivity-preserving technique to construct positive numerical approximations. Following the idea in [48], we restrict ourselves to the P 1-LDG formulation. For simplicity, we consider Euler forward time discretization and assume the numerical solutions at time level n are positive: unh > 0, vhn > 0. Then we will construct positive numerical approximations at time level n + 1. To do so, we will firstly prove that the cell average at time level n + 1 is positive. Next, we can use a simple limiter to modify the DG polynomial and make it to be positive without changing its cell average. Throughout this section, we use oij for the numerical approximation oh in Kij , and the cell average is o¯ij . 21
Let us consider how to construct positive uh first, and the equation satisfied by the cell averages is u¯n+1 = Hijc (u, r, s) + Hijd (u, p, q), ij where Hijc (u, r, s)
∆t 1 = u¯nij − 2 ∆x∆y
Hijd (u, p, q)
1 ∆t = u¯nij + 2 ∆x∆y
Z Jj
cri− 1 ,j u cri+ 1 ,j − u 2
2
! Z dy + u csi,j+ 1 − u csi,j− 1 dx , 2
Ii
2
! Z Z qˆi,j+ 1 − qˆi,j− 1 dx . pˆi+ 1 ,j − pˆi− 1 ,j dy + Jj
2
2
Ii
2
2
For simplicity, if not otherwise stated, we will omit the subscript ij in Hijc and Hijd . To analyze H c , we approximate the integral by a 2-point Gauss quadrature. The Gauss quadrature points on Ii and Jj are denoted by pxi
=
n
xβi
: β = 1, 2
o
and
pyj
n o β = yj : β = 1, 2 ,
respectively. Also, we denote wβ as the corresponding weights on the interval − 21 , 21 . Let λ1 =
∆t ∆x
and λ2 =
∆t , ∆y
then H c becomes
X X 1 n wβ u cr i− 1 ,β − u wβ u csβ,j− 1 − u cr i+ 1 ,β + λ2 csβ,j+ 1 , H (u, r, s) = uij + λ1 2 2 2 2 2 2
2
β=1
β=1
c
where u cri− 1 ,β is the numerical flux at (xi− 1 , yjβ ). Likewise for the other fluxes. As the general 2
2
treatment, we rewrite the cell average on the right hand side as
1X 1X − − + + u + u = wβ u + w u , = β i+ 21 ,β β,j+ 21 i− 21 ,β β,j− 21 2 β=1 2 β=1 2
unij
2
where u− = u− (y β ) is a point value of the numerical approximation in the Gauss i− 1 ,β i− 1 ,j j 2
2
quadrature. Likewise for the other point values. Define µ = λ1 + λ2 , then 2 X 1 + − c ui− 1 ,β + ui+ 1 ,β + u cri+ 1 ,β H (u, r, s) = λ1 wβ cri− 1 ,β − u 2 2 2 2 4µ β=1 2 X 1 + − uβ,j− 1 + uβ,j+ 1 + u csβ,j− 1 − u csβ,j+ 1 . + λ2 wβ 2 2 2 2 4µ β=1 22
We need to suitably choose the parameter α in the Lax-Friedrichs flux, and the result is given below. Lemma 4.1. We can choose α>
max 1 ≤ i ≤ Nx − 1 1 ≤ j ≤ Ny − 1 β = 1, 2
o n − + − + , −r , s , −s , ri+ 1 ,β i+ 1 ,β β,j+ 1 β,j+ 1 2
2
2
2
then un+1 > 0 under a CFL condition j ∆t 1 ∆t A ≤ , + ∆x ∆y 2 where A=α−
min 1 ≤ i ≤ Nx − 1 1 ≤ j ≤ Ny − 1 β = 1, 2
(4.1)
o n + − + − ≥ 0, ri+ , −r , s 1 1 1 , −s 1 ,β i+ ,β β,j+ β,j+ 2
2
2
2
Proof. It is easy to check, for i = 2, 3, · · · , Nx , 1 + cri− 1 ,β u 1 +u 2 4µ i− 2 ,β 1 − 1 + − + + + − + u r − α(u − u ) = ui− 1 ,β ri− ui− 1 ,β + 1 1 1 1 1 ,β i− 2 ,β i− 2 ,β i− 2 ,β i− 2 ,β 2 2 2 4µ 2 1 1 1 − + = u− + − α u+ . α + ri− + ri− 1 1 1 ,β i− ,β ,β i− 21 ,β 2 2 2 2 2 2µ − Therefore, we can choose α+ri− > 0 and 1 ,β 2
+ 1 −α +ri− 1 2µ ,β 2
> 0 to obtain
1 + +c uri− 1 ,β u 4µ i− 21 ,β 2
0. If i = 1, then by (2.12), 1 + 1 + cri− 1 ,β = ui− 1 ,β + u u 1 > 0. 2 2 4µ 4µ i− 2 ,β
Also, for i = 1, 2, · · · , Ny − 1,
1 − cri+ 1 ,β u 1 −u 2 4µ i+ 2 ,β 1 − 1 − − + + + − + u r − α(u − u ) ui+ 1 ,β ri+ ui+ 1 ,β − = 1 ,β i+ 21 ,β i+ 12 ,β i+ 21 ,β i+ 21 ,β 2 2 2 4µ 2 1 1 1 + − + = α − ri+ − ri+ + − α u− . u 1 1 1 ,β ,β i+ 12 ,β i+ 2 ,β 2 2 2 2 2µ 23
>
+ > 0 and We can choose α − ri+ 1 ,β 2
1 2µ
− − α > 0 such that − ri+ 1 ,β 2
1 − u 4µ i+ 21 ,β
i = Ny , then by (2.12), 1 − 1 − ui+ 1 ,β − u u 1 > 0. cri+ 1 ,β = 2 2 4µ 4µ i+ 2 ,β
> 0 and Similarly, for j = 2, 3, · · · , Ny we require α + s− β,j− 1 2
1 + +c usβ,j− 1 u 4µ β,j− 21 2
such that
> 0. For j = 1, 2, · · · , Ny −1, we need
1 − u 4µ β,j+ 21
1 2µ
α−s+ β,j+ 21
−u cri+ 1 ,β > 0. If 2
+ s+ − α > 0 to obtain β,j− 1 2
> 0 and
1 −α −s− 2µ β,j+ 21
>0
−u csβ,j+ 1 > 0. 2
Now, we proceed to analyze H d (p, q). Following the same analysis in [48] with some minor changes, we can show the following lemma. Lemma 4.2. Suppose unij > 0 for all i and j, then H d (u, p, q) > 0 under a CFL condition ∆t 1 ∆t + ≤ . 2 2 ∆x ∆y 20 Proof. Define H d (u, p, q) =
Z
+
1 H 2 (u, q), ∆x
where
Z
Λ1 + − ui− 1 ,j + ui+ 1 ,j − λ1 pbi− 1 ,j − pbi+ 1 ,j dy = 2 2 2 2 Jj 4Λ Z Λ2 + − d ui,j− 1 + ui,j+ 1 − λ2 qbi,j− 1 − qbi,j+ 1 dx, H2 (u, q) = 2 2 2 2 Ii 4Λ
H1d (u, p)
with Λ1 =
1 H 1 (u, p) ∆y
∆t , ∆x2
Λ2 =
∆t ∆y 2
and Λ = Λ1 + Λ2 . It is easy to check that
Z
1 + − + 4ui+ 1 ,j − 3ui+ 1 ,j − ui− 1 ,j dy, i = 1, · · · , Nx − 1, = 2 2 2 Jj Jj ∆x Z Z 1 − + + − dy = − 3u − u dx, j = 1, · · · , Ny − 1. 4u qi,j+ 1 i,j+ 21 i,j− 21 i,j+ 21 2 Ii ∆y Ii p− dy i+ 12 ,j
Plug the flux (2.7) into (4.3) to obtain 1. i = 1 H1d (u, p) 24
(4.2)
(4.3) (4.4)
Λ1 + Λ1 − − ui+ 1 ,j + u 1 + λ1 pi+ 1 ,j dy = 2 2 4Λ i− 2 ,j Jj 4Λ Z Λ1 Λ1 + + − = − Λ1 ui− 1 ,j + − 3Λ1 ui+ 1 ,j + 4Λ1 ui+ 1 ,j dy 2 2 2 4Λ 4Λ Jj Z
> 0. 2. 2 ≤ i ≤ Nx − 1
H1d (u, p) Z Λ1 + Λ1 − − − = ui+ 1 ,j + u 1 + λ1 pi+ 1 ,j − pi− 1 ,j dy 2 2 2 4Λ i− 2 ,j Jj 4Λ Z Λ1 Λ1 − + + − + − 5Λ1 ui− 1 ,j + − 3Λ1 ui+ 1 ,j + 4Λ1 ui+ 1 ,j dy = Λ1 ui− 3 ,j + 3Λ1 ui− 1 ,j + 2 2 2 2 2 4Λ 4Λ Jj
> 0.
3. i = Nx H1d (u, p) Z Λ1 − Λ1 + − = ui+ 1 ,j + u 1 − λ1 pi− 1 ,j dy 2 2 4Λ i− 2 ,j Jj 4Λ Z Λ1 Λ1 − − + + − 4Λ1 ui− 1 ,j + u 1 dy = Λ1 ui− 3 ,j + 3Λ1 ui− 1 ,j + 2 2 2 4Λ 4Λ i+ 2 ,j Jj
> 0.
We can analyze H2d in a similar way, so we skip the proof. Based on the above two lemmas, we have the following theorem. Theorem 4.1. Suppose unij > 0 for all i and j, then un+1 >0 ij under the CFL conditions (4.1) and (4.2). Now, let us proceed to analyze v, and the equation satisfied by the numerical cell averages is v n+1 = ij
1 n v ij + H d (v, r, s) + ∆t unij − v nij 2 25
d
= H (v, r, s) +
1 − ∆t v nij + ∆tunij . 2
Applying Lemma 4.2, it is easy to prove the following theorem. Theorem 4.2. Suppose vijn > 0 for all i and j, then v n+1 >0 ij provided ∆t ∆t 1 + ≤ , 2 2 ∆x ∆y 20 and 1 ∆t ≤ . 2 Based on Theorems 4.1 and 4.2, the numerical cell averages we obtained are positive. However, the numerical solutions unij and vijn might still be negative. Hence, we have to modify the numerical solutions while keeping the cell averages untouched. For simplicity, we discuss the modification of unij only and the procedure is given in the following steps. • Set up a small number ε = 10−13 . • If unij > ε, then proceed to the following steps. Otherwise, unij is identified as the approximation to vacuum, and we will take u enij = unij as the numerical solution and
skip the following steps.
• Modify the density first: Compute bij = If bij < ε, then take u enij as with
min unij (x, y).
(x,y)∈Kij
u enij = unij + θij unij − unij , unij − ε , θij = n uij − bij
and use u enij as the new numerical density unij . 26
Following [42], we can show the L1 -stability of the numerical scheme with the positivitypreserving limiter. Since unh is positive, we have kunh kL1
=
Z
Ω
unh (x)dx
=
Z
Ω
u0h (x)dx = ku0h kL1 ,
where kukL1 is the standard L1 -norm of u on Ω. This implies the L1 -stability of the scheme.
4.1
High order time discretizations
All the previous analyses are based on first-order Euler forward time discretization. We can also use strong stability preserving (SSP) high-order time discretizations to solve the ODE system wt = Lw. More details of these time discretizations can be found in [36, 35, 20]. In this paper, we use the third order SSP Runge-Kutta method [36] w(1) = wn + ∆tL(wn ), 3 n w + 4 1 = wn + 3
w(2) = wn+1
1 w(1) + ∆tL(w(1) ) , 4 2 w(2) + ∆tL(w(2) ) , 3
(4.5)
and the third order SSP multi-step method [35] w
n+1
16 n 12 11 n−3 n n−3 = w + ∆tL(w ) . (w + 3∆tL(w )) + 27 27 11
(4.6)
Since an SSP time discretization is a convex combination of Euler forward, by using the limiter mentioned in Section 4, the numerical solutions obtained from the full scheme are also positive.
5
Numerical experiments
In this section, we present numerical examples in two space dimensions to verify the theoretical analysis and the positivity-preserving property of the proposed method. If not otherwise stated, we consider the third-order RKDG method (k = 2).
27
Example 5.1. We solve the following problem on the domain Ω = [0, 2π] × [0, 2π] ut − div(∇u − u∇v) = fu (x, y, t), vt − △v = u − v + fv (x, y, t),
x ∈ Ω, t > 0 x ∈ Ω, t > 0,
(5.1)
subject to boundary condition (1.2). We choose suitable fu and fv such that the exact solutions are u(x, y, t) = exp(−t)(cos(x) + cos(y)),
v(x, y, t) = exp(−t)(cos(x) + cos(y)).
We can easily check that the source terms are fu (x, y, t) = −2 exp(−2t)(cos2 (x) + cos(x) cos(y) + cos2 (y) − 1),
fv (x, y, t) = 0.
In Table 5.1, we present the numerical results for the proposed method without the bound preserving limiter. From the table, we can observe optimal rate of convergence with Vhk finite element spaces.
N 10 20 40 80
k=1 k=2 2 L error order L error order 2.69e-1 – 1.32e-2 – 7.16e-2 1.91 1.66e-3 2.99 1.74e-2 2.04 2.07e-4 3.00 4.25e-3 2.03 2.59e-5 3.00 2
Table 5.1: Example 5.1: accuracy test at T = 0.1 for the second- and third-order LDG methods without the positivity-preserving limiter. k is the degree of polynomial and Nx = Ny = N.
Example 5.2. In this example, we solve (1.1) with the following initial condition. u0 = 840 exp(−84(x2 + y 2)),
v0 = 420 exp(−42(x2 + y 2 )).
We use second-order positivity-preserving LDG methods to solve and the numerical approximations at time t is given in figure 5.1. In [15], the authors demonstrated that the blow-up time should be approximately t = 1.21 × 10−4 . However, we can continue our numerical simulation to t = 2 × 10−4 , and the numerical approximation is given in figure 5.2. We also solve 28
Figure 5.1: Example 5.2: Numerical approximations of u at t = 6 × 10−5 (left) and t = 1.2 × 10−4 with positivity-preserving limiter for P 1 polynomials and N = 160. .
Figure 5.2: Example 5.2: Numerical approximations of u at t = 2.0 × 10−4 with positivitypreserving limiter for P 1 polynomials and N = 160. .
29
the problem without the positivity-preserving limiter. We find that at about t = 8 × 10−5 the numerical scheme yields negative u. Moreover, we also test the numerical blow-up time. For simplicity, we take Nx = Ny = N, and compute the L2 -norm numerical approximations at time t with N × N cells, defined as S(N, t). Define tb (N) = inf{t : S(2N, t) ∗ 1.05% ≥ S(N, t)}
(5.2)
as the numerical blow-up time. We anticipate that as we refine the mesh, the numerical blow-up time will converge to the exact value. However, due to the computational cost, we have to apply adaptive methods and refine the mesh in the vicinity of the blow-up point, and this work will be considered in the future. To verify our anticipation, we consider the
Figure 5.3: L2 -norm of the numerical approximation for Example 5.2 with different N’s. . following function 1 exp f (x, t) = 1−t
−x2 2(t − 1)2
It is easy to see that lim f (x, t) → δ(x)
t→1−
30
.
(5.3)
in the sense of distribution. We consider the interval [-1,1] and divided in into N uniform cells. Denote uh to be the L2 -projection of f (x, t) in each cell and S(N, t) to be the L2 -norm of uh over [-1,1] with different t′ s. We also compute the numerical blow-up time tb (N) by (5.2) and the results are given in Table 5.2. We can clearly see that as we refine the meshes, tb (N) converges to 1, the exact blow-up time, in first order accuracy. N Blow-up time Error
10 -
20 40 80 160 320 0.671 0.836 0.918 0.959 0.980 0.329 0.164 0.092 0.041 0.020
Table 5.2: Numerical blow-up time with different mesh sizes for (5.3).
6
Conclusion
In this paper, we develop LDG methods to the KS chemotaxis model. We improve the result given by [15], and give the optimal error estimate under special finite element spaces. A special positivity-preserving limiter is constructed to obtain physically relevant numerical approximations. Moreover, we also prove the L1 -stability of the LDG scheme with the limiter. Numerical experiments are given to demonstrate the good performance of the LDG scheme and the estimate of the numerical blow-up time. In the future, we plan to apply adaptive methods to calculate the blow-up time more accurately.
7
Acknowledgments
Xingjie Helen Li would like to thank the Shanghai Centre for Mathematics Science (SCMS), Fudan University, for support during her visit.
References [1] F. Bassi and S. Rebay, A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations, Journal of Compu31
tational Physics, 131 (1997), 267-279. [2] A. Chertock and A. Kurganov, A second-order positivity preserving central-upwind scheme for chemotaxis and haptotaxis models, Numerische Mathematik, 111 (2008), 169-205. [3] S. Childress and J. Percus, Nonlinear aspects of chemotaxis, Mathematical Biosciences, 56 (1981), 217-237. [4] P. Ciarlet, Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [5] B. Cockburn and B. Dong, An analysis of the minimal dissipation local discontinuous Galerkin method for convection-diffusion problems, Journal of Scientific Computing, 32 (2007), 233-262. [6] B. Cockburn, S. Hou and C.-W. Shu, The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case, Mathematics of Computation, 54 (1990), 545-581. [7] B. Cockburn, G. Kanschat, I. Perugia, and D. Schotzau, Superconvergence of the local discontinuous galerkin method for elliptic problems on cartesian grids, SIAM Journal on Numerical Analysis, 39 (2001), 264285. [8] B. Cockburn, S.-Y. Lin and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one-dimensional systems, Journal of Computational Physics, 84 (1989), 90-113. [9] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework, Mathematics of Computation, 52 (1989), 411-435.
32
[10] B. Cockburn and C.-W. Shu, The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems, Journal of Computational Physics, 141 (1998), 199-224. [11] B. Cockburn and C.-W. Shu, The local discontinuous Galerkin method for time dependent convection-diffusion systems, SIAM Journal on Numerical Analysis, 35 (1998), 2440-2463. [12] Y. Epshteyn, Discontinuous Galerkin methods for the chemotaxis and haptotaxis models, Journal of Computational and Applied Mathematics, 224 (2009), 168-181. [13] Y. Epshteyn, Upwind-difference potentials method for Patlak-Keller-Segel chemotaxis model, Journal of Scientific Computing, 53 (2012), 689-713. [14] Y. Epshteyn and A. Izmirlioglu, Fully discrete analysis of a discontinuous finite element method for the Keller-Segel Chemotaxis model, Journal of Scientific Computing, 40 (2009), 211-256. [15] Y. Epshteyn and A. Kurganov, New interior penalty discontinuous Galerkin methods for the Keller-Segel Chemotaxis model, SIAM Journal on Numerical Analysis, 47 (2008), 368-408. [16] F. Filbet, A finite volume scheme for the Patlak-Keller-Segel chemotaxis model, Numerische Mathematik, 104 (2006), 457-488. [17] F. Filbet, P. Lauren¸cot and B. Perthame, Derivation of hyperbolic models for chemosensitive movement, Journal of Mathematical Biology, 50 (2005), 189-207. [18] H. Gajewski and K. Zacharias, Global behaviour of a reaction-diffusion system modelling chemotaxis, Mathematische Nachrichten, 195 (1998), 77-114. [19] I.M. Gelfand, Some questions of analysis and differential equations, American Mathematical Society Translations, 26 (1963), 201-219. 33
[20] S. Gottlieb, C.-W. Shu and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Review, 43 (2001), 89-112. [21] L. Guo and Y. Yang, Positivity-preserving high-order local discontinuous Galerkin method for parabolic equations with blow-up solutions, Journal of Computational Physics, 289 (2015), 181-195. [22] J. Hakovec and C. Schmeiser, Stochastic particle approximation for measure valued solutions of the 2d Keller-Segel system, Journal of Statistical Physics, 135 (2009), 133-151. [23] M.A. Herrero, E. Medina and J.J.L. Velazquez, Finite-time aggregation into a single point in a reaction-diffusion system, Nonlinearity, 10 (1997), 1739-1754. [24] D. Horstman, From 1970 until now: The Keller-Segel model in chemotaxis and its consequences I, Jahresbericht DMV, 105 (2003), 103-165. [25] D. Horstmann, From 1970 until now: The Keller-Segel model in chemotaxis and its consequences II, Jahresbericht DMV, 106 (2004), 51-69. [26] A.E. Hurd and D.H. Sattinger, Questions of existence and uniqueness for hyperbolic equations with discontinuous coefficients, Transactions of the American Mathematical Society, 132 (1968), 159-174. [27] E.F. Keller and L.A. Segel, Initiation on slime mold aggregation viewed as instability, Journal of Theoretical Biology, 26 (1970), 399-415. [28] E.F. Keller and L.A. Segel, Traveling bands of chemotactic bacteria, a theoretical analysis, Journal of Theoretical Biology, 30 (1971), 235-248. [29] A. Marrocco, 2D simulation of chemotaxis bacteria aggregation, ESAIM: Mathematical Modelling and Numerical Analysis (M 2 AN), 37 (2003), 617-630.
34
[30] X. Meng, C.-W. Shu, Q. Zhang and B. Wu, Superconvergence of discontinuous Galerkin methods for scalar nonlinear conservation laws in one space dimension, SIAM Journal on Numerical Analysis, 50 (2012), 2336-2356. [31] C. Patlak, Random walk with persistence and external bias, The Bulletin of Mathematical Biophysics, 15 (1953), 311-338. [32] T. Qin, C.-W. Shu and Y. Yang, Bound-preserving discontinuous Galerkin methods for relativistic hydrodynamics, submitted to Journal of Computational Physics. [33] W.H. Reed and T.R. Hill, Triangular mesh methods for the Neutron transport equation, Los Alamos Scientific Laboratory Report LA-UR-73-479, Los Alamos, NM, 1973. [34] N. Saito, Conservative upwind finite-element method for a simplified Keller-Segel system modelling chemotaxis, IMA Journal of Numerical Analysis, 27 (2007), 332-365. [35] C.-W. Shu, Total-variation-diminishing time discretizations, SIAM Journal on Scientific and Statistical Computing, 9 (1988), 1073-1084. [36] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shockcapturing schemes, Journal of Computational Physics, 77 (1988), 439-471. [37] R. Strehl, A. Sokolov, D. Kuzmin, D. Horstmann, and S. Turek, A positivity-preserving finite element method for chemotaxis problems in 3D, Journal of Computational and Applied Mathematics, 239 (2013), 290-303. [38] R. Tyson, L.J. Stern and R.J. LeVeque, Fractional step methods applied to a chemotaxis model, Journal of Mathematical Biology, 41 (2000), 455-475. [39] H. Wang, C.-W. Shu and Q. Zhang, Stability and error estimates of local discontinuous Galerkin methods with implicit-explicit time-marching for advection-diffusion problems, SIAM Journal on Numerical Analysis, 53 (2015), 206-227.
35
[40] H. Wang, S. Wang, C.-W. Shu and Q. Zhang, Local discontinuous Galerkin methods with implicit-explicit time-marching for multi-dimensional convection-diffusion problems, ESAIM: Mathematical Modelling and Numerical Analysis (M 2 AN), to appear. [41] Y. Yang and C.-W. Shu, Discontinuous Galerkin method for hyperbolic equations involving δ-singularities: negative-order norm error estimates and applications, Numerische Mathematik, 124, (2013), 753-781. [42] Y. Yang, D. Wei and C.-W. Shu, Discontinuous Galerkin method for Krause’s consensus models and pressureless Euler equations, Journal of Computational Physics, 252 (2013), 109-127. [43] Q. Zhang and C.-W. Shu, Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws, SIAM Journal on Numerical Analysis, 42 (2004), 641-666. [44] Q. Zhang and C.-W. Shu, Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws, SIAM Journal on Numerical Analysis, 48 (2010), 1038-1063. [45] X. Zhang and C.-W. Shu, On maximum-principle-satisfying high order schemes for scalar conservation laws, Journal of Computational Physics, 229 (2010), 3091-3120. [46] X. Zhang and C.-W. Shu, On positivity preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes, Journal of Computational Physics, 229 (2010), 8918-8934. [47] X. Zhang and C.-W. Shu, Positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations with source terms, Journal of Computational Physics, 230 (2011), 1238-1248.
36
[48] Y. Zhang, X. Zhang and C.-W. Shu, Maximum-principle-satisfying second order discontinuous Galerkin schemes for convection-diffusion equations on triangular meshes, Journal of Computational Physics, 234 (2013), 295-316. [49] X. Zhao, Y. Yang and C. Syler, A positivity-preserving semi-implicit discontinuous Galerkin scheme for solving extended magnetohydrodynamics equations, Journal of Computational Physics, 278 (2014), 400-415.
37
Keywords: Local discontinuous Galerkin method, Keller-Segel chemotaxis model, positivity preserving, error estimate, Neumann boundary condition, blow-up, L1 stability 1
Department of Mathematics and Statistics, University of North Carolina at Charlotte, Charlotte, NC 28223. E-mail: [email protected] 2 Division of Applied Mathematics, Brown University, Providence, RI 02912. E-mail: [email protected]. Research supported by ARO grant W911NF-15-1-0226 and NSF grant DMS1418750. 3 Department of Mathematical Sciences, Michigan Technological University, Houghton, MI 49931. E-mail: [email protected]
1
1
Introduction
In this paper, we study the Keller-Segel (KS) chemotaxis model in two space dimensions [31, 27] and focus on the following common formulation [3], ut − div(∇u − χu∇v) = 0, vt − △v = u − v,
x ∈ Ω, t > 0, x ∈ Ω, t > 0,
(1.1)
where Ω is assumed to be a convex, bounded and open set in R2 . Chemotaxis is the highly nonlinear term which indicates movements by cells in reaction to a chemical substance, where cells approach chemically favorable environments and avoid unpleasant ones. In (1.1), u and v denote the densities of cells and the chemical concentration, respectively. The chemotactic sensitivity function χ is assumed to be a positive constant. For simplicity, we take χ ≡ 1. The boundary conditions are set to be ∇u · n = ∇v · n = 0, where n is the outer normal of the boundary ∂Ω. With this boundary condition, constant during the time evolution and the system is thus isolated.
(1.2) R
Ω
u is a
The exact solutions of the KS chemotaxis model are always positive. Moreover, the model exhibits blow-up patterns with certain initial conditions [18, 17]. When the blow-up patterns occur, the density u of cells will strengthen in the neighborhood of isolated points, and these regions become sharper and eventually result in finite time point-wise blow-up. Mathematical proofs for spherically symmetric solutions in a ball have been investigated in [23]. The blow-up point in this case is the center of the ball, and it is proved to be the only possible singularity. For nonsymmetric cases, blow-up results have also been studied in [24, 25], and the blow-up points will be located on the boundary. More theoretical works can be found in [28, 18]. It is difficult to construct numerical schemes for (1.1), and most of the previous works are for the following simplified system ut − div(∇u − χu∇v) = 0, −△v = u − v, 2
x ∈ Ω, t > 0, x ∈ Ω, t > 0,
(See, for example, [38, 29, 16, 22] and the references therein). Recently, there are some significant works designed to solve (1.1) directly [15, 37]. In [37], the authors constructed implicit second-order positivity preserving finite-volume schemes in three-dimensional space, and their technique requires solving a large linear system of equations coupling together all grid points at each stage of the two stage TR-BDF2 method when updating the diffusion terms at each time step. In [15], the authors applied the interior penalty discontinuous Galerkin (IPDG) method on rectangular meshes to obtain suboptimal rate of convergence, and the finite element space is assumed to be piecewise polynomials of degree k ≥ 2. Other related works in this direction include [14, 12, 13] and the positivity-preserving property was demonstrated by numerical experiments only. Besides the above, in [2], the authors constructed a second order positivity-preserving scheme to a revised system by differentiating (1.1) with respect to x and y. Hence, the scheme was not designed to solve (1.1) directly, similar work can also be found in [34]. Recently, Zhang and Shu introduced positivity-preserving limiter to hyperbolic equations in [45, 46, 47]. Subsequently, the idea was extended to parabolic problems in [48]. In this paper, we will apply the positivity-preserving limiter introduced in [48] to construct second-order positivity-preserving local discontinuous Galerkin (LDG) schemes to obtain physically relevant numerical approximations. The method we plan to use preserves the positivity of the numerical solutions, and can be applied to unstructured meshes. The DG method was first introduced in 1973 by Reed and Hill [33] in the framework of neutron linear transport. Subsequently, Cockburn et al. developed Runge-Kutta discontinuous Galerkin (RKDG) methods for hyperbolic conservation laws in a series of papers [9, 6, 8, 10]. In [11], Cockburn and Shu introduced the LDG method to solve the convectiondiffusion equations. Their idea was motivated by Bassi and Rebay [1], where the compressible Navier-Stokes equations were successfully solved. Recently, the DG methods were applied to linear hyperbolic equations with δ-singularities [41] to obtain high-order approximations under suitable negative-order norms. Subsequently, the methods have also been applied to
3
nonlinear hyperbolic equations with δ-singularities [42, 49]. Combined with special limiters, the schemes were proved to be L1 stable [42, 32]. Recently, the idea has been extended to parabolic equations with blow-up solutions by using the LDG method [21]. In this paper, we follow the same direction and employ the LDG method to capture the blow-up phenomenon. In the LDG method, we introduce auxiliary variables p = ux and q = uy . Numerical experiments will be given to demonstrate the optimal rate of convergence. However, in this problem the approximations of p and q are discontinuous across the cell interfaces and we can hardly obtain error estimates if we analyze the convection and diffusion terms separately. To explain this point, let us consider the following hyperbolic equation ut + (a(x)u)x = 0, where a(x) is discontinuous at x = x0 . In [19, 26], the authors studied such a problem and defined Q=
a(x0 + b) − a(x0 ) . b
If Q is bounded from below for all b, then the solution exists, but may not be unique. If Q is bounded from above for all b, we can guarantee the uniqueness, but the solution may not exist. To bridge this gap, most of the previous works for similar problems are based on mixed finite element methods (see e.g. [29]). In this paper, we consider a new technique for error analysis following the idea in [39, 40]. Moreover, we also apply the positivity-preserving limiter to guarantee positivity of the numerical approximation. With this special technique, the L1 norm of the numerical approximation is a constant during the time evolution due to the mass conservation. However, the L2 norm might still be unbounded when blow-up patterns occur. We thus introduce a new idea to capture the numerical blow-up time based on the L2 norm of the numerical approximations. To our best knowledge, this is the first paper studying the numerical blow-up time with a concrete theoretical background. The organization of this paper is as follows. In Section 2, we construct the LDG scheme. In Section 3, we give the error estimates based on two different finite element spaces. In 4
Section 4, we discuss the positivity-preserving technique, the foundation of limiters and high order time discretizations. Numerical experiments are given in Section 5. Finally, we will end in Section 6 with concluding remarks and remarks for future work.
2
LDG scheme
In this section, we define the finite element spaces and proceed to construct the LDG scheme. In this paper, we consider the computational domain to be Ω = [0, 1] × [0, 1] and use rectangular meshes for simplicity. The analyses for triangular meshes can be extended with some minor changes, so we will skip this part. Let Kij = (xi− 1 , xi+ 1 ) × (yj− 1 , yj+ 1 ), i = 2
2
2
2
1, . . . , Nx , j = 1, . . . , Ny , be a partition of the computational domain Ω, and the mesh sizes in the x and y directions are given as ∆xi = xi+ 1 −xi− 1 and ∆yj = yj+ 1 −yj− 1 , respectively. For 2
2
2
2
simplicity, in this paper, we assume uniform meshes and denote h = ∆xi = ∆yj . However, this assumption is not essential in the proof. Moreover, we also denote Ii = [xi− 1 , xi+ 1 ] and 2
2
Jj = [yj− 1 , yj+ 1 ] to be the cells in x and y directions. 2
2
We define the finite element spaces Vhk as n o Vhk = z : z Kij ∈ P k (Kij ) ,
where P k (Kij ) denotes the set of polynomials of degree up to k in cell Kij . To construct the LDG scheme, we introduce the axillary variables p = ux , q = uy , r = vx and s = vy , then (1.1) can be written as ut = −(ru)x − (su)y + px + qy , p = ux , q = uy , vt = rx + sy + u − v, r = vx , s = vy . 5
For simplicity, we use notations o to represent the exact solutions of the six variables u, p, q, v, r, s and oh for the corresponding numerical approximations. Then the LDG scheme is to find uh , ph , qh ∈ Vhk1 and vh , rh , sh ∈ Vhk2 , such that for any test functions w u , w p, w q ∈ Vhk1 and w v , w r , w s ∈ Vhk2 u
(uh t , w )ij =
(rh uh , wxu )ij
+
Z
Jj
Z
u + rd h uh (w )i− 1 ,j dy 2
−
Z
Z
Jj
u − rd h uh (w )i+ 1 ,j dy 2
u − u + sd + (sh uh , wyu )ij + sd h uh (w )i,j+ 1 dx h uh (w )i,j− 1 dx − 2 2 Ii Z Ii Z − (ph , wxu )ij − pbh (w u )+ dy + pbh (w u )− dy i− 12 ,j i+ 21 ,j Jj Jj Z Z u u + qbh (w u )− dx, − (qh , wy )ij − qbh (w )i,j− 1 dx + i,j+ 21 2 Ii Ii Z Z p p p + uh (w p )− c dy, (ph , w )ij = −(uh , wx )ij − uh (w )i− 1 ,j dy + c i+ 21 ,j 2 Jj Jj Z Z q q q + uh (w q )− c dx, (qh , w )ij = −(uh , wy )ij − uh (w )i,j− 1 dx + c i,j+ 21 2 Ii Ii Z Z v v v + rbh (w v )− dy rbh (w )i− 1 ,j dy + (vht , w )ij = −(rh , wx )ij − i+ 21 ,j 2 Jj Jj Z Z v v + − (sh , wy )ij − sbh (w )i,j− 1 dx + sbh (w v )− dx i,j+ 1 2
Ii
2
where (u, v)ij = for
− wi+ 1 , ,j 2
+ wi,j− 1 2
R
Kij
and
(2.2) (2.3)
2
Ii
+ (uh − vh , w v )ij , Z Z r + r r vbh (w r )− vbh (w )i− 1 ,j dy + dy, (rh , w )ij = −(vh , wx )ij − i+ 12 ,j 2 Jj Jj Z Z s s s + (sh , w )ij = −(vh , wy )ij − vbh (w )i,j− 1 dx + vbh (w s )− dx, i,j+ 1 Ii
(2.1)
(2.4) (2.5) (2.6)
2
Ii
+ is the trace of w at x = xi− 1 in cell Kij . Likewise uvdxdy, and wi− 1 ,j 2
2
− wi,j+ 1. 2
We also denote rd h uh , sd h uh and obh to be the numerical fluxes.
In this paper, we use the alternating fluxes for the diffusion terms and Lax-Friedrichs fluxes for the convection terms. For simplicity, if we consider the generic trace at the interior cell interfaces, the subscript will be omitted. Due to the Neumann boundary condition (1.2), we take
6
1. For j = 1, 2, · · · , Ny , + , i = 0, u 1 ,0 h i+ 2 ,j − uh + , i = 1, · · · , Nx − 1, u ch i+ 1 ,j , pbh i+ 1 ,j = 1 , ph 1 i+ ,j 2 2 2 i+ 2 ,j uh − 1 , 0 , i = Nx , i+ ,j
(2.7)
2
+ , 0 , i = 0, v h i+ 1 ,j 2 − vh + , i = 1, · · · , Nx − 1, vbh i+ 1 ,j , rbh i+ 1 ,j = 1 , rh 1 i+ ,j 2 2 2 i+ 2 ,j vh − 1 , 0 , i = Nx . i+ ,j
(2.8)
2
2. For i = 1, 2, · · · , Nx
+ , 0 , j = 0, u h 1 i+ 2 ,j − uh + , j = 1, · · · , Ny − 1, uh i,j+ 1 , qbh i,j+ 1 = c 1 , qh 1 i,j+ 2 2 2 i,j+ 2 uh − 1 , 0 , j = Ny , i,j+
(2.9)
2
+ v , 0 , j = 0, h i+ 1 ,j 2 , sh − , j = 1, · · · , Ny − 1, vh + vbh i,j+ 1 , sbh i,j+ 1 = 1 i,j+ 21 2 2 i,j+ 2 vh − 1 , 0 , j = Ny . i,j+
(2.10)
2
For the convection term, we consider Lax-Friedrichs flux rd h uh =
1 + + + − rh uh + rh− u− h − α(uh − uh ) , 2
(2.11)
where α is chosen by the positivity-preserving technique. Due to the mass conservation, we also take rd h uh Nx + 1 ,j = 0, h uh 1 ,j = rd 2
2
j = 1, · · · , Ny .
(2.12)
We also define sd h uh is a similar way.
Also let [w] = w + − w − denote the jump of w at the interface of each cell and (u, v) =
Ny Nx X X
(u, v)ij ,
i=1 j=1
then (2.1)-(2.6) can be written as (uht , w u ) = Lcx (rh , uh , w u ) + Lcy (sh , uh , w u ) − Ldx (ph , w u ) − Ldy (qh , w u ), 7
(2.13)
(ph , w p ) = −Qx (uh , w p ),
(2.14)
(qh , w q ) = −Qy (uh , w q ),
(2.15)
(vht , w v ) = −Ldx (rh , w v ) − Ldy (sh , w v ) + (uh − vh , w v ),
(2.16)
(rh , w r ) = −Qx (vh , w r ),
(2.17)
(sh , w s ) = −Qy (vh , w s ),
(2.18)
where Lcx (r, u, w)
Ny Z y Nx −1 X 1 j+ 2 1 X r + u+ + r − u− − α[u] [w]i+ 1 ,j dy, = (ru, wx ) + 2 2 i=1 j=1 y 1
(2.19)
j− 2
Lcy (s, u, w) = (su, wy ) + Ldx (p, w) = (p, wx ) +
1 2
y −1 Z x Nx N X X i+ 1
2
i=1 j=1
2
Q (u, w) = (u, wx ) +
yj− 1
2
xi− 1
q − [w]i,j+ 1 dx,
(2.22)
2
2
u [w]i+ 1 ,j dy + 2
Ny Z y N x −1 X X j+ 1 2
yj− 1
xi− 1
y −1 Z Nx N X X
i=1 j=1
yj+ 1 2
u− [w]Nx + 1 ,j dy (2.23) 2
yj− 1
2
− [u]wi+ 1 dy, ,j
(2.24)
2
2
2
i=1 j=0
Ny Z X j=1
2
y −1 Z x Nx N X X i+ 1
= −(uy , w) −
+
yj− 1
i=1 j=1
Q (u, w) = (u, wy ) +
(2.21)
2
2
i=0 j=1
y
p− [w]i+ 1 ,j dy,
Ny Z y N x −1 X X j+ 1
= −(ux , w) −
(2.20)
2
2
y −1 Z x Nx N X X i+ 1
i=1 j=1
x
2
Ny Z y 1 N x −1 X X j+ i=1 j=1
Ldy (q, w) = (q, wy ) +
xi− 1
s+ u+ + s− u− − α[u] [w]i,j+ 1 dx,
+
u [w]i,j+ 1 dx + 2
Nx Z X i=1
2
xi+ 1
2
xi+ 1
2
xi− 1
u− [w]i,Ny + 1 dx (2.25) 2
2
− [u]wi,j+ 1 dx.
(2.26)
2
xi− 1 2
+ + Here, for simplicity, we take w −1 ,j = wN = wi,−1 = wi,N 1 1 = 0. It is easy to check the x + ,j y+ 2
2
2
2
following identities by integration by parts on each cell: for any functions u and w, Ldx (u, w) + Qx (w, u) = 0,
(2.27)
Ldy (u, w) + Qy (w, u) = 0.
(2.28)
8
3
Error estimates
In this section, we analyze the error between the numerical and exact solutions. We first introduce the norms, construct the error equations and state the main theorem. Then we proceed to the proof. The whole procedure can be divided into several steps. We first construct the special projections that will be used in this section. Then we give an a priori error estimate and its verification. Finally, we will consider another finite element space and the corresponding error estimate. Let us define the norms first.
3.1
Norms
In this subsection ,we define several norms that will be used throughout the paper. Denote kuk0,Kij to be the standard L2 norm of u in cell Kij . For any natural number ℓ, we consider the norm of the Sobolev space H ℓ (Kij ), defined by ) 12 (
X
∂ α+β u 2
. kukℓ,Kij =
∂ α x∂ β y 0,Kij 0≤α+β≤ℓ
Moreover, we define the norms on the whole computational domain as ! 12 Ny Nx X X kukℓ = . kuk2ℓ,Kij i=1 j=1
For convenience, if we consider the standard L2 norm, then the corresponding subscript will be omitted. Let Γij be the edges of Kij , and we define Z Z 2 − 2 2 − 2 2 + (u+ dx. kukΓij = (ui− 1 ,j ) + (ui+ 1 ,j ) dy + 1 ) + (u 1) i,j− i,j+ Jj
2
2
Ii
2
2
Moreover, we also define
kuk2Γh
=
Ny Nx X X
kuk2Γij .
i=1 j=1
Finally, we define the standard L∞ norm of u in Kij as kuk∞,Kij , and define the L∞ norm on the whole computational domain as kuk∞ = max kuk∞,Kij . i,j
9
3.2
Error equations and the main theorem
In this subsection, we proceed to construct the error equations. Denote the error between the exact and numerical solutions to be eo = o − oh , then we have the equations satisfied by the errors as (eut , w u ) =Lcx (r, u, w u) − Lcx (rh , uh , w u ) + Lcy (s, u, w u) − Lyc (sh , uh , w u ) − Ldx (ep , w u ) − Ldy (eq , w u ),
(3.1)
(ep , w p ) = − Qx (eu , w p ),
(3.2)
(eq , w q ) = − Qy (eu , w q ),
(3.3)
(ev t , w v ) = − Ldx (er , w v ) − Ldy (es , w v ) + (eu − ev , w v ),
(3.4)
(er , w r ) = − Qx (ev , w r ),
(3.5)
(es , w s ) = − Qy (ev , w s).
(3.6)
Now, we can state our main theorem. Theorem 3.1. Suppose u and v are smooth functions, r and s are uniformly bounded for t ≤ T . The numerical approximations uh , ph , qh ∈ Vhk1 and vh , sh , rh ∈ Vhk2 . The initial discretization is given as the standard L2 -projection. If we take k1 ≥ 1 and k2 ≥ 2, then we have ku − uh k + kv − vh k ≤ Chmin{k1 +1,k2 } , where the positive constant C does not depend on h.
3.3
Preliminaries and projections
In this subsection, we study the basic properties of the finite element space. Let us start with the classical inverse properties. Lemma 3.1. Assuming u ∈ Vhk , there exists a positive constant C independent of h and u such that hkuk∞,Kij + h1/2 kukΓij ≤ CkukKij . 10
In this paper, we consider several special projections. We define the elliptic projections Pu, Pp, Pq ∈ Vhk1 such that for any w u , w p, w q ∈ Vhk1 Ldx (p, w u ) + Ldy (q, w u) =Ldx (Pp, w u ) + Ldy (Pq, w u ),
(3.7)
(Pp, w p ) + (Pq, w q ) = − Qx (Pu, w p ) − Qy (Pu, w q ).
(3.8)
Similarly, we also define Pv, Pr, Ps ∈ Vhk2 such that for any w v , w r , w s ∈ Vhk2 Ldx (r, w v ) + Ldy (s, w v ) =Ldx (Pr, w v ) + Ldy (Ps, w v ),
(3.9)
(Pr, w r ) + (Ps, w s ) = − Qx (Pv, w r ) − Qy (Pv, w s ).
(3.10)
The existence of the elliptic projections will be given in the following lemma [5], Lemma 3.2. Suppose
R
Pudxdy = Ω
R
udxdy and Ω
R
Pvdxdy = Ω
R
Ω
vdxdy, then Po are
uniquely determined by (3.7)-(3.10). Moreover, we have the following estimates kp − Ppk + kq − Pqk ≤Chk1 , ku − Puk ≤Chk1 +1 , kr − Prk + ks − Psk ≤Chk2 , kv − Pvk ≤Chk2 +1 ,
(3.11) (3.12) (3.13) (3.14)
where C is independent of h. In addition, we also define the standard L2 projection P by (P u, v)ij = (u, v)ij ,
∀v ∈ P k (Kij ).
The L2 projection satisfies the following property [4]. Lemma 3.3. Suppose u ∈ C k+1 (Ω). Then there exists a positive constant C independent of h and u such that ku − P uk + hku − P uk∞ + h1/2 ku − P ukΓh ≤ Chk+1 kukk+1. 11
Let us finish this section by proving the following lemma. Lemma 3.4. Let u ∈ C k+1(Ω) and Πu ∈ Vhk . Suppose ku − Πuk ≤ Chκ for some positive constant C and κ ≤ k + 1. Then hku − Πuk∞ + h1/2 ku − ΠukΓh ≤ Chκ , where the positive constant C does not depend on h. Proof. Actually, hku − Πuk∞ + h1/2 ku − ΠukΓh ≤ hkP u − Πuk∞ + hku − P uk∞ + h1/2 kP u − ΠukΓh + h1/2 ku − P ukΓh ≤ CkP u − Πuk + Chk+1 + CkP u − Πuk + Chk+1 ≤ Cku − Πuk + Cku − P uk + Chk+1 ≤ Chκ , where the first and third steps follow from triangle inequality, the second step requires Lemma 3.1 and Lemma 3.3, and the last step we use Lemma 3.3.
3.4
A priori error estimate
In this subsection, we would like to make an a priori error estimate that ku − uh k ≤ h,
(3.15)
which further implies ku − uh k∞ ≤ C by Lemma 3.4. Moreover, if u is bounded, then kuh k∞ ≤ C. This assumption, which will be verified in Section 3.6, is used for the estimate of the convection terms. This idea has been used to obtain the error estimates for nonlinear hyperbolic equations [43, 44, 30]. With this assumption, we can proceed to the main proof of the theorem.
12
3.5
Proof of Theorem 3.1
In this subsection, we give the proof of Theorem 3.1. As the general treatment of the finite element methods, we divide the error as eo = ηo − ξo , where ηo = o − Po and ξo = oh − Po . Clearly, ξo is chosen from the desired finite element space, and following [39, 40] we have Lemma 3.5. k(ξu )x k2 + k(ξu )y k2 + h−1 k[ξu ]k2Γh ≤ C(kξp k2 + kξq k2 ). Proof. By (2.14), (2.15) and (3.8), we have (ξp , w p ) + (ξq , w q ) = −Qx (ξu , w p ) − Qy (ξu , w q ) We take w p (x, y) =
xi+ 1 −x 2
∆xi
(3.16)
(ξu )x , (x, y) ∈ Kij and w q = 0, then by (2.24) and (2.14) it is
easy to show |||(ξu )x |||2Kij where
xi+ 1 − x 2 (ξu )x ≤ k(ξu )x kKij kξp kKij . = ξp , ∆xi ij |||u|||2Kij
=
Z
xi+ 1 − x 2
∆xi
Kij
u2 dxdy.
2
Clearly, ||| · |||Kij is a weighted L -norm in P k1 (Kij ), and is equivalent to the standard L2 -norm, i.e. for any u ∈ P k1 (Kij ), kukKij ≤ C|||u|||Kij , where the constant C does not depend on Kij . Therefore, k(ξu )x k2Kij ≤ C|||(ξu )x |||2Kij ≤ Ck(ξu )x kKij kξp kKij , which further yields k(ξu )x kKij ≤ Ckξp kKij . Similarly, we can also prove k(ξu )y k ≤ Ckξq k 13
following the same line. Next, we will take care of interfacial terms. We take w p (x, y) = ξu (x, y) − ξu+ (xi+ 1 , y), 2
(x, y) ∈ Kij and w q = 0 in (3.16) to obtain in cell Kij , Z
Jj
[ξu (xi+ 1 , y)]2 dy = (ξp − (ξu )x , w p )ij ≤ Ckw p kKij kξp kKij . 2
It is easy to see that p
w (x, y) = −[ξu (xi+ 1 , y)] − 2
Z
x+ 12
(ξu )x dx. x
Therefore, by Minkowski’s inequality, Z
[ξu (xi+ 1 , y)] dy ≤ Ckξp kKij [ξu (xi+ 1 , y)] + ∆xi k(ξu )x kKij 2 2 Kij Jj !1/2 Z 1/2 ≤ Ckξp kKij ∆xi [ξu (xi+ 1 , y)]2dy + ∆xi kξp kKij 2
1 ≤ 2
which further yields
Similarly, we have
Jj
Z
Z
Jj
Ii
[ξu (xi+ 1 , y)]2dy + C∆xi kξp k2Kij 2
[ξu (xi+ 1 , y)]2 dy ≤ C∆xi kξp k2Kij .
(3.17)
[ξu (x, yj+ 1 )]2 dx ≤ C∆yj kξq k2Kij .
(3.18)
Jj
Z
2
2
2
Hence, h−1 k[ξu ]k2Γh ≤ C(kξp k2 + kξq k2 ).
With the above lemma, we can proceed to the proof of the main theorem. We take w o = ξo in (3.1)-(3.6) to obtain 1 dkξu k2 + kξp k2 + kξq k2 = T1 − T2 − T3 , 2 dt 1 dkξv k2 + kξr k2 + kξs k2 = T4 , 2 dt
14
(3.19) (3.20)
where T1 = ((ηu )t , ξu ), T2 = (ru − rh uh , (ξu )x ) + hru − rd h uh , [ξu ]ix , T3 = (su − sh uh , (ξu )y ) + hsu − sd h uh , [ξu ]iy , T4 = ((ηv )t , ξv ) − (ηu − ξu − ηv + ξv , ξv ),
with hu, vix =
Ny Z N x −1 X X i=1 j=1
Jj
uvi+ 1 ,j dy 2
and hu, viy =
y −1 Z Nx N X X
i=1 j=1
Furthermore, we also define
kukΓxh = hu, uix
Ii
uvi,j+ 1 dx. 2
and kukΓyh = hu, uiy .
Then by (3.18) and (3.17), we have k[ξu ]kΓxh ≤ Ch1/2 kξp k and k[ξu ]kΓyh ≤ Ch1/2 kξq k.
(3.21)
Now we give the estimates of Ti′ s. By Cauchy-Schwarz inequality and Lemma 3.2 T1 ≤ k(ηu )t kkξu k ≤ Chk1 +1 kξu k.
(3.22)
The estimates of T2 and T3 are similar, and we consider T2 only. T2 = (r(u − uh ) + (r − rh )uh , (ξu )x ) + hru − rd h uh , [ξu ]ix , = T21 + T22 + T23 ,
where T21 = (r(u − uh ) + (r − rh )uh , (ξu )x ) T22 =
1 − − h2ru − rh+ u+ h − rh uh , [ξu ]ix 2
T23 = −hα[ηu − ξu ], [ξu ]ix , Using the fact that krk∞ ≤ C and the a priori error estimate kuh k∞ ≤ C, we have T21 ≤ C (ku − uh k + kr − rh k) k(ξu )x k 15
(3.23)
≤ C(kηu k + kξu k + kηr k + kξr k) kξp k ≤ C(hk1 +1 + kξu k + hk2 + kξr k) kξp k,
(3.24)
where the first step is based on Cauchy-Schwarz inequality, the second step follows from Lemma 3.5 and triangle inequality, and in the last one we use Lemma 3.2. The estimate of T22 also requires the boundedness of r and uh , T22 =
1 + + − − − hr(u − u+ h ) + (r − rh )uh + r(u − uh ) + (r − rh )uh , [ξu ]ix 2
≤ C (ku − uh kΓh + kr − rh kΓh ) k[ξu ]kΓxh ≤ Ch1/2 (kηu kΓh + kξu kΓh + kηr kΓh + kξr kΓh ) kξp k ≤ C(hk1 +1 + kξu k + hk2 + kξr k) kξp k.
(3.25)
Here in the second step we use Cauchy-Schwarz inequality, the third step follows from triangle inequality and (3.21), the last one requires Lemma 3.4 and Lemma 3.2. Now we proceed to the estimate of T23 , T23 ≤ C(k[ηu ]kΓh + k[ξu ]kΓxh ) k[ξu ]kΓxh ≤ Ch1/2 (k[ηu ]kΓh + k[ξu ]kΓh ) kξp k ≤ C(hk1 +1 + kξu k) kξp k,
(3.26)
where the first step follows from Cauchy-Schwarz inequality, the second step is based on (3.21), the third one requires Lemma 3.4 and Lemma 3.2. Plug (3.24), (3.25) and (3.26) into (3.23) to obtain T2 ≤ (Chk1 +1 + Chk2 + Ckξu k + Ckξr k) kξp k.
(3.27)
Following the same analysis of T2 , we obtain a similar estimate of T3 : T3 ≤ (Chk1 +1 + Chk2 + Ckξu k + Ckξs k) kξq k.
(3.28)
The estimate of T4 is simple, we only use Cauchy-Schwartz inequality and Lemma 3.2, T4 ≤ (k(ηv )t k + kηu k + kξu k + kηv k + kξv k) kξv k 16
≤ C hmin(k1 ,k2 )+1 + kξu k + kξv k kξv k
(3.29)
Plug (3.22), (3.27) and (3.28) into (3.19) to obtain 1 dkξu k2 + kξp k2 + kξq k2 ≤ Chk1 +1 kξu k + C(hmin(k1 +1,k2 ) + kξu k)(kξp k + kξq k) 2 dt +Ckξr kkξp k + Ckξs kkξq k ≤ C h2 min(k1 +1,k2 ) + kξu k2 + kξs k2 + kξr k2 + kξp k2 + kξq k2 , which further implies 1 dkξu k2 ≤ C h2 min(k1 +1,k2 ) + kξu k2 + kξs k2 + kξr k2 . 2 dt
(3.30)
Moreover, plug (3.29) into (3.20) to yield 1 dkξv k2 + kξr k2 + kξs k2 ≤ C hmin(k1 ,k2 )+1 + kξu k + kξv k kξv k. 2 dt
(3.31)
Combining (3.30) and (3.31), we can find a constant γ such that 1 dkξu k2 +γ 2 dt
1 dkξv k2 + kξr k2 + kξs k2 2 dt
≤ C h2 min(k1 +1,k2 ) + kξu k2 + kξs k2 + kξr k2 + γC hmin(k1 ,k2 )+1 + kξu k + kξv k kξv k Take γ to be sufficiently large, we have dkξu k2 dkξv k2 +γ ≤ Ch2 min(k1 +1,k2 ) + Ckξu k2 + Ckξv k2 . dt dt Therefore, by the Gronwall’s inequality and use L2 -projection as the initial discretization, kξu k2 + kξv k2 ≤ Ch2 min(k1 +1,k2 ) , which further yield ku − uh k + kv − vh k ≤ Chmin(k1 +1,k2 ) by Lemma 3.2.
17
3.6
Verification of the a priori error estimate
In this section, we proceed to verify the a priori error estimate (3.15). Actually, (3.15) is true at t = 0. Denote t⋆ = inf{t|k(u − uh )(t)k > h}. If t⋆ < T , then at t = t⋆ , we have h = ku − uh k ≤ Chmin(k1 +1,k2 ) , which is a contradiction if k1 ≥ 1, k2 ≥ 2 and h is small. Therefore, we have ku − uh k ≤ Chmin(k1 +1,k2 ) . Now we finish the whole proof.
3.7
Qk polynomials
In this section, we consider another finite element space and give the error estimates. If not otherwise stated, we follow the same notations used before. We first introduce the finite element space Whk defined as n o Whk = z : z Kij ∈ Qk (Kij ) ,
where Qk (Kij ) denotes the set of tensor product polynomials of degree up to k in cell Kij . In this section, we will take vh , rh , sh ∈ Whk and uh , ph , qh ∈ Vhk to obtain optimal error estimate. We first define several special projections. We define P+ into Whk which is, for each cell Kij , Z k−1 (P+ u − u, v)ij = 0, ∀v ∈ Q (Kij ), (P+ u − u)(xi− 1 , y)v(y)dy = 0, ∀v ∈ P k−1(Jj ), 2 Jj Z (P+ u − u)(xi− 1 , yj− 1 ) = 0, (P+ u − u)(x, yj− 1 )v(x)dx = 0, ∀v ∈ P k−1(Ii ). (3.32) 2
2
2
Ii
Moreover, we also define Πxk and Πyk into Whk which are, for each cell Kij , Z x k (Πk u − u, vx )ij = 0, ∀v ∈ Q (Kij ), (Πxk u − u)(xi+ 1 , y)v(y)dy = 0, ∀v ∈ P k (Jj ), 2 J Zj (3.33) y y k k (Πk u − u, vy )ij = 0, ∀v ∈ Q (Kij ), (Πk u − u)(x, yj+ 1 )v(x)dx = 0, ∀v ∈ P (Ii ). 2
Ii
Following the same analysis before, we define ηv = v − P+ v,
ηr = r − Πxk r,
ηs = s − Πyk s,
ξv = vh − P+ v,
ξr = rh − Πxk r,
ξs = sh − Πyk s.
and
Following the proof of Lemma 3.5, we obtain another similar one 18
Lemma 3.6. Suppose ξv , ξr and ξs are defined above, we have k(ξv )x k ≤ Ckξr k,
k(ξv )y k ≤ Ckξs k,
h−1 k[ξv ]k2Γxh ≤ Ckξr k2 ,
h−1 k[ξv ]k2Γy ≤ Ckξs k2 . h
Moreover, we will use the following lemma [4] Lemma 3.7. Suppose v is sufficiently smooth, then we have kηv k + hk∇ηv k + h1/2 kηv kΓh ≤Chk+1 ,
(3.34)
kηr k + hk∇ηr k + h1/2 kηr kΓh ≤Chk+1 ,
(3.35)
kηs k + hk∇ηs k + h1/2 kηs kΓh ≤Chk+1 .
(3.36)
Finally, we also need the following superconvergence property [7]. Lemma 3.8. For any v ∈ Whk , we have |Qx (ηu , v)| ≤ Chk+1 kukk+2kvk,
|Qy (ηu , v)| ≤ Chk+1 kukk+2kvk.
Now we can state the main theorem. Theorem 3.2. Suppose the exact solutions u, v are smooth and the derivatives r, s are uniformly bounded. The LDG scheme is defined as (2.1)-(2.6) with uh , qh , qh ∈ Vhk and vh , rh , sh ∈ Whk for k ≥ 1. Moreover, the initial discretization is also given as the L2 projection. Then we have ku − uh k + kv − vh k ≤ Chk+1 . Proof. In the proof, we skip the a posterior error estimate, its verification in Section 3.5 and most of the derivation steps, since they can be obtained from Section 3.5 with some minor changes. We take w o = ξo in (3.1)-(3.6) to obtain 1 dkξu k2 + kξp k2 + kξq k2 = T1 − T2 − T3 , 2 dt 1 dkξv k2 + kξr k2 + kξs k2 = T4 + T5 , 2 dt
19
(3.37) (3.38)
where T1 = ((ηu )t , ξu ), T2 = (ru − rh uh , (ξu )x ) + hru − rd h uh , [ξu ]ix , T3 = (su − sh uh , (ξu )y ) + hsu − sd h uh , [ξu ]iy , T4 = ((ηv )t , ξv ) − (ηu − ξu − ηv + ξv , ξv ), T5 = Qx (ηu , ξv ) + Qy (ηu , ξv ) Following the analyses in Section 3.5, we obtain the estimates of Ti′ s. T1 ≤ Chk1 +1 kξu k.
(3.39)
T2 ≤ (Chk+1 + Ckξu k + Ckξr k) kξp k
(3.40)
T3 ≤ (Chk+1 + Ckξu k + Ckξs k) kξq k T4 ≤ C hk+1 + kξu k + kξv k kξv k
(3.41) (3.42)
Finally, the estimate of T5 follows from Lemma 3.8 directly, T5 ≤ Chk+1 (kξr k + kξs k).
(3.43)
Plug (3.39), (3.40) and (3.41) into (3.37) to obtain 1 dkξu k2 + kξp k2 + kξq k2 ≤ Chk+1 kξu k + C(hk+1 + kξu k)(kξpk + kξq k) 2 dt +Ckξr kkξp k + Ckξs kkξq k
which further implies
≤ C hk+1 + kξu k2 + kξs k2 + kξr k2 + kξp k2 + kξq k2 , 1 dkξu k2 ≤ C h2k+2 + kξu k2 + kξs k2 + kξr k2 . 2 dt
Moreover, plug (3.42) and (3.43) into (3.38) to yield
1 dkξv k2 + kξr k2 + kξs k2 ≤ C hk+1 + kξu k + kξv k kξv k + Chk+1(kξr k + kξs k) 2 dt 20
(3.44)
which further implies
1 1 ≤ C hk+1 + kξu k + kξv k kξv k + Ch2k+2 + kξr k2 + kξs k2 , 2 2
dkξv k2 + kξr k2 + kξs k2 ≤ C h2k+2 + kξu k2 + kξv k2 dt
(3.45)
Combining (3.44) and (3.45), we can find a constant γ such that dkξu k2 dkξv k2 2 2 +γ + kξr k + kξs k dt dt ≤ C h2k+2 + kξu k2 + kξs k2 + kξr k2 + γC h2k+2 + kξu k2 + kξv k2 Take γ to be sufficiently large, then
dkξu k2 dkξv k2 +γ ≤ Ch2k+2 + Ckξu k2 + Ckξv k2 . dt dt Therefore, by the Gronwall’s inequality and L2 initial discretization kξu k2 + kξv k2 ≤ Ch2k+2 , which further yield ku − uh k + kv − vh k ≤ Chk+1 by Lemma 3.2 and Lemma 3.7.
4
Positivity preserving property
In this subsection, we apply the positivity-preserving technique to construct positive numerical approximations. Following the idea in [48], we restrict ourselves to the P 1-LDG formulation. For simplicity, we consider Euler forward time discretization and assume the numerical solutions at time level n are positive: unh > 0, vhn > 0. Then we will construct positive numerical approximations at time level n + 1. To do so, we will firstly prove that the cell average at time level n + 1 is positive. Next, we can use a simple limiter to modify the DG polynomial and make it to be positive without changing its cell average. Throughout this section, we use oij for the numerical approximation oh in Kij , and the cell average is o¯ij . 21
Let us consider how to construct positive uh first, and the equation satisfied by the cell averages is u¯n+1 = Hijc (u, r, s) + Hijd (u, p, q), ij where Hijc (u, r, s)
∆t 1 = u¯nij − 2 ∆x∆y
Hijd (u, p, q)
1 ∆t = u¯nij + 2 ∆x∆y
Z Jj
cri− 1 ,j u cri+ 1 ,j − u 2
2
! Z dy + u csi,j+ 1 − u csi,j− 1 dx , 2
Ii
2
! Z Z qˆi,j+ 1 − qˆi,j− 1 dx . pˆi+ 1 ,j − pˆi− 1 ,j dy + Jj
2
2
Ii
2
2
For simplicity, if not otherwise stated, we will omit the subscript ij in Hijc and Hijd . To analyze H c , we approximate the integral by a 2-point Gauss quadrature. The Gauss quadrature points on Ii and Jj are denoted by pxi
=
n
xβi
: β = 1, 2
o
and
pyj
n o β = yj : β = 1, 2 ,
respectively. Also, we denote wβ as the corresponding weights on the interval − 21 , 21 . Let λ1 =
∆t ∆x
and λ2 =
∆t , ∆y
then H c becomes
X X 1 n wβ u cr i− 1 ,β − u wβ u csβ,j− 1 − u cr i+ 1 ,β + λ2 csβ,j+ 1 , H (u, r, s) = uij + λ1 2 2 2 2 2 2
2
β=1
β=1
c
where u cri− 1 ,β is the numerical flux at (xi− 1 , yjβ ). Likewise for the other fluxes. As the general 2
2
treatment, we rewrite the cell average on the right hand side as
1X 1X − − + + u + u = wβ u + w u , = β i+ 21 ,β β,j+ 21 i− 21 ,β β,j− 21 2 β=1 2 β=1 2
unij
2
where u− = u− (y β ) is a point value of the numerical approximation in the Gauss i− 1 ,β i− 1 ,j j 2
2
quadrature. Likewise for the other point values. Define µ = λ1 + λ2 , then 2 X 1 + − c ui− 1 ,β + ui+ 1 ,β + u cri+ 1 ,β H (u, r, s) = λ1 wβ cri− 1 ,β − u 2 2 2 2 4µ β=1 2 X 1 + − uβ,j− 1 + uβ,j+ 1 + u csβ,j− 1 − u csβ,j+ 1 . + λ2 wβ 2 2 2 2 4µ β=1 22
We need to suitably choose the parameter α in the Lax-Friedrichs flux, and the result is given below. Lemma 4.1. We can choose α>
max 1 ≤ i ≤ Nx − 1 1 ≤ j ≤ Ny − 1 β = 1, 2
o n − + − + , −r , s , −s , ri+ 1 ,β i+ 1 ,β β,j+ 1 β,j+ 1 2
2
2
2
then un+1 > 0 under a CFL condition j ∆t 1 ∆t A ≤ , + ∆x ∆y 2 where A=α−
min 1 ≤ i ≤ Nx − 1 1 ≤ j ≤ Ny − 1 β = 1, 2
(4.1)
o n + − + − ≥ 0, ri+ , −r , s 1 1 1 , −s 1 ,β i+ ,β β,j+ β,j+ 2
2
2
2
Proof. It is easy to check, for i = 2, 3, · · · , Nx , 1 + cri− 1 ,β u 1 +u 2 4µ i− 2 ,β 1 − 1 + − + + + − + u r − α(u − u ) = ui− 1 ,β ri− ui− 1 ,β + 1 1 1 1 1 ,β i− 2 ,β i− 2 ,β i− 2 ,β i− 2 ,β 2 2 2 4µ 2 1 1 1 − + = u− + − α u+ . α + ri− + ri− 1 1 1 ,β i− ,β ,β i− 21 ,β 2 2 2 2 2 2µ − Therefore, we can choose α+ri− > 0 and 1 ,β 2
+ 1 −α +ri− 1 2µ ,β 2
> 0 to obtain
1 + +c uri− 1 ,β u 4µ i− 21 ,β 2
0. If i = 1, then by (2.12), 1 + 1 + cri− 1 ,β = ui− 1 ,β + u u 1 > 0. 2 2 4µ 4µ i− 2 ,β
Also, for i = 1, 2, · · · , Ny − 1,
1 − cri+ 1 ,β u 1 −u 2 4µ i+ 2 ,β 1 − 1 − − + + + − + u r − α(u − u ) ui+ 1 ,β ri+ ui+ 1 ,β − = 1 ,β i+ 21 ,β i+ 12 ,β i+ 21 ,β i+ 21 ,β 2 2 2 4µ 2 1 1 1 + − + = α − ri+ − ri+ + − α u− . u 1 1 1 ,β ,β i+ 12 ,β i+ 2 ,β 2 2 2 2 2µ 23
>
+ > 0 and We can choose α − ri+ 1 ,β 2
1 2µ
− − α > 0 such that − ri+ 1 ,β 2
1 − u 4µ i+ 21 ,β
i = Ny , then by (2.12), 1 − 1 − ui+ 1 ,β − u u 1 > 0. cri+ 1 ,β = 2 2 4µ 4µ i+ 2 ,β
> 0 and Similarly, for j = 2, 3, · · · , Ny we require α + s− β,j− 1 2
1 + +c usβ,j− 1 u 4µ β,j− 21 2
such that
> 0. For j = 1, 2, · · · , Ny −1, we need
1 − u 4µ β,j+ 21
1 2µ
α−s+ β,j+ 21
−u cri+ 1 ,β > 0. If 2
+ s+ − α > 0 to obtain β,j− 1 2
> 0 and
1 −α −s− 2µ β,j+ 21
>0
−u csβ,j+ 1 > 0. 2
Now, we proceed to analyze H d (p, q). Following the same analysis in [48] with some minor changes, we can show the following lemma. Lemma 4.2. Suppose unij > 0 for all i and j, then H d (u, p, q) > 0 under a CFL condition ∆t 1 ∆t + ≤ . 2 2 ∆x ∆y 20 Proof. Define H d (u, p, q) =
Z
+
1 H 2 (u, q), ∆x
where
Z
Λ1 + − ui− 1 ,j + ui+ 1 ,j − λ1 pbi− 1 ,j − pbi+ 1 ,j dy = 2 2 2 2 Jj 4Λ Z Λ2 + − d ui,j− 1 + ui,j+ 1 − λ2 qbi,j− 1 − qbi,j+ 1 dx, H2 (u, q) = 2 2 2 2 Ii 4Λ
H1d (u, p)
with Λ1 =
1 H 1 (u, p) ∆y
∆t , ∆x2
Λ2 =
∆t ∆y 2
and Λ = Λ1 + Λ2 . It is easy to check that
Z
1 + − + 4ui+ 1 ,j − 3ui+ 1 ,j − ui− 1 ,j dy, i = 1, · · · , Nx − 1, = 2 2 2 Jj Jj ∆x Z Z 1 − + + − dy = − 3u − u dx, j = 1, · · · , Ny − 1. 4u qi,j+ 1 i,j+ 21 i,j− 21 i,j+ 21 2 Ii ∆y Ii p− dy i+ 12 ,j
Plug the flux (2.7) into (4.3) to obtain 1. i = 1 H1d (u, p) 24
(4.2)
(4.3) (4.4)
Λ1 + Λ1 − − ui+ 1 ,j + u 1 + λ1 pi+ 1 ,j dy = 2 2 4Λ i− 2 ,j Jj 4Λ Z Λ1 Λ1 + + − = − Λ1 ui− 1 ,j + − 3Λ1 ui+ 1 ,j + 4Λ1 ui+ 1 ,j dy 2 2 2 4Λ 4Λ Jj Z
> 0. 2. 2 ≤ i ≤ Nx − 1
H1d (u, p) Z Λ1 + Λ1 − − − = ui+ 1 ,j + u 1 + λ1 pi+ 1 ,j − pi− 1 ,j dy 2 2 2 4Λ i− 2 ,j Jj 4Λ Z Λ1 Λ1 − + + − + − 5Λ1 ui− 1 ,j + − 3Λ1 ui+ 1 ,j + 4Λ1 ui+ 1 ,j dy = Λ1 ui− 3 ,j + 3Λ1 ui− 1 ,j + 2 2 2 2 2 4Λ 4Λ Jj
> 0.
3. i = Nx H1d (u, p) Z Λ1 − Λ1 + − = ui+ 1 ,j + u 1 − λ1 pi− 1 ,j dy 2 2 4Λ i− 2 ,j Jj 4Λ Z Λ1 Λ1 − − + + − 4Λ1 ui− 1 ,j + u 1 dy = Λ1 ui− 3 ,j + 3Λ1 ui− 1 ,j + 2 2 2 4Λ 4Λ i+ 2 ,j Jj
> 0.
We can analyze H2d in a similar way, so we skip the proof. Based on the above two lemmas, we have the following theorem. Theorem 4.1. Suppose unij > 0 for all i and j, then un+1 >0 ij under the CFL conditions (4.1) and (4.2). Now, let us proceed to analyze v, and the equation satisfied by the numerical cell averages is v n+1 = ij
1 n v ij + H d (v, r, s) + ∆t unij − v nij 2 25
d
= H (v, r, s) +
1 − ∆t v nij + ∆tunij . 2
Applying Lemma 4.2, it is easy to prove the following theorem. Theorem 4.2. Suppose vijn > 0 for all i and j, then v n+1 >0 ij provided ∆t ∆t 1 + ≤ , 2 2 ∆x ∆y 20 and 1 ∆t ≤ . 2 Based on Theorems 4.1 and 4.2, the numerical cell averages we obtained are positive. However, the numerical solutions unij and vijn might still be negative. Hence, we have to modify the numerical solutions while keeping the cell averages untouched. For simplicity, we discuss the modification of unij only and the procedure is given in the following steps. • Set up a small number ε = 10−13 . • If unij > ε, then proceed to the following steps. Otherwise, unij is identified as the approximation to vacuum, and we will take u enij = unij as the numerical solution and
skip the following steps.
• Modify the density first: Compute bij = If bij < ε, then take u enij as with
min unij (x, y).
(x,y)∈Kij
u enij = unij + θij unij − unij , unij − ε , θij = n uij − bij
and use u enij as the new numerical density unij . 26
Following [42], we can show the L1 -stability of the numerical scheme with the positivitypreserving limiter. Since unh is positive, we have kunh kL1
=
Z
Ω
unh (x)dx
=
Z
Ω
u0h (x)dx = ku0h kL1 ,
where kukL1 is the standard L1 -norm of u on Ω. This implies the L1 -stability of the scheme.
4.1
High order time discretizations
All the previous analyses are based on first-order Euler forward time discretization. We can also use strong stability preserving (SSP) high-order time discretizations to solve the ODE system wt = Lw. More details of these time discretizations can be found in [36, 35, 20]. In this paper, we use the third order SSP Runge-Kutta method [36] w(1) = wn + ∆tL(wn ), 3 n w + 4 1 = wn + 3
w(2) = wn+1
1 w(1) + ∆tL(w(1) ) , 4 2 w(2) + ∆tL(w(2) ) , 3
(4.5)
and the third order SSP multi-step method [35] w
n+1
16 n 12 11 n−3 n n−3 = w + ∆tL(w ) . (w + 3∆tL(w )) + 27 27 11
(4.6)
Since an SSP time discretization is a convex combination of Euler forward, by using the limiter mentioned in Section 4, the numerical solutions obtained from the full scheme are also positive.
5
Numerical experiments
In this section, we present numerical examples in two space dimensions to verify the theoretical analysis and the positivity-preserving property of the proposed method. If not otherwise stated, we consider the third-order RKDG method (k = 2).
27
Example 5.1. We solve the following problem on the domain Ω = [0, 2π] × [0, 2π] ut − div(∇u − u∇v) = fu (x, y, t), vt − △v = u − v + fv (x, y, t),
x ∈ Ω, t > 0 x ∈ Ω, t > 0,
(5.1)
subject to boundary condition (1.2). We choose suitable fu and fv such that the exact solutions are u(x, y, t) = exp(−t)(cos(x) + cos(y)),
v(x, y, t) = exp(−t)(cos(x) + cos(y)).
We can easily check that the source terms are fu (x, y, t) = −2 exp(−2t)(cos2 (x) + cos(x) cos(y) + cos2 (y) − 1),
fv (x, y, t) = 0.
In Table 5.1, we present the numerical results for the proposed method without the bound preserving limiter. From the table, we can observe optimal rate of convergence with Vhk finite element spaces.
N 10 20 40 80
k=1 k=2 2 L error order L error order 2.69e-1 – 1.32e-2 – 7.16e-2 1.91 1.66e-3 2.99 1.74e-2 2.04 2.07e-4 3.00 4.25e-3 2.03 2.59e-5 3.00 2
Table 5.1: Example 5.1: accuracy test at T = 0.1 for the second- and third-order LDG methods without the positivity-preserving limiter. k is the degree of polynomial and Nx = Ny = N.
Example 5.2. In this example, we solve (1.1) with the following initial condition. u0 = 840 exp(−84(x2 + y 2)),
v0 = 420 exp(−42(x2 + y 2 )).
We use second-order positivity-preserving LDG methods to solve and the numerical approximations at time t is given in figure 5.1. In [15], the authors demonstrated that the blow-up time should be approximately t = 1.21 × 10−4 . However, we can continue our numerical simulation to t = 2 × 10−4 , and the numerical approximation is given in figure 5.2. We also solve 28
Figure 5.1: Example 5.2: Numerical approximations of u at t = 6 × 10−5 (left) and t = 1.2 × 10−4 with positivity-preserving limiter for P 1 polynomials and N = 160. .
Figure 5.2: Example 5.2: Numerical approximations of u at t = 2.0 × 10−4 with positivitypreserving limiter for P 1 polynomials and N = 160. .
29
the problem without the positivity-preserving limiter. We find that at about t = 8 × 10−5 the numerical scheme yields negative u. Moreover, we also test the numerical blow-up time. For simplicity, we take Nx = Ny = N, and compute the L2 -norm numerical approximations at time t with N × N cells, defined as S(N, t). Define tb (N) = inf{t : S(2N, t) ∗ 1.05% ≥ S(N, t)}
(5.2)
as the numerical blow-up time. We anticipate that as we refine the mesh, the numerical blow-up time will converge to the exact value. However, due to the computational cost, we have to apply adaptive methods and refine the mesh in the vicinity of the blow-up point, and this work will be considered in the future. To verify our anticipation, we consider the
Figure 5.3: L2 -norm of the numerical approximation for Example 5.2 with different N’s. . following function 1 exp f (x, t) = 1−t
−x2 2(t − 1)2
It is easy to see that lim f (x, t) → δ(x)
t→1−
30
.
(5.3)
in the sense of distribution. We consider the interval [-1,1] and divided in into N uniform cells. Denote uh to be the L2 -projection of f (x, t) in each cell and S(N, t) to be the L2 -norm of uh over [-1,1] with different t′ s. We also compute the numerical blow-up time tb (N) by (5.2) and the results are given in Table 5.2. We can clearly see that as we refine the meshes, tb (N) converges to 1, the exact blow-up time, in first order accuracy. N Blow-up time Error
10 -
20 40 80 160 320 0.671 0.836 0.918 0.959 0.980 0.329 0.164 0.092 0.041 0.020
Table 5.2: Numerical blow-up time with different mesh sizes for (5.3).
6
Conclusion
In this paper, we develop LDG methods to the KS chemotaxis model. We improve the result given by [15], and give the optimal error estimate under special finite element spaces. A special positivity-preserving limiter is constructed to obtain physically relevant numerical approximations. Moreover, we also prove the L1 -stability of the LDG scheme with the limiter. Numerical experiments are given to demonstrate the good performance of the LDG scheme and the estimate of the numerical blow-up time. In the future, we plan to apply adaptive methods to calculate the blow-up time more accurately.
7
Acknowledgments
Xingjie Helen Li would like to thank the Shanghai Centre for Mathematics Science (SCMS), Fudan University, for support during her visit.
References [1] F. Bassi and S. Rebay, A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations, Journal of Compu31
tational Physics, 131 (1997), 267-279. [2] A. Chertock and A. Kurganov, A second-order positivity preserving central-upwind scheme for chemotaxis and haptotaxis models, Numerische Mathematik, 111 (2008), 169-205. [3] S. Childress and J. Percus, Nonlinear aspects of chemotaxis, Mathematical Biosciences, 56 (1981), 217-237. [4] P. Ciarlet, Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [5] B. Cockburn and B. Dong, An analysis of the minimal dissipation local discontinuous Galerkin method for convection-diffusion problems, Journal of Scientific Computing, 32 (2007), 233-262. [6] B. Cockburn, S. Hou and C.-W. Shu, The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case, Mathematics of Computation, 54 (1990), 545-581. [7] B. Cockburn, G. Kanschat, I. Perugia, and D. Schotzau, Superconvergence of the local discontinuous galerkin method for elliptic problems on cartesian grids, SIAM Journal on Numerical Analysis, 39 (2001), 264285. [8] B. Cockburn, S.-Y. Lin and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one-dimensional systems, Journal of Computational Physics, 84 (1989), 90-113. [9] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework, Mathematics of Computation, 52 (1989), 411-435.
32
[10] B. Cockburn and C.-W. Shu, The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems, Journal of Computational Physics, 141 (1998), 199-224. [11] B. Cockburn and C.-W. Shu, The local discontinuous Galerkin method for time dependent convection-diffusion systems, SIAM Journal on Numerical Analysis, 35 (1998), 2440-2463. [12] Y. Epshteyn, Discontinuous Galerkin methods for the chemotaxis and haptotaxis models, Journal of Computational and Applied Mathematics, 224 (2009), 168-181. [13] Y. Epshteyn, Upwind-difference potentials method for Patlak-Keller-Segel chemotaxis model, Journal of Scientific Computing, 53 (2012), 689-713. [14] Y. Epshteyn and A. Izmirlioglu, Fully discrete analysis of a discontinuous finite element method for the Keller-Segel Chemotaxis model, Journal of Scientific Computing, 40 (2009), 211-256. [15] Y. Epshteyn and A. Kurganov, New interior penalty discontinuous Galerkin methods for the Keller-Segel Chemotaxis model, SIAM Journal on Numerical Analysis, 47 (2008), 368-408. [16] F. Filbet, A finite volume scheme for the Patlak-Keller-Segel chemotaxis model, Numerische Mathematik, 104 (2006), 457-488. [17] F. Filbet, P. Lauren¸cot and B. Perthame, Derivation of hyperbolic models for chemosensitive movement, Journal of Mathematical Biology, 50 (2005), 189-207. [18] H. Gajewski and K. Zacharias, Global behaviour of a reaction-diffusion system modelling chemotaxis, Mathematische Nachrichten, 195 (1998), 77-114. [19] I.M. Gelfand, Some questions of analysis and differential equations, American Mathematical Society Translations, 26 (1963), 201-219. 33
[20] S. Gottlieb, C.-W. Shu and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Review, 43 (2001), 89-112. [21] L. Guo and Y. Yang, Positivity-preserving high-order local discontinuous Galerkin method for parabolic equations with blow-up solutions, Journal of Computational Physics, 289 (2015), 181-195. [22] J. Hakovec and C. Schmeiser, Stochastic particle approximation for measure valued solutions of the 2d Keller-Segel system, Journal of Statistical Physics, 135 (2009), 133-151. [23] M.A. Herrero, E. Medina and J.J.L. Velazquez, Finite-time aggregation into a single point in a reaction-diffusion system, Nonlinearity, 10 (1997), 1739-1754. [24] D. Horstman, From 1970 until now: The Keller-Segel model in chemotaxis and its consequences I, Jahresbericht DMV, 105 (2003), 103-165. [25] D. Horstmann, From 1970 until now: The Keller-Segel model in chemotaxis and its consequences II, Jahresbericht DMV, 106 (2004), 51-69. [26] A.E. Hurd and D.H. Sattinger, Questions of existence and uniqueness for hyperbolic equations with discontinuous coefficients, Transactions of the American Mathematical Society, 132 (1968), 159-174. [27] E.F. Keller and L.A. Segel, Initiation on slime mold aggregation viewed as instability, Journal of Theoretical Biology, 26 (1970), 399-415. [28] E.F. Keller and L.A. Segel, Traveling bands of chemotactic bacteria, a theoretical analysis, Journal of Theoretical Biology, 30 (1971), 235-248. [29] A. Marrocco, 2D simulation of chemotaxis bacteria aggregation, ESAIM: Mathematical Modelling and Numerical Analysis (M 2 AN), 37 (2003), 617-630.
34
[30] X. Meng, C.-W. Shu, Q. Zhang and B. Wu, Superconvergence of discontinuous Galerkin methods for scalar nonlinear conservation laws in one space dimension, SIAM Journal on Numerical Analysis, 50 (2012), 2336-2356. [31] C. Patlak, Random walk with persistence and external bias, The Bulletin of Mathematical Biophysics, 15 (1953), 311-338. [32] T. Qin, C.-W. Shu and Y. Yang, Bound-preserving discontinuous Galerkin methods for relativistic hydrodynamics, submitted to Journal of Computational Physics. [33] W.H. Reed and T.R. Hill, Triangular mesh methods for the Neutron transport equation, Los Alamos Scientific Laboratory Report LA-UR-73-479, Los Alamos, NM, 1973. [34] N. Saito, Conservative upwind finite-element method for a simplified Keller-Segel system modelling chemotaxis, IMA Journal of Numerical Analysis, 27 (2007), 332-365. [35] C.-W. Shu, Total-variation-diminishing time discretizations, SIAM Journal on Scientific and Statistical Computing, 9 (1988), 1073-1084. [36] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shockcapturing schemes, Journal of Computational Physics, 77 (1988), 439-471. [37] R. Strehl, A. Sokolov, D. Kuzmin, D. Horstmann, and S. Turek, A positivity-preserving finite element method for chemotaxis problems in 3D, Journal of Computational and Applied Mathematics, 239 (2013), 290-303. [38] R. Tyson, L.J. Stern and R.J. LeVeque, Fractional step methods applied to a chemotaxis model, Journal of Mathematical Biology, 41 (2000), 455-475. [39] H. Wang, C.-W. Shu and Q. Zhang, Stability and error estimates of local discontinuous Galerkin methods with implicit-explicit time-marching for advection-diffusion problems, SIAM Journal on Numerical Analysis, 53 (2015), 206-227.
35
[40] H. Wang, S. Wang, C.-W. Shu and Q. Zhang, Local discontinuous Galerkin methods with implicit-explicit time-marching for multi-dimensional convection-diffusion problems, ESAIM: Mathematical Modelling and Numerical Analysis (M 2 AN), to appear. [41] Y. Yang and C.-W. Shu, Discontinuous Galerkin method for hyperbolic equations involving δ-singularities: negative-order norm error estimates and applications, Numerische Mathematik, 124, (2013), 753-781. [42] Y. Yang, D. Wei and C.-W. Shu, Discontinuous Galerkin method for Krause’s consensus models and pressureless Euler equations, Journal of Computational Physics, 252 (2013), 109-127. [43] Q. Zhang and C.-W. Shu, Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws, SIAM Journal on Numerical Analysis, 42 (2004), 641-666. [44] Q. Zhang and C.-W. Shu, Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws, SIAM Journal on Numerical Analysis, 48 (2010), 1038-1063. [45] X. Zhang and C.-W. Shu, On maximum-principle-satisfying high order schemes for scalar conservation laws, Journal of Computational Physics, 229 (2010), 3091-3120. [46] X. Zhang and C.-W. Shu, On positivity preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes, Journal of Computational Physics, 229 (2010), 8918-8934. [47] X. Zhang and C.-W. Shu, Positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations with source terms, Journal of Computational Physics, 230 (2011), 1238-1248.
36
[48] Y. Zhang, X. Zhang and C.-W. Shu, Maximum-principle-satisfying second order discontinuous Galerkin schemes for convection-diffusion equations on triangular meshes, Journal of Computational Physics, 234 (2013), 295-316. [49] X. Zhao, Y. Yang and C. Syler, A positivity-preserving semi-implicit discontinuous Galerkin scheme for solving extended magnetohydrodynamics equations, Journal of Computational Physics, 278 (2014), 400-415.
37
