A Simplified Weak Galerkin Finite Element Method: Algorithm and Error Estimates

# A Simplified Weak Galerkin Finite Element Method: Algorithm and Error Estimates

Yujie Liu School of Data and Computer Science, Sun Yat-sen University, Guangzhou, 510275, China (liuyujie5@mail.sysu.edu.cn). The research of Liu was partially supported by Guangdong Provincial Natural Science Foundation (No. 2017A030310285), Shandong Provincial natural Science Foundation (No. ZR2016AB15) and Youthful Teacher Foster Plan Of Sun Yat-Sen University (No. 171gpy118),    Junping Wang Division of Mathematical Sciences, National Science Foundation, Alexandria, VA 22314 (jwang@nsf.gov). The research of Wang was supported by the NSF IR/D program, while working at National Science Foundation. However, any opinion, finding, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation.
###### Abstract

In this article a simplified weak Galerkin finite element method is developed for the Dirichlet boundary value problem of convection-diffusion-reaction equations. The simplified weak Galerkin method utilizes only the degrees of freedom on the boundary of each element and, hence, has significantly reduced computational complexity over the regular weak Galerkin finite element method. A stability and some optimal order error estimates in the and norms are established for the corresponding numerical solutions. Numerical results are presented to verify the theory error estimates and a superconvergence phenomena on rectangular partitions.

c

onvection-diffusion-reaction equations, simplified weak Galerkin, finite element methods, error estimates.

{AMS}

Primary, 65N30, 65N15; Secondary, 35J50

## 1 Introduction

This paper is concerned with the development of a simplified formulation for the weak Galerkin finite element method for second order elliptic equations. For simplicity, consider the model problem that seeks an unknown function satisfying

 (1) −∇⋅(α∇u)+β⋅∇u+cu = fin Ω (2) u = gon ∂Ω

where is a bounded polytopal domain in with boundary , is the diffusion coefficient, is the convection, and is the reaction coefficient in relevant applications. We assume that is sufficient smooth, , and is piecewise smooth with respect to a partition of the domain. For well-posedness of the problem (1)-(2), we assume , , and

 (3) c−12∇⋅β≥0,α(x)≥α0∀x∈Ω

for a constant .

The model problem (1)-(2) arises from many scientific applications such as fluid flow in porous media. Mostly importantly, this model problem has served, and still serves, the scientific computing community as a testbed in the search and design of new and efficient computational algorithms for partial differential equations. The classical Galerkin finite element method (see, e.g., [10, 28, 16]) is particularly a numerical technique originated from the study of elliptic problems closed related to (1)-(2) or its variations. In the last three decades, various finite element methods using discontinuous trial and test functions, including discontinuous Galerkin (DG) methods and weak Galerkin (WG) methods, have been developed for numerical solutions of partial differential equations. These developments were often tested over testbed problems such as (1)-(2) before they were generalized or applied to more complex problems in science and engineering. The DG method, also known as the interior penalty method in different contexts, was originated in early 70s of the last century for a numerical study of model problems such as (1)-(2); see, e.g., [3, 14, 25, 38] for early incubations and [1, 13, 17, 27] for a detailed discussion and recent developments.

The weak Galerkin finite element method is a recently developed discretization framework for partial differential equations [36, 37, 24, 34]. With new concepts referred to as weak differential operators (e.g., weak gradient, weak curl, weak Laplacian etc.) and weak continuity through the use of various stabilizers, the method allows the use of totally discontinuous functions and provides stable numerical schemes that are parameter-independent or free of locking [33]. For the convection-diffusion-reaction equation (1)-(2), the recent work in the context of weak Galerkin includes the algorithm developed and analyzed in [9], the one in [20] for singularly perturbed problems, and an earlier one in [39]. The WG finite element method has been rapidly developed and applied to several different types of problems, including second order elliptic problems, the Stokes and Navier-Stokes equations, the biharmonic and elasticity equations, div-curl systems and the Maxwell’s equations, etc. The latest development of the WG methods is the prime-dual formulation for problems that are either nonsymmetric or do not have variational forms friendly for numerical use. Details on the new developments can be found in [30] for second order elliptic equations in nondivergence form, [31] for the Fokker-Planck equation, and [32] for elliptic Cauchy problems.

The typical WG method for the model problem (1)-(2) seeks weak finite element approximations satisfying and

 (4) S(uh,v)+(α∇wuh,∇wv)+(β⋅∇wuh,v0)+(cu0,v0)=(f,v0)

for all test functions satisfying , where is an interpolation of the Dirichlet boundary data, is the discrete weak gradient operator, and is a properly selected stabilizer that gives weak continuities for the numerical solutions. The numerical solution consists of two components: the approximation on each element and the approximation on the boundary of each element. To reduce the computational complexity, some hybridized formulations have been introduced in [22, 29] for the method when applied to the diffusion equation and the biharmonic equation through the elimination of the degrees of freedom associated with the unknown function locally on each element. In the superconvergence study for WG [18] on rectangular elements, this hybridized formulation was further simplified in the description of the numerical algorithm, yielding a simplified weak Galerkin (SWG) finite element scheme for the diffusion equation. In our further investigation of the SWG to the convection-diffusion-reaction equation (1), we came to the conclusion that SWG represents a new discretization scheme that is different from the usual WG through a simple elimination of the unknown . As a result, we believe that a systematic study of the SWG for the convection-diffusion-reaction problem (1)-(2) should be conducted for its stability and convergence. This paper is in response to this observation and shall provide a mathematical theory for the stability and the convergence of the simplified weak Galerkin finite element method for the model problem (1)-(2). We believe that the result of this paper can be extended to other types of modeling equations.

The paper is organized as follows: In Section 2, we shall describe the simplified weak Galerkin finite element method for (1)-(2) on general polygonal partitions. In Section 3, we shall present a computational formula for the element stiffness matrices and the element load vectors from SWG. In Section 4, we provide a mathematical theory for the stability and well-posedness of the SWG scheme. Sections 5 and 6 are devoted to a discussion of the error estimates in a discrete and the norm for the numerical solutions. Finally, in Section 7, we present some numerical results to demonstrate the efficiency and accuracy of the SWG method.

Throughout the rest of the paper, we assume and shall use the standard notations for Sobolev spaces and norms [10, 16]. For any open set , and denote the norm and inner-product in the Sobolev space consisting of square integrable partial derivatives up to order . When or , we shall drop the corresponding subscripts in the norm and inner-product notation.

## 2 Algorithm on Polymesh

Assume that the domain is of polygonal type and is partitioned into non-overlap polygons that are shape regular. For each , denote by its diameter and by the number of edges. For each edge , denote by the midpoints and the outward normal direction of (see Fig. 1 for an illuatration). The meshsize of is defined as .

Let be a piecewise constant function defined on the boundary of , i.e.,

 vb|ei=vb,i,

with being a constant. We define the weak gradient of on by:

 (5) ∇wvb:=1|T|N∑i=1vb,i|ei|ni,

where is the length of the edge and is the area of the element . It is not hard to see that the weak gradient satisfies the following equation:

 (6) (∇wvb,ϕ)T=⟨vb,ϕ⋅n⟩∂T

for all constant vector . Here and in what follows of the paper, stands for the usual inner product in .

Denote by the space of piecewise constant functions on . The global finite element space is constructed by patching together all the local elements through single values on interior edges. The subspace of consisting of functions with vanishing boundary value is denoted as .

We use the conventional notation of for the space of polynomials of degree on . For each , we associate it with a linear extension in , denoted as , satisfying

 (7) N∑i=1(s(vb)(Mi)−vb,i)ϕ(Mi)|ei|=0,∀ϕ∈P1(T).

It is easy to see that is well defined by (7), and its computation is local and straightforward. In fact, can be viewed as an extension of from to through a least-squares fitting.

On each element , we introduce the following bilinear forms:

 (8) aT(ub,vb) := (α∇wub,∇wvb)T, (9) bT(ub,vb) := (β⋅∇wub,s(vb))T, (10) cT(ub,vb) := (cs(ub),s(vb))T.

For simplicity, we set

 (11) BT(ub,vb):=aT(ub,vb)+bT(ub,vb)+cT(ub,vb)

for . We further introduce the stabilizer

 (12) ST(ub,vb):=h−1N∑i=1(s(ub)(Mi)−ub,i)(s(vb)(Mi)−vb,i)|ei|=h−1⟨Qbs(ub)−ub,Qbs(vb)−vb⟩∂T,

where is the projection operator onto ; namely is the average of on each edge. In particular, is well-defined and takes the average of the Dirichlet data on each boundary edge.

###### SWG Algorithm 2.1

The simplified weak Galerkin (SWG) scheme for the elliptic equation (1)-(2) seeks satisfying on and

 (13) A(ub,vb)=(f,s(vb))∀vb∈W0h(Th),

where ,

 S(ub,vb) = ∑T∈ThST(ub,vb), B(ub,vb) = ∑T∈ThBT(ub,vb)

are bilinear forms in and is a linear form in .

## 3 Element Stiffness Matrices

The simplified weak Galerkin finite element method (13) is user-friendly in computer implementation. In this section, we present a formula for the computation of the element stiffness matrices and the element load vector on general polygonal elements.

{theorem}

Let be a polygonal element of sides. Denote by the vector representation of given by . Then, the element stiffness matrix and the element load vector for the SWG scheme (13) are given in a block matrix form as follows:

 (14) (κh−1AT+B+R+C)Xub≅F,

where the block components in (14) are given by:

• ,

• , with ,

• , with ,

• , with ,

• , with ,

• and ,

• and are given by

 M=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣1x1−xTy1−yT1x2−xTy2−yT⋮⋮⋮1xN−xTyN−yT⎤⎥ ⎥ ⎥ ⎥ ⎥⎦N×3,E=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣|e1||e2|⋱|eN|⎤⎥ ⎥ ⎥ ⎥ ⎥⎦N×N.

Here is any point on the plane (e.g., the center of as a specific case), is the midpoint of , is the length of edge , is the unit outward normal vector on , and is the area of the element .

From (13), the element stiffness matrix on consists of two sub-matrices corresponding to the following forms:

 ST(ub,vb) and BT(ub,vb).

The bilinear form is composed of three bilinear forms given by (11). The rest of this section is devoted to a computation of the element stiffness matrices for each of the bilinear forms involved.

### 3.1 The stiffness matrix for ST(⋅,⋅)

For the element stiffness matrix corresponding to , the key is to compute and which can be accomplished through its definition (7); readers are referred to [21] for a detailed derivation. Specifically, let be the center of T (or any point on the plane), the extension can be represented as follows:

 s(ub)=γ0+γ1(x−xT)+γ2(y−yT),

where

 (15) ⎡⎢⎣γ0γ1γ2⎤⎥⎦=(MTEM)−1MTE⎡⎢ ⎢ ⎢ ⎢ ⎢⎣ub,1ub,2⋮ub,N⎤⎥ ⎥ ⎥ ⎥ ⎥⎦.

From and (15), we have

 (16) ⎡⎢ ⎢ ⎢ ⎢ ⎢⎣s(ub)(M1)s(ub)(M2)⋮s(ub)(MN)⎤⎥ ⎥ ⎥ ⎥ ⎥⎦=M⎡⎢⎣γ0γ1γ2⎤⎥⎦=M(MTEM)−1MTE⎡⎢ ⎢ ⎢ ⎢ ⎢⎣ub,1ub,2⋮ub,N⎤⎥ ⎥ ⎥ ⎥ ⎥⎦.

Let be the basis function corresponding to the edge of :

 vb={1, on ej,0, otherwise.

Then the coefficient for is given by

 ⎡⎢⎣~γ0~γ1~γ2⎤⎥⎦=(MTEM)−1MTE⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣vb,1⋮vb,j⋮vb,N⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦=(MTEM)−1MTE⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0⋮1⋮0⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦≜⎡⎢⎣d1,jd2,jd3,j⎤⎥⎦.

It follows that

 (17)

where is the identity matrix of size .

### 3.2 The stiffness matrix for aT(⋅,⋅)

For a computation of the element stiffness matrix corresponding to the bilinear form , we have from the weak gradient formula (5) that

 (α∇wub,∇wvb)T = (α1|T|N∑j=1ub,jnj|ej|,1|T|N∑i=1vb,ini|ei|)T = N∑i,j=1(α1|T|ub,jnj|ej|,1|T|vb,ini|ei|)T = N∑i,j=1|ej||ei||T|2(αnj,ni)Tub,jvb,i, = N∑i,j=1bi,jub,jvb,i,

which leads to the block matrix in the element stiffness matrix.

### 3.3 The stiffness matrix for bT(⋅,⋅)

Recall that the bilinear form is given by

 bT(ub,vb)=(β⋅∇wub,s(vb))T.

Note that the extension has the following representation:

 s(vb)=γ0+γ1(x−xT)+γ2(y−yT),

where

 (18) ⎡⎢⎣γ0γ1γ2⎤⎥⎦=(MTEM)−1MTE⎡⎢ ⎢ ⎢ ⎢ ⎢⎣vb,1vb,2⋮vb,N⎤⎥ ⎥ ⎥ ⎥ ⎥⎦.

Thus, with , we have from the weak gradient formula (5) that

 (19) (β⋅∇wub,s(vb))T=1|T|N∑i,j=1(β⋅nj,d1,i+d2,i(x−xT)+d3,i(y−yT))T|ej|ub,jvb,i=1|T|N∑i,j=1∫Tβ⋅nj(d1,i+d2,i(x−xT)+d3,i(y−yT))dT|ej|ub,jvb,i.

For simplicity, we introduce the following functions:

 (20) ζi(x,y)=d1,i+d2,i(x−xT)+d3,i(y−yT),i=1,…,N.

Then, the equation (19) indicates that the element stiffness matrix corresponding to the bilinear form is given by

 R={rij}N×N,  rij=|ej||T|∫Tβ⋅njζidT.

### 3.4 The stiffness matrix for cT(⋅,⋅)

Recall that the bilinear form is given by

 cT(ub,vb)=(cs(ub),s(vb))T.

Thus, the element stiffness matrix corresponding to has the following formula:

 C={cij}N×N,cij=∫Tc(x,y)ζjζidT,

where is the function defined in (20).

### 3.5 The element load vector

Finally, the element load vector can be obtained from

 (f,s(vb))T = ∫Tfs(vb)dT = ∫Tf(x,y)(d1,i+d2,i(x−xT)+d3,i(y−yT))dT = ∫Tf(x,y)ζi(x,y)dT

for .

## 4 Stability and Well-Posedness

The SWG scheme (13) can be derived from the classical weak Galerkin finite element method [36, 24, 37] by eliminating the degrees of freedom associated with the interior of each element when and . But for the general case of and , the SWG finite element method (13) is different from the weak Galerkin schemes in existing literature. It is thus necessary to provide a mathematical theory for the stability and well-posedness of the numerical scheme (13).

{lemma}

Let be a shape-regular polygonal partition of the domain . There exists a constant such that

 (21) ∥∇s(vb)∥2T ≤ (22) ∥vb−s(vb)∥20,∂T ≤ Ch(∥∇wvb∥2T+h−1∥vb−Qbs(vb)∥2∂T).

Moreover, the following Poincaré-type estimate holds true:

 (23) ∥s(vb)∥2 ≤
{proof}

From the formula (6) for the weak gradient, we have for any constant vector that

 (∇wvb,ϕ)T=⟨vb,ϕ⋅n⟩∂T=⟨vb−s(vb),ϕ⋅n⟩∂T+⟨s(vb),ϕ⋅n⟩∂T=⟨vb−Qbs(vb),ϕ⋅n⟩∂T+(∇s(vb),ϕ)T,

which gives

 (∇s(vb),ϕ)T=(∇wvb,ϕ)T−⟨vb−Qbs(vb),ϕ⋅n⟩∂T.

Hence, by letting we arrive at

 ∥∇s(vb)∥2T≤C(∥∇wvb∥2T+h−1∥vb−Qbs(vb)∥2∂T),

which verifies (21).

Next, from the usual error estimate for the projection operator and the estimate (21), we have

 ∥s(vb)−Qbs(vb)∥2∂T≤Ch2∥∇s(vb)∥2∂T≤Ch∥∇s(vb)∥2T≤C(h∥∇wvb∥2T+∥vb−Qbs(vb)∥2∂T).

It follows that

 (24) ∥vb−s(vb)∥0,∂T≤∥vb−Qbs(vb)∥0,∂T+∥s(vb)−Qbs(vb)∥0,∂T≤C(h∥∇wvb∥2T+∥vb−Qbs(vb)∥2∂T)1/2,

which verifies the estimate (22).

To derive the inequality (23), we note the following discrete Poincaré inequality:

 ∥s(vb)∥2≤C∑T∈Th(∥∇s(vb)∥2T+h−1T∥s(vb)−vb∥2∂T).

Combining the above estimate with (21) and (21) gives rise to the desired inequality (23). This completes the proof of the lemma.

{lemma}

On each element , the following identity holds true:

 (25) bT(vb,vb)=12⟨vb,vbβ⋅n⟩∂T−12((∇⋅β)s(vb),s(vb))T−12⟨vb−s(vb),(vb−s(vb))β⋅n⟩∂T+⟨vb−s(vb),¯¯¯¯¯¯¯¯¯¯¯¯¯¯s(vb)β⋅n−s(vb)β⋅n⟩∂T,

where is the average of on the element .

{proof}

From the formula (5), we have

 (26) bT(vb,vb)=(β⋅∇wvb,s(vb))T=(∇wvb,s(vb)β)T=(∇wvb,¯¯¯¯¯¯¯¯¯¯¯¯¯¯s(vb)β)T=⟨vb,¯¯¯¯¯¯¯¯¯¯¯¯¯¯s(vb)β⋅n⟩∂T=⟨vb−s(vb),¯¯¯¯¯¯¯¯¯¯¯¯¯¯s(vb)β⋅n⟩∂T+⟨s(vb),¯¯¯¯¯¯¯¯¯¯¯¯¯¯s(vb)β⋅n⟩∂T.

Note that

 ⟨s(vb),¯¯¯¯¯¯¯¯¯¯¯¯¯¯s(vb)β⋅n⟩∂T=(∇s(vb),¯¯¯¯¯¯¯¯¯¯¯¯¯¯s(vb)β)T=(∇s(vb),s(vb)β)T=12⟨s(vb),s(vb)β⋅n⟩∂T−12((∇⋅β)s(vb),s(vb))T.

Substituting the above identity into (26) yields

 (27) bT(vb,vb)=⟨vb−s(vb),¯¯¯¯¯¯¯¯¯¯¯¯¯¯s(vb)β⋅n⟩∂T+12⟨s(vb),s(vb)β⋅n⟩∂T  −12((∇⋅β)s(vb),s(vb))T=⟨vb−s(vb),¯¯¯¯¯¯¯¯¯¯¯¯¯¯s(vb)β⋅n−s(vb)β⋅n⟩∂T+⟨vb,s(vb)β⋅n⟩∂T  −12⟨s(vb),s(vb)β⋅n⟩∂T−12((∇⋅β)s(vb),s(vb))T=⟨vb−s(vb),¯¯¯¯¯¯¯¯¯¯¯¯¯¯s(vb)β⋅n−s(vb)β⋅n⟩∂T −12⟨vb−s(vb),(vb−s(vb))β⋅n⟩∂T +12⟨vb,vbβ⋅n⟩∂T−12((∇⋅β)s(vb),s(vb))T,

which leads to the identify (25).

In the finite element space , we introduce the following semi-norm:

 (28) |||vb|||2:=∑T∈Th(κST(vb,vb)+aT(vb,vb))

We claim that defines a norm in the closed subspace . It suffices to show that for any satisfying . In fact, if , then from (28) we have

 κ∑TST(vb,vb)+∑T(α∇wvb,∇wvb)T=0.

It follows that on each element

 (29) ∇wvb=0,(vb−s(vb))(Mi)=0

for . Thus,

so that has constant value on each element . By using (29) we see that on each edge, which, together with the fact that on , leads to in .

{lemma}

For the model problem (1), assume that and the condition (3) is satisfied. Then, the bilinear form is bounded and coercive in the finite element space ; i.e., there exist constants and such that

 (30) |κS(vb,wb)+B(vb,wb)| ≤ M|||vb||||||wb|||∀vb,wb∈W0h(Th), (31) κS(vb,vb)+B(vb,vb) ≥ Λ|||vb|||2∀vb∈W0h(Th),

provided that the meshsize of is sufficiently small.

{proof}

Recall that for any we have

 (32) B(vb,wb) =∑T∈Th(aT(vb,wb)+bT(vb,wb)+cT(vb,wb)),S(vb,wb) =∑T∈ThST(vb,wb).

The boundedness estimate (30) is then straightforward from the usual Cauchy-Schwarz and the inequality (23). We shall focus on the derivation of the coercivity inequality (31) in the rest of the proof.

In comparison with (28), the key to the coercivity inequality (31) is to derive an estimate of the following type:

 (33) ∑T∈Th(bT(vb,vb)+cT(vb,vb))≥η−ε(h)|||vb|||2,

where and is a parameter satisfying