# Sub-optimal convergence of discontinuous Galerkin methods with central fluxes for even degree polynomial approximations

In this paper, we theoretically and numerically verify that the discontinuous Galerkin (DG) methods with central fluxes on non-uniform meshes have sub-optimal convergence properties when measured in the L^2-norm for even degree polynomial approximations. On uniform meshes, the optimal error estimates are provided for arbitrary number of cells in one and multi-dimensions, improving previous results. The theoretical findings are found to be sharp and consistent with numerical results.

There are no comments yet.

## Authors

• 141 publications
• 16 publications
• 3 publications
12/21/2020

### A p-robust polygonal discontinuous Galerkin method with minus one stabilization

We introduce a new stabilization for discontinuous Galerkin methods for ...
08/19/2019

### Discontinuous Galerkin approximations in computational mechanics: hybridization, exact geometry and degree adaptivity

Discontinuous Galerkin (DG) discretizations with exact representation of...
06/07/2018

### Optimal Design of Process Flexibility for General Production Systems

Process flexibility is widely adopted as an effective strategy for respo...
01/09/2021

### On the suboptimality of the p-version discontinuous Galerkin methods for first order hyperbolic problems

We address the issue of the suboptimality in the p-version discontinuous...
12/30/2021

### DPG methods for a fourth-order div problem

We study a fourth-order div problem and its approximation by the discont...
10/16/2019

### Convergence in the incompressible limit of new discontinuous Galerkin methods with general quadrilateral and hexahedral elements

Standard low-order finite elements, which perform well for problems invo...
10/14/2020

### An adaptive multigrid solver for DPG methods with applications in linear acoustics and electromagnetics

We propose an adaptive multigrid preconditioning technology for solving ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

A fundamental form of energy transmission is wave propagation, which arises in many fields of science, engineering and industry, such as petroleum engineering, geoscience, telecommunication, and the defense industry (see [7, 11]). It is important for these applications to study efficient and accurate numerical methods to solve wave propagation problems. Experience reveals that energy-conserving numerical methods, which conserve the discrete approximation of energy, are favorable, because they are able to maintain the phase and shape of the waves more accurately, especially for long-time simulation.

Various numerical approximations of wave problems modeled by linear hyperbolic systems can be found in the literature. Here, we will focus on the classical Runge-Kutta DG method of Cockburn and Shu [6]. There are several approaches to obtain an optimal, energy conserving DG method. Chung and Engquist [4] presented an optimal, energy conserving DG method for the acoustic wave equation on staggered grids. Chou et al. [3] proposed an optimal energy conserving DG using alternating fluxes for the second order wave equation. More recently, Fu and Shu [8] developed an optimal energy conserving DG method by introducing an auxiliary zero function.

As is well known, the simplest energy conserving DG method for hyperbolic equations is the one using central fluxes. However, it has sub-optimal convergence of order when measured in the -norm when piece-wise polynomials of an odd degree are used; see, e.g. [14]. When is even, we usually observe higher convergence rates than th order for a general regular non-uniform meshes, such as random perturbation over an uniform mesh, see section 4. In fact, many papers have mentioned that the optimal convergence rates can be observed when even degree polynomials are used; see for example [1, 2, 14]. In this paper, we provide a counter example to show that the scheme only has sub-optimal error accuracy of order for a regular non-uniform mesh, when is even. We refer to the work of Guzmán and Rivière [10] in which they constructed a special mesh sequence to produce the sup-optimal accuracy for the non-symmetric DG methods for elliptic problems when

is odd. For uniform meshes, the classical DG scheme with the central flux does have the optimal convergence rate

, observed in the numerical experiments and proved theoretically under the condition that the number of cells in the mesh is [1, 14]. In this paper, we provide a new proof which is available for arbitrary number of cells and dimensions. We have used the shifting technique [12, 13] to construct the special local projection to obtain the optimal error estimate on uniform meshes. We also numerically find the superconvergence phenomenon for the cell averages and numerical fluxes.

The outline of the paper is as follows. In section 2, we review the DG scheme for hyperbolic equations with central fluxes and give the error estimates for the semi-discrete version in one dimension. We extend our analysis to multi-dimensions in section 3. In section 4, we give numerical examples to show the sub-optimal convergence for non-uniform meshes and optimal convergence for uniform meshes in both one and two-dimensional cases. Finally, we give concluding remarks in section 5. Some of the technical proof of the lemmas and propositions is included in the Appendix A.

## 2 One dimensional problems

We consider the following one dimensional linear hyperbolic equation

 {ut+ux=0, x∈[0,1], t≥0u(x,0)=u0(x), x∈[0,1], (2.1)

with periodic boundary condition. We first introduce the usual notations of the DG method. For a given interval and the index set , the usual DG mesh is defined as:

 0=x12

We denote

 Ij=(xj−12,xj+12),xj=12(xj−12+xj+12),hj=xj+12−xj−12,j∈ZN. (2.3)

We also assume the mesh is regular, i.e., the ratio between the maximum and minimum mesh sizes shall stay bounded during mesh refinements. That means there exists a positive constant , such that,

 1σh≤hj≤σh,h=1N,∀j∈ZN. (2.4)

We define the approximation space as

 Vkh={vh:(vh)|Ij∈Pk(Ij),j=1,…,N}. (2.5)

Here denotes the set of all polynomials of degree at most on . We first introduce some standard Sobolev space notations. For any integer , denote the standard Sobolev spaces on the sub-domain equipped with the norm and the semi-norm . If , we set , and and we omit the index , when .

The semi-discrete DG scheme is to seek such that for all ,

 ((uh)t,vh)j+aj(uh,vh)=0,∀j∈ZN, (2.6)

where

 aj(uh,vh)=−(uh,(vh)x)j+^uhv−h|j+12−^uhv+h|j−12, (2.7)

where , and denote the left and right limits of at the point , respectively, and are the numerical fluxes. Here, we consider the central flux,

 ^uh={uh}=12(u−h+u+h). (2.8)

The initial datum is obtained by the standard projection,

 (u0−Pu0,vh)j=0,∀vh∈Pk(Ij). (2.9)

Thus, we have,

 ∥u0−uh(⋅,0)∥≲hk+1∥u∥k+1 (2.10)

Here and below, an unmarked norm denotes the norm, and denotes that can be bounded by multiplied by a constant independent of the mesh size . As mentioned earlier, we have the following energy-conserving results [14].

###### Theorem 2.1.

Suppose is the solution of DG scheme (2.6), then it satisfies

 ddt∥uh∥2=0. (2.11)

Next we consider the error estimate, first we recall the following basic facts [5]. For any function ,

 (i) ∥(wh)x∥≲h−1∥wh∥, (ii) ∥wh∥Γh≲h−12∥wh∥. (2.12)

where denotes the set of boundary points of all elements , and the norm . In order to obtain the optimal error estimate for the case of uniform meshes, we need to use the shifting technique [12, 13] to construct a special projection , which is defined as follows. For any given function and each ,

 ∫IjP⋆hw(x)dx=∫Ijw(x)dx (2.13) ˜Ph(P⋆hw;v)j=˜Ph(w;v)j∀v∈Pk(Ij) (2.14)

where is defined as

 ˜Ph(w;v)j=−(w,vx)j+w(x−j+12)+w(x+j−12)2(v(x−j+12)−v(x+j−12)) (2.15)

Note that the projection is a local projection, so we only consider the projection defined on the reference interval . We have the following lemma to establish the fact that the projection is well defined.

###### Lemma 2.1.

When is even, the projection defined by (2.13) on the interval exists and is unique for any function , and the projection is bounded in the norm, i.e.,

 ∥P⋆hw∥∞≤C(k)∥w∥∞ (2.16)

where is a constant that only depends on but is independent of .

###### Proof.

We provide the proof of this lemma in the appendix; see section A.1.

###### Remark 2.1.

The projection is only well defined when is even. In fact, when is odd, for example , we can take , which satisfies

 ∫1−1wI(x)dx=0, (2.17) ˜Ph(wI;v)=−∫1−1wI(x)vxdx+wI(1)+wI(−1)2(v(1)−v(−1))=0,∀v∈P1([−1,1]). (2.18)

It means that there exists a nonzero function , where . This implies that is not unique.

###### Remark 2.2.

In fact, the projection has an equivalent definition as follows,

 ∫xj+12xj−12P⋆hwvdx=∫xj+12xj−12wvdx,∀v∈Pk−1(Ij), (2.19) 12(P⋆hw(x−j+12)+P⋆hw(x+j−12))=12(w(xj+12)+w(xj−12)) (2.20)

As a direct corollary of lemma 2.1, the standard approximation theory [5] implies, for a smooth function ,

 ∥P⋆hw(x)−w(x)∥+h12∥P⋆hw(x)−w(x)∥Γh≲hk+1∥w∥k+1. (2.21)

We also have the following properties of the projection ,

###### Lemma 2.2.

Suppose that . Let ,. If , then we have following relationship:

 (x−h)k+1−uj−1(x−h)=xk+1−uj(x)=(x+h)k+1−uj+1(x+h),∀x∈Ij. (2.22)

where means that the projection of is defined on the subinterval , and , refer to the projection of on the element and respectively, since implies and .

###### Proof.

The proof of this lemma is by the same arguments as in [12, 13], so we omit it here.

By this lemma, we also have the following superconvergence results.

###### Proposition 2.1.

For any given the index , suppose that is a th degree polynomial function in . If , we have

 aj(P⋆hu,vh)=aj(u,vh)∀vh∈Pk(Ij), (2.23)

where is defined by (2.7)

Then we can state the main theorem of this paper.

###### Theorem 2.2.

Suppose is the numerical solution of the DG scheme (2.6) for equation (2.1) with a smooth initial condition , and is the exact solution of (2.1), then the approximation satisfies the following error estimate:

 ∥u(⋅,T)−uh(⋅,T)∥≲hk, (2.24)

where is the degree of the piecewise polynomials in the finite element spaces . Furthermore, when is even and the mesh is uniform, we have the optimal error estimate:

 ∥u(⋅,T)−uh(⋅,T)∥≲hk+1. (2.25)
###### Proof.

Obviously, the exact solution of (2.1) also satisfies

 (ut,vh)+aj(u,vh)=0,∀vh∈Vh. (2.26)

Subtracting (2.6) from (2.26), we obtain the error equation

 ((u−uh)t,vh)j+aj(u−uh,vh)=0,∀vh∈Vh. (2.27)

We denote

 ξ=uh−P⋆u;η=u−P⋆u, (2.28)

where is some projection. From the error equation (2.27), and taking , we have

 (ξt,ξ)j+aj(ξ,ξ)=(ηt,ξ)j+aj(η,ξ) (2.29)

For the nonuniform mesh case, the sub-optimal error estimate can be easily obtained by using the standard projection . We take as the standard projection , then we have,

 ∥u−Pu∥+h12∥u−Pu∥Γh≲hk+1∥u∥k+1 (2.30)

For the left-hand side of (2.29), we can use the stability result (2.1) to obtain

 12ddt∥ξ∥2 =−N∑j=1{η}j+12[ξ]j+12 ≲hk∥ξ∥∥u∥k+1 (2.31)

where the last inequality is from (2.30) and (ii) of (2). Thus, by using Gronwall’s inequality and (2.10), we have,

 ∥ξ∥≲hk∥u∥k+1 (2.32)

The triangle inequality implies our designed results for the general non-uniform mesh case.

For the case of uniform meshes, when is even, we take as which is defined in (2.13). Let be the Taylor expansion polynomial of order of over , i.e., , . Let denote the remainder term, i.e., . Recalling the Bramble-Hilbert lemma in [5], we have

 ∥rju∥L∞(Dj)≲hk+32|u|k+2,Dj. (2.33)

Thus, using Proposition 2.1, we have

 aj(η,ξ) =aj(ujI−P⋆hujI,ξ)+aj(rju−P⋆hrju,ξ) =aj(rju−P⋆hrju,ξ)

By using the property of the projection (2.16) and (2.33), and the inverse inequality in (2) for , we have

 ∑jaj(η,ξ)≲h2k+2∥u∥k+2+∥ξ∥2 (2.34)

Therefore, form (2.29), (2.21) and the stability result (2.1), we have

 12ddt∥ξ∥2≲h2k+2∥u∥k+2+∥ξ∥2 (2.35)

This together with the approximation results (2.21) and the initial datum (2.10), implies the desired error estimate (2.25). ∎

We summarize the theoretical findings and numerical findings in Table 2.1.

From Table 2.1, we can see that our theoretical findings are sharp and consistent with the numerical results. We emphasize that when is even, in order to produce the sub-optimal accuracy, we have designed a special regular mesh sequence which is motivated by [10], see section 4.

## 3 Multi-dimensional problems

In this section, we consider the semidiscrete DG method with central fluxes for multidimensional linear hyperbolic equations. Without loss of generality, we only study the two dimensional problem; all the arguments we present in our analysis depends on the tensor product structure of the mesh and the finite element space and can be easily extended to the more general cases

. Hence, we consider the following two-dimensional problem

 {ut+ux+uy=0, (x,y,t)∈Ω×(0,T],u(x,y,0)=u0(x,y), (x,y)∈Ω. (3.1)

again with periodic boundary conditions. Without loss of generality, we assume . We use the regular Cartesian mesh, , , . We denote , and . Let , where denotes the space of tensor-product polynomials of degrees at most in each variable defined on .

The semidiscrete DG scheme with central fluxes is as follows. We seek , such that for all test functions , and all ,

 ∫Ki,j(uh)tvdxdy= ∫Ki,j(uhvx+uhvy)dxdy −∫yj+12yj−12(^uh(xi+12,y)v(x−i+12,y)−^uh(xi−12,y)v(x+i−12,y))dy −∫xi+12xi−12(~uh(x,yj+12)v(x,y−j+12)−~uh(x,yj−12)v(x,y+j−12))dx (3.2) =:bi,j(uh,v) (3.3)

where

 ^uh(xi+12,y)=uh(x+i+12,y)+uh(x−i+12,y)2;~uh(x,yj+12)=uh(x,y+j+12)+uh(x,y−j+12)2 (3.4)

For the initial data, we take , where is the projection into , and we have

 ∥u0−Pu0∥≲hk+1∥u0∥k+1. (3.5)

We also have the following energy conservative property

###### Proposition 3.1.

The numerical solution of (3.2) satisfies

 12ddt∥uh∥2=0. (3.6)

### 3.1 A priori error estimates

Let us now state our main result as a theorem, whose proof will be provided in the next subsection.

###### Theorem 3.1.

Suppose is the numerical solution of the DG scheme (3.2) for equation (3.1) with a smooth initial condition , and is the exact solution of (3.1), then the approximation satisfies the following error estimate:

 ∥u(x,y,T)−uh(x,y,T)∥≲hk, (3.7)

where is the degree of the piecewise tensor-product polynomials in the finite element spaces . Furthermore, when is even and the mesh is uniform, we have the optimal error estimate,

 ∥u(x,y,T)−uh(x,y,T)∥≲hk+1. (3.8)
###### Remark 3.1.

We note that the finite element space , where denotes the space of polynomials of degrees at most defined on , can also be taken as the approximation space. But it only has the sub-optimal accuracy of order in the numerical examples, see section 4. Thus, here we only consider the tensor product space.

By the same arguments as in the one dimensional problem, we also have the error equation

 ∫Ki,j(u−uh)tvdxdy−bi,j(u−uh,v)=0,∀v∈Wh,∀i,j. (3.9)

### 3.2 Proof of the error estimates

We divide the proof of Theorem 3.1 into several steps. First, for non-uniform meshes, the proof of the sub-optimal error estimate is straightforward. We just need to use the standard

projection and follow the standard error estimates of DG methods which is the same as in the one dimensional case. Thus next we only consider the uniform mesh case. In order to prove the optimal error estimate when

is even, we need to construct the special local projection . In addition, the optimal approximation properties of are derived. The superconvergence results of the special projections would be given in the subsection 3.2.2. Finally, we finish the proof of Theorem 3.1 in subsection 3.2.3.

#### 3.2.1 The special projection Π⋆h

Since our finite element space consists of piecewise polynomials, we use the tensor product technique to construct the 2D projection. We define as the following projection into . For each ,

 ∫Ki,jΠ⋆hw(x,y)v(x,y)dxdy=∫Ki,jw(x,y)v(x,y)dxdy,∀v∈Qk−1(Ki,j). (3.10a) ∫IiΠ⋆hw(x,y−j+12)+Π⋆hw(x,y+j−12)2φ(x)dx=∫Iiw(x,y−j+12)+w(x,y+j−12)2φ(x)dx,∀φ(x)∈Pk−1(Ii) (3.10b) ∫JjΠ⋆hw(x−i+12,y)+Π⋆hw(x+i−12,y)2φ(y)dy=∫Jjw(x−i+12,y)+w(x+i−12,y)2φ(y)dy,∀φ(y)∈Pk−1(Jj) (3.10c) 14(Π⋆hw(x−i+12,y−j+12)+Π⋆hw(x−i+12,y+j−12)+Π⋆hw(x+i−12,y−j+12)+Π⋆hw(x+i−12,y+j−12)) =14(w(x−i+12,y−j+12)+w(x−i+12,y+j−12)+ w(x+i−12,y−j+12)+w(x+i−12,y+j−12)) (3.10d)

Again, since the projection is local, we only consider the projection defined on the reference cell . We establish the existence and uniqueness of the projection when is even in the following lemma

###### Lemma 3.1.

When is even, the projection defined by (3.10) on the cell exists and is unique for any function , and the projection is bounded in the norm, i.e.

 ∥Π⋆hw∥∞≤C(k)∥w∥∞, (3.11)

where is a constant that only depends on but is independent of .

###### Proof.

The proof of this lemma is given in the Appendix; see section A.2.

Since the projection is a -th degree polynomial preserving local projection, standard approximation theory [5] implies, for a smooth function ,

 ∥w−Π⋆hw∥L2(Ki,j)≲hk+1∥w∥k+1,Ki,j (3.12)

For the two dimensional space, for any , the following inequalities hold,

 ∥∂xωh∥≲h−1∥ωh∥,∥ωh∥L2(∂Ki,j)≲h−1/2∥ωh∥,∥ωh∥∞≲h−1∥ωh∥, (3.13)

where is the boundary of cell .

###### Remark 3.2.

By similar arguments as in the one dimensional problem, we note that the projection is not well defined when is odd.

#### 3.2.2 Properties of the projection Π⋆h

Since the projection is local, we have the following lemma:

###### Lemma 3.2.

Assume that or . Let . If and , then , we have following relationship:

 u(x−hx,y)−ui−1,j(x−hx,y)=u(x,y)−ui,j(x,y)=u(x+hx)−ui+1,j(x+hx,y) =u(x,y+hy)−ui,j+1(x,y+hy)=u(x,y−hy)−ui,j−1(x,y−hy). (3.14)

Similar to the one dimensional case, we also have the following superconvergence result.

###### Proposition 3.2.

For a given index , suppose that is a th degree polynomial function in , where . If and , then we have

 (3.15)

where is defined by (3.3).

###### Proof.

We provide the proof of this Proposition in the Appendix; see section A.3.

#### 3.2.3 Proof of Theorem 3.1

Let

 ξ=uh−Π⋆hu;η=u−Π⋆hu. (3.16)

From (3.9), we obtain

 ∫Ki,j(ξ)tvdxdy−bi,j(ξ,v)=∫Ki,j(η)tvdxdy−bi,j(η,v),∀v∈Qk(Ki,j). (3.17)

Take , for the left hand side of (3.17), we use the stability result Proposition 3.1 to obtain

 ∑i,j∫Ki,j(ξ)tξdxdy−bi,j(ξ,ξ)=12ddt∥ξ∥2. (3.18)

For each element , we consider the Taylor expansion of around :

 u=Tu+Ru, (3.19)

where

 Tu=k+1∑l=0l∑m=01(l−m)!m!∂lu(xi,yj)∂xl−m∂ym(x−xi)l−m(y−yj)m,
 Ru=(k+2)k+2∑m=0(x−xi)k+2−m(y−yj)m(k+2−m)!m!∫10(1−s)∂k+2u(xsi,ysj)∂xk+2−m∂ym ds.

with , . Clearly, . By the linearity of the projection, and from (3.15), we then get

 bi,j(η,v) =bi,j(Tu−Π⋆hTu,v)+bi,j(Ru−Π⋆hRu,v) =bi,j(Ru−Π⋆hRu,v) (3.20)

Again recalling the Bramble-Hilbert lemma in [5] , we have

 ∥Ru∥L∞(Di,j)≤Chk+1|u|Hk+2(Di,j). (3.21)

Thus, this together with the standard approximate proposition of the projection (3.12), and the inverse inequality in (3.13) for , we have

 ∑i,jbi,j(η,ξ)≲h2k+2∥u∥2k+2+∥ξ∥2 (3.22)

From (3.18), (3.22) and (3.17), we have

 12ddt∥ξ∥2≲h2k+2∥u∥2k+2+∥ξ∥2 (3.23)

This together with the approximation results (3.12) and the initial discretization (3.5), implies the desired error estimate (3.8).

To end this section, we summarize our theoretical findings and numerical findings for the 2D problem in Table 3.1. Again our theoretical proof is sharp and consistent with the numerical results.

## 4 Numerical examples

In this section, we present some numerical examples to verify our theoretical findings. In our numerical experiments, we present the , , and errors, respectively. They are defined by

 E2= ∥u−uh∥, (4.1) EA= ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩(1N∑Nj=1(1hj∫Ij(u−uh)dx)2)12, for one% dimension(1NxNy