# An L^p- Primal-Dual Weak Galerkin Method for Convection-Diffusion Equations

In this article, the authors present a new L^p- primal-dual weak Galerkin method (L^p-PDWG) for convection-diffusion equations with p>1. The existence and uniqueness of the numerical solution is discussed, and an optimal-order error estimate is derived in the L^q-norm for the primal variable, where 1/p+1/q=1. Furthermore, error estimates are established for the numerical approximation of the dual variable in the standard W^m,p norm, 0≤ m≤ 2. Numerical results are presented to demonstrate the efficiency and accuracy of the proposed L^p-PDWG method.

## Authors

• 8 publications
• 18 publications
• 14 publications
01/19/2020

10/30/2019

### Primal-dual weak Galerkin finite element methods for linear convection equations in non-divergence form

A new primal-dual weak Galerkin (PD-WG) finite element method was develo...
03/22/2021

### Analysis of backward Euler primal DPG methods

We analyse backward Euler time stepping schemes for the primal DPG formu...
06/06/2021

### An L^p-weak Galerkin method for second order elliptic equations in non-divergence form

This article presents a new primal-dual weak Galerkin method for second ...
03/26/2020

08/02/2019

02/09/2019

### An Optimal-Storage Approach to Semidefinite Programming using Approximate Complementarity

This paper develops a new storage-optimal algorithm that provably solves...
##### 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, the authors are concerned with the development of an - primal-dual weak Galerkin (-PDWG) finite element method for second order elliptic boundary value problems that seek such that

 (1.1) −Δu+∇⋅(βu)=f,in Ω,u=g,on ∂Ω,

where () is an open bounded and connected domain with piecewise smooth Lipschitz boundary . By introducing the space , we obtain the following weak form for the model problem (1.1) following an use of the usual integration by parts: Find such that

 (1.2) (u,Δσ+β⋅∇σ)=⟨g,∇σ⋅n⟩∂Ω−(f,σ),∀σ∈V.

Assume the model problem (1.1) has one and only one solution for any given and in appropriate spaces.

The weak Galerkin (WG) method was first introduced by Wang and Ye in [18] for second order elliptic equations, where weak gradient and its discrete weak gradient were constructed to replace the standard gradient and its discrete gradient. Later, the authors in [8] developed a primal-dual weak Galerkin (PDWG) finite element method for the second order elliptic problem in non-divergence form, where a discrete weak Hessian operator in the weak formulation of the model PDEs was designed. This PDWG algorithm can be characterized as a constrained optimization problem with constraints given by the weak formulation of the model PDEs weakly defined on each element. In the past several years, many theoretical a priori error estimates for weak Galerkin methods have been established in and discrete norms. Readers are referred to [12, 13, 14, 18, 19, 6, 16, 7, 15, 9, 10, 11] for an incomplete list of references.

The purpose of this paper is to present an - primal-dual weak Galerkin method for the problem (1.1), and establish a general and theory for the numerical method. To our best knowledge, there is one existing result in the error estimate for the mixed finite element method developed by Duran [3] for second order elliptic problems in , but no results in are known for the weak Galerkin finite element methods in the literature. Different from the method in [8], our numerical scheme is based on the weak formulation (1.2) together with a weak version of the dual operator applied to the test functions. The new PDWG method can be characterized as a constrained optimization problem with constraints that satisfy the PDE weakly on each element, which extends the idea of minimization problem in [8] to a more general setting.

To conclude this section, we point out that our theory for the primal-dual weak Gelerkin finite element method is based on the assumption that the solution to the following adjoint problem

 (1.3) −Δφ−β⋅∇φ= χ, in Ω,φ= 0, on ∂Ω,

is -regular in the sense that it has a unique solution in and the solution satisfies

 ∥φ∥2,p≤C∥χ∥0,p.

Under this assumption, we shall derive an optimal order error estimate in the standard norm for the primal variable and the standard norms for the dual variable. Numerical experiments demonstrate that our error estimate for the primal variable is optimal; i.e., the error bound is sharp.

The rest of this paper is organized as follows. In Section 2, we briefly review the weak differential operators and their discrete versions. In Section 3, the primal-dual weak Galerkin scheme is introduced for the model problem (1.1) based on theory. Section 4 is devoted to the establishment of the solution existence, uniqueness and stability. In Section 5, we derive an error equations for our numerical methods, which is of essential importance in our later error estimates. Section 6 and Section 7 establish error estimates for the primal variable in norms and for the dual variable in norms, respectively. Finally, a series of numerical examples are presented in Section 8 to verify the mathematical convergence theory.

## 2 Weak Differential Operators

The Laplacian and the gradient are the principle differential operators used in the weak formulation (1.2) for the second order elliptic model problem (1.1). This section gives a brief discussion of the weak Laplacian and gradient operators as well as their discrete analogies [5].

Let be a polygonal or polyhedral domain with boundary . A weak function on refers to a triplet such that , , and . The first and second components and can be identified as the value of in the interior and on the boundary of . The third component is meant to represent the value of on the boundary of the element . Note that and might be totally independent of the trace of and on , respectively. Denote by the space of all scalar-valued weak functions on ; i.e.,

 (2.1) W(T)={σ={σ0,σb,σn}:σ0∈Lp(T),σb∈Lp(∂T),σn∈Lp(∂T)}.

The weak Laplacian operator, denoted by , is defined as a linear functional in such that

 (Δwσ,w)T=(σ0,Δw)T−⟨σb,∇w⋅n⟩∂T+⟨σn,w⟩∂T,

for all .

Denote by the space of all polynomials on with degree no more than . A discrete analogy of for is defined as the unique polynomial satisfying

 (2.2) (Δw,r,Tσ,w)T=(σ0,Δw)T−⟨σb,∇w⋅n⟩∂T+⟨σn,w⟩∂T,∀w∈Pr(T).

For smooth such that , we have from the integration by parts

 (2.3) (Δw,r,Tσ,w)T=(Δσ0,w)T+⟨σ0−σb,∇w⋅n⟩∂T−⟨∇σ0⋅n−σn,w⟩∂T,

for all . Similarly, the discrete weak gradient operator is defined as the unique polynomial satisfying

 (2.4) (∇w,r,Tσ,φ)T=−(σ0,∇⋅φ)T+⟨σb,φ⋅n⟩∂T,∀φ∈[Pr(T)]d.

When , the following identify holds true:

 (2.5) (∇w,r,Tσ,φ)T=(∇σ0,φ)T+⟨σb−σ0,φ⋅n⟩∂T,

for all .

## 3 Numerical Algorithm

Denote by a partition of the domain into polygons in 2D or polyhedra in 3D which is shape regular in the sense described in [18]. Denote by the set of all edges or flat faces in and the set of all interior edges or flat faces. Denote by the meshsize of and the meshsize for the partition .

For any given integer , denote by the local discrete space of the weak functions defined by

 Wk(T)={{σ0,σb,σn}:σ0∈Pk(T),σb∈Pk(e),σn∈Pk−1(e),e⊂∂T}.

Patching over all the elements through a common value of and on the interior interface yields a weak finite element space :

 (3.1) Wh={{σ0,σb,σn}:{σ0,σb,σn}|T∈Wk(T),∀T∈Th}.

Note that has two values and on each interior interface as seen from the two elements and , and they must satisfy . Denote by the subspace of with homogeneous boundary condition; i.e.,

 W0h={{σ0,σb,σn}:{σ0,σb,σn}|T∈Wh,σb|∂Ω=0, ∀e∈∂T,T∈Th}.

Denote by the finite element space consisting of piecewise polynomials of degree where ; i.e.,

 (3.2) Mh={w:w|T∈Ps(T),∀T∈Th}.

We emphasize that both the weak gradient and the weak Laplacian operators are defined by using piecewise polynomials of degree . For purely diffusive equations, one may assume the value of .

For simplicity of notation and without confusion, for any , denote by and the discrete weak Laplacian and discrete weak gradient computed by (2.2) and (2.4) on each element , respectively; i.e.,

 (Δwσ)|T=Δw,s,T(σ|T),(∇wσ)|T=∇w,s,T(σ|T), s=k−1.

For any and , we introduce the following forms

 (3.3) s(λ,σ)= ∑T∈ThsT(λ,σ), (3.4) b(u,λ)= ∑T∈ThbT(u,λ),

where

 sT(λ,σ)=h1−2pT∫∂T|λ0−λb|p−1sgn(λ0−λb)(σ0−σb)ds+h1−pT∫∂T|∇λ0⋅n−λn|p−1sgn(∇λ0⋅n−λn)(∇σ0⋅n−σn)dsbT(u,λ)=(u,−β⋅∇wλ−Δwλ)T.

The numerical scheme for the second order elliptic model problem (1.1) based on the variational formulation (1.2) can be stated as follows:

###### Primal-Dual Weak Galerkin Algorithm 3.1

Find , such that

 (3.5) s(λh,σ)+b(uh,σ) = (f,σ0)−⟨g,σn⟩∂Ω,∀σ∈W0h, (3.6) b(v,λh) = 0,∀v∈Mh.

In next section, we shall study the solution existence and uniqueness for the primal-dual weak Galerkin finite element algorithm (3.5)-(3.6). For simplicity of analysis, we assume constant value for the convection term on each element in the rest of the paper.

## 4 Solution Existence and Uniqueness

Denote by the projection operator onto for each element . For each edge or face , denote by and the projection operators onto and , respectively. For any , define by the projection onto the weak finite element space such that on each element ,

 Qhw={Q0w,Qbw,Qn(∇w⋅n)}.

Denote by the projection onto the finite element space .

[5] The projection operators and satisfy the following commutative properties:

 (4.1) Δw(Qhw) = Qh(Δw),w∈W2,p(T), (4.2) ∇w(Qhv) = Qh(∇v),v∈W1,p(T).

To show the existence of solutions, we consider the functional

 J(σ,v):=1ps(σ,σ)+b(v,σ)−(F,σ),

where and . If is the solution of (3.5)-(3.6), then we have from (3.6)

 J(λh,v)=1ps(λh,λh)−(F,λh)=J(λh,uh),  ∀v∈Mh.

On the other hand, the equation (3.5) indicates that for all , where is the Gateaux partial derivative at in the direction of . It follows that is a global minimizer of the functional ; i.e.,

 J(λh,uh)≤J(σ,uh),∀σ∈W0h.

Consequently,

 (4.3) J(λh,v)≤J(λh,uh)≤J(σ,uh),∀v∈Mh,σ∈W0h.

The above inequality implies that the solution is a saddle point of the functional . Thus, (3.5)-(3.6) can be formulated as the following min-max problem: Find and such that

 (λh,uh)=argminσ∈W0hmaxv∈MhJ(σ,v).

As a convex minimization problem, the above problem has a solution so that there must be a solution satisfying (3.5)-(3.6).

The rest of this section is devoted to a discussion of the uniqueness of the numerical solution .

The numerical scheme (3.5)-(3.6) has one and only one solution in the finite element space . Let and be two solutions of (3.5)-(3.6). Denote

 ϵh=λ(1)h−λ(2)h={ϵ0,ϵb,ϵn},  eh=u(1)h−u(2)h.

For any constants , we choose in (3.5) and use (3.6) to obtain

 s(λ(1)h,θ1λ(1)h+θ2λ(2)h)−s(λ(2)h,θ1λ(1)h+θ2λ(2)h)=0.

In particular, by taking , we have

 (4.4) s(λ(1)h,λ(1)h)=s(λ(2)h,λ(1)h),  s(λ(2)h,λ(2)h)=s(λ(1)h,λ(2)h),

which yields, together with Young’s inequality , that

 s(λ(1)h,λ(1)h)≤s(λ(2)h,λ(2)h)q+s(λ(1)h,λ(1)h)p,  s(λ(2)h,λ(2)h)≤s(λ(1)h,λ(1)h)q+s(λ(2)h,λ(2)h)p,

which yields

 (4.5) s(λ(1)h,λ(1)h)=s(λ(2)h,λ(2)h).

On the other hand, for any two real numbers , there holds

 ∣∣∣A+B2∣∣∣p≤(|A|p+|B|p)/2,

and the equality holds true if and only if . It follows that

 (4.6) s(λ(1)h+λ(2)h2,λ(1)h+λ(2)h2)≤12(s(λ(1)h,λ(1)h)+s(λ(2)h,λ(2)h))=s(λ(1)h,λ(1)h).

By (4.4)-(4.5) and Young’s inequality,

 s(λ(1)h,λ(1)h) = 12(s(λ(1)h,λ(1)h)+s(λ(1)h,λ(2)h))=s(λ(1)h,λ(1)h+λ(2)h2) ≤ 1qs(λ(1)h,λ(1)h)+1ps(λ(1)h+λ(2)h2,λ(1)h+λ(2)h2),

which indicates that

 s(λ(1)h,λ(1)h)≤s(λ(1)h+λ(2)h2,λ(1)h+λ(2)h2).

In light of (4.6), we easily obtain that

 s(λ(1)h+λ(2)h2,λ(1)h+λ(2)h2)=s(λ(1)h,λ(1)h)=s(λ(2)h,λ(2)h).

The above equality holds true if and only if

 λ(1)0−λ(1)b = λ(2)0−λ(2)b, on ∂T, ∇λ(1)0⋅n−λ(1)n = ∇λ(2)0⋅n−λ(2)n, on ∂T,

or equivalently,

 (4.7) ϵ0 = ϵb, on ∂T, (4.8) ∇ϵ0⋅n = ϵn, on ∂T.

Let and be two solutions of (3.5)-(3.6). We have from (3.6) that . Using (2.3) and (2.5), we have

 0=b(v,ϵh)=∑T∈Th(v,−β⋅∇wϵh−Δwϵh)T=∑T∈Th−(∇ϵ0,βv)T−⟨ϵb−ϵ0,βv⋅n⟩∂T−(Δϵ0,v)−⟨ϵ0−ϵb,∇v⋅n⟩∂T+⟨∇ϵ0⋅n−ϵn,v⟩∂T=∑T∈Th(−β⋅∇ϵ0−Δϵ0,v)T,

where we used (4.7) -(4.8), which gives on each . Together with (4.7) -(4.8), we arrive at in , with the boundary condition on due to the fact that (4.7) and . Therefore, we obtain in . Furthermore, we have and , which leads to .

We next show . To this end, using and the equation (3.5) we obtain

 b(eh,σ)=s(λ(1)h,σ)−s(λ(2)h,σ)+b(eh,σ)=0,∀σ∈W0h,

which, together with the definition of the weak Laplacian and the weak gradient , yields

 0=b(eh,σ) = ∑T∈Th(eh,−β⋅∇wσ−Δwσ)T = ∑T∈Th(∇⋅(βeh)−Δeh,σ0)T +∑T∈Th⟨σb,∇eh⋅n⟩∂T−⟨σn,eh⟩∂T−⟨eh,β⋅nσb⟩∂T = ∑T∈Th(∇⋅(βeh)−Δeh,σ0)T +∑e∈Eh∫e([[∇eh−βeh]]⋅neσb−[[eh]]σn)ds

for all , where is the assigned outward normal direction to and is the jump on the edge . In particular, by taking , , and we obtain on each

 −△eh+∇⋅(βeh)=0, % in T, [[eh]]=0,  [[∇eh−βeh]]⋅ne=0, on ∂T.

Consequently, from the solution uniqueness for (1.1) we have

 eh≡0,  or equivalently  u(1)h=u(2)h.

This completes the proof.

## 5 Error Equation

Let and be the exact solution of (1.1) and its numerical solution arising from the PDWG scheme (3.5)-(3.6), respectively. Note that the exact solution of the Lagrangian multiplier is . Denote two error functions by

 (5.1) eh =Qhu−uh, (5.2) ϵh =Qhλ−λh=−λh.

Let and be the exact solution of (1.1) and its numerical solution arising from PDWG scheme (3.5)-(3.6). The error functions and satisfy the following equations:

 (5.3) s(ϵh,σ)+b(eh,σ) = lu(σ),∀σ∈W0h, (5.4) b(v,ϵh) = 0,∀v∈Mh.

Here

 (5.5) lu(σ)=∑T∈Th⟨u−Qhu,σn−∇σ0⋅n⟩∂T+⟨∇u⋅n−∇Qhu⋅n,σ0−σb⟩∂T+∑T∈Th⟨(u−Qhu)β⋅n,σb−σ0⟩∂T.

First, from (5.2) and (3.6) we may readily derive (5.4). Next, by using (3.4) for and choosing in (2.3) and (2.5), we have

 (5.6) b(Qhu,σ)=∑T∈Th(Qhu,−β⋅∇wσ−Δwσ)T=∑T∈Th(Qhu,−β⋅∇σ0−Δσ0)T−⟨Qhuβ⋅n,σb−σ0⟩∂T+⟨∇Qhu⋅n,σb−σ0⟩∂T+⟨Qhu,∇σ0⋅n−σn⟩∂T=∑T∈Th(u,−β⋅∇σ0−Δσ0)T−⟨Qhuβ⋅n,σb−σ0⟩∂T+⟨∇Qhu⋅n,σb−σ0⟩∂T+⟨Qhu,∇σ0⋅n−σn⟩∂T.

Now applying the usual integration by parts to the integrals on yields

 (5.7) b(Qhu,σ)=∑T∈Th(−Δu+∇⋅(βu),σ0)T+⟨(u−Qhu)β⋅n,σb−σ0⟩∂T+⟨∇(Qhu−u)⋅n,σb−σ0⟩∂T+⟨(Qhu−u),(∇σ0⋅n−σn⟩∂T−⟨g,σn⟩∂Ω= (f,σ0)T−⟨g,σn⟩∂Ω+∑T∈Th⟨∇(Qhu−u)⋅n,σb−σ0⟩∂T+⟨(Qhu−u),(∇σ0⋅n−σn⟩∂T+∑T∈Th⟨(u−Qhu)β⋅n,σb−σ0⟩∂T,

where we have used and on . From , we have . Subtracting (3.5) from (5.7) yields the error equation (5.3). This completes the proof of the lemma.

The equations (5.3)-(5.4) are called error equations for the primal-dual WG finite element scheme (3.5)-(3.6).

## 6 Lq-Error Estimate for the Primal Variable uh

Recall that is a shape-regular finite element partition of the domain . For any and with , the following trace inequality holds true:

 (6.1) ∥w∥qLq(∂T)≤Ch−1T(∥w∥qLq(T)+hqT∥∇w∥qLq(T)).

Let and . Let and be the exact solution of the second order elliptic model problem (1.1) and the numerical solution arising from PDWG scheme (3.5)-(3.6). The following error estimate holds true:

 (6.2) s(λh,λh)≤Chqk∥∇ku∥qLq(Ω).

Recall that . By letting in (5.3), we have from (5.4) and (5.5) that

 (6.3) s(λh,λh)=∑T∈Th⟨u−Qhu,∇λ0⋅n−λn⟩∂T+⟨∇u⋅n−∇Qhu⋅n,λb−λ0⟩∂T+∑T∈Th⟨(u−Qhu)β⋅n,λ0−λb⟩∂T.

For the first term on the right-hand side of (6.3), we use the Cauchy-Schwarz inequality to obtain

 (6.4) ∣∣∑T∈Th⟨Qhu−u,λn−∇λ0⋅n⟩∂T∣∣≤(∑T∈Th∥u−Qhu∥qLq(∂T))1q(∑T∈Th∥λn−∇λ0⋅n∥pLp(∂T))1p.

For the term