A LOCAL DISCONTINUOUS GALERKIN METHOD ... - Semantic Scholar
A LOCAL DISCONTINUOUS GALERKIN METHOD FOR THE BURGERS–POISSON EQUATION HAILIANG LIU AND NATTAPOL PLOYMAKLAM Abstract. In this work, we design, analyze and test a local discontinuous Galerkin method for solving the Burgers–Poisson equation. This model, proposed by Whitham [Linear and Nonlinear Waves, John Wiley & Sons, New York, 1974] as a simplified model for shallow water waves, admits conservation of both momentum and energy as two invariants. The proposed numerical method is high order accurate and preserves two invariants, hence producing solutions with satisfying long time behavior. The L2 -stability of the scheme for general solutions is a consequence of the energy preserving property. The optimal order of accuracy for polynomial elements of even degree is proven. The numerical tests for two types of solutions are provided to illustrate both accuracy and capability of the method.
1. Introduction In this paper, we are interested in numerical approximations to the Burgers-Poisson (BP) equation of the form � 2 � u (1a) ut + − φ = 0, 2 x (1b) φxx − φ = u.
The subscript t (or x, respectively) denotes the differentiation with respect to time variable t (or x). System (1) can be rewritten as a nonlocal equation (2)
ut + uux + [G ∗ u]x = 0
with the kernel G(x) = 12 e−|x| . This nonlocal model was found as a simplified shallow water model by Whitham [30] to approximate the model with a singular kernel � � � 1 tanh k 1/2 ikx G(x) = e dk. 2π R k
As for the BP equation, analysis of traveling waves in [14] showed that the equation features wave breaking in finite time, and the authors also proved the global existence of weak entropy solutions. In this work, we develop a local discontinuous Galerkin method (LDG) to solve this nonlinear BP equation. Our proposed scheme is high order accurate, and preserves two invariants of momentum and energy, hence producing solutions with satisfying long time behavior. The L2 stability of the scheme for general solutions is a consequence of the energy preserving property. In the context of water waves, one of the best known local models is probably the Korteweg de Vries (KdV)-equation, ut + uux + uxxx = 0. This equation possess soliton solutions’ coherent structures that interact nonlinearly among themselves, then reemerge, retaining their identity and showing particle-like scattering behavior. Date: September 07, 2013. 2000 Mathematics Subject Classification. 65M60, 65M12, 35Q53. Key words and phrases. Discontinuous Galerkin method, Burgers-Poisson equation, conservation, stability. 1
2
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
In shallow water wave theory, the nonlinear shallow water equations which neglect dispersion altogether lead to the finite time wave breaking. On the other hand the third order derivative term in the KdV equation will prevent this ever happening in its solutions. In reality, some water waves appear to break, if the wave height is above certain threshold. Therefore in [30] Whitham raised an intriguing question: what kind of mathematical equation can describe waves with breaking? He suggested equation (2) with the above two kernels; many competing models have since been suggested to capture one aspect or another of the classical water-wave problem, see e.g. [15, 2, 3, 19, 13, 18, 26, 12, 25]. One common feature of these models is the associated global invariants, infinitely many or finitely many. The BP equation has the following two invariants: � � (3a) u(0, x) dx = u(t, x) dx =: E1 (t), � � 2 (3b) u (0, x) dx = u2 (t, x) dx =: E2 (t). It is desirable to design stable and high order accurate numerical schemes which preserve two invariants for solving the BP equation. It is believed that numerical methods preserving more invariants are advantageous: besides the high accuracy of numerical solutions, an invariant preserving scheme can preserve good stability properties after long-time numerical integration. Much more effort has been devoted in this topic for different integrable PDEs recently [16, 28, 27, 5]. The goal of this paper is to develop a discontinuous Galerkin (DG) method to preserve both momentum and energy at the discrete setting. The DG method is a class of finite element methods using completely discontinuous piecewise-polynomials for the numerical solution and the test functions. It was first designed and has been successful for solving first order PDEs such as nonlinear conservation laws [29, 8, 9, 7, 11]. The local discontinuous Galerkin (LDG) method is an extension of the DG method for solving higher order PDEs. It was first designed for convection-diffusion equations [10], and has been extended to other higher order wave equations, including the KdV equation [35, 33, 22, 34] and the Camassa-Holm equation [31], see also the recent review paper [32] on the LDG methods for higher order PDEs. The idea of the LDG method is to rewrite higher order equations into a first order system, and then apply the DG method on the system. In contrast, the direct discontinuous Galerkin methods, proposed in [23, 24] primarily for diffusion equations, aimed at directly solving higher order PDEs by the DG discretization, see e.g., [1, 36] for energy preserving DG methods for KdV type equations, and [21] for the Degasperis-Procesi equation. In this work we propose an LDG method based on formulation (1), for which the second equation was rewritten into a first order system for applying the LDG discretization. In the algorithm we update the solution in two steps: (1) given u, obtain φ by solving (1b) with the LDG method; (2) with the obtained φ, update u by solving (1a) with a standard DG method using a conservative numerical flux so that the resulting scheme preserves two integrals E1 and E2 in smooth region. As for error estimates, we define a global projection dictated by the selected numerical flux and obtain the needed projection error, following the strategy of error estimates carried out in [20] for the DDG method to solve convection-diffusion equations. Through careful estimates using this global projection, we obtain the optimal order of accuracy for polynomial elements of even degree. This is confirmed by the numerical tests with k = 2, 4. Numerical tests also show that for k odd, only k-th order of accuracy is observed. Such an optimal error estimate only for k =even was also shown in [1], and numerically observed in [1, 36] for KdV type equations. The
LDG FOR BP EQUATION
3
main feature of the scheme presented in this work is its capability to produce wave solutions with satisfying long time behavior. The rest of the paper is organized as follows. In section 2, we formulate our LDG method with a class of numerical fluxes. In section 3, we show that the LDG method for solving the Poisson equation (1b) is well defined and stable, and the method is shown to conserve both momentum and energy for the given numerical fluxes. In section 4, we obtain the optimal order of error between the numerical solution and smooth solutions for the conservative scheme when using polynomial elements of even degree. Finally, in section 5, we present numerical examples to illustrate the capacity of the LDG scheme to preserve two invariants after a long-time simulation. 2. The discontinuous Galerkin method 2.1. LDG formulation. We develop a local discontinuous Galerkin (LDG) method for the BP equation subject to initial data u0 (x) and periodic boundary conditions. Let us partition the interval I = [0, L] into 0 = x1/2 , x3/2 , . . . , xN +1/2 = L to get N equal subintervals and denote � � each cell by Ij = [xj−1/2 , xj+1/2 ], j = 1, . . . , N . The center of the cell is xj = 12 xj−1/2 + xj+1/2 . The piecewise polynomial space Vhk is defined as the space of polynomials of degree up to k in each cell Ij , that is, Vhk = {v : v|Ij ∈ P k (Ij ), j = 1, 2, . . . , N }.
(4)
Note that functions in Vhk are allowed to have discontinuities across the interfaces. The solution of the DG method is denoted by uh , which belongs to the finite element space Vhk . We denote − the limit values of uh at xj+1/2 from the right and from the left by (uh )+ j+1/2 and (uh )j+1/2 respectively. Let ω be a piecewise smooth function, its jump across the cell interface be denoted + − by [ω] := ω + − ω − , and its average at the cell interface, ω +ω , be denote by {ω}. 2 To define the LDG method, we introduce an auxiliary variable p = φx and rewrite (1a)-(1b) as follows: � 2� u (5a) ut + − p = 0, 2 x (5b) p − φx = 0, px − φ − u = 0.
(5c)
Then, the scheme is defined as follows: find uh , ph , φh ∈ Vhk such that �
(6a) (6b)
�
� �2 � u u2h h � (uh )t ρ dx − ρx dx + ρ ∂Ij − ph ρ dx = 0, 2 Ij Ij 2 Ij � � � �h γ �∂I = 0, ph γ dx + φh γx dx − φ j −
(6c) (6d)
�
Ij
Ij
Ij
� ph qx dx + p�h q �∂Ij − �
�
Ij
(φh + uh )q dx = 0,
Ij
(uh − u)|t=0 v dx = 0,
for all test functions ρ, γ, q, v in the finite element space Vhk . The choice for numerical fluxes �2 , φ � , p� is given by u h
h
h
4
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
�2 = 1 �(u+ )2 + u+ u− + (u− )2 � , u h h h h h 3 + − �h = θφ + (1 − θ)φ , φ h h
(7a) (7b)
− p�h = (1 − θ)p+ h + θph ,
(7c)
where θ ∈ [1/2, 1]. Here, the numerical fluxes at the endpoints of I can be defined using − − + + 2 U1/2 := UN +1/2 and UN +1/2 := U1/2 where U represents uh , φh , or ph . The resulting LDG scheme (6) with the fluxes (7) is called LDG-C. �2 is needed in order to capture the entropy For discontinuous solutions, an entropy flux for u h solution. One well-known choice is the Lax-Friedrich flux of the form �2 = 1 �(u− )2 + (u+ )2 − σ(u+ − u− )� , σ = 2 max |u|, (8) u h h h h h + 2 u∈[u− h ,uh ] with which the resulting LDG scheme is called LDG-D. In practice, one may adopt an adaptive numerical flux � � 1 (u+ )2 + u+ u− + (u− )2 if |u+ − u− | < 10−2 h h h h � 2 3 (9) uh = 1 �(u− )2 + (u+ )2 − 2σ(u+ − u− )� otherwise. h h h h 2
Here, 10−2 may vary as long as it can serve as a shock detector. Notation. We use � · �m,Ω as the H m -norm over domain Ω, and | · |m,Ω as its semi-norm. For m = 0, we simply use � · �Ω to denote the L2 -norm over domain Ω. We also use the notation � · �∞,Ω to denote the L∞ norm over domain Ω. The domain Ω could be a computational cell Ij or a master domain Iˆ := [−1, 1]. If Ω is the whole domain, we do not specify the domain unless necessary. 3. Analytical properties of the scheme 3.1. Existence, uniqueness, and stability. In this section, we prove the existence, uniqueness, and stability of ph , φh obtained from (6b)-(6c) with numerical fluxes (7b)-(7c), given uh . Lemma 3.1. The numerical scheme (6b)-(6c) with the numerical flux (7b)-(7c) satisfies �ph �2 + �φh �2 ≤ �uh �2 .
(10)
Proof. We choose γ = ph and q = φh . Then (6b)-(6c) gives � � � 2 �h ph �∂I = 0, ph dx + φh (ph )x dx − φ j Ij
−
�
Ij
Ij
� ph (φh )x dx + p�h φ �∂Ij −
Subtracting the two equations above gives � � 2 � ph + φ2h dx Ij
=−
�
Ij
� �h ph �∂I − φh (ph )x dx + φ j
�
Ij
�
Ij
φ2h dx
=
�
uh φh dx. Ij
� ph (φh )x dx + p�h φ �∂Ij −
�
uh φh dx. Ij
LDG FOR BP EQUATION
5
Take summation over j and use the periodic boundary condition to get � � � � �� � 2 � 2 � ph + φh dx = − φh (ph )x dx − φh [ph ] − ph (φh )x dx 1 I
I
−
� j
j+ 2
j
(p�h [φ])j+ 1 − 2
�
uh φh dx.
I
The first four terms on the right-hand side can be simplified to � � �� � − (φh ph )x dx − φh [ph ] + p�h [φ] I
��
=
j
I
j
�h [ph ] − p�h [φ] [φh ph ] − φ
�
j+ 12
j+ 12
= 0,
because of the choice of numerical fluxes (7b)-(7c). Therefore, we have that � � � 2 � 1 1 2 ph + φh dx ≤ |uh φh | dx ≤ �uh �2 + �φh �2 , 2 2 I I
�
which proves (10).
Remark: The inequality (10) implies the stabiliy of ph and φh since their L2 norms are bounded by �uh �. Moreover, (6b)-(6c) can be written in terms of a finite dimensional invertible linear mapping. Therefore, the inequality (10) shows that if uh = 0, then so are ph and φh . This proves the uniqueness and existence of ph and φh for any given uh . 3.2. Discrete conservation laws. In this section, we look at the properties of the numerical solution uh that are analogous to (3a)-(3b). Theorem 3.2. The following relations hold for all t > 0: � L � L (11) uh (t, x)dx = uh (0, x)dx, (12)
�
0
L
0
u2h (t, x)dx
=
�
0
L
0
u2h (0, x)dx
+
� j
(2θ − 1)
�
0
t�
[φh ]2 + [ph ]2
�
j+ 12
dτ .
Hence the scheme is conservative for θ = 1/2, and the scheme is energy stable for 1/2 < θ ≤ 1. Proof. Because (6) holds for any test function in Vhk , we choose ρ = 1 and γ = 1 in (6a)-(6b), respectively, to obtain �
(uh )t dx + Ij
�2 � � u h � �� ∂Ij − φh ∂Ij = 0. 2
Take summation over all j and use the periodic boundary condition, we have d dt This proves (11).
�
uh dx = I
�2 � �u � h � �� ∂Ij + φh ∂Ij = 0, 2 j
6
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
Next, we choose the test functions ρ = uh , γ = −φh , and q = −ph in (6a)-(6c) to obtain � � � �2 � u u2h h � (13) (uh )t uh dx − (uh )x dx + uh ∂Ij − ph uh dx = 0, 2 Ij Ij 2 Ij � � � �h φh �∂I = 0, (14) − ph φh dx − φh (φh )x dx + φ j �
(15)
Ij
Ij
Ij
� ph (ph )x dx − p�h ph �∂Ij +
�
(φh + uh )ph dx = 0,
Ij
Integrating some terms out and adding the above three relations together, we get � � � � � � � �2 u u3h �� φ2h �� p2h �� h � (uh )t uh dx + uh − ∂Ij + φh φh − ∂Ij − p�h ph − ∂Ij = 0. 2 6 2 2 Ij
Summing the terms above for all j = 1, . . . , N and using the periodic boundary condition, we get � � � � � �2 2�� 2�� � u u3h �� φ p 1d 2 h h � �∂I + p�h ph − h �∂I uh dx = − uh − ∂Ij − φh φh − j j 2 dt 2 6 2 2 j � � � � �2 � u 1 3 1 2 h �h [φh ] − [φ ] = [uh ] − [uh ] + φ 2 6 2 h j+ 1 1 j 2 j+ 2 � � 1 − p�h [ph ] − [p2h ] 2 j+ 12 � �� �h − {φh } [φh ] + ({ph } − p�h ) [ph ] = φ j
=
�� j
1 θ− 2
�
which proves (12).
�
[φh ]2 + [ph ]2
�
� 4. Error estimations
In this section we estimate the error from approximating ph , φh , uh . We proceed by defining a global projection with some established properties to prove that the errors from approximating ph and φh can be controlled by the errors from approximating uh . Then, we show that the error from approximating uh is of optimal order. 4.1. The global projection. Let ω be a piecewise smooth function on I = [0, L], we define the projection Qθ such that it satisfies the following properties: � � (16a) (Qθ ω) v dx = ωv dx, ∀v ∈ P k−1 , j = 1, · · · , N, Ij
(16b) where
Ij
� Q �j+1/2 , θ ω j+ 1 = ω 2
j = 1, . . . , N,
V� := θV + + (1 − θ)V − .
LDG FOR BP EQUATION
7
For j = N , we use the periodic extension to define (Qθ ω)+ N +1/2 , in order to be consistent with the numerical flux defined in (7). We first show that the projection Qθ w is well defined. Lemma 4.1. The projection Qθ satisfying (16) is uniquely defined for either θ �= with k even and N odd.
1 2
or θ =
1 2
Proof. Let {ψl }kl=0 be a set of orthogonal Legendre polynomials on [−1, 1] of degree up to k. We can write the projection Qθ of ω ∈ H k+1 (I) on each cell Ij as k
(Qθ ω)(xj +
� j h ξ) = al ψl (ξ), 2
ξ ∈ [−1, 1].
l=0
With v = ψi , the condition (16a) gives � 2i + 1 1 h j (17) ai = ω(xj + ξ)ψi (ξ) dξ, i = 0, . . . , k − 1, 2 2 −1 �1 2 2 where we have used −1 ψi (ξ) dξ = 2i+1 . j It remains to determine ak for j = 1, . . . , N . Since ψl (±1) = (±1)l , the condition (16b) gives � k � � k � � j+1 � j (18) θ al (−1)l + (1 − θ) al = ω ˆ (xj+ 1 ), j = 1, . . . , N. l=0
2
l=0
Because ω is periodic, we require that k �
+1 aN ψl (ξ) = l
l=0
k �
a1l ψl (ξ),
l=0
∀ξ ∈ [−1, 1].
which allows us to write the system (18) as 1 1 − θ (−1)k θ 0 ··· 0 ak b1 2 0 1 − θ (−1)k θ · · · 0 a k b2 (19) .. .. · .. = .. , .. .. .. . . . . . . . k bN aN (−1) θ 0 0 ··· 1 − θ k �� � �� � k−1 j+1 k−1 j where bj = ω ˆ (xj+ 1 )−θ (−1)l −(1−θ) l=0 al l=0 al . The determinant of the coefficient 2
matrix A above is given by (1 − θ)N + (−1)N +1+kN θN , which is non-zero for all θ �= 12 . When θ = 12 , the determinant is non-zero whenever N is odd and k is even. This shows that the projection Qθ is unique whenever θ �= 12 . If θ = 12 , Qθ is unique when k is even and N is odd. � Lemma 4.2. We have the following projection error (20) where C∗ depends on k ≥ 1 and θ.
�Qθ ω − ω� ≤ C∗ hk+1 ,
Proof. The proof is carried out in two steps. Step 1. We first establish the following inequality (21)
�(Qθ − I)ω�2 ≤ Ch
N � j=1
�˜ ω j �2k+1,Iˆ,
8
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
where h ˆ ξ), ξ ∈ [−1, 1] = I. 2 By the Cauchy-Schwarz inequality, we have from (17) that for j = 1, . . . , N , ω ˜ j (ξ) := ω(xj +
(22)
2i + 1 j 2 �w ˜ �0,Iˆ, 2
|aji |2 ≤
(23)
i = 0, . . . , k − 1.
Hence N � k−1 � j=0 i=0
From (19) of the form
|aji |2 ≤ (k − 1/2)
N � j=1
�˜ ω j �20,Iˆ.
ak = A−1 b, T T where ak = [a1k , . . . , aN k ] and b = [b1 , . . . , bN ] , it follows that N � j=1
(ajk )2 = bT (A−1 )T A−1 b ≤ C ≤C ≤C
N � j=1
N � j=1
�
(˜ ω j (1))2 +
N �
(bj )2
j=1
k−1 �
(aj+1 )2 + l
l=0
k−1 �
(ajl )2
l=0
�
�˜ ω j �21,Iˆ,
where we have used the Sobolev inequality |˜ ω j |∞,Iˆ ≤ C�˜ ω j �1,Iˆ. Hence, 2
�Qθ ω� =
N � j=1
�Qθ ω�20,Ij
� k−1 � N � h� j 2 h j = (al ) �ψl �20,Iˆ + (ak )2 �ψk �20,Iˆ 2 2 j=1 l=0 � � N k−1 N � � � j 2 j 2 ≤h |ai | + (ak ) ≤ Ch �˜ ω j �21,Iˆ. j=1
i=0
j=1
Step 2. For any v ∈ Vhk (I), we have that Qθ v = v. Therefore, �Qθ ω − ω�2 = inf �(Qθ − I)(ω − v)�2 v∈Vhk
≤ Ch ≤ Ch which proves (20).
N �
inf
�˜ ω j − vˆ�2k+1,Iˆ
j=1
vˆ∈P k [−1,1]
N �
|˜ ω j |2k+1,Iˆ = Ch2k+2 |ω|2k+1 ,
j=1
�
LDG FOR BP EQUATION
9
Lemma 4.3. For k ≥ 1 the following inequality holds, �2 N � � � � − 2k+1 � � (24) . �(ω − Qθ ω)(xj+ 1 )� ≤ Ch 2
j=1
The constant C depends on k and θ.
Proof. On each interval Ij , using the orthogonality relation (16a), we have j
ω(x)|Ij := ω ˜ (ξ) = j
� Qθ ω(x)|Ij := (Q θ ω) (ξ) =
∞ � l=0 k=1 �
ωlj ψl (ξ), ωlj ψl (ξ) + αkj ψk (ξ).
l=0
Hence, by ψl (1) = 1, we have (25)
�∞ � � �2 � � �2 � � � � j − �(ω − Qθ ω)(x 1 )� ≤ 2 � ωl � + 2|αkj |2 . � j+ 2 � � � l=k
To control the first term on the right-hand side of (25), we consider the following expression j
(26)
∂ξ ω ˜ (ξ) =
∞ �
βlj ψl (ξ).
l=0
Following the idea in [4], we integrate (26) with respect to ξ to get � ξ ∞ � j j j ω ˜ (ξ) = ω ˜ (−1) + βl ψl (ν)dν. −1
l=0
Using the property of Legendre polynomials � ξ 1 ψi (ν)dν = (ψi+1 (ξ) − ψi−1 (ξ)) , 2i + 1 −1 we can write ω ˜ j (ξ) = ω ˜ j (−1) +
�
β0j
β1j − 3
�
ψ0 (ξ) +
∞ � l=1
�
j j βl−1 βl+1 − 2l − 1 2l + 3
Therefore, ωij =
�
j j βi−1 βi+1 − 2i − 1 2i + 3
�
,
i ≥ 1.
Thus, � j j β β k−1 k ωlj = + , 2k − 1 2k + 1 l=k � j+1 � ∞ � βk−1 βkj+1 j+1 l k ωl (−1) = (−1) − . 2k − 1 2k + 1 ∞ �
l=k
�
�
ψl (ξ).
10
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
These ensure the following estimates � �2 � � j j ∞ �� � 2(βk−1 )2 2(βk+1 )2 1 1 � j� (27a) ωl � ≤ + ≤ �∂ξ ω ˜ j �20,Iˆ � � � 2k − 1 2k − 1 2k + 1 2k − 1 l=k � �2 ∞ �� � 1 � j+1 l� (27b) ωl (−1) � ≤ �∂ξ ω ˜ j+1 �20,Iˆ. � � � 2k − 1 l=k
The second term on the right-hand side of (25) is determined by (16b), i.e., �∞ � �∞ � � j+1 � j j+1 j k l (28) θαk (−1) + (1 − θ)αk = θ ωl (−1) + (1 − θ) ωl , l=k
l=k
where we have used ω ˆ (xj+1/2 ) = θω ˜ j+1 (−1) + (1 − θ)˜ ω j (1). We then have from (27) and (28) that � �2 � �2 N N ∞ ∞ �� � �� � � � � � � �� |αkj |2 ≤ C ωlj+1 (−1)l � + � ωlj � � � � � j=1
j=1
≤C ≤C
N � j=1
N � j=1
l=k
l=k
�∂ξ ω ˜ j+1 �20,Iˆ + �∂ξ ω ˜ j �20,Iˆ �∂ξ ω ˜ j �20,Iˆ.
Insertion of these estimates back into (25) yields �2 N � N � � � � �(ω − Qθ ω)(x− 1 )� ≤ C �∂ξ ω ˜ j �20,Iˆ. � j+ 2 � j=1
j=1
Recall that Qθ v = v for any v ∈ P k , we proceed
�2 N � N � � � � − �(ω − Qθ ω)(x 1 )� ≤ C inf �∂ξ ω ˜ j − ∂ξ v˜�20,Iˆ � j+ 2 � v˜∈P k j=1
j=1
=C
≤C
N �
inf
�∂ξ ω ˜ j − p˜�2k,Iˆ
j=1
p˜∈P k−1
N �
|∂ξ ω ˜ j |2k,Iˆ
j=1
� �2k+2 � �−1 h h = C|ω|2k+1 ≤ C|ω|2k+1 h2k+1 . 2 2 � We will use the error estimates obtained in Lemmas 4.2-4.3 to estimate the error of the computed solution. We also need to use the following well-known inverse properties ( see e.g.,
LDG FOR BP EQUATION
11
[6]). �∂x w� ≤ Ch−1 �w�,
(29a)
�w�Γh ≤ Ch−1/2 �w�,
(29b)
�w�∞ ≤ Ch−1/2 �w�,
(29c)
for any w ∈ Vhk . The constant C is independent of w and h. The subscript Γh denotes the set of the boundary points of all the cells Ij . 4.2. Approximating p and φ. Lemma 4.4. Let (u, p, φ) be the exact solution of the system (5). Let (uh , ph , φh ) be obtained from (6) with the choice of fluxes (7). Let θ ∈ [0, 1] be such that both Qθ and Q1−θ are uniquely defined, then the following inequality holds for all t > 0. � � (30) �Q1−θ p − ph �2 + �Qθ φ − φh �2 ≤ 2 �Q1−θ p − p�2 + �Qθ φ − φ�2 + �u − uh �2 . Proof. Since the scheme with fluxes (7) is consistent, (6b)-(6c) also hold for (u, p, φ). In other words, � � � (31a) pγ + φγx − φγ �∂Ij = 0, Ij
−
(31b)
�
Ij
Ij
� pqx + pq �∂Ij −
�
(φ + u)q = 0,
Ij
Subtracting (6b)-(6c) from (31a)-(31b), we get the error equations � � � � � �h γ �∂I = 0, (32a) (p − ph )γ + (φ − φh )γx − φ − φ j Ij
−
(32b)
�
Ij
Ij
� (p − ph )qx + (p − p�h ) q �∂Ij −
�
Ij
(φ − φh )q =
�
Ij
(u − uh )q.
� � Define �p = Q1−θ p − p, wp = Q1−θ p − ph , ��p = Q �p = Q 1−θ p − p, and w 1−θ p − p�h . (Similar definition can be given for �φ , wφ , ��φ , and w �φ associated with Qθ .) We then choose γ = wp , q = wφ and take the summation of (32) over j to get �� �� � (wp − �p )wp + (wφ − �φ )(wp )x + (w �φ − ��φ ) [wp ]j+ 1 = 0, j
−
Ij
�� j
=
�
I
j
Ij
(wp − �p )(wφ )x −
(u − uh )wφ .
Ij
j
� j
(� wp − ��p ) [wφ ]j+ 1 − 2
2
�� j
Ij
(wφ − �φ )wφ
Take the difference of both equations, we get � � � � � 2 2 wp + wφ = �p wp + �φ wφ + �1 + �2 + �3 − (u − uh )wφ , I
I
I
I
I
12
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
where �1 = − �2 =
j
�� j
�3 =
��
�� j
Ij
(wφ (wp )x + wp (wφ )x ) −
�φ (wp )x + Ij
� j
�p (wφ )x + Ij
� j
� j
��φ [wp ]j+ 1 ,
w �φ [wp ]j+ 1 − 2
� j
w �p [wφ ]j+ 1 , 2
2
��p [wφ ]j+ 1 . 2
First, note that wφ (wp )x + wp (wφ )x = (wp wφ )x , so � � � �1 = [wp wφ ]j+ 1 − w �φ [wp ]j+ 1 − w �p [wφ ]j+ 1 = 0, 2
j
2
j
2
j
with the choice of numerical fluxes (7b)-(7c). As for �2 , the property (16a) of Qθ gives �� �� �φ (wp )x = (Qθ φ − φ)(wp )x = 0, Ij
j
j
Ij
since (wp )x is in P k−1 . On the other hand, the propery (16b) of Qθ gives � � �� � ��φ [wp ]j+ 1 = Q θ φ − φ [wp ]j+ 1 = 0. j
2
j
2
Similarly, the term �3 vanishes by the properties (16) of Q1−θ . Using the Young’s inequality together with the results from �1 , �2 , and �3 , we get (33)
1 1 1 �wp �2 + �wφ �2 ≤ ��p �2 + ��φ �2 + �u − uh �2 , 2 2 2
which proves (30).
�
4.2.1. Approximating u. Theorem 4.5. If k is even, then the numerical solution, uh , obtained from the scheme (6) and the numerical fluxes (7) with θ = 12 satisfies (34)
�u − uh � ≤ Chk+1 .
The constant C is independent of the mesh size. Proof. Since the scheme (6) with fluxes (7) is consistent, (6a) also holds for (u, p, φ). In other words, � � � �2 � u2 u � ut ρ − ρx + ρ ∂Ij − (35) pρ = 0. 2 Ij Ij 2 Ij
Define w = Q1/2 u − uh and � = Q1/2 u − u. We have that u − uh = w − �. Subtracting (6a) from (35) and choose ρ = w, we get � � � � � � � 2 � �2 � u2h u u2 u h � wt w = �t w + − wx − − w ∂Ij + ep w, 2 2 2 2 Ij Ij Ij Ij where ep = p − ph .
LDG FOR BP EQUATION
13
Take summation over all j and introduce {uh }2 /2 into the third term on the right-hand side to get � � � � � u2 u2 � � � u2 {uh }2 � h wt w = �t w + − wx + − [w]j+ 1 2 2 2 2 2 I I Ij j j � � � �2 � {uh }2 u h + − [w]j+ 1 + ep w. 2 2 2 I j
Using the identity A2 /2 − B 2 /2 = A(A − B) − (A − B)2 /2, we get � � � �� 1� wt w = �t w + u (u − uh ) wx − (u − uh )2 wx 2 I I Ij Ij j
+
� j
+
� j
j
u (u − {uh }) [w]j+ 1 − 2
�
�2 {uh }2 u − h 2 2
�
1� 2
[w]j+ 1 + 2
j
�
(u − {uh })2 [w]j+ 1 2
ep w. I
Let {w} = {Q1/2 u} − {uh } and {�} = {Q1/2 u} − u, we can write � � � (36) wt w = �t w + τ1 + τ2 + τ3 + τ4 + τ5 + ep w, I
I
I
where
τ1 =
�� j
τ2 = −
uwwx + Ij
�� j
�
u{w}[w]j+ 1 ,
j
Ij
��
u�wx −
� j
2
u{�}[w]j+ 1 , 2
1� {w}2 [w]j+ 1 , 2 2 I j j j � �� � 1� 1� τ4 = w�wx − �2 wx + {w}{�}[w]j+ 1 − {�}2 [w]j+ 1 , 2 2 2 2 Ij Ij j j j j � � �2 � {uh }2 u h τ5 = − [w]j+ 1 . 2 2 2 τ3 = −
1 2
w 2 wx −
j
Note that
� � w2 τ1 = u + u{w}[w]j+ 1 2 2 x Ij j j � 2� � � � w2 � �� w =− u − ux + u{w}[w]j+ 1 2 2 j+ 1 2 Ij j j j 2 � � � � w2 1 =− ux ≤ �ux �∞ �w�2 . 2 2 Ij ��
j
�
14
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
As for τ2 , we write u(x) = u(xj ) + u� (x∗ )(x − xj ) for all x ∈ Ij where x∗ is between x and xj . Therefore, � � � � � �� � � � � � � � � � � � � � ∗ �− � � u�wx − u{�}[w]j+ 1 � ≤ �− u(xj ) �wx − u (x ) �(x − xj )wx �� + |0| � 2 Ij Ij � j Ij � � j � j j � � � � � � �� � �� � ∗ � � � � � ≤� u(xj ) �wx � + u (x ) |�hwx | Ij Ij � j � j � � 1 ≤ |0| + �ux �∞ ���2 + h2 �wx �2 2 ≤ C(h2k+2 + �w�2 ), because of the projection properties (16), the inverse property (29a), and Lemma 4.2. For τ3 , we can show that � � � � 3� 1� w 1� 1 � w3 1� 2 − − {w} [w]j+ 1 = − {w}2 [w]j+ 1 2 2 2 3 x 2 2 3 2 Ij j
j
j
j
1 � 3 = [w]j+ 1 ≤ Ch−1 �w�∞ �w�2 ≤ Ch−3/2 �w�3 , 2 24 j
by the inverse properties (29b)-(29c). From the inverse properties (29a) and (29c), it follows that �wx �∞ ≤ Ch−3/2 �w�, with which we are able to estimate terms in τ4 by using Lemma 4.2, �� w�wx ≤ C�w�����wx �∞ ≤ Chk+1−3/2 �w�2 = Chk−1/2 �w�2 , j
Ij
� 1� 2
j
1 �2 wx ≤ �wx �∞ ���2 ≤ Ch2k+1/2 �w�. 2 Ij
As for the remaining terms in τ4 , � 1� {w}{�}[w]j+ 1 − {�}2 [w]j+ 1 = 0 2 2 2 j
j
because of the projection property (16b). �2 /2 = −[u ]2 /24, and [u ] = [u − u ] = [w] − [�], we Finally, using the fact that {uh }2 /2 − u h h h h have � � �2 � {uh }2 u � 1 h τ5 = − [w]j+ 1 = − [uh ]2 [w]j+ 1 , 2 2 2 2 24 j
j
=
� j
−
1 1 1 [w]3j+ 1 + [�][w]2j+ 1 − [�]2 [w]j+ 1 2 2 2 24 12 24
≤ C�w�∞ (�w�2Γh + ���Γh �w�Γh + ���2Γh )
≤ Ch−1/2 �w�(h−1 �w�2 + hk+1/2 h−1/2 �w� + h2k+1 ) ≤ C(h2k+2 + �w�2 + h−3/2 �w�3 ),
LDG FOR BP EQUATION
15
where we have used the inverse properties (29b) and (29c). The results from τ1 to τ5 and (36) give � � d (37) �w�2 ≤ C h2k+2 + �w�2 + h−3/2 �w�3 . dt With this, we can write �w(t, ·)�2 ≤ CG(t),
(38) where (39)
G(t) = h h2k+2
and
F (η) =
�
This gives G(0) =
(41)
η 1
+
�
t 0
�w(τ, ·)�2 + h−3/2 �w(τ, ·)�3 dτ .
� � G� (t) ≤ C∗ G(t) + h−3/2 G(t)3/2 .
(40) Now, define
2k+2
1 � dξ = −3/2 ξ+h G(0)ξ 3/2
�
η 1
1 ξ+
hk−1/2 ξ 3/2
dξ.
For k ≥ 1, we have that F � (η) is uniformly (with respect to h) positive and bounded above by 1 for all η > 1. Thus, there exists C˜ = F −1 (C∗ T ) for given T > 0. Integrate (40) to get � � G ˜ (42) F ≤ C∗ T = F (C). G(0) G ˜ Using this and (38), we prove (34) as desired. Hence ≤ C. G(0) � 4.3. Time discretization. We partition the time interval [0, T ] into M equal subintervals with boundaries {tn }, n = 0, 1, 2, . . . , M . Set ∆t = T /M as the time step. In order to preserve both mass and energy at the fully discrete level, we may use the Crank-Nicolson time discretization to find un+1 = 2u∗h − unh , h where u∗h is determined by �� �2 � � � � u∗h u∗h − unh (u∗h )2 2 ρ− ρx + ρ �∂Ij − p∗h ρ = 0, ∆t 2 2 I Ij Ij � j � � �∗ γ �∂I = 0, (43) p∗h γ + φ∗h γx − φ j h Ij
−
�
Ij
Ij
p∗h qx
� + p�∗h q �∂Ij −
�
Ij
(φ∗h + u∗h ) q = 0.
Indeed, this time discretization has the desired and provable properties. Theorem 4.6. The fully-discrete scheme (43) gives solution unh that satisfies � L � L n (44) uh dx = u0h dx, (45) for all 0 ≤ n ≤ M .
�
0
L
0
(unh )2 dx =
�
0
0
L�
u0h
�2
dx,
16
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
Proof. Take the test functions ρ = 1, γ = 1 in (43), adding together, to get �� �2 � � u∗h �� un+1 − unh h �∗ � + ∂Ij − φh ∂Ij = 0, ∆t 2 Ij
which upon summation over j proves (44). Next, we choose the test functions ρ = u∗h , γ = −φ∗h , and q = −p∗h so that �� �2 � � � u∗h un+1 − unh ∗ (u∗h )2 ∗ h uh − (uh )x + u∗h �∂Ij ∆t 2 2 Ij Ij � − p∗h u∗h = 0, �
− Ij
�
Ij
p∗h φ∗h −
p∗h (p∗h )x
�
Ij
Ij
� �∗ φ∗ �∂I = 0, φ∗h (φ∗h )x + φ j h h
� − p�∗h p∗h �∂Ij +
�
Ij
(φ∗h + u∗h )p∗h = 0.
Summation of the above three equations over j gives � � � � � � un+1 2 − (un )2 � L un+1 2 − (un )2 h h h h = = 0, 2∆t 2∆t Ij 0 j
which leads to (45).
�
Note that the above time discretization is fully nonlinear and requires the costly iteration solver. In practice, one would prefer to use some explicit solver with high order accuracy for time discretization. In our numerical simulation we choose to use the TVD third-order RungeKutta method [17] to solve the ODE system of the form a˙ = L(a): a(1) = an + ∆tL(an ), 3 1 (46) a(2) = an + a(1) + 4 4 1 2 an+1 = an + a(2) + 3 3 n n where a is the coefficient vector of uh .
1 ∆tL(a(1) ), 4 2 ∆tL(a(2) ), 3
5. Numerical Tests It is known [14] that one steady solution of system (1a)-(1b) is given by � 4 � −|x|/2 (47) U1 (x) = e −1 . 3 The system also has a steady periodic solution of the form � � � � 4 cosh x2 � � −1 , (48) U2 (x) = 3 cosh p2
for −p < x < p and by periodic continuation with period 2p. Because the system (1a)-(1b) is Galilean invariant, a family of traveling-wave solutions (1a)-(1b) may be obtained from the steady solutions as (49)
u(t, x) = U (x − u0 t) + u0 ,
LDG FOR BP EQUATION
17
where U is the steady state solution (47) or (48). We will use both steady and traveling wave solutions to test our scheme. Example 1. (Accuracy test) We run the semi-discrete scheme (6) and the numerical flux (7) with θ = 12 , 1, along with the third order Runge-Kutta method (46) on the steady state problem which has (48) as its exact solution. The results for k = 1, 2, 3, 4 are given in the two tables below. Here, we use ∆t = 0.001, final time tmax = 2, and p = 2. The norms of the error were computed by using the sixteen-point Guass quadrature rule. k N 1 10 20 40 80 2 10 20 40 80 3 10 20 40 80 4 10 20 40 80 k N 1
5 10 20 40 80 2 5 10 20 40 80 3 5 10 20 40 80 4 5 10 20 40 80
L1 5.5907e-02 2.8223e-02 1.4137e-02 7.0694e-03 1.6478e-03 2.0461e-04 2.5457e-05 3.1778e-06 2.5013e-04 3.1133e-05 3.8845e-06 4.8514e-07 1.2182e-06 3.7401e-08 1.1545e-09 3.5841e-11 L1 1.3616e-01 5.9425e-02 2.7403e-02 1.3131e-02 6.4241e-03 2.5044e-03 2.2204e-04 1.8941e-05 1.5556e-06 1.3670e-07 2.0147e-04 1.9490e-05 2.0788e-06 2.3547e-07 2.7898e-08 4.0456e-06 9.1182e-08 2.0673e-09 4.5009e-11 1.0746e-12
order 0.9862 0.9974 0.9998 3.0096 3.0067 3.0019 3.0062 3.0026 3.0013 5.0256 5.0178 5.0095 order 1.1962 1.1167 1.0614 1.0314 3.4956 3.5513 3.6059 3.5084 3.3697 3.2289 3.1421 3.0773 5.4715 5.4630 5.5214 5.3883
θ = 1/2 L2 order 3.3360e-02 1.6765e-02 0.9926 8.3909e-03 0.9986 4.1959e-03 0.9999 1.0558e-03 1.3282e-04 2.9908 1.6615e-05 2.9990 2.0768e-06 3.0000 1.7960e-04 2.2535e-05 2.9945 2.8200e-06 2.9984 3.5262e-07 2.9995 9.6189e-07 3.0251e-08 4.9908 9.4598e-10 4.9990 2.9535e-11 5.0013
L∞ 4.8122e-02 2.6199e-02 1.3877e-02 7.1781e-03 2.6115e-03 3.8008e-04 5.3250e-05 7.2587e-06 6.4039e-04 9.9254e-05 1.4970e-05 2.2034e-06 4.5036e-06 1.8170e-07 7.1522e-09 2.7574e-10
θ=1 L2 order 8.2323e-02 3.6667e-02 1.1668 1.7300e-02 1.0837 8.4184e-03 1.0391 4.1554e-03 1.0186 1.4860e-03 1.4181e-04 3.3894 1.3206e-05 3.4246 1.2308e-06 3.4235 1.1673e-07 3.3984 1.2239e-04 1.2763e-05 3.2614 1.4660e-06 3.1221 1.7711e-07 3.0492 2.1807e-08 3.0218 2.5910e-06 5.7074e-08 5.5045 1.3796e-09 5.3705 3.1612e-11 5.4476 7.7906e-13 5.3426
L∞ 1.1532e-01 5.2580e-02 2.3941e-02 1.1365e-02 5.6409e-03 2.4239e-03 3.1615e-04 3.8885e-05 4.7536e-06 5.8323e-07 1.8214e-04 2.8368e-05 4.6833e-06 7.2428e-07 1.0497e-07 7.9819e-06 1.3149e-07 4.7551e-09 1.7362e-10 5.7253e-12
order 0.8772 0.9168 0.9510 2.7805 2.8355 2.8750 2.6897 2.7290 2.7643 4.6314 4.6670 4.6970 order 1.1330 1.1350 1.0749 1.0106 2.9386 3.0233 3.0321 3.0269 2.6827 2.5987 2.6929 2.7866 5.9238 4.7893 4.7754 4.9225
18
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
The results show that the optimal order of accuracy is achieved only when k =even, which is consistent with our theoretical result on the optimal error estimates for k = even. Also such an observation seems unaffected by the choice of θ ∈ [0, 1], though we only display results for θ = 1/2 and θ = 1. Example 2. (Energy-preserving test) We compare the performance of the LDG-C and LDG-D schemes on the traveling wave version of (47), with u0 = 1/100. We use θ = 12 , 1 along with the third order Runge-Kutta method (46). The relative L2 energy �uh (t, x) − uh (0, x)� of the computed solution for k = 4 are shown in Figure 1. Here, we use ∆t = 0.001, tmax = 100. The L2 norms were computed by using the sixteen-point Gauss quadrature rule. In next two examples, we use polynomial elements of degree k = 2 along with the TVBM limiter introduced in [9]. Here, we use the mesh size h = 1/16. Example 3. We test the conservative scheme for initial data 0 ≤ x ≤ 1, 0.5 −x + 1.5, 1 ≤ x ≤ 4, u0 (x) = −2.5, x ≥ 4.
This initial data has a downward ramp of height 3 and the constant states lying symmetric with respect to u = −1, the solution is expected to converge to a stationary solution. In Figure 2, we observe a stable pattern formation as analyzed in [14]. In our experiment we use a modified initial data in C 2 , which agrees with the original data everywhere except for near x = 1, 4, so that we can apply directly the TVBM limiter introduced in [9]. Our goal is to observe the stable wave pattern, so the choice of modification is not essential. Example 4. We consider another initial data of the form x ≤ 8, −0.5 15.5 − 2x, 8 ≤ x ≤ 8.5, u0 (x) = −1.5, x ≥ 8.5. This example with smaller jump has no stable stationary solution to converge. We plot the computed solution at the different times in Figure 3, from which we can see that dispersive effects with oscillations propagate to the left of the ramp as analyzed in [14]. Acknowledgments This research was partially supported by the National Science Foundation under grant DMS1312636 and the KI-Net research network. References [1] J. L. Bona, H. Chen, O. Karakashian, and Y. Xing. Conservative, discontinuous Galerkin methods for the generalized Korteweg-de Vries equation. Math. Comp., 82(283):1401–1432, 2013. [2] R. Camassa and D. D. Holm. An integrable shallow water equation with peaked solitons. Phys. Rev. Lett., 71(11):1661–1664, 1993. [3] R. Camassa, D. D. Holm, and J. M. Hyman. A new integrable shallow water equation. Advances in Applied Mechanics, 31(31):1–33, 1994. [4] P. Castillo, B. Cockburn, D. Sch¨ otzau, and C. Schwab. Optimal a priori error estimates for the hp-version of the local discontinuous Galerkin method for convection-diffusion problems. Math. Comp., 71(238):455–478, 2002. [5] E. Celledoni, V. Grimm, R. I. McLachlan, D. I. McLaren, D. O’Neale, B. Owren, and G. R. W. Quispel. Preserving energy resp. dissipation in numerical PDEs using the “average vector field” method. J. Comput. Phys., 231(20):6770–6789, 2012.
LDG FOR BP EQUATION
19
Figure 1. Example 2: the evolution of the relative L2 energy over long-time simulation. [6] P. G. Ciarlet. The Finite Element Method for Elliptic Problems. North-Holland, New York, 1978. [7] 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. Math. Comp., 54(190):545–581, 1990. [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. J. Comput. Phys., 84(1):90–113, 1989. [9] B. Cockburn and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework. Math. Comp., 52(186):411–435, 1989. [10] B. Cockburn and C.-W. Shu. The local discontinuous Galerkin method for time-dependent convectiondiffusion systems. SIAM J. Numer. Anal., 35(6):2440–2463 (electronic), 1998. [11] B. Cockburn and C.-W. Shu. The Runge-Kutta discontinuous Galerkin method for conservation laws. V. Multidimensional systems. J. Comput. Phys., 141(2):199–224, 1998. [12] A. Constantin and D. Lannes. The hydrodynamical relevance of the Camassa-Holm and Degasperis-Procesi equations. Arch. Ration. Mech. Anal., 192(1):165–186, 2009. [13] A. Degasperis, D. D. Kholm, and A. N. I. Khon. A new integrable equation with peakon solutions. Teoret. Mat. Fiz., 133(2):170–183, 2002. [14] K. Fellner and C. Schmeiser. Burgers-Poisson: a nonlinear dispersive model equation. SIAM J. Appl. Math., 64(5):1509–1525 (electronic), 2004. [15] B. Fuchssteiner and A. S. Fokas. Symplectic structures, their B¨ acklund transformations and hereditary symmetries. Phys. D, 4(1):47–66, 1981/82. [16] D. Furihata and M. Mori. General derivation of finite difference schemes by means of a discrete variation. Transactions-Japan Society for Industrial and Applied Mathematics, 8:317–340, 1998. [17] S. Gottlieb and C.-W. Shu. Total variation diminishing Runge-Kutta schemes. Math. Comp., 67(221):73–85, 1998. [18] D. D. Holm and M. F. Staley. Wave structure and nonlinear balances in a family of evolutionary PDEs. SIAM J. Appl. Dyn. Syst., 2(3):323–380 (electronic), 2003.
20
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
Figure 2. Example 3: the computed solution at t = 0, 10, 100. [19] R. S. Johnson and Camassa-Holm. Korteweg-de Vries and related models for water waves. J. Fluid Mech., 455:63–82, 2002. [20] H. Liu. Optimal error estimates of the direct discontinuous Galerkin method for convection-diffusion equations. Math. Comp., 2013. in review. [21] H. Liu, Y. Huang, and N. Yi. A direct discontinuous Galerkin method for the Degasperis–Procesi equation. Methods Appl. Anal., 2013. in review. [22] H. Liu and J. Yan. A local discontinuous Galerkin method for the Korteweg–de Vries equation with boundary effect. J. Comput. Phys., 215(1):197–218, 2006. [23] H. Liu and J. Yan. The direct discontinuous Galerkin (DDG) methods for diffusion problems. SIAM J. Numer. Anal., 47(1):675–698, 2008/09. [24] H. Liu and J. Yan. The direct discontinuous Galerkin (DDG) method for diffusion with interface corrections. Commun. Comput. Phys., 8(3):541–564, 2010. [25] H. Liu and Z. Yin. Global regularity, and wave breaking phenomena in a class of nonlocal dispersive equations. In Nonlinear partial differential equations and hyperbolic wave phenomena, volume 526 of Contemp. Math., pages 273–294. Amer. Math. Soc., Providence, RI, 2010. [26] Y. Liu and Z. Yin. Global existence and blow-up phenomena for the Degasperis-Procesi equation. Comm. Math. Phys., 267(3):801–820, 2006. [27] T. Matsuo. Dissipative/conservative Galerkin method using discrete partial derivatives for nonlinear evolution equations. J. Comput. Appl. Math., 218(2):506–521, 2008. [28] R. I. McLachlan, G. R. W. Quispel, and N. Robidoux. Geometric integration using discrete gradients. R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci., 357(1754):1021–1045, 1999. [29] W. H. Reed and T. R. Hill. Triangular mesh methods for the neutron transport equation. Los Alamos Report LA-UR-73-479, 1973. [30] G. B. Whitham. Linear and Nonlinear Waves. John Wiley & Sons, New York, 1974. [31] Y. Xu and C.-W. Shu. A local discontinuous Galerkin method for the Camassa-Holm equation. SIAM J. Numer. Anal., 46(4):1998–2021, 2008.
LDG FOR BP EQUATION
21
Figure 3. Example 4: the computed solution at t = 0, 5, 20. [32] Y. Xu and C.-W. Shu. Local discontinuous Galerkin methods for high-order time-dependent partial differential equations. Commun. Comput. Phys., 7(1):1–46, 2010. [33] Y. Xu and C.W. Shu. Local discontinuous Galerkin methods for two classes of two-dimensional nonlinear wave equations. Physica D: Nonlinear phenomena, 208(1-2):21–58, 2005. [34] Y. Xu and C.W. Shu. Error estimates of the semi-discrete local discontinuous Galerkin method for nonlinear convection-diffusion and KdV equations. Comput. Methods Appl. Mech. Engrg., 196(37-40):3805–3822, 2007. [35] J. Yan and C.W. Shu. A local discontinuous Galerkin method for KdV type equations. SIAM J. Numer. Anal., 40:769–791, 2003. [36] N. Yi, Y. Huang, and H. Liu. A direct discontinuous Galerkin method for the generalized Korteweg-de Vries equation: Energy conservation and boundary effect. J. Comput. Phys., 242:351–366, 2013. Iowa State University, Mathematics Department, Ames, IA 50011 E-mail address: [email protected] Iowa State University, Mathematics Department, Ames, IA 50011 E-mail address: [email protected]
1. Introduction In this paper, we are interested in numerical approximations to the Burgers-Poisson (BP) equation of the form � 2 � u (1a) ut + − φ = 0, 2 x (1b) φxx − φ = u.
The subscript t (or x, respectively) denotes the differentiation with respect to time variable t (or x). System (1) can be rewritten as a nonlocal equation (2)
ut + uux + [G ∗ u]x = 0
with the kernel G(x) = 12 e−|x| . This nonlocal model was found as a simplified shallow water model by Whitham [30] to approximate the model with a singular kernel � � � 1 tanh k 1/2 ikx G(x) = e dk. 2π R k
As for the BP equation, analysis of traveling waves in [14] showed that the equation features wave breaking in finite time, and the authors also proved the global existence of weak entropy solutions. In this work, we develop a local discontinuous Galerkin method (LDG) to solve this nonlinear BP equation. Our proposed scheme is high order accurate, and preserves two invariants of momentum and energy, hence producing solutions with satisfying long time behavior. The L2 stability of the scheme for general solutions is a consequence of the energy preserving property. In the context of water waves, one of the best known local models is probably the Korteweg de Vries (KdV)-equation, ut + uux + uxxx = 0. This equation possess soliton solutions’ coherent structures that interact nonlinearly among themselves, then reemerge, retaining their identity and showing particle-like scattering behavior. Date: September 07, 2013. 2000 Mathematics Subject Classification. 65M60, 65M12, 35Q53. Key words and phrases. Discontinuous Galerkin method, Burgers-Poisson equation, conservation, stability. 1
2
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
In shallow water wave theory, the nonlinear shallow water equations which neglect dispersion altogether lead to the finite time wave breaking. On the other hand the third order derivative term in the KdV equation will prevent this ever happening in its solutions. In reality, some water waves appear to break, if the wave height is above certain threshold. Therefore in [30] Whitham raised an intriguing question: what kind of mathematical equation can describe waves with breaking? He suggested equation (2) with the above two kernels; many competing models have since been suggested to capture one aspect or another of the classical water-wave problem, see e.g. [15, 2, 3, 19, 13, 18, 26, 12, 25]. One common feature of these models is the associated global invariants, infinitely many or finitely many. The BP equation has the following two invariants: � � (3a) u(0, x) dx = u(t, x) dx =: E1 (t), � � 2 (3b) u (0, x) dx = u2 (t, x) dx =: E2 (t). It is desirable to design stable and high order accurate numerical schemes which preserve two invariants for solving the BP equation. It is believed that numerical methods preserving more invariants are advantageous: besides the high accuracy of numerical solutions, an invariant preserving scheme can preserve good stability properties after long-time numerical integration. Much more effort has been devoted in this topic for different integrable PDEs recently [16, 28, 27, 5]. The goal of this paper is to develop a discontinuous Galerkin (DG) method to preserve both momentum and energy at the discrete setting. The DG method is a class of finite element methods using completely discontinuous piecewise-polynomials for the numerical solution and the test functions. It was first designed and has been successful for solving first order PDEs such as nonlinear conservation laws [29, 8, 9, 7, 11]. The local discontinuous Galerkin (LDG) method is an extension of the DG method for solving higher order PDEs. It was first designed for convection-diffusion equations [10], and has been extended to other higher order wave equations, including the KdV equation [35, 33, 22, 34] and the Camassa-Holm equation [31], see also the recent review paper [32] on the LDG methods for higher order PDEs. The idea of the LDG method is to rewrite higher order equations into a first order system, and then apply the DG method on the system. In contrast, the direct discontinuous Galerkin methods, proposed in [23, 24] primarily for diffusion equations, aimed at directly solving higher order PDEs by the DG discretization, see e.g., [1, 36] for energy preserving DG methods for KdV type equations, and [21] for the Degasperis-Procesi equation. In this work we propose an LDG method based on formulation (1), for which the second equation was rewritten into a first order system for applying the LDG discretization. In the algorithm we update the solution in two steps: (1) given u, obtain φ by solving (1b) with the LDG method; (2) with the obtained φ, update u by solving (1a) with a standard DG method using a conservative numerical flux so that the resulting scheme preserves two integrals E1 and E2 in smooth region. As for error estimates, we define a global projection dictated by the selected numerical flux and obtain the needed projection error, following the strategy of error estimates carried out in [20] for the DDG method to solve convection-diffusion equations. Through careful estimates using this global projection, we obtain the optimal order of accuracy for polynomial elements of even degree. This is confirmed by the numerical tests with k = 2, 4. Numerical tests also show that for k odd, only k-th order of accuracy is observed. Such an optimal error estimate only for k =even was also shown in [1], and numerically observed in [1, 36] for KdV type equations. The
LDG FOR BP EQUATION
3
main feature of the scheme presented in this work is its capability to produce wave solutions with satisfying long time behavior. The rest of the paper is organized as follows. In section 2, we formulate our LDG method with a class of numerical fluxes. In section 3, we show that the LDG method for solving the Poisson equation (1b) is well defined and stable, and the method is shown to conserve both momentum and energy for the given numerical fluxes. In section 4, we obtain the optimal order of error between the numerical solution and smooth solutions for the conservative scheme when using polynomial elements of even degree. Finally, in section 5, we present numerical examples to illustrate the capacity of the LDG scheme to preserve two invariants after a long-time simulation. 2. The discontinuous Galerkin method 2.1. LDG formulation. We develop a local discontinuous Galerkin (LDG) method for the BP equation subject to initial data u0 (x) and periodic boundary conditions. Let us partition the interval I = [0, L] into 0 = x1/2 , x3/2 , . . . , xN +1/2 = L to get N equal subintervals and denote � � each cell by Ij = [xj−1/2 , xj+1/2 ], j = 1, . . . , N . The center of the cell is xj = 12 xj−1/2 + xj+1/2 . The piecewise polynomial space Vhk is defined as the space of polynomials of degree up to k in each cell Ij , that is, Vhk = {v : v|Ij ∈ P k (Ij ), j = 1, 2, . . . , N }.
(4)
Note that functions in Vhk are allowed to have discontinuities across the interfaces. The solution of the DG method is denoted by uh , which belongs to the finite element space Vhk . We denote − the limit values of uh at xj+1/2 from the right and from the left by (uh )+ j+1/2 and (uh )j+1/2 respectively. Let ω be a piecewise smooth function, its jump across the cell interface be denoted + − by [ω] := ω + − ω − , and its average at the cell interface, ω +ω , be denote by {ω}. 2 To define the LDG method, we introduce an auxiliary variable p = φx and rewrite (1a)-(1b) as follows: � 2� u (5a) ut + − p = 0, 2 x (5b) p − φx = 0, px − φ − u = 0.
(5c)
Then, the scheme is defined as follows: find uh , ph , φh ∈ Vhk such that �
(6a) (6b)
�
� �2 � u u2h h � (uh )t ρ dx − ρx dx + ρ ∂Ij − ph ρ dx = 0, 2 Ij Ij 2 Ij � � � �h γ �∂I = 0, ph γ dx + φh γx dx − φ j −
(6c) (6d)
�
Ij
Ij
Ij
� ph qx dx + p�h q �∂Ij − �
�
Ij
(φh + uh )q dx = 0,
Ij
(uh − u)|t=0 v dx = 0,
for all test functions ρ, γ, q, v in the finite element space Vhk . The choice for numerical fluxes �2 , φ � , p� is given by u h
h
h
4
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
�2 = 1 �(u+ )2 + u+ u− + (u− )2 � , u h h h h h 3 + − �h = θφ + (1 − θ)φ , φ h h
(7a) (7b)
− p�h = (1 − θ)p+ h + θph ,
(7c)
where θ ∈ [1/2, 1]. Here, the numerical fluxes at the endpoints of I can be defined using − − + + 2 U1/2 := UN +1/2 and UN +1/2 := U1/2 where U represents uh , φh , or ph . The resulting LDG scheme (6) with the fluxes (7) is called LDG-C. �2 is needed in order to capture the entropy For discontinuous solutions, an entropy flux for u h solution. One well-known choice is the Lax-Friedrich flux of the form �2 = 1 �(u− )2 + (u+ )2 − σ(u+ − u− )� , σ = 2 max |u|, (8) u h h h h h + 2 u∈[u− h ,uh ] with which the resulting LDG scheme is called LDG-D. In practice, one may adopt an adaptive numerical flux � � 1 (u+ )2 + u+ u− + (u− )2 if |u+ − u− | < 10−2 h h h h � 2 3 (9) uh = 1 �(u− )2 + (u+ )2 − 2σ(u+ − u− )� otherwise. h h h h 2
Here, 10−2 may vary as long as it can serve as a shock detector. Notation. We use � · �m,Ω as the H m -norm over domain Ω, and | · |m,Ω as its semi-norm. For m = 0, we simply use � · �Ω to denote the L2 -norm over domain Ω. We also use the notation � · �∞,Ω to denote the L∞ norm over domain Ω. The domain Ω could be a computational cell Ij or a master domain Iˆ := [−1, 1]. If Ω is the whole domain, we do not specify the domain unless necessary. 3. Analytical properties of the scheme 3.1. Existence, uniqueness, and stability. In this section, we prove the existence, uniqueness, and stability of ph , φh obtained from (6b)-(6c) with numerical fluxes (7b)-(7c), given uh . Lemma 3.1. The numerical scheme (6b)-(6c) with the numerical flux (7b)-(7c) satisfies �ph �2 + �φh �2 ≤ �uh �2 .
(10)
Proof. We choose γ = ph and q = φh . Then (6b)-(6c) gives � � � 2 �h ph �∂I = 0, ph dx + φh (ph )x dx − φ j Ij
−
�
Ij
Ij
� ph (φh )x dx + p�h φ �∂Ij −
Subtracting the two equations above gives � � 2 � ph + φ2h dx Ij
=−
�
Ij
� �h ph �∂I − φh (ph )x dx + φ j
�
Ij
�
Ij
φ2h dx
=
�
uh φh dx. Ij
� ph (φh )x dx + p�h φ �∂Ij −
�
uh φh dx. Ij
LDG FOR BP EQUATION
5
Take summation over j and use the periodic boundary condition to get � � � � �� � 2 � 2 � ph + φh dx = − φh (ph )x dx − φh [ph ] − ph (φh )x dx 1 I
I
−
� j
j+ 2
j
(p�h [φ])j+ 1 − 2
�
uh φh dx.
I
The first four terms on the right-hand side can be simplified to � � �� � − (φh ph )x dx − φh [ph ] + p�h [φ] I
��
=
j
I
j
�h [ph ] − p�h [φ] [φh ph ] − φ
�
j+ 12
j+ 12
= 0,
because of the choice of numerical fluxes (7b)-(7c). Therefore, we have that � � � 2 � 1 1 2 ph + φh dx ≤ |uh φh | dx ≤ �uh �2 + �φh �2 , 2 2 I I
�
which proves (10).
Remark: The inequality (10) implies the stabiliy of ph and φh since their L2 norms are bounded by �uh �. Moreover, (6b)-(6c) can be written in terms of a finite dimensional invertible linear mapping. Therefore, the inequality (10) shows that if uh = 0, then so are ph and φh . This proves the uniqueness and existence of ph and φh for any given uh . 3.2. Discrete conservation laws. In this section, we look at the properties of the numerical solution uh that are analogous to (3a)-(3b). Theorem 3.2. The following relations hold for all t > 0: � L � L (11) uh (t, x)dx = uh (0, x)dx, (12)
�
0
L
0
u2h (t, x)dx
=
�
0
L
0
u2h (0, x)dx
+
� j
(2θ − 1)
�
0
t�
[φh ]2 + [ph ]2
�
j+ 12
dτ .
Hence the scheme is conservative for θ = 1/2, and the scheme is energy stable for 1/2 < θ ≤ 1. Proof. Because (6) holds for any test function in Vhk , we choose ρ = 1 and γ = 1 in (6a)-(6b), respectively, to obtain �
(uh )t dx + Ij
�2 � � u h � �� ∂Ij − φh ∂Ij = 0. 2
Take summation over all j and use the periodic boundary condition, we have d dt This proves (11).
�
uh dx = I
�2 � �u � h � �� ∂Ij + φh ∂Ij = 0, 2 j
6
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
Next, we choose the test functions ρ = uh , γ = −φh , and q = −ph in (6a)-(6c) to obtain � � � �2 � u u2h h � (13) (uh )t uh dx − (uh )x dx + uh ∂Ij − ph uh dx = 0, 2 Ij Ij 2 Ij � � � �h φh �∂I = 0, (14) − ph φh dx − φh (φh )x dx + φ j �
(15)
Ij
Ij
Ij
� ph (ph )x dx − p�h ph �∂Ij +
�
(φh + uh )ph dx = 0,
Ij
Integrating some terms out and adding the above three relations together, we get � � � � � � � �2 u u3h �� φ2h �� p2h �� h � (uh )t uh dx + uh − ∂Ij + φh φh − ∂Ij − p�h ph − ∂Ij = 0. 2 6 2 2 Ij
Summing the terms above for all j = 1, . . . , N and using the periodic boundary condition, we get � � � � � �2 2�� 2�� � u u3h �� φ p 1d 2 h h � �∂I + p�h ph − h �∂I uh dx = − uh − ∂Ij − φh φh − j j 2 dt 2 6 2 2 j � � � � �2 � u 1 3 1 2 h �h [φh ] − [φ ] = [uh ] − [uh ] + φ 2 6 2 h j+ 1 1 j 2 j+ 2 � � 1 − p�h [ph ] − [p2h ] 2 j+ 12 � �� �h − {φh } [φh ] + ({ph } − p�h ) [ph ] = φ j
=
�� j
1 θ− 2
�
which proves (12).
�
[φh ]2 + [ph ]2
�
� 4. Error estimations
In this section we estimate the error from approximating ph , φh , uh . We proceed by defining a global projection with some established properties to prove that the errors from approximating ph and φh can be controlled by the errors from approximating uh . Then, we show that the error from approximating uh is of optimal order. 4.1. The global projection. Let ω be a piecewise smooth function on I = [0, L], we define the projection Qθ such that it satisfies the following properties: � � (16a) (Qθ ω) v dx = ωv dx, ∀v ∈ P k−1 , j = 1, · · · , N, Ij
(16b) where
Ij
� Q �j+1/2 , θ ω j+ 1 = ω 2
j = 1, . . . , N,
V� := θV + + (1 − θ)V − .
LDG FOR BP EQUATION
7
For j = N , we use the periodic extension to define (Qθ ω)+ N +1/2 , in order to be consistent with the numerical flux defined in (7). We first show that the projection Qθ w is well defined. Lemma 4.1. The projection Qθ satisfying (16) is uniquely defined for either θ �= with k even and N odd.
1 2
or θ =
1 2
Proof. Let {ψl }kl=0 be a set of orthogonal Legendre polynomials on [−1, 1] of degree up to k. We can write the projection Qθ of ω ∈ H k+1 (I) on each cell Ij as k
(Qθ ω)(xj +
� j h ξ) = al ψl (ξ), 2
ξ ∈ [−1, 1].
l=0
With v = ψi , the condition (16a) gives � 2i + 1 1 h j (17) ai = ω(xj + ξ)ψi (ξ) dξ, i = 0, . . . , k − 1, 2 2 −1 �1 2 2 where we have used −1 ψi (ξ) dξ = 2i+1 . j It remains to determine ak for j = 1, . . . , N . Since ψl (±1) = (±1)l , the condition (16b) gives � k � � k � � j+1 � j (18) θ al (−1)l + (1 − θ) al = ω ˆ (xj+ 1 ), j = 1, . . . , N. l=0
2
l=0
Because ω is periodic, we require that k �
+1 aN ψl (ξ) = l
l=0
k �
a1l ψl (ξ),
l=0
∀ξ ∈ [−1, 1].
which allows us to write the system (18) as 1 1 − θ (−1)k θ 0 ··· 0 ak b1 2 0 1 − θ (−1)k θ · · · 0 a k b2 (19) .. .. · .. = .. , .. .. .. . . . . . . . k bN aN (−1) θ 0 0 ··· 1 − θ k �� � �� � k−1 j+1 k−1 j where bj = ω ˆ (xj+ 1 )−θ (−1)l −(1−θ) l=0 al l=0 al . The determinant of the coefficient 2
matrix A above is given by (1 − θ)N + (−1)N +1+kN θN , which is non-zero for all θ �= 12 . When θ = 12 , the determinant is non-zero whenever N is odd and k is even. This shows that the projection Qθ is unique whenever θ �= 12 . If θ = 12 , Qθ is unique when k is even and N is odd. � Lemma 4.2. We have the following projection error (20) where C∗ depends on k ≥ 1 and θ.
�Qθ ω − ω� ≤ C∗ hk+1 ,
Proof. The proof is carried out in two steps. Step 1. We first establish the following inequality (21)
�(Qθ − I)ω�2 ≤ Ch
N � j=1
�˜ ω j �2k+1,Iˆ,
8
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
where h ˆ ξ), ξ ∈ [−1, 1] = I. 2 By the Cauchy-Schwarz inequality, we have from (17) that for j = 1, . . . , N , ω ˜ j (ξ) := ω(xj +
(22)
2i + 1 j 2 �w ˜ �0,Iˆ, 2
|aji |2 ≤
(23)
i = 0, . . . , k − 1.
Hence N � k−1 � j=0 i=0
From (19) of the form
|aji |2 ≤ (k − 1/2)
N � j=1
�˜ ω j �20,Iˆ.
ak = A−1 b, T T where ak = [a1k , . . . , aN k ] and b = [b1 , . . . , bN ] , it follows that N � j=1
(ajk )2 = bT (A−1 )T A−1 b ≤ C ≤C ≤C
N � j=1
N � j=1
�
(˜ ω j (1))2 +
N �
(bj )2
j=1
k−1 �
(aj+1 )2 + l
l=0
k−1 �
(ajl )2
l=0
�
�˜ ω j �21,Iˆ,
where we have used the Sobolev inequality |˜ ω j |∞,Iˆ ≤ C�˜ ω j �1,Iˆ. Hence, 2
�Qθ ω� =
N � j=1
�Qθ ω�20,Ij
� k−1 � N � h� j 2 h j = (al ) �ψl �20,Iˆ + (ak )2 �ψk �20,Iˆ 2 2 j=1 l=0 � � N k−1 N � � � j 2 j 2 ≤h |ai | + (ak ) ≤ Ch �˜ ω j �21,Iˆ. j=1
i=0
j=1
Step 2. For any v ∈ Vhk (I), we have that Qθ v = v. Therefore, �Qθ ω − ω�2 = inf �(Qθ − I)(ω − v)�2 v∈Vhk
≤ Ch ≤ Ch which proves (20).
N �
inf
�˜ ω j − vˆ�2k+1,Iˆ
j=1
vˆ∈P k [−1,1]
N �
|˜ ω j |2k+1,Iˆ = Ch2k+2 |ω|2k+1 ,
j=1
�
LDG FOR BP EQUATION
9
Lemma 4.3. For k ≥ 1 the following inequality holds, �2 N � � � � − 2k+1 � � (24) . �(ω − Qθ ω)(xj+ 1 )� ≤ Ch 2
j=1
The constant C depends on k and θ.
Proof. On each interval Ij , using the orthogonality relation (16a), we have j
ω(x)|Ij := ω ˜ (ξ) = j
� Qθ ω(x)|Ij := (Q θ ω) (ξ) =
∞ � l=0 k=1 �
ωlj ψl (ξ), ωlj ψl (ξ) + αkj ψk (ξ).
l=0
Hence, by ψl (1) = 1, we have (25)
�∞ � � �2 � � �2 � � � � j − �(ω − Qθ ω)(x 1 )� ≤ 2 � ωl � + 2|αkj |2 . � j+ 2 � � � l=k
To control the first term on the right-hand side of (25), we consider the following expression j
(26)
∂ξ ω ˜ (ξ) =
∞ �
βlj ψl (ξ).
l=0
Following the idea in [4], we integrate (26) with respect to ξ to get � ξ ∞ � j j j ω ˜ (ξ) = ω ˜ (−1) + βl ψl (ν)dν. −1
l=0
Using the property of Legendre polynomials � ξ 1 ψi (ν)dν = (ψi+1 (ξ) − ψi−1 (ξ)) , 2i + 1 −1 we can write ω ˜ j (ξ) = ω ˜ j (−1) +
�
β0j
β1j − 3
�
ψ0 (ξ) +
∞ � l=1
�
j j βl−1 βl+1 − 2l − 1 2l + 3
Therefore, ωij =
�
j j βi−1 βi+1 − 2i − 1 2i + 3
�
,
i ≥ 1.
Thus, � j j β β k−1 k ωlj = + , 2k − 1 2k + 1 l=k � j+1 � ∞ � βk−1 βkj+1 j+1 l k ωl (−1) = (−1) − . 2k − 1 2k + 1 ∞ �
l=k
�
�
ψl (ξ).
10
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
These ensure the following estimates � �2 � � j j ∞ �� � 2(βk−1 )2 2(βk+1 )2 1 1 � j� (27a) ωl � ≤ + ≤ �∂ξ ω ˜ j �20,Iˆ � � � 2k − 1 2k − 1 2k + 1 2k − 1 l=k � �2 ∞ �� � 1 � j+1 l� (27b) ωl (−1) � ≤ �∂ξ ω ˜ j+1 �20,Iˆ. � � � 2k − 1 l=k
The second term on the right-hand side of (25) is determined by (16b), i.e., �∞ � �∞ � � j+1 � j j+1 j k l (28) θαk (−1) + (1 − θ)αk = θ ωl (−1) + (1 − θ) ωl , l=k
l=k
where we have used ω ˆ (xj+1/2 ) = θω ˜ j+1 (−1) + (1 − θ)˜ ω j (1). We then have from (27) and (28) that � �2 � �2 N N ∞ ∞ �� � �� � � � � � � �� |αkj |2 ≤ C ωlj+1 (−1)l � + � ωlj � � � � � j=1
j=1
≤C ≤C
N � j=1
N � j=1
l=k
l=k
�∂ξ ω ˜ j+1 �20,Iˆ + �∂ξ ω ˜ j �20,Iˆ �∂ξ ω ˜ j �20,Iˆ.
Insertion of these estimates back into (25) yields �2 N � N � � � � �(ω − Qθ ω)(x− 1 )� ≤ C �∂ξ ω ˜ j �20,Iˆ. � j+ 2 � j=1
j=1
Recall that Qθ v = v for any v ∈ P k , we proceed
�2 N � N � � � � − �(ω − Qθ ω)(x 1 )� ≤ C inf �∂ξ ω ˜ j − ∂ξ v˜�20,Iˆ � j+ 2 � v˜∈P k j=1
j=1
=C
≤C
N �
inf
�∂ξ ω ˜ j − p˜�2k,Iˆ
j=1
p˜∈P k−1
N �
|∂ξ ω ˜ j |2k,Iˆ
j=1
� �2k+2 � �−1 h h = C|ω|2k+1 ≤ C|ω|2k+1 h2k+1 . 2 2 � We will use the error estimates obtained in Lemmas 4.2-4.3 to estimate the error of the computed solution. We also need to use the following well-known inverse properties ( see e.g.,
LDG FOR BP EQUATION
11
[6]). �∂x w� ≤ Ch−1 �w�,
(29a)
�w�Γh ≤ Ch−1/2 �w�,
(29b)
�w�∞ ≤ Ch−1/2 �w�,
(29c)
for any w ∈ Vhk . The constant C is independent of w and h. The subscript Γh denotes the set of the boundary points of all the cells Ij . 4.2. Approximating p and φ. Lemma 4.4. Let (u, p, φ) be the exact solution of the system (5). Let (uh , ph , φh ) be obtained from (6) with the choice of fluxes (7). Let θ ∈ [0, 1] be such that both Qθ and Q1−θ are uniquely defined, then the following inequality holds for all t > 0. � � (30) �Q1−θ p − ph �2 + �Qθ φ − φh �2 ≤ 2 �Q1−θ p − p�2 + �Qθ φ − φ�2 + �u − uh �2 . Proof. Since the scheme with fluxes (7) is consistent, (6b)-(6c) also hold for (u, p, φ). In other words, � � � (31a) pγ + φγx − φγ �∂Ij = 0, Ij
−
(31b)
�
Ij
Ij
� pqx + pq �∂Ij −
�
(φ + u)q = 0,
Ij
Subtracting (6b)-(6c) from (31a)-(31b), we get the error equations � � � � � �h γ �∂I = 0, (32a) (p − ph )γ + (φ − φh )γx − φ − φ j Ij
−
(32b)
�
Ij
Ij
� (p − ph )qx + (p − p�h ) q �∂Ij −
�
Ij
(φ − φh )q =
�
Ij
(u − uh )q.
� � Define �p = Q1−θ p − p, wp = Q1−θ p − ph , ��p = Q �p = Q 1−θ p − p, and w 1−θ p − p�h . (Similar definition can be given for �φ , wφ , ��φ , and w �φ associated with Qθ .) We then choose γ = wp , q = wφ and take the summation of (32) over j to get �� �� � (wp − �p )wp + (wφ − �φ )(wp )x + (w �φ − ��φ ) [wp ]j+ 1 = 0, j
−
Ij
�� j
=
�
I
j
Ij
(wp − �p )(wφ )x −
(u − uh )wφ .
Ij
j
� j
(� wp − ��p ) [wφ ]j+ 1 − 2
2
�� j
Ij
(wφ − �φ )wφ
Take the difference of both equations, we get � � � � � 2 2 wp + wφ = �p wp + �φ wφ + �1 + �2 + �3 − (u − uh )wφ , I
I
I
I
I
12
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
where �1 = − �2 =
j
�� j
�3 =
��
�� j
Ij
(wφ (wp )x + wp (wφ )x ) −
�φ (wp )x + Ij
� j
�p (wφ )x + Ij
� j
� j
��φ [wp ]j+ 1 ,
w �φ [wp ]j+ 1 − 2
� j
w �p [wφ ]j+ 1 , 2
2
��p [wφ ]j+ 1 . 2
First, note that wφ (wp )x + wp (wφ )x = (wp wφ )x , so � � � �1 = [wp wφ ]j+ 1 − w �φ [wp ]j+ 1 − w �p [wφ ]j+ 1 = 0, 2
j
2
j
2
j
with the choice of numerical fluxes (7b)-(7c). As for �2 , the property (16a) of Qθ gives �� �� �φ (wp )x = (Qθ φ − φ)(wp )x = 0, Ij
j
j
Ij
since (wp )x is in P k−1 . On the other hand, the propery (16b) of Qθ gives � � �� � ��φ [wp ]j+ 1 = Q θ φ − φ [wp ]j+ 1 = 0. j
2
j
2
Similarly, the term �3 vanishes by the properties (16) of Q1−θ . Using the Young’s inequality together with the results from �1 , �2 , and �3 , we get (33)
1 1 1 �wp �2 + �wφ �2 ≤ ��p �2 + ��φ �2 + �u − uh �2 , 2 2 2
which proves (30).
�
4.2.1. Approximating u. Theorem 4.5. If k is even, then the numerical solution, uh , obtained from the scheme (6) and the numerical fluxes (7) with θ = 12 satisfies (34)
�u − uh � ≤ Chk+1 .
The constant C is independent of the mesh size. Proof. Since the scheme (6) with fluxes (7) is consistent, (6a) also holds for (u, p, φ). In other words, � � � �2 � u2 u � ut ρ − ρx + ρ ∂Ij − (35) pρ = 0. 2 Ij Ij 2 Ij
Define w = Q1/2 u − uh and � = Q1/2 u − u. We have that u − uh = w − �. Subtracting (6a) from (35) and choose ρ = w, we get � � � � � � � 2 � �2 � u2h u u2 u h � wt w = �t w + − wx − − w ∂Ij + ep w, 2 2 2 2 Ij Ij Ij Ij where ep = p − ph .
LDG FOR BP EQUATION
13
Take summation over all j and introduce {uh }2 /2 into the third term on the right-hand side to get � � � � � u2 u2 � � � u2 {uh }2 � h wt w = �t w + − wx + − [w]j+ 1 2 2 2 2 2 I I Ij j j � � � �2 � {uh }2 u h + − [w]j+ 1 + ep w. 2 2 2 I j
Using the identity A2 /2 − B 2 /2 = A(A − B) − (A − B)2 /2, we get � � � �� 1� wt w = �t w + u (u − uh ) wx − (u − uh )2 wx 2 I I Ij Ij j
+
� j
+
� j
j
u (u − {uh }) [w]j+ 1 − 2
�
�2 {uh }2 u − h 2 2
�
1� 2
[w]j+ 1 + 2
j
�
(u − {uh })2 [w]j+ 1 2
ep w. I
Let {w} = {Q1/2 u} − {uh } and {�} = {Q1/2 u} − u, we can write � � � (36) wt w = �t w + τ1 + τ2 + τ3 + τ4 + τ5 + ep w, I
I
I
where
τ1 =
�� j
τ2 = −
uwwx + Ij
�� j
�
u{w}[w]j+ 1 ,
j
Ij
��
u�wx −
� j
2
u{�}[w]j+ 1 , 2
1� {w}2 [w]j+ 1 , 2 2 I j j j � �� � 1� 1� τ4 = w�wx − �2 wx + {w}{�}[w]j+ 1 − {�}2 [w]j+ 1 , 2 2 2 2 Ij Ij j j j j � � �2 � {uh }2 u h τ5 = − [w]j+ 1 . 2 2 2 τ3 = −
1 2
w 2 wx −
j
Note that
� � w2 τ1 = u + u{w}[w]j+ 1 2 2 x Ij j j � 2� � � � w2 � �� w =− u − ux + u{w}[w]j+ 1 2 2 j+ 1 2 Ij j j j 2 � � � � w2 1 =− ux ≤ �ux �∞ �w�2 . 2 2 Ij ��
j
�
14
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
As for τ2 , we write u(x) = u(xj ) + u� (x∗ )(x − xj ) for all x ∈ Ij where x∗ is between x and xj . Therefore, � � � � � �� � � � � � � � � � � � � � ∗ �− � � u�wx − u{�}[w]j+ 1 � ≤ �− u(xj ) �wx − u (x ) �(x − xj )wx �� + |0| � 2 Ij Ij � j Ij � � j � j j � � � � � � �� � �� � ∗ � � � � � ≤� u(xj ) �wx � + u (x ) |�hwx | Ij Ij � j � j � � 1 ≤ |0| + �ux �∞ ���2 + h2 �wx �2 2 ≤ C(h2k+2 + �w�2 ), because of the projection properties (16), the inverse property (29a), and Lemma 4.2. For τ3 , we can show that � � � � 3� 1� w 1� 1 � w3 1� 2 − − {w} [w]j+ 1 = − {w}2 [w]j+ 1 2 2 2 3 x 2 2 3 2 Ij j
j
j
j
1 � 3 = [w]j+ 1 ≤ Ch−1 �w�∞ �w�2 ≤ Ch−3/2 �w�3 , 2 24 j
by the inverse properties (29b)-(29c). From the inverse properties (29a) and (29c), it follows that �wx �∞ ≤ Ch−3/2 �w�, with which we are able to estimate terms in τ4 by using Lemma 4.2, �� w�wx ≤ C�w�����wx �∞ ≤ Chk+1−3/2 �w�2 = Chk−1/2 �w�2 , j
Ij
� 1� 2
j
1 �2 wx ≤ �wx �∞ ���2 ≤ Ch2k+1/2 �w�. 2 Ij
As for the remaining terms in τ4 , � 1� {w}{�}[w]j+ 1 − {�}2 [w]j+ 1 = 0 2 2 2 j
j
because of the projection property (16b). �2 /2 = −[u ]2 /24, and [u ] = [u − u ] = [w] − [�], we Finally, using the fact that {uh }2 /2 − u h h h h have � � �2 � {uh }2 u � 1 h τ5 = − [w]j+ 1 = − [uh ]2 [w]j+ 1 , 2 2 2 2 24 j
j
=
� j
−
1 1 1 [w]3j+ 1 + [�][w]2j+ 1 − [�]2 [w]j+ 1 2 2 2 24 12 24
≤ C�w�∞ (�w�2Γh + ���Γh �w�Γh + ���2Γh )
≤ Ch−1/2 �w�(h−1 �w�2 + hk+1/2 h−1/2 �w� + h2k+1 ) ≤ C(h2k+2 + �w�2 + h−3/2 �w�3 ),
LDG FOR BP EQUATION
15
where we have used the inverse properties (29b) and (29c). The results from τ1 to τ5 and (36) give � � d (37) �w�2 ≤ C h2k+2 + �w�2 + h−3/2 �w�3 . dt With this, we can write �w(t, ·)�2 ≤ CG(t),
(38) where (39)
G(t) = h h2k+2
and
F (η) =
�
This gives G(0) =
(41)
η 1
+
�
t 0
�w(τ, ·)�2 + h−3/2 �w(τ, ·)�3 dτ .
� � G� (t) ≤ C∗ G(t) + h−3/2 G(t)3/2 .
(40) Now, define
2k+2
1 � dξ = −3/2 ξ+h G(0)ξ 3/2
�
η 1
1 ξ+
hk−1/2 ξ 3/2
dξ.
For k ≥ 1, we have that F � (η) is uniformly (with respect to h) positive and bounded above by 1 for all η > 1. Thus, there exists C˜ = F −1 (C∗ T ) for given T > 0. Integrate (40) to get � � G ˜ (42) F ≤ C∗ T = F (C). G(0) G ˜ Using this and (38), we prove (34) as desired. Hence ≤ C. G(0) � 4.3. Time discretization. We partition the time interval [0, T ] into M equal subintervals with boundaries {tn }, n = 0, 1, 2, . . . , M . Set ∆t = T /M as the time step. In order to preserve both mass and energy at the fully discrete level, we may use the Crank-Nicolson time discretization to find un+1 = 2u∗h − unh , h where u∗h is determined by �� �2 � � � � u∗h u∗h − unh (u∗h )2 2 ρ− ρx + ρ �∂Ij − p∗h ρ = 0, ∆t 2 2 I Ij Ij � j � � �∗ γ �∂I = 0, (43) p∗h γ + φ∗h γx − φ j h Ij
−
�
Ij
Ij
p∗h qx
� + p�∗h q �∂Ij −
�
Ij
(φ∗h + u∗h ) q = 0.
Indeed, this time discretization has the desired and provable properties. Theorem 4.6. The fully-discrete scheme (43) gives solution unh that satisfies � L � L n (44) uh dx = u0h dx, (45) for all 0 ≤ n ≤ M .
�
0
L
0
(unh )2 dx =
�
0
0
L�
u0h
�2
dx,
16
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
Proof. Take the test functions ρ = 1, γ = 1 in (43), adding together, to get �� �2 � � u∗h �� un+1 − unh h �∗ � + ∂Ij − φh ∂Ij = 0, ∆t 2 Ij
which upon summation over j proves (44). Next, we choose the test functions ρ = u∗h , γ = −φ∗h , and q = −p∗h so that �� �2 � � � u∗h un+1 − unh ∗ (u∗h )2 ∗ h uh − (uh )x + u∗h �∂Ij ∆t 2 2 Ij Ij � − p∗h u∗h = 0, �
− Ij
�
Ij
p∗h φ∗h −
p∗h (p∗h )x
�
Ij
Ij
� �∗ φ∗ �∂I = 0, φ∗h (φ∗h )x + φ j h h
� − p�∗h p∗h �∂Ij +
�
Ij
(φ∗h + u∗h )p∗h = 0.
Summation of the above three equations over j gives � � � � � � un+1 2 − (un )2 � L un+1 2 − (un )2 h h h h = = 0, 2∆t 2∆t Ij 0 j
which leads to (45).
�
Note that the above time discretization is fully nonlinear and requires the costly iteration solver. In practice, one would prefer to use some explicit solver with high order accuracy for time discretization. In our numerical simulation we choose to use the TVD third-order RungeKutta method [17] to solve the ODE system of the form a˙ = L(a): a(1) = an + ∆tL(an ), 3 1 (46) a(2) = an + a(1) + 4 4 1 2 an+1 = an + a(2) + 3 3 n n where a is the coefficient vector of uh .
1 ∆tL(a(1) ), 4 2 ∆tL(a(2) ), 3
5. Numerical Tests It is known [14] that one steady solution of system (1a)-(1b) is given by � 4 � −|x|/2 (47) U1 (x) = e −1 . 3 The system also has a steady periodic solution of the form � � � � 4 cosh x2 � � −1 , (48) U2 (x) = 3 cosh p2
for −p < x < p and by periodic continuation with period 2p. Because the system (1a)-(1b) is Galilean invariant, a family of traveling-wave solutions (1a)-(1b) may be obtained from the steady solutions as (49)
u(t, x) = U (x − u0 t) + u0 ,
LDG FOR BP EQUATION
17
where U is the steady state solution (47) or (48). We will use both steady and traveling wave solutions to test our scheme. Example 1. (Accuracy test) We run the semi-discrete scheme (6) and the numerical flux (7) with θ = 12 , 1, along with the third order Runge-Kutta method (46) on the steady state problem which has (48) as its exact solution. The results for k = 1, 2, 3, 4 are given in the two tables below. Here, we use ∆t = 0.001, final time tmax = 2, and p = 2. The norms of the error were computed by using the sixteen-point Guass quadrature rule. k N 1 10 20 40 80 2 10 20 40 80 3 10 20 40 80 4 10 20 40 80 k N 1
5 10 20 40 80 2 5 10 20 40 80 3 5 10 20 40 80 4 5 10 20 40 80
L1 5.5907e-02 2.8223e-02 1.4137e-02 7.0694e-03 1.6478e-03 2.0461e-04 2.5457e-05 3.1778e-06 2.5013e-04 3.1133e-05 3.8845e-06 4.8514e-07 1.2182e-06 3.7401e-08 1.1545e-09 3.5841e-11 L1 1.3616e-01 5.9425e-02 2.7403e-02 1.3131e-02 6.4241e-03 2.5044e-03 2.2204e-04 1.8941e-05 1.5556e-06 1.3670e-07 2.0147e-04 1.9490e-05 2.0788e-06 2.3547e-07 2.7898e-08 4.0456e-06 9.1182e-08 2.0673e-09 4.5009e-11 1.0746e-12
order 0.9862 0.9974 0.9998 3.0096 3.0067 3.0019 3.0062 3.0026 3.0013 5.0256 5.0178 5.0095 order 1.1962 1.1167 1.0614 1.0314 3.4956 3.5513 3.6059 3.5084 3.3697 3.2289 3.1421 3.0773 5.4715 5.4630 5.5214 5.3883
θ = 1/2 L2 order 3.3360e-02 1.6765e-02 0.9926 8.3909e-03 0.9986 4.1959e-03 0.9999 1.0558e-03 1.3282e-04 2.9908 1.6615e-05 2.9990 2.0768e-06 3.0000 1.7960e-04 2.2535e-05 2.9945 2.8200e-06 2.9984 3.5262e-07 2.9995 9.6189e-07 3.0251e-08 4.9908 9.4598e-10 4.9990 2.9535e-11 5.0013
L∞ 4.8122e-02 2.6199e-02 1.3877e-02 7.1781e-03 2.6115e-03 3.8008e-04 5.3250e-05 7.2587e-06 6.4039e-04 9.9254e-05 1.4970e-05 2.2034e-06 4.5036e-06 1.8170e-07 7.1522e-09 2.7574e-10
θ=1 L2 order 8.2323e-02 3.6667e-02 1.1668 1.7300e-02 1.0837 8.4184e-03 1.0391 4.1554e-03 1.0186 1.4860e-03 1.4181e-04 3.3894 1.3206e-05 3.4246 1.2308e-06 3.4235 1.1673e-07 3.3984 1.2239e-04 1.2763e-05 3.2614 1.4660e-06 3.1221 1.7711e-07 3.0492 2.1807e-08 3.0218 2.5910e-06 5.7074e-08 5.5045 1.3796e-09 5.3705 3.1612e-11 5.4476 7.7906e-13 5.3426
L∞ 1.1532e-01 5.2580e-02 2.3941e-02 1.1365e-02 5.6409e-03 2.4239e-03 3.1615e-04 3.8885e-05 4.7536e-06 5.8323e-07 1.8214e-04 2.8368e-05 4.6833e-06 7.2428e-07 1.0497e-07 7.9819e-06 1.3149e-07 4.7551e-09 1.7362e-10 5.7253e-12
order 0.8772 0.9168 0.9510 2.7805 2.8355 2.8750 2.6897 2.7290 2.7643 4.6314 4.6670 4.6970 order 1.1330 1.1350 1.0749 1.0106 2.9386 3.0233 3.0321 3.0269 2.6827 2.5987 2.6929 2.7866 5.9238 4.7893 4.7754 4.9225
18
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
The results show that the optimal order of accuracy is achieved only when k =even, which is consistent with our theoretical result on the optimal error estimates for k = even. Also such an observation seems unaffected by the choice of θ ∈ [0, 1], though we only display results for θ = 1/2 and θ = 1. Example 2. (Energy-preserving test) We compare the performance of the LDG-C and LDG-D schemes on the traveling wave version of (47), with u0 = 1/100. We use θ = 12 , 1 along with the third order Runge-Kutta method (46). The relative L2 energy �uh (t, x) − uh (0, x)� of the computed solution for k = 4 are shown in Figure 1. Here, we use ∆t = 0.001, tmax = 100. The L2 norms were computed by using the sixteen-point Gauss quadrature rule. In next two examples, we use polynomial elements of degree k = 2 along with the TVBM limiter introduced in [9]. Here, we use the mesh size h = 1/16. Example 3. We test the conservative scheme for initial data 0 ≤ x ≤ 1, 0.5 −x + 1.5, 1 ≤ x ≤ 4, u0 (x) = −2.5, x ≥ 4.
This initial data has a downward ramp of height 3 and the constant states lying symmetric with respect to u = −1, the solution is expected to converge to a stationary solution. In Figure 2, we observe a stable pattern formation as analyzed in [14]. In our experiment we use a modified initial data in C 2 , which agrees with the original data everywhere except for near x = 1, 4, so that we can apply directly the TVBM limiter introduced in [9]. Our goal is to observe the stable wave pattern, so the choice of modification is not essential. Example 4. We consider another initial data of the form x ≤ 8, −0.5 15.5 − 2x, 8 ≤ x ≤ 8.5, u0 (x) = −1.5, x ≥ 8.5. This example with smaller jump has no stable stationary solution to converge. We plot the computed solution at the different times in Figure 3, from which we can see that dispersive effects with oscillations propagate to the left of the ramp as analyzed in [14]. Acknowledgments This research was partially supported by the National Science Foundation under grant DMS1312636 and the KI-Net research network. References [1] J. L. Bona, H. Chen, O. Karakashian, and Y. Xing. Conservative, discontinuous Galerkin methods for the generalized Korteweg-de Vries equation. Math. Comp., 82(283):1401–1432, 2013. [2] R. Camassa and D. D. Holm. An integrable shallow water equation with peaked solitons. Phys. Rev. Lett., 71(11):1661–1664, 1993. [3] R. Camassa, D. D. Holm, and J. M. Hyman. A new integrable shallow water equation. Advances in Applied Mechanics, 31(31):1–33, 1994. [4] P. Castillo, B. Cockburn, D. Sch¨ otzau, and C. Schwab. Optimal a priori error estimates for the hp-version of the local discontinuous Galerkin method for convection-diffusion problems. Math. Comp., 71(238):455–478, 2002. [5] E. Celledoni, V. Grimm, R. I. McLachlan, D. I. McLaren, D. O’Neale, B. Owren, and G. R. W. Quispel. Preserving energy resp. dissipation in numerical PDEs using the “average vector field” method. J. Comput. Phys., 231(20):6770–6789, 2012.
LDG FOR BP EQUATION
19
Figure 1. Example 2: the evolution of the relative L2 energy over long-time simulation. [6] P. G. Ciarlet. The Finite Element Method for Elliptic Problems. North-Holland, New York, 1978. [7] 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. Math. Comp., 54(190):545–581, 1990. [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. J. Comput. Phys., 84(1):90–113, 1989. [9] B. Cockburn and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework. Math. Comp., 52(186):411–435, 1989. [10] B. Cockburn and C.-W. Shu. The local discontinuous Galerkin method for time-dependent convectiondiffusion systems. SIAM J. Numer. Anal., 35(6):2440–2463 (electronic), 1998. [11] B. Cockburn and C.-W. Shu. The Runge-Kutta discontinuous Galerkin method for conservation laws. V. Multidimensional systems. J. Comput. Phys., 141(2):199–224, 1998. [12] A. Constantin and D. Lannes. The hydrodynamical relevance of the Camassa-Holm and Degasperis-Procesi equations. Arch. Ration. Mech. Anal., 192(1):165–186, 2009. [13] A. Degasperis, D. D. Kholm, and A. N. I. Khon. A new integrable equation with peakon solutions. Teoret. Mat. Fiz., 133(2):170–183, 2002. [14] K. Fellner and C. Schmeiser. Burgers-Poisson: a nonlinear dispersive model equation. SIAM J. Appl. Math., 64(5):1509–1525 (electronic), 2004. [15] B. Fuchssteiner and A. S. Fokas. Symplectic structures, their B¨ acklund transformations and hereditary symmetries. Phys. D, 4(1):47–66, 1981/82. [16] D. Furihata and M. Mori. General derivation of finite difference schemes by means of a discrete variation. Transactions-Japan Society for Industrial and Applied Mathematics, 8:317–340, 1998. [17] S. Gottlieb and C.-W. Shu. Total variation diminishing Runge-Kutta schemes. Math. Comp., 67(221):73–85, 1998. [18] D. D. Holm and M. F. Staley. Wave structure and nonlinear balances in a family of evolutionary PDEs. SIAM J. Appl. Dyn. Syst., 2(3):323–380 (electronic), 2003.
20
HAILIANG LIU AND NATTAPOL PLOYMAKLAM
Figure 2. Example 3: the computed solution at t = 0, 10, 100. [19] R. S. Johnson and Camassa-Holm. Korteweg-de Vries and related models for water waves. J. Fluid Mech., 455:63–82, 2002. [20] H. Liu. Optimal error estimates of the direct discontinuous Galerkin method for convection-diffusion equations. Math. Comp., 2013. in review. [21] H. Liu, Y. Huang, and N. Yi. A direct discontinuous Galerkin method for the Degasperis–Procesi equation. Methods Appl. Anal., 2013. in review. [22] H. Liu and J. Yan. A local discontinuous Galerkin method for the Korteweg–de Vries equation with boundary effect. J. Comput. Phys., 215(1):197–218, 2006. [23] H. Liu and J. Yan. The direct discontinuous Galerkin (DDG) methods for diffusion problems. SIAM J. Numer. Anal., 47(1):675–698, 2008/09. [24] H. Liu and J. Yan. The direct discontinuous Galerkin (DDG) method for diffusion with interface corrections. Commun. Comput. Phys., 8(3):541–564, 2010. [25] H. Liu and Z. Yin. Global regularity, and wave breaking phenomena in a class of nonlocal dispersive equations. In Nonlinear partial differential equations and hyperbolic wave phenomena, volume 526 of Contemp. Math., pages 273–294. Amer. Math. Soc., Providence, RI, 2010. [26] Y. Liu and Z. Yin. Global existence and blow-up phenomena for the Degasperis-Procesi equation. Comm. Math. Phys., 267(3):801–820, 2006. [27] T. Matsuo. Dissipative/conservative Galerkin method using discrete partial derivatives for nonlinear evolution equations. J. Comput. Appl. Math., 218(2):506–521, 2008. [28] R. I. McLachlan, G. R. W. Quispel, and N. Robidoux. Geometric integration using discrete gradients. R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci., 357(1754):1021–1045, 1999. [29] W. H. Reed and T. R. Hill. Triangular mesh methods for the neutron transport equation. Los Alamos Report LA-UR-73-479, 1973. [30] G. B. Whitham. Linear and Nonlinear Waves. John Wiley & Sons, New York, 1974. [31] Y. Xu and C.-W. Shu. A local discontinuous Galerkin method for the Camassa-Holm equation. SIAM J. Numer. Anal., 46(4):1998–2021, 2008.
LDG FOR BP EQUATION
21
Figure 3. Example 4: the computed solution at t = 0, 5, 20. [32] Y. Xu and C.-W. Shu. Local discontinuous Galerkin methods for high-order time-dependent partial differential equations. Commun. Comput. Phys., 7(1):1–46, 2010. [33] Y. Xu and C.W. Shu. Local discontinuous Galerkin methods for two classes of two-dimensional nonlinear wave equations. Physica D: Nonlinear phenomena, 208(1-2):21–58, 2005. [34] Y. Xu and C.W. Shu. Error estimates of the semi-discrete local discontinuous Galerkin method for nonlinear convection-diffusion and KdV equations. Comput. Methods Appl. Mech. Engrg., 196(37-40):3805–3822, 2007. [35] J. Yan and C.W. Shu. A local discontinuous Galerkin method for KdV type equations. SIAM J. Numer. Anal., 40:769–791, 2003. [36] N. Yi, Y. Huang, and H. Liu. A direct discontinuous Galerkin method for the generalized Korteweg-de Vries equation: Energy conservation and boundary effect. J. Comput. Phys., 242:351–366, 2013. Iowa State University, Mathematics Department, Ames, IA 50011 E-mail address: [email protected] Iowa State University, Mathematics Department, Ames, IA 50011 E-mail address: [email protected]
