A Weak Galerkin Finite Element Scheme for the Biharmonic Equations by Using Polynomials of Reduced Order

# A Weak Galerkin Finite Element Scheme for the Biharmonic Equations by Using Polynomials of Reduced Order

## Abstract

A new weak Galerkin (WG) finite element method for solving the biharmonic equation in two or three dimensional spaces by using polynomials of reduced order is introduced and analyzed. The WG method is on the use of weak functions and their weak derivatives defined as distributions. Weak functions and weak derivatives can be approximated by polynomials with various degrees. Different combination of polynomial spaces leads to different WG finite element methods, which makes WG methods highly flexible and efficient in practical computation. This paper explores the possibility of optimal combination of polynomial spaces that minimize the number of unknowns in the numerical scheme, yet without compromising the accuracy of the numerical approximation. Error estimates of optimal order are established for the corresponding WG approximations in both a discrete norm and the standard norm. In addition, the paper also presents some numerical experiments to demonstrate the power of the WG method. The numerical results show a great promise of the robustness, reliability, flexibility and accuracy of the WG method.

w

eak Galerkin finite element methods, weak Laplacian, biharmonic equation, polyhedral meshes.

{AMS}

Primary, 65N30, 65N15, 65N12, 74N20; Secondary, 35B45, 35J50, 35J35

## 1 Introduction

This paper will concern with approximating the solution of the biharmonic equation

 (1.1) Δ2u = f,in Ω,

with clamped boundary conditions

 (1.2) u = g,on ∂Ω, (1.3) ∂u∂n = ϕ,on ∂Ω,

where is the Laplacian operator, is a bounded polygonal or polyhedral domain in for and n denotes the outward unit normal vector along . We assume that are given, sufficiently smooth functions.

This problem mainly arises in fluid dynamics where the stream functions of incompressible flows are sought and elasticity theory, in which the deflection of a thin plate of the clamped plate bending problem is sought [27, 36, 39].

Due to the significance of the biharmonic problem, a large number of methods for discretizing (1.1) - (1.3) have been proposed. These methods include dealing with the biharmonic operator directly, such as discretizing (1.1)-(1.3) on a uniform grid using a 13-point or 25-point direct approximation of the fourth order differential operator [9, 25]; mixed methods, that is, splitting the biharmonic equation into two coupled Poisson equations [1, 5, 18, 4, 13, 16, 19, 20, 21, 28, 26, 6, 7]. Also there are some other approaches to the biharmonic problems, like the conformal mapping methods [12, 37], integral equations , orthogonal spline collocation method  and the fast multipole methods , etc.

Among these methods, finite element methods are one of the most widely used technique, which is based on variational formulations of the equations considered. In fact, the biharmonic equation is also one of the most important applicable problems of the finite element methods, cf. [17, 23, 2, 44, 14, 15]. The Galerkin methods, discretizing the corresponding variational form of (1.1) is given by seeking satisfying

 u|∂Ω=g,∂u∂n|∂Ω=ϕ

such that

 (1.4) (Δu,Δv)=(f,v),∀v∈H20(Ω),

where is the subspace of consisting of functions with vanishing value and normal derivative on .

Standard finite element methods for solving (1.1) - (1.3) based on the variational form (1.4) with conforming finite element require rather sophisticated finite elements such as the 21-degrees-of-freedom of Argyris (see ) or nonconforming elements of Hermite type. Since the complexity in the construction for the finite element with high continuous elements, conforming element are seldom used in practice for the biharmonic problem. To avoid using of -elements, besides the mixed methods, an alternative approach, nonconforming and discontinuous Galerkin finite element methods have been developed for solving the biharmonic equation over the last several decades. Morley element  is a well known nonconforming element for the biharmonic equation for its simplicity. A interior penalty method was developed in [10, 22]. In , a hp-version interior penalty discontinuous Galerkin method was presented for the biharmonic equation.

Recently a new class of finite element methods, called weak Galerkin(WG) finite element methods were developed for the biharmonic equation for its highly flexible and robust properties. The WG method refers to a numerical scheme for partial differential equations in which differential operators are approximated by weak forms as distributions over a set of generalized functions. This thought was first proposed in  for a model second order elliptic problem, and this method was further developed in [42, 32, 43]. In , a weak Galerkin method for the biharmonic equation was derived by using discontinuous functions of piecewise polynomials on general partitions of polygons or polyhedra of arbitrary shape. After that, in order to reduce the number of unknowns, a WG method  was proposed and analyzed. However, due to the continuity limitation, the WG scheme only works for the traditional finite partitions, while not arbitrary polygonal or polyhedral girds as allowed in .

In order to realize the aim that reducing the unknown numbers and suit for general partitions of polygons or polyhedra of arbitrary shape at the same time, in this paper we construct a reduction WG scheme based on the use of a discrete weak Laplacian plus a new stabilization that is also parameter free. The goal of this paper is to specify all the details for the reduction WG method for the biharmonic equations and present the numerical analysis by presenting a mathematical convergence theory.

An outline of the paper is as follows. In the remainder of the introduction we shall introduce some preliminaries and notations for Sobolev spaces. In Section 2 is devoted to the definitions of weak functions and weak derivatives. The WG finite element schemes for the biharmonic equation (1.1)-(1.3) are presented in Section 3. In Section 4, we establish an optimal order error estimates for the WG finite element approximation in an equivalent discrete norm. In Section 5, we shall drive an error estimate for the WG finite element method in the standard norm. Section 6 contains the numerical results of the WG method. The theoretical results are illustrated by these numerical examples. Finally, we present some technical estimates for quantities related to the local projections into various finite element spaces and some approximation properties which are useful in the convergence analysis in Appendix A.

Now let us define some notations. Let be any open bounded domain with Lipschitz continuous boundary in . We use the standard definition for the Sobloev space and their associated inner products , norms , and seminorms for any .

The space coincides with , for which the norm and the inner product are denoted by and , respectively. When , we shall drop the subscript in the norm and in the inner product notation.

The space is defined as the set of vector-valued functions on which, together with their divergence, are square integrable; i.e.,

 H(div;D)={v:v∈[L2(D)]d,∇⋅v∈L2(D)}.

The norm in is defined by

 ∥v∥H(div;D)=(∥v∥2D+∥∇⋅v% ∥2D)12.

## 2 Weak Laplacain and Discrete Weak Laplacian

For the biharmonic equation (1.1), the underlying differential operator is the Laplacian . Thus, we shall first introduce a weak version for the Laplacian operator defined on a class of discontinuous functions as distributions .

Let be any polygonal or polyhedral domain with boundary . A weak function on the region refers to a function such that , , and , where n is the outward unit normal vector along . Denote by the space of all weak functions on , that is,

 (2.1) W(K)={v={v0,vb,vg}:v0∈L2(K),vb,vg⋅n∈L2(∂K)}.

Recall that, for any , the weak Laplacian of is defined as a linear functional in the dual space of whose action on each is given by

 (2.2) (Δwv,φ)K=(v0,Δφ)K−⟨vb,∇φ⋅n⟩∂K+⟨vg⋅n,φ⟩∂K,

where stands for the -inner product in and is the inner product in .

The Sobolev space can be embedded into the space by an inclusion map defined as follows

 iW(ϕ)={ϕ|K,ϕ|∂K,(∇ϕ⋅n% )n|∂K},ϕ∈H2(K).

With the help of the inclusion map , the Sobolev space can be viewed as a subspace of by identifying each with .

Analogously, a weak function is said to be in if it can be identified with a function through the above inclusion map. Here the first components can be seen as the value of in the interior and the second component represents the value of on . Denote by , then the third component represents . Obviously, . Note that if , then and may not necessarily be related to the trace of and on , respectively.

For , from integration by parts we have

 (Δwv,φ)K = (v,Δφ)K−⟨v,∇φ⋅n⟩∂K+⟨∇v⋅n,φ⟩∂K = (v0,Δφ)K−⟨vb,∇φ⋅%n⟩∂K+⟨vg⋅n,φ⟩∂K.

Thus the weak Laplacian is identical with the strong Laplacian, i.e.,

 ΔwiW(v)=Δv

for smooth functions in .

For numerical implementation purpose, we define a discrete version of the weak Laplacain operator by approximating in polynomial subspaces of the dual of . To this end, for any non-negative integer , let be the set of polynomials on with degree no more than .

{definition}

() A discrete weak Laplacian operator, denoted by , is defined as the unique polynomial satisfying

 (2.3) (Δw,r,Kv,φ)K=(v0,Δφ)K−⟨vb,∇φ⋅n⟩∂K+⟨% \rm\bf vn⋅n,φ⟩∂K,∀φ∈Pr(K).

From the integration by parts, we have

 (v0,Δφ)K=(Δv0,φ)K+⟨v0,∇φ⋅n⟩∂K−⟨∇v0⋅n,φ⟩∂K.

Substituting the above identity into (2.3) yields

 (2.4) (Δw,r,Kv,φ)K−(Δv0,φ)K=⟨v0−vb,∇φ⋅n⟩∂K−⟨(∇v0−vg)⋅n,φ⟩∂K,

for all .

## 3 Weak Galerkin Finite Element Scheme

Let be a partition of the domain into polygons in 2D or polyhedra in 3D. Assume that is shape regular in the sense as defined in . Denote by the set of all edges or flat faces in , and let be the set of all interior edges or flat faces.

Since represents , then is naturally dependent on n. To ensure a single valued function on , we introduce a set of normal directions on as follows

 (3.1) Nh={ne:ne is unit and normal to e, e∈Eh}.

For any given integer , denote by the discrete weak function space given by

 (3.2) Wk(T)={{v0,vb,vnne}:v0∈Pk(T),vb,vn∈Pk−1(e),e⊂∂T}.

By patching over all the elements through a common value on the interface , we arrive at a weak finite element space given by

 Vh={{v0,vb,vnne}:{v0,vb,vnne}∣∣T∈Wk(T),∀T∈Th}.

Denote by the subspace of constituting discrete weak functions with vanishing traces; i.e.,

 V0h={{v0,vb,vnne}:{v0,vb,vnne}∈Vh,vb|e=0,vn|e=0,e∈∂T∩∂Ω}.

Denote by the trace of on from the component . It is obvious that consists of piecewise polynomials of degree . Similarly, denote by the trace of from the component of as piecewise polynomials of degree . Denote by the discrete weak Laplacian operator on the finite element space computed by using (2.3) on each element for , that is,

 (3.3) (Δw,k−2v)|T=Δw,k−2,T(v|T)∀v∈Vh.

For simplicity, we shall drop the subscript in the notation for the discrete weak Laplacian operator. We also introduce the following notation

 (Δwv,Δww)h=∑T∈Th(Δwv,Δww)T.

For each element , denote by the projection onto , . For each edge/face , denote by the projection onto . Now for any , we shall combine these two projections together to define a projection into the finite element space such that on the element

 Qhu={Q0u,Qbu,(Qb(∇u⋅ne))ne}.
{theorem}

Let be the local projection onto . Then the following commutative diagram holds true on each element :

 (3.4) ΔwQhu=QhΔu,∀u∈H2(T).
{proof}

For any , from the definition of the discrete weak Laplacian and the projection

 (ΔwQhu,ϕ)T = (Q0u,Δϕ)T−⟨Qbu,∇ϕ⋅n⟩∂T+⟨Qb(∇u⋅ne)ne⋅n,ϕ⟩∂T = (u,Δϕ)T−⟨u,∇ϕ⋅n⟩∂T+⟨∇u⋅n,ϕ⟩∂T = (Δu,ϕ)T=(QhΔu,ϕ),

which implies (3.4).

The commutative property (3.4) indicates that the discrete weak Laplacian of the projection of is a good approximation of the Laplacian of in the classical sense. This is a good property of the discrete weak Laplacian in application to algorithm and analysis.

For any and in , we introduce a bilinear form as follows

 s(uh,v) = ∑T∈Thh−1T⟨∇u0⋅ne−un,∇v0⋅ne−vn⟩∂T +∑T∈Thh−3T⟨Qbu0−ub,Qbv0−vb⟩∂T.
###### Weak Galerkin Algorithm 1

Find satisfying and on and the following equation:

 (3.5) (Δwuh,Δwv)h+s(uh,v)=(f,v0), ∀v={v0,vb,vnne}∈V0h.
{lemma}

For any , let be given by

 (3.6) |||v|||2=(Δwv,Δwv)h+s(v,v).

Then, defines a norm in the linear space .

{proof}

For simplicity, we shall only prove the positivity property for . Assume that for some . It follows that on each element T, and on each edge . We claim that holds true locally on each element T. To this end, for any we use and the identify (2.4) to obtain

 0 = (Δwv,φ)T = (Δv0,φ)T+⟨Qbv0−vb,∇φ⋅n⟩∂T+⟨vnne⋅n−∇v0⋅n,φ⟩∂T = (Δv0,φ)T,

where we have used the fact that and

 vnne⋅n−∇v0⋅n=±(vn−∇v0⋅ne)=0

in the last equality. The identity (3) implies that holds true locally on each element .

Next, we claim that also holds true locally on each element . For this purpose, for any , we utilize the Gauss formula to obtain

 (3.8) (∇v0,∇ϕ)T = −(Δv0,ϕ)T+⟨∇v0⋅n,ϕ⟩∂T=⟨∇v0⋅n,ϕ⟩∂T.

By letting on each element and summing over all we obtain

 (3.9) ∑T∈Th(∇v0,∇v0)T=∑T∈Th⟨∇v0⋅n,v0⟩∂T.

For two elements , which share as a common edge, denote the values of in the interior of , respectively. It follows from on edge and the fact that

 ⟨∇v10⋅nT1,v10⟩e+⟨∇v20⋅nT2,v20⟩e = ±⟨vn,v10−v20⟩e=±⟨vn,Qbv10−Qbv20⟩e=0,

where denote the outward unit normal vectors on according to elements , respectively. This, together with on the boundary edge implies

 ∑T∈Th⟨∇v0⋅n,v0⟩∂T=0.

It follows from equation (3.9) that on each element . Thus, locally on each element and is then continuous across each interior edge as

 v0|e=Qbv0=vb.

The boundary condition of then implies that on , which completes the proof.

{lemma}

The weak Galerkin finite element scheme (3.5) has a unique solution.

{proof}

Assume and are two solutions of the WG finite element scheme (3.5). It is obvious that the difference is a finite element function in satisfying

 (3.10) (Δwρh,Δwv)h+s(ρh,v)=0, ∀v∈V0h.

By letting in above equation (3.10) we obtain the following indentity

 (Δwρh,Δwρh)h+s(ρh,ρh)=0.

It follows from Lemma 3 that , which shows that . This completes the proof.

## 4 An Error Estimate

The goal of this section is to establish an error estimate for the WG-FEM solution arising from (3.5).

First of all, let us derive an error equation for the WG finite element solution obtained from (3.5). This error equation is critical in convergence analysis.

{lemma}

Let and be the solution of (1.1)-(1.3) and (3.5), respectively. Denote by

 eh=Qhu−uh

the error function between the projection of and its weak Galerkin finite element solution. Then the error function satisfies the following equation

 (4.1) (Δωeh,Δωv)h+s(eh,v)=ℓu(v)

for all . Here

 ℓu(v) = ∑T∈Th⟨Δu−QhΔu,∇v0⋅n−vnne⋅n⟩∂T −∑T∈Th⟨∇(Δu−QhΔu)⋅n,v0−vb⟩∂T+s(Qhu,v).
{proof}

Using (2.4) with we obtain

 (ΔωQhu,Δωv)T = (Δv0,QhΔu)T+⟨v0−vb,∇(QhΔu)⋅n⟩∂T−⟨(∇v0−vnne)⋅n,QhΔu⟩∂T = (Δu,Δv0)T+⟨v0−vb,∇(QhΔu)⋅n⟩∂T−⟨(∇v0−vnne)⋅n,QhΔu⟩∂T,

which implies that

 (Δu,Δv0)T = (ΔωQhu,Δωv)T−⟨v0−vb,∇(QhΔu)⋅n⟩∂T

Next, it follows from the integration by parts that

 (Δu,Δv0)T=(Δ2u,v0)T+⟨Δu,∇v0⋅n⟩∂T−⟨∇(Δu)⋅n,v0⟩∂T.

By summing over all and then using the identity we arrive at

 ∑T∈Th(Δu,Δv0)T = (f,v0)+∑T∈Th⟨Δu,∇v0⋅n−vnne⋅n⟩∂T −∑T∈Th⟨∇(Δu)⋅n,v0−vb⟩∂T,

where we have used the fact that and vanish on the boundary of the domain. Combining the above equation with (4) yields

 (ΔωQhu,Δωv)h = (f,v0)+∑T∈Th⟨Δu−QhΔu,(∇v0−vnne)⋅n⟩∂T −∑T∈Th⟨∇(Δu−QhΔu)⋅n,v0−vb⟩∂T.

Adding to both sides of the above equation gives

 (ΔωQhu,Δωv)h+s(Qhu,v) = (f,v0)+∑T∈Th⟨Δu−QhΔu,(∇v0−vnne)⋅n⟩∂T −∑T∈Th⟨∇(Δu−QhΔu)⋅n,v0−vb⟩∂T+s(Qhu,v).

Subtracting (3.5) from (4) leads to the following error equation

 (Δωeh,Δωv)h+s(eh,v) = ∑T∈Th⟨Δu−QhΔu,(∇v0−vnne)⋅n⟩∂T −∑T∈Th⟨∇(Δu−QhΔu)⋅n,v0−vb⟩∂T+s(Qhu,v)

for all . This completes the derivation of (4.1).

The following Theorem presents an optimal order error estimate for the error function in the trip-bar norm. We believe this tripe-bar norm provides a discrete analogue of the usual -norm.

{theorem}

Let be the weak Galerkin finite element solution arising from (3.5) with finite element functions of order . Assume that the exact solution of (1.1)-(1.3) is sufficiently regular such that . Then, there exists a constant such that

 (4.5) |||uh−Qhu|||≤Chk−1 ∥u∥k+2.

The above estimate is of optimal order in terms of the meshsize , but not in the regularity assumption on the exact solution of the biharmonic equation.

{proof}

By letting in the error equation (4.1), we have

 (4.6)

where

 ℓ(eh) = ∑T∈Th⟨Δu−QhΔu,(∇e0−enne)⋅n⟩∂T −∑T∈Th⟨∇(Δu−QhΔu)⋅n,e0−eb⟩∂T +∑T∈Thh−1T⟨∇Q0u⋅ne−Qb(∇u⋅ne),∇e0⋅ne−en⟩∂T +∑T∈Thh−3T⟨QbQ0u−Qbu,Qbe0−eb⟩∂T.

The rest of the proof shall estimate each of the terms on the right-hand side of (4). For the first term, we use the Cauchy-Schwarz inequality and the estimates (A.5) and (A.6) in Lemma A.1 (see Appendix A) with to obtain

 ∣∣ ∣∣∑T∈Th⟨Δu−QhΔu,(∇e0−enne)⋅n⟩∂T∣∣ ∣∣ ≤ ⎛⎝∑T∈ThhT∥Δu−QhΔu∥2∂T⎞⎠12⎛⎝∑T∈Thh−1T∥∇e0⋅ne−en∥2∂T⎞⎠12 ≤ Chk−1∥u∥k+1|||eh|||.

For the second term, using Lemma A.1, A.2 and A.2 we obtain

 ∣∣ ∣∣∑T∈Th⟨∇(Δu−QhΔu)⋅n,e0−eb⟩∂T∣∣ ∣∣ ≤ ∣∣ ∣∣∑T∈Th⟨∇(Δu−QhΔu)⋅n,Qbe0−eb⟩∂T∣∣ ∣∣ +∣∣ ∣∣∑T∈Th⟨∇(Δu−QhΔu)⋅n,e0−Qbe0⟩∂T∣∣ ∣∣ = ∣∣ ∣∣∑T∈Th⟨∇(Δu−QhΔu)⋅n,Qbe0−eb⟩∂T∣∣ ∣∣ +∣∣ ∣∣∑T∈Th⟨(∇(Δu)−Qb(∇(Δu)))⋅n,e0−Qbe0⟩∂T∣∣ ∣∣ ≤ ⎛⎝∑T∈Thh3T∥∇(Δu−QhΔu)∥2∂T⎞⎠12⋅⎛⎝∑T∈Thh−3T∥Qbe0−eb∥2∂T⎞⎠12 +⎛⎝∑T∈Th∥∇(Δu)−Qb(∇(Δu))∥2∂T⎞⎠12⋅⎛⎝∑T∈Th∥e0−Qbe0∥2∂T⎞⎠12 ≤ Chk−1 ∥u∥k+2|||eh|||,

where the -norm of is used because the estimate in Lemma A.2 is not optimal in terms of the mesh parameter .

The third and fourth terms can be estimated by using the Cauchy-Schwarz inequality and the estimates (A.7) and (A.8) in Lemma A.1 as follows

 (4.10) ∣∣ ∣∣∑T∈Thh−1T⟨∇Q0u⋅ne−Qb(∇u⋅ne),∇e0⋅ne−en⟩∂T∣∣ ∣∣≤Chk−1∥u∥k+1|||eh|||

and

 (4.11) ∣∣ ∣∣∑T∈Thh−3T⟨QbQ0u−Qbu,Qbe0−eb⟩∂T∣∣ ∣∣≤Chk−1∥u∥k+1|||eh|||.

Substituting (4)-(4.11) into (4.6) gives

 |||eh|||2≤Chk−1 ∥u∥k+2|||eh|||,

which implies (4.5) and hence completes the proof.

## 5 Error Estimates in L2

In this section, we shall establish some error estimates for all three components of the error function in the standard norm.

First of all, let us derive an error estimate for the first component of the error function by applying the usual duality argument in the finite element analysis. To this end, we consider the problem of seeking such that

 (5.1) Δ2φ = e0,in Ω, φ = 0, on ∂Ω, ∂φ∂n = 0, on ∂Ω.

Assume that the dual problem has the regularity property in the sense that the solution function and there exists a constant such that

 (5.2) ∥φ∥4≤C∥e0∥.
{theorem}

Let be the weak Galerkin finite element solution arising from (3.5) with finite element functions of order . Let . Assume that the exact solution of (1.1)-(1.3) is sufficiently regular such that and the dual problem (5.1) has the regularity. Then, there exists a constant such that

 (5.3) ∥u0−Q0u∥≤Chk+k0−2∥u∥k+1,

which means we have a sub-optimal order of convergence for and optimal order of convergence for .

{proof}

Testing (5.1) by error function and then using the integration by parts gives

 ∥e0∥2 = (Δ2φ,e0) = ∑T∈T