A mixed finite element method for a biharmonic problem with weakly imposed Dirichlet boundary condition

Bishnu P. Lamichhane School of Information and Physical Sciences, University of Newcastle, NSW 2308, Callaghan, [email protected]

Abstract

We consider a mixed finite element method for a biharmonic equation with clamped boundary conditions based on biorthogonal systems with weakly imposed Dirichlet boundary condition. We show that the weak imposition of the boundary condition arising from a natural minimisation formulation allows to get an optimal a priori error estimate for the finite element scheme improving the existing error estimate for such a formulation without weakly imposed Dirichlet boundary condition. We also briefly outline the algebraic formulation arising from the finite element method.

keywords:

Biharmonic problem, mixed finite elements, biorthogonal system, weak Dirichlet boundary condition, Nitsche approach

AMS:

65N30, 65N15

1 Introduction

Thin plates and beams, strain gradient elasticity, phase separation of a binary mixture and fluid flow problems are often modelled by fourth order elliptic and parabolic problems [7, 11, 15, 30]. This difficulty of constructing $H^{2}$ - conforming finite element spaces is avoided either by using a discontinuous Galerkin method as in [11, 6, 30] or by using a mixed formulation as in [9, 8, 12, 7, 13, 2, 26, 10].

In this paper, we start with a mixed finite method due to Ciarlet and Raviart [9, 8, 7] using different spaces for the stream function and vorticity for a fourth order problem with clamped boundary conditions. The great advantage of this formulation is that it allows the use of the standard $H^{1}$ -conforming finite element method. Working with this formulation for clamped boundary conditions the a priori error estimate is sub-optimal [9, 7, 27, 15, 10, 18, 31], where the finite element method of order $k$ converges with $h^{k-\frac{1}{2}}$ in the energy norm. The strong imposition of the Dirichlet boundary condition is the main reason for the sub-optimal convergence rate. In order to get an optimal estimate, we impose the Dirichlet boundary condition weakly using a Nitsche type approach. This leads to an optimal order of convergence improving the existing a priori error estimate for the biharmonic problem with clamped boundary conditions. As in [18] we work with discrete spaces having local basis functions satisfying the condition of biorthogonality for the discretisation of the stream function and vorticity. This yields a very efficient finite element method to approximate the solution of a fourth order problem. While the standard symmetric Nitsche apporach requires a penalty parameter [28], our approach does not require a penalty parameter.

The structure of the rest of the paper is organised as follows. In the rest of this section, we briefly recall a mixed formulation for a biharmonic equation with clamped boundary conditions and extend the formulation to include non-homogeneous clamped boundary conditions. Section 2 is devoted for the numerical analysis of the approach. We give an algebraic formulation of the finite element scheme in Section 3. Finally, we draw a conclusion in the last section.

1.1 Mixed formulation

We now derive a mixed formulation of a fourth order problem. We first briefly recall a mixed formulation of the biharmonic problem with homogeneous clamped boundary conditions.

Homogeneous clamped boundary conditions

Let $\Omega\subset\mbox{$\mathbb{R}$}^{2}$ be a bounded convex domain with polygonal boundary $\Gamma=\partial\Omega$ and outward pointing normal $n$ on $\Gamma$ . We consider the biharmonic equation

\Delta^{2}u=f\quad\text{in}\quad\Omega

(1.1)

with clamped boundary conditions

u=\frac{\partial u}{\partial\mbox{\boldmath{$n$}}}=0\quad\text{on}\quad\Gamma.

(1.2)

Following the same approach as in [9, 7, 18] we recast the biharmonic problem as a minimisation problem with a constraint and then reformulate the problem as a three-field formulation. The main idea here is to include the weak form of the Dirichlet boundary condition. We note that the main difficulty to get optimal error estimates using simplicial Lagrange finite element methods for the biharmonic problem is the imposition of the Dirichlet boundary condition on the boundary in the strong sense, which induces a loss of accuracy in the error estimates. To rectify this we propose to impose the Dirichlet boundary condition weakly using a minimisation formulation or equivalently Nitsche approach. In contrast to other Nitsche approaches, we do not require a penalty parameter in our formulation.

We use usual notations for Sobolev spaces as [23, 1, 16, 5]. We consider the following variational form of the biharmonic problem

J(u)=\inf_{v\in H^{2}_{0}(\Omega)}J(v),

(1.3)

with

J(v)=\frac{1}{2}\int_{\Omega}|\Delta v|^{2}\,dx-\int_{\Omega}f\,v\,dx.

(1.4)

Let $H^{*}(\Omega)$ be the dual space of $H^{1}(\Omega)$ . We now introduce a new unknown $\phi=\Delta u$ and write a weak form of this equation as

\int_{\Omega}\phi\mu\,dx-\langle u,\Delta\mu\rangle=0,\quad\mu\in Q,

where $\langle u,\Delta\mu\rangle$ is the duality pairing between the spaces $H^{1}(\Omega)$ and its dual $H^{*}(\Omega)$ , and

Q=\{v\in H^{1}(\Omega):\,\int_{\Omega}v\,dx=0\}.

This is a right choice for the Lagrange multiplier space as

\int_{\Omega}\phi\,dx=0.

Let $V=H^{1}(\Omega)\times L^{2}(\Omega)$ . The variational problem (1.3) can be recast as the minimization problem [7]

\mathcal{J}(u,\phi)=\inf_{(v,\psi)\in\mathcal{V}}\mathcal{J}(v,\psi),

(1.5)

where

	$\displaystyle\mathcal{J}(v,\psi)$	$\displaystyle=$	$\displaystyle\frac{1}{2}\int_{\Omega}\|\psi\|^{2}\,dx+\frac{1}{2}\\|v\\|^{2}_{\frac{1}{2},\Gamma}-\int_{\Omega}f\,v\,dx,$
	$\displaystyle\mathcal{V}$	$\displaystyle=$	$\displaystyle\{(v,\psi)\in V:\;\int_{\Omega}\psi\,q\,dx-\langle v,\Delta q\rangle=0,\;q\in Q\}.$

In the following, the $H^{\frac{1}{2}}(\Gamma)$ inner product is denoted by $\langle\cdot,\cdot\rangle_{\frac{1}{2},\Gamma}$ and $H^{\frac{1}{2}}$ -norm by $\|\cdot\|^{2}_{\frac{1}{2},\Gamma}$ . The dual space of $H^{\frac{1}{2}}(\Gamma)$ is denoted by $H^{-\frac{1}{2}}(\Gamma)$ .

Non-homogeneous boundary conditions

In the following, we consider the biharmonic problem (1.1) with non-homogeneous clamped boundary conditions with $g_{D}\in H^{\frac{1}{2}}(\Gamma),\;g_{N}\in H^{-\frac{1}{2}}(\Gamma)$ . These boundary conditions are as follows:

u=g_{D}\quad\text{and}\quad\frac{\partial u}{\partial\mbox{\boldmath{$n$}}}=g_{N}\quad\text{on}\quad\Gamma.

(1.6)

Then, we have the minimisation problem (1.5) with

	$\displaystyle\mathcal{J}(v,\psi)$	$\displaystyle=$	$\displaystyle\frac{1}{2}\int_{\Omega}\|\psi\|^{2}\,dx+\frac{1}{2}\\|v-g_{D}\\|^{2}_{\frac{1}{2},\Gamma}-\int_{\Omega}f\,v\,dx,$
	$\displaystyle\mathcal{W}$	$\displaystyle=$	$\displaystyle\{(v,\psi)\in V:\;\int_{\Omega}\psi\,q\,dx-\langle v,\Delta q\rangle=\langle g_{N},q\rangle_{\Gamma}-\langle\frac{\partial q}{\partial\mbox{\boldmath{$n$}}},g_{D}\rangle_{\Gamma},\;q\in Q\},$

where $\langle\cdot,\cdot\rangle_{\Gamma}$ is the duality pairing between the spaces $H^{\frac{1}{2}}(\Gamma)$ and its dual $H^{-\frac{1}{2}}(\Gamma)$

Remark 1.

Here, the normal derivative of an $H^{1}$ -function is a generalised normal derivative as defined in [25, 24]. Lemma 4.3 of [24] gives the following bound for the normal derivative of $q\in H^{1}(\Omega)$ (see also [25])

\|\frac{\partial q}{\partial\mbox{\boldmath{$n$}}}\|_{-\frac{1}{2},\Gamma}\leq C(\|q\|_{1,\Omega}+\|\Delta q\|_{H^{*}(\Omega)}).

The problem (1.5) can be recast as a saddle point formulation [18, 9, 7, 10]. The saddle point problem is: Given $\ell\in H^{-1}(\Omega)$ , find $((u,\phi),p)\in V\times Q$ such that

\begin{array}[]{llccc}a((u,\phi),(v,\psi))+&b((v,\psi),p)&=&\ell(v),&\quad(v,\psi)\in V,\\ b((u,\phi),q)&&=&g(q),&\quad q\in Q,\end{array}

(1.7)

where

	$\displaystyle a((u,\phi),(v,\psi))=\int_{\Omega}\phi\psi\,dx+\langle u,v\rangle_{\frac{1}{2},\Gamma},$		(1.8)
	$\displaystyle\ell(v)=\int_{\Omega}fv\,dx+\langle g_{D},v\rangle_{\frac{1}{2},\Gamma},\quad b((v,\psi),q)=\int_{\Omega}\psi\,q\,dx-\langle v,\Delta q\rangle,$
	$\displaystyle\text{and}\quad g(q)=\langle g_{N},q\rangle_{\Gamma}-\langle\frac{\partial q}{\partial\mbox{\boldmath{$n$}}},g_{D}\rangle_{\Gamma}.$

Consistency

Let $u\in H^{2}(\Omega)$ be the solution of the biharmonic problem (1.1) with the non-homogeneous boundary conditions (1.6). Let $\phi=\Delta u$ and $p=-\phi$ . An integration by parts can be performed to show that they satisfy the saddle point equations (1.7).

Remark 2 (Existence and uniqueness of the solution).

There is a difficulty in proving the coercivity of the bilinear form $a(\cdot,\cdot)$ in the saddle point problem (1.7) as the standard trace theorem [16] does not work for the generalised normal derivative [24, 25]. However, there is no problem for defining the standard normal derivative for a function $q_{h}$ in the standard finite element space, see the next section. Therefore, we do not analyse the existence and uniqueness of the saddle point problem (1.7), but rather focus on its discrete counterpart in the following section.

2 Finite element discretizations

We consider a quasi-uniform and shape-regular triangulation $\mathcal{T}_{h}$ of the polygonal domain $\Omega$ with the global mesh-size $h$ , where $\mathcal{T}_{h}$ consists of triangles or parallelograms. Let $\mathcal{C}_{h}$ be the collection of boundary edges of the triangulation of $\Omega$ . We use $h_{K}$ and $h_{e}$ to denote the sizes of the elements in $\mathcal{T}_{h}$ and $\mathcal{C}_{h}$ , respectively. Let $S_{h}\subset H^{1}(\Omega)$ be a standard Lagrange finite element space of order $k\in\mbox{$\mathbb{N}$}$ , and $M_{h}\subset L^{2}(\Omega)$ be another piecewise polynomial space. We also set $V_{h}=S_{h}\times M_{h}$ . We have a well-known approximation result for every $u\in H^{k+1}(\Omega)$ [3]: there exists a function $u_{h}\in S_{h}$ such that

h\|u-u_{h}\|_{1,\Omega}+\|u-u_{h}\|_{0,\Omega}\leq Ch^{k+1}\|u\|_{k+1,\Omega}.

In the following, we use a generic constant $C$ , which takes different values in different occurrences but is always independent of the mesh-size. We impose the following assumptions on $M_{h}$ .

Assumption 3.

We assume that there is a constant $C>0$ independent of the mesh-size such that

\displaystyle\|q_{h}\|_{0,\Omega}\leq C\sup_{\phi_{h}\in S_{h}}\frac{\int_{\Omega}\phi_{h}q_{h}\,dx}{\|\phi_{h}\|_{0,\Omega}},\quad q_{h}\in M_{h},

(2.1)

Assumption 4.

The space $M_{h}$ has the approximation property:

\inf_{\lambda_{h}\in M_{h}}\|\phi-\lambda_{h}\|_{0,\Omega}\leq Ch^{k}|\phi|_{k,\Omega},\quad\phi\in H^{k}(\Omega).

(2.2)

We use

Q_{h}=\{v_{h}\in S_{h}:\,\int_{\Omega}v_{h}\,dx=0\}

to approximate the Lagrange multiplier space $Q$ . Our analysis is based on the following mesh-dependent inner product and the norm induced by this inner product on the boundary of $\Omega$ for $s\in[-1,1]$ [28]:

\langle v,\,w\rangle_{s,h}=\sum_{e\in\mathcal{C}_{h}}\frac{1}{h_{e}^{2s}}\int_{e}v\,w\,d\sigma,\quad v,w\in L^{2}(\Omega).

(2.3)

We will use the mesh-dependent norm for $v_{h}\in S_{h}$ ,

\|v_{h}\|_{1,h}^{2}=\|v_{h}\|_{1,\Omega}^{2}+\|v_{h}\|_{\frac{1}{2},h}^{2},

where $\|\cdot\|_{\frac{1}{2},h}$ is the norm induced by the inner product (2.3). In fact,

\|u_{h}\|_{\frac{1}{2},h}^{2}=\sum_{e\in\mathcal{C}_{h}}\frac{1}{h_{e}}\int_{e}u_{h}^{2}\,d\sigma.

With the definition of $\|\cdot\|_{s,h}$ -norm we have the following Cauchy-Schwarz type inequality for the inner product $\langle\cdot,\cdot\rangle_{\frac{1}{2},h}$ [3.13 of [28]]:

\displaystyle\langle v,\,w\rangle_{\frac{1}{2},h}\leq\|v\|_{\frac{1}{2},h}\|w\|_{-\frac{1}{2},h},\quad v\in H^{1}(\Omega),\;w\in L^{2}(\Omega).

(2.4)

The discrete biharmonic problem is given as a saddle point problem: given $f\in H^{-1}(\Omega)$ , $g_{D}\in H^{\frac{1}{2}}(\Gamma),\;g_{N}\in H^{-\frac{1}{2}}(\Gamma)$ , find $((u_{h},\phi_{h}),p_{h})\in V_{h}\times S_{h}$ such that

\begin{array}[]{llccc}a_{h}((u_{h},\phi_{h}),(v_{h},\psi_{h}))+&b_{h}((v_{h},\psi_{h}),p_{h})&=&\ell_{h}(v_{h}),&\quad(v_{h},\psi_{h})\in V_{h},\\ b_{h}((u_{h},\phi_{h}),q_{h})&&=&g_{h}(q_{h}),&\quad q_{h}\in Q_{h},\end{array}

(2.5)

where

	$\displaystyle a_{h}((u_{h},\phi_{h}),(v_{h},\psi_{h}))$	$\displaystyle=$	$\displaystyle\int_{\Omega}\phi_{h}\psi_{h}\,dx+\langle u_{h},\,v_{h}\rangle_{\frac{1}{2},h},\;b_{h}((v_{h},\psi_{h}),q_{h})=\int_{\Omega}\psi_{h}\,q_{h}\,dx-\langle v_{h},\Delta_{h}q_{h}\rangle$
	$\displaystyle\ell_{h}(v_{h})$	$\displaystyle=$	$\displaystyle\int_{\Omega}fv_{h}\,dx+\langle g_{D},\,v_{h}\rangle_{\frac{1}{2},h}\quad\text{and}\quad g_{h}(q_{h})=\langle g_{N},q_{h}\rangle_{\Gamma}-\int_{\Gamma}\frac{\partial q_{h}}{\partial\mbox{\boldmath{$n$}}}g_{D}\,d\sigma,$

where for $q\in H^{\frac{3}{2}+\epsilon}(\Omega)$ with $\epsilon>0$ , $\Delta_{h}:H^{\frac{3}{2}+\epsilon}(\Omega)\to M_{h}$ is defined as

\langle v_{h},\Delta_{h}q\rangle=-\int_{\Omega}\nabla v_{h}\cdot\nabla q\,dx+\int_{\Gamma}\frac{\partial q}{\partial\mbox{\boldmath{$n$}}}v_{h}\,d\sigma,\quad v_{h}\in S_{h}.

We note that $\Delta_{h}q$ is well-defined due to Assumption 3.

In order to analyse the finite element problem we introduce the mesh-dependent graph norm on $V_{h}$ defined as

\displaystyle\|(v_{h},\psi_{h})\|_{a}=\sqrt{\|\psi_{h}\|^{2}_{0,\Omega}+\|v_{h}\|_{1,h}^{2}}

(2.6)

and the following mesh-dependent norm for the Lagrange multiplier $q_{h}\in Q_{h}$ defined as

\|q_{h}\|^{2}_{Q_{h}}=\|q_{h}\|^{2}_{0,\Omega}+\|\Delta_{h}q_{h}\|_{-1,h}^{2},

where

\|\Delta_{h}q_{h}\|_{-1,h}=\sup_{v_{h}\in S_{h}}\frac{\langle\Delta_{h}q_{h},v_{h}\rangle}{\|v_{h}\|_{1,h}}.

We can see that the continuity of the bilinear form $a_{h}(\cdot,\cdot)$ and linear forms $\ell_{h}(\cdot)$ and $g_{h}(\cdot)$ follows from the Cauchy-Schwarz and trace inequalities [14]. The continuity of the bilinear form $b_{h}(\cdot,\cdot)$ follows from

\|w_{h}\|_{1,h}\|\Delta_{h}q_{h}\|_{-1,h}=\|w_{h}\|_{1,h}\sup_{v_{h}\in S_{h}}\frac{\langle\Delta_{h}q_{h},v_{h}\rangle}{\|v_{h}\|_{1,h}}\geq\left|\langle\Delta_{h}q_{h},w_{h}\rangle\right|,\;w_{h}\in S_{h},\,q_{h}\in Q_{h}.

Thus

|b_{h}((w_{h},\psi_{h}),q_{h})|\leq\|\psi_{h}\|_{0,\Omega}\|q_{h}\|_{0,\Omega}+\|w_{h}\|_{1,h}\|\Delta_{h}q_{h}\|_{-1,h}.

We now show the inf-sup condition for the bilinear form $b_{h}(\cdot,\cdot)$ . We need to show the existence of a mesh-independent constant $C$ such that

\sup_{(v_{h},\psi_{h})\in V_{h}}\frac{b_{h}((v_{h},\psi_{h}),q_{h})}{\|(v_{h},\psi_{h})\|_{a}}\geq C\|q_{h}\|_{Q_{h}}.

(2.7)

First we set $v_{h}=0$ on the left hand side of the above inequality and use (3) to obtain

\sup_{(v_{h},\psi_{h})\in V_{h}}\frac{b_{h}((v_{h},\psi_{h}),q_{h})}{\|(v_{h},\psi_{h})\|_{a}}\geq\sup_{\psi_{h}\in M_{h}}\frac{\int_{\Omega}q_{h}\psi_{h}}{\|\psi_{h}\|_{0,\Omega}}\geq C\|q_{h}\|_{0,\Omega}.

In the second step, we set $\psi_{h}=0$ on the left hand side of the inequality (2.7) and use the definition of the norm $\|\cdot\|_{-1,h}$ to obtain

\sup_{(v_{h},\psi_{h})\in V_{h}}\frac{b_{h}((v_{h},\psi_{h}),q_{h})}{\|(v_{h},\psi_{h})\|_{a}}\geq\sup_{v_{h}\in S_{h}}\frac{\langle v_{h},\Delta_{h}q_{h}\rangle}{\|v_{h}\|_{1,h}}=\|\Delta_{h}q_{h}\|_{-1,h}.

Now we turn our attention to prove the coercivity of the bilinear form $a_{h}(\cdot,\cdot)$ on the kernel space $\mathcal{V}_{h}$ defined as

\mathcal{V}_{h}=\{(v_{h},\psi_{h})\in V_{h}:\;\int_{\Omega}\psi_{h}\,q_{h}\,dx-\langle\Delta_{h}q_{h},v_{h}\rangle=0,\;q_{h}\in Q_{h}\}.

(2.8)

First, we note that

a_{h}((v_{h},\psi_{h}),(v_{h},\psi_{h}))=\|\psi_{h}\|^{2}_{0,\Omega}+\|v_{h}\|^{2}_{\frac{1}{2},h}.

If $(v_{h},\psi_{h})\in\mathcal{V}_{h}$ , we have

\int_{\Omega}\left(\psi_{h}\,q_{h}+\nabla v_{h}\cdot\nabla q_{h}\right)\,dx=\int_{\Gamma}\frac{\partial q_{h}}{\partial\mbox{\boldmath{$n$}}}v_{h}\,d\sigma,\quad q_{h}\in Q_{h}.

(2.9)

Let

q_{h}=v_{h}-\frac{1}{|\Omega|}\int_{\Omega}v_{h}\,dx\in Q_{h}.

Then we have

\frac{\partial q_{h}}{\partial\mbox{\boldmath{$n$}}}=\frac{\partial v_{h}}{\partial\mbox{\boldmath{$n$}}}\quad\text{and}\quad\nabla q_{h}=\nabla v_{h}.

Hence for $(v_{h},\psi_{h})\in\mathcal{V}_{h}$ , using this $q_{h}$ in (2.9), we obtain

\|\nabla v_{h}\|_{0,\Omega}^{2}=\int_{\Gamma}\frac{\partial v_{h}}{\partial\mbox{\boldmath{$n$}}}v_{h}\,d\sigma-\int_{\Omega}\psi_{h}\,\left(v_{h}-\frac{1}{|\Omega|}\int_{\Omega}v_{h}\,dx\right)\,dx.

(2.10)

We now apply the Cauchy-Schwarz type inequality for the boundary integral of the first term on the right of the above equation

\left|\int_{\Gamma}\frac{\partial v_{h}}{\partial\mbox{\boldmath{$n$}}}v_{h}\,d\sigma\right|\leq\left\|\frac{\partial v_{h}}{\partial\mbox{\boldmath{$n$}}}\right\|_{-\frac{1}{2},h}\left\|v_{h}\right\|_{\frac{1}{2},h},

so that (2.10) yields

\|\nabla v_{h}\|_{0,\Omega}^{2}\leq\left\|\frac{\partial v_{h}}{\partial\mbox{\boldmath{$n$}}}\right\|_{-\frac{1}{2},h}\left\|v_{h}\right\|_{\frac{1}{2},h}+\|\psi_{h}\|_{0,\Omega}\left\|v_{h}-\frac{1}{|\Omega|}\int_{\Omega}v_{h}\,dx\right\|_{0,\Omega}.

(2.11)

In terms of the following trace inequality [(4) of [14]]

\left\|\frac{\partial v_{h}}{\partial\mbox{\boldmath{$n$}}}\right\|_{-\frac{1}{2},h}\leq C\|\nabla v_{h}\|_{0,\Omega},

and Poincaré-Friedrichs inequality

\left\|v_{h}-\frac{1}{|\Omega|}\int_{\Omega}v_{h}\,dx\right\|_{0,\Omega}\leq C\|\nabla v_{h}\|_{0,\Omega},

we get from (2.11)

\|\nabla v_{h}\|_{0,\Omega}^{2}\leq C\left(\|\nabla v_{h}\|_{0,\Omega}\left\|v_{h}\right\|_{\frac{1}{2},h}+\|\psi_{h}\|_{0,\Omega}\|\nabla v_{h}\|_{0,\Omega}\right).

Hence we have

\|\nabla v_{h}\|_{0,\Omega}\leq C(\|\psi_{h}\|_{0,\Omega}+\left\|v_{h}\right\|_{\frac{1}{2},h}).

Moreover, we have a mesh-independent constant $C$ such that [4]

\|v_{h}\|_{0,\Omega}\leq C(\|\nabla v_{h}\|_{0,\Omega}+\left\|v_{h}\right\|_{\frac{1}{2},h}).

Thus we have the following lemma for the coercivity of the bilinear form $a_{h}(\cdot,\cdot)$ on $\mathcal{V}_{h}$ .

Lemma 5.

There exists $\alpha_{0}>0$ independent of the mesh-size $h$ such that

a_{h}((v_{h},\psi_{h}),(v_{h},\psi_{h}))\geq\alpha_{0}(\|v_{h}\|^{2}_{1,h}+\|\psi_{h}\|^{2}_{0,\Omega}),\;(v_{h},\psi_{h})\in\mathcal{V}_{h}.

Hence we have obtained the well-posedness of the saddle point problem (2.5).

Lemma 6.

The saddle point problem (2.5) has a unique solution $((u_{h},\phi_{h}),p_{h})\in V_{h}\times S_{h}$ .

We use the following lemma to prove the a priori error estimate for the discrete solution [18].

Lemma 7.

Let $u$ be the solution of the biharmonic equation (1.1) with non-homogeneous boundary condition (1.6), and $\phi=\Delta u$ as well as $p=-\phi$ . Let $p\in H^{k+1}(\Omega)$ . Let $((u_{h},\phi_{h}),p_{h})\in V_{h}\times Q_{h}$ be the solution of the discrete problem (2.5). Then there exists a constant $C>0$ independent of the mesh-size $h$ so that

\|(u-u_{h},\phi-\phi_{h})\|_{a}\leq C\left(\inf_{(w_{h},\xi_{h})\in\mathcal{W}_{h}}\|(u-w_{h},\phi-\xi_{h})\|_{a}+h^{k}\|p\|_{k+1,\Omega}\right),

(2.12)

where

\mathcal{W}_{h}=\{(w_{h},\xi_{h})\in V_{h}|\,\int_{\Omega}\xi_{h}q_{h}\,dx-\langle\Delta_{h}q_{h},w_{h}\rangle=\langle g_{N},q_{h}\rangle_{\Gamma}-\langle\frac{\partial q_{h}}{\partial\mbox{\boldmath{$n$}}},g_{D}\rangle_{\Gamma},\;q_{h}\in Q_{h}\}.

Proof. Let $(w_{h},\xi_{h})\in\mathcal{W}_{h}$ . Then $(w_{h},\xi_{h})$ satisfies

\int_{\Omega}\xi_{h}q_{h}\,dx-\langle\Delta_{h}q_{h},w_{h}\rangle=\langle g_{N},q_{h}\rangle_{\Gamma}-\langle\frac{\partial q_{h}}{\partial\mbox{\boldmath{$n$}}},g_{D}\rangle_{\Gamma},\;q_{h}\in Q_{h}.

Thus (2.5) implies $(u_{h}-w_{h},\phi_{h}-\xi_{h})\in\mathcal{V}_{h}$ , and hence coercivity of $a_{h}(\cdot,\cdot)$ on $\mathcal{V}_{h}$ yields

\alpha_{0}\|(u_{h}-w_{h},\phi_{h}-\xi_{h})\|_{a}\leq\sup_{(v_{h},\psi_{h})\in\mathcal{V}_{h}}\frac{a_{h}((u_{h}-w_{h},\phi_{h}-\xi_{h}),(v_{h},\psi_{h}))}{\|(v_{h},\psi_{h})\|_{a}}.

Since from (2.5) and (1.7) $a_{h}((u-u_{h},\phi-\phi_{h}),(v_{h},\psi_{h}))+b_{h}((v_{h},\psi_{h}),p)=0$ for all $(v_{h},\psi_{h})\in\mathcal{V}_{h}$ , we have

	$\displaystyle a_{h}((u_{h}-w_{h},\phi_{h}-\xi_{h}),(v_{h},\psi_{h}))$	$\displaystyle=$	$\displaystyle a_{h}((u-w_{h},\phi-\xi_{h}),(v_{h},\psi_{h}))+a_{h}((u_{h}-u,\phi_{h}-\phi),(v_{h},\psi_{h}))$
		$\displaystyle=$	$\displaystyle a_{h}((u-w_{h},\phi-\xi_{h}),(v_{h},\psi_{h}))+b_{h}((v_{h},\psi_{h}),p).$

Let $\tilde{p}_{h}\in Q_{h}$ be a finite element interpolant for $p$ . Using the fact that

b_{h}((v_{h},\psi_{h}),p)=\int_{\Omega}\psi_{h}p\,dx+\int_{\Omega}\nabla p\cdot\nabla v_{h}\,dx-\langle\frac{\partial p}{\partial\mbox{\boldmath{$n$}}},v_{h}\rangle_{\Gamma},\;\;\text{and}\;\;(v_{h},\psi_{h})\in\mathcal{V}_{h},

we get

b_{h}((v_{h},\psi_{h}),p)=b_{h}((v_{h},\psi_{h}),p-\tilde{p}_{h})=\int_{\Omega}\psi_{h}(p-\tilde{p}_{h})\,dx+\int_{\Omega}\nabla(p-\tilde{p}_{h})\cdot\nabla v_{h}\,dx-\langle\frac{\partial(p-\tilde{p}_{h})}{\partial\mbox{\boldmath{$n$}}},v_{h}\rangle_{\Gamma}.

We note that the interpolant $\tilde{p}_{h}$ satisfies [29, Lemma 2.3]

\left|\langle\frac{\partial(p-\tilde{p}_{h})}{\partial\mbox{\boldmath{$n$}}},v_{h}\rangle_{\Gamma}\right|\leq h^{k}\|p\|_{k+1,\Omega}\|v_{h}\|_{\frac{1}{2},h}.

And hence

|b_{h}((v_{h},\psi_{h}),p)|\leq Ch^{k}\|p\|_{k+1,\Omega}\,\|(v_{h},\psi_{h})\|_{a}.

Thus

	$\displaystyle\alpha_{0}\\|(u_{h}-w_{h},\phi_{h}-\xi_{h})\\|_{a}$	$\displaystyle\leq$	$\displaystyle\sup_{(v_{h},\psi_{h})\in\mathcal{V}_{h}}\frac{a_{h}((u-w_{h},\phi-\xi_{h}),(v_{h},\psi_{h}))}{\\|(v_{h},\psi_{h})\\|_{a}}+Ch^{k}\\|p\\|_{k+1,\Omega}$
		$\displaystyle\leq$	$\displaystyle\\|(u-w_{h},\phi-\xi_{h})\\|_{a}+Ch^{k}\\|p\\|_{k+1,\Omega},$

where we have used the fact that the continuity constant of the bilinear form $a(\cdot,\cdot)$ is 1. Finally, a triangle inequality yields the estimate (2.12):

	$\displaystyle\\|(u-u_{h},\phi-\phi_{h})\\|_{a}$	$\displaystyle\leq$	$\displaystyle\\|(u-w_{h},\phi-\xi_{h})\\|_{a}+\\|(w_{h}-u_{h},\xi_{h}-\phi_{h})\\|_{a}$
		$\displaystyle\leq$	$\displaystyle\left(1+\frac{1}{\alpha_{0}}\right)\\|(u-w_{h},\phi-\xi_{h})\\|_{a}+\frac{C}{\alpha_{0}}h^{k}\\|p\\|_{k+1,\Omega}.$

Theorem 8.

Let $u$ be the solution of the biharmonic equation (1.1) with non-homogeneous boundary condition (1.6), and $\phi=\Delta u$ as well as $p=-\phi$ . Let $((u_{h},\phi_{h}),p_{h})\in V_{h}\times Q_{h}$ be the solution of the discrete saddle point problem (2.5). Let $u\in H^{k+1}(\Omega)\cap H^{1}_{0}(\Omega)$ , $\phi\in H^{k}(\Omega),\,p\in H^{k+1}(\Omega)$ , and Assumptions (3) and (4) are satisfied. Then there exists a constant $C>0$ independent of the mesh-size $h$ so that

\|(u-u_{h},\phi-\phi_{h})\|_{a}\leq Ch^{k}\left(\|u\|_{k+1,\Omega}+|\phi|_{k,\Omega}+\|p\|_{k+1,\Omega}\right).

(2.13)

Proof. Let $\Pi_{h}:L^{2}(\Omega)\to M_{h}$ and $\Pi^{*}_{h}:L^{2}(\Omega)\to S_{h}$ be two projections defined by

\int_{\Omega}\Pi_{h}v\,q_{h}\,dx=\int_{\Omega}v\,q_{h}\,dx,\;q_{h}\in S_{h},\quad\text{and}

\int_{\Omega}\Pi^{*}_{h}v\,\eta_{h}\,dx=\int_{\Omega}v\,\eta_{h}\,dx,\;\eta_{h}\in M_{h}.

These projectors are well-defined by Assumption 3. Moreover, using Assumptions 3 and 4 we have [20]

\|\Pi_{h}v\|_{0,\Omega}\leq C\|v\|_{0,\Omega},\;\text{and}\;\|\Pi_{h}w-w\|_{0,\Omega}\leq Ch^{k}\|w\|_{k,\Omega}\;\text{for}\;v\in L^{2}(\Omega),\;\text{and}\;w\in H^{k}(\Omega).

(2.14)

Similarly, for $v\in L^{2}(\Omega)$ and $w\in H^{1}(\Omega)$ , we have [20]

\|\Pi^{*}_{h}v\|_{0,\Omega}\leq C\|v\|_{0,\Omega},\quad\text{and}\quad\|\Pi^{*}_{h}w\|_{1,\Omega}\leq C\|w\|_{1,\Omega}.

(2.15)

We also have for $r=\{0,1\}$ and $w\in H^{k+1}(\Omega)$

\|\Pi^{*}_{h}w-w\|_{r,\Omega}\leq Ch^{k+1-r}\|w\|_{k+1,\Omega}.

(2.16)

Moreover, for $w\in H^{k+1}(\Omega)$ , for the projector $\Pi_{h}^{*}$ , we have [Lemma 1 of [28]]

\|w-\Pi^{*}_{h}w\|_{1,h}\leq Ch^{k}\|w\|_{k+1,\Omega}.

(2.17)

For the exact solution $\phi=\Delta u$ , we get

\int_{\Omega}\phi\,q_{h}\,dx-\langle\Delta_{h}q_{h},u\rangle=\langle\frac{\partial q_{h}}{\partial\mbox{\boldmath{$n$}}},g_{D}\rangle_{\Gamma}+\langle g_{N},q_{h}\rangle_{\Gamma},\;q_{h}\in Q_{h}.

(2.18)

Since $\Delta_{h}q_{h}\in M_{h}$ , we have

\langle\Delta_{h}q_{h},\Pi^{*}_{h}u\rangle=\int_{\Omega}\Delta_{h}q_{h}\,u\,dx.

Thus we have

\int_{\Omega}\Pi_{h}\phi\,q_{h}\,dx-\langle\Delta_{h}q_{h},\Pi^{*}_{h}u\rangle=\langle\frac{\partial q_{h}}{\partial\mbox{\boldmath{$n$}}},g_{D}\rangle_{\Gamma}+\langle g_{N},q_{h}\rangle_{\Gamma},\;q_{h}\in Q_{h}.

(2.19)

Hence we have obtained that $(\Pi^{*}_{h}u,\Pi_{h}\phi)\in\mathcal{W}_{h}$ , and

\|(u-\Pi^{*}_{h}u,\phi-\Pi_{h}\phi)\|_{a}\leq Ch^{k}\left(\|u\|_{k+1,\Omega}+|\phi|_{k,\Omega}\right).

The proof now follows from Lemma 7.

Remark 9.

The existing error estimate approaches require an extra regularity of the solution $u$ [22, 18]. The energy error estimate in [18, 10] is sub-optimal even with the extra regularity, whereas the error estimate in [22] is optimal but the approach works only on rectangular meshes with a special structure.

3 Algebraic formulation

To obtain an efficient numerical scheme in which all the auxiliary variables (the vorticity $\phi_{h}$ and the Lagrange multiplier $p_{h}$ ) can be statically condensed out from the system, we construct a biorthogonal system for the sets of basis functions of $Q_{h}$ and $M_{h}$ . Let $\{\varphi_{1},\cdots,\varphi_{n}\}$ be a finite element basis for the space $Q_{h}$ . A finite element basis $\{\mu_{1},\cdots,\mu_{n}\}$ for the space $M_{h}$ with $\mathop{\operator@font supp}\nolimits\mu_{i}=\mathop{\operator@font supp}\nolimits\varphi_{i}$ , $1\leq i\leq n$ , is constructed in such a way that the basis functions of $Q_{h}$ and $M_{h}$ satisfy a condition of biorthogonality relation

\displaystyle\int_{\Omega}\mu_{i}\ \varphi_{j}\,dx=c_{j}\delta_{ij},\;c_{j}\neq 0,\;1\leq i,j\leq n,

(3.1)

where $n:=\dim M_{h}=\dim Q_{h}$ , $\delta_{ij}$ is the Kronecker symbol, and $c_{j}$ a scaling factor proportional to the area $|\mathop{\operator@font supp}\nolimits\phi_{j}|$ . The basis functions of $M_{h}$ are constructed in a reference element and they satisfy (3), (4) and (3.1) [19, 17, 21].

Let $\vec{u}$ , $\vec{\phi}$ and $\vec{p}$ be the vector representations of the solution $u_{h}$ , $\phi_{h}$ and $p_{h}$ , respectively. Let $\mathtt{A}\vec{u}$ , $\mathtt{M}\vec{\phi}$ and $\mathtt{D}\vec{\phi}$ be algebraic representations of the bilinear forms $\int_{\Omega}u_{h}\Delta_{h}q_{h}\,dx$ , $\int_{\Omega}\phi_{h}\psi_{h}\,dx$ and $\int_{\Omega}\phi_{h}q_{h}\,dx$ , respectively, where $u_{h}\in S_{h}$ , $q_{h}\in Q_{h}$ , $\phi_{h},\psi_{h}\in M_{h}$ . We also denote the algebraic representation of the bilinear form $\langle u_{h},v_{h}\rangle_{\frac{1}{2},h}$ by $\mathtt{B}_{\Gamma}\vec{u}$ . Although the bilinear form $\langle u_{h},v_{h}\rangle_{\frac{1}{2},h}$ is restricted to the boundary $\Gamma$ of the domain $\Omega$ , $\mathtt{B}_{\Gamma}$ is the extended form of the algebraic representation so that the number of columns of the matrix $\mathtt{B}_{\Gamma}$ is equal to the number of components in $\vec{u}$ , where entries of the matrix $\mathtt{B}_{\Gamma}$ corresponding to interior nodes of the mesh are all set to zero. Then the algebraic formulation of the saddle point problem (2.5) is given by

\left[\begin{array}[]{cccc}\mathtt{B}_{\Gamma}&0&-\mathtt{A}^{T}\\ 0&\mathtt{M}&\mathtt{D}\\ -\mathtt{A}&\mathtt{D}&0\end{array}\right]\left[\begin{array}[]{ccc}\vec{u}\\ \vec{\phi}\\ \vec{p}\end{array}\right]=\left[\begin{array}[]{ccc}\vec{f}\\ 0\\ \vec{g}\end{array}\right],

(3.2)

where $\vec{f}$ is the vector associated with the linear form $\ell_{h}(v_{h})$ , and $\vec{g}$ is the vector representation of $g_{h}(q_{h})$ . Since the matrix $\mathtt{D}$ is diagonal, we can do the static condensation of unknowns $\vec{\phi}$ and $\vec{p}$ and arrive at the following linear system based on the unknown $\vec{u}$ associated only with the stream function:

\left(\mathtt{M}_{\Gamma}+\mathtt{A}^{T}\mathtt{D}^{-1}\mathtt{M}\mathtt{D}^{-1}\mathtt{A}\right)\vec{u}=(\vec{f}-(\mathtt{A}^{T}\mathtt{D}^{-1}\mathtt{M}\mathtt{D}^{-1})\vec{g}).

(3.3)

Since the inverse of the matrix $\mathtt{D}$ is diagonal, the system matrix in (3.3) is sparse. It is important to have the system matrix to have sparse structure if an iterative solver is to be applied. The vector corresponding to the vorticity $\vec{\phi}$ and the Lagrange multiplier $\vec{p}$ can be computed by simply inverting the diagonal matrix using the second and third blocks of (3.2).

4 Conclusion

We have proposed a finite element formulation for the biharmonic equation with clamped boundary conditions leading to an optimal convergence rate improving the existing a priori error estimate in the energy norm. The main idea is to impose the Dirichlet boundary condition weakly using the Nitsche technique. The new formulation also allows to use a biorthogonal system that gives an efficient finite element approach. In contrast to other Nitsche approaches, we do not require a penalty parameter in our formulation.

Acknowledgement

Part of this work was completed during my visit to the Indian Institute of Technology, Mumbai in 2023. I gratefully acknowledge their hospitality. I especially thank my host Prof. Neela Nataraj so much for her wonderful hospitality and kindness during my stay. I also thank Dr Devika Shylaja for carefully reading an earlier version of this manuscript and providing many constructive comments.

References

[1] R. Adams, Sobolev Spaces, Academic Press New York, 1975.
[2] I. Babuška, J. Osborn, and J. Pitkäranta, Analysis of mixed methods using mesh dependent norms, Mathematics of Computation, 35 (1980), pp. 1039–1062.
[3] D. Braess, Finite Elements. Theory, fast solver, and applications in solid mechanics, Cambridge University Press, Second Edition, 2001.
[4] S. Brenner, Poincaré–Friedrichs inequalities for piecewise ${H}^{1}$ functions, SIAM Journal on Numerical Analysis, 41 (2003), pp. 306–324.
[5] S. Brenner and L. Scott, The Mathematical Theory of Finite Element Methods, Springer–Verlag, New York, 1994.
[6] S. Brenner and L.-Y. Sung, $C^{0}$ interior penalty methods for fourth order elliptic boundary value problems on polygonal domains, Journal of Scientific Computing, 22-23 (2005), pp. 83 – 118.
[7] P. Ciarlet, The finite element method for elliptic problems, North Holland, Amsterdam, 1978.
[8] P. Ciarlet and R. Glowinski, Dual iterative techniques for solving a finite element approximation of the biharmonic euation, Computer Methods in Applied Mechanics and Engineering, 5 (1975), pp. 277–295.
[9] P. Ciarlet and P.-A. Raviart, A mixed finite element method for the biharmonic equation, in Symposium on Mathematical Aspects of Finite Elements in Partial Differential Equations, C. D. Boor, ed., New York, 1974, Academic Press, pp. 125–143.
[10] C. Davini and I. Pitacco, An uncontrained mixed method for the biharmonic problem, SIAM Journal on Numerical Analysis, 38 (2001), pp. 820–836.
[11] G. Engel, K. Garikipati, T. Hughes, M. Larson, L. Mazzei, and R. Taylor, Continuous/discontinuous finite element approximations of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity, Computer Methods in Applied Mechanics and Engineering, 191 (2002), pp. 3669–3750.
[12] R. Falk, Approximation of the biharmonic equation by a mixed finite element method, SIAM Journal on Numerical Analysis, 15 (1978), pp. 556–567.
[13] R. Falk and J. Osborn, Error estimates for mixed methods, RAIRO Anal. Numér., 14 (1980), pp. 249–277.
[14] J. Freund and R. Stenberg, On weakly imposed boundary conditions for second order problems, in Proceedings of the Ninth International Conference on Finite Elements in Fluids, M. Morandi Cecchi, ed., Venice, 1995, pp. 327–326.
[15] V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer-Verlag, Berlin, 1986.
[16] P. Grisvard, Elliptic Problems in Nonsmooth Domains, vol. 24 of Monographs and Studies in Mathematics, Pitman (Advanced Publishing Program), Boston, MA, 1985.
[17] B. Lamichhane, Higher Order Mortar Finite Elements with Dual Lagrange Multiplier Spaces and Applications, PhD thesis, Universität Stuttgart, 2006.
[18] , A mixed finite element method for the biharmonic problem using biorthogonal or quasi-biorthogonal systems, Journal of Scientific Computing, 46 (2011), pp. 379–396.
[19] B. Lamichhane, R. Stevenson, and B. Wohlmuth, Higher order mortar finite element methods in 3D with dual Lagrange multiplier bases, Numerische Mathematik, 102 (2005), pp. 93–121.
[20] B. Lamichhane and B. Wohlmuth, Mortar finite elements for interface problems, Computing, 72 (2004), pp. 333–348.
[21] , Biorthogonal bases with local support and approximation properties, Mathematics of Computation, 76 (2007), pp. 233–249.
[22] J. Li, Full-order convergence of a mixed finite element method for fourth-order elliptic equations, Journal of Mathematical Analysis and Applications, 230 (1999), pp. 329–349.
[23] J.-L. Lions and E. Magenes, Non-homogeneous Boundary Value Problems and Applications. Vol. I, Springer-Verlag, New York, 1972. Translated from the French by P. Kenneth, Die Grundlehren der mathematischen Wissenschaften, Band 181.
[24] W. C. H. McLean, Strongly elliptic systems and boundary integral equations, Cambridge university press, 2000.
[25] S. E. Mikhailov, Traces, extensions and co-normal derivatives for elliptic systems on lipschitz domains, Journal of mathematical analysis and applications, 378 (2011), pp. 324–342.
[26] P. Monk, A mixed finite element method for the biharmonic equation, SIAM Journal on Numerical Analysis, 24 (1987), pp. 737–749.
[27] R. Scholz, A mixed method for 4th order problems using linear finite elements, RAIRO Anal. Numér., 12 (1978), pp. 85–90.
[28] R. Stenberg, On some techniques for approximating boundary conditions in the finite element method, Journal of Computational and Applied Mathematics, 63 (1995), pp. 139–148.
[29] V. Thomée, Galerkin finite element methods for parabolic problems, vol. 25, Springer Science & Business Media, 2007.
[30] G. Wells, E. Kuhl, and K. Garikipati, A discontinuous Galerkin method for the Cahn-Hilliard equation, Journal of Computational Physics, 218 (2006), pp. 860 – 877.
[31] W. Zulehner, The Ciarlet–Raviart method for biharmonic problems on general polygonal domains: Mapping properties and preconditioning, SIAM Journal on Numerical Analysis, 53 (2015), pp. 984–1004.

	$\displaystyle\alpha_{0}\\|(u_{h}-w_{h},\phi_{h}-\xi_{h})\\|_{a}$	$\displaystyle\leq$	$\displaystyle\sup_{(v_{h},\psi_{h})\in\mathcal{V}_{h}}\frac{a_{h}((u-w_{h},\phi-\xi_{h}),(v_{h},\psi_{h}))}{\\|(v_{h},\psi_{h})\\|_{a}}+Ch^{k}\\|p\\|_{k+1,\Omega}$
		$\displaystyle\leq$	$\displaystyle\\|(u-w_{h},\phi-\xi_{h})\\|_{a}+Ch^{k}\\|p\\|_{k+1,\Omega},$

	$\displaystyle\\|(u-u_{h},\phi-\phi_{h})\\|_{a}$	$\displaystyle\leq$	$\displaystyle\\|(u-w_{h},\phi-\xi_{h})\\|_{a}+\\|(w_{h}-u_{h},\xi_{h}-\phi_{h})\\|_{a}$
		$\displaystyle\leq$	$\displaystyle\left(1+\frac{1}{\alpha_{0}}\right)\\|(u-w_{h},\phi-\xi_{h})\\|_{a}+\frac{C}{\alpha_{0}}h^{k}\\|p\\|_{k+1,\Omega}.$