# A $C^{0}$ finite element method for the biharmonic problem with Dirichlet
boundary conditions in a polygonal domain

Hengguang Li†, Charuka D. Wickramasinghe†, Peimeng Yin∗ † Department of
Mathematics, Wayne State University, Detroit, MI 48202, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>∗ Multiscale Methods and Dynamics Group, Computer Science and
Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee
37831, USA<EMAIL_ADDRESS>
###### Abstract.

In this paper, we study the biharmonic equation with Dirichlet boundary
conditions in a polygonal domain. In particular, we propose a method that
effectively decouples the fourth-order problem into a system of two Poison
equations and one Stokes equation, or a system of one Stokes equation and one
Poisson equation. It is shown that the solution of each system is equivalent
to that of the original fourth-order problem on both convex and non-convex
polygonal domains. Two finite element algorithms are in turn proposed to solve
the decoupled systems. In addition, we show the regularity of the solutions in
each decoupled system in both the Sobolev space and the weighted Sobolev
space, and we derive the optimal error estimates for the numerical solutions
on both quasi-uniform meshes and graded meshes. Numerical test results are
presented to justify the theoretical findings.

###### Key words and phrases:

Biharmonic equation, reentrant corner, Stokes equation, Taylor-Hood method,
optimal error estimates.

###### 1991 Mathematics Subject Classification:

65N12, 65N30, 35J40

∗ Corresponding author.

##  1\. Introduction

We are interested in $C^{0}$ finite element method for the biharmonic problem

$\displaystyle\Delta^{2}\phi=f\quad{\rm{in}}\
\Omega,\qquad\quad\phi=0\quad{\rm{and}}\quad\partial_{\mathbf{n}}\phi=0\quad{\rm{on}}\
\partial\Omega,$ (1.1)

where $\Omega\subset\mathbb{R}^{2}$ is a polygonal domain, $\mathbf{n}$ is the
outward normal derivation, and $f$ is a given function. The boundary
conditions in (1.1) are referred to as the homogeneous Dirichlet boundary
conditions or clamped boundary conditions [22, 41] that occur for example in
fluid mechanics [25] and linear elasticity [17]. Equation (1.1) is a fourth-
order elliptic equation, there are three classical approaches to discretizing
the biharmonic equation in the literature. The first type is conforming finite
element methods, such as the Argyris finite element method [3], which requires
globally $C^{1}$ finite element spaces. The second type is the nonconforming
method, such as Morley finite element methods [37], which generally involve
low order polynomials. The third type is the mixed finite element methods that
require only Lagrange finite element spaces.

To obtain the mixed finite element methods, one first decomposes the fourth-
order equation into a lower order system, then applies the finite element
method to the system. Typically, problem (1.1) is decoupled into two Poisson
problems by introducing one intermediate function $\psi=-\Delta\phi$.
Different from biharmonic equation with Navior boundary conditions ($u=\Delta
u=0$ on $\partial\Omega$) that allows one to obtain two Poisson equations that
are completely decoupled [35], applying such decomposition to (1.1) leads to
two Poisson equation with either undertermined or overdetermined boundary
condition,

$\displaystyle\left\\{\begin{array}[]{ll}-\Delta\psi=f\quad{\rm{in}}\
\Omega,\\\ \text{no data}\quad{\rm{on}}\
\partial\Omega;\end{array}\right.\qquad\qquad{\rm{and}}\qquad\qquad\left\\{\begin{array}[]{ll}-\Delta\phi=\psi\quad{\rm{in}}\
\Omega,\\\ \phi=\partial_{\mathbf{n}}\phi=0\quad{\rm{on}}\
\partial\Omega.\end{array}\right.$

To overcome this difficulty, Ciarlet and Raviart [18] introduced a mixed
finite element method for problem (1.1) in a smooth domain by introducing a
conditioned function space. Later, Monk improved the mixed finite element
method in [36] to allow the intermediate function $\psi$ in $H^{1}(\Omega)$.
Their main difference lies in the smoothness of the intermediate function
$\psi$, namely $\psi$ is in different Sobolev spaces, but their solutions
could be equivalent under some smoothness assumptions [21]. However, these two
mixed finite element methods can only be applied to smooth domains and convex
polygonal domains. For non-convex polygonal domains, the solution of the
Ciarlet-Raviart mixed method is equivalent to the weak solution of (1.1), but
the corresponding mixed finite element solutions won’t converge to the exact
solution due to the low regularity [18, 21], while the Monk mixed method could
result in spurious solutions as discussed in [21], in which an augmented Monk
method was also introduced to remove the spurious solutions.

There is rich literature on mixed finite element methods for the biharmonic
problem (1.1). For example, the work in [20] studied a mixed method, and in
[45] preconditioning techniques were investigated for the discrete system,
both of which are based on the Ciarlet-Raviart formulation [18]. In addition,
more mixed variational formulations can be found in [44, 23].

We concern with the efficient numerical method for (1.1) in a polygonal
domain, especially, the non-convex polygonal domain. Motivated by the fact
that a Stokes problem could be reduced to a biharmonic problem on convex
polygonal domains or smooth domains [17, 15, 16, 5, 23], we resort to
decompose problem (1.1) in the following way. For given source terms $f$ of
the biharmonic problem (1.1), we first identify a source term of the Stokes
problem, though there are infinity many choices for the source term of the
Stokes problem, the velocity of the Stokes problem is uniquely determined by
$f$. Then we decompose the biharmonic problem (1.1) into a system of two
Poisson problems and one Stokes problem (P-S-P sequentially), or a system of
one Stokes problem and one Poisson problem (S-P sequentially), and we prove
that the solution of each decoupled system, namely the solution of the last
Poisson problem, is equivalent to the solution of the biharmonic problem (1.1)
in both convex and non-convex polygonal domains.

Many fast solvers or finite element methods are available for the involved
Stokes problem and the Poisson problems, see [17] for a review of the finite
element methods. For each decoupled system, we provide a finite element
algorithm, in which the Mini element method [4] or the Taylor-Hood element
method [25, 40, 11] is used to solve the Stokes problem, and the usual finite
element method [17, 12] is used for the Poisson problems. For each algorithm,
we carry out the error analysis for both the finite element approximations of
the biharmonic problem and the Taylor-Hood element approximations of the
involved Stokes problem.

For the algorithm involving only one Stokes problem and one Poisson problem
(S-P), the error estimates for the Stokes problem are standard. The error
estimate in $H^{1}$ norm for the biharmonic problem has an extra bound by the
$L^{2}$ error for the velocity of the Stokes problem, and it has a convergence
rate $h^{\min\\{k,\alpha_{0}+1,2\alpha_{0}\\}}$, where $k$ is the degree of
polynomials for finite element approximations of the biharmonic problem, and
$\alpha_{0}$ depending on the largest interior angle $\omega$ is given by
(2.2), which is visualized in Figure 2. The error estimate in $L^{2}$ norm has
an extra bound depending on $\omega$ by either the $L^{2}$ or the $H^{-1}$
error for the velocity of the Stokes problem, and it has a convergence rate
$h^{\min\\{k+1,\alpha_{0}+2,2\alpha_{0}\\}}$.

For the algorithm involving two Poisson problems and one Stokes problem
(P-S-P), the error estimates for the Poisson problem (the first Poisson
problem) sharing the same source term as the biharmonic equation are standard.
The $H^{1}$ error estimate for the velocity of the Stokes problem has an extra
bound by the $L^{2}$ error for the first Possion problem. The error estimates
in $L^{2}$ and $H^{-1}$ norms essentially use the property that the finite
element approximations for the source term of the Stokes problem and the
source term itself are invariant in inner product with any functions in curl
of finite element space (see Lemma 3.5 below). Finally, the $H^{1}$ error for
the biharmonic problem has a convergence rate
$h^{\min\\{k,\beta_{0}+2,\alpha_{0}+1,2\alpha_{0}\\}}$ with
$\beta_{0}=\frac{\pi}{\omega}$, and the error in $L^{2}$ norm has a
convergence rate $h^{\min\\{k+1,\beta_{0}+3,\alpha_{0}+2,2\alpha_{0}\\}}$. For
$k\leq 4$, we find $\beta_{0}+2$ in $H^{1}$ convergence rate and $\beta_{0}+3$
in $L^{2}$ convergence rate cannot achieve the minimum so the convergence
rates are the same as those in the previous algorithm. Moreover, two
algorithms give the same convergence rates in non-convex polygonal domains for
any $k\geq 1$.

In addition, we derive the regularity estimates for each proposed system in a
class of Kondratiev-type weighted Sobolev space. Based on the regularity
results, we in turn propose graded mesh refinement algorithms, such that the
associated finite element algorithms recover the optimal convergence rates for
both the finite element approximations of the biharmonic problem, and Mini
element approximation or Taylor-Hood approximations of the involved Stokes
problem. We shall point out that the thresholds of the grading parameter for
the graded meshes could depend on the norm of the error and the regularity of
the solution to a specific equation.

The rest of the paper is organized as follows. In Section 2, based on the
general regularity theory for second-order elliptic equations and the Stokes
equation, we review the weak solutions of all involved equations and show the
equivalence of the solution of the proposed system to that of the original
problem. In Section 3, we propose two finite element algorithms and obtain
error estimates on quasi-uniform meshes. In Section 4, we introduce a weighted
Sobolev space and derive regularity estimates for the solution. We also
present the graded mesh algorithm and provide optimal error estimates on
graded meshes. We report numerical test results in Section 5 to validate the
theory.

Throughout the paper, the generic constant $C>0$ in our estimates may be
different at different occurrences. It will depend on the computational
domain, but not on the functions involved nor on the mesh level in the finite
element algorithms.

##  2\. The biharmonic problem and its decoupled formulation

###  2.1. Well-posdedness and regularity of the solution

Denote by $H^{m}(\Omega)$, $m\geq 0$, the Sobolev space that consists of
functions whose $i$th ($0\leq i\leq m$) derivatives are square integrable.
Denote by $H^{1}_{0}(\Omega)\subset H^{1}(\Omega)$ the subspace consisting of
functions with zero trace on the boundary $\partial\Omega$. Let
$L^{2}(\Omega):=H^{0}(\Omega)$, we denote the norm $\|\cdot\|_{L^{2}(\Omega)}$
by $\|\cdot\|$ when there is no ambiguity about the underlying domain. For
$s>0$, let $s=m+t$, where $m\in\mathbb{Z}_{\geq 0}$ and $0<t<1$. Recall that
for $D\subseteq\mathbb{R}^{d}$, the fractional order Sobolev space $H^{s}(D)$
consists of distributions $v$ in $D$ satisfying

$\|v\|^{2}_{H^{s}(D)}:=\|v\|^{2}_{H^{m}(D)}+\sum_{|\alpha|=m}\int_{D}\int_{D}\frac{|\partial^{\alpha}v(x)-\partial^{\alpha}v(y)|^{2}}{|x-y|^{d+2t}}dxdy<\infty,$

where $\alpha=(\alpha_{1},\ldots,\alpha_{d})\in\mathbb{Z}^{d}_{\geq 0}$ is a
multi-index such that
$\partial^{\alpha}=\partial_{x_{1}}^{\alpha_{1}}\cdots\partial^{\alpha_{d}}_{x_{d}}$
and $|\alpha|=\sum_{i=1}^{d}\alpha_{i}$. We denote by $H_{0}^{s}(D)$ the
closure of $C_{0}^{\infty}(D)$ in $H^{s}(D)$, and $H^{-s}(D)$ the dual space
of $H_{0}^{s}(D)$. The notation $[\cdot]^{2}$ represents the vector space. For
example, $\mathbf{v}=(v_{1},v_{2})^{T}\in[H_{0}^{1}(\Omega)]^{2}$ represents
$v_{i}\in H_{0}^{1}(\Omega)$, $i=1,2$, where $T$ is the transposition of a
matrix or a vector. For $\mathbf{v}=(v_{1},v_{2})^{T}$, we denote ($\text{curl
}\mathbf{v}):=\frac{\partial v_{2}}{\partial x_{1}}-\frac{\partial
v_{1}}{\partial x_{2}}$. For a scalar function $\psi$, we denote
$(\textbf{curl }{\psi}):=({\psi}_{x_{2}},-{\psi}_{x_{1}})^{T}$.

By applying Green’s formulas, the variational formulation for biharmonic
problem (1.1) can be written as:

$\displaystyle a(\phi,\psi):=\int_{\Omega}\Delta\phi\Delta\psi
dx=\int_{\Omega}f\psi dx=(f,\psi)\quad\forall\psi\in H_{0}^{2}(\Omega).$ (2.1)

For a function $\psi\in H_{0}^{2}(\Omega)$, applying the Poincaré-type
inequality [27] twice, it follows

$a(\psi,\psi)=\|\Delta\psi\|^{2}=|\psi|^{2}_{H^{2}(\Omega)}\geq
C\|\psi\|^{2}_{H^{2}(\Omega)}.$

Thus, for any $f\in H^{-2}(\Omega)$, we have by the Lax-Milgram Theorem that
(2.1) admits a unique weak solution $\phi\in H_{0}^{2}(\Omega).$

$\theta=0$$\theta=\omega$$\Omega$$\bullet$$Q$ Figure 1. Domain $\Omega$
containing one reentrant corner.

The regularity of the solution $\phi$ depends on the given data $f$ and the
domain geometry [1, 9]. In order to decouple (1.1), we assume that the
polygonal domain $\Omega$ consists of $N$ vertices $Q_{i}$, $i=1,\cdots,N$,
and the corresponding interior angles are $\omega_{i}\in(0,2\pi)$. The largest
interior angle $\omega=\max_{i}\omega_{i}\in[\frac{\pi}{3},2\pi)$ associated
with the vertex $Q$. A sketch of the domain is given in Figure 1. We set
$z_{j}$, $j=1,2,\cdots,n$ the solutions of the following characteristic
equation corresponding to the the biharmonic problem (1.1) (see, e.g., [27]),

$\sin^{2}(z\omega)=z^{2}\sin^{2}(\omega),$

then there exists a threshold

$\alpha_{0}:=\min\\{\rm{Re}(z_{j}),\ j=1,2,\cdots,n\\}>\frac{1}{2},$ (2.2)

such that when $0\leq\alpha<\alpha_{0}$, the biharmonic problem (1.1) holds
the regularity estimate [29, 6, 10]

$\|\phi\|_{H^{2+\alpha}(\Omega)}\leq C\|f\|_{H^{-2+\alpha}(\Omega)}.$

The sketches of the threshold $\alpha_{0}$ in terms of the interior angle
$\omega$ is shown in Figure 2a. In Figure 2a, as a comparison we also show
$\beta_{0}=\frac{\pi}{\omega}$, which is the threshold of the characteristic
equation for the Poisson equation with homogeneous Dirichlet boundary
condition in the same polygonal domain. In Figure 2b, we show the difference
$\alpha_{0}-\beta_{0}$ in term of the interior angle $\omega$, from which we
find $\beta_{0}+1>\alpha_{0}$, when the largest interior angle $\omega$
satisfies $\frac{\pi}{3}<\omega^{\prime}<\omega<\pi$ for some
$\omega^{\prime}$. In Table 1, we present some numerical values of
$\alpha_{0}$ in terms of different interior angles $\omega$.

Figure 2. (a) $\alpha_{0},\beta_{0}$ in terms of $\omega$; (b) $\alpha_{0}-\beta_{0}$ in terms of $\omega$.  Table 1. Some numerical values of $\alpha_{0}$ in terms of different interior angles $\omega$. $\omega(<\pi)$ | $\alpha_{0}\approx$ | $\omega(>\pi)$ | $\alpha_{0}\approx$
---|---|---|---
$\frac{\pi}{3}$ | 4.059329012151345 | $\frac{7\pi}{6}$ | 0.751974545407645
$\frac{\pi}{2}$ | 2.739593356324596 | $\frac{6\pi}{5}$ | 0.717799308407060
$\frac{2\pi}{3}$ | 2.094139108847751 | $\frac{5\pi}{4}$ | 0.673583432221468
$\frac{3\pi}{4}$ | 1.885371778114173 | $\frac{4\pi}{3}$ | 0.615731059491289
$\frac{5\pi}{6}$ | 1.533859976323978 | $\frac{3\pi}{2}$ | 0.544483736993940
$\frac{11\pi}{12}$ | 1.200631594651580 | $\frac{7\pi}{4}$ | 0.505009699452470

For $\alpha_{0}$ and $\beta_{0}$, we have the following result from [29,
Theorem 7.1.1].

######  Lemma 2.1.

If $\omega\in(0,\pi)$, it follows that $\alpha_{0}$ in (2.2) satisfies

$\beta_{0}<\alpha_{0}<2\beta_{0},$ (2.3)

and if $\omega\in(\pi,2\pi)$, it follows

$\frac{1}{2}<\alpha_{0}<\beta_{0},$ (2.4)

where $\beta_{0}=\frac{\pi}{\omega}$.

By Lemma 2.1, when $\omega<\pi$, it follows

$\frac{1}{2}<\frac{\beta_{0}}{\alpha_{0}}<1,$ (2.5)

and when $\omega>\pi$, it follows

$\frac{\beta_{0}}{\alpha_{0}}>1.$ (2.6)

###  2.2. Decoupled formulation of the biharmonic problem

It is known that solving high order problems numerically, such as (1.1), is
much harder than solving the lower order problem. To decouple (1.1) into lower
order problems, we first introduce a steady-state Stokes problem

$\displaystyle-\Delta\mathbf{u}+\nabla p=$
$\displaystyle\mathbf{F}\quad\text{in }\Omega,$ (2.7) $\displaystyle\text{div
}\mathbf{u}=$ $\displaystyle 0\quad\text{in }\Omega,$
$\displaystyle\mathbf{u}=$ $\displaystyle 0\quad\text{on }\partial\Omega,$

where $\mathbf{u}=(u_{1},u_{2})^{T}$ is the velocity field of an
incompressible fluid motion, $p$ is the associated pressure. The source term
$\mathbf{F}=(f_{1},f_{2})^{T}$ satisfies

$\text{curl }\mathbf{F}=\frac{\partial f_{2}}{\partial x_{1}}-\frac{\partial
f_{1}}{\partial x_{2}}=f,$ (2.8)

for $f$ given in (1.1).

The weak formulation of the Stokes equations (2.7) is to find
$\mathbf{u}\in[H_{0}^{1}(\Omega)]^{2}$ and $p\in L_{0}^{2}(\Omega)$ such that

$\displaystyle(\nabla\mathbf{u},\nabla\mathbf{v})-(\text{div }\mathbf{v},p)=$
$\displaystyle\langle\mathbf{F},\mathbf{v}\rangle\quad\forall\mathbf{v}\in[H_{0}^{1}(\Omega)]^{2},$
(2.9) $\displaystyle-(\text{div }\mathbf{u},q)=$ $\displaystyle 0\quad\forall
q\in L_{0}^{2}(\Omega),$

where

$L_{0}^{2}(\Omega)=\\{q\in L^{2}(\Omega),\int_{\Omega}qdx=0\\}.$

For the bilinear forms in weak formulation (2.9), we have the following
Ladyzhenskaya-Babuska-Breezi (LBB) or inf-sup conditions,

$\displaystyle\inf_{q\in
L^{2}_{0}(\Omega)}\sup_{\mathbf{v}\in[H_{0}^{1}(\Omega)]^{2}}\frac{-(\text{div
}\mathbf{v},q)}{\|\mathbf{v}\|_{[H_{0}^{1}(\Omega)]^{2}}\|q\|}\geq\gamma_{1}>0,$
(2.10)
$\displaystyle\inf_{\mathbf{u}\in[H_{0}^{1}(\Omega)]^{2}}\sup_{\mathbf{v}\in[H_{0}^{1}(\Omega)]^{2}}\frac{(\nabla\mathbf{u},\nabla\mathbf{v})}{\|\mathbf{u}\|_{[H_{0}^{1}(\Omega)]^{2}}\|\mathbf{v}\|_{[H_{0}^{1}(\Omega)]^{2}}}\geq\gamma_{2}>0,$

and the boundedness

$\displaystyle(\nabla\mathbf{u},\nabla\mathbf{v})\leq
C_{1}\|\mathbf{u}\|_{[H_{0}^{1}(\Omega)]^{2}}\|\mathbf{v}\|_{[H_{0}^{1}(\Omega)]^{2}},$
(2.11) $\displaystyle(\text{div }\mathbf{v},q)\leq
C_{2}\|\mathbf{v}\|_{[H_{0}^{1}(\Omega)]^{2}}\|q\|,$

where $\gamma_{1},\gamma_{2},C_{1},C_{2}$ are constants.

Given that $\mathbf{F}\in[H^{-1}(\Omega)]^{2}$, under conditions (2.10) and
(2.11), the weak formulation (2.9) admits a unique solution
$(\mathbf{u},p)\in[H_{0}^{1}(\Omega)]^{2}\times L_{0}^{2}(\Omega)$ (see, e.g.
[30, 42, 24]). Moreover, if $\mathbf{F}\in[H^{-1+\alpha}(\Omega)]^{2}$ for
$\alpha<\alpha_{0}$, the Stokes problem holds the regularity estimate [8, 27,
38],

$\|\mathbf{u}\|_{[H^{1+\alpha}(\Omega)]^{2}}+\|p\|_{H^{\alpha}(\Omega)}\leq
C\|\mathbf{F}\|_{[H^{-1+\alpha}(\Omega)]^{2}}.$ (2.12)

Next, we introduce

$\mathbf{H}(\text{curl};\Omega):=\\{\mathbf{F}\in[L^{2}(\Omega)]^{2}:\text{curl
}\mathbf{F}\in L^{2}(\Omega)\\}.$

For given $f$, we introduce a Poisson problem

$-\Delta w=f\quad\text{in }\Omega,\quad w=0\quad{\rm{on}}\ \partial\Omega.$
(2.13)

Then we have the following results.

######  Lemma 2.2.

For $f\in H^{-1}(\Omega)$, assume that $w$ is the solution of (2.13), then it
follows

$\mathbf{F}=\textbf{curl
}w\in\mathbf{H}(\text{curl};\Omega)\subset[L^{2}(\Omega)]^{2}$ (2.14)

satisfies (2.8) and

$\|\mathbf{F}\|_{[L^{2}(\Omega)]^{2}}\leq C\|f\|_{H^{-1}(\Omega)}.$ (2.15)

###### Proof.

It is easy to verify that $\mathbf{F}$ satisfies (2.8), i.e.,

$\text{curl }\mathbf{F}=\text{curl }(\textbf{curl }w)=-\Delta w=f.$

The Poisson problem (2.13) admits a unique $w\in H_{0}^{1}(\Omega)$, which
satisfies

$\|w\|_{H^{1}(\Omega)}\leq C\|f\|_{H^{-1}(\Omega)}.$

Note that $\|\textbf{curl }w\|_{[L^{2}(\Omega)]^{2}}=|w|_{H^{1}(\Omega)}$, so
we have

$\|\textbf{curl }w\|_{[L^{2}(\Omega)]^{2}}\leq C\|f\|_{H^{-1}(\Omega)}.$
(2.16)

Thus, the estimate (2.15) holds. ∎

In a polygonal domain $\Omega$, if $f\in H^{l}(\Omega)$ for $l\geq-1$, the
regularity estimate [26, 27] for the Poisson problem (2.13) gives

$\|w\|_{H^{\min\\{1+\beta,l+2\\}}(\Omega)}\leq\|f\|_{H^{l}(\Omega)},$ (2.17)

where $\beta<\beta_{0}=\frac{\pi}{\omega}$ with $\omega$ being the largest
interior angles of $\Omega$. So for $\mathbf{F}$ obtained from (2.14), we have
$\mathbf{F}\in[H^{\min\\{\beta,l+1\\}}(\Omega)]^{2}$.

Since no boundary data is enforced to (2.8), so $\mathbf{F}$ obtained through
(2.8) is not unique. Assume that $\mathbf{F}_{0}\in[L^{2}(\Omega)]^{2}$ is a
solution of (2.8), then it follows that

$\mathbf{F}=\mathbf{F}_{0}+\nabla q\quad\forall q\in H^{1}(\Omega),$

is also a solution of (2.8) in $[L^{2}(\Omega)]^{2}$, since for $q\in
H^{1}(\Omega)$, we have $(\text{curl }\nabla q)\equiv 0.$

In addition to obtaining $\mathbf{F}$ by Lemma 2.2, we provide another way to
obtain $\mathbf{F}$.

######  Lemma 2.3.

Assume that $f\in L^{2}(\Omega)$.
(i) For any fixed $x_{2}$, if $f(\xi,x_{2})$ is integrable on $[c_{1},x_{1}]$
for some constant $c_{1}$, then

$\mathbf{F}=\left[0,\ \int_{c_{1}}^{x_{1}}f(\xi,x_{2})d\xi\right]^{T}$ (2.18)

satisfies (2.8).
(ii) Similarly, for any fixed $x_{1}$, if $f(x_{1},\zeta)$ is integrable on
$[c_{2},x_{2}]$ for some constant $c_{2}$, then

$\mathbf{F}=\left[-\int_{c_{2}}^{x_{2}}f(x_{1},\zeta)d\zeta,\ 0\right]^{T}$
(2.19)

also satisfies (2.8).
(iii) If both $f(\xi,x_{2})$ and $f(x_{1},\zeta)$ are integrable, then for any
constant $\eta$,

$\mathbf{F}=\left[-\eta\int_{c_{2}}^{x_{2}}f(x_{1},\zeta)d\zeta,\
(1-\eta)\int_{c_{1}}^{x_{1}}f(\xi,x_{2})d\xi\right]^{T}$ (2.20)

also satisfies (2.8).

It is obvious that $\mathbf{F}$ obtained from Lemma 2.3 satisfies
$\mathbf{F}\in[L^{2}(\Omega)]^{2}$. For all these
$\mathbf{F}\in[L^{2}(\Omega)]^{2}$ satisfying (2.8), we have the following
result.

######  Theorem 2.4.

Assume that $\mathbf{F}_{l}\in[L^{2}(\Omega)]^{2}$, $l=1,2$ both satisfy
(2.8). Let $(\mathbf{u}_{l},p_{l})$ be solutions of (2.7) or (2.9)
corresponding to $\mathbf{F}_{l}$, then it follows that

$\displaystyle\mathbf{u}_{1}=$ $\displaystyle\mathbf{u}_{2}\quad\text{in
}[H_{0}^{1}(\Omega)]^{2}\cap[H^{1+\alpha}(\Omega)]^{2},$ (2.21) $\displaystyle
p_{1}=$ $\displaystyle p_{2}+q\quad\text{in }L^{2}_{0}(\Omega)\cap
H^{\alpha}(\Omega),$

where $q\in L_{0}^{2}(\Omega)\cap H^{1}(\Omega)$ satisfies $\nabla
q=\mathbf{F}_{1}-\mathbf{F}_{2}$.

###### Proof.

We take
$\bar{\mathbf{F}}=\mathbf{F}_{1}-\mathbf{F}_{2}\in[L^{2}(\Omega)]^{2}$, then
by Helmholtz decomposition [24], there exist a stream-function $\psi$ and a
potential-function $q\in H^{1}(\Omega)$ uniquely up to a constant such that

$\bar{\mathbf{F}}=\nabla q+\textbf{curl }\psi,$ (2.22)

and

$(\bar{\mathbf{F}}-\nabla q)\cdot\mathbf{n}=(\textbf{curl
}\psi)\cdot\mathbf{n}=0\quad\text{in }H^{-\frac{1}{2}}(\partial\Omega).$
(2.23)

From (2.23), we have

$\frac{\partial\psi}{\partial\mathbf{\tau}}=(\textbf{curl
}\psi)\cdot\mathbf{n}=0\quad\text{in }H^{-\frac{1}{2}}(\partial\Omega),$

where $\tau$ is the unit tangential vector on $\partial\Omega$, thus we have

$\psi=C_{0}\quad\text{in }H^{\frac{1}{2}}(\partial\Omega),$ (2.24)

where $C_{0}$ is a constant. Take curl  on (2.22), we have

$-\Delta\psi=\text{curl }(\textbf{curl }\psi)=\text{curl }\bar{\mathbf{F}}=0,$
(2.25)

where the last equality is based on the fact that
$\mathbf{F}_{1},\mathbf{F}_{2}$ satisfy (2.8). By the Lax-Milgram Theorem, the
Poisson equation (2.25) with the boundary condition (2.24) admits a unique
solution $\psi=C_{0}$ in $H^{1}(\Omega)$. Therefore, the decomposition (2.22)
is equivalent to

$\bar{\mathbf{F}}=\nabla q.$ (2.26)

Let $\bar{\mathbf{u}}=\mathbf{u}_{1}-\mathbf{u}_{2}$ and
$\bar{p}=p_{1}-p_{2}$, then $(\bar{\mathbf{u}},\bar{p})$ satisfies

$\displaystyle-\Delta\bar{\mathbf{u}}+\nabla(\bar{p}-q)=$ $\displaystyle
0\quad\text{in }\Omega,$ (2.27) $\displaystyle\text{div }\bar{\mathbf{u}}=$
$\displaystyle 0\quad\text{in }\Omega,$ $\displaystyle\bar{\mathbf{u}}=$
$\displaystyle 0\quad\text{on }\partial\Omega.$

By the regularity of the Stokes problem (2.27), the conclusion holds. ∎

######  Lemma 2.5.

Assume that the source term $\mathbf{F}\in[L^{2}(\Omega)]^{2}$ of the Stokes
problem (2.7) is any vector function determined by $f\in H^{-1}(\Omega)$
satisfying (2.8), then (2.7) admits a unique solution
$\mathbf{u}\in[H^{1+\alpha}(\Omega)]^{2}$ and satisfies

$\|\mathbf{u}\|_{[H^{1+\alpha}(\Omega)]^{2}}\leq C\|f\|_{H^{-1}(\Omega)}.$
(2.28)

###### Proof.

Given $f\in H^{-1}(\Omega)$, we can always find a vector function
$\mathbf{F}_{0}\in[L^{2}(\Omega)]^{2}$ following Lemma 2.2 such that the
corresponding Stokes problem (2.7) admits a unique solution
$\mathbf{u}_{0}\in[H^{1+\alpha}(\Omega)]^{2}$ satisfying

$\|\mathbf{u}_{0}\|_{[H^{1+\alpha}(\Omega)]^{2}}\leq
C\|\mathbf{F}_{0}\|_{[H^{-1+\alpha}(\Omega)]^{2}}\leq
C\|\mathbf{F}_{0}\|_{[L^{2}(\Omega)]^{2}}\leq C\|f\|_{H^{-1}(\Omega)}.$ (2.29)

For any source term $\mathbf{F}$ also satisfying (2.14), it follows by Theorem
2.4 that the corresponding solution $\mathbf{u}=\mathbf{u}_{0}$, so the
conclusion holds. ∎

To show the connection of the solution $\mathbf{u}$ to the Stokes problem
(2.7) with the biharmonic problem (1.1), we introduce the following result
from [25, Theorem 3.1].

######  Lemma 2.6.

A function $\mathbf{v}\in[H^{m}(\Omega)]^{2}$ for integer $m\geq 0$ satisfies

$\text{div
}\mathbf{v}=0,\quad\langle\mathbf{v}\cdot\mathbf{n},1\rangle|_{\partial\Omega}=0,$

then there exists a stream function $\psi\in H^{m+1}(\Omega)$ uniquely up to
an additive constant satisfying

$\mathbf{v}=\textbf{curl }\psi.$

Since $\mathbf{u}\in[H_{0}^{1}(\Omega)]^{2}\cap[H^{1+\alpha}(\Omega)]^{2}$ and
$\text{div }\mathbf{u}=0$, so we have by Lemma 2.6 that there exists
$\bar{\phi}\in H^{2}(\Omega)$ uniquely up to an additive constant satisfying

$(u_{1},u_{2})^{T}=\mathbf{u}=\textbf{curl
}\bar{\phi}=(\bar{\phi}_{x_{2}},-\bar{\phi}_{x_{1}})^{T},$ (2.30)

which further implies $|\nabla\bar{\phi}|\in H^{1+\alpha}(\Omega)$, thus we
have

$\bar{\phi}\in H^{2+\alpha}(\Omega).$ (2.31)

######  Lemma 2.7.

There exists a unique

$\bar{\phi}\in H_{0}^{2}(\Omega)\cap H^{2+\alpha}(\Omega).$ (2.32)

satisfying (2.30).

###### Proof.

By calculation,

$\bar{\phi}_{\mathbf{\tau}}=\textbf{curl
}\bar{\phi}\cdot\mathbf{n}=\mathbf{u}\cdot\mathbf{n}=0,$

where $\mathbf{\tau}$ is the unit tangent to $\partial\Omega$, thus it follows

$\bar{\phi}=\text{constant}\quad\text{on }\partial\Omega.$

Without loss of generality, we can take

$\bar{\phi}=0\quad\text{on }\partial\Omega.$ (2.33)

From (2.30), we also have

$\nabla\bar{\phi}=(\bar{\phi}_{x_{1}},\bar{\phi}_{x_{2}})^{T}=(-u_{2},u_{1})^{T}=\mathbf{0}\quad\text{on
}\partial\Omega.$ (2.34)

Thus, the conclusion (2.32) follows from (2.33), (2.34) and (2.31). ∎

Instead of solving for $\bar{\phi}\in H_{0}^{2}(\Omega)\cap
H^{2+\alpha}(\Omega)$ from (2.30) directly, we apply the operator curl on
(2.30) to obtain the following Poisson problem

$-\Delta\bar{\phi}=\text{curl }\mathbf{u}\quad\text{in
}\Omega\quad\bar{\phi}=0\quad{\rm{on}}\ \partial\Omega.$ (2.35)

The weak formulation of (2.35) is to find $\bar{\phi}\in H_{0}^{1}(\Omega)$,
such that

$(\nabla\bar{\phi},\nabla\psi)=(\text{curl
}\mathbf{u},\psi)\quad\forall\psi\in H_{0}^{1}(\Omega).$ (2.36)

Since $(\text{curl }\mathbf{u})\in L^{2}(\Omega)$, so we have by the Lax-
Milgram Theorem that (2.36) admits a unique solution $\bar{\phi}\in
H_{0}^{1}(\Omega)$.

######  Lemma 2.8.

The Poisson problem (2.35) admits a unique solution $\bar{\phi}\in
H_{0}^{2}(\Omega)\cap H^{2+\alpha}(\Omega)$.

###### Proof.

Since $\bar{\phi}\in H_{0}^{2}(\Omega)\cap H^{2+\alpha}(\Omega)\subset
H_{0}^{1}(\Omega)$ is a solution of (2.30), so it is also a solution of the
Poisson problem (2.35). By the uniqueness of the solution of (2.35) in
$H_{0}^{1}(\Omega)$, the conclusion holds. ∎

######  Lemma 2.9.

The solution $\bar{\phi}\in H_{0}^{2}(\Omega)\cap H^{2+\alpha}(\Omega)$
obtained through (2.30) or the Poisson problem (2.35) satisfies the biharmonic
problem

$\Delta^{2}\bar{\phi}=\text{curl }\mathbf{F}=f\quad\text{in
}\Omega,\qquad\quad\bar{\phi}=0\quad{\rm{and}}\quad\partial_{\mathbf{n}}\bar{\phi}=0\quad{\rm{on}}\
\partial\Omega.$ (2.37)

###### Proof.

Following (2.30), we replace $\mathbf{u}$ by $\textbf{curl }\bar{\phi}$ in
(2.7) and obtain

$\displaystyle-\Delta(\bar{\phi}_{x_{2}})+p_{x_{1}}=$ $\displaystyle
f_{1}\quad\text{in }\Omega,$ (2.38a)
$\displaystyle-\Delta(-\bar{\phi}_{x_{1}})+p_{x_{2}}=$ $\displaystyle
f_{2}\quad\text{in }\Omega.$ (2.38b)

Applying differential operators $-\frac{\partial}{\partial x_{2}}$ and
$\frac{\partial}{\partial x_{1}}$ to (2.38a) and (2.38b), respectively, and
taking the summation lead to the conclusion. ∎

From Lemma 2.9, we find that $\bar{\phi}$ in (2.37) satisfies exactly the same
problem as $\phi$ in (1.1) in the following sense,

$\phi=\bar{\phi}\quad\text{in }H_{0}^{2}(\Omega)\cap H^{2+\alpha}(\Omega).$
(2.39)

From now on, we will use $\phi$ to replace the notation $\bar{\phi}$. Thus the
Poisson problem (2.35) is equivalent to

$-\Delta\phi=\text{curl }\mathbf{u}\quad\text{in
}\Omega,\quad\phi=0\quad{\rm{on}}\ \partial\Omega.$ (2.40)

The weak formulation (2.41) is equivalent to $\phi\in H_{0}^{1}(\Omega)$
satisfying

$(\nabla\phi,\nabla\psi)=(\text{curl }\mathbf{u},\psi)\quad\forall\psi\in
H_{0}^{1}(\Omega).$ (2.41)

By the regularity of the Poisson problem (2.40) and Lemma 2.5, we have that

$\|\phi\|_{H^{2+\alpha}(\Omega)}\leq C\|\text{curl
}\mathbf{u}\|_{H^{\alpha}(\Omega)}\leq
C\|\nabla\mathbf{u}\|_{[H^{\alpha}(\Omega)]^{2}}\leq
C\|\mathbf{u}\|_{[H^{1+\alpha}(\Omega)]^{2}}\leq C\|f\|_{H^{-1}(\Omega)}.$
(2.42)

In summary, we can obtain the solution $\phi$ of the biharmonic problem (1.1)
by solving the lower order problems in the following steps,

  1. (a)

Choose an appropriate $\mathbf{F}$ by Lemma 2.2 or Lemma 2.3;

  2. (b)

Solve $\mathbf{u}$ from the Stokes problem (2.7);

  3. (c)

Solve $\phi$ from the Poisson problem (2.40).

##  3\. The finite element method and error estimates

In this section, we propose a linear $C^{0}$ finite element method for solving
the biharmonic problem (1.1) based on the results in the previous section.

###  3.1. The finite element algorithm

Let $\mathcal{T}_{n}$ be a triangulation of $\Omega$ with shape-regular
triangles and let $\mathcal{P}_{k}(\mathcal{T}_{n})$ be the $C^{0}$ Lagrange
finite element space associated with $\mathcal{T}_{n}$,

$\mathcal{P}_{k}(\mathcal{T}_{n}):=\\{v\in C^{0}(\Omega):\ v|_{T}\in P_{k},\
\forall T\in\mathcal{T}_{n}\\},$ (3.1)

where $P_{k}$ is the space of polynomials of degree no more than $k$. Further,
we introduce the following specific $C^{0}$ Lagrange finite element spaces
associated with $\mathcal{T}_{n}$,

$\displaystyle V_{n}^{k}:=$ $\displaystyle\mathcal{P}_{k}(\mathcal{T}_{n})\cap
H_{0}^{1}(\Omega),$ (3.2) $\displaystyle S_{n}^{k}:=$
$\displaystyle\mathcal{P}_{k}(\mathcal{T}_{n})\cap L_{0}^{2}(\Omega),$

and the bubble function space

$B_{n}^{3}:=\\{v\in C^{0}(\Omega):\
v|_{T}\in\text{span}\\{\lambda_{1}\lambda_{2}\lambda_{3}\\},\ \forall
T\in\mathcal{T}_{n}\\},$

where $\lambda_{i}$, $i=1,2,3$ are the barycentric coordinates on $T$.

We define the finite element solution of the biharmonic problem (1.1) by
utilizing the decomposition in the previous section as follows.

######  Algorithm 3.1.

For any $f\in H^{-1}(\Omega)$ and $k\geq 1$, we consider the following steps.

  * •

Step 1. Find $w_{n}\in V_{n}^{k}$ of the Poisson equation

$(\nabla w_{n},\nabla\psi)=(f,\psi)\quad\forall\psi\in V_{n}^{k},$ (3.3)

then take $\mathbf{F}_{n}=\textbf{curl }w_{n}$.

  * •

Step 2. If $k=1$, we find the Mini element approximation $\mathbf{u}_{n}\times
p_{n}\in[V_{n}^{1}\oplus B_{n}^{3}]^{2}\times S_{n}^{1}$ of the Stokes
equation

$\displaystyle(\nabla\mathbf{u}_{n},\nabla\mathbf{v})-(p_{n},\text{div
}\mathbf{v})=$
$\displaystyle\langle\mathbf{F}_{n},\mathbf{v}\rangle\quad\forall\mathbf{v}\in[V_{n}^{1}\oplus
B_{n}^{3}]^{2},$ (3.4) $\displaystyle-(\text{div }\mathbf{u}_{n},q)=$
$\displaystyle 0\quad\forall q\in S_{n}^{1}.$

If $k\geq 2$, we find the Taylor-Hood element solution $\mathbf{u}_{n}\times
p_{n}\in[V_{n}^{k}]^{2}\times S_{n}^{k-1}$ of the Stokes equation

$\displaystyle(\nabla\mathbf{u}_{n},\nabla\mathbf{v})-(p_{n},\text{div
}\mathbf{v})=$
$\displaystyle\langle\mathbf{F}_{n},\mathbf{v}\rangle\quad\forall\mathbf{v}\in[V_{n}^{k}]^{2},$
(3.5) $\displaystyle-(\text{div }\mathbf{u}_{n},q)=$ $\displaystyle
0\quad\forall q\in S_{n}^{k-1}.$

  * •

Step 3. Find the finite element solution $\phi_{n}\in V_{n}^{k}$ of the
Poisson equation

$(\nabla\phi_{n},\nabla\psi)=(\text{curl
}\mathbf{u}_{n},\psi)\quad\forall\psi\in V_{n}^{k}.$ (3.6)

For some source term $f$, Algorithm 3.1 could be updated as follows.

######  Algorithm 3.2.

If $f\in L^{2}(\Omega)$ satisfies the condition of Lemma 2.3, we do the
following updates.

  * •

Step 1. This step is the same as Algorithm 3.1 Step 2 with
$\langle\mathbf{F}_{n},\mathbf{v}\rangle$ replaced by
$\langle\mathbf{F},\mathbf{v}\rangle$, where $\mathbf{F}$ obtained following
Lemma 2.3.

  * •

Step 2. The same as Algorithm 3.1 Step 3.

The finite element approximations for the Poisson problems in both Algorithm
3.1 and Algorithm 3.2 are well defined by the Lax-Milgram Theorem. We take the
Mini element method [4] and the Taylor-Hood element method [43, 14] for
solving the Stokes problem, other methods could also be used. The Mini element
approximations or the Taylor-Hood element approximations are well defined, if
i) the bilinear forms in Mini element method satisfy the following LBB
condition,

$\displaystyle\inf_{q\in S_{n}^{1}}\sup_{\mathbf{v}\in[V_{n}^{1}\oplus
B_{n}^{3}]^{2}}\frac{-(\text{div
}\mathbf{v},q)}{\|\mathbf{v}\|_{[H_{0}^{1}(\Omega)]^{2}}\|q\|}\geq\tilde{\gamma}_{1}>0,$
(3.7a) $\displaystyle\inf_{\mathbf{u}\in[V_{n}^{1}\oplus
B_{n}^{3}]^{2}}\sup_{\mathbf{v}\in[V_{n}^{1}\oplus
B_{n}^{3}]^{2}}\frac{(\nabla\mathbf{u},\nabla\mathbf{v})}{\|\mathbf{u}\|_{[H_{0}^{1}(\Omega)]^{2}}\|\mathbf{v}\|_{[H_{0}^{1}(\Omega)]^{2}}}\geq\tilde{\gamma}_{2}>0,$
(3.7b)

and these in Taylor-Hood method satisfies the following LBB condition,

$\displaystyle\inf_{q\in
S_{n}^{k-1}}\sup_{\mathbf{v}\in[V_{n}^{k}]^{2}}\frac{-(\text{div
}\mathbf{v},q)}{\|\mathbf{v}\|_{[H_{0}^{1}(\Omega)]^{2}}\|q\|}\geq\tilde{\gamma}_{1}>0,$
(3.8a)
$\displaystyle\inf_{\mathbf{u}\in[V_{n}^{k}]^{2}}\sup_{\mathbf{v}\in[V_{n}^{k}]^{2}}\frac{(\nabla\mathbf{u},\nabla\mathbf{v})}{\|\mathbf{u}\|_{[H_{0}^{1}(\Omega)]^{2}}\|\mathbf{v}\|_{[H_{0}^{1}(\Omega)]^{2}}}\geq\tilde{\gamma}_{2}>0,$
(3.8b)

where $\tilde{\gamma}_{1},\tilde{\gamma}_{2}$ are some constants; ii) the
bilinear forms are bounded.

######  Remark 3.3.

The Algorithm 3.1 is similar to method in [15, 16, 23] for fourth order
problems in smooth domains or convex polygonal domains. Error estimates were
derived in [16] for $P_{k}$ element approximations in a convex polygonal
domain by assuming that the solutions are smooth enough. In [23], error
estimates were given for $P_{1}$ element approximations based on the
regularity assumption $|\nabla\mathbf{u}|,|\mathbf{F}|,p,|\nabla\phi|\in
H^{s+1}(\Omega)$ for $-1<s\leq 0$. In this work, we carry out the error
analysis of Algorithm 3.1 and Algorithm 3.2 for biharmonic problem (1.1) in
both convex and non-convex polygonal domains for $P_{k}$ polynomials based on
the regularity estimates obtained in Section 2.

###  3.2. Optimal error estimates on quasi-uniform meshes

Suppose that the mesh $\mathcal{T}_{n}$ consists of quasi-uniform triangles
with size $h$. The interpolation error estimate on $\mathcal{T}_{n}$ (see
e.g., [17]) for any $v\in H^{\sigma}(\Omega)$, $\sigma>1$,

$\|v-v_{I}\|_{H^{\tau}(\Omega)}\leq
Ch^{\sigma-\tau}\|v\|_{H^{\sigma}(\Omega)},$ (3.9)

where $\tau=0,1$ and $v_{I}\in V_{n}^{k}$ represents the nodal interpolation
of $v$.

To make the analysis simple and clear, we we assume that $f\in
H^{\max\\{\alpha_{0},\beta_{0}\\}-1}(\Omega)\cap L^{2}(\Omega)$, where
$\alpha_{0}$ given in (2.2), and $\beta_{0}=\frac{\pi}{\omega}$. Since
Algorithm 3.2 involves only one Stokes equation and one Poisson equation, we
first give the error estimate of Algorithm 3.2.

For $f\in H^{\max\\{\alpha_{0},\beta_{0}\\}-1}(\Omega)\cap L^{2}(\Omega)$, if
$\mathbf{F}$ is given by Lemma 2.3, we have
$\mathbf{u}\in[H^{1+\alpha}(\Omega)]^{2},p\in H^{\alpha}(\Omega)$. Note that
the bilinear forms in the Mini element method ($k=1$) or Taylor-Hood element
method ($k\geq 2$) satisfying the LBB condition (3.7) or (3.8) on quasi-
uniform meshes [4, 43, 14], then the standard arguments for error estimate
(see e.g., [25, 40, 11]) give the following error estimate.

######  Lemma 3.4.

Let $(\mathbf{u},p)$ be the solution of the Stokes problem (2.9), and
$(\mathbf{u}_{n},p_{n})$ be the Mini element solution ($k=1$) or Taylor-Hood
element solution ($k\geq 2$) in Algorithm 3.2 on quasi-uniform meshes, then it
follows

$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{1}(\Omega)]^{2}}+\|p-p_{n}\|\leq
Ch^{\min\\{k,\alpha\\}},$ (3.10a)
$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\leq
Ch^{\min\\{k+1,\alpha+1,2\alpha\\}},$ (3.10b)
$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{-1}(\Omega)]^{2}}\leq
Ch^{\min\\{2k,k+2,k+\alpha,\alpha+2,2\alpha\\}}.$ (3.10c)

If the largest interior angle $\omega<\pi$, it follows
$\min\\{\alpha+1,2\alpha\\}=\alpha+1$, and if $\omega>\pi$, we have
$\min\\{\alpha+1,2\alpha\\}=2\alpha$.

If $\mathbf{F}$ is given by Lemma 2.2, then for the Poisson equations (2.13)
in a polygonal domain with $f\in
H^{\max\\{\alpha_{0},\beta_{0}\\}-1}(\Omega)\cap L^{2}(\Omega)$, the
regularity estimate gives $v\in H^{1+\beta}(\Omega)$ for
$\beta<\beta_{0}=\frac{\pi}{\omega}$ (see e.g., [26, 27]), it implies that

$\mathbf{F}=\textbf{curl
}w\in[H^{\beta}(\Omega)]^{2}\subset[H^{\alpha-1}(\Omega)]^{2}\cap[H^{\beta}(\Omega)]^{2},$

where $\alpha<\alpha_{0}$. Therefore, we also have
$\mathbf{u}\in[H^{\alpha+1}(\Omega)]^{2}\cap[H^{\beta+2}(\Omega)]^{2},p\in
H^{\alpha}(\Omega)\cap H^{\beta+1}(\Omega)$, and $\phi\in
H^{\alpha+2}(\Omega)\cap H^{\beta+3}(\Omega)$. For the finite element
approximations $w_{n}$ in (3.3), the standard error estimate [17] yields

$\|w-w_{n}\|_{H^{1}(\Omega)}\leq Ch^{\min\\{k,\beta\\}},\quad\|w-w_{n}\|\leq
Ch^{\min\\{k+1,\beta+1,2\beta\\}},$ (3.11)

which implies that

$\displaystyle\|\mathbf{F}-\mathbf{F}_{n}\|_{[L^{2}(\Omega)]^{2}}=$
$\displaystyle\|\textbf{curl }w-\textbf{curl
}w_{n}\|_{[L^{2}(\Omega)]^{2}}\leq\|w-w_{n}\|_{H^{1}(\Omega)}\leq
Ch^{\min\\{k,\beta\\}},$ (3.12)
$\displaystyle\|\mathbf{F}-\mathbf{F}_{n}\|_{[H^{-1}(\Omega)]^{2}}=$
$\displaystyle\|\textbf{curl }w-\textbf{curl
}w_{n}\|_{[H^{-1}(\Omega)]^{2}}\leq C\|w-w_{n}\|\leq
Ch^{\min\\{k+1,\beta+1,2\beta\\}}.$

For $\mathbf{F}_{n}$ in Algorithm 3.1, we further have the following result.

######  Lemma 3.5.

If $\mathbf{F}_{n}$ is given by $\mathbf{F}_{n}=\textbf{curl }w_{n}$ in Step 1
of Algorithm 3.1, and $\mathbf{F}=\textbf{curl }w$ is given in (2.14), then it
follows

$\langle\mathbf{F}-\mathbf{F}_{n},\textbf{curl
}\psi\rangle=0\quad\forall\psi\in V_{n}^{k}.$ (3.13)

###### Proof.

Subtract (3.3) from the weak formulation of (2.13), then we have the Galerkin
orthogonality,

$\displaystyle(\nabla(w-w_{n}),\nabla\psi)=(w-w_{n})_{x_{1}}\psi_{x_{1}}+(w-w_{n})_{x_{2}}\psi_{x_{2}}=0\quad\forall\psi\in
V_{n}^{k},$ (3.14)

which implies that

$\displaystyle\langle\mathbf{F}-\mathbf{F}_{n},\textbf{curl
}\psi\rangle=\langle\textbf{curl }(w-w_{n}),\textbf{curl
}\psi\rangle=(w-w_{n})_{x_{2}}\psi_{x_{2}}+(w-w_{n})_{x_{1}}\psi_{x_{1}}=0.$
(3.15)

∎

Next, we consider the error estimates of Taylor-Hood element approximations.
Subtract (3.5) from (2.9), we have the following equality,

$\displaystyle(\nabla(\mathbf{u}-\mathbf{u}_{n}),\nabla\mathbf{v})-(\text{div
}\mathbf{v},p-p_{n})=$
$\displaystyle\langle\mathbf{F}-\mathbf{F}_{n},\mathbf{v}\rangle\quad\forall\mathbf{v}\in[V_{n}^{k}]^{2},$
(3.16a) $\displaystyle-(\text{div }(\mathbf{u}-\mathbf{u}_{n}),q)=$
$\displaystyle 0\quad\forall q\in S_{n}^{k-1}.$ (3.16b)

We introduce the adjoint problem of the the Stokes equations (2.7),

$\displaystyle-\Delta\mathbf{r}+\nabla s=$
$\displaystyle\mathbf{g}\quad\text{in }\Omega,$ (3.17) $\displaystyle\text{div
}\mathbf{r}=$ $\displaystyle 0\quad\text{in }\Omega,$
$\displaystyle\mathbf{r}=$ $\displaystyle 0\quad\text{on }\partial\Omega,$

where $\mathbf{g}\in[H_{0}^{l}(\Omega)]^{2}$ for some $l=0,1$. Here, the
notation $H_{0}^{0}(\Omega):=H^{0}(\Omega)=L^{2}(\Omega)$. The weak
formulation of (3.17) is to find $\mathbf{r}\in[H_{0}^{1}(\Omega)]^{2}$ and
$s\in L_{0}^{2}(\Omega)$ such that

$\displaystyle(\nabla\mathbf{r},\nabla\mathbf{v})-(\text{div }\mathbf{v},s)=$
$\displaystyle\langle\mathbf{g},\mathbf{v}\rangle\quad\forall\mathbf{v}\in[H_{0}^{1}(\Omega)]^{2},$
(3.18a) $\displaystyle-(\text{div }\mathbf{r},q)=$ $\displaystyle
0\quad\forall q\in L_{0}^{2}(\Omega).$ (3.18b)

We have the following regularity result,

$\|\mathbf{r}\|_{[H^{1+\min\\{\alpha,l+1\\}}(\Omega)]^{2}}+\|s\|_{H^{\min\\{\alpha,l+1\\}}(\Omega)}\leq
C\|\mathbf{g}\|_{[H^{\min\\{\alpha,l+1\\}-1}(\Omega)]^{2}}\leq
C\|\mathbf{g}\|_{[H^{l}(\Omega)]^{2}},$ (3.19)

where $\alpha<\alpha_{0}$.

Note that $\mathbf{r}\in[H^{1+\min\\{\alpha,l+1\\}}(\Omega)]^{2}$ satisfying
(3.19) and (3.17), we have by Lemma 2.6 that there exists $\psi\in
H^{2+\min\\{\alpha,l+1\\}}(\Omega)\cap H_{0}^{1}(\Omega)$ such that

$\mathbf{r}=\textbf{curl }\psi.$ (3.20)

We also have that $\|\psi\|_{H^{2+\min\\{\alpha,l+1\\}}(\Omega)}\leq
C\|\mathbf{r}\|_{[H^{1+\min\\{\alpha,l+1\\}}(\Omega)]^{2}}$.

Let $(\mathbf{r}_{n},s_{n})$ be the Taylor-Hood solution of (3.17), we have

$\displaystyle(\nabla\mathbf{r}_{n},\nabla\mathbf{v})-(\text{div
}\mathbf{v},s_{n})=$
$\displaystyle\langle\mathbf{g},\mathbf{v}\rangle\quad\forall\mathbf{v}\in[V_{n}^{k}]^{2},$
(3.21a) $\displaystyle-(\text{div }\mathbf{r}_{n},q)=$ $\displaystyle
0\quad\forall q\in S_{n}^{k-1}.$ (3.21b)

Subtracting (3.21) from (3.18), we have the Galerkin orthogonality,

$\displaystyle(\nabla(\mathbf{r}-\mathbf{r}_{n}),\nabla\mathbf{v})-(\text{div
}\mathbf{v},s-s_{n})=$ $\displaystyle
0\quad\forall\mathbf{v}\in[V_{n}^{k}]^{2},$ (3.22a) $\displaystyle-(\text{div
}(\mathbf{r}-\mathbf{r}_{n}),q)=$ $\displaystyle 0\quad\forall q\in
S_{n}^{k-1}.$ (3.22b)

######  Lemma 3.6.

Let $(\mathbf{u},p)$ be the solution of the Stokes problem (2.9), and
$(\mathbf{u}_{n},p_{n})$ be the Mini element solution ($k=1$) or Taylor-Hood
element solution ($k\geq 2$) in Algorithm 3.1 on quasi-uniform meshes, then it
follows the error estimates

$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{1}(\Omega)]^{2}}+\|p-p_{n}\|\leq
Ch^{\min\\{k,\alpha,\beta+1\\}},$ (3.23a)
$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\leq
Ch^{\min\\{k+1,\alpha+1,\beta+2,2\alpha\\}},$ (3.23b)
$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{-1}(\Omega)]^{2}}\leq
Ch^{\min\\{2k,k+2,k+\beta,\alpha+2,\beta+3,2\alpha\\}}.$ (3.23c)

###### Proof.

We will only present the proof of error estimates of the Taylor-Hood element
approximations, the error estimates of the Mini element approximations can be
proved similarly. By Lemma A.1, we have

$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{1}(\Omega)]^{2}}+\|p-p_{n}\|\leq$
$\displaystyle
C\left(\inf_{\mathbf{v}\in[V_{n}^{k}]^{2}}\|\mathbf{u}-\mathbf{v}\|_{[H^{1}(\Omega)]^{2}}+\inf_{q\in
S_{n}^{k-1}}\|p-q\|+\|\mathbf{F}-\mathbf{F}_{n}\|_{[H^{-1}(\Omega)]^{2}}\right)$
$\displaystyle\leq$ $\displaystyle
Ch^{\min\\{k,\alpha,\beta+1\\}}+Ch^{\min\\{k+1,\beta+1,2\beta\\}}\leq
Ch^{\min\\{k,\alpha,\beta+1\\}}.$

We take $\mathbf{v}=\mathbf{u}-\mathbf{u}_{n}$, $q=p-p_{n}$ in (3.18), then we
have

$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{-l}(\Omega)]^{2}}=\sup_{g\in[H_{0}^{l}(\Omega)]^{2}}\frac{\langle\mathbf{g},\mathbf{u}-\mathbf{u}_{n}\rangle}{\|\mathbf{g}\|_{[H^{l}(\Omega)]^{2}}},$

where $l=0,1$. Subtracting (3.21) with $\mathbf{v}=\mathbf{u}_{n}$ from (3.18)
with $\mathbf{v}=\mathbf{u}$, we have

$\displaystyle\langle\mathbf{g},\mathbf{u}-\mathbf{u}_{n}\rangle=$
$\displaystyle(\nabla\mathbf{r},\nabla\mathbf{u})-(\nabla\mathbf{r}_{n},\nabla\mathbf{u}_{n})-(\text{div
}\mathbf{u},s)+(\text{div }\mathbf{u}_{n},s_{n})$ $\displaystyle=$
$\displaystyle(\nabla(\mathbf{u}-\mathbf{u}_{n}),\nabla\mathbf{r})-(\text{div
}(\mathbf{u}-\mathbf{u}_{n}),s)+(\nabla(\mathbf{r}-\mathbf{r}_{n}),\nabla\mathbf{u}_{n})-(\text{div
}\mathbf{u}_{n},s-s_{n})$ $\displaystyle=$
$\displaystyle(\nabla(\mathbf{u}-\mathbf{u}_{n}),\nabla\mathbf{r})-(\text{div
}(\mathbf{u}-\mathbf{u}_{n}),s)$

where we have used (3.22a) with $\mathbf{v}=\mathbf{u}_{n}$ in the last
equality. Subtracting (3.16a) with $\mathbf{v}=\mathbf{r}_{n}$, (3.16b) with
$q=s_{n}$, and (3.18b) with $q=-(p-p_{n})$ from the equation above,
respectively, we have

$\displaystyle\langle\mathbf{g},\mathbf{u}-\mathbf{u}_{n}\rangle=$
$\displaystyle(\nabla(\mathbf{u}-\mathbf{u}_{n}),\nabla(\mathbf{r}-\mathbf{r}_{n}))-(\text{div
}(\mathbf{r}-\mathbf{r}_{n}),p-p_{n})$ (3.24) $\displaystyle-(\text{div
}(\mathbf{u}-\mathbf{u}_{n}),s-s_{n})+\langle\mathbf{F}-\mathbf{F}_{n},\mathbf{r}_{n}\rangle$
$\displaystyle:=$ $\displaystyle T_{1}+T_{2}+T_{3}+T_{4},$

For $\mathbf{r}-\mathbf{r}_{n}$ and $s-s_{n}$, we have the following estimate

$\displaystyle\|\mathbf{r}-\mathbf{r}_{n}\|_{[H^{1}(\Omega)]^{2}}+\|s-s_{n}\|\leq
C\left(\inf_{\mathbf{r}_{I}\in[V_{n}^{k}(\Omega)]^{2}}\|\mathbf{r}-\mathbf{r}_{I}\|_{[H^{1}(\Omega)]^{2}}+\inf_{s_{I}\in
S_{n}^{k-1}}\|s-s_{I}\|\right)$ (3.25) $\displaystyle\leq
Ch^{\min\\{k,\alpha,l+1\\}}(\|\mathbf{r}\|_{[H^{\min\\{k+1,\alpha+1,l+2\\}}(\Omega)]^{2}}+\|s\|_{H^{\min\\{k,\alpha,l+1\\}}(\Omega)}).$

Then we have the following estimates for terms $T_{i}$, $i=1,\cdots,4$ in
(3.24).

$\displaystyle|T_{1}|\leq$
$\displaystyle|\mathbf{u}-\mathbf{u}_{n}|_{[H^{1}(\Omega)]^{2}}|\mathbf{r}-\mathbf{r}_{n}|_{[H^{1}(\Omega)]^{2}}\leq
Ch^{\min\\{k,\alpha,\beta+1\\}+\min\\{k,\alpha,l+1\\}}\|\mathbf{r}\|_{[H^{1+\min\\{k,\alpha,l+1\\}}(\Omega)]^{2}}$
$\displaystyle=$ $\displaystyle
Ch^{\min\\{2k,k+l+1,k+\alpha,\alpha+l+1,k+\beta+1,\alpha+\beta+1,\beta+l+2,2\alpha\\}}\|\mathbf{r}\|_{[H^{1+\min\\{k,\alpha,l+1\\}}(\Omega)]^{2}}.$
$\displaystyle|T_{2}|\leq\|p-p_{n}\||\mathbf{r}-\mathbf{r}_{n}|_{[H^{1}(\Omega)]^{2}}\leq
Ch^{\min\\{2k,k+l+1,k+\alpha,\alpha+l+1,k+\beta+1,\alpha+\beta+1,\beta+l+2,2\alpha\\}}\|\mathbf{r}\|_{[H^{1+\min\\{k,\alpha,l+1\\}}(\Omega)]^{2}}.$
$\displaystyle|T_{3}|\leq|\mathbf{u}-\mathbf{u}_{n}|_{[H^{1}(\Omega)]^{2}}|\|s-s_{n}\|\leq
Ch^{\min\\{2k,k+l+1,k+\alpha,\alpha+l+1,,k+\beta+1,\alpha+\beta+1,\beta+l+2,2\alpha\\}}\|s\|_{H^{\min\\{k,\alpha,l+1\\}}(\Omega)}.$

For $T_{4}$, we have

$\displaystyle
T_{4}=\langle\mathbf{F}-\mathbf{F}_{n},\mathbf{r}_{n}-\mathbf{r}\rangle+\langle\mathbf{F}-\mathbf{F}_{n},\mathbf{r}\rangle:=T_{41}+T_{42}.$
$\displaystyle|T_{41}|\leq$
$\displaystyle\|\mathbf{F}-\mathbf{F}_{n}\|_{[H^{-1}(\Omega)]^{2}}\|\mathbf{r}-\mathbf{r}_{n}\|_{[H^{1}(\Omega)]^{2}}\leq
Ch^{\min\\{k+1,\beta+1,2\beta\\}+\min\\{k,\alpha,l+1\\}}\|\mathbf{r}\|_{[H^{1+\min\\{\alpha,l+1\\}}(\Omega)]^{2}}$
$\displaystyle=$ $\displaystyle
Ch^{\min\\{2k+1,k+l+2,k+\beta+1,\beta+l+2,\alpha+\beta+1,2\beta+\alpha\\}}\|\mathbf{r}\|_{[H^{1+\min\\{\alpha,l+1\\}}(\Omega)]^{2}}.$

By (3.20) and Lemma 3.5, we have

$\displaystyle|T_{42}|=$
$\displaystyle|\langle\mathbf{F}-\mathbf{F}_{n},\textbf{curl
}\psi\rangle|=|\langle\mathbf{F}-\mathbf{F}_{n},\textbf{curl
}(\psi-\psi_{I})\rangle|\leq\|\mathbf{F}-\mathbf{F}_{n}\|_{[L^{2}(\Omega)]^{2}}\|\textbf{curl
}(\psi-\psi_{I})\|_{[L^{2}(\Omega)]^{2}}$ $\displaystyle\leq$
$\displaystyle\|\mathbf{F}-\mathbf{F}_{n}\|_{[L^{2}(\Omega)]^{2}}\|\psi-\psi_{I}\|_{H^{1}(\Omega)}\leq
Ch^{\min\\{k,\beta\\}+\min\\{k,\alpha+1,l+2\\}}\|\mathbf{r}\|_{[H^{1+\min\\{\alpha,l+1\\}}(\Omega)]^{2}}$
$\displaystyle\leq$ $\displaystyle
Ch^{\min\\{2k,k+l+2,k+\beta,\beta+l+2,\alpha+\beta+1\\}}\|\mathbf{r}\|_{[H^{1+\min\\{\alpha,l+1\\}}(\Omega)]^{2}}.$

where $\psi_{I}$ is the nodal interpolation of $\psi$. It can be verified that

$|T_{4}|\leq|T_{41}|+|T_{42}|\leq
Ch^{\min\\{2k,k+l+2,k+\beta,\beta+l+2,\alpha+\beta+1,\alpha+2\beta\\}}\|\mathbf{r}\|_{[H^{1+\min\\{\alpha,l+1\\}}(\Omega)]^{2}}.$

By the regularity (3.19) and the summation of estimates $|T_{i}|$,
$i=1,\cdots,4$, the error estimate
$\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{-l}(\Omega)]^{2}}$ holds.

∎

######  Theorem 3.7.

Let $\phi_{n}\in V_{n}^{k}$ be the solution of finite element solution of
(3.6) from Algorithm 3.2 or Algorithm 3.1, and $\phi$ is the solution of the
biharmonic problem (1.1), if $\phi_{n}$ is the solution of Algorithm 3.2, then
we have

$\|\phi-\phi_{n}\|_{H^{1}(\Omega)}\leq Ch^{\min\\{k,\alpha+1,2\alpha\\}};$
(3.26)

if it is the solution of of Algorithm 3.1, then we have

$\|\phi-\phi_{n}\|_{H^{1}(\Omega)}\leq
Ch^{\min\\{k,\beta+2,\alpha+1,2\alpha\\}}.$ (3.27)

###### Proof.

Subtracting (3.6) from (2.41) gives

$(\nabla(\phi-\phi_{n}),\nabla\psi)=(\text{curl
}(\mathbf{u}-\mathbf{u}_{n}),\psi)\quad\forall\psi\in V_{n}^{k}.$

Denote by $\phi_{I}\in V_{n}^{k}$ the nodal interpolation of $\phi$. Set
$\epsilon=\phi_{I}-\phi$, $e=\phi_{I}-\phi_{n}$ and take $\psi=e$, then we
have

$(\nabla e,\nabla e)=(\nabla\epsilon,\nabla e)+(\text{curl
}(\mathbf{u}-\mathbf{u}_{n}),e)=(\nabla\epsilon,\nabla
e)+(\mathbf{u}-\mathbf{u}_{n},\textbf{curl }e),$

which gives

$\displaystyle\|e\|^{2}_{H^{1}(\Omega)}\leq$
$\displaystyle\|\epsilon\|_{H^{1}(\Omega)}\|e\|_{H^{1}(\Omega)}+\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\|\textbf{curl
}e\|_{[L^{2}(\Omega)]^{2}}$ (3.28) $\displaystyle\leq$ $\displaystyle
C(\|\epsilon\|_{H^{1}(\Omega)}+\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}})\|e\|_{H^{1}(\Omega)},$

By the triangle inequality, we have

$\displaystyle\|\phi-\phi_{n}\|_{H^{1}(\Omega)}\leq$
$\displaystyle\|\epsilon\|_{H^{1}(\Omega)}+\|e\|_{H^{1}(\Omega)}\leq
C\left(\|\epsilon\|_{H^{1}(\Omega)}+\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\right)$
(3.29)

Recall that $\phi\in H^{2+\alpha}(\Omega)$, so it follows

$\|\epsilon\|_{H^{1}(\Omega)}\leq Ch^{\min\\{k,1+\alpha\\}},$

which together with (3.10) for Algorithm 3.2 and (3.23) for Algorithm 3.1
leads to the conclusion. ∎

######  Theorem 3.8.

Let $\phi_{n}$ be the solution of finite element solution of (3.6) from
Algorithm 3.2 or Algorithm 3.1, and $\phi$ be the solution of the biharmonic
problem (1.1). If $\phi_{n}$ is the solution of Algorithm 3.2, then we have

$\|\phi-\phi_{n}\|\leq Ch^{\min\\{k+1,\alpha+2,2\alpha\\}};$ (3.30)

if it is the solution of Algorithm 3.1, then we have

$\|\phi-\phi_{n}\|\leq Ch^{\min\\{k+1,\beta+3,\alpha+2,2\alpha\\}}.$ (3.31)

###### Proof.

Consider the Poisson problem

$-\Delta v=\phi-\phi_{n}\text{ in }\Omega\quad v=0\text{ on }\partial\Omega.$
(3.32)

Then we have

$\|\phi-\phi_{n}\|^{2}=(\nabla(\phi-\phi_{n}),\nabla v).$ (3.33)

By Subtracting (3.6) from (2.41), it follows

$(\nabla(\phi-\phi_{n}),\nabla\psi)=(\text{curl
}(\mathbf{u}-\mathbf{u}_{n}),\psi)\quad\forall\psi\in V_{n}^{k}.$ (3.34)

Set $\psi=v_{I}\in V_{n}^{k}$ the nodal interpolation of $v$ and subtract
(3.34) from (3.32), we have

$\displaystyle\|\phi-\phi_{n}\|^{2}=$
$\displaystyle(\nabla(\phi-\phi_{n}),\nabla(v-v_{I}))+(\text{curl
}(\mathbf{u}-\mathbf{u}_{n}),v_{I}),$ $\displaystyle=$
$\displaystyle(\nabla(\phi-\phi_{n}),\nabla(v-v_{I}))+(\text{curl
}(\mathbf{u}-\mathbf{u}_{n}),v_{I}-v)+(\text{curl
}(\mathbf{u}-\mathbf{u}_{n}),v),$ $\displaystyle=$
$\displaystyle(\nabla(\phi-\phi_{n}),\nabla(v-v_{I}))+(\mathbf{u}-\mathbf{u}_{n},\textbf{curl
}(v_{I}-v))+(\mathbf{u}-\mathbf{u}_{n},\textbf{curl }v),$ $\displaystyle\leq$
$\displaystyle\|\phi-\phi_{n}\|_{H^{1}(\Omega)}\|v-v_{I}\|_{H^{1}(\Omega)}+\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\|v-v_{I}\|_{H^{1}(\Omega)}$
$\displaystyle+\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{-\min\\{\lfloor\beta\rfloor,1\\}}(\Omega)]^{2}}\|\textbf{curl
}v\|_{H^{\min\\{\lfloor\beta\rfloor,1\\}}(\Omega)},$

where $\lfloor\cdot\rfloor$ represents the floor function.

The regularity result [26, 27] of the Poisson problem (3.33) gives

$\|v\|_{H^{\min\\{1+\beta,2\\}}(\Omega)}\leq
C\|\phi-\phi_{n}\|_{H^{\min\\{\beta-1,0\\}}(\Omega)}\leq C\|\phi-\phi_{n}\|.,$
(3.35)

where $\beta<\frac{\pi}{\omega}$. From (3.9), we have

$\|v-v_{I}\|_{H^{1}(\Omega)}\leq
Ch^{\min\\{\beta,1\\}}\|v\|_{H^{\min\\{1+\beta,2\\}}(\Omega)}.$

For Algorithm 3.2, we have the following result by (3.10). Since
$\beta<\frac{\pi}{\omega}$, so if $\omega>\pi$ we have
$\lfloor\beta\rfloor=0$, and

$\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{-\min\\{\lfloor\beta\rfloor,1\\}}(\Omega)]^{2}}=\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\leq
Ch^{2\alpha}.$ (3.36)

and if $\omega<\pi$, we have $\lfloor\beta\rfloor=1$, and

$\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{-\min\\{\lfloor\beta\rfloor,1\\}}(\Omega)]^{2}}\leq
Ch^{\min\\{2k,k+2,k+\alpha,\alpha+2,2\alpha\\}}.$ (3.37)

For $\omega\in(0,2\pi)\setminus\\{\pi\\}$, (3.36) and (3.37) imply that

$\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{-\min\\{\lfloor\beta\rfloor,1\\}}(\Omega)]^{2}}\leq
Ch^{\min\\{2k,k+2,k+\alpha,\alpha+2,2\alpha\\}}.$ (3.38)

Thus, we have by (3.26), (3.10), and (3.38),

$\displaystyle\|\phi-\phi_{n}\|^{2}\leq$ $\displaystyle
Ch^{\min\\{k+1,\alpha+2,k+\beta,2\alpha+\beta\\}}\|v\|_{H^{\min\\{1+\beta,2\\}}(\Omega)}+Ch^{\min\\{2k,k+2,k+\alpha,\alpha+2,2\alpha\\}}\|v\|_{H^{\min\\{1+\beta,2\\}}(\Omega)}$
(3.39) $\displaystyle\leq$ $\displaystyle
Ch^{\min\\{k+1,\alpha+2,2\alpha\\}}\|v\|_{H^{\min\\{1+\beta,2\\}}(\Omega)}.$

By (3.35) and (3.39), the estimate (3.30) holds.

Similarly, for Algorithm 3.1 we have

$\displaystyle\|\phi-\phi_{n}\|^{2}\leq$ $\displaystyle
Ch^{\min\\{k+1,\alpha+2,\beta+3,k+\beta,2\alpha+\beta\\}}\|v\|_{H^{\min\\{1+\beta,2\\}}(\Omega)}$
(3.40)
$\displaystyle+Ch^{\min\\{2k,k+2,k+\beta,\alpha+2,\beta+3,2\alpha\\}}\|v\|_{H^{\min\\{1+\beta,2\\}}(\Omega)}$
$\displaystyle\leq$ $\displaystyle
Ch^{\min\\{k+1,\alpha+2,\beta+3,2\alpha\\}}\|v\|_{H^{\min\\{1+\beta,2\\}}(\Omega)}.$

By (3.35) and (3.40), the estimate (3.31) holds. ∎

######  Remark 3.9.

For results in Theorem 3.7 and Theorem 3.8, we have the following comparisons.

  * •

From Figure 2(b), we find $\alpha>\beta+1$ only when the largest interior
angle $\omega$ of the domain $\Omega$ is close to $\frac{\pi}{3}$, so there
exists $\omega_{0}\in(\frac{\pi}{3},\frac{\pi}{2})$ such that when
$\omega>\omega_{0}$ the Algorithm 3.1 and the Algorithm 3.2 give the same
convergence rates. In particular, these two algorithms give the same
convergence rates on non-convex domains.

  * •

For $k\leq 4$, from Figure 2 we can find $\beta+2$ in estimate (3.27) and
$\beta+3$ in estimate (3.31) cannot achieve the minimum, so these two
algorithms also give the same convergence rates.

##  4\. Optimal error estimates on graded meshes

To improve the convergence rate, we consider the Algorithm 3.2 on graded
meshes. We start with the regularity in weighted Sobolev space.

###  4.1. Weighted Sobolev space

Recall that $Q_{i}$, $i=1,\cdots,N$ are the vertices of domain $\Omega$. Let
$r_{i}=r_{i}(x,Q_{i})$ be the distance from $x$ to $Q_{i}$ and let

$\displaystyle\rho(x)=\Pi_{1\leq i\leq N}r_{i}(x,Q_{i}).$ (4.1)

Let $\mathbf{a}=(a_{1},\cdots,a_{i},\cdots,a_{N})$ be a vector with $i$th
component associated with $Q_{i}$. We denote
$t+\mathbf{a}=(t+a_{1},\cdots,t+a_{N})$, so we have

$\rho(x)^{(t+\mathbf{a})}=\Pi_{1\leq i\leq
N}r_{i}^{(t+\mathbf{a})}(x,Q_{i})=\Pi_{1\leq i\leq
N}r_{i}^{t}(x,Q_{i})\Pi_{1\leq i\leq
N}r_{i}^{\mathbf{a}}(x,Q_{i})=\rho(x)^{t}\rho(x)^{\mathbf{a}}.$

Then, we introduce the Kondratiev-type weighted Sobolev spaces for the
analysis of the Stokes problem (2.9) and the Poisson problem (2.40).

######  Definition 4.1.

(Weighted Sobolev spaces) For $a\in\mathbb{R}$, $m\geq 0$, and
$G\subset\Omega$, we define the weighted Sobolev space

$\mathcal{K}_{\mathbf{a}}^{m}(G):=\\{v|\
\rho^{|\nu|-\mathbf{a}}\partial^{\nu}v\in L^{2}(G),\forall\ |\nu|\leq m\\},$

where the multi-index $\nu=(\nu_{1},\nu_{2})\in\mathbb{Z}^{2}_{\geq 0}$,
$|\nu|=\nu_{1}+\nu_{2}$, and
$\partial^{\nu}=\partial_{x}^{\nu_{1}}\partial_{y}^{\nu_{2}}$. The
$\mathcal{K}_{\mathbf{a}}^{m}(G)$ norm for $v$ is defined by

$\|v\|_{\mathcal{K}_{\mathbf{a}}^{m}(G)}=\big{(}\sum_{|\nu|\leq
m}\iint_{G}|\rho^{|\nu|-\mathbf{a}}\partial^{\alpha}v|^{2}dxdy\big{)}^{\frac{1}{2}}.$

######  Remark 4.2.

According to Definition 4.1, in the region that is away from the corners, the
weighted space $\mathcal{K}^{m}_{\mathbf{a}}$ is equivalent to the Sobolev
space $H^{m}$. In the neighborhood of $Q_{i}$, the space
$\mathcal{K}^{m}_{\mathbf{a}}(B_{i})$ is the equivalent to the Kondratiev
space [19, 26, 28],

$\mathcal{K}_{a_{i}}^{m}(B_{i}):=\\{v|\
r_{i}^{|\nu|-a_{i}}\partial^{\alpha}v\in L^{2}(B_{i}),\forall\ |\nu|\leq
m\\},$

where $B_{i}\subset\Omega$ represents the neighborhood of $Q_{i}$ satisfying
$B_{i}\cap B_{j}=\emptyset$ for $i\not=j$.

###  4.2. Graded meshes

We now present the construction of graded meshes to improve the convergence
rate of the numerical approximation from Algorithm 3.2.

######  Algorithm 4.3.

(Graded refinements) Let $\mathcal{T}$ be a triangulation of $\Omega$ with
shape-regular triangles. Recall that $Q_{i}$, $i=1,\cdots,N$ are the vertices
of $\Omega$. Let ${AB}$ be an edge in the triangulation $\mathcal{T}$ with $A$
and $B$ as the endpoints. Then, in a graded refinement, a new node $D$ on $AB$
is produced according to the following conditions:

  * 1.

(Neither $A$ nor $B$ coincides with $Q_{i}$.) We choose $D$ as the midpoint
($|AD|=|BD|$).

  * 2.

($A$ coincides with $Q_{i}$.) We choose $r$ such that
$|AD|=\kappa_{Q_{i}}|AB|$, where $\kappa_{Q_{i}}\in(0,0.5)$ is a parameter
that will be specified later. See Figure 3 for example.

Then, the graded refinement, denoted by $\kappa(\mathcal{T})$, proceeds as
follows. For each triangle $T\in\mathcal{T}$, a new node is generated on each
edge of $T$ as described above. Then, $T$ is decomposed into four small
triangles by connecting these new nodes (Figure 4). Given an initial mesh
$\mathcal{T}_{0}$ satisfying the condition above, the associated family of
graded meshes $\\{\mathcal{T}_{n},\ n\geq 0\\}$ is defined recursively
$\mathcal{T}_{n+1}=\kappa(\mathcal{T}_{n})$.

$B$$A$$D$$B$$A$$D$ Figure 3. The new node on an edge $AB$. (left): $A\neq
Q_{i}$ and $B\neq Q_{i}$ (midpoint); (right) $A=Q_{i}$
($|AD|=\kappa_{Q_{i}}|AB|$, $\kappa_{Q_{i}}<0.5$).
$x_{2}$$x_{0}$$x_{1}$$x_{2}$$x_{0}$$x_{1}$$x_{12}$$x_{01}$$x_{02}$$x_{2}$$x_{0}$$x_{1}$$x_{12}$$x_{01}$$x_{02}$$x_{2}$$x_{0}$$x_{1}$$x_{12}$$x_{01}$$x_{02}$
Figure 4. Refinement of a triangle $\triangle x_{0}x_{1}x_{2}$. First row:
(left – right): the initial triangle and the midpoint refinement; second row:
two consecutive graded refinements toward $x_{0}=Q_{i}$,
($\kappa_{Q_{i}}<0.5$).

Given a grading parameter $\kappa_{Q_{i}}$, Algorithm 4.3 produces smaller
elements near $Q_{i}$ for a better approximation of the singular solution. It
is an explicit construction of graded meshes based on recursive refinements.
See also [2, 7, 31, 33] and references therein for more discussions on the
graded mesh.

Note that after $n$ refinements, the number of triangles in the mesh
$\mathcal{T}_{n}$ is $O(4^{n})$, so we denote the “mesh size” of
$\mathcal{T}_{n}$ by

$h=2^{-n}.$ (4.2)

In Algorithm 4.3, we choose the parameter $\kappa_{Q_{i}}$ for each vertex
$Q_{i}$ as follows. Recall that and $\alpha_{0}^{i}$ is the solution of (2.2)
with $\omega$ being replaced by the interior angle $\omega_{i}$ at $Q_{i}$.
Given the degree of polynomials $k,k^{\prime}$ in Algorithm 3.2, we choose

$\kappa_{Q_{i}}=2^{-\frac{\theta}{a_{i}}}\left(\leq\frac{1}{2}\right),$ (4.3)

where $a_{i}>0$ and $\theta$ could be any possible constants satisfying

$a_{i}\leq\theta\leq\min\\{k,m\\}.$ (4.4)

In (4.4), if we take $a_{i}=\theta$, the grading parameter
$\kappa_{Q_{i}}=\frac{1}{2}$.

###  4.3. Interpolation error estimates on graded meshes

######  Lemma 4.4.

Let $T_{(0)}\in\mathcal{T}_{0}$ be an initial triangle of the triangulation
$\mathcal{T}_{n}$ in Algorithm 4.3 with grading parameters $\kappa_{Q_{i}}$
given by (4.3). For $m\geq 1,k\geq 1$, we denote $v_{I}\in V_{n}^{k}$ (resp.
$q_{I}\in S_{n}^{k-1}$) the nodal interpolation of
$v\in\mathcal{K}_{\mathbf{a}+1}^{m+1}(\Omega)$ (resp.
$q\in\mathcal{K}_{\mathbf{a}}^{m}(\Omega)$). If $\bar{T}_{(0)}$ does not
contain any vertices $Q_{i}$, $i=1,\cdots,N$, then

$\|v-v_{I}\|_{H^{1}(T_{(0)})}\leq
Ch^{\min\\{k,m\\}},\quad\|q-q_{I}\|_{L^{2}(T_{(0)})}\leq Ch^{\min\\{k,m\\}},$

where $h=2^{-n}$.

###### Proof.

If $\bar{T}_{0}$ does not contain any vertices $Q_{i}$ of the domain $\Omega$,
we have $v\in\mathcal{K}_{\mathbf{a}+1}^{m+1}(\Omega)\subset H^{m+1}(T_{(0)})$
(see Remark 4.2) and the mesh on $T_{(0)}$ is quasi-uniform (Algorithm 4.3)
with size $O(2^{-n})$. Therefore, based on the standard interpolation error
estimate, we have

$\displaystyle\|v-v_{I}\|_{H^{1}(T_{(0)})}\leq
Ch^{\min\\{k,m\\}}\|v\|_{H^{m+1}(T_{(0)})}.$ (4.5)

Note that $q\in\mathcal{K}_{\mathbf{a}}^{m}(\Omega)\subset H^{m}(\Omega)$, we
can similar obtain the estimate for $\|q-q_{I}\|$. ∎

We now study the interpolation error in the neighborhood $Q_{i}$,
$i=1,\cdots,N$. In the rest of this subsection, we assume
$T_{(0)}\in\mathcal{T}_{0}$ is an initial triangle such that the $i$th vertex
$Q_{i}$ is a vertex of $T_{(0)}$. We first define mesh layers on $T_{(0)}$
which are collections of triangles in $\mathcal{T}_{n}$.

######  Definition 4.5.

(Mesh layers) Let $T_{(t)}\subset T_{(0)}$ be the triangle in
$\mathcal{T}_{t}$, $0\leq t\leq n$, that is attached to the singular vertex
$Q_{i}$ of $T_{(0)}$. For $0\leq t<n$, we define the $t$th mesh layer of
$\mathcal{T}_{n}$ on $T_{(0)}$ to be the region $L_{t}:=T_{(t)}\setminus
T_{(t+1)}$; and for $t=n$, the $n$th layer is $L_{n}:=T_{(n)}$. See Figure 5
for example.

$Q_{i}$$Q_{i}$$L_{0}$$L_{1}$$Q_{i}$$L_{0}$$L_{1}$$L_{2}$ Figure 5. Mesh layers
(left – right): the initial triangle $T_{(0)}$ with a vertex $Q_{i}$; two
layers after one refinement; three layers after two refinements.

######  Remark 4.6.

The triangles in $\mathcal{T}_{n}$ constitute $n$ mesh layers on $T_{(0)}$.
According to Algorithm 4.3 and the choice of grading parameters
$\kappa_{Q_{i}}$ given by (4.3), the mesh size in the $t$th layer $L_{t}$ is

$O(\kappa_{Q_{i}}^{t}2^{t-n}).$ (4.6)

Meanwhile, the weight function $\rho$ in (4.1) satisfies

$\displaystyle\rho=O(\kappa_{Q_{i}}^{t})\ \ \ {\rm{in\ }}L_{t}\ (0\leq
t<n)\qquad{\rm{and}}\qquad\rho\leq C\kappa_{Q_{i}}^{n}\ \ \ {\rm{in\ }}L_{n}.$
(4.7)

Although the mesh size varies in different layers, the triangles in
$\mathcal{T}_{n}$ are shape regular. In addition, using the local Cartesian
coordinates such that $Q$ is the origin, the mapping

$\displaystyle\mathbf{B}_{t}=\begin{pmatrix}\kappa_{Q_{i}}^{-t}&0\\\
0&\kappa_{Q_{i}}^{-t}\\\ \end{pmatrix}\qquad 0\leq t\leq n,$ (4.8)

is a bijection between $L_{t}$ and $L_{0}$ for $0\leq t<n$ and a bijection
between $L_{n}$ and $T_{(0)}$. We call $L_{0}$ (resp. $T_{(0)}$) the reference
region associated to $L_{t}$ for $0\leq t<n$ (resp. $L_{n}$).

With the mapping (4.8), we have that for any point $(x,y)\in L_{t}$, $0\leq
t\leq n$, the image point $(\hat{x},\hat{y}):=\mathbf{B}_{t}(x,y)$ is in its
reference region. We introduce the following dilation result.

######  Lemma 4.7.

For $0\leq t\leq n$, given a function $v(x,y)\in\mathcal{K}_{a}^{l}(L_{t})$,
the function $\hat{v}(\hat{x},\hat{y}):=v(x,y)$ belongs to
$\mathcal{K}_{a}^{l}(\hat{L})$, where
$(\hat{x},\hat{y}):=\mathbf{B}_{t}(x,y)$, $\hat{L}=L_{0}$ for $0\leq t<n$, and
$\hat{L}=T_{(0)}$ for $t=n$. Then, it follows

$\|\hat{v}(\hat{x},\hat{y})\|_{\mathcal{K}_{a}^{l}(\hat{L})}=\kappa_{Q_{i}}^{t(a-1)}\|v(x,y)\|_{\mathcal{K}_{a}^{l}(L_{i})}.$

###### Proof.

The proof can be found in [34, Lemma 4.5]. ∎

We then derive the interpolation error estimate in each layer.

######  Lemma 4.8.

For $k\geq 1,m\geq 1$, set $\kappa_{Q_{i}}$ in (4.3) with $\theta$ satisfying
(4.4) for the graded mesh on $T_{(0)}$. Let $h:=2^{-n}$, then in the $t$th
layer $L_{t}$ on $T_{(0)}$, $0\leq t<n$,
(i) if $v_{I}\in V_{n}^{k}$ be the nodal interpolation of
$v\in\mathcal{K}_{\mathbf{a}+1}^{m+1}(\Omega)$, it follows

$|v-v_{I}|_{H^{1}(L_{t})}\leq
Ch^{\theta}\|v\|_{\mathcal{K}_{a_{i}+1}^{m+1}(L_{t})};$ (4.9)

(ii) if $q_{I}\in V_{n}^{k-1}$ be the nodal interpolation of
$q\in\mathcal{K}_{\mathbf{a}}^{m}(\Omega)$, it follows

$\|q-q_{I}\|_{L^{2}(L_{t})}\leq
Ch^{\theta}\|q\|_{\mathcal{K}_{a_{i}}^{m}(L_{t})}.$ (4.10)

###### Proof.

For $L_{t}$ associated with $Q_{i}$, $0\leq t<n$, the space
$\mathcal{K}_{a_{i}+1}^{m+1}(L_{t})$ (resp. $\mathcal{K}_{a_{i}}^{m}(L_{t})$)
is equivalent to $H^{m+1}(L_{t})$ (resp. $H^{m}(L_{t})$). Therefore, $v$
(resp. $q$) is a continuous function in $L_{t}$. For any point $(x,y)\in
L_{t}$, let $(\hat{x},\hat{y})=\mathbf{B}_{t}(x,y)\in L_{0}$. For $v(x,y)$
(resp. $q(x,y)$) in $L_{t}$, we define $\hat{v}(\hat{x},\hat{y}):=v(x,y)$
(resp. $\hat{q}(\hat{x},\hat{y}):=q(x,y)$) in $L_{0}$.

(i). Using the standard interpolation error estimate, the scaling argument,
the estimate in (4.6), and the mapping in (4.8), we have

$\displaystyle|v-v_{I}|_{H^{1}(L_{t})}$ $\displaystyle=$
$\displaystyle|\hat{v}-\hat{v}_{I}|_{H^{1}(L_{0})}\leq
C2^{(t-n)\mu}\|\hat{v}\|_{\mathcal{K}_{a_{i}+1}^{m+1}(L_{0})}\leq
C2^{(t-n)\mu}\kappa_{Q_{i}}^{a_{i}t}\|v\|_{\mathcal{K}_{a_{i}+1}^{m+1}(L_{t})},$

where we have used Lemma 4.7 in the last inequality. Since
$\kappa_{Q_{i}}=2^{-\frac{\theta}{a_{i}}}$, so we have
$\kappa_{Q_{i}}^{a_{i}t}=2^{-\theta t}$. Set $\mu=\min\\{k,m\\}$, by
$\theta\leq\mu$ from (4.4) and $t<n$, we have $2^{(n-t)(\theta-\mu)}<2^{0}=1$.
Therefore, we have the estimate

$\displaystyle|v-v_{I}|_{H^{1}(L_{t})}$ $\displaystyle\leq$ $\displaystyle
C2^{(t-n)\mu-\theta
t}\|v\|_{\mathcal{K}_{a_{i}+1}^{m+1}(L_{t})}=C2^{-n\theta}2^{(n-t)(\theta-\mu)}\|v\|_{\mathcal{K}_{a_{i}+1}^{m+1}(L_{t})}$
$\displaystyle\leq$ $\displaystyle
C2^{-n\theta}\|v\|_{\mathcal{K}_{a_{i}+1}^{m+1}(L_{t})}\leq
Ch^{\theta}\|v\|_{\mathcal{K}_{a_{i}+1}^{m+1}(L_{t})}.$

(ii). We can show that

$\displaystyle\|q-q_{I}\|_{L^{2}(L_{t})}$ $\displaystyle=$
$\displaystyle\kappa_{Q_{i}}^{t}\|\hat{q}-\hat{q}_{I}\|_{L^{2}(L_{0})}\leq
C\kappa_{Q_{i}}^{t}2^{(t-n)\mu}\|\hat{q}\|_{\mathcal{K}_{a_{i}}^{m}(L_{0})}\leq
C2^{(t-n)\mu}\kappa_{Q_{i}}^{a_{i}t}\|q\|_{\mathcal{K}_{a_{i}}^{m}(L_{t})},$

where again we used Lemma 4.7 in the last inequality. Using the similar
argument as in (i), the estimate (4.10) holds. ∎

Before deriving the interpolation error estimate in the last layer $L_{n}$ on
$T_{(0)}$, we first introduce the following results.

######  Lemma 4.9.

For $\forall v\in\mathcal{K}_{a}^{l}(L_{n})$, if $0\leq l^{\prime}\leq l$ and
$a^{\prime}\leq a$, then it follows

$\|v\|_{\mathcal{K}_{a^{\prime}}^{l^{\prime}}(L_{n})}\leq
C\kappa_{Q_{i}}^{n(a-a^{\prime})}\|v\|_{\mathcal{K}_{a}^{l}(L_{n})}.$ (4.11)

###### Proof.

This is a direct application of [33, Lemma 2.6] under condition (4.7) on
$L_{n}$. ∎

######  Lemma 4.10.

For $\forall v\in\mathcal{K}_{a}^{l}(L_{n})$ , if $a\geq l$, then it follows
that

$\|v\|_{H^{l}(L_{n})}\leq
C\kappa_{Q_{i}}^{n(a-l)}\|v\|_{\mathcal{K}_{a}^{l}(L^{n})}.$ (4.12)

###### Proof.

This is a direct application of [33, Lemma 2.8] under condition (4.7) on
$L_{n}$. ∎

######  Lemma 4.11.

For $k\geq 1,m\geq 1$, set $\kappa_{Q_{i}}$ in (4.3) with $\theta$ satisfying
(4.4) for the graded mesh on $T_{(0)}$. Let $h:=2^{-n}$, then in the $n$th
layer $L_{n}$ on $T_{(0)}$ for $n$ sufficiently large,
(i) if $v_{I}\in V_{n}^{k}$ be the nodal interpolation of
$v\in\mathcal{K}_{\mathbf{a}+1}^{m+1}(\Omega)$, it follows

$|v-v_{I}|_{H^{1}(L_{n})}\leq
Ch^{\theta}\|v\|_{\mathcal{K}^{m+1}_{{a_{i}}+1}(L_{n})};$ (4.13)

(ii) if $q_{I}\in V_{n}^{k-1}$ be the nodal interpolation of
$q\in\mathcal{K}_{\mathbf{a}}^{m}(\Omega)$, it follows

$\|q-q_{I}\|_{L^{2}(L_{n})}\leq
Ch^{\theta}\|q\|_{\mathcal{K}_{a_{i}}^{m}(L_{n})}.$ (4.14)

###### Proof.

Recall the mapping $\mathbf{B}_{n}$ in (4.8). For any point $(x,y)\in L_{n}$,
let $(\hat{x},\hat{y})=\mathbf{B}_{n}(x,y)\in T_{(0)}$.

(i). Let $\eta:T_{(0)}\rightarrow[0,1]$ be a smooth function that is equal to
$0$ in a neighborhood of $Q_{i}$, but is equal to 1 at all the other nodal
points in $\mathcal{T}_{0}$. For a function $v(x,y)$ in $L_{n}$, we define
$\hat{v}(\hat{x},\hat{y}):=v(x,y)$ in $T_{(0)}$. We take $w=\eta\hat{v}$ in
$T_{(0)}$. Consequently, we have for $l\geq 0$

$\displaystyle\|w\|^{2}_{\mathcal{K}^{l}_{1}(T_{(0)})}$ $\displaystyle=$
$\displaystyle\|\eta\hat{v}\|^{2}_{\mathcal{K}^{l}_{1}(T_{(0)})}\leq
C\|\hat{v}\|^{2}_{\mathcal{K}^{l}_{1}(T_{(0)})},$ (4.15)

where $C$ depends on $l$ and the smooth function $\eta$. Moreover, the
condition $\hat{v}\in\mathcal{K}_{a^{i}+1}^{m+1}(T_{(0)})$ with and $m\geq 2$
implies $\hat{v}(Q)=0$ (see, e.g., [33, Lemma 4.7]). Let $w_{\hat{I}}$ be the
nodal interpolation of $w$ associated with the mesh $\mathcal{T}_{0}$ on
$T_{(0)}$. Therefore, by the definition of $w$, we have

$\displaystyle w_{\hat{I}}=\hat{v}_{\hat{I}}=\widehat{v_{I}}\quad{\rm{in}}\
T_{(0)}.$ (4.16)

Note that the $\mathcal{K}^{l}_{1}$ norm and the $H^{l}$ norm are equivalent
for $w$ on $T_{(0)}$, since $w=0$ in the neighborhood of the vertex $Q_{i}$.
Let $r$ be the distance from $(x,y)$ to $Q_{i}$, and $\hat{r}$ be the distance
from $(\hat{x},\hat{y})$ to $Q_{i}$. Then, by the definition of the weighted
space, the scaling argument, (4.15), (4.16), and (4.7), we have

$\displaystyle|v-v_{I}|_{H^{1}(L_{n})}^{2}$ $\displaystyle\leq$ $\displaystyle
C\|v-v_{I}\|_{\mathcal{K}^{1}_{{1}}(L_{n})}^{2}\leq C\sum_{|\nu|\leq
1}\|r(x,y)^{|\nu|-1}\partial^{\nu}(v-v_{I})\|_{L^{2}(L_{n})}^{2}$
$\displaystyle\leq$ $\displaystyle C\sum_{|\nu|\leq
1}\|\hat{r}(\hat{x},\hat{y})^{|\nu|-1}\partial^{\nu}(\hat{v}-\widehat{v_{I}})\|_{L^{2}(T_{(0)})}^{2}\leq
C\|\hat{v}-w+w-\widehat{v_{I}}\|_{\mathcal{K}^{1}_{1}(T_{(0)})}^{2}$
$\displaystyle\leq$ $\displaystyle
C\big{(}\|\hat{v}-w\|^{2}_{\mathcal{K}^{1}_{1}(T_{(0)})}+\|w-\widehat{v_{I}}\|^{2}_{\mathcal{K}^{1}_{1}(T_{(0)})}\big{)}=C\big{(}\|\hat{v}-w\|^{2}_{\mathcal{K}^{1}_{1}(T_{(0)})}+\|w-w_{\hat{I}}\|^{2}_{\mathcal{K}^{1}_{1}(T_{(0)})}\big{)}$
$\displaystyle\leq$ $\displaystyle
C\big{(}\|\hat{v}\|^{2}_{\mathcal{K}^{1}_{1}(T_{(0)})}+\|w\|^{2}_{\mathcal{K}^{m+1}_{1}(T_{(0)})}\big{)}\leq
C\big{(}\|\hat{v}\|^{2}_{\mathcal{K}^{1}_{1}(T_{(0)})}+\|\hat{v}\|^{2}_{\mathcal{K}^{m+1}_{1}(T_{(0)})}\big{)}$
$\displaystyle=$ $\displaystyle
C\big{(}\|v\|^{2}_{\mathcal{K}^{1}_{1}(L_{n})}+\|v\|^{2}_{\mathcal{K}^{m+1}_{1}(L_{n})}\big{)}\leq
C\kappa_{Q_{i}}^{2na_{i}}\|v\|_{\mathcal{K}^{m+1}_{{a_{i}}+1}(L_{n})}^{2}$
$\displaystyle\leq$ $\displaystyle
C2^{-2n\theta}\|v\|_{\mathcal{K}^{m+1}_{{a_{i}}+1}(L_{n})}^{2}\leq
Ch^{2\theta}\|v\|_{\mathcal{K}^{m+1}_{{a_{i}}+1}(L_{n})}^{2},$

where the ninth and tenth relationships are based on Lemma 4.7 and Lemma 4.9,
respectively. This completes the proof of (4.9).

(ii). Since $q\in L^{2}(\Omega)$, we have that the interpolation operator is
$L^{2}$ stable

$\|q_{I}\|_{L^{2}(L_{n})}\leq C\|q\|_{L^{2}(L_{n})}.$ (4.17)

Thus, by (4.17) and (4.12), we have

$\displaystyle\|q-q_{I}\|_{L^{2}(L_{n})}\leq C\|q\|_{L^{2}(L_{n})}\leq
C\kappa_{Q_{i}}^{na_{i}}\|q\|_{\mathcal{K}_{a_{i}}^{m}(L_{n})}\leq
C2^{-n\theta}\|q\|_{\mathcal{K}_{a_{i}}^{m}(L_{n})}\leq
Ch^{\theta}\|q\|_{\mathcal{K}_{a_{i}}^{m}(L_{n})}.$

∎

######  Theorem 4.12.

Let $\mathcal{T}_{0}$ be an initial triangle of the triangulation
$\mathcal{T}_{n}$ in Algorithm 4.3 with grading parameters $\kappa_{Q_{i}}$ in
(4.3). For $k\geq 1,m\geq 1$, if $v_{I}\in V_{n}^{k}$ (resp. $q_{I}\in
V_{n}^{k-1}$) be the nodal interpolation of
$v\in\mathcal{K}_{\mathbf{a}+1}^{m+1}(\Omega)$ (resp.
$q\in\mathcal{K}_{\mathbf{a}}^{m}(\Omega)$). Then, it follows the following
interpolation error

$\|v-v_{I}\|_{H^{1}(\Omega)}\leq
Ch^{\theta}\|v\|_{\mathcal{K}_{\mathbf{a}+1}^{m+1}(\Omega)},\quad\|q-q_{I}\|\leq
Ch^{\theta}\|q\|_{\mathcal{K}_{\mathbf{a}}^{m}(\Omega)},$ (4.18)

where $h:=2^{-n}$, and $\theta$ satisfying (4.4).

###### Proof.

By summing the estimates in Lemma 4.4, Lemma 4.8, and Lemma 4.11, we have

$\displaystyle\|v-v_{I}\|^{2}_{H^{1}(\Omega)}=$
$\displaystyle\sum_{T_{(0)}\in\mathcal{T}_{0}}\|v-v_{I}\|^{2}_{H^{1}(T_{(0)})}\leq
Ch^{2\theta}\|v\|^{2}_{\mathcal{K}_{\mathbf{a}+1}^{m+1}(\Omega)},$
$\displaystyle\|q-q_{I}\|^{2}=$
$\displaystyle\sum_{T_{(0)}\in\mathcal{T}_{0}}\|v-v_{I}\|^{2}_{L^{2}(T_{(0)})}\leq
Ch^{2\theta}\|q\|^{2}_{\mathcal{K}_{\mathbf{a}}^{m}(\Omega)}.$

∎

###  4.4. Regularity in weighted Sobolev space

Let $\alpha_{0}^{i}$ the solution of (2.2) with $\omega$ being replaced by the
interior angle $\omega_{i}$ at $Q_{i}$. We denote the vector
$\mbox{\boldmath$\alpha$}_{0}=(\alpha_{0}^{1},\cdots,\alpha_{0}^{i},\cdots,\alpha_{0}^{N})$.
By Lemma 2.2, for $f\in H^{-1}(\Omega)$, there exists
$\mathbf{F}\in[\mathcal{K}_{\mathbf{a}-1}^{0}(\Omega)]^{2}$ with
$a_{i}<\min_{i}\\{\alpha_{0}^{i}\\}$ satisfying (2.8) and

$\|\mathbf{F}\|_{[\mathcal{K}_{\mathbf{a}-1}^{0}(\Omega)]^{2}}\leq
C\|f\|_{H^{-1}(\Omega)}.$ (4.19)

If $f\in\mathcal{K}_{\mathbf{a}-2}^{m-2}(\Omega)\cap H^{-1}(\Omega)$ with
$m\geq 1$ and $0\leq\mathbf{a}<\mbox{\boldmath$\alpha$}_{0}$, then we can find
$\mathbf{F}\in[\mathcal{K}_{\mathbf{a}-1}^{m-1}(\Omega)]^{2}$ satisfying (2.8)
and

$\|\mathbf{F}\|_{[\mathcal{K}_{\mathbf{a}-1}^{m-1}(\Omega)]^{2}}\leq
C\|f\|_{\mathcal{K}_{\mathbf{a}-2}^{m-2}(\Omega)}.$ (4.20)

For the Stokes problem (2.7), we have the following regularity estimate in
weighted Sobolev space [8].

######  Lemma 4.13.

Let $(\mathbf{u},p)\in[H_{0}^{1}(\Omega)]^{2}\times L_{0}^{2}(\Omega)$ be the
solution of the Stokes problem (2.7). For $m\geq 1$ and
$0\leq\mathbf{a}\leq\mbox{\boldmath$\alpha$}_{0}$, if
$\mathbf{F}\in[\mathcal{K}_{\mathbf{a}-1}^{m-1}(\Omega)]^{2}$, then it follows

$\|\mathbf{u}\|_{[\mathcal{K}_{\mathbf{a}+1}^{m+1}(\Omega)]^{2}}+\|p\|_{\mathcal{K}_{\mathbf{a}}^{m}(\Omega)}\leq
C\|\mathbf{F}\|_{[\mathcal{K}_{\mathbf{a}-1}^{m-1}(\Omega)]^{2}}.$ (4.21)

We then have the following result.

######  Lemma 4.14.

Given $f\in\mathcal{K}_{\mathbf{a}-2}^{m-2}(\Omega)\cap H^{-1}(\Omega)$ for
$0\leq\mathbf{a}<\mbox{\boldmath$\alpha$}_{0}$ and $m\geq 1$, let $\phi\in
H_{0}^{2}(\Omega)$ be the solution of the Possible problem (2.40), then it
follows

$\|\phi\|_{\mathcal{K}_{\mathbf{a}+2}^{m+2}(\Omega)}\leq
C\|f\|_{\mathcal{K}_{\mathbf{a}-2}^{m-2}(\Omega)}.$ (4.22)

###### Proof.

By Lemma 4.13, we have
$\mathbf{u}\in[\mathcal{K}_{\mathbf{a}+1}^{m+1}(\Omega)]^{2}\cap[H_{0}^{1}(\Omega)]^{2}$,
thus we have $(\text{curl }\mathbf{u})\in\mathcal{K}_{\mathbf{a}}^{m}(\Omega)$
and

$\|\text{curl }\mathbf{u}\|_{\mathcal{K}_{\mathbf{a}}^{m}(\Omega)}\leq
C\|\nabla\mathbf{u}\|_{[\mathcal{K}_{\mathbf{a}}^{m}(\Omega)]^{2}}\leq
C\|\mathbf{u}\|_{[\mathcal{K}_{\mathbf{a}+1}^{m+1}(\Omega)]^{2}}\leq
C\|\mathbf{F}\|_{[\mathcal{K}_{\mathbf{a}-1}^{m-1}(\Omega)]^{2}}.$ (4.23)

By the regularity estimate [26, 31] for the Poisson problem (2.40), we have

$\|\phi\|_{\mathcal{K}_{\mathbf{a}+2}^{m+2}(\Omega)}\leq\|\text{curl
}\mathbf{u}\|_{\mathcal{K}_{\mathbf{a}}^{m}(\Omega)}.$ (4.24)

The conclusion holds by combining (4.20), (4.23) and (4.24). ∎

###  4.5. Optimal error estimates on graded meshes

To better observe the threshold of grading parameter $\kappa_{Q_{i}}$ in
obtaining the optimal convergence rates, we always assume $1\leq k\leq m$ in
the following discussions, otherwise, we just replace $k$ by $\min\\{k,m\\}$.
In this section, we assume that
$f\in\mathcal{K}_{\mathbf{a}-1}^{m-1}(\Omega)\cap\mathcal{K}_{\mathbf{b}-1}^{m-1}(\Omega)$
with $0<\mathbf{a}<\mbox{\boldmath$\alpha$}_{0}$, and
$0<\mathbf{b}<\mbox{\boldmath$\beta$}_{0}$, where
$\mbox{\boldmath$\beta$}_{0}=(\frac{\pi}{\omega_{1}},\cdots,\frac{\pi}{\omega_{N}})$.

If $\mathbf{F}$ is given by Lemma 2.3, then we have
$\mathbf{F}\in[\mathcal{K}_{\mathbf{a}-1}^{m-1}(\Omega)]^{2}$, and the
regularities in Lemma 4.13 and Lemma 4.14 hold.

If $\mathbf{F}$ is given by Lemma 2.2, by the regularity estimate [7] for the
Poisson problem (2.13) on weighted Sobolev space, it follows that

$\|w\|_{\mathcal{K}_{\mathbf{b}+1}^{m+1}(\Omega)}\leq
C\|f\|_{\mathcal{K}_{\mathbf{b}-1}^{m-1}(\Omega)},$ (4.25)

which implies $\mathbf{F}=\textbf{curl
}w\in[\mathcal{K}_{\mathbf{b}}^{m}(\Omega)]^{2}\subset[\mathcal{K}_{\mathbf{b}}^{m-1}(\Omega)]^{2}$.
Then the solution of the Stokes problem (2.7) satisfies

$\|\mathbf{u}\|_{[\mathcal{K}_{\mathbf{c}+1}^{m+1}(\Omega)]^{2}}+\|p\|_{\mathcal{K}_{\mathbf{c}}^{m}(\Omega)}\leq
C\|\mathbf{F}\|_{[\mathcal{K}_{{\mathbf{c}}-1}^{m-1}(\Omega)]^{2}}\leq
C\|\mathbf{F}\|_{[\mathcal{K}_{{\mathbf{b}}-1}^{m-1}(\Omega)]^{2}},$ (4.26)

where $\mathbf{c}=(c_{1},\cdots,c_{N})$ with $c_{i}=\min\\{b_{i}+1,a_{i}\\}$.
Thus, the solution of the Poisson problem (2.40) satisfies
$\phi\in\mathcal{K}_{\mathbf{c}+2}^{m+2}(\Omega)$.

Since the bilinear functional in (3.3) is coercive and continuous on
$V_{n}^{k}$, so we have by Céa’s Theorem,

$\|w-w_{n}\|_{H^{1}(\Omega)}\leq C\inf_{v\in
V_{n}^{k}}\|w-v\|_{H^{1}(\Omega)}.$ (4.27)

Recall that $\alpha_{0}=\min_{i}\\{\alpha_{0}^{i}\\}$ given by (2.2), and
$\beta_{0}=\min_{i}\\{\beta_{0}^{i}\\}=\frac{\pi}{\omega}$ are the thresholds
corresponding to the largest interior angle $\omega$, then we have the
following result.

######  Lemma 4.15.

Set the grading parameters
$\kappa_{Q_{i}}=2^{-\frac{\theta}{a_{i}}}\left(=2^{-\frac{\theta^{\prime}}{b_{i}}}\right)$
with $0<a_{i}<\alpha_{0}^{i}$, $0<b_{i}<\beta_{0}^{i}$, $\theta$ being any
constant satisfying $a_{i}\leq\theta\leq k$, and
$\theta^{\prime}=\min\left\\{\frac{\beta_{0}}{\alpha_{0}}\max\\{\theta,\alpha_{0}\\},k\right\\}$
satisfying $b_{i}\leq\theta^{\prime}\leq k$. Let $w_{n}\in V_{n}^{k}$ be the
solution of finite element solution of (3.3), and $w$ is the solution of the
Poisson problem (2.13), then it follows

$\|w-w_{n}\|_{H^{1}(\Omega)}\leq Ch^{\theta^{\prime}},\quad\|w-w_{n}\|\leq
Ch^{\min\left\\{2\theta^{\prime},\theta^{\prime}+1\right\\}},$ (4.28)

where $h:=2^{-n}$.

###### Proof.

By (4.27) and the interpolation error estimates in Lemma 4.12 under the
regularity result in (4.25) and
$\kappa_{Q_{i}}=2^{-\frac{\theta}{a_{i}}}=2^{-\frac{\theta^{\prime}}{b_{i}}}$,
we have the estimate

$\|w-w_{n}\|_{H^{1}(\Omega)}\leq C\|w-w_{I}\|_{H^{1}(\Omega)}\leq
Ch^{\theta^{\prime}}.$

Consider the Poisson problem

$-\Delta v=w-w_{n}\text{ in }\Omega,\quad v=0\text{ on }\partial\Omega.$
(4.29)

Then we have

$\|w-w_{n}\|^{2}=(\nabla(w-w_{n}),\nabla v).$ (4.30)

Subtract (3.3) from (2.13), we have the Galerkin orthogonality,

$(\nabla(w-w_{n}),\nabla\psi)=0\quad\forall\psi\in V_{n}^{k}.$ (4.31)

Setting $\psi=v_{I}\in V_{n}^{k}$ the nodal interpolation of $v$ and subtract
(4.31) from (4.30), we have

$\displaystyle\|w-w_{n}\|^{2}=(\nabla(w-w_{n}),\nabla(v-v_{I}))\leq\|w-w_{n}\|_{H^{1}(\Omega)}\|v-v_{I}\|_{H^{1}(\Omega)}.$
(4.32)

Similarly, the solution $v\in\mathcal{K}^{2}_{\mathbf{b}^{\prime}+1}(\Omega)$
satisfies the regularity estimate

$\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}\leq
C\|w-w_{n}\|_{K^{0}_{\mathbf{b}^{\prime}-1}(\Omega)}\leq C\|w-w_{n}\|,$ (4.33)

where the $i$th entry of $\mathbf{b}^{\prime}$ satisfying
$b^{\prime}_{i}=\min\left\\{b_{i},1\right\\}$. By Lemma 4.12 with grading
parameter
$\kappa_{Q_{i}}=2^{-\frac{\theta^{\prime}}{b_{i}}}(=2^{-\frac{\theta}{a_{i}}})$
again, we have the interpolation error

$\|v-v_{I}\|_{H^{1}(\Omega)}\leq
Ch^{\min\\{\theta^{\prime},1\\}}\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}.$
(4.34)

The $L^{2}$ error estimate in (4.28) can be obtained by combining (4.32),
(4.33), and (4.34). ∎

Thus, we have the following result,

$\displaystyle\|\mathbf{F}-\mathbf{F}_{n}\|_{[L^{2}(\Omega)]^{2}}=$
$\displaystyle\|\textbf{curl }v-\textbf{curl
}v_{n}\|_{[L^{2}(\Omega)]^{2}}\leq\|w-w_{n}\|_{H^{1}(\Omega)}\leq
Ch^{\theta^{\prime}},$ (4.35)
$\displaystyle\|\mathbf{F}-\mathbf{F}_{n}\|_{[H^{-1}(\Omega)]^{2}}=$
$\displaystyle\|\textbf{curl }w-\textbf{curl
}w_{n}\|_{[H^{-1}(\Omega)]^{2}}\leq C\|w-w_{n}\|\leq
Ch^{\min\\{\theta^{\prime}+1,2\theta^{\prime}\\}}.$

Now we have the following error estimate of the Mini element approximation
($k=1$) or Taylor-Hood method ($k\geq 2$) in Algorithm 3.2 on graded meshes
for the Stokes problem (2.7).

######  Lemma 4.16.

The bilinear forms in both Mini element method and the Taylor-Hood method on
graded meshes satisfies the LBB or inf-sup condition (3.7) or (3.8).

###### Proof.

For given $\kappa=\min_{i}\\{\kappa_{Q_{i}}\\}$, Algorithm 4.3 implies that
there exists a constant $\sigma(\kappa)>0$ such that

$h_{T}\leq\sigma(\kappa)\rho_{T},\quad\forall T\in\mathcal{T}_{n},$ (4.36)

where $h_{T}$ is the diameter of $T$, and $\rho_{T}$ is the maximum diameter
of all circles contained in $T$. Under condition (4.36) of the graded mesh,
the conclusion follows from [39, Theorem 3.1]. ∎

######  Theorem 4.17.

Set the grading parameters $\kappa_{Q_{i}}=2^{-\frac{\theta}{a_{i}}}$ with
$0<a_{i}<\alpha_{0}^{i}$ and $\theta$ being any constants satisfying
$a_{i}\leq\theta\leq k$. Let $(\mathbf{u},p)$ be the solution of the Stokes
problem (2.9), and $(\mathbf{u}_{n},p_{n})$ be the Mini element solution
($k=1$) or Taylor-Hood element solution ($k\geq 2$) on graded meshes
$\mathcal{T}_{n}$. If $(\mathbf{u}_{n},p_{n})$ is the solution of in Algorithm
3.2, then it follows

$\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{1}(\Omega)]^{2}}+\|p-p_{n}\|\leq
Ch^{\theta};$ (4.37)

if it is the solution of in Algorithm 3.1, then it follows

$\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{1}(\Omega)]^{2}}+\|p-p_{n}\|\leq
Ch^{\min\\{\theta,\theta^{\prime}+1\\}},$ (4.38)

where $\theta^{\prime}$ is given in Lemma 4.15.

###### Proof.

For Algorithm 3.2, by Corollary A.2 and the interpolation error estimates in
Lemma 4.12 under the regularity result in Lemma 4.13, the estimate (4.37)
holds.

For Algorithm 3.1, by Lemma A.1 with the estimate (4.35) and the interpolation
error estimates in Lemma 4.12 under the regularity result (4.26), it follows

$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{1}(\Omega)]^{2}}+\|p-p_{n}\|\leq
Ch^{\theta}+Ch^{\min\\{\theta^{\prime}+1,2\theta^{\prime}\\}}\leq
Ch^{\min\\{\theta,\theta^{\prime}+1\\}}.$

Here, we have used the fact that if $\omega>\pi$,
$\theta\leq\theta^{\prime}<2\theta^{\prime}$. Note that
$\theta^{\prime}=\min\left\\{\frac{\beta_{0}}{\alpha_{0}}\max\\{\theta,\alpha_{0}\\},k\right\\}$,
so if $\theta^{\prime}=k$, then $\theta\leq k=\theta^{\prime}$; otherwise

$\theta^{\prime}=\frac{\beta_{0}}{\alpha_{0}}\max\\{\alpha_{0},\theta\\}\geq\frac{\beta_{0}}{\alpha_{0}}\theta>\theta,$

where we have used (2.6).

If $\omega<\pi$, by taking
$1<b_{i}=\frac{\beta_{0}}{\alpha_{0}}a_{i}<\beta_{0}\leq\beta_{0}^{i}$, it
follows $\theta^{\prime}=\frac{b_{i}}{a_{i}}\theta\geq b_{i}>1$, so that
$\theta^{\prime}+1<2\theta^{\prime}$. Thus, the estimate (4.38) holds. ∎

In weighted Sobolev space, the regularity result for (3.19) with $l=0,1$ has
the form

$\|\mathbf{r}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}+1}^{l+2}(\Omega)]^{2}}+\|s\|_{\mathcal{K}_{\mathbf{a}^{\prime}}^{l+1}(\Omega)}\leq
C\|\mathbf{g}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}-1}^{l}(\Omega)]^{2}}\leq
C\|\mathbf{g}\|_{[\mathcal{K}_{\mathbf{b}}^{l}(\Omega)]^{2}},$ (4.39)

where $0<\mathbf{a}^{\prime}=\min\\{\mathbf{a},l+1\\}$ with
$0<\mathbf{a}<\mbox{\boldmath$\alpha$}_{0}$ and
$0<\mathbf{b}<\mbox{\boldmath$\beta$}_{0}$ . Then we have the following
result.

######  Theorem 4.18.

Set the grading parameters $\kappa_{Q_{i}}=2^{-\frac{\theta}{a_{i}}}$ with
$0<a_{i}<\alpha_{0}^{i}$ and $\theta$ being any constants satisfying
$a_{i}\leq\theta\leq k$. Let $(\mathbf{u},p)$ be the solution of the Stokes
problem (2.9), and $(\mathbf{u}_{n},p_{n})$ be the Mini element solution
($k=1$) or Taylor-Hood element solution ($k\geq 2$) in Algorithm 3.1 on graded
meshes $\mathcal{T}_{n}$. If $(\mathbf{u}_{n},p_{n})$ is the solution of in
Algorithm 3.2, then it follows

$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\leq$
$\displaystyle Ch^{\min\\{2\theta,\theta+1\\}},$ (4.40)
$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[(\mathcal{K}_{\mathbf{b}}^{1}(\Omega))^{*}]^{2}}\leq$
$\displaystyle Ch^{\min\\{2\theta,\theta+2\\}};$

if it is the solution of in Algorithm 3.2, then it follows

$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\leq$
$\displaystyle Ch^{\min\\{2\theta,\theta+1,\theta^{\prime}+2\\}},$ (4.41)
$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[(\mathcal{K}_{\mathbf{b}}^{1}(\Omega))^{*}]^{2}}\leq$
$\displaystyle
Ch^{\min\\{2\theta,\theta+2,k+\theta^{\prime},\theta+\theta^{\prime}+1,\theta^{\prime}+3\\}},$

where $\theta^{\prime}$ is given in Lemma 4.15. Here, $(\cdot)^{*}$ represents
the dual space.

###### Proof.

We only prove (4.41) for Taylor-Hood method, all other cases can be proved
similarly. Similar to Lemma 3.6, we take
$\mathbf{v}=\mathbf{u}-\mathbf{u}_{n}$, $q=p-p_{n}$ in (3.18), then we have

$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[(\mathcal{K}_{\mathbf{b}}^{l}(\Omega))^{*}]^{2}}=\sup_{g\in(\mathcal{K}_{\mathbf{b}}^{l}(\Omega))]^{2}}\frac{\langle\mathbf{g},\mathbf{u}-\mathbf{u}_{n}\rangle}{\|\mathbf{g}\|_{[(\mathcal{K}_{\mathbf{b}}^{l}(\Omega))]^{2}}},$

where $l=0,1$. Let $(\mathbf{r}_{n},s_{n})$ be the Taylor-Hood solution of
(3.17), then it follows

$\displaystyle\langle\mathbf{g},\mathbf{u}-\mathbf{u}_{n}\rangle=$
$\displaystyle T_{1}+T_{2}+T_{3}+T_{4},$

where $T_{i}$, $i=1,\cdots,4$ have the same expressions as those in Lemma 3.6.
For $\mathbf{r}_{n}$ and $s_{n}$, we have the following estimate in the
weighted Sobolev space

$\displaystyle\|\mathbf{r}-\mathbf{r}_{n}\|_{[H^{1}(\Omega)]^{2}}+\|s-s_{n}\|\leq
C\left(\inf_{\mathbf{r}_{I}\in[V_{n}^{k}(\Omega)]^{2}}\|\mathbf{r}-\mathbf{r}_{I}\|_{[H^{1}(\Omega)]^{2}}+\inf_{s_{I}\in
S_{n}^{k-1}}\|s-s_{I}\|\right)$ (4.42) $\displaystyle\leq
Ch^{\min\\{\theta,l+1\\}}(\|\mathbf{r}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}+1}^{l+2}(\Omega)]^{2}}+\|s\|_{\mathcal{K}_{\mathbf{a}^{\prime}}^{l+1}(\Omega)}).$

Here we have

$\displaystyle|T_{j}|\leq$ $\displaystyle
Ch^{\min\\{\theta,\theta^{\prime}+1\\}+\min\\{\theta,l+1\\}}(\|\mathbf{r}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}+1}^{l+2}(\Omega)]^{2}}+\|s\|_{\mathcal{K}_{\mathbf{a}^{\prime}}^{l+1}(\Omega)})$
$\displaystyle=$ $\displaystyle
Ch^{\min\\{2\theta,\theta+l+1,\theta+\theta^{\prime}+1,\theta^{\prime}+l+2\\}}(\|\mathbf{r}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}+1}^{l+2}(\Omega)]^{2}}+\|s\|_{\mathcal{K}_{\mathbf{a}^{\prime}}^{l+1}(\Omega)}),$

where $j=1,2,3$.

Note that

$\displaystyle
T_{4}=T_{41}+T_{42}=\langle\mathbf{F}-\mathbf{F}_{n},\mathbf{r}_{n}-\mathbf{r}\rangle+\langle\mathbf{F}-\mathbf{F}_{n},\mathbf{r}\rangle.$

By (4.35) and (4.42), we have

$\displaystyle|T_{41}|\leq$
$\displaystyle\|\mathbf{F}-\mathbf{F}_{n}\|_{[H^{-1}(\Omega)]^{2}}\|\mathbf{r}-\mathbf{r}_{n}\|_{[H^{1}(\Omega)]^{2}}\leq
Ch^{\min\\{\theta^{\prime}+1,2\theta^{\prime}\\}+\min\\{\theta,l+1\\}}\|\mathbf{r}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}+1}^{l+2}(\Omega)]^{2}}$
$\displaystyle=$ $\displaystyle
Ch^{\min\\{\theta+\theta^{\prime}+1,\theta^{\prime}+l+2,2\theta^{\prime}+l+1,\theta+2\theta^{\prime}\\}}\|\mathbf{r}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}+1}^{l+2}(\Omega)]^{2}}.$

By Theorem 4.12 for $\psi\in\mathcal{K}_{\mathbf{a}^{\prime}+2}^{l+3}(\Omega)$
satisfying (3.20) with
$\kappa_{Q_{i}}=2^{-\frac{\theta}{a_{i}}}=2^{-\frac{\theta_{1}}{1+a^{\prime}_{i}}}$,
we have

$\|\psi-\psi_{I}\|_{H^{1}(\Omega)}\leq
Ch^{\min\\{k,\theta_{1},l+2\\}}\|\mathbf{r}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}+1}^{l+2}(\Omega)]^{2}},$
(4.43)

where
$\theta_{1}=(1+\frac{1}{a^{\prime}_{i}})\theta\geq(1+\frac{1}{a_{i}})\theta\geq\theta+1$
and $\psi_{I}$ is the nodal interpolation of $\psi$.

By (4.35) and (4.43), we have

$\displaystyle|T_{42}|\leq$
$\displaystyle\|\mathbf{F}-\mathbf{F}_{n}\|_{[L^{2}(\Omega)]^{2}}\|\psi-\psi_{I}\|_{H^{1}(\Omega)}\leq
Ch^{\theta^{\prime}+\min\\{k,\theta+\theta^{\prime}_{1},l+2\\}}\|\mathbf{r}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}+1}^{l+2}(\Omega)]^{2}}$
$\displaystyle\leq$ $\displaystyle
Ch^{\min\\{k+\theta^{\prime},\theta^{\prime}+l+2,\theta+\theta^{\prime}+\theta^{\prime}_{1}\\}}\|\mathbf{r}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}+1}^{l+2}(\Omega)]^{2}},$

where $\theta^{\prime}_{1}=\min_{i}\\{\frac{1}{a^{\prime}_{i}}\theta\\}\geq
1$. It can be verified that

$|T_{4}|\leq|T_{41}|+|T_{42}|\leq
Ch^{\min\\{k+\theta^{\prime},\theta+\theta^{\prime}+1,\theta^{\prime}+l+2,2\theta^{\prime}+l+1,\theta+2\theta^{\prime}\\}}\|\mathbf{r}\|_{[\mathcal{K}_{\mathbf{a}^{\prime}+1}^{l+2}(\Omega)]^{2}}.$

By the regularity (4.39) and the summation of estimates $|T_{i}|$,
$i=1,\cdots,4$, and $\theta<2\theta^{\prime}$, we have the error estimate

$\|\mathbf{u}-\mathbf{u}_{n}\|_{[(\mathcal{K}_{\mathbf{b}}^{1}(\Omega))^{*}]^{2}}\leq
Ch^{\min\\{2\theta,\theta+l+1,k+\theta^{\prime},\theta+\theta^{\prime}+1,\theta^{\prime}+l+2\\}}.$
(4.44)

Recall that $k\geq 1$, $\theta\leq k$, and when $\omega>\pi$, we have
$\theta<\theta^{\prime}$, then it follows

$\theta+1\leq k+\theta^{\prime},$ (4.45)

and when $\omega<\pi$, we have $\theta^{\prime}>1$, so the inequality (4.45)
still holds. The estimates in (4.41) follows from (4.44) with the fact (4.45).
∎

######  Remark 4.19.

By Theorem 4.17 and Theorem 4.18, we can find that if we take

$\theta=k$ (4.46)

in the grading parameter $\kappa_{Q_{i}}$, then we can obtain the optimal
convergence rate for the Stokes approximations in Algorithm 3.2,

$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{1}(\Omega)]^{2}}+\|p-p_{n}\|\leq
Ch^{k},$ (4.47a)
$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\leq
Ch^{k+1}.$ (4.47b)

For Algorithm 3.1, it follows

$\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{1}(\Omega)]^{2}}+\|p-p_{n}\|\leq
Ch^{\min\left\\{\frac{\beta_{0}}{\alpha_{0}}\max\\{k,\alpha_{0}\\}+1,k\right\\}},$
(4.48a) $\displaystyle\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\leq
Ch^{\min\left\\{k+1,\frac{\beta_{0}}{\alpha_{0}}\max\\{k,\alpha_{0}\\}+2\right\\}}.$
(4.48b)

However, to obtain the optimal convergence rate for the biharmonic
approximation, the convergence rates of Mini element or Taylor-Hood element
approximations don’t have to be optimal. Therefore, we shall figure out the
admissible parameters $\theta$ such that the convergence rate of the
biharmonic approximation is optimal.

######  Theorem 4.20.

Set the grading parameters $\kappa_{Q_{i}}=2^{-\frac{\theta}{a_{i}}}$ with
$0<a_{i}<\alpha_{0}^{i}$ and

$\theta=\max\\{k-1,a^{\prime}_{i}\\},$ (4.49)

where $a^{\prime}_{i}=\min\\{\alpha_{0},a_{i}\\}\leq\alpha_{0}$ for
$\alpha_{0}$ given by (2.2). Let $\phi_{n}\in V_{n}^{k}$ be the solution of
finite element solution of (3.6), and $\phi$ is the solution of the biharmonic
problem (1.1). If $\phi_{n}$ is the solution in Algorithm 3.2, then it follows

$\|\phi-\phi_{n}\|_{H^{1}(\Omega)}\leq Ch^{k};$ (4.50)

it is the solution in Algorithm 3.1, then it follows

$\|\phi-\phi_{n}\|_{H^{1}(\Omega)}\leq
Ch^{\min\left\\{k,\max\left\\{\frac{\beta_{0}}{\alpha_{0}}(k-1)+2,\beta_{0}+2\right\\}\right\\}}.$
(4.51)

###### Proof.

Denote by $\phi_{I}\in V_{n}^{k}$ the nodal interpolation of $\phi$. Similar
to Theorem 3.7, we have

$\displaystyle\|\phi-\phi_{n}\|_{H^{1}(\Omega)}\leq$ $\displaystyle
C\left(\|\phi-\phi_{I}\|_{H^{1}(\Omega)}+\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\right)$
(4.52)

Recall that
$\phi\in\mathcal{K}_{\mathbf{a}+2}^{m+2}(\Omega)=\mathcal{K}_{(\mathbf{a}+1)+1}^{(m+1)+1}(\Omega)$
with $m\geq k$, so by Lemma 4.12 with grading parameter
$\kappa_{Q_{i}}=2^{-\frac{\theta_{1}}{1+a_{i}}}(=2^{-\frac{\theta}{a_{i}}})$
with $\theta$ given in (4.49) and
$\theta_{1}=\frac{1+a_{i}}{a_{i}}\theta=\theta+\frac{1}{a_{i}}\theta\geq\theta+1\geq
k$, we have

$\|\phi-\phi_{I}\|_{H^{1}(\Omega)}\leq Ch^{\min\\{k,\theta_{1}\\}}=Ch^{k}.$
(4.53)

For $\theta$ given in (4.49), Theorem 4.18 indicates for Algorithm 3.2,

$\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\leq
Ch^{\min\left\\{\max\\{2\alpha_{0},2(k-1)\\},\max\\{\alpha_{0}+1,k\\}\right\\}}.$
(4.54)

Plugging (4.53) and (4.54) into (4.52), the estimate (4.50) holds.

For Algorithm 3.1, we have

$\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\leq
Ch^{\min\left\\{\max\\{2\alpha_{0},2(k-1)\\},\max\\{\alpha_{0}+1,k\\},\frac{\beta_{0}}{\alpha_{0}}\max\\{k-1,\alpha_{0}\\}+2\right\\}}.$
(4.55)

By plugging (4.53) and (4.55) into (4.52), it follows the estimate (4.50). ∎

######  Theorem 4.21.

Set the grading parameters $\kappa_{Q_{i}}=2^{-\frac{\theta}{a_{i}}}$ with
$0<a_{i}<\alpha_{0}^{i}$ and $\theta$ given by

$\theta=\max\left\\{k-1,\frac{k+1}{2},a^{\prime}_{i}\right\\},$ (4.56)

where $a^{\prime}_{i}=\min\\{\alpha_{0},a_{i}\\}\leq\alpha_{0}$ for
$\alpha_{0}$ given by (2.2). Let $\phi_{n}$ be the solution of finite element
solution of (3.6), and $\phi$ be the solution of the biharmonic problem (1.1).
If $\phi_{n}$ is the solution in Algorithm 3.2, then it follows

$\|\phi-\phi_{n}\|\leq Ch^{k+1}.$ (4.57)

it is the solution in Algorithm 3.1, then it follows

$\|\phi-\phi_{n}\|\leq
Ch^{\min\left\\{k+1,\max\left\\{\frac{\beta_{0}}{\alpha_{0}}(k-1)+2,\beta_{0}+2\right\\}+1\right\\}}.$
(4.58)

###### Proof.

Set $\psi=v_{I}\in V_{n}^{k}$ the nodal interpolation of $v$ of the Poisson
problem (3.33). Similar to Theorem 3.8, we have

$\displaystyle\|\phi-\phi_{n}\|^{2}\leq$
$\displaystyle\|\phi-\phi_{n}\|_{H^{1}(\Omega)}\|v-v_{I}\|_{H^{1}(\Omega)}+\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\|v-v_{I}\|_{H^{1}(\Omega)}$
(4.59)
$\displaystyle+\|\mathbf{u}-\mathbf{u}_{n}\|_{[(\mathcal{K}_{\mathbf{b}}^{1}(\Omega))^{*}]^{2}}\|\textbf{curl
}v\|_{[\mathcal{K}_{\mathbf{b}}^{1}(\Omega)]^{2}}:=T_{1}+T_{2}+T_{3}.$

Based on the results in [7], the solution
$v\in\mathcal{K}^{2}_{\mathbf{b}^{\prime}+1}(\Omega)$ satisfies the regularity
estimate

$\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}\leq C\|\phi-\phi_{n}\|,$ (4.60)

where the $i$th entry of $\mathbf{b}^{\prime}$ is given by
$b^{\prime}_{i}=\min\left\\{b_{i},1\right\\}$ with
$b_{i}<\frac{\pi}{\omega_{i}}$. If $\omega>\pi$, we have

$\theta^{\prime}\geq\theta\geq\frac{k+1}{2}\geq 1,$

so it follows the interpolation error

$\|v-v_{I}\|_{H^{1}(\Omega)}\leq
Ch^{\min\\{\theta^{\prime},1\\}}\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}=Ch\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}.$
(4.61)

If $\omega<\pi$, the interpolation error (4.61) is obvious since $v\in
H^{2}(\Omega)$.

For Algorithm 3.2, we have the following estimate for each $T_{i}$, $i=1,2,3$.
By Theorem 4.20 and (4.61), it follows

$T_{1}=\|\phi-\phi_{n}\|_{H^{1}(\Omega)}\|v-v_{I}\|_{H^{1}(\Omega)}\leq
Ch^{k+1}\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}.$

By Theorem 4.18 and (4.61), it follows

$T_{2}=\|\mathbf{u}-\mathbf{u}_{n}\|_{[L^{2}(\Omega)]^{2}}\|v-v_{I}\|_{H^{1}(\Omega)}\leq
Ch^{k+2}\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}.$

Again, by Theorem 4.18, it follows

$T_{3}\leq
C\|\mathbf{u}-\mathbf{u}_{n}\|_{[(\mathcal{K}_{\mathbf{b}^{\prime}}^{1}(\Omega))^{*}]^{2}}\|v\|_{\mathcal{K}_{\mathbf{b}^{\prime}+1}^{2}(\Omega)}\leq
Ch^{k+1}\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}.$

Thus, the regularity estimate (4.60) and the summation of $T_{i}$, $i=1,2,3$
give the estimate (4.57).

For Algorithm 3.1, we have the following estimate for $T_{i}$. By Theorem 4.20
and (4.61), it follows

$T_{1}\leq
Ch^{\min\left\\{k+1,\max\left\\{\frac{\beta_{0}}{\alpha_{0}}(k-1)+2,\beta_{0}+2\right\\}+1\right\\}}\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}.$

By Theorem 4.18 and (4.61), it follows

$T_{2}\leq
Ch^{\min\left\\{k+2,\max\left\\{\frac{\beta_{0}}{\alpha_{0}}(k-1)+2,\beta_{0}+2\right\\}+2\right\\}}\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}.$

Again, it follows by Theorem 4.18,

$T_{3}\leq
Ch^{\min\left\\{k+1,\max\left\\{\frac{\beta_{0}}{\alpha_{0}}(k-1)+2,\beta_{0}+2\right\\}+1\right\\}}\|v\|_{K^{2}_{\mathbf{b}^{\prime}+1}(\Omega)}.$

Thus, the regularity estimate (4.60) and the summation of $T_{i}$, $i=1,2,3$
again give the estimate (4.58).

∎

######  Remark 4.22.

For the results in Theorem 4.20 and Theorem 4.21, we have the following facts.

  * •

If $k=1$ in (4.49), then $\theta=a_{i}$, $i=1,\cdots,N$ gives
$\kappa_{Q_{i}}=\frac{1}{2}$, which indicates the mesh is exactly the uniform
mesh.

  * •

With the grading parameters given, the finite element approximations
$\phi_{n}$ from Algorithm 3.1 achieve the optimal convergence rates at least
for $k=1,2,3$ in both $H^{1}$ and $L^{2}$ norm. Moreover, the optimal
convergence rates can be obtained for any $k\geq 1$ when $\omega>\pi$ by
(2.6).

  * •

With the grading parameters given, the finite element approximations
$\phi_{n}$ from Algorithm 3.2 can achieve optimal convergence rate for any
$k\geq 1$.

  * •

Based on the regularity (4.26) for the involved Stokes problem in Algorithm
3.1, if we take the grading parameter
$\kappa_{Q_{i}}=2^{-\frac{\theta}{c_{i}}}\left(\leq
2^{-\frac{\theta}{a_{i}}}\right)$, then the corresponding error estimates
(4.28), (4.38), (4.41) and (4.51) still hold.

  * •

By (2.40) and $\phi\in H_{0}^{2}(\Omega)$, we have that

$\|\phi\|_{H^{2}(\Omega)}\leq C\|\mathbf{u}\|_{[H^{1}(\Omega)]^{2}},$

so the error $\|\mathbf{u}-\mathbf{u}_{n}\|_{[H^{1}(\Omega)]^{2}}$ is an
estimate of the solution error to $\phi$ in $H^{2}$ norm.

##  5\. Numerical illustrations

In this section, we present numerical tests to validate our theoretical
predictions for the proposed finite element algorithm for solving the
biharmonic problem (1.1). If an exact solution (or vector) $v$ is unknown, we
use the following numerical convergence rate

$\displaystyle{\mathcal{R}}=\log_{2}\frac{|v_{j}-v_{j-1}|_{[H^{l}(\Omega)]^{l^{\prime}}}}{|v_{j+1}-v_{j}|_{[H^{l}(\Omega)]^{l^{\prime}}}},$
(5.1)

$l=0,1$ as an indicator of the actual convergence rate [32]. Here $v_{j}$
denotes the finite element solution on the mesh $\mathcal{T}_{j}$ obtained
after $j$ refinements of the initial triangulation $\mathcal{T}_{0}$. For
scalar functions, we take $l^{\prime}=1$, otherwise, $l^{\prime}=2$. So if
$v_{j}=w_{j},\phi_{j},p_{j}$, we take $l^{\prime}=1$; if
$v_{j}=\mathbf{u}_{j}$, we take $l^{\prime}=2$.

To test the performance of Algorithm 3.1 and Algorithm 3.2 for solving the
biharmonic problem (1.1), we shall use the $H^{2}$-conforming Argyris finite
element approximation [3] as a reference solution $\phi_{R}$, which is
computed on the same mesh as that for Algorithm 3.2. Since the solution of the
$H^{2}$-conforming finite element method converges to the exact solution
$\phi$ regardless of the convexity of the domain as the mesh is refined, so we
can use $\phi_{R}$ as a good approximation of the exact solution $\phi$.

Since the convergence rate of the finite element approximation $w_{j}$ of the
Poisson equation in Algorithm 3.1 has been well investigated in many papers
(see e.g., [35, 34]), so we will not report the convergence rates of $w_{j}$
in the following numerical tests.

######  Example 5.1.

We solve the bibarmonic problem (1.1) with $f=1$ using Algorithm 3.2 based on
polynomials with $k=2$ on uniform meshes obtained by the midpoint refinements.
The source term of involved Stokes problem (2.7) in Algorithm 3.2 is taken as
$\mathbf{F}=(0,x)^{T}$, which satisfies (2.8).

Test case 1. We first consider this problem in a square domain
$\Omega=(-1,1)^{2}$ with the initial mesh given in Figure 6a. The errors in
$L^{\infty}$ norm between the finite element solution $\phi_{j}$ and the
reference solution $\phi_{R}$ are given in Table 2. The finite element
solution and its difference with the reference solution are shown in 6b and
6c, respectively. These results indicate that the solution of Algorithm 3.2
converges to the exact solution. In table 3, we report the $H^{1}$ and $L^{2}$
convergence rates of the finite element solutions $\phi_{j}$, respectively.
The results indicate that optimal convergence rates are obtained for finite
element solutions of the biharmonic problem.

In addition, the Taylor-Hood element approximations $\mathbf{u}_{7}$ and
$p_{7}$ for the involved Stokes problem are shown in Figure 6d-6f. In Table 3,
we also report the $H^{1}$ and/or $L^{2}$ convergence rates of the Taylor-Hood
element approximations $\mathbf{u}_{j}$ and $p_{j}$, the results imply that
optimal convergence rates are obtained for Taylor-Hood element approximations
of the Stokes problem.

These results are consistent with our expectation in Lemma 3.4, Theorem 3.7,
and Theorem 3.8 for the biharmonic problem (1.1) and the involved Stokes
problem (2.7) in a convex domain.

Table 2. The $L^{\infty}$ error $\|\phi_{R}-\phi_{j}\|_{L^{\infty}(\Omega)}$ in the square domain on quasi-uniform meshes. | $j=3$ | $j=4$ | $j=5$ | $j=6$ | $j=7$ | $j=8$
---|---|---|---|---|---|---
$k=2$ | 6.66250e-05 | 9.10404e-06 | 1.17798e-06 | 1.49996e-07 | 1.89341e-08 | 2.37885e-09

Figure 6. The Square domain (Example 5.1): (a) the initial mesh; (b) $P_{2}$ finite element approximation $\phi_{7}$; (c) the difference $|\phi_{R}-\phi_{7}|$; (d) Taylor-Hood element approximation $u_{1}$ of $\mathbf{u}_{7}$; (e) Taylor-Hood element approximation $u_{2}$ of $\mathbf{u}_{7}$; (d) Taylor-Hood element approximation $p_{7}$.  Table 3. Numerical convergence rates in the square domain on uniform meshes. |  |  $H^{1}$ rate of $\phi_{j}$ |  $L^{2}$ rate of $\phi_{j}$ |  $H^{1}$ rate of $\mathbf{u}_{j}$ |  $L^{2}$ rate of $\mathbf{u}_{j}$ | $L^{2}$ rate of $p_{j}$
---|---|---|---|---|---|---
$k=2$ | $j=4$ | 1.96 | 3.00 | 1.99 | 3.01 | 2.03
$j=5$ | 1.99 | 3.00 | 2.00 | 3.02 | 2.02
$j=6$ | 2.00 | 3.00 | 2.00 | 3.01 | 2.01
$j=7$ | 2.00 | 3.00 | 2.00 | 3.00 | 2.00
$j=8$ | 2.00 | 3.00 | 2.00 | 3.00 | 2.00

Test case 2. We then consider this problem in an L-shaped domain
$\Omega=\Omega_{0}\setminus\Omega_{1}$ with $\Omega_{0}=(-1,1)^{2}$ and
$\Omega_{1}=(0,1)\times(-1,0)$ based on the initial mesh given in Figure 7a.
The error $\|\phi_{R}-\phi_{j}\|_{L^{\infty}(\Omega)}$ is given in Table 4.
The finite element solution and its difference with the reference solution are
shown in Figure 7b and 7c, respectively. These results indicate that the
solutions of Algorithm 3.2 converge to the exact solution in a nonconvex
polygonal domain.

Table 4. The $L^{\infty}$ error $\|\phi_{R}-\phi_{j}\|_{L^{\infty}(\Omega)}$ in the L-shaped domain on quasi-uniform meshes. | $j=3$ | $j=4$ | $j=5$ | $j=6$ | $j=7$ | $j=8$
---|---|---|---|---|---|---
$k=2$ | 8.74987e-04 | 3.94122e-04 | 1.77980e-04 | 8.26205e-05 | 3.86434e-05 | 1.81330e-05

Figure 7. The L-shaped domain (Example 5.1): (a) the initial mesh; (b) $P_{2}$
finite element approximation $\phi_{7}$; (c) the difference
$|\phi_{R}-\phi_{7}|$.

######  Example 5.2.

We test Example 5.1 Test case 2 again using Algorithm 3.1 and Algorithm 3.2.
For Algorithm 3.2, we will consider a different source term
${\mathbf{F}}=(-y,0)^{T}$ for the involved Stokes problem, and the new source
term also satisfies (2.8), namely, $\text{curl }\mathbf{F}=f=1$. Here we
denote the finite element approximations in Example 5.1 Test case 2 (namely,
Algorithm 3.2 with $\mathbf{F}=(0,x)^{T}$) by $\tilde{\phi}_{j}$ for the
biharmonic problem and by $\tilde{\mathbf{u}}_{j}$ and $\tilde{p}_{j}$ for the
involved Stokes problem.

We compare the finite element approximations with those obtained in Example
5.1 Test case 2 on uniform meshes. The errors of the finite element
approximations of the biharmonic problem are shown in Table 5, and these for
the Taylor-Hood element approximations of the involved Stokes problem are
reported in Table 6. From these results, we find that the errors of finite
element approximations ${\phi}_{j}$ and Taylor-Hood element approximations
${\mathbf{u}}_{j}$ of both Algorithm 3.1 and Algorithm 3.2 are converging to
the solutions of solutions in Example 5.1 Test case 2. Therefore, the finite
element approximations $\phi_{j}$ from both Algorithm 3.1 and Algorithm 3.2
converge to the exact solution of the biharmonic problem (1.1). All these
results also indicate that the finite element approximations $\phi_{j}$ and
the Taylor-Hood element approximations $\mathbf{u}_{j}$ are independent of the
choice of the source term $\mathbf{F}$ as long as it satisfies (2.8), and
these approximations are uniquely determined by $f$ in (1.1). The results are
consistent with the theoretical results in Lemma 2.4 and Corollary 2.5. We
also find that the errors between $\tilde{p}_{j}$ and $p_{j}$ are not
converging, which implies that $p_{j}$ depends on the specific $\mathbf{F}$ as
indicated in Lemma 2.4.

Table 5. Errors with the reference solution and convergence rates in the L-shaped domain. | Algorithm 3.1 | Algorithm 3.2
---|---|---
| $\|\phi_{j}-\tilde{\phi}_{j}\|_{H^{1}(\Omega)}$ | $\|\phi_{j}-\tilde{\phi}_{j}\|$ | $\|\phi_{j}-\tilde{\phi}_{j}\|_{H^{1}(\Omega)}$ | $\|\phi_{j}-\tilde{\phi}_{j}\|$
$j=5$ | 1.73930e-06 | 3.83797e-07 | 7.97013e-09 | 3.72412e-10
$j=6$ | 3.73791e-07 | 8.38856e-08 | 5.88353e-10 | 2.10894e-11
$j=7$ | 8.05609e-08 | 1.82423e-08 | 4.56271e-11 | 1.24448e-12
$j=8$ | 1.73829e-08 | 3.95457e-09 | 3.70076e-12 | 7.54393e-14
$j=9$ | 3.75263e-09 | 8.55733e-10 | 3.10543e-13 | 4.80570e-15
$j=10$ | 8.10299e-10 | 1.84999e-10 | 2.81240e-14 | $--$
Table 6. Errors with the reference solution and convergence rates in the L-shaped domain. | Algorithm 3.1 | Algorithm 3.2
---|---|---
| $\|\mathbf{u}_{j}-\tilde{\mathbf{u}}_{j}\|_{[H^{1}(\Omega)]^{2}}$ | $\|\mathbf{u}_{j}-\tilde{\mathbf{u}}_{j}\|_{[L^{2}(\Omega)]^{2}}$ | $\|p_{j}-\tilde{p}_{j}\|$ | $\|\mathbf{u}_{j}-\tilde{\mathbf{u}}_{j}\|_{[H^{1}(\Omega)]^{2}}$ | $\|\mathbf{u}_{j}-\tilde{\mathbf{u}}_{j}\|_{[L^{2}(\Omega)]^{2}}$ | $\|p_{j}-\tilde{p}_{j}\|$
$j=5$ | 7.99655e-05 | 1.36408e-06 | 4.19809e-01 | 1.10073e-05 | 7.49412e-08 | 5.59017e-01
$j=6$ | 2.48630e-05 | 2.79879e-07 | 4.19758e-01 | 1.95571e-06 | 6.67031e-09 | 5.59017e-01
$j=7$ | 7.79196e-06 | 5.87681e-08 | 4.19736e-01 | 3.46595e-07 | 5.91602e-10 | 5.59017e-01
$j=8$ | 2.44954e-06 | 1.24998e-08 | 4.19727e-01 | 6.13467e-08 | 5.23793e-11 | 5.59017e-01
$j=9$ | 7.70977e-07 | 2.67752e-09 | 4.19723e-01 | 1.08514e-08 | 4.63360e-12 | 5.59017e-01
$j=10$ | 2.42772e-07 | 5.75794e-10 | 4.19722e-01 | 1.91888e-09 | 4.09765e-13 | 5.59017e-01

######  Example 5.3.

We solve the problem in Example 5.1 Test Case 2 using both Algorithm 3.1 and
Algorithm 3.2 with polynomials $k=1,2$ on a sequence of graded meshes
(including uniform mesh). The initial mesh and the graded mesh after 2 mesh
refinements are shown in Figure 7a and Figure 8a, respectively.

Figure 8. The L-shaped domain (Example 5.3): (a) the graded mesh after two
mesh refinements with $\kappa=0.2$; (b) Taylor-Hood element approximation
$u_{1}$ of $\mathbf{u}_{7}$; (e) Taylor-Hood element approximation $u_{2}$ of
$\mathbf{u}_{7}$; (d) Taylor-Hood element approximation $p_{7}$.

In Table 7 (Algorithm 3.2) and Table 8 (Algorithm 3.1), we show the numerical
convergence rates of finite element approximation $\phi_{j}$ in both $H^{1}$
and $L^{2}$ norm to the solution of the biharmonic problem. We note that the
convergence rates from both Algorithm 3.2 and Algorithm 3.1 are almost the
same. From Table Table 7 and Table 8, we find on uniform meshes ($\kappa=0.5$)
that the $H^{1}$ convergence rate of the $P_{1}$ finite element approximation
is optimal with $\mathcal{R}=1$, and that of the $P_{2}$ finite element
approximations is suboptimal with $\mathcal{R}\approx 1.10$. Both of them are
consistent the theoretical result in Theorem 3.7 in an L-shaped domain, that
is
$\mathcal{R}_{\text{exact}}=\min\\{k,\alpha_{0}+1,2\alpha_{0}\\}\approx\min\\{k,1.09\\}$
for Algorithm 3.2, and
$\mathcal{R}_{\text{exact}}=\min\\{k,\beta_{0}+2,\alpha_{0}+1,2\alpha_{0}\\}\approx\min\\{k,1.09\\}$
for Algorithm 3.1, where $\alpha_{0}$ is given in Table 1 with
$\omega=\frac{3\pi}{2}$, and $\beta_{0}=\frac{\pi}{\omega}$. We also find that
the convergence rates of $P_{1}$ finite element approximations are optimal
with $\mathcal{R}=1$ on graded meshes with $\kappa<0.5$, and that of the
$P_{2}$ finite element approximations are optimal with $\mathcal{R}=2$ on
graded meshes with $\kappa\leq 0.3$, which are close to the theoretical result
in Theorem 4.20, namely, the optimal convergence rate can be obtained when
$\kappa<2^{-\frac{\alpha_{0}}{\alpha_{0}}}=0.5$ for $P_{1}$ finite element
approximations, and $\kappa<2^{-\frac{1}{\alpha_{0}}}\approx 0.28$ for $P_{2}$
finite element approximations.

Again from Table Table 7 and Table 8, the $L^{2}$ convergence rates of both
$P_{1}$ and $P_{2}$ finite element approximations on uniform meshes are
suboptimal with $\mathcal{R}=1.14$ and $\mathcal{R}\approx 1.09$, which are
consistent with the theoretical result
$\mathcal{R}_{\text{exact}}\approx\min\\{k+1,1.09\\}$ for both Algorithm 3.2
and Algorithm 3.1 in Theorem 3.8. On graded meshes, the convergence rates of
$P_{1}$ finite element approximations are optimal with $\kappa\leq 0.2$, and
those of the $P_{2}$ finite element approximations are optimal with
$\kappa\leq 0.1$, which are consistent with the theoretical result in Theorem
4.21, namely, the optimal convergence rate $\mathcal{R}=2$ can be obtained
when $\kappa<2^{-\frac{1}{\alpha_{0}}}\approx 0.28$ for $P_{1}$ finite element
approximations, and $\mathcal{R}=3$ can be achieved when
$\kappa<2^{-\frac{1.5}{\alpha_{0}}}\approx 0.15$ for $P_{2}$ finite element
approximations.

Table 7. Convergence history of finite element approximation $\phi_{j}$ of the biharmonic problem from Algorithm 3.2 in the L-shaped domain. |  |  $H_{1}$ rate of $\phi_{j}$ | $L_{2}$ rate of $\phi_{j}$
---|---|---|---
| $\kappa$ | 0.05 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.05 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5
$k=1$ | $j=6$ | $--$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | $--$ | 2.00 | 1.99 | 1.96 | 1.84 | 1.52
$j=7$ | $--$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | $--$ | 2.00 | 2.00 | 1.96 | 1.79 | 1.37
$j=8$ | $--$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | $--$ | 2.00 | 2.00 | 1.96 | 1.72 | 1.26
$j=9$ | $--$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | $--$ | 2.00 | 2.00 | 1.95 | 1.66 | 1.18
$j=10$ | $--$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | $--$ | 2.00 | 2.00 | 1.95 | 1.60 | 1.14
$k=2$ | $j=6$ | 1.98 | 1.99 | 1.99 | 1.99 | 1.89 | 1.37 | 3.05 | 3.02 | 2.99 | 1.94 | 1.43 | 1.08
$j=7$ | 1.99 | 2.00 | 2.00 | 2.00 | 1.81 | 1.23 | 3.03 | 3.01 | 2.98 | 1.90 | 1.43 | 1.08
$j=8$ | 2.00 | 2.00 | 2.00 | 2.00 | 1.71 | 1.15 | 3.02 | 3.00 | 2.96 | 1.89 | 1.44 | 1.08
$j=9$ | 2.00 | 2.00 | 2.00 | 2.00 | 1.61 | 1.12 | 3.01 | 3.00 | 2.92 | 1.89 | 1.44 | 1.08
Table 8. Convergence history of finite element approximation $\phi_{j}$ of the biharmonic problem from Algorithm 3.1 in the L-shaped domain. |  |  $H_{1}$ rate of $\phi_{j}$ | $L_{2}$ rate of $\phi_{j}$
---|---|---|---
| $\kappa$ | 0.05 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.05 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5
$k=1$ | $j=6$ | $--$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | $--$ | 1.99 | 1.99 | 1.96 | 1.85 | 1.54
$j=7$ | $--$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | $--$ | 2.00 | 2.00 | 1.96 | 1.80 | 1.38
$j=8$ | $--$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | $--$ | 2.00 | 2.00 | 1.96 | 1.73 | 1.26
$j=9$ | $--$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | $--$ | 2.00 | 2.00 | 1.96 | 1.67 | 1.19
$j=10$ | $--$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | $--$ | 2.00 | 2.00 | 1.95 | 1.61 | 1.14
$k=2$ | $j=6$ | 1.98 | 1.99 | 1.99 | 1.99 | 1.89 | 1.35 | 3.04 | 3.02 | 2.99 | 1.94 | 1.41 | 1.04
$j=7$ | 1.99 | 2.00 | 2.00 | 2.00 | 1.81 | 1.22 | 3.03 | 3.01 | 2.98 | 1.90 | 1.43 | 1.06
$j=8$ | 2.00 | 2.00 | 2.00 | 2.00 | 1.71 | 1.15 | 3.01 | 3.00 | 2.96 | 1.89 | 1.43 | 1.07
$j=9$ | 2.00 | 2.00 | 2.00 | 2.00 | 1.61 | 1.12 | 3.01 | 3.00 | 2.92 | 1.89 | 1.44 | 1.08

The Taylor-Hood element approximations $\mathbf{u}_{7}$ and $p_{7}$ based on
Algorithm 3.2 on uniform meshes for the involved Stokes problem are shown in
Figure 8b-8d. In Table 9-12, we display numerical convergence rates of the
Mini element approximations and the Taylor-Hood approximations from both
Algorithm 3.2 and Algorithm 3.1 for the involved Stokes problem. The $H^{1}$
convergence rates of $\mathbf{u}_{j}$ and the $L^{2}$ convergence rates of
$p_{j}$ with $k=1,2$ are suboptimal on uniform meshes with convergence rates
$\mathcal{R}\approx 0.54$, which are consistent with the theoretical result
$\mathcal{R}_{\text{exact}}=\alpha_{0}\approx 0.54$ in Lemma 3.4 and Lemma 3.6
in an L-shaped domain. On graded meshes, the convergence rates are optimal
with $\mathcal{R}=1$ when $\kappa\leq 0.2$ for $k=1$, and $\mathcal{R}=2$ when
$\kappa\leq 0.05$ for $k=2$, these are consistent with the results in Theorem
4.17 and Remark 4.19, namely, the optimal convergence rate can be achieved
when $\kappa<2^{-\frac{1}{\alpha_{0}}}\approx 0.28$ for $k=1$, and
$\kappa<2^{-\frac{2}{\alpha_{0}}}\approx 0.08$ for $k=2$.

The $L^{2}$ convergence rates of $\mathbf{u}_{j}$ are suboptimal on uniform
meshes with convergence rates $\mathcal{R}\approx 1.13$ and
$\mathcal{R}\approx 1.12$, which are consistent with the theoretical result
$\mathcal{R}=2\alpha_{0}\approx 1.09$ in Lemma 3.4 and Lemma 3.6 in an
L-shaped domain. On graded meshes, the convergence rates are optimal with
$\mathcal{R}=2$ for $\kappa\leq 0.2$, and with $\mathcal{R}=3$ for $\kappa\leq
0.1$, which are close to the theoretical results in Theorem 4.17 and Remark
4.19, namely, the optimal convergence rate can be achieved when
$\kappa<2^{-\frac{1}{\alpha_{0}}}\approx 0.28$ for $k=1$, and
$\kappa<2^{-\frac{2}{\alpha_{0}}}\approx 0.08$ for $k=2$.

Table 9. Convergence history of the Mini element approximations ($k=1$) of Stokes problem from Algorithm 3.2 in the L-shaped domain. |  $H_{1}$ rate of $\mathbf{u}_{j}$ |  $L^{2}$ rate of $\mathbf{u}_{j}$ | $L_{2}$ rate of $p_{j}$
---|---|---|---
$\kappa$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5
$j=6$ | 1.00 | 1.00 | 0.98 | 0.88 | 0.69 | 2.01 | 2.01 | 1.97 | 1.80 | 1.40 | 1.23 | 1.35 | 1.15 | 0.77 | 0.57
$j=7$ | 1.00 | 1.00 | 0.98 | 0.85 | 0.64 | 2.01 | 2.01 | 1.96 | 1.71 | 1.28 | 1.10 | 1.35 | 1.07 | 0.74 | 0.56
$j=8$ | 1.00 | 1.00 | 0.97 | 0.82 | 0.60 | 2.01 | 2.00 | 1.95 | 1.62 | 1.20 | 1.05 | 1.35 | 1.01 | 0.73 | 0.55
$j=9$ | 1.00 | 1.00 | 0.97 | 0.80 | 0.58 | 2.01 | 2.00 | 1.94 | 1.55 | 1.16 | 1.15 | 1.34 | 0.98 | 0.72 | 0.55
$j=10$ | 1.00 | 1.00 | 0.97 | 0.78 | 0.56 | 2.00 | 2.00 | 1.94 | 1.51 | 1.13 | 1.29 | 1.33 | 0.96 | 0.72 | 0.55
Table 10. Convergence history of the Taylor-Hood element approximations ($k=2$) of Stokes problem from Algorithm 3.2 in the L-shaped domain. |  |  $H_{1}$ rate | $L_{2}$ rate
---|---|---|---
| $\kappa$ | 0.05 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.05 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5
$\mathbf{u}_{j}$ | $j=6$ | 1.82 | 1.83 | 1.37 | 0.96 | 0.72 | 0.55 | 2.97 | 3.00 | 2.99 | 2.04 | 1.57 | 1.32
$j=7$ | 1.90 | 1.84 | 1.31 | 0.95 | 0.72 | 0.54 | 3.03 | 3.01 | 2.99 | 1.93 | 1.50 | 1.24
$j=8$ | 1.95 | 1.83 | 1.28 | 0.95 | 0.72 | 0.54 | 3.03 | 3.01 | 2.98 | 1.90 | 1.47 | 1.19
$j=9$ | 1.97 | 1.83 | 1.27 | 0.95 | 0.72 | 0.54 | 3.01 | 3.00 | 2.95 | 1.89 | 1.45 | 1.15
$j=10$ | 1.98 | 1.82 | 1.27 | 0.95 | 0.72 | 0.54 | 3.01 | 3.00 | 2.91 | 1.89 | 1.45 | 1.12
$p_{j}$ | $j=6$ |  | 1.81 | 1.82 | 1.31 | 0.95 | 0.72 | 0.55
$j=7$ |  | 1.91 | 1.82 | 1.28 | 0.95 | 0.72 | 0.55
$j=8$ |  | 1.97 | 1.82 | 1.27 | 0.95 | 0.72 | 0.54
$j=9$ |  | 1.99 | 1.82 | 1.27 | 0.95 | 0.72 | 0.54
$j=10$ |  | 2.00 | 1.81 | 1.27 | 0.95 | 0.72 | 0.54
Table 11. Convergence history of the Mini element approximations ($k=1$) of Stokes problem from Algorithm 3.1 in the L-shaped domain. |  $H_{1}$ rate of $\mathbf{u}_{j}$ |  $L^{2}$ rate of $\mathbf{u}_{j}$ | $L_{2}$ rate of $p_{j}$
---|---|---|---
$\kappa$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5
$j=6$ | 1.00 | 1.00 | 0.98 | 0.88 | 0.69 | 2.01 | 2.00 | 1.97 | 1.81 | 1.41 | 1.22 | 1.35 | 1.15 | 0.76 | 0.54
$j=7$ | 1.00 | 1.00 | 0.98 | 0.85 | 0.64 | 2.01 | 2.00 | 1.97 | 1.73 | 1.29 | 1.10 | 1.35 | 1.07 | 0.74 | 0.54
$j=8$ | 1.00 | 1.00 | 0.97 | 0.82 | 0.60 | 2.01 | 2.00 | 1.96 | 1.64 | 1.21 | 1.04 | 1.35 | 1.01 | 0.73 | 0.54
$j=9$ | 1.00 | 1.00 | 0.97 | 0.80 | 0.58 | 2.00 | 2.00 | 1.95 | 1.57 | 1.16 | 1.15 | 1.34 | 0.98 | 0.72 | 0.54
$j=10$ | 1.00 | 1.00 | 0.97 | 0.78 | 0.56 | 2.00 | 2.00 | 1.94 | 1.52 | 1.13 | 1.29 | 1.35 | 0.96 | 0.72 | 0.54
Table 12. Convergence history of the Taylor-Hood element approximations ($k=2$) of Stokes problem from Algorithm 3.1 in the L-shaped domain. |  |  $H_{1}$ rate | $L_{2}$ rate
---|---|---|---
| $\kappa$ | 0.05 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.05 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5
$\mathbf{u}_{j}$ | $j=6$ | 1.82 | 1.83 | 1.37 | 0.96 | 0.72 | 0.53 | 2.97 | 2.99 | 2.98 | 2.03 | 1.56 | 1.29
$j=7$ | 1.90 | 1.84 | 1.31 | 0.95 | 0.72 | 0.54 | 3.03 | 3.01 | 2.99 | 1.93 | 1.50 | 1.23
$j=8$ | 1.95 | 1.83 | 1.28 | 0.95 | 0.72 | 0.54 | 3.03 | 3.01 | 2.97 | 1.90 | 1.47 | 1.18
$j=9$ | 1.97 | 1.83 | 1.27 | 0.95 | 0.72 | 0.54 | 3.01 | 3.00 | 2.95 | 1.89 | 1.45 | 1.14
$j=10$ | 1.98 | 1.82 | 1.27 | 0.95 | 0.72 | 0.54 | 3.00 | 2.98 | 2.93 | 1.89 | 1.45 | 1.12
$p_{j}$ | $j=6$ |  | 1.75 | 1.71 | 1.27 | 0.94 | 0.71 | 0.53
$j=7$ |  | 1.89 | 1.75 | 1.27 | 0.95 | 0.72 | 0.53
$j=8$ |  | 1.96 | 1.77 | 1.27 | 0.95 | 0.72 | 0.54
$j=9$ |  | 1.99 | 1.78 | 1.26 | 0.95 | 0.72 | 0.54
$j=10$ |  | 2.00 | 1.79 | 1.26 | 0.95 | 0.72 | 0.54

######  Example 5.4.

We consider the biharmonic problem (1.1) with $f=1$ in the polygonal domain
$\Omega$ (see Figure 9a) enclosed by line segments $Q_{1}Q_{2}$, $Q_{2}Q_{3}$,
$Q_{3}Q_{4}$, and $Q_{4}Q_{1}$ with $Q_{1}(0,0)$,
$Q_{2}(\frac{2}{c_{1}/c_{2}+1},-\frac{2c_{1}}{c_{1}/c_{2}+1})$, $Q_{3}(2,0)$,
and $Q_{4}(\frac{2}{c_{1}/c_{2}+1},\frac{2c_{1}}{c_{1}/c_{2}+1})$, here
$c_{1}=\tan\left(\frac{11\pi}{24}\right)$, and
$c_{2}=\tan\left(\frac{13\pi}{72}\right)$. The largest interior angle
$\omega=\angle Q_{2}Q_{1}Q_{4}=\frac{11\pi}{12}$. We solve this problem by
Algorithm 3.1 and Algorithm 3.2 on both graded meshes and uniform meshes based
on polynomials with $k=2$. We take $\mathbf{F}=(0,x)^{T}$ as the source term
of involved Stokes problem (2.7) in Algorithm 3.2. The finite element
approximations $w_{7}$ and $\phi_{7}$ from Algorithm 3.1 are shown in Figure
9b and Figure 9c, respectively. We show both $H^{1}$ and $L^{2}$ convergence
rates of the finite element approximation $\phi_{j}$ from both Algorithm 3.1
and Algorithm 3.2 in Table 13. We find the the convergence rates are almost
the same for these two algorithms.

For $H^{1}$ convergence rate of $\phi_{j}$ from both algorithms, we find that
the convergence rates on uniform meshes are $\mathcal{R}=2.00$, which is
consistent with the expected convergence rate
$\mathcal{R}_{\text{exact}}=\min\\{k,\alpha_{0}+1,2\alpha_{0}\\}=2$ in Theorem
3.7, where $\alpha_{0}\approx 1.20$ corresponding to the interior angle
$\omega=\frac{11\pi}{12}$ shown in Table 1. On graded meshes, the optimal
convergence rate $\mathcal{R}=2.0$ can be observed on graded meshes for
$\kappa<0.5$, which is consistent with the theoretical result
$\kappa<2^{-\frac{\alpha_{0}}{\alpha_{0}}}=0.5$ in Theorem 4.20. For $L^{2}$
convergence rate, we find that on uniform meshes the convergence rates
$\mathcal{R}\approx 2.98$, which is a little bit larger than the expected rate
$\mathcal{R}_{\text{exact}}\approx 2\alpha_{0}=2.40$. On graded meshes, we
find that the optimal convergence rates $\mathcal{R}=3$ can be obtained when
$\kappa\leq 0.4$, which is consistent with the theoretical result
$\kappa<2^{-\frac{1.5}{\alpha_{0}}}\approx 0.42$ in Theorem 4.21.

Figure 9. A convex domain with $\omega=\frac{11\pi}{12}$ (Example 5.4): (a) the initial mesh; (b) $P_{2}$ finite element approximation $w_{7}$; (c) $P_{2}$ finite element approximation $\phi_{7}$; (d) Taylor-Hood element approximation $u_{1}$ of $\mathbf{u}_{7}$; (e) Taylor-Hood element approximation $u_{2}$ of $\mathbf{u}_{7}$; (f) Taylor-Hood element approximation $p_{7}$. Table 13. Convergence history of finite element approximation of the biharmonic problem with $\omega=\frac{11\pi}{12}$. |  |  $H_{1}$ rate of $\phi_{j}$ | $L_{2}$ rate of $\phi_{j}$
---|---|---|---
| $\kappa$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5
Algorithm 3.2 | $j=6$ | 1.99 | 1.99 | 1.99 | 2.00 | 1.99 | 3.02 | 3.01 | 3.00 | 3.00 | 2.994
$j=7$ | 2.00 | 2.00 | 2.00 | 2.00 | 1.99 | 3.00 | 3.00 | 3.00 | 3.00 | 3.995
$j=8$ | 2.00 | 2.00 | 2.00 | 2.00 | 2.00 | 3.00 | 3.00 | 3.00 | 3.00 | 2.994
$j=9$ | 2.00 | 2.00 | 2.00 | 2.00 | 2.00 | 3.00 | 3.00 | 3.00 | 3.00 | 2.990
$j=10$ | 2.00 | 2.00 | 2.00 | 2.00 | 2.00 | 3.00 | 3.00 | 3.00 | 3.00 | 2.992
Algorithm 3.1 | $j=6$ | 1.99 | 1.99 | 1.99 | 2.00 | 1.99 | 3.02 | 3.01 | 3.00 | 3.00 | 2.994
$j=7$ | 2.00 | 2.00 | 2.00 | 2.00 | 1.99 | 3.01 | 3.00 | 3.00 | 3.00 | 2.995
$j=8$ | 2.00 | 2.00 | 2.00 | 2.00 | 2.00 | 3.00 | 3.00 | 3.00 | 3.00 | 2.994
$j=9$ | 2.00 | 2.00 | 2.00 | 2.00 | 2.00 | 3.00 | 3.00 | 3.00 | 3.00 | 2.990
$j=10$ | 2.00 | 2.00 | 2.00 | 2.00 | 2.00 | 3.00 | 3.00 | 3.00 | 3.00 | 2.984

We show the Taylor-Hood element approximations $\mathbf{u}_{7}$ and $p_{7}$
from Algorithm 3.1 on uniform meshes in Figure 9d-9f. We also report the
convergence rate of the Taylor-Hood element approximations of the involved
Stokes problem in Table 14 (Algorithm 3.2) and Table 15 (Algorithm 3.1). The
$H^{1}$ convergence rate of $\mathbf{u}_{j}$ and the $L^{2}$ convergence rate
of $p_{j}$ are converging with the rate $\mathcal{R}=1.20$ on uniform meshes,
and it is consistent with the theoretical result
$\mathcal{R}_{\text{exact}}=\alpha_{0}\approx 1.20$ in Lemma 3.4 and Lemma
3.6. The optimal convergence rate $\mathcal{R}=2$ can be observed for
$\kappa\leq 0.3$, which is consistent with the theoretical result
$\kappa<2^{-\frac{2}{\alpha_{0}}}\approx 0.31$ in Theorem 4.17 and Remark
4.19. The $L^{2}$ convergence rate of $\mathbf{u}_{j}$ on uniform meshes is
$\mathcal{R}=2.20$, which is consistent with the expected rate
$\mathcal{R}_{\text{exact}}=\alpha_{0}+1\approx 2.20$. On graded meshes, the
optimal convergence rate $\mathcal{R}=3$ is observed with $\kappa\leq 0.3$,
which is also consistent with the theoretical result
$\kappa<2^{-\frac{2}{\alpha_{0}}}\approx 0.31$ in Theorem 4.18 and Remark
4.19.

Table 14. Convergence history of the Taylor-Hood element approximations ($k=2$) of Stokes problem from Algorithm 3.2 with $\omega=\frac{11\pi}{12}$. |  $H_{1}$ rate of $\mathbf{u}_{j}$ |  $L^{2}$ rate of $\mathbf{u}_{j}$ | $L_{2}$ rate of $p_{j}$
---|---|---|---
$\kappa$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5
$j=6$ | 2.00 | 2.00 | 2.00 | 1.88 | 1.32 | 3.01 | 3.00 | 3.01 | 2.99 | 2.26 | 2.01 | 2.00 | 2.00 | 1.94 | 1.42
$j=7$ | 2.00 | 2.00 | 2.00 | 1.83 | 1.24 | 3.00 | 3.00 | 3.00 | 2.98 | 2.22 | 2.00 | 2.00 | 2.00 | 1.90 | 1.29
$j=8$ | 2.00 | 2.00 | 2.00 | 1.77 | 1.21 | 3.00 | 3.00 | 3.00 | 2.98 | 2.21 | 2.00 | 2.00 | 2.00 | 1.85 | 1.23
$j=9$ | 2.00 | 2.00 | 2.00 | 1.71 | 1.20 | 3.00 | 3.00 | 3.00 | 2.98 | 2.20 | 2.00 | 2.00 | 2.00 | 1.79 | 1.21
$j=10$ | 2.00 | 2.00 | 2.00 | 1.67 | 1.20 | 3.00 | 3.00 | 3.00 | 2.97 | 2.20 | 2.00 | 2.00 | 2.00 | 1.73 | 1.20
Table 15. Convergence history of the Taylor-Hood element approximations ($k=2$) of Stokes problem from Algorithm 3.1 with $\omega=\frac{11\pi}{12}$. |  $H_{1}$ rate of $\mathbf{u}_{j}$ |  $L^{2}$ rate of $\mathbf{u}_{j}$ | $L_{2}$ rate of $p_{j}$
---|---|---|---
$\kappa$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5
$j=6$ | 2.00 | 2.00 | 1.99 | 1.88 | 1.31 | 3.01 | 3.00 | 3.00 | 2.99 | 2.26 | 2.04 | 2.02 | 2.00 | 1.75 | 1.20
$j=7$ | 2.00 | 2.00 | 2.00 | 1.83 | 1.24 | 3.00 | 3.00 | 3.00 | 2.98 | 2.22 | 2.02 | 2.01 | 2.00 | 1.69 | 1.19
$j=8$ | 2.00 | 2.00 | 2.00 | 1.77 | 1.21 | 3.00 | 3.00 | 3.00 | 2.98 | 2.20 | 2.01 | 2.00 | 1.99 | 1.65 | 1.19
$j=9$ | 2.00 | 2.00 | 2.00 | 1.71 | 1.20 | 3.00 | 3.00 | 3.00 | 2.98 | 2.20 | 2.00 | 2.00 | 1.99 | 1.63 | 1.19
$j=10$ | 2.00 | 2.00 | 2.00 | 1.67 | 1.20 | 3.00 | 3.00 | 3.00 | 2.97 | 2.20 | 2.00 | 2.00 | 2.00 | 1.61 | 1.20

######  Example 5.5.

In this example, we compare the CPU time and the memory usage of the proposed
finite element algorithms (Algorithm 3.2 and Algorithm 3.1) with those of the
$H^{2}$-conforming Argyris finite element method by solving the biharmonic
problem (1.1) in Example 5.1 on the same meshes. The results of the CPU time
comparison (in seconds) are shown in Table 16. The results of the memory usage
comparison (in GB) are shown in Table 17. All results are tested on MATLAB
R2021a in Linux with 16 GB memory and Intel® Core${}^{\text{TM}}$ i7-6600U
processors.

Table 16. The CPU time (in seconds) of Algorithm 3.2, Aglorithm 3.1, and the Argyris finite element method. (“$--$” represents running out of memory) method$\backslash j$ | Example 5.1 Test case 1 | Example 5.1 Test case 2
---|---|---
$j=4$ | $j=5$ | $j=6$ | $j=7$ | $j=8$ | $j=4$ | $j=5$ | $j=6$ | $j=7$ | $j=8$
Argyris FEM | 7.09 | 28.57 | 118.10 | 492.43 | $--$ | 5.27 | 21.62 | 87.70 | 367.76 | $--$
Algorithm 3.2(k=1) | 0.08 | 0.39 | 2.82 | 15.92 | 90.94 | 0.05 | 0.24 | 1.13 | 9.59 | 58.84
Algorithm 3.1(k=1) | 0.10 | 0.42 | 3.01 | 16.09 | 96.83 | 0.08 | 0.30 | 1.36 | 11.37 | 66.22
Algorithm 3.2(k=2) | 0.88 | 2.29 | 9.05 | 37.72 | 150.40 | 0.60 | 1.90 | 7.03 | 28.21 | 120.95
Algorithm 3.1(k=2) | 0.91 | 2.74 | 10.79 | 43.29 | 181.42 | 0.77 | 2.27 | 8.25 | 34.25 | 147.97
Table 17. The memory usage (in GB) of Algorithm 3.2, Aglorithm 3.1, and the Argyris finite element method. (“$--$” represents running out of memory) method$\backslash j$ | Example 5.1 Test case 1 | Example 5.1 Test case 2
---|---|---
$j=4$ | $j=5$ | $j=6$ | $j=7$ | $j=8$ | $j=4$ | $j=5$ | $j=6$ | $j=7$ | $j=8$
Argyris FEM | 0.06 | 0.27 | 1.00 | 5.10 | $--$ | 0.04 | 0.13 | 0.54 | 3.00 | $--$
Algorithm 3.2(k=1) |  $<$0.01 | 0.01 | 0.02 | 0.41 | 3.03 |  $<$0.01 | 0.01 | 0.02 | 0.27 | 2.20
Algorithm 3.1(k=1) |  $<$0.01 | 0.01 | 0.02 | 0.46 | 2.98 |  $<$0.01 | 0.01 | 0.02 | 0.27 | 2.19
Algorithm 3.2(k=2) | 0.01 | 0.05 | 0.23 | 1.07 | 4.14 |  $<$0.01 | 0.04 | 0.20 | 0.80 | 3.02
Algorithm 3.1(k=2) | 0.01 | 0.04 | 0.21 | 0.93 | 4.19 |  $<$0.01 | 0.03 | 0.19 | 0.80 | 3.14

From Table 16, we find that Algorithm 3.2 and Algorithm 3.1 are much faster
than the Argyris finite element method due to the availability of fast Stokes
solvers and Poisson solvers. Moreover, Algorithm 3.2 is faster than Algorithm
3.1, since Algorithm 3.1 has one extra Poisson problem to compute. The results
in Table 17 indicate that Algorithm 3.2 and Algorithm 3.1 use much less memory
compared with the Argyris finite element method. We notice that there are not
too many differences in memory usage for both Algorithm 3.2 and Algorithm 3.1,
which is because both algorithms are serial, and the Stokes solver accounts
for the maximum memory usage.

##  Appendix A Error bounds of finite element approximations to the Stokes
problem

In this section, we only prove the error bound of the Taylor-Hood
approximations to the Stokes problem, the error bound of the the Mini element
approximations can be proved similarly. We introduce the following operators.

$\displaystyle\mathcal{B}=$ $\displaystyle-\text{div}\ :\
[H_{0}^{1}(\Omega)]^{2}\rightarrow(L_{0}^{2}(\Omega))^{\prime}=L_{0}^{2}(\Omega),\quad\langle\mathcal{B}\mathbf{v},q\rangle=-(\text{div
}\mathbf{v},q),$ (A.1) $\displaystyle\mathcal{B}^{\prime}=$
$\displaystyle\nabla\ :\
L_{0}^{2}(\Omega)\rightarrow[H^{-1}(\Omega)]^{2},\quad\langle\mathbf{v},\mathcal{B}^{\prime}q\rangle=-(\text{div
}\mathbf{v},q).$

We denote $\mathcal{B}_{n}:[V_{n}^{k}]^{2}\rightarrow(S_{n}^{k-1})^{\prime}$
be the discrete counterpart of the operator $\mathcal{B}$ and it satisfies

$\langle\mathcal{B}_{n}\mathbf{v},q\rangle=\langle\mathcal{B}\mathbf{v},q\rangle=-(\text{div
}\mathbf{v},q)\quad\forall(\mathbf{v},q)\in[V_{n}^{k}]^{2}\times S_{n}^{k-1}.$

The nullspace of $\mathcal{B}_{n}$ is given by