A conforming DG method for the biharmonic equation on polytopal meshes

TL;DR

This paper introduces a simple, conforming discontinuous Galerkin method for the biharmonic equation on polytopal meshes, providing optimal error estimates and confirming convergence through numerical experiments.

Contribution

It presents a novel, easy-to-implement conforming DG method for biharmonic problems with proven optimal error estimates.

Findings

Optimal order error estimates in discrete $H^2$ norm.

Sub-optimal $L^2$ error estimates for lowest order elements.

Numerical results confirm theoretical convergence.

Abstract

A conforming discontinuous Galerkin finite element method is introduced for solving the biharmonic equation. This method, by its name, uses discontinuous approximations and keeps simple formulation of the conforming finite element method at the same time. The ultra simple formulation of the method will reduce programming complexity in practice. Optimal order error estimates in a discrete $H^2$ norm is established for the corresponding finite element solutions. Error estimates in the $L^2$ norm are also derived with a sub-optimal order of convergence for the lowest order element and an optimal order of convergence for all high order of elements. Numerical results are presented to confirm the theory of convergence.

Tables

Table 1. Table 1: The error and the order of convergence for ( 45 ) on triangular grids (Figure 1 )

level	${\Vert u_h-u\Vert }_0$	rate	${\left|u_h-u\right|}_{1,h}$	rate	$\Vert \left|u_h-u\right|\Vert$	rate
	by the $P_2$ conforming DG finite element
4	0.3653E-03	2.0	0.3281E-02	1.9	0.1229E+01	0.9
5	0.9566E-04	1.9	0.8733E-03	1.9	0.6312E+00	1.0
6	0.2480E-04	1.9	0.2268E-03	1.9	0.3199E+00	1.0
	by the $P_3$ conforming DG finite element
2	0.2291E-03	4.4	0.3275E-02	3.1	0.1612E+00	2.0
3	0.1143E-04	4.3	0.3889E-03	3.1	0.4577E-01	1.8
4	0.7148E-06	4.0	0.4743E-04	3.0	0.1243E-01	1.9

Table 2. Table 2: Error profiles and convergence rates for ( 45 ) on polygonal grids (Figure 2 )

level	${\Vert u_h-u\Vert }_0$	rate	${\left|u_h-u\right|}_{1,h}$	rate	$\Vert \left|u_h-u\right|\Vert$	rate
	by the $P_2$ conforming DG finite element
4	0.3171E-03	1.9	0.4537E-02	1.9	0.3286E+01	1.0
5	0.8671E-04	1.9	0.1175E-02	1.9	0.1647E+01	1.0
6	0.2428E-04	1.8	0.3023E-03	2.0	0.8243E+00	1.0
	by the $P_3$ conforming DG finite element
1	0.3402E-02	0.0	0.3868E-01	0.0	0.3702E+01	0.0
2	0.2027E-03	4.1	0.4895E-02	3.0	0.9408E+00	2.0
3	0.1476E-04	3.8	0.6244E-03	3.0	0.2368E+00	2.0

Equations

$${\Delta }^{2}u$$

$$u$$

$$\frac{\partial u}{\partial n}$$

$$(\Delta u,\Delta v)=(f,v)\forall v\in H_0^2(\Omega ).$$

$$(\Delta u_h,\Delta v)=(f,v)\forall v\in V_h.$$

$$(\Delta u_h,\Delta v{)}_{{\mathcal{T}}_h}$$

$$({\Delta }_w{u}_h,{\Delta }_{w}v)=(f,v),$$

$$(v,w{)}_{{\mathcal{T}}_h}$$

$$⟨v,w{⟩}_{\partial {\mathcal{T}}_h}$$

$$V_h=\left\{v\in L^2(\Omega ):v{|}_T\in P_k(T)T\in {\mathcal{T}}_h\right\}.$$

$$[v]=v{|}_{T_1}{\mathbf{n}}_1+v{|}_{T_2}{\mathbf{n}}_2,[\mathbf{q}]=\mathbf{q}{|}_{T_1}\cdot {\mathbf{n}}_1+\mathbf{q}{|}_{T_2}\cdot {\mathbf{n}}_2,$$

$$\{v\}=\frac{1}{2}(v{|}_{T_1}+v{|}_{T_2})\{\mathbf{q}\}=\frac{1}{2}(\mathbf{q}{|}_{T_1}+\mathbf{q}{|}_{T_2}).$$

$$\{v\}=0,\{\mathbf{q}\}=0,[v]=v\mathbf{n},[\mathbf{q}]=\mathbf{q}\cdot \mathbf{n}.$$

$$({\Delta }_w{u}_h,{\Delta }_{w}v{)}_{{\mathcal{T}}_h}=(f,v)\forall v\in V_h.$$

$$({\Delta }_{w}v,\phi {)}_T=(v,\Delta \phi {)}_T-⟨\{v\},\nabla \phi \cdot \mathbf{n}{⟩}_{\partial T}+⟨\{\nabla v\}\cdot \mathbf{n},\phi {⟩}_{\partial T},\forall \phi \in P_j(T).$$

$${\Delta }_w\varphi ={\mathbb{Q}}_h(\Delta \varphi ),$$

$$({\Delta }_w\varphi ,\tau {)}_T$$

$${\Delta }_w\varphi ={\mathbb{Q}}_h(\Delta \varphi ).$$

$$|||v{|||}^2=({\Delta }_{w}v,{\Delta }_{w}v{)}_{{\mathcal{T}}_h}.$$

$$\Vert v{\Vert }_{2,h}={\left(\sum _{T\in {\mathcal{T}}_h}\left(\Vert \Delta v{\Vert }_T^2+h_T^{-3}\Vert [v]{\Vert }_{\partial T}^2+h_T^{-1}\Vert [\nabla v]{\Vert }_{\partial T}^2\right)\right)}^{\frac{1}{2}}.$$

$$\Vert \phi {\Vert }_e^2\le C\left(h_T^{-1}\Vert \phi {\Vert }_T^2+h_T\Vert \nabla \phi {\Vert }_T^2\right).$$

$$C_1\Vert v{\Vert }_{2,h}\le |||v|||\le C_2\Vert v{\Vert }_{2,h}.$$

$$({\Delta }_{w}v,\phi {)}_T$$

$$\Vert {\Delta }_{w}v{\Vert }_T^2$$

$$\Vert v-\{v\}{\Vert }_e=\Vert [v]{\Vert }_e\mathrm{if}e\subset \partial \Omega ,\Vert v-\{v\}{\Vert }_e=\frac{1}{2}\Vert [v]{\Vert }_e\mathrm{if}e\in {\mathcal{E}}_h^0.$$

$$\begin{array}{rl}& \Vert (\nabla v-\{\nabla v\})\cdot \mathbf{n}{\Vert }_e=\Vert [\nabla v]{\Vert }_e\mathrm{if}e\subset \partial \Omega ,\\ & |(\nabla v-\{\nabla v\})\cdot \mathbf{n}{\Vert }_e=\frac{1}{2}\Vert [\nabla v]{\Vert }_e\mathrm{if}e\in {\mathcal{E}}_h^0.\end{array}$$

$$\Vert {\Delta }_{w}v{\Vert }_T^2$$

$$\Vert {\Delta }_{w}v{\Vert }_T\le C\left(\Vert \Delta v{\Vert }_T+h_T^{-3/2}\Vert [v]{\Vert }_{\partial T}+h_T^{-1/2}\Vert [\nabla v]{\Vert }_{\partial T}\right),$$

$$|||v|||\le C_2\Vert v{\Vert }_{2,h}.$$

$$\sum h_T^{-3}\Vert [v]{\Vert }_{\partial T}^2\le C|||v|||.$$

$$({\Delta }_{w}v,\phi {)}_T$$

$$(\Delta v,{\phi }_0{)}_T=0,⟨(\{\nabla v\}-\nabla v)\cdot \mathbf{n},{\phi }_0{⟩}_{\partial T}=0,$$

Full text

1 Introduction

We consider the biharmonic equation of the form

[TABLE]

where $\Omega$ is a bounded polytopal domain in $\mathbb{R}^{d}$.

The weak formulation of the boundary value problem (2) and (3) is seeking $u\in H^{2}_{0}(\Omega)$ satisfying

[TABLE]

The $H^{2}$ conforming finite element method for the problem (1)-(3) keeps the same simple form as in (4): find $u_{h}\in V_{h}\subset H^{2}_{0}(\Omega)$ such that

[TABLE]

However, it is known that $H^{2}$-conforming methods require $C^{1}$-continuous piecewise polynomials on a simplicial meshes, which imposes difficulty in practical computation. Due to the complexity in the construction of $C^{1}$-continuous elements, $H^{2}$-conforming finite element methods are rarely used in practice for solving the biharmonic equation.

An approach of avoiding construction of $H^{2}$-conforming elements is to use discontinuous approximations. Due to the flexibility of discontinuous Galerkin (DG) finite element methods in element constructions and in mesh generations, many finite element methods have been developed using totally discontinuous polynomials. Here we are only interested in interior penalty discontinuous Galerkin (IPDG) methods since the proposed the method shares the same finite element spaces with IPDG method. For the biharmonic equation, interior penalty discontinuous Galerkin finite element methods have been studied in [1, 2, 3, 4, 5, 8]. One obvious disadvantage of discontinuous finite element methods is their rather complicated formulations which are often necessary to guarantee well posedness and convergence of the methods. For example, the symmetric IPDG method for the biharmonic equation with homogenous boundary conditions [1, 2] has the following formulation:

[TABLE]

where $\sigma$ and $\tau$ are two parameters that need to be tuned.

The purpose of this work is to introduce a conforming DG finite element method for the biharmonic equation which has the following ultra simple formulation without any stabilizing/penalty terms and other mixed terms of lower dimension integrations in (6):

[TABLE]

where $\Delta_{w}$ is called weak Laplacian, an approximation of $\Delta$. The formulation (7) can be viewed as a counterpart of (5) for discontinuous approximations. The conforming DG method was first introduced in [12, 13] for second order elliptic equations, which, by name, means the method using the finite element spaces of DG methods and the simple formulations of conforming methods. This new finite element method shares the same finite element space with the IPDG methods but having much simpler formulation. This simple formulation can be obtained by defining weak Laplacian $\Delta_{w}$ appropriately. The idea here is to raise the degree of polynomials used to compute weak Laplacian $\Delta_{w}$. Using higher degree polynomials in computation of weak Laplacian will not change the size, neither the global sparsity of the stiffness matrix. Optimal order error estimates in a discrete $H^{2}$ for $k\geq 2$ and in $L^{2}$ norm for $k>2$ are established for the corresponding finite element solutions. Numerical results are provided to confirm the theories.

2 A Conforming DG Finite Element Method

Let ${\mathcal{T}}_{h}$ be a partition of the domain $\Omega$ consisting of polygons in two dimension or polyhedra in three dimension satisfying a set of conditions defined in [9] and additional conditions specified in [14]. Denote by ${\cal E}_{h}$ the set of all edges or flat faces in ${\cal T}_{h}$, and let ${\cal E}_{h}^{0}={\cal E}_{h}\backslash\partial\Omega$ be the set of all interior edges or flat faces.

For simplicity, we adopt the following notations,

[TABLE]

Let $P_{k}(K)$ consist all the polynomials degree less or equal to $k$ defined on $K$.

We define a finite element space $V_{h}$ for $k\geq 2$ as follows

[TABLE]

Let $T_{1}$ and $T_{2}$ be two polygons/polyhedrons sharing $e$ if $e\in{\mathcal{E}}_{h}^{0}$. Let $v$ and ${\bf q}$ be scalar and vector valued functions, the jumps $[v]$ and $[{\bf q}]$ are defined as

[TABLE]

and the averages $\{v\}$ and $\{{\bf q}\}$ are defined as

[TABLE]

If $e$ is on $\partial\Omega$, then

[TABLE]

The new conforming DG finite element method for the biharmonic equation (1)-(3) is defined as follows.

Weak Galerkin Algorithm 1**.**

A numerical approximation for (1)-(3) can be obtained by seeking $u_{h}\in V_{h}$ satisfying the following equation:

[TABLE]

Next we will discuss how to compute the weak Laplacian $\Delta_{w}u_{h}$ and $\Delta_{w}v$ in (12). The concept of weak derivative was first introduced in [10, 9] for weak functions in weak Galerkin methods and was modified in [7, 11]. A weak Laplacian operator, denoted by $\Delta_{w}$, is defined as the unique polynomial $\Delta_{w}v\in P_{j}(T)$ for $j>k$ that satisfies the following equation

[TABLE]

Lemma 1**.**

Let $\phi\in H^{2}(\Omega)$, then on any $T\in{\mathcal{T}}_{h}$,

[TABLE]

where $\mathbb{Q}_{h}$ is a locally defined $L^{2}$ projections onto $P_{j}(T)$ on each element $T\in{\mathcal{T}}_{h}$.

Proof.

It is not hard to see that for any $\tau\in P_{j}(T)$ we have

[TABLE]

which implies

[TABLE]

We complete the proof. ∎

3 Well Posedness

First we define a semi-norm ${|||}\cdot{|||}$ as

[TABLE]

Then we introduce a discrete $H^{2}$ norm as follows:

[TABLE]

The following lemma indicates that the two norms $\|\cdot\|_{2,h}$ and ${|||}\cdot{|||}$ are equivalent.

First we need the following trace inequality. For any function $\varphi\in H^{1}(T)$, the trace inequality holds true (see [9] for details):

[TABLE]

Lemma 2**.**

There exist two positive constants $C_{1}$ and $C_{2}$ such that for any $v\in V_{h}$, we have

[TABLE]

Proof.

For any $v\in V_{h}$, it follows from the definition of weak Laplacian (13) and integration by parts that

[TABLE]

By letting $\varphi=\Delta_{w}v$ in (20) we arrive at

[TABLE]

It is easy to see that the following equations hold true for $v$ on $T$ with $e\subset{\partial T}$,

[TABLE]

and

[TABLE]

From the trace inequality (18), (21)-(22) and the inverse inequality we have

[TABLE]

which implies

[TABLE]

and consequently

[TABLE]

Next we will prove

[TABLE]

It follows from (20) that for any $\varphi\in P_{j}(T)$,

[TABLE]

By Lemma 3.1 in [14], there exist a $\varphi_{0}$ such that for $e\subset{\partial T}$,

[TABLE]

and

[TABLE]

Letting $\varphi=\varphi_{0}$ in (23) yields

[TABLE]

which implies

[TABLE]

Taking the summation of the above equation over $T\in{\mathcal{T}}_{h}$ and using (21), one has

[TABLE]

Similarly, by Lemma 3.2 in [14], we can have

[TABLE]

Finally, by letting $\varphi=\Delta v$ in (23) we arrive at

[TABLE]

Using the trace inequality (18), the inverse inequality and (24)-(25), one has

[TABLE]

which gives

[TABLE]

We complete the proof.

∎

Lemma 3**.**

The finite element scheme (12) has a unique solution.

Proof.

It suffices to show that the solution of (12) is trivial if $f=g=\phi=0$. It follows that

[TABLE]

Then the norm equivalence (19) implies $\|u_{h}\|_{2,h}=0$, i.e.

[TABLE]

Therefore, $u_{h}$ is a smooth harmonic function on $\Omega$ and $u_{h}=0$ on $\partial\Omega$. Thus we have $u_{h}=0$, which completes the proof. ∎

4 An Error Equation

Let $e_{h}=u-u_{h}$. Next we derive an error equation that $e_{h}$ satisfies.

Lemma 4**.**

For any $v\in V_{h}$, we have

[TABLE]

where

[TABLE]

Proof.

Testing (1) by $v\in V_{h}$ and using the fact that $\sum_{T\in{\mathcal{T}}_{h}}\langle\nabla(\Delta u)\cdot{\bf n},\{v\}\rangle_{\partial T}=0$ and $\sum_{T\in{\mathcal{T}}_{h}}\langle\Delta u,\{\nabla v\}\cdot{\bf n}\rangle_{\partial T}=0$ and integration by parts, we arrive at

[TABLE]

Next we investigate the term $(\Delta u,\Delta v)_{{\mathcal{T}}_{h}}$ in the above equation. Using (14), integration by parts and the definition of weak Laplacian (13), we have

[TABLE]

Combining the above two equations gives

[TABLE]

which implies that

[TABLE]

The error equation follows from subtracting (12) from the above equation,

[TABLE]

We have proved the lemma. ∎

5 An Error Estimate in $H^{2}$

We start this section by defining some approximation operator. Let $Q_{h}$ be the element-wise defined $L^{2}$ projection onto $P_{k}(T)$ on each element $T$.

Lemma 5**.**

Let $k\geq 2$ and $w\in H^{\max\{k+1,4\}}(\Omega)$. There exists a constant $C$ such that the following estimates hold true:

[TABLE]

Here $\delta_{i,j}$ is the usual Kronecker’s delta with value $1$ when $i=j$ and value [math] otherwise.

The above lemma can be proved by using the trace inequality (18) and the definition of $\mathbb{Q}_{h}$. The proof can also be found in [6].

Lemma 6**.**

Let $w\in H^{\max\{k+1,4\}}(\Omega)$, and $v\in V_{h}$. There exists a constant $C$ such that

[TABLE]

Proof.

Using the Cauchy-Schwartz inequality, (29), (30), (21), (22) and (19), we have

[TABLE]

and

[TABLE]

We complete the proof. ∎

Lemma 7**.**

Let $w\in H^{\max\{k+1,4\}}(\Omega)$, then

[TABLE]

Proof.

For any $T\in{\mathcal{T}}_{h}$, it follows from (13), integration by parts, (18) and inverse inequality that for $w\in P_{j}(T)$,

[TABLE]

which implies

[TABLE]

Taking the summation over $T\in{\mathcal{T}}_{h}$, we have proved the lemma. ∎

Theorem 8**.**

Let $u_{h}\in V_{h}$ be the finite element solution arising from (12). Assume that the exact solution $u\in H^{\max\{k+1,4\}}(\Omega)$. Then, there exists a constant $C$ such that

[TABLE]

Proof.

Let $\epsilon_{h}=Q_{h}u-u_{h}$. Then it is straightforward to obtain

[TABLE]

We will bound the term $(\Delta_{w}e_{h},\Delta_{w}\epsilon_{h})_{{\mathcal{T}}_{h}}$ on right hand side of (5) first. Letting $v=\epsilon_{h}\in V_{h}$ in (26) and using (31)-(32) and (35), we have

[TABLE]

To bound the second term on right hand side of (5), we have by (35),

[TABLE]

Combining the estimates (38) and (39) with (5), we arrive

[TABLE]

which completes the proof. ∎

6 Error Estimates in $L^{2}$ Norm

In this section, we will obtain an error bound for the finite element solution $u_{h}$ in $L^{2}$ norm.

The dual problem considered has the following form,

[TABLE]

Assume that the $H^{4}$ regularity holds,

[TABLE]

Theorem 9**.**

Let $u_{h}\in V_{h}$ be the finite element solution arising from (12). Assume that the exact solution $u\in H^{k+1}(\Omega)$ and (43) holds true. Then, there exists a constant $C$ such that

[TABLE]

Proof.

Testing (40) by $e_{h}$ and using the fact that $\sum_{T\in{\mathcal{T}}_{h}}\langle\nabla(\Delta w)\cdot{\bf n},\{e_{h}\}\rangle_{\partial T}=0$ and $\sum_{T\in{\mathcal{T}}_{h}}\langle\Delta w,\{\nabla e_{h}\}\cdot{\bf n}\rangle_{\partial T}=0$ and integration by parts, we arrive at

[TABLE]

It follows from integration by parts, the definition of weak Laplacian (13) and (14),

[TABLE]

Combining the two equations above implies

[TABLE]

By simple manipulation and (26), we have

[TABLE]

Combining the two equations above implies

[TABLE]

Next, we will estimate all the terms on the right hand side of the above equation. Using the Cauchy-Schwartz inequality, (29)-(30), (21) and (33), we have

[TABLE]

Similarly, by the Cauchy-Schwartz inequality, (29)-(30), (22) and (34), we have

[TABLE]

The estimates (36) and (35) give

[TABLE]

To estimate $I_{4}$, we need to bound $\|\Delta e_{h}\|_{T}$. By (19), (35), (36) and the definition of $Q_{h}$, we have

[TABLE]

The above estimate and the definition of $\mathbb{Q}_{h}$ imply

[TABLE]

Using the Cauchy-Schwartz inequality, (21), (29), (19), (35) and (36), we have

[TABLE]

Similarly, we obtain

[TABLE]

Combining all the estimates above yields

[TABLE]

It follows from the above inequality and the regularity assumption (43).

[TABLE]

We have completed the proof. ∎

7 Numerical Experiments

We solve the following biharmonic equation by $P_{k}$ conforming DG finite element methods:

[TABLE]

with the boundary conditions $u=g_{1}$ and $\nabla u\cdot{\mathbf{n}}=g_{2}$ on $\partial\Omega$. We choose $f$, $g_{1}$ and $g_{2}$ so that the exact solution is

[TABLE]

In the first computation, the first three levels of grids are plotted in Figure 1. The error and the order of convergence for the method are listed in Tables 1. Here on triangular grids, we compute the weak Laplacian $\Delta_{w}v$ by $P_{k+2}$ polynomials. The numerical results confirm the convergence theory.

In the next computation, we use a family of polygonal grids (with pentagons and of 8-side polygons) shown in Figure 2. We let the polynomial degree $j=k+3$ for the weak Laplacian on such polygonal meshes. The rate of convergence is listed in Table 2. The convergence history confirms the theory.