Discontinuous Galerkin method for Krause's ... - Brown University
Discontinuous Galerkin method for Krause’s consensus models and pressureless Euler equations1 Yang Yang2 , Dongming Wei3 and Chi-Wang Shu4 Abstract In this paper, we apply discontinuous Galerkin (DG) methods to solve two model equations: Krause’s consensus models and pressureless Euler equations. These two models are used to describe the collisions of particles, and the distributions can be identified as density functions. If the particles are placed at a single point, then the density function turns out to be a δ-function and is difficult to be well approximated numerically. In this paper, we use DG method to approximate such a singularity and demonstrate the good performance of the scheme. Since the density functions are always positive, we apply a positivity-preserving limiter to them. Moreover, for pressureless Euler equations, the velocity satisfies the maximum principle. We also construct special limiters to fulfill this requirement.
Keywords: δ-singularities, discontinuous Galerkin method, Krause’s model, pressureless gas, positivity preserving, symmetry preserving, maximum principle. 1
Research supported by DOE grant DE-FG02-08ER25863 and NSF grant DMS-1112700. Division of Applied Mathematics, Brown University, Providence, RI 02912. E-mail: yang [email protected] 3 Division of Applied Mathematics, Brown University, Providence, RI 02912. Current address: PNC Financial Services Group, Inc., New York, New York. E-mail:[email protected] 4 Division of Applied Mathematics, Brown University, Providence, RI 02912. E-mail: [email protected] 2
1
1
Introduction
In this paper, we apply discontinuous Galerkin (DG) methods to solve hyperbolic conservation law ut + f(u)x = 0, u(x, 0) = u0 (x),
(x, t) ∈ R × (0, T ], x ∈ R,
(1.1)
and its two dimensional version, where the exact solution u(x, t) contains δ-singularities. Such problems appear often in applications and are difficult to be well approximated numerically. Most of the previous techniques are designed to modify the singularities with smooth kernels in some narrow region, and hence may smear such singularities severely, leading to large errors in the approximation. On the other hand, the DG methods rely on weak formulations and can solve such problems directly, leading to very accurate results. The first work on approximating δ-singularities with DG method is [24], in which the authors proved the superconvergence results for linear hyperbolic equations with singular initial data and singular source terms. Moreover, several numerical examples were also given in [24] to demonstrate the advantages of the DG scheme in approximating δ-singularities for both linear and nonlinear hyperbolic equations. In this paper, we extend the work in [24] and consider applications to two model equations: Krause’s consensus models and pressureless Euler equations. The DG method, first introduced by Reed and Hill in 1973 [19], was extended to solve scalar linear hyperbolic equations by Johnson and Pitk¨aranta [18]. Later, Cockburn et al. studied Runge-Kutta discontinuous Galerkin (RKDG) methods for hyperbolic conservation laws [14, 12, 13, 15]. Recently, in [25], genuinely maximum-principle-satisfying high order DG schemes for scalar equations and two-dimensional incompressible flows in vorticitystreamfunction formulation have been constructed. Subsequently, positivity-preserving high order DG schemes for compressible Euler equations were given in [26]. In this paper, we extend the ideas in [25, 26] to construct bound-preserving high order DG schemes for the Krause’s consensus models and pressureless Euler equations.
2
For the first example, we discuss the following Krause’s consensus model equation ρt + Fx = 0, x ∈ [0, 1], t > 0, ρ(x, 0) = ρ0 (x), t > 0,
(1.2)
where ρ is the density function, which is always positive. The flux F is given by F (x, t) = v(x, t)ρ(x, t), and the velocity v is defined by v(x, t) =
Z
0
1
(y − x)ξ(y − x)ρ(y, t)dy,
where 0 ≤ ξ(x) ≤ 1 is supported on a ball centered at zero with radius R. In [7], Canuto et al. investigated the discretized version of the PDE and proved that when the time t tends to infinity, the density function ρ will converge to some isolated δ-singularities, and the distance between any two of these δ-singularities cannot be less than R. Some computational results are shown in [7] based on a first order finite volume method. For two dimensions, if the initial density is rotationally invariant, the limit density should also be rotationally invariant, and hence can only be a single delta located at the center. However, direct computations on rectangle meshes yield more than one delta singularity for sufficiently small R. This is because the meshes are not invariant under rotation. In this paper, following the idea in [9, 10], we construct a special mesh to obtain symmetry and convergence to physically relevant solutions. Computational results are given to demonstrate the advantages of high order DG schemes. For the second example, we discuss the pressureless Euler equation wt + f(w)x = 0, t > 0, x ∈ R, m ρ , , f(w) = w= ρu2 m
(1.3)
with m = ρu, where ρ is the density function and u is the velocity. Pressureless Euler equations in one space dimension have been analyzed at the theoretical level intensively, e.g. [2, 3, 6, 8, 16]. Some numerical methods have also been studied by several authors [4, 1, 11, 5]. 3
In [4], only first and second order numerical schemes were considered. Except for those in [4], no other methods seem to have been designed to solve the equation (1.3) directly. In [5], the authors added an artificial viscosity and built a diffusive scheme. In [11], the authors applied the sticky particle methods to the equation and showed the approximation satisfies the original system within a certain residual. In [1], the authors introduced a new variable and added one more equation to the system, leading to more computational cost. In this paper, we consider high order DG scheme and approximate the equation without modification. Physically, the density ρ is positive and the velocity u satisfies a maximum principle. We extend the idea in [26] and construct suitable limiters to fulfill these two requirements while maintaining high order accuracy. Moreover, numerical evidences also demonstrate that the scheme is good for approximations in the presence of vacuum. Finally, our scheme works well in two dimensions. To the authors’ knowledge, few works in the literature focus on the computational aspects of two dimensional equations, although some theoretical results can be found in [20, 21]. It appears that complete existence and uniqueness results are not available. Therefore, our scheme should be a good tool to study two-dimensional pressureless Euler equations and other similar equations. The organization of this paper is as follows. In section 2, we present some preliminaries, including the implementation of the DG scheme, the foundation of the limiters, and high order time discretizations. In sections 3 and 4, we discuss the two models in detail, and numerical evidences are given to demonstrate the advantages of the DG methods. Finally, we will end in section 5 with some concluding remarks and remarks for future works.
2 2.1
Preliminaries The DG scheme
First, we divide the computational domain Ω = [0, 1] into N cells 0 = x 1 < x 3 < · · · < xN + 1 = 1, 2
2
2
4
and denote
Ij = xj− 1 , xj+ 1 2
2
as the cells. For simplicity, we consider uniform mesh in this paper, and denote ∆x as the length of each cell. Next, we define Vh = v : each of its components vi |Ij ∈ P k (Ij ), j = 1, · · · , N
as the finite element space, where P k (Ij ) denotes the space of polynomials in Ij of degree at most k. The DG scheme for (1.1) is the following: find uh ∈ Vh , such that for any vh ∈ Vh ((uh )t , vh )j = (f(uh ), (vh )x )j + ˆ fj− 1 vh+ |j− 1 − ˆfj+ 1 vh− |j+ 1 , 2
where (w, v)j =
R
Ij
xj+ 1 . Likewise for 2
2
2
2
(2.1)
) denotes the left limit of the vector vh at wvdx, and vh− |j+ 1 = vh (x− j+ 1 2
vh+ .
2
Moreover, the numerical flux ˆ f is a single valued vector defined at
the cell interfaces and in general depends on the values of the numerical solution uh from both sides of the interfaces ˆ fi+ 1 = ˆ f (uh (x− ), uh (x+ )). i+ 1 i+ 1 2
2
2
In general, for scalar equations, we use monotone fluxes.
2.2
Limiters
In this subsection, we use Euler-forward for time discretization and briefly discuss the construction of bound-preserving limiters [27]. We denote unj and unj to be the numerical solution and its cell average at time level n in cell Ij . For simplicity, throughout the paper, if we consider generic numerical solution on the whole computational domain Ω, then the subscript j will be omitted. Suppose the exact solution of equation (1.1) is in some convex set G, we are interested in constructing numerical solutions which are also in G. The whole precess can be divided into three steps.
5
In the first step, we consider a first order scheme n n n n n ˆ ˆ un+1 = u + λ f u , u − f u , u j j−1 j j j+1 j 1 n 1 n n n n n ˆ ˆ u + 2λf uj−1 , uj + u − 2λf uj , uj+1 = 2 j 2 j 1 1 = H1 unj−1 , unj , 2λ + H2 unj , unj+1 , 2λ , 2 2
where
H1 (u, v, c) = v + c ˆ f(u, v),
H2 (u, v, c) = u − c ˆ f(u, v).
Here unj = unj is a constant in each cell Ij , and λ =
∆t ∆x
(2.2)
(2.3)
is the ratio of time and space mesh
sizes. For many two-point first order numerical fluxes, we can prove the following property. Property 2.1 Suppose G is a convex set and u, v ∈ G, then there exists a positive constant C⋆ , such that, for any 0 < c < C⋆ , we have H1 (u, v, c) , H2 (u, v, c) ∈ G. Based on the above property, we can easily obtain that, under the CFL condition λ
0}. Clearly, G is a convex set. We start with the following first order scheme n n n n n 1 h(ρ , ρ 1 h(ρ , ρ ) − v ) , ρn+1 = ρ + λ v j−1 j j+ j j+1 j j− j 2
(3.1)
2
where h(·, ·) is a numerical flux, and ρnj = ρnj is the numerical approximation to the exact solution in cell Ij at time level n, with ρnj being its cell average. Moreover, vj− 1 is the 2
numerical velocity at the interface xj− 1 , given by 2
vj− 1 = 2
N Z X i=1
Ii
(y − xj− 1 )ξ(y − xj− 1 )ρni (y)dy. 2
2
For this model, we use an upwind flux, i.e. vh(u, w) =
vu v ≥ 0, vw v < 0,
and define H1 and H2 as H1 (u, w, c) = w + cvh(u, w),
H2 (u, w, c) = u − cvh(u, w).
Then the scheme (3.1) can be written as 1 1 ρn+1 = H1 (ρnj−1 , ρnj , 2λ) + H2 (ρnj , ρnj+1 , 2λ). j 2 2 Based on the above notations, we can prove the following lemma. Lemma 3.1 Suppose ρn > 0, then under the CFL condition 1 max |vj− 1 |λ < , 2 j 2 we have ρn+1 > 0. Proof: We consider H1 (ρnj−1 , ρnj , 2λ) first. If vj− 1 < 0, then 2
H1 (ρnj−1 , ρnj , 2λ) = ρnj + 2λvj− 1 h(ρnj−1 , ρnj ) = (1 + 2λvj− 1 )ρnj > 0. 2
2
8
On the other hand, if vj− 1 ≥ 0, then 2
H1 (ρnj−1 , ρnj , 2λ) = ρnj + 2λvj− 1 h(ρnj−1 , ρnj ) = ρnj + 2λvj− 1 ρnj−1 > 0. 2
2
Similarly, we have H2 (ρnj , ρnj+1, 2λ) = ρnj − 2λvj+ 1 h(ρnj , ρnj+1 ) > 0. 2
Therefore, 1 1 ρn+1 = H1 (ρnj−1 , ρnj , 2λ) + H2 (ρnj , ρnj+1, 2λ) > 0. j 2 2 Now, we consider high order schemes. The analysis in section 2.2 implies the following theorem. Theorem 3.1 Suppose the DG solution ρn > 0, then under the CFL condition max |vj− 1 |λ < α0 , j
2
we have ρ¯n+1 > 0. Based on the above theorem, we can modify the density ρnj in the following steps. • Set up a small number ε = min 10−13 , ρ¯nj .
• Compute mj = mini ρnj (xji ), where {xji } are the Gauss-Lobatto points in cell Ij . • If mj < ε, then we take θ=
ρnj − ε , ρnj − mj
and use ρenj = ρnj + θ(ρnj − ρnj )
(3.2)
as the DG approximation in cell Ij at time level n.
Based on the above steps, the numerical density is always positive. Therefore n
kρ kL1 (Ω) =
Z
n
ρ (x)dx =
Ω
Z
Ω
9
ρ0 (x)dx = kρ0 kL1 (Ω) ,
(3.3)
where kukL1 (Ω) is the standard L1 norm of u on Ω. Clearly, (3.3) implies the L1 stability for the DG scheme. Moreover, we can also derive a sufficient CFL condition which does not depend on the numerical velocity in Theorem 3.1. Actually, Notice the fact that vj− 1 = 2
N Z X i=1
Ii
(y − xj− 1 )ξ(y − xj− 1 )ρni (y)dy ≤ Rkρ0 kL1 (Ω) , 2
2
the sufficient CFL condition is λ≤
α0 . Rkρ0 kL1 (Ω)
(3.4)
To sum up, we have the following theorem. Theorem 3.2 Under the CFL condition (3.4), the DG scheme with the positivity-preserving limiter for equation (1.2) is L1 stable and the density function is always positive.
3.2
Numerical experiments
In this subsection, some numerical examples will be given to demonstrate the good performance of the DG scheme. Example 1. We consider the following problem ρt + (vρ)x = 0, ρ(x, 0) = ρ0 (x),
x ∈ [0, 1], t > 0, t > 0,
(3.5)
where the velocity v is defined by v(x, t) =
Z
x+R x−R
(y − x)ρ(y, t)dy.
We apply the positivity-preserving limiter and use P 0 and P 1 polynomials. Moreover, we use the third order SSP Runge-Kutta discretization in time [23] with ∆t = 0.1∆x. Figure 3.1 shows the numerical approximations of ρ(x) at t = 1000, with N = 400, ρ0 = 1, and R = 0.02. We can observe 22 δ-singularities in each panel, and the distance between any two adjacent singularities is greater than R. The algorithm is quite stable in this simulation. Moreover, the P 1 solution in the right panel is more accurate than the P 0 one in the left panel, since the heights of the δ-singularities are almost doubled. 10
30
30
25
25
20
20
15
15
10
10
5
5
0
0.2
0.4
X
0.6
0
0.8
0.2
0.4
X
0.6
0.8
Figure 3.1: Numerical density for (3.5) at t = 1000 with N = 400 when using P 0 (left) and P 1 (right) polynomials. Example 2. We consider the model problem in two dimensions. x ∈ [−1, 1]2 , t > 0, t > 0,
ρt + div(vρ) = 0, ρ(x, 0) = ρ0 (x),
(3.6)
where the velocity v is defined by v(x, t) =
Z
BR (x)
(y − x)ρ(y, t)dy.
In this example, we take R = 0.1 and ρ0 (x) =
1 r < 0.5, 0 r > 0.5,
where r = kxk is the Euclidean norm of x. In [7], the authors demonstrated that the exact solution should be a single delta placed at the origin. However, by using rectangle meshes, we observe more than one delta singularity for R sufficiently small. This is because the meshes are not invariant under rotation. To tackle this problem, we follow the same ideas in [9, 10], and construct a special mesh, namely equal-angle-zoned mesh. The structure of the mesh is given in figure 3.2. By using such a special mesh, the limit density given in figure 3.3 is a single delta placed at the origin.
11
Figure 3.2: Equal-angle-zoned mesh.
Figure 3.3: Numerical density ρ for (3.6) at t = 2000 with N = 200 when using P 0 polynomials.
12
4
Pressureless Euler equations
In this section, we apply DG methods to pressureless Euler equations.
4.1
Numerical schemes in one dimension
We study equation (1.3) in more detail. Physically, the density is positive and the velocity satisfies the maximum principle. Therefore, we define ρ : ρ > 0, aρ ≤ m ≤ bρ , G= w= m where a = min u0 (x),
b = max u0 (x),
(4.1)
with u0 being the initial velocity. Clearly, G is a convex set. As mentioned in section 2.2, we start with the following first order scheme, n n wjn+1 = wjn + λ h(wj−1 , wjn ) − h(wjn , wj+1 ) ,
where h(·, ·) is a numerical flux and wjn = ρnj , mnj
(4.2)
T
is the numerical approximation to the T as its cell exact solution in cell Ij at time level n. Moreover, we define wnj = ρnj , mnj
average. Clearly, for a first order scheme, wjn = wnj in (4.2). For simplicity, we use unj for
mn j ρn j
as the numerical velocity throughout this section. In this problem, we consider the Godunov flux [4]. Suppose at the cell interface x = xj− 1 we have two numerical approximations 2
T
T
wℓ = (ρℓ , mℓ ) and wr = (ρr , mr ) from left and right respectively. Then the Godunov flux is given as (mℓ , ρℓ u2ℓ ) (0, 0) (mr , ρr u2r ) ρu2 j− 1 ) = (h (wℓ , wr ))T = (c ρuj− 1 , d 2 2 (mℓ , ρℓ u2ℓ ) (m , ρ u2 ) mℓr+mrr r 2 ( 2 , ρℓ uℓ = ρr u2r )
where
mℓ , uℓ = ρℓ
mr ur = , and v = ρr 13
√
uℓ uℓ uℓ uℓ uℓ uℓ
> 0, ur ≤ 0, ur ≤ 0, ur > 0, ur > 0, ur > 0, ur
√ ρℓ uℓ + ρr ur . √ √ ρℓ + ρr
> 0, > 0, ≤ 0, ≤ 0, v > 0, ≤ 0, v < 0, ≤ 0, v = 0,
(4.3)
For this problem, H1 and H2 are taken to be H1 (u, v, c) = v + c h (u, v) ,
H2 (u, v, c) = u − c h (u, v) .
(4.4)
Clearly, the scheme (4.2) can be written as 1 1 n n , wjn , 2λ + H2 wjn , wj+1 , 2λ . wjn+1 = H1 wj−1 2 2
Before proceeding to the theoretical results for the scheme above, we would like to introduce the following lemma, whose proof is trivial and is omitted. Lemma 4.1 Suppose {xi } are positive real numbers, and a ≤ yi ≤ b, ∀i, then Pn xi yi a ≤ Pi=1 ≤ b. n i=1 xi
We will use Lemma 4.1 to prove the following lemma.
Lemma 4.2 Suppose wn ∈ G, then under the CFL condition λ
0 and a ≤ uˇ ≤ b. To do so, we have to determine what {xi } and {yi } should be in Lemma 4.1, by testing the different choices for the numerical flux in (4.3). For simplicity, we define u bj− 1 = 2
d2 1 ρu j− 2
ρc uj− 1 2
, if ρc uj− 1 6= 0. 2
• If ρc uj− 1 = mnj−1 , then u bj− 1 = unj−1 > 0. We take x1 = ρnj , y1 = unj and x2 = 2λmnj−1 , 2
y2 =
2
unj−1.
bj− 1 = unj ≤ 0. We take x1 = ρnj + 2λmnj , y1 = unj . • If ρc uj− 1 = mnj , then u 2
2
14
d2 1 = (mn un + mn un )/2. We combine the two • If ρc uj− 1 = (mnj−1 + mnj )/2, then ρu j−1 j−1 j j j− 2
2
situations above, and take x1 = ρnj + λmnj , y1 = unj and x2 = λmnj−1 , y2 = unj−1.
ρu2 j− 1 = 0. We take x1 = ρnj , y1 = unj . • If ρc uj− 1 = 0, then d 2
2
Clearly, in each case, xi > 0, therefore ρˇ > 0. Moreover, we observe a ≤ yi ≤ b, then by Lemma 4.1, we have a ≤ uˇ ≤ b. Now, we consider high order schemes. By the same analysis as in section 2.2, we have the following theorem. Theorem 4.1 Suppose wn ∈ G, then under the CFL condition λ
ε, then proceed to the following steps. Otherwise, ρnj is identified as the e jn = wnj approximation to vacuum, and the velocity is undefined. Therefore, we take w
as the numerical solution and skip the following steps.
• Modify the density first: Compute mj = mini ρnj (xji ), where {xji } are the Gauss-Lobatto points in cell Ij , and get ρenj by (3.2). Then use ρenj as the new numerical density ρnj . 15
• Modify the velocity: Define qji = wjn (xji ) in cell Ij . If qji ∈ Gε , then take θij = 1. Otherwise, take θij
n
w − sj j i
, = n
wj − qj i
sji
where k · k is the Euclidean norm, and
is the intersection point of the straight line
s(t) = (1 − t)wnj + tqji ,
0 ≤ t ≤ 1,
and the surface ∂Gε . Define θj = mini=0,··· ,m θij , and use e jn = wnj + θj (wjn − wnj ), w
as the DG approximation in cell Ij .
4.2
Numerical schemes in two dimensions
We extend our work to two dimensions and consider the following equation wt + f(w)x + g(w)y = 0,
ρ w = m , n
with
t > 0, (x, y) ∈ R2 ,
m f(w) = ρu2 , ρuv
m = ρu,
(4.6)
n g(w) = ρuv , ρv 2
n = ρv,
where ρ is the density function and (u, v) is the velocity field. We define ρ G = w = m : ρ > 0, m2 + n2 ≤ S 2 ρ2 , n
where
S > 0, and S 2 = max u2 (x, y, 0) + v 2 (x, y, 0) x,y
with (u, v)(x, y, 0) being the initial velocity field. Clearly, G is a convex set.
i h For simplicity, we use uniform rectangular meshes. The cell is defined as Iij = xi− 1 , xi+ 1 × 2 2 h i yj− 1 , yj+ 1 , and the mesh sizes in x and y directions are denoted as ∆x and ∆y respectively. 2
2
16
At time level n, we approximate the exact solution with a vector of polynomials of degree k, n wij = (ρnij , mnij , nnij )T , and define the cell average wnij = (ρnij , mnij , nnij )T . Moreover, we denote + − + − wi− (y), wi,j− (x) as the traces of w on the four edges in cell Iij re1 (y), w 1 (x), w ,j i+ 1 ,j i,j+ 1 2
2
2
2
spectively. More details can be found in [26]. In this subsection, we always use (unij , vijn ) for mn
( ρnij , ij
nn ij ) ρn ij
as the numerical velocity field in cell Iij at time level n, and define a1 = maxij |unij |
and a2 = maxij |vijn |. For simplicity, if we consider a generic numerical solution on the whole computational domain at time level n, then the subscript ij will be omitted. In this section, we only consider high order schemes, and the one satisfied by the cell averages can be written as wn+1 ij
Z y 1 j+ 2 ∆t + − + − = + h1 wi− 1 ,j (y), wi− 1 ,j (y) − h1 wi+ 1 ,j (y), wi+ 1 ,j (y) dy 2 2 2 2 ∆x∆y y 1 j− 2 Z x 1 i+ 2 ∆t + − + − − h w dx, (4.7) h2 wi,j− + 1 (x) 1 (x), w 1 (x) 1 (x), w 2 i,j− 2 i,j+ 2 i,j+ 2 2 ∆x∆y x 1 wnij
i− 2
where h1 (·, ·) and h2 (·, ·) are one-dimensional numerical fluxes. For this problem, we also use the Godunov flux. Suppose (x, y) = (xi− 1 , y0 ) is a point on the vertical cell interface, 2
at which we have two numerical approximations wℓ = (ρℓ , mℓ , nℓ )T and wr = (ρr , mr , nr )T from left and right respectively. Then the Godunov flux (h1 (wℓ , wr ))T can be written as (mℓ , ρℓ u2ℓ , ρℓ uℓ vℓ ) uℓ > 0, ur > 0, uℓ ≤ 0, ur > 0, (0, 0, 0) 2 (m , ρ u , ρ u v ) uℓ ≤ 0, ur ≤ 0, r r r r r r ρc u, d ρu2 , ρd uv = 2 (mℓ , ρℓ uℓ , ρℓ uℓ vℓ ) uℓ > 0, ur ≤ 0, v > 0, 2 (m , ρ u , ρ u v ) uℓ > 0, ur ≤ 0, v < 0, 1 r r r r r 2r 2 (m + m , ρ u + ρ u , m v + m v ) u ℓ r ℓ ℓ r r ℓ ℓ r r ℓ > 0, ur ≤ 0, v = 0, 2 where
√ √ ρℓ uℓ + ρr ur mℓ nℓ mr nℓ (uℓ , vℓ ) = , (ur , vr ) = , and v = √ , , . √ ρℓ ρℓ ρr ρℓ ρℓ + ρr T c2 The numerical flux h2 = ρc v, ρd uv, ρv can be defined in a similar way on the horizontal
cell interfaces.
For accuracy, we use L-point Gauss quadratures with L ≥ k + 1 to approximate the integrals in (4.7). More details of this requirement can be found in [12]. The Gauss quadrature 17
i h i h points on xi− 1 , xi+ 1 and yj− 1 , yj+ 1 are denoted by 2
2
pxi
2
=
n
2
o n o y β : β = 1, · · · , L and pj = yj : β = 1, · · · , L ,
xβi
respectively. Also, we denote wβ as the corresponding weights on the interval − 21 , 12 .
Different from the notations in previous sections, we use
pˆxi = {ˆ xαi : α = 0, · · · , M} and pˆyj = yˆjα : α = 0, · · · , M
h i h i as the Gauss-Lobatto points on xi− 1 , xi+ 1 and yj− 1 , yj+ 1 respectively. Also, we denote 2 2 2 2 1 1 wˆα as the corresponding weights on the interval − 2 , 2 . ∆t ∆x
and λ2 =
wn+1 ij
wnij
Let λ1 =
=
∆t , ∆y
+ λ1
then the numerical scheme (4.7) becomes
L X β=1
+ λ2
L X β=1
h i − + + − wβ h1 wi− 1 ,β , wi− 1 ,β − h1 wi+ 1 ,β , wi+ 1 ,β 2
2
2
2
i h − + + − − h w , , w , w wβ h2 wβ,j− 1 2 β,j+ 1 β,j− 1 β,j+ 1 2
2
2
(4.8)
2
β − − = wi− where wi− 1 (yj ) is a point value in the Gauss quadrature. Likewise for the other 1 ,β ,j 2
2
point values. As the general treatment, we rewrite the cell average on the right hand side as wnij
=
M X L X
1 wˆα wβ wαβ
=
M X L X
2 wˆα wβ wβα ,
α=0 β=1
α=0 β=1
1 2 n n where wαβ and wβα denote wij (ˆ xαi , yjβ ) and wij (xβi , yˆjα) respectively. We extend the defini-
tions of H1 and H2 in (4.4) to two-dimensional problems and define H11 (u, v, c) = v + c h1 (u, v) ,
H12 (u, v, c) = u − c h1 (u, v) ,
H21 (u, v, c) = v + c h2 (u, v) ,
H22 (u, v, c) = u − c h2 (u, v) .
Let µ = a1 λ1 + a2 λ2 , then scheme (4.8) can be written as wn+1 = C1 ij
L X
wβ
α=1
β=1
+ C2
L X β=1
M −1 X
wβ
M −1 X α=1
1 wˆα wαβ
! + − + − , wi− , µ1 + wˆM H12 wi+ , wi+ , µ1 + wˆ0 H11 wi− 1 1 1 1 ,β ,β ,β ,β 2
2
2
2
! − − + + 2 , µ2 + wˆM H22 wβ,i+ , µ2 , wˆα wβα + wˆ0 H21 wβ,j− 1,w 1,w β,j− 1 β,j+ 1 2
18
2
2
2
where C1 =
a2 λ2 µ µ a1 λ1 , C2 = , µ1 = , µ2 = . µ µ a1 wˆ0 a2 wˆ0
Now, we can state the main theorem. Theorem 4.2 Suppose wn ∈ G, then under the CFL condition ∆t ∆t a1 + a2 ≤ w ˆ0 , ∆x ∆y we have wn+1 ∈ G. Proof: For simplicity, we only prove H11
H11
+ − , wi− , µ1 wi− 1 1 ,β ,β 2 2
− + wi− , wi− , µ1 1 1 ,β ,β 2 2
= (ˇ ρ, m, ˇ n ˇ )T , uˇ =
∈ G, ∀β, and define n ˇ m ˇ , vˇ = . ρˇ ρˇ
Following the same analysis as in Lemma 4.2, we have ρˇ > 0. Therefore, we need only prove uˇ2 + vˇ2 ≤ S 2 . By the assumption, we have − wi− 1 ,β 2
=
ρ− , m− , n− i− 12 ,β i− 21 ,β i− 21 ,β
+ − , wi− Denote h1 wi− 1 1 ,β ,β 2
2
T
∈ G and
+ wi− 1 ,β 2
T + + + = ρi− 1 ,β , mi− 1 ,β , ni− 1 ,β ∈ G. 2
2
2
T d 2 = ρc ui− 1 ,β , ρu i− 1 ,β , ρd uv i− 1 ,β as the corresponding numerical 2
2
2
T
flux, and for any unit vector n = (n1 , n2 ) , define w ˇ = uˇn1 + vˇn2 . Then + d 2 1 n +ρ 1 n + n n + µ ρu d uv n m+ 1 i− 2 ,β 1 i− 2 ,β 2 i− 12 ,β 2 i− 21 ,β 1 wˇ = + ρi− 1 ,β + µ1 ρc ui− 1 ,β 2
2
=
where + wi− 1 ,β 2
=
ρ+ w+ i− 21 ,β i− 21 ,β ρ+ i− 21 ,β
bi− 1 ,β + µ1 ρc ui− 1 ,β w 2
+ µ1 ρc ui− 1 ,β
,
2
n + n+ n m+ i− 1 ,β 2 i− 1 ,β 1 2
2
2
ρ+ i− 12 ,β
,
w b
i− 21 ,β
=
d2 1 n1 + ρd ρu uv i− 1 ,β n2 i− ,β 2
2
ρc u
.
i− 21 ,β
+ We can easily show that |wi− | ≤ S and |w bi− 1 ,β | ≤ S. Following the same lines as the proof 1 ,β 2
2
of Lemma 4.2, we have |w| ˇ ≤ S. Especially, choosing n to be parallel with (ˇ u, vˇ), we have uˇ2 + vˇ2 ≤ S 2 , completing the proof. 19
Remark 4.1 Since a1 ≤ S and a2 ≤ S, another sufficient CFL condition in Theorem 4.2 is ∆t ∆x
+
∆t ∆y
≤
w ˆ0 . S
n Based on the above theorem, we can modify the numerical solution wij keeping the cell
average untouched. Due to the rounding error, we define ρ ε 2 2 2 2 m G = w= : ρ ≥ ε, m + n ≤ (S + ε) ρ , n ρ ∂Gε = w = m : ρ ≥ ε, m2 + n2 = (S + ε)2 ρ2 . n n Then the modification of wij is given in the following steps.
• Set up a small number ε = 10−13 . • If ρnij > ε, then proceed to the following steps. Otherwise, ρnij is identified as the n e ij approximation to vacuum, and the velocity is undefined. Therefore, we take w = wnij
as the numerical solution and skip the following steps.
o n • Modify the density first: Compute mij = minαβ ρnij (b xαi , yjβ ), ρnij (xβi , ybjα) . If mij < ε, then take ρenij as
ρenij = ρnij + θij ρnij − ρnij ,
with
θij =
ρnij − ε , ρnij − mij
and use ρenij as the new numerical density ρnij .
1 2 1 • Modify the velocity: Consider wαβ and wβα in cell Iij respectively. If wαβ ∈ Gε , then 1 take θαβ = 1. Otherwise, take 1 θαβ
n
wij − s1 αβ
=
wn − w1 , ij
αβ
where k · k is the Euclidean norm, and s1αβ is the intersection point of the straight line 1 s1 (t) = (1 − t)wnij + twαβ ,
20
0 ≤ t ≤ 1,
2 2 and the surface ∂Gε . Similarly, we can define θβα in the same way for wβα . Finally, we
use n n e ij − wnij ), w = wnij + θ(wij
as the DG approximation in cell Iij .
4.3
1 2 , θβα , θ = min θαβ α,β
Numerical experiments
In this subsection, we provide numerical experiments to demonstrate the good performance of the DG scheme for solving pressureless Euler equations. In all the numerical simulations, if not otherwise stated, we use third order schemes and take N = 100. 4.3.1
One space dimension
We consider the problem in one space dimension and solve equation (1.3) with different initial conditions. Example 1. We consider the following initial data ρ0 (x) = sin(x) + 2, u0 (x) = sin(x) + 2,
(4.9)
with periodic boundary condition. Clearly, the exact solution is u(x, t) = u0 (x0 ),
ρ(x, t) =
ρ0 (x0 ) , 1 + u′0 (x0 )
where x0 is given implicitly by x0 + tu0 (x0 ) = x. We use the third order SSP multi-step method in time [22] with ∆t = 0.01∆x2 , and test the example by using P k polynomials with k = 1, 2, 3 on uniform meshes. Table 4.1, shows the L2 -norm of the error at t = 0.1. We observe (k + 0.5)-th order convergence. Example 2. We consider the following initial condition ρ0 (x) =
1 x < 0, 0.25 x > 0, 21
u0 (x) =
1 x < 0, 0 x > 0.
(4.10)
Table 4.1: L2 -norm of the error between the numerical density and the exact density for initial condition (4.9). k=1 N error 20 1.41E-02 40 4.18E-03 80 1.30E-03 160 4.24E-04 320 1.51E-04
order 1.76 1.68 1.62 1.49
k=2 error 6.84E-04 1.04E-04 1.55E-05 2.41E-06 3.80E-07
k=3 order error order 3.40e-5 2.72 2.82e-6 3.59 2.74 2.26e-7 3.64 2.69 1.83e-8 3.62 2.67 1.49e-9 3.63
Clearly, the exact solution is (ρ(x, t), u(x, t)) = and at x =
2t , 3
(1, 1) x < 2t/3, (0.25, 0) x > 2t/3,
the density should be a δ-function. Figure 4.1 shows the numerical density
1 15
0.8
0.6
u
10
0.4 5 0.2
0
0 -0.4
-0.2
0
x
0.2
0.4
-0.4
-0.2
0
x
0.2
0.4
Figure 4.1: Numerical density (left) and velocity (right) at t = 0.5 with P 1 polynomials for initial condition (4.10). and velocity at t = 0.5 with P 1 polynomials. From the figure, we observe the numerical solution capture the profile of the exact solution quite well. Example 3. We consider the following initial condition −0.5 x < −0.5, 0.4 −0.5 < x < 0, ρ0 (x) = 0.5, u0 (x) = 0.4 − x 0 < x < 0.8, −0.4 x > 0.8, 22
(4.11)
and the exact solution for t < 1 (ρ(x, t), u(x, t)) =
is (0.5, −0.5) (0, undefined) (0.5, 0.4) 0.5 0.4−x ( 1−t , 1−t ) (0.5, −0.4)
x < −0.5 − 0.5t, −0.5 − 0.5t < x < −0.5 + 0.4t, −0.5 + 0.4t < x < 0.4t, 0.4t < x < 0.8 − 0.4t, x > 0.8 − 0.4t.
Figure 4.2 shows the numerical density and velocity at t = 0.5. From the figure, we can
0.4 1 0.2
0.8
0
u
0.6
0.4
-0.2
0.2 -0.4 0 -1
-0.5
0
x
0.5
1
-1
-0.5
0
x
0.5
1
Figure 4.2: Numerical density (left) and velocity (right) at t = 0.5 for initial condition (4.11). The solid line shows the exact solution while the symbols show the numerical solution. observe some local oscillations near the singularities. This is not surprising as we have not used any limiters other than the bound-preserving ones for the DG scheme. Example 4. We consider the following initial condition ρ0 (x) = 0.5,
u0 (x) =
−0.5 x < 0, 0.4 x > 0,
(4.12)
and the exact solution is x < −0.5t, (0.5, −0.5) (0, undefined) −0.5t < x < 0.4t, (ρ(x, t), u(x, t)) = (0.5, 0.4) x > 0.4t.
Figure 4.3 shows the numerical density and velocity at t = 0.5. From the figure, we can observe some local oscillations near the singularities.
23
0.5
0.4
0.4
0.2
0.3
u
0
0.2 -0.2 0.1 -0.4 0 -0.4
-0.2
0
0.2
x
0.4
-0.4
-0.2
0
0.2
x
0.4
Figure 4.3: Numerical density (left) and velocity (right) at t = 0.5 for initial condition (4.12). 4.3.2
Two dimensions
We consider the problem in two dimensions and solve equation (4.6) with different initial conditions. Example 1. We consider the following initial condition ρ(x, y, 0) = ρ0 (x + y) = exp(sin(x + y)), u(x, y, 0) = u0 (x + y) = 31 (cos(x + y) + 2),
(4.13)
v(x, y, 0) = v0 (x + y) = 31 (sin(x + y) + 2). The exact solution is u(x, y, t) = u0 (z0 ),
v(x, y, t) = v0 (z0 ),
ρ(x, y, t) =
1+
ρ0 (z0 ) ′ u0 (z0 ) +
v0′ (z0 )
,
where z0 is given implicitly by z0 + t(u0 (z0 ) + v0 (z0 )) = x + y. We use the third order SSP multi-step method in time [22] with ∆t = 0.01∆x3/2 , and test the example by using P k polynomials with k = 1, 2, 3. Table 4.2 shows the L2 -norm of the error at t = 0.1. From the table, we again observe about (k + 0.5)-th order convergence. 24
Table 4.2: L2 -norm of the error between the numerical density and the exact density for initial condition (4.13). k=1 N error 10 0.512 20 0.176 40 6.48E-02 80 2.32E-02 160 9.08E-03
order 1.54 1.44 1.48 1.35
k=2 error 0.107 3.12E-02 8.52E-03 1.39E-03 1.92E-04
k=3 order error 3.42E-02 1.78 3.57E-03 1.87 4.86E-04 2.62 3.97E-05 2.86 3.65E-06
order 3.26 2.88 3.61 3.45
Example 2. We consider the following initial condition ρ(x, y, 0) =
1 , 100
(u, v)(x, y, 0) = (−
1 1 cos θ, − sin θ), 10 10
(4.14)
where θ is the polar angle. Since all the particles are moving towards the origin, the density
0.4
y
0.2
0
-0.2
-0.4 -0.4
-0.2
0
x
0.2
0.4
Figure 4.4: Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (4.14). function at t > 0 should be a single delta at the origin. Different from example 2 in section 3.2, we can observe only one delta located at the origin by using rectangle mesh and the result is given in figure 4.4.
25
Example 3. We consider the following initial condition (−0.25, −0.25) 1 (0.25, −0.25) ρ(x, y, 0) = , (u, v)(x, y, 0) = (0.25, 0.25) 10 (−0.25, 0.25)
x > 0, y x < 0, y x < 0, y x > 0, y
> 0, > 0, < 0, < 0.
(4.15)
Figure 4.5 shows the numerical density and velocity field at t = 0.5. From the figure, we
1
0.8
y
0.6
0.4
0.2
0
0
0.2
0.4
x
0.6
0.8
1
Figure 4.5: Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (4.15). can observe δ-singularities located at the origin and two axes. Example 4. We consider the following initial condition 1 (cos θ, sin θ) r < 0.3, ρ(x, y, 0) = , (u, v)(x, y, 0) = (4.16) 1 1 (− 2 cos θ, − 2 sin θ) r > 0.3, 100 p x2 + y 2 and θ is the polar angle. Figure 4.6 shows the numerical density where r = (contour plot) and velocity field at t = 0.5. From the figure, we can observe δ-shocks located on a circle and vacuum inside. Example 5. We consider the following initial condition (0.3, 0.4) (−0.4, 0.3) ρ(x, y, 0) = 0.5, (u, v)(x, y, 0) = (−0.3, −0.4) (0.4, −0.3) 26
x > 0, y x < 0, y x < 0, y x > 0, y
> 0, > 0, < 0, < 0.
(4.17)
0.4
y
0.2
0
-0.2
-0.4 -0.4
-0.2
0
0.2
x
0.4
Figure 4.6: Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (4.16).
1
0.8
y
0.6
0.4
0.2
0
0
0.2
0.4
x
0.6
0.8
1
Figure 4.7: Numerical density (left) and velocity field (right) at t = 0.4 with N = 50 for initial condition (4.17).
27
Figure 4.7 shows the numerical density (contour plot) and velocity field with N = 50 at t = 0.4. From the figure, we can observe that the numerical solution approximates the vacuum quite well.
5
Concluding remarks
In this paper, we apply discontinuous Galerkin (DG) method to solve hyperbolic conservation laws involving δ-singularities. We study Krause’s consensus models and pressureless Euler equations to demonstrate the stability and high resolution of the DG approximations. Moreover, numerical experiments show that the scheme is also good for approximations in the presence of vacuum. In future work we will extend DG methods to other equations involving δ-singularities in wider areas of applications.
References [1] C. Berthon, M. Breuss and M.-O. Titeux, A relaxation scheme for the approximation of the pressureless Euler equations, Numerical Methods for Partial Differential Equations, 22, (2006), 484-505. [2] F. Bouchut, On zero pressure gas dynamics, Advances in Kinetic Theory and Computing, Would Scientific, River Edge, NJ, (1994), 171-190. [3] F. Bouchut and F. James Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness, Communications in Partial Differential Equations, 24, (1999), 2173-2189. [4] F. Bouchut, S. Jin and X. Li, Numerical approximations of pressureless and isothermal gas dynamics, SIAM Journal on Numerical Analysis, 41, (2003), 135-158. [5] L. Boudin and J. Mathiaud, A numerical scheme for the one-dimensional pressureless gases system, Numerical Methods for Partial Differential Equations, 28, (2006), 17291746. 28
[6] Y. Brenier and E. Grenier, Sticky particles and scalar conservation laws, SIAM Journal on Numerical Analysis, 35, (1998), 2317-2328. [7] C. Canuto, F. Fagnani and P. Tilli, An Eulerian approach to the analysis of Krause’s consensus models, SIAM Journal on Control and Optimization, 50, (2012), 243-265. [8] G.-Q. Chen and H. Liu, Formation of δ-shocks and vacuum states in the vanishing pressure limit of solutions to the Euler equations for isentropic fluids, SIAM Journal on Mathematical Analysis, 34, (2003), 925-938. [9] J. Cheng and C.-W. Shu, A cell-centered Lagrangian scheme with the preservation of symmetry and conservation properties for compressible fluid flows in two-dimensional cylindrical geometry, Journal of Computational Physics, 229, (2010), 7191-7206. [10] J. Cheng and C.-W. Shu, Improvement on spherical symmetry in two-dimensional cylindrical coordinates for a class of control volume Lagrangian schemes, Communications in Computational Physics, 11, (2012), 1144-1168. [11] A. Chertock, A. Kurganov and Y. Rykov, A new sticky particle method for pressureless gas dynamics, SIAM Journal on Numerical Analysis, 45, (2007), 2408-2441. [12] 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. [13] 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. [14] 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. 29
[15] 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. [16] W. E., Yu. G. Rykov and Ya. G. Sinai, Generalized variational principles, global weak solutions and behavior with random initial data for systems of conservation laws arising in adhesion particle dynamics, Communications in Mathematical Physics, 177 (1996), 349-380. [17] S. Gottlieb, C.-W. Shu and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Review, 43 (2001), 89-112. [18] C. Johnson and J. Pitk¨aranta, An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation, Mathematics of Computation, 46 (1986), 1-26. [19] 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. [20] Yu. G. Rykov, Propagation of shock wave type singularities in equations of twodimensional zero-pressure gas dynamics, Mathematical Notes, 66 (1999), 628-635. [21] Yu. G. Rykov, On the nonhamiltonian character of shocks in 2-D pressureless gas, Bollettino della Unione Matematica Italiana. Serie VIII. Sezione B. Articoli di Ricerca Matematica, 5 (2002), 55-78. [22] C.-W. Shu, Total-variation-diminishing time discretizations, SIAM Journal on Scientific and Statistical Computing, 9 (1988), 1073-1084. [23] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shockcapturing schemes, Journal of Computational Physics, 77 (1988), 439-471.
30
[24] Y. Yang and C.-W. Shu, Discontinuous Galerkin method for hyperbolic equations involving δ-singularities: negative-order norm error estimates and applications, Numerische Mathematik, to appear. [25] 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. [26] 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. [27] X. Zhang and C.-W. Shu, Maximum-principle-satisfying and positivity-preserving high order schemes for conservation laws: Survey and new developments, Proceedings of the Royal Society A, 467 (2011), 2752-2776.
31
Keywords: δ-singularities, discontinuous Galerkin method, Krause’s model, pressureless gas, positivity preserving, symmetry preserving, maximum principle. 1
Research supported by DOE grant DE-FG02-08ER25863 and NSF grant DMS-1112700. Division of Applied Mathematics, Brown University, Providence, RI 02912. E-mail: yang [email protected] 3 Division of Applied Mathematics, Brown University, Providence, RI 02912. Current address: PNC Financial Services Group, Inc., New York, New York. E-mail:[email protected] 4 Division of Applied Mathematics, Brown University, Providence, RI 02912. E-mail: [email protected] 2
1
1
Introduction
In this paper, we apply discontinuous Galerkin (DG) methods to solve hyperbolic conservation law ut + f(u)x = 0, u(x, 0) = u0 (x),
(x, t) ∈ R × (0, T ], x ∈ R,
(1.1)
and its two dimensional version, where the exact solution u(x, t) contains δ-singularities. Such problems appear often in applications and are difficult to be well approximated numerically. Most of the previous techniques are designed to modify the singularities with smooth kernels in some narrow region, and hence may smear such singularities severely, leading to large errors in the approximation. On the other hand, the DG methods rely on weak formulations and can solve such problems directly, leading to very accurate results. The first work on approximating δ-singularities with DG method is [24], in which the authors proved the superconvergence results for linear hyperbolic equations with singular initial data and singular source terms. Moreover, several numerical examples were also given in [24] to demonstrate the advantages of the DG scheme in approximating δ-singularities for both linear and nonlinear hyperbolic equations. In this paper, we extend the work in [24] and consider applications to two model equations: Krause’s consensus models and pressureless Euler equations. The DG method, first introduced by Reed and Hill in 1973 [19], was extended to solve scalar linear hyperbolic equations by Johnson and Pitk¨aranta [18]. Later, Cockburn et al. studied Runge-Kutta discontinuous Galerkin (RKDG) methods for hyperbolic conservation laws [14, 12, 13, 15]. Recently, in [25], genuinely maximum-principle-satisfying high order DG schemes for scalar equations and two-dimensional incompressible flows in vorticitystreamfunction formulation have been constructed. Subsequently, positivity-preserving high order DG schemes for compressible Euler equations were given in [26]. In this paper, we extend the ideas in [25, 26] to construct bound-preserving high order DG schemes for the Krause’s consensus models and pressureless Euler equations.
2
For the first example, we discuss the following Krause’s consensus model equation ρt + Fx = 0, x ∈ [0, 1], t > 0, ρ(x, 0) = ρ0 (x), t > 0,
(1.2)
where ρ is the density function, which is always positive. The flux F is given by F (x, t) = v(x, t)ρ(x, t), and the velocity v is defined by v(x, t) =
Z
0
1
(y − x)ξ(y − x)ρ(y, t)dy,
where 0 ≤ ξ(x) ≤ 1 is supported on a ball centered at zero with radius R. In [7], Canuto et al. investigated the discretized version of the PDE and proved that when the time t tends to infinity, the density function ρ will converge to some isolated δ-singularities, and the distance between any two of these δ-singularities cannot be less than R. Some computational results are shown in [7] based on a first order finite volume method. For two dimensions, if the initial density is rotationally invariant, the limit density should also be rotationally invariant, and hence can only be a single delta located at the center. However, direct computations on rectangle meshes yield more than one delta singularity for sufficiently small R. This is because the meshes are not invariant under rotation. In this paper, following the idea in [9, 10], we construct a special mesh to obtain symmetry and convergence to physically relevant solutions. Computational results are given to demonstrate the advantages of high order DG schemes. For the second example, we discuss the pressureless Euler equation wt + f(w)x = 0, t > 0, x ∈ R, m ρ , , f(w) = w= ρu2 m
(1.3)
with m = ρu, where ρ is the density function and u is the velocity. Pressureless Euler equations in one space dimension have been analyzed at the theoretical level intensively, e.g. [2, 3, 6, 8, 16]. Some numerical methods have also been studied by several authors [4, 1, 11, 5]. 3
In [4], only first and second order numerical schemes were considered. Except for those in [4], no other methods seem to have been designed to solve the equation (1.3) directly. In [5], the authors added an artificial viscosity and built a diffusive scheme. In [11], the authors applied the sticky particle methods to the equation and showed the approximation satisfies the original system within a certain residual. In [1], the authors introduced a new variable and added one more equation to the system, leading to more computational cost. In this paper, we consider high order DG scheme and approximate the equation without modification. Physically, the density ρ is positive and the velocity u satisfies a maximum principle. We extend the idea in [26] and construct suitable limiters to fulfill these two requirements while maintaining high order accuracy. Moreover, numerical evidences also demonstrate that the scheme is good for approximations in the presence of vacuum. Finally, our scheme works well in two dimensions. To the authors’ knowledge, few works in the literature focus on the computational aspects of two dimensional equations, although some theoretical results can be found in [20, 21]. It appears that complete existence and uniqueness results are not available. Therefore, our scheme should be a good tool to study two-dimensional pressureless Euler equations and other similar equations. The organization of this paper is as follows. In section 2, we present some preliminaries, including the implementation of the DG scheme, the foundation of the limiters, and high order time discretizations. In sections 3 and 4, we discuss the two models in detail, and numerical evidences are given to demonstrate the advantages of the DG methods. Finally, we will end in section 5 with some concluding remarks and remarks for future works.
2 2.1
Preliminaries The DG scheme
First, we divide the computational domain Ω = [0, 1] into N cells 0 = x 1 < x 3 < · · · < xN + 1 = 1, 2
2
2
4
and denote
Ij = xj− 1 , xj+ 1 2
2
as the cells. For simplicity, we consider uniform mesh in this paper, and denote ∆x as the length of each cell. Next, we define Vh = v : each of its components vi |Ij ∈ P k (Ij ), j = 1, · · · , N
as the finite element space, where P k (Ij ) denotes the space of polynomials in Ij of degree at most k. The DG scheme for (1.1) is the following: find uh ∈ Vh , such that for any vh ∈ Vh ((uh )t , vh )j = (f(uh ), (vh )x )j + ˆ fj− 1 vh+ |j− 1 − ˆfj+ 1 vh− |j+ 1 , 2
where (w, v)j =
R
Ij
xj+ 1 . Likewise for 2
2
2
2
(2.1)
) denotes the left limit of the vector vh at wvdx, and vh− |j+ 1 = vh (x− j+ 1 2
vh+ .
2
Moreover, the numerical flux ˆ f is a single valued vector defined at
the cell interfaces and in general depends on the values of the numerical solution uh from both sides of the interfaces ˆ fi+ 1 = ˆ f (uh (x− ), uh (x+ )). i+ 1 i+ 1 2
2
2
In general, for scalar equations, we use monotone fluxes.
2.2
Limiters
In this subsection, we use Euler-forward for time discretization and briefly discuss the construction of bound-preserving limiters [27]. We denote unj and unj to be the numerical solution and its cell average at time level n in cell Ij . For simplicity, throughout the paper, if we consider generic numerical solution on the whole computational domain Ω, then the subscript j will be omitted. Suppose the exact solution of equation (1.1) is in some convex set G, we are interested in constructing numerical solutions which are also in G. The whole precess can be divided into three steps.
5
In the first step, we consider a first order scheme n n n n n ˆ ˆ un+1 = u + λ f u , u − f u , u j j−1 j j j+1 j 1 n 1 n n n n n ˆ ˆ u + 2λf uj−1 , uj + u − 2λf uj , uj+1 = 2 j 2 j 1 1 = H1 unj−1 , unj , 2λ + H2 unj , unj+1 , 2λ , 2 2
where
H1 (u, v, c) = v + c ˆ f(u, v),
H2 (u, v, c) = u − c ˆ f(u, v).
Here unj = unj is a constant in each cell Ij , and λ =
∆t ∆x
(2.2)
(2.3)
is the ratio of time and space mesh
sizes. For many two-point first order numerical fluxes, we can prove the following property. Property 2.1 Suppose G is a convex set and u, v ∈ G, then there exists a positive constant C⋆ , such that, for any 0 < c < C⋆ , we have H1 (u, v, c) , H2 (u, v, c) ∈ G. Based on the above property, we can easily obtain that, under the CFL condition λ
0}. Clearly, G is a convex set. We start with the following first order scheme n n n n n 1 h(ρ , ρ 1 h(ρ , ρ ) − v ) , ρn+1 = ρ + λ v j−1 j j+ j j+1 j j− j 2
(3.1)
2
where h(·, ·) is a numerical flux, and ρnj = ρnj is the numerical approximation to the exact solution in cell Ij at time level n, with ρnj being its cell average. Moreover, vj− 1 is the 2
numerical velocity at the interface xj− 1 , given by 2
vj− 1 = 2
N Z X i=1
Ii
(y − xj− 1 )ξ(y − xj− 1 )ρni (y)dy. 2
2
For this model, we use an upwind flux, i.e. vh(u, w) =
vu v ≥ 0, vw v < 0,
and define H1 and H2 as H1 (u, w, c) = w + cvh(u, w),
H2 (u, w, c) = u − cvh(u, w).
Then the scheme (3.1) can be written as 1 1 ρn+1 = H1 (ρnj−1 , ρnj , 2λ) + H2 (ρnj , ρnj+1 , 2λ). j 2 2 Based on the above notations, we can prove the following lemma. Lemma 3.1 Suppose ρn > 0, then under the CFL condition 1 max |vj− 1 |λ < , 2 j 2 we have ρn+1 > 0. Proof: We consider H1 (ρnj−1 , ρnj , 2λ) first. If vj− 1 < 0, then 2
H1 (ρnj−1 , ρnj , 2λ) = ρnj + 2λvj− 1 h(ρnj−1 , ρnj ) = (1 + 2λvj− 1 )ρnj > 0. 2
2
8
On the other hand, if vj− 1 ≥ 0, then 2
H1 (ρnj−1 , ρnj , 2λ) = ρnj + 2λvj− 1 h(ρnj−1 , ρnj ) = ρnj + 2λvj− 1 ρnj−1 > 0. 2
2
Similarly, we have H2 (ρnj , ρnj+1, 2λ) = ρnj − 2λvj+ 1 h(ρnj , ρnj+1 ) > 0. 2
Therefore, 1 1 ρn+1 = H1 (ρnj−1 , ρnj , 2λ) + H2 (ρnj , ρnj+1, 2λ) > 0. j 2 2 Now, we consider high order schemes. The analysis in section 2.2 implies the following theorem. Theorem 3.1 Suppose the DG solution ρn > 0, then under the CFL condition max |vj− 1 |λ < α0 , j
2
we have ρ¯n+1 > 0. Based on the above theorem, we can modify the density ρnj in the following steps. • Set up a small number ε = min 10−13 , ρ¯nj .
• Compute mj = mini ρnj (xji ), where {xji } are the Gauss-Lobatto points in cell Ij . • If mj < ε, then we take θ=
ρnj − ε , ρnj − mj
and use ρenj = ρnj + θ(ρnj − ρnj )
(3.2)
as the DG approximation in cell Ij at time level n.
Based on the above steps, the numerical density is always positive. Therefore n
kρ kL1 (Ω) =
Z
n
ρ (x)dx =
Ω
Z
Ω
9
ρ0 (x)dx = kρ0 kL1 (Ω) ,
(3.3)
where kukL1 (Ω) is the standard L1 norm of u on Ω. Clearly, (3.3) implies the L1 stability for the DG scheme. Moreover, we can also derive a sufficient CFL condition which does not depend on the numerical velocity in Theorem 3.1. Actually, Notice the fact that vj− 1 = 2
N Z X i=1
Ii
(y − xj− 1 )ξ(y − xj− 1 )ρni (y)dy ≤ Rkρ0 kL1 (Ω) , 2
2
the sufficient CFL condition is λ≤
α0 . Rkρ0 kL1 (Ω)
(3.4)
To sum up, we have the following theorem. Theorem 3.2 Under the CFL condition (3.4), the DG scheme with the positivity-preserving limiter for equation (1.2) is L1 stable and the density function is always positive.
3.2
Numerical experiments
In this subsection, some numerical examples will be given to demonstrate the good performance of the DG scheme. Example 1. We consider the following problem ρt + (vρ)x = 0, ρ(x, 0) = ρ0 (x),
x ∈ [0, 1], t > 0, t > 0,
(3.5)
where the velocity v is defined by v(x, t) =
Z
x+R x−R
(y − x)ρ(y, t)dy.
We apply the positivity-preserving limiter and use P 0 and P 1 polynomials. Moreover, we use the third order SSP Runge-Kutta discretization in time [23] with ∆t = 0.1∆x. Figure 3.1 shows the numerical approximations of ρ(x) at t = 1000, with N = 400, ρ0 = 1, and R = 0.02. We can observe 22 δ-singularities in each panel, and the distance between any two adjacent singularities is greater than R. The algorithm is quite stable in this simulation. Moreover, the P 1 solution in the right panel is more accurate than the P 0 one in the left panel, since the heights of the δ-singularities are almost doubled. 10
30
30
25
25
20
20
15
15
10
10
5
5
0
0.2
0.4
X
0.6
0
0.8
0.2
0.4
X
0.6
0.8
Figure 3.1: Numerical density for (3.5) at t = 1000 with N = 400 when using P 0 (left) and P 1 (right) polynomials. Example 2. We consider the model problem in two dimensions. x ∈ [−1, 1]2 , t > 0, t > 0,
ρt + div(vρ) = 0, ρ(x, 0) = ρ0 (x),
(3.6)
where the velocity v is defined by v(x, t) =
Z
BR (x)
(y − x)ρ(y, t)dy.
In this example, we take R = 0.1 and ρ0 (x) =
1 r < 0.5, 0 r > 0.5,
where r = kxk is the Euclidean norm of x. In [7], the authors demonstrated that the exact solution should be a single delta placed at the origin. However, by using rectangle meshes, we observe more than one delta singularity for R sufficiently small. This is because the meshes are not invariant under rotation. To tackle this problem, we follow the same ideas in [9, 10], and construct a special mesh, namely equal-angle-zoned mesh. The structure of the mesh is given in figure 3.2. By using such a special mesh, the limit density given in figure 3.3 is a single delta placed at the origin.
11
Figure 3.2: Equal-angle-zoned mesh.
Figure 3.3: Numerical density ρ for (3.6) at t = 2000 with N = 200 when using P 0 polynomials.
12
4
Pressureless Euler equations
In this section, we apply DG methods to pressureless Euler equations.
4.1
Numerical schemes in one dimension
We study equation (1.3) in more detail. Physically, the density is positive and the velocity satisfies the maximum principle. Therefore, we define ρ : ρ > 0, aρ ≤ m ≤ bρ , G= w= m where a = min u0 (x),
b = max u0 (x),
(4.1)
with u0 being the initial velocity. Clearly, G is a convex set. As mentioned in section 2.2, we start with the following first order scheme, n n wjn+1 = wjn + λ h(wj−1 , wjn ) − h(wjn , wj+1 ) ,
where h(·, ·) is a numerical flux and wjn = ρnj , mnj
(4.2)
T
is the numerical approximation to the T as its cell exact solution in cell Ij at time level n. Moreover, we define wnj = ρnj , mnj
average. Clearly, for a first order scheme, wjn = wnj in (4.2). For simplicity, we use unj for
mn j ρn j
as the numerical velocity throughout this section. In this problem, we consider the Godunov flux [4]. Suppose at the cell interface x = xj− 1 we have two numerical approximations 2
T
T
wℓ = (ρℓ , mℓ ) and wr = (ρr , mr ) from left and right respectively. Then the Godunov flux is given as (mℓ , ρℓ u2ℓ ) (0, 0) (mr , ρr u2r ) ρu2 j− 1 ) = (h (wℓ , wr ))T = (c ρuj− 1 , d 2 2 (mℓ , ρℓ u2ℓ ) (m , ρ u2 ) mℓr+mrr r 2 ( 2 , ρℓ uℓ = ρr u2r )
where
mℓ , uℓ = ρℓ
mr ur = , and v = ρr 13
√
uℓ uℓ uℓ uℓ uℓ uℓ
> 0, ur ≤ 0, ur ≤ 0, ur > 0, ur > 0, ur > 0, ur
√ ρℓ uℓ + ρr ur . √ √ ρℓ + ρr
> 0, > 0, ≤ 0, ≤ 0, v > 0, ≤ 0, v < 0, ≤ 0, v = 0,
(4.3)
For this problem, H1 and H2 are taken to be H1 (u, v, c) = v + c h (u, v) ,
H2 (u, v, c) = u − c h (u, v) .
(4.4)
Clearly, the scheme (4.2) can be written as 1 1 n n , wjn , 2λ + H2 wjn , wj+1 , 2λ . wjn+1 = H1 wj−1 2 2
Before proceeding to the theoretical results for the scheme above, we would like to introduce the following lemma, whose proof is trivial and is omitted. Lemma 4.1 Suppose {xi } are positive real numbers, and a ≤ yi ≤ b, ∀i, then Pn xi yi a ≤ Pi=1 ≤ b. n i=1 xi
We will use Lemma 4.1 to prove the following lemma.
Lemma 4.2 Suppose wn ∈ G, then under the CFL condition λ
0 and a ≤ uˇ ≤ b. To do so, we have to determine what {xi } and {yi } should be in Lemma 4.1, by testing the different choices for the numerical flux in (4.3). For simplicity, we define u bj− 1 = 2
d2 1 ρu j− 2
ρc uj− 1 2
, if ρc uj− 1 6= 0. 2
• If ρc uj− 1 = mnj−1 , then u bj− 1 = unj−1 > 0. We take x1 = ρnj , y1 = unj and x2 = 2λmnj−1 , 2
y2 =
2
unj−1.
bj− 1 = unj ≤ 0. We take x1 = ρnj + 2λmnj , y1 = unj . • If ρc uj− 1 = mnj , then u 2
2
14
d2 1 = (mn un + mn un )/2. We combine the two • If ρc uj− 1 = (mnj−1 + mnj )/2, then ρu j−1 j−1 j j j− 2
2
situations above, and take x1 = ρnj + λmnj , y1 = unj and x2 = λmnj−1 , y2 = unj−1.
ρu2 j− 1 = 0. We take x1 = ρnj , y1 = unj . • If ρc uj− 1 = 0, then d 2
2
Clearly, in each case, xi > 0, therefore ρˇ > 0. Moreover, we observe a ≤ yi ≤ b, then by Lemma 4.1, we have a ≤ uˇ ≤ b. Now, we consider high order schemes. By the same analysis as in section 2.2, we have the following theorem. Theorem 4.1 Suppose wn ∈ G, then under the CFL condition λ
ε, then proceed to the following steps. Otherwise, ρnj is identified as the e jn = wnj approximation to vacuum, and the velocity is undefined. Therefore, we take w
as the numerical solution and skip the following steps.
• Modify the density first: Compute mj = mini ρnj (xji ), where {xji } are the Gauss-Lobatto points in cell Ij , and get ρenj by (3.2). Then use ρenj as the new numerical density ρnj . 15
• Modify the velocity: Define qji = wjn (xji ) in cell Ij . If qji ∈ Gε , then take θij = 1. Otherwise, take θij
n
w − sj j i
, = n
wj − qj i
sji
where k · k is the Euclidean norm, and
is the intersection point of the straight line
s(t) = (1 − t)wnj + tqji ,
0 ≤ t ≤ 1,
and the surface ∂Gε . Define θj = mini=0,··· ,m θij , and use e jn = wnj + θj (wjn − wnj ), w
as the DG approximation in cell Ij .
4.2
Numerical schemes in two dimensions
We extend our work to two dimensions and consider the following equation wt + f(w)x + g(w)y = 0,
ρ w = m , n
with
t > 0, (x, y) ∈ R2 ,
m f(w) = ρu2 , ρuv
m = ρu,
(4.6)
n g(w) = ρuv , ρv 2
n = ρv,
where ρ is the density function and (u, v) is the velocity field. We define ρ G = w = m : ρ > 0, m2 + n2 ≤ S 2 ρ2 , n
where
S > 0, and S 2 = max u2 (x, y, 0) + v 2 (x, y, 0) x,y
with (u, v)(x, y, 0) being the initial velocity field. Clearly, G is a convex set.
i h For simplicity, we use uniform rectangular meshes. The cell is defined as Iij = xi− 1 , xi+ 1 × 2 2 h i yj− 1 , yj+ 1 , and the mesh sizes in x and y directions are denoted as ∆x and ∆y respectively. 2
2
16
At time level n, we approximate the exact solution with a vector of polynomials of degree k, n wij = (ρnij , mnij , nnij )T , and define the cell average wnij = (ρnij , mnij , nnij )T . Moreover, we denote + − + − wi− (y), wi,j− (x) as the traces of w on the four edges in cell Iij re1 (y), w 1 (x), w ,j i+ 1 ,j i,j+ 1 2
2
2
2
spectively. More details can be found in [26]. In this subsection, we always use (unij , vijn ) for mn
( ρnij , ij
nn ij ) ρn ij
as the numerical velocity field in cell Iij at time level n, and define a1 = maxij |unij |
and a2 = maxij |vijn |. For simplicity, if we consider a generic numerical solution on the whole computational domain at time level n, then the subscript ij will be omitted. In this section, we only consider high order schemes, and the one satisfied by the cell averages can be written as wn+1 ij
Z y 1 j+ 2 ∆t + − + − = + h1 wi− 1 ,j (y), wi− 1 ,j (y) − h1 wi+ 1 ,j (y), wi+ 1 ,j (y) dy 2 2 2 2 ∆x∆y y 1 j− 2 Z x 1 i+ 2 ∆t + − + − − h w dx, (4.7) h2 wi,j− + 1 (x) 1 (x), w 1 (x) 1 (x), w 2 i,j− 2 i,j+ 2 i,j+ 2 2 ∆x∆y x 1 wnij
i− 2
where h1 (·, ·) and h2 (·, ·) are one-dimensional numerical fluxes. For this problem, we also use the Godunov flux. Suppose (x, y) = (xi− 1 , y0 ) is a point on the vertical cell interface, 2
at which we have two numerical approximations wℓ = (ρℓ , mℓ , nℓ )T and wr = (ρr , mr , nr )T from left and right respectively. Then the Godunov flux (h1 (wℓ , wr ))T can be written as (mℓ , ρℓ u2ℓ , ρℓ uℓ vℓ ) uℓ > 0, ur > 0, uℓ ≤ 0, ur > 0, (0, 0, 0) 2 (m , ρ u , ρ u v ) uℓ ≤ 0, ur ≤ 0, r r r r r r ρc u, d ρu2 , ρd uv = 2 (mℓ , ρℓ uℓ , ρℓ uℓ vℓ ) uℓ > 0, ur ≤ 0, v > 0, 2 (m , ρ u , ρ u v ) uℓ > 0, ur ≤ 0, v < 0, 1 r r r r r 2r 2 (m + m , ρ u + ρ u , m v + m v ) u ℓ r ℓ ℓ r r ℓ ℓ r r ℓ > 0, ur ≤ 0, v = 0, 2 where
√ √ ρℓ uℓ + ρr ur mℓ nℓ mr nℓ (uℓ , vℓ ) = , (ur , vr ) = , and v = √ , , . √ ρℓ ρℓ ρr ρℓ ρℓ + ρr T c2 The numerical flux h2 = ρc v, ρd uv, ρv can be defined in a similar way on the horizontal
cell interfaces.
For accuracy, we use L-point Gauss quadratures with L ≥ k + 1 to approximate the integrals in (4.7). More details of this requirement can be found in [12]. The Gauss quadrature 17
i h i h points on xi− 1 , xi+ 1 and yj− 1 , yj+ 1 are denoted by 2
2
pxi
2
=
n
2
o n o y β : β = 1, · · · , L and pj = yj : β = 1, · · · , L ,
xβi
respectively. Also, we denote wβ as the corresponding weights on the interval − 21 , 12 .
Different from the notations in previous sections, we use
pˆxi = {ˆ xαi : α = 0, · · · , M} and pˆyj = yˆjα : α = 0, · · · , M
h i h i as the Gauss-Lobatto points on xi− 1 , xi+ 1 and yj− 1 , yj+ 1 respectively. Also, we denote 2 2 2 2 1 1 wˆα as the corresponding weights on the interval − 2 , 2 . ∆t ∆x
and λ2 =
wn+1 ij
wnij
Let λ1 =
=
∆t , ∆y
+ λ1
then the numerical scheme (4.7) becomes
L X β=1
+ λ2
L X β=1
h i − + + − wβ h1 wi− 1 ,β , wi− 1 ,β − h1 wi+ 1 ,β , wi+ 1 ,β 2
2
2
2
i h − + + − − h w , , w , w wβ h2 wβ,j− 1 2 β,j+ 1 β,j− 1 β,j+ 1 2
2
2
(4.8)
2
β − − = wi− where wi− 1 (yj ) is a point value in the Gauss quadrature. Likewise for the other 1 ,β ,j 2
2
point values. As the general treatment, we rewrite the cell average on the right hand side as wnij
=
M X L X
1 wˆα wβ wαβ
=
M X L X
2 wˆα wβ wβα ,
α=0 β=1
α=0 β=1
1 2 n n where wαβ and wβα denote wij (ˆ xαi , yjβ ) and wij (xβi , yˆjα) respectively. We extend the defini-
tions of H1 and H2 in (4.4) to two-dimensional problems and define H11 (u, v, c) = v + c h1 (u, v) ,
H12 (u, v, c) = u − c h1 (u, v) ,
H21 (u, v, c) = v + c h2 (u, v) ,
H22 (u, v, c) = u − c h2 (u, v) .
Let µ = a1 λ1 + a2 λ2 , then scheme (4.8) can be written as wn+1 = C1 ij
L X
wβ
α=1
β=1
+ C2
L X β=1
M −1 X
wβ
M −1 X α=1
1 wˆα wαβ
! + − + − , wi− , µ1 + wˆM H12 wi+ , wi+ , µ1 + wˆ0 H11 wi− 1 1 1 1 ,β ,β ,β ,β 2
2
2
2
! − − + + 2 , µ2 + wˆM H22 wβ,i+ , µ2 , wˆα wβα + wˆ0 H21 wβ,j− 1,w 1,w β,j− 1 β,j+ 1 2
18
2
2
2
where C1 =
a2 λ2 µ µ a1 λ1 , C2 = , µ1 = , µ2 = . µ µ a1 wˆ0 a2 wˆ0
Now, we can state the main theorem. Theorem 4.2 Suppose wn ∈ G, then under the CFL condition ∆t ∆t a1 + a2 ≤ w ˆ0 , ∆x ∆y we have wn+1 ∈ G. Proof: For simplicity, we only prove H11
H11
+ − , wi− , µ1 wi− 1 1 ,β ,β 2 2
− + wi− , wi− , µ1 1 1 ,β ,β 2 2
= (ˇ ρ, m, ˇ n ˇ )T , uˇ =
∈ G, ∀β, and define n ˇ m ˇ , vˇ = . ρˇ ρˇ
Following the same analysis as in Lemma 4.2, we have ρˇ > 0. Therefore, we need only prove uˇ2 + vˇ2 ≤ S 2 . By the assumption, we have − wi− 1 ,β 2
=
ρ− , m− , n− i− 12 ,β i− 21 ,β i− 21 ,β
+ − , wi− Denote h1 wi− 1 1 ,β ,β 2
2
T
∈ G and
+ wi− 1 ,β 2
T + + + = ρi− 1 ,β , mi− 1 ,β , ni− 1 ,β ∈ G. 2
2
2
T d 2 = ρc ui− 1 ,β , ρu i− 1 ,β , ρd uv i− 1 ,β as the corresponding numerical 2
2
2
T
flux, and for any unit vector n = (n1 , n2 ) , define w ˇ = uˇn1 + vˇn2 . Then + d 2 1 n +ρ 1 n + n n + µ ρu d uv n m+ 1 i− 2 ,β 1 i− 2 ,β 2 i− 12 ,β 2 i− 21 ,β 1 wˇ = + ρi− 1 ,β + µ1 ρc ui− 1 ,β 2
2
=
where + wi− 1 ,β 2
=
ρ+ w+ i− 21 ,β i− 21 ,β ρ+ i− 21 ,β
bi− 1 ,β + µ1 ρc ui− 1 ,β w 2
+ µ1 ρc ui− 1 ,β
,
2
n + n+ n m+ i− 1 ,β 2 i− 1 ,β 1 2
2
2
ρ+ i− 12 ,β
,
w b
i− 21 ,β
=
d2 1 n1 + ρd ρu uv i− 1 ,β n2 i− ,β 2
2
ρc u
.
i− 21 ,β
+ We can easily show that |wi− | ≤ S and |w bi− 1 ,β | ≤ S. Following the same lines as the proof 1 ,β 2
2
of Lemma 4.2, we have |w| ˇ ≤ S. Especially, choosing n to be parallel with (ˇ u, vˇ), we have uˇ2 + vˇ2 ≤ S 2 , completing the proof. 19
Remark 4.1 Since a1 ≤ S and a2 ≤ S, another sufficient CFL condition in Theorem 4.2 is ∆t ∆x
+
∆t ∆y
≤
w ˆ0 . S
n Based on the above theorem, we can modify the numerical solution wij keeping the cell
average untouched. Due to the rounding error, we define ρ ε 2 2 2 2 m G = w= : ρ ≥ ε, m + n ≤ (S + ε) ρ , n ρ ∂Gε = w = m : ρ ≥ ε, m2 + n2 = (S + ε)2 ρ2 . n n Then the modification of wij is given in the following steps.
• Set up a small number ε = 10−13 . • If ρnij > ε, then proceed to the following steps. Otherwise, ρnij is identified as the n e ij approximation to vacuum, and the velocity is undefined. Therefore, we take w = wnij
as the numerical solution and skip the following steps.
o n • Modify the density first: Compute mij = minαβ ρnij (b xαi , yjβ ), ρnij (xβi , ybjα) . If mij < ε, then take ρenij as
ρenij = ρnij + θij ρnij − ρnij ,
with
θij =
ρnij − ε , ρnij − mij
and use ρenij as the new numerical density ρnij .
1 2 1 • Modify the velocity: Consider wαβ and wβα in cell Iij respectively. If wαβ ∈ Gε , then 1 take θαβ = 1. Otherwise, take 1 θαβ
n
wij − s1 αβ
=
wn − w1 , ij
αβ
where k · k is the Euclidean norm, and s1αβ is the intersection point of the straight line 1 s1 (t) = (1 − t)wnij + twαβ ,
20
0 ≤ t ≤ 1,
2 2 and the surface ∂Gε . Similarly, we can define θβα in the same way for wβα . Finally, we
use n n e ij − wnij ), w = wnij + θ(wij
as the DG approximation in cell Iij .
4.3
1 2 , θβα , θ = min θαβ α,β
Numerical experiments
In this subsection, we provide numerical experiments to demonstrate the good performance of the DG scheme for solving pressureless Euler equations. In all the numerical simulations, if not otherwise stated, we use third order schemes and take N = 100. 4.3.1
One space dimension
We consider the problem in one space dimension and solve equation (1.3) with different initial conditions. Example 1. We consider the following initial data ρ0 (x) = sin(x) + 2, u0 (x) = sin(x) + 2,
(4.9)
with periodic boundary condition. Clearly, the exact solution is u(x, t) = u0 (x0 ),
ρ(x, t) =
ρ0 (x0 ) , 1 + u′0 (x0 )
where x0 is given implicitly by x0 + tu0 (x0 ) = x. We use the third order SSP multi-step method in time [22] with ∆t = 0.01∆x2 , and test the example by using P k polynomials with k = 1, 2, 3 on uniform meshes. Table 4.1, shows the L2 -norm of the error at t = 0.1. We observe (k + 0.5)-th order convergence. Example 2. We consider the following initial condition ρ0 (x) =
1 x < 0, 0.25 x > 0, 21
u0 (x) =
1 x < 0, 0 x > 0.
(4.10)
Table 4.1: L2 -norm of the error between the numerical density and the exact density for initial condition (4.9). k=1 N error 20 1.41E-02 40 4.18E-03 80 1.30E-03 160 4.24E-04 320 1.51E-04
order 1.76 1.68 1.62 1.49
k=2 error 6.84E-04 1.04E-04 1.55E-05 2.41E-06 3.80E-07
k=3 order error order 3.40e-5 2.72 2.82e-6 3.59 2.74 2.26e-7 3.64 2.69 1.83e-8 3.62 2.67 1.49e-9 3.63
Clearly, the exact solution is (ρ(x, t), u(x, t)) = and at x =
2t , 3
(1, 1) x < 2t/3, (0.25, 0) x > 2t/3,
the density should be a δ-function. Figure 4.1 shows the numerical density
1 15
0.8
0.6
u
10
0.4 5 0.2
0
0 -0.4
-0.2
0
x
0.2
0.4
-0.4
-0.2
0
x
0.2
0.4
Figure 4.1: Numerical density (left) and velocity (right) at t = 0.5 with P 1 polynomials for initial condition (4.10). and velocity at t = 0.5 with P 1 polynomials. From the figure, we observe the numerical solution capture the profile of the exact solution quite well. Example 3. We consider the following initial condition −0.5 x < −0.5, 0.4 −0.5 < x < 0, ρ0 (x) = 0.5, u0 (x) = 0.4 − x 0 < x < 0.8, −0.4 x > 0.8, 22
(4.11)
and the exact solution for t < 1 (ρ(x, t), u(x, t)) =
is (0.5, −0.5) (0, undefined) (0.5, 0.4) 0.5 0.4−x ( 1−t , 1−t ) (0.5, −0.4)
x < −0.5 − 0.5t, −0.5 − 0.5t < x < −0.5 + 0.4t, −0.5 + 0.4t < x < 0.4t, 0.4t < x < 0.8 − 0.4t, x > 0.8 − 0.4t.
Figure 4.2 shows the numerical density and velocity at t = 0.5. From the figure, we can
0.4 1 0.2
0.8
0
u
0.6
0.4
-0.2
0.2 -0.4 0 -1
-0.5
0
x
0.5
1
-1
-0.5
0
x
0.5
1
Figure 4.2: Numerical density (left) and velocity (right) at t = 0.5 for initial condition (4.11). The solid line shows the exact solution while the symbols show the numerical solution. observe some local oscillations near the singularities. This is not surprising as we have not used any limiters other than the bound-preserving ones for the DG scheme. Example 4. We consider the following initial condition ρ0 (x) = 0.5,
u0 (x) =
−0.5 x < 0, 0.4 x > 0,
(4.12)
and the exact solution is x < −0.5t, (0.5, −0.5) (0, undefined) −0.5t < x < 0.4t, (ρ(x, t), u(x, t)) = (0.5, 0.4) x > 0.4t.
Figure 4.3 shows the numerical density and velocity at t = 0.5. From the figure, we can observe some local oscillations near the singularities.
23
0.5
0.4
0.4
0.2
0.3
u
0
0.2 -0.2 0.1 -0.4 0 -0.4
-0.2
0
0.2
x
0.4
-0.4
-0.2
0
0.2
x
0.4
Figure 4.3: Numerical density (left) and velocity (right) at t = 0.5 for initial condition (4.12). 4.3.2
Two dimensions
We consider the problem in two dimensions and solve equation (4.6) with different initial conditions. Example 1. We consider the following initial condition ρ(x, y, 0) = ρ0 (x + y) = exp(sin(x + y)), u(x, y, 0) = u0 (x + y) = 31 (cos(x + y) + 2),
(4.13)
v(x, y, 0) = v0 (x + y) = 31 (sin(x + y) + 2). The exact solution is u(x, y, t) = u0 (z0 ),
v(x, y, t) = v0 (z0 ),
ρ(x, y, t) =
1+
ρ0 (z0 ) ′ u0 (z0 ) +
v0′ (z0 )
,
where z0 is given implicitly by z0 + t(u0 (z0 ) + v0 (z0 )) = x + y. We use the third order SSP multi-step method in time [22] with ∆t = 0.01∆x3/2 , and test the example by using P k polynomials with k = 1, 2, 3. Table 4.2 shows the L2 -norm of the error at t = 0.1. From the table, we again observe about (k + 0.5)-th order convergence. 24
Table 4.2: L2 -norm of the error between the numerical density and the exact density for initial condition (4.13). k=1 N error 10 0.512 20 0.176 40 6.48E-02 80 2.32E-02 160 9.08E-03
order 1.54 1.44 1.48 1.35
k=2 error 0.107 3.12E-02 8.52E-03 1.39E-03 1.92E-04
k=3 order error 3.42E-02 1.78 3.57E-03 1.87 4.86E-04 2.62 3.97E-05 2.86 3.65E-06
order 3.26 2.88 3.61 3.45
Example 2. We consider the following initial condition ρ(x, y, 0) =
1 , 100
(u, v)(x, y, 0) = (−
1 1 cos θ, − sin θ), 10 10
(4.14)
where θ is the polar angle. Since all the particles are moving towards the origin, the density
0.4
y
0.2
0
-0.2
-0.4 -0.4
-0.2
0
x
0.2
0.4
Figure 4.4: Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (4.14). function at t > 0 should be a single delta at the origin. Different from example 2 in section 3.2, we can observe only one delta located at the origin by using rectangle mesh and the result is given in figure 4.4.
25
Example 3. We consider the following initial condition (−0.25, −0.25) 1 (0.25, −0.25) ρ(x, y, 0) = , (u, v)(x, y, 0) = (0.25, 0.25) 10 (−0.25, 0.25)
x > 0, y x < 0, y x < 0, y x > 0, y
> 0, > 0, < 0, < 0.
(4.15)
Figure 4.5 shows the numerical density and velocity field at t = 0.5. From the figure, we
1
0.8
y
0.6
0.4
0.2
0
0
0.2
0.4
x
0.6
0.8
1
Figure 4.5: Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (4.15). can observe δ-singularities located at the origin and two axes. Example 4. We consider the following initial condition 1 (cos θ, sin θ) r < 0.3, ρ(x, y, 0) = , (u, v)(x, y, 0) = (4.16) 1 1 (− 2 cos θ, − 2 sin θ) r > 0.3, 100 p x2 + y 2 and θ is the polar angle. Figure 4.6 shows the numerical density where r = (contour plot) and velocity field at t = 0.5. From the figure, we can observe δ-shocks located on a circle and vacuum inside. Example 5. We consider the following initial condition (0.3, 0.4) (−0.4, 0.3) ρ(x, y, 0) = 0.5, (u, v)(x, y, 0) = (−0.3, −0.4) (0.4, −0.3) 26
x > 0, y x < 0, y x < 0, y x > 0, y
> 0, > 0, < 0, < 0.
(4.17)
0.4
y
0.2
0
-0.2
-0.4 -0.4
-0.2
0
0.2
x
0.4
Figure 4.6: Numerical density (left) and velocity field (right) at t = 0.5 for initial condition (4.16).
1
0.8
y
0.6
0.4
0.2
0
0
0.2
0.4
x
0.6
0.8
1
Figure 4.7: Numerical density (left) and velocity field (right) at t = 0.4 with N = 50 for initial condition (4.17).
27
Figure 4.7 shows the numerical density (contour plot) and velocity field with N = 50 at t = 0.4. From the figure, we can observe that the numerical solution approximates the vacuum quite well.
5
Concluding remarks
In this paper, we apply discontinuous Galerkin (DG) method to solve hyperbolic conservation laws involving δ-singularities. We study Krause’s consensus models and pressureless Euler equations to demonstrate the stability and high resolution of the DG approximations. Moreover, numerical experiments show that the scheme is also good for approximations in the presence of vacuum. In future work we will extend DG methods to other equations involving δ-singularities in wider areas of applications.
References [1] C. Berthon, M. Breuss and M.-O. Titeux, A relaxation scheme for the approximation of the pressureless Euler equations, Numerical Methods for Partial Differential Equations, 22, (2006), 484-505. [2] F. Bouchut, On zero pressure gas dynamics, Advances in Kinetic Theory and Computing, Would Scientific, River Edge, NJ, (1994), 171-190. [3] F. Bouchut and F. James Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness, Communications in Partial Differential Equations, 24, (1999), 2173-2189. [4] F. Bouchut, S. Jin and X. Li, Numerical approximations of pressureless and isothermal gas dynamics, SIAM Journal on Numerical Analysis, 41, (2003), 135-158. [5] L. Boudin and J. Mathiaud, A numerical scheme for the one-dimensional pressureless gases system, Numerical Methods for Partial Differential Equations, 28, (2006), 17291746. 28
[6] Y. Brenier and E. Grenier, Sticky particles and scalar conservation laws, SIAM Journal on Numerical Analysis, 35, (1998), 2317-2328. [7] C. Canuto, F. Fagnani and P. Tilli, An Eulerian approach to the analysis of Krause’s consensus models, SIAM Journal on Control and Optimization, 50, (2012), 243-265. [8] G.-Q. Chen and H. Liu, Formation of δ-shocks and vacuum states in the vanishing pressure limit of solutions to the Euler equations for isentropic fluids, SIAM Journal on Mathematical Analysis, 34, (2003), 925-938. [9] J. Cheng and C.-W. Shu, A cell-centered Lagrangian scheme with the preservation of symmetry and conservation properties for compressible fluid flows in two-dimensional cylindrical geometry, Journal of Computational Physics, 229, (2010), 7191-7206. [10] J. Cheng and C.-W. Shu, Improvement on spherical symmetry in two-dimensional cylindrical coordinates for a class of control volume Lagrangian schemes, Communications in Computational Physics, 11, (2012), 1144-1168. [11] A. Chertock, A. Kurganov and Y. Rykov, A new sticky particle method for pressureless gas dynamics, SIAM Journal on Numerical Analysis, 45, (2007), 2408-2441. [12] 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. [13] 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. [14] 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. 29
[15] 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. [16] W. E., Yu. G. Rykov and Ya. G. Sinai, Generalized variational principles, global weak solutions and behavior with random initial data for systems of conservation laws arising in adhesion particle dynamics, Communications in Mathematical Physics, 177 (1996), 349-380. [17] S. Gottlieb, C.-W. Shu and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Review, 43 (2001), 89-112. [18] C. Johnson and J. Pitk¨aranta, An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation, Mathematics of Computation, 46 (1986), 1-26. [19] 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. [20] Yu. G. Rykov, Propagation of shock wave type singularities in equations of twodimensional zero-pressure gas dynamics, Mathematical Notes, 66 (1999), 628-635. [21] Yu. G. Rykov, On the nonhamiltonian character of shocks in 2-D pressureless gas, Bollettino della Unione Matematica Italiana. Serie VIII. Sezione B. Articoli di Ricerca Matematica, 5 (2002), 55-78. [22] C.-W. Shu, Total-variation-diminishing time discretizations, SIAM Journal on Scientific and Statistical Computing, 9 (1988), 1073-1084. [23] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shockcapturing schemes, Journal of Computational Physics, 77 (1988), 439-471.
30
[24] Y. Yang and C.-W. Shu, Discontinuous Galerkin method for hyperbolic equations involving δ-singularities: negative-order norm error estimates and applications, Numerische Mathematik, to appear. [25] 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. [26] 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. [27] X. Zhang and C.-W. Shu, Maximum-principle-satisfying and positivity-preserving high order schemes for conservation laws: Survey and new developments, Proceedings of the Royal Society A, 467 (2011), 2752-2776.
31
