    # Superconvergence of local discontinuous Galerkin methods with generalized alternating fluxes for 1D linear convection-diffusion equations

This paper investigates superconvergence properties of the local discontinuous Galerkin methods with generalized alternating fluxes for one-dimensional linear convection-diffusion equations. By the technique of constructing some special correction functions, we prove the (2k+1)th order superconvergence for the cell averages, and the numerical traces in the discrete L^2 norm. In addition, superconvergence of order k+2 and k+1 are obtained for the error and its derivative at generalized Radau points. All theoretical findings are confirmed by numerical experiments.

## Authors

##### 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

In this paper, we consider the local discontinuous Galerkin (LDG) methods for one-dimensional linear convection-diffusion equations

 ut+ux−uxx=0,       (x,t)∈[0,2π]×(0,T], (1.1a) u(x,0)=u0(x),           x∈R, (1.1b)

where is sufficiently smooth. We will consider the periodic boundary condition , the mixed boundary condition and the Dirichlet boundary condition . We study the superconvergence property concerning Radau points, cell averages, supercloseness of the LDG method with generalized alternating numerical fluxes, including the case for which the parameters involved in the numerical fluxes for the prime variable regarding the convection part and the diffusion part are independently chosen for solving (1.1).

As an extension of discontinuous Galerkin (DG) method for solving first order hyperbolic equations, the LDG method was proposed by Cockburn and Shu  in the framework of solving second-order convection-diffusion equations. The idea of the LDG methods is to rewrite the original equation with high spatial derivatives as a first order system so that the DG method can be applied. Remark that in addition to stability issue, the local solvability of auxiliary variables introduced should also be guaranteed when choosing numerical fluxes.

Being a deeper insight of DG methods, superconvergence has been investigated basically measured in the discrete norm for Radau points as well as cell averages, in the norm for the error between the numerical solution and a particular projection of the exact solution (supercloseness), and in weak negative-order norm for enhancing accuracy. For example, by virtue of the duality argument in combination with the standard optimal a priori

error estimates, Cockburn et al.

 proved that the post-processed error is of order superconvergence in the norm for linear hyperbolic systems and Ji et al.  demonstrated that the smoothness-increasing accuracy-conserving (SIAC) filter can be extended to the multidimensional linear convection-diffusion equation in order to obtain th order superconvergence, where or . Here and blow, denotes the polynomial degree of the discontinuous finite element space. Later, to efficiently compute multi-dimensional problems, the line filter, namely the one-dimensional kernel is designed via rotation in , and a rigorous proof of the post-processed errors is also given. For arbitrary non-uniform regular meshes, by establishing the relation of the numerical solution and auxiliary variable as well as its time derivative, superconvergence of order is proved for linear convection-diffusion equations . For supercloseness results concerning high order equations, see, for example, [20, 24]. Note that aforementioned supercloseness results are not sharp. In view of this, Yang and Shu  adopted the dual argument to study the sharp superconvergence of the LDG method for one-dimensional linear parabolic equations, and improved superconvergence results of order were obtained in terms of supercloseness and Radau points.

Recently, motivated by the successful applications of correction functions to finite element methods and finite volume methods for elliptic equations , Cao et al. [8, 7, 6, 5] studied superconvergence properties of DG and LDG methods for linear hyperbolic and parabolic equations. Specifically, they offered a novel proof to derive the th or th order superconvergence rate for the cell average and numerical fluxes, which will lead to the sharp

th order superconvergence for supercloseness as well as the function errors at downwind-biased points. Note that these superconvergent results are based on a supercloseness property of the DG solution to an interpolation function consisting of the difference between a standard Gauss–Radau projection of the exact solution and a carefully designed correction function. It is worth pointing out that a suitable correction is introduced to balance the difference between the projection error for the inner product term and for the DG operator term, and in standard optimal error estimates when a Gauss–Radau projection is used, the projection error involved in the DG operator term is exactly zero. This indicates that the standard Gauss–Radau projection is not the best choice for superconvergence analysis. The superconvergence of the direct DG (DDG) method for the one-dimensional linear convection-diffusion equation was studied in

. We would like to remark that all the work mentioned above are focused on purely upwind and alternating numerical fluxes.

In order to obtain flexible numerical dissipation with potential applications to nonlinear systems, the upwind-biased flux was proposed in , which is a linear combination of the numerical solution from both sides of interfaces. Stability and optimal error estimates were obtained by constructing and analyzing some suitable global projections with emphasize on analysis of some circulant matrices. Note that the design of global projections was similar to those in the work for the Burgers–Poisson equation . Moreover, Cheng et al.  studied the LDG methods for the linear convection-diffusion equations when the generalized alternating fluxes were used, and they obtained the optimal norm error estimate in a unified setting, especially when numerical fluxes with different weights are considered. In , Cao et al. investigated the superconvergence of DG methods based on upwind-biased fluxes for one-dimensional linear hyperbolic equations. More recently, Frean and Ryan  proved that the use of SIAC filters was still able to extract the superconvergence information and obtain a globally smooth and superconvergent solution of order for linear hyperbolic equations based on upwind-biased fluxes. Moreover, the -fluxes, which were introduced as linear combinations of average and jumps of the solution as well as the auxiliary variables at cell interfaces, has been a hot research topic in recent years [1, 19, 12].

In current paper, we aim at analyzing the superconvergence properties of LDG methods using generalized alternating numerical fluxes for the convection-diffusion equations. The contribution of this paper is to consider the more flexible generalized alternating fluxes. The critical step in deriving superconvergence is to construct special interpolation functions for both variables (the exact solution and the auxiliary variable ) with the aid of some suitable correction functions, essentially following . Taking into account the stability result, we use special projections to eliminate or control the troublesome jump terms involving projection errors; see, e.g. Lemma 3.1. To be more precise, we will establish the superconvergence between the LDG solution and special interpolation functions as well as , where and are correction functions to be specified later, with the main technicality being the construction and analysis of some suitable projections tailored to the very choice of the numerical fluxes. By a rigorous mathematical proof, we prove a superconvergence rate of for the errors of numerical traces and for the cell averages, and for the DG error at generalized Radau points.

The paper is organized as follows. In section 2, we present the LDG method with generalized alternating fluxes. In section 3, we construct special functions to correct the error between the LDG solution and the standard Gauss–Radau projections of the exact solution. Section 4 is the main body of the paper, in which we show and prove some superconvergence phenomena for cell averages and generalized Radau points for periodic boundary conditions. Other boundary cases including the mixed boundary condition and Dirichlet boundary condition will be considered in section 5, and the choice of numerical initial discretization is also given. In section 6, we present some numerical experiments that confirm the sharpness of our theoretical results. We will end in section 7 with concluding remarks and some possible future work.

## 2 The LDG scheme

In this section, we present the LDG scheme with generalized alternating fluxes for the linear convection-diffusion equation (1.1). As usual, we divide the computational domain into cells

 0=x12

For any positive integer , we define and denote

as the cell centers and cells, respectively. Let be the length of the cell for , and . We assume the partition is quasi-uniform in the sense that there exists a constant independent of such that . Define the finite element space

 Vkh={v∈L2(Ω):v|Ij∈Pk(Ij), ∀j∈ZN},

where is the space of polynomials on of degree at most . We use

 ¯uj+12=12(u+j+12+u−j+12),       [u]j+12=u+j+12−u−j+12

to denote the mean and jump of the function at each element boundary point , and the weighted average is denoted by

 u(θ)j+12=θu−j+12+~θu+j+12,  ~θ=1−θ,

where and are the traces from the right and left cells, respectively.

Throughout this paper, we employ to denote the standard Sobolev space on equipped with the norm with . For simplicity, we set with or . The subscript will be omitted when , and can be written as when .

In order to construct the LDG scheme, we first introduce an auxiliary variable , then the problem (1.1) can be written into a first order system

 ut+(u−q)x=0,   q−ux=0, (2.1)

where is the physical flux and is the so-called prime variable. The LDG scheme is thus to find such that for all test functions

 (uht,v)j−(uh−qh,vx)j+(~uh−^qh)v−|j+12−(~uh−^qh)v+|j−12=0, (2.2a) (qh,ψ)j+(uh,ψx)j−^uhψ−|j+12+^uhψ+|j−12=0. (2.2b)

Here , and are numerical fluxes. We use the generalized alternating numerical fluxes related to arbitrary parameters and as in . That is,

 (~uh−^qh,^uh)=(u(λ)h−q(~θ)h,u(θ)h). (2.3)

Remark that the parameters in the numerical flux regarding the convection part and diffusion part can be chosen independently, and to ensure stability the weight should satisfy .

For simplicity, we introduce the notation pertaining to the DG operator

 H1(u,q;v)=N∑j=1H1j(u,q;v),      H2(u;ψ)=N∑j=1H2j(u;ψ),

where

 H1j(u,q;v) =(q−u,vx)j−(^q−~u)v−|j+12+(^q−~u)v+|j−12, H2j(u;ψ) =(u,ψx)j−^uψ−|j+12+^uψ+|j−12.

Thus, by Galerkin orthogonality, the cell error equation can be written as

 (eut,v)j+(eq,ψ)j+H1j(eu,eq;v)+H2j(eu;ψ)=0,   ∀v,ψ∈Vkh, (2.4)

where .

For optimal error estimates of the LDG scheme using the generalized numerical fluxes (2.3) solving convection-diffusion equations with periodic boundary conditions, a globally defined projection together with is usually needed. For , the generalized Gauss–Radau projection is defined as the element of that satisfies

 ∫Ij(Pθz−z)vhdx=0,     ∀vh∈Pk−1(Ij), (2.5a) (Pθz)(θ)j+12=(z(θ))j+12,      ∀j∈ZN. (2.5b)

It has been shown in [22, 23] that the projection is well defined for , and for some restrictions on the mesh as well as polynomial degree are needed to guarantee existence and optimal approximation property of the projection . Note that when the parameter is taken as or , the projection reduces to the standard local Gauss–Radau projection or as defined in . Besides, the projection satisfies the following optimal approximation property [22, 23]

 ∥z−Pθz∥Ij+h12∥z−Pθz∥∞,Ij≤Chk+32∥z∥k+1,∞, (2.6)

where is independent of and .

To obtain the superconvergence results, the following lemma is useful in describing correction functions.  Suppose is an circulant matrix with the first row and the last row , where

. Then, for any vectors

satisfying , there holds

 |xj|≲∥b∥∞,∀j∈ZN.

## 3 Correction functions

In what follows, we will present the construction of correction functions. The cases for the weights of the prime variable in (2.3) being the same or different are discussed in the following two subsections.

### 3.1 The case with λ=θ in (2.3)

When in (2.3), to construct special interpolation functions by modifying generalized Gauss–Radau projections with correction functions so that they are superclose to the LDG solution , we start by denoting as the standard Legendre polynomial of degree on the interval , and assume that the function has the following Legendre expansion. That is, on each ,

 v(x,t)=∞∑m=0vj,m(t)Lj,m(x),      vj,m=2m+1hj(v,Lj,m)j.

By the definition of in (2.5a), we can rewrite into the following form

 Pθv=k∑m=0vj,m(t)Lj,m(x)+¯vj,k(t)Lj,k(x),

where can be determined by with

 v−Pθv=−¯vj,k(t)Lj,k(x)+∞∑m=k+1vj,m(t)Lj,m(x). (3.1)

It follows from the orthogonality of Legendre polynomials and (2.6) that,

 |¯vj,k|≲2k+1hj|(v−Pθv,Lj,k)j|≲hk+1∥v∥k+1,∞.

Following , to balance projection errors for the inner product term and the DG operator term, we define an integral operator by

 D−1xu(x)=1¯hj∫xxj−12u(τ)dτ,     τ∈Ij,

where . Obviously, . Moreover, by the properties of Legendre polynomials, we have

 D−1xLj,k(x)=12k+1(Lj,k+1−Lj,k−1)(x). (3.2)

To clearly see how to cancel terms involving projection errors with the goal of obtaining superconvergence, we split the error into two parts:

 eu=u−uh=u−uℓI+uℓI−uh≜ηu+ξu, eq=q−qh=q−qℓI+qℓI−qh≜ηq+ξq.

Then error equation (2.4) becomes

 (ξut,v)j+(ξq,ψ)j+H1j(ξu,ξq;v)+H2j(ξu;ψ) = −(ηut,v)j−(ηq,ψ)j−H1j(ηu,ηq;v)−H2j(ηu;ψ).

For periodic boundary conditions, by choosing and summing over all , we have

 12ddt∥ξu∥2+∥ξq∥2+(λ−12)N∑j=1[ξu]2j+12 = −(ηut,ξu)−(ηq,ξq)−H1(ηu,ηq;ξu)−H2(ηu;ξq). (3.3)

From the equation (3.3), we can see that in order to obtain the supercloseness properties between the numerical solution and interpolation function , we need to obtain a sharp superconvergent bound for the right-hand term, essentially using the switch of the time derivative and spatial derivative through the integral operator in combination with integration by parts; see Lemma 3.1 below. Next, we show how to construct interpolation functions and estimate the right-hand side of (3.3).

To construct the interpolation functions , we define a series of functions as follows:

 (wu,i−¯hjD−1xwq,i−1,v)j=0,              (w(θ)u,i)j+12=0, (3.4a) (wq,i−wu,i−¯hjD−1x∂twu,i−1,v)j=0,   (w(~θ)q,i)j+12=0, (3.4b)

where , and

 wu,0=u−Pθu,      wq,0=q−P~θq.

The functions defined in (3.4) have the following properties

 ∥∂twu,i∥∞≲hk+i+1∥u∥k+i+3,∞,  ∥wq,i∥∞≲hk+i+1∥u∥k+i+2,∞, (3.5a) (wu,i,v)j=0, (wq,i,v)j=0, ∀v∈Pk−i−1(Ij). (3.5b)
###### Proof.

The proof of this lemma is based on deriving the following expression of and in each element , which can be obtained by induction. It reads,

 wu,i∣∣Ij=k∑m=k−iβji,mLj,m(x),    wq,i∣∣Ij=k∑m=k−iγji,mLj,m(x),   i∈Zk. (3.6)

Step 1: When , by taking with in (3.4a) and using (3.2) together with the orthogonality property of Legendre polynomials, we obtain

 (wu,1−¯hjD−1xwq,0,v)=(βj1,k−1Lj,k−1−¯qj,k2k+1¯hjLj,k−1,v)=0.

Obviously, , where is the coefficient of the Legendre expansion for ; see (3.1) with replaced by and replaced by . Using the fact that we have

 θβj1,k+(−1)k(1−θ)βj+11,k=(−1)k(1−θ)βj+11,k−1−θβj1,k−1. (3.7)

Then the linear system (3.7) can be written in the matrix-vector form

 Aβ1,k=b,

where is an circulant matrix, and

 β1,k=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝β11,kβ21,k⋮βN1,k⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,     b=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−θβ11,k−1+(−1)k(1−θ)β21,k−1−θβ21,k−1+(−1)k(1−θ)β31,k−1⋮−θβN1,k−1+(−1)k(1−θ)β11,k−1⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

It is easy to compute the determinant of the in the form

 |A|=θN(1−pN),    p=(−1)k(θ−1)θ,

and for the matrix is always invertible. Therefore, the linear system (3.7) has the unique solution. Moreover, by Lemma 2, we have

 |βj1,k|≲max1⩽ℓ⩽N|bℓ|≲hk+2∥u∥k+2,∞,  ∀j∈ZN.

Thus,

 ∥∂twu,1∥∞,Ij=∥∂t(βj1,k−1Lj,k−1+βj1,kLj,k)∥∞,Ij≲hj|∂t¯qj,k|≲hk+2∥u∥k+4,∞.

Similarly, when choosing in (3.4b), we obtain

 wq,1∣∣Ij=k∑m=k−1γj1,mLj,m,

where

 γj1,k−1=βj1,k−1+∂t¯uj,k2k+1¯hj,

and is the solution of the following linear system

 ~Aγ1,k=~b,

with being an circulant matrix, and

 γ1,k=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝γ11,kγ21,k⋮γN1,k⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,     ~b=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−~θγ11,k−1+(−1)k(1−~θ)γ21,k−1−~θγ21,k−1+(−1)k(1−~θ)γ31,k−1⋮−~θγN1,k−1+(−1)k(1−~θ)γ11,k−1⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

Consequently, the estimate of in (3.5a) follows by using Lemma 2 and optimal approximation property (2.6). Moreover, (3.5b) is a trivial consequence of the expression (3.6) when the orthogonality property of Legendre polynomials is taken into account.

Step 2: Suppose that (3.5a) and (3.6) are valid for all and we want to prove it still holds for . From equation (3.4a), we can get

 (wu,i+1−¯hjD−1x(k∑m=k−iγji,mLj,m),v)j=0,  ∀v∈Pk−1(Ij).

In order to get the superconvergent bounds of , we need to find out the expression of coefficient . After a direct calculation, there hold

 βji+1,k−i−1=−γji,k−i¯hj2(k−i)+1,  βji+1,k−i=−γji,k−i+1¯hj2(k−i)+3, βji+1,m=¯hj(γji,m−12m−1−γji,m+12m+3),  m=k−i+1,⋯,k−1.

Moreover, by the fact that , we get

 θ(βji+1,k−i−1+⋯+βji+1,k)+(1−θ)(−1)k−i−1βj+1i+1,k−i−1+⋯+(1−θ)(−1)kβj+1i+1,k=0.

Again, we can write the above linear system into the matrix-vector form

 Aβi+1,k=c,

and when , we arrive at the unique existence of the system. Consequently, it follows from Lemma 2 that

 ∥∂twu,i+1∥∞,Ij ≲k∑m=k−i−1|∂tβji+1,m|≲hk∑m=k−i|∂tγji,m| ≲h∥∂twq,i∥∞≲hk+i+2∥∂tq∥k+i+1,∞≲hk+i+2∥u∥k+i+4,∞.

Analogously, the other estimate of (3.5a) can be obtained, and the orthogonality property in (3.5b) is a trivial consequence of expression of and in (3.6) with replaced by . This finishes the proof of Lemma 3.1. ∎

We are now ready to define the correction functions as follows. For any positive integer , we define in each element

 Wℓu=ℓ∑i=1wu,i,     Wℓq=ℓ∑i=1wq,i, (3.8)

and the special interpolation functions are

 uℓI=Pθu−Wℓu,     qℓI=P~θq−Wℓq. (3.9)

Suppose is the solution of (1.1), and are defined by (3.9), then for , we have

 |((u−uℓI)t,v)j−(Wℓu,vx)j+(Wℓq,vx)j|≲hk+ℓ+1∥u∥k+ℓ+3,∞∥v∥1,Ij, (3.10a) |(q−qℓI,ψ)j+(Wℓu,ψx)j|≲hk+ℓ+1∥u∥k+ℓ+2,∞∥ψ∥1,Ij. (3.10b)
###### Proof.

By the orthogonality property of , we have

 D−1xwu,i(x−j+12)=1¯hj(wu,i,1)j=0=D−1xwu,i(x+j−12), D−1xwq,i(x−j+12)=1¯hj(wq,i,1)j=0=D−1xwq,i(x+j−12).

It follows from integration by parts that

 (∂twu,i,v)j=−¯hj(D−1x∂twu,i,vx)j=−(wq,i+1−wu,i+1,vx), v∈Vkh, i∈Zk−1, (wq,i,v)j=−¯hj(D−1xwq,i,vx)j=−(wu,i+1,vx),                  v∈Vkh, i∈Zk−1.

Then

 ((u−uℓI)t,v)j−(Wℓu,vx)j+(Wℓq,vx)j = ((u−Pθu)t,v)j+ℓ∑i=1[(∂twu,i,v)j+(wq,i−wu,i,vx)j] = (∂twu,ℓ,v)j.

Similarly, there holds

 (q−qℓI,ψ)j+(Wℓu,ψx)j=(wq,ℓ,ψ)j.

By (3.5a), we can get the desired result (3.10). ∎

### 3.2 The case with λ≠θ in (2.3)

When parameters and in (2.3) pertaining to convection and diffusion terms are chosen differently, a pair of suitable interpolation functions in possession of supercloseness property should be constructed, which are based on a combination of modified projections and new correction functions. To do that, let us first recall a new modified projection. That is,

 Πh(u,q)=(Pθu,P∗~θq),

in which has been given in (2.5a), and depends on both and satisfying