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


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.


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


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

1 Introduction

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


with clamped boundary conditions


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 [30], orthogonal spline collocation method [8] and the fast multipole methods [24], 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

such that


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 [3]) 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 [29] is a well known nonconforming element for the biharmonic equation for its simplicity. A interior penalty method was developed in [10, 22]. In [31], 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 [41] for a model second order elliptic problem, and this method was further developed in [42, 32, 43]. In [34], 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 [35] 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 [34].

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.,

The norm in is defined by

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 [34].

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,


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


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

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

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

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 .


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


From the integration by parts, we have

Substituting the above identity into (2.3) yields


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 [42]. 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


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


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

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

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,


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

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


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


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

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

Weak Galerkin Algorithm 1

Find satisfying and on and the following equation:


For any , let be given by


Then, defines a norm in the linear space .


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

where we have used the fact that and

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


By letting on each element and summing over all we obtain


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

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

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

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


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


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


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

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.


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

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


for all . Here


Using (2.4) with we obtain

which implies that

Next, it follows from the integration by parts that

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

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

Adding to both sides of the above equation gives

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

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.


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


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.


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



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

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

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




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

which implies (4.5) and hence completes the proof.

5 Error Estimates in

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


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


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


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


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