The a posteriori error estimates of the Ciarlet-Raviart mixed finite element method for the biharmonic eigenvalue problem

1.   Introduction

The biharmonic equation/eigenvalue problem is a fundamental model in mathematics and physics, and many numerical methods for these problems have been developed. Among these methods, the Ciarlet-Raviart mixed finite element method ^[1] is popular and classical, and it has been applied to the biharmonic equation (see ^{[2,3,4,5,6,7]}, etc.), the biharmonic eigenvalue problem (see ^{[8,9,10,11,12]}, etc.), and the transmission eigenvalue problem which has a similar structure with the biharmonic eigenvalue problem (see ^[13,14], etc.).

In practical calculations, in order to obtain high-precision approximations, a posteriori error estimation and adaptive algorithms have been widely applied (such as those in introductory textbooks ^[15,16] and review article ^[17]). For the biharmonic eigenvalue problem, Li and Yang ^[18] gave $C^0$IPG adaptive algorithms. Under the condition that the eigenfunctions $u$ and $v = \Delta u$ have the same regularity, Wang et al. ^[10] proposed a mixed discontinuous Galerkin (denoted as DG mixed) approximation scheme, and got the residual-based a posteriori error estimator of the approximate eigenpair. Feng et al. ^[19] proposed the reliable residual-based a posteriori error estimator of the approximate eigenvalue under the condition that the eigenfunction $u$ and $v = \Delta u$ have different regularity. This paper aims to study the a posteriori error estimation and adaptive algorithms of the Ciarlet-Raviart mixed conforming finite element method (denoted as the C-R mixed method) for the biharmonic eigenvalue problem. Discontinuous Galerkin methods are also effective methods for solving the biharmonic eigenvalue problem (see ^[10,19]) and they have advantages for irregular regions as they preserve local conservative properties and allow hanging nodes in the mesh adaption. But, on the same adaptive mesh without hanging nodes, the C-R mixed method has much fewer degrees of freedom than the DG mixed method. For the biharmonic eigenvalue problem on convex polygons, the C-R mixed method is simple and efficient. However, we have not seen literature on the a posteriori error analysis of this method.

As we know, the finite element method and its error estimates for an eigenvalue problem are based on the finite element method and its error estimates for the corresponding source problem. For the biharmonic equations, Charbonneau et al. ^[20] explored the residual-based a posteriori error estimate of the C-R mixed method, and Gudi ^[21] further studied the a posteriori error estimate under the condition that there are no quasi-uniformity assumptions on the triangulation.

In this paper, we extend the a posteriori error analysis of the biharmonic equation in ^[21] to the eigenvalue problem, prove the reliability and efficiency of the estimator of the approximate eigenfunction, use the error identity (2.15) to study the a posteriori error estimates of the approximate eigenvalues, and analyze the reliability of the error estimator of the approximate eigenvalues. We also implement adaptive computation. Numerical experiments indicate that our method is efficient and can get an approximate solution with high accuracy.

The organization of this paper is as follows. In the next section, we introduce the biharmonic eigenvalue problem and its C-R mixed approximation. In Section 3, we discuss the a posteriori error estimates. Finally, we present some numerical experiments to validate our theoretical results.

In this paper, $C$ represents a generic positive constant independent of the mesh size $h$, which may not be the same constant in different places. For simplicity, we use the symbol $a\lesssim b$ to mean that $a\leq C b$.

2.   Preliminaries

Consider the biharmonic eigenvalue problem

where $\Omega \subset \mathbb{R}^{2}$ is a bounded convex polygonal domain with boundary $\partial \Omega$, and $\nu$ is the unit outward normal to $\partial\Omega$.

Let $v = \Delta u$. We can rewrite the forth-order problem (2.1) as a system of second-order problems:

Multiplying the first and the second equations of (2.2) by test functions $\psi$ and $\varphi$, respectively, integrating by parts and using the boundary conditions, we can obtain the following C-R mixed variational form of (2.1): find $(\lambda, u, v)\in \mathbb{R}\times H^1_0(\Omega) \times H^1(\Omega)$ such that $\|u\|_{0} = 1$ and

where the bilinear forms are defined as follows:

In this paper, we assume $D\subseteq \Omega$. Let $H^{\rho}(D)$ denote the standard Sobolev space on $D$ with norm $\|\cdot\|_{\rho, D}$, seminorm $|\cdot|_{\rho, D}$, and $H^{0}(D) = L^{2}(D)$. When $D = \Omega$, $\|\cdot\|_{\rho, \Omega}$ and $|\cdot|_{\rho, \Omega}$ are simply denoted by $\|\cdot\|_{\rho}$ and $|\cdot|_{\rho}$, respectively. Let $H^{\rho}(\partial D)$ denote the Sobolev space on $\partial D$ with norm $\|\cdot\|_{\rho, \partial D}$ and seminorm $|\cdot|_{\rho, \partial D}$.

Assume that $\mathcal{J}_{h} = \{\kappa\}$ is a family of regular triangulation of $\Omega$ (see ^[2]). Let $h_{\kappa}$ be the diameter of $\kappa$ and $h = \max\{h_{\kappa}:\kappa\in\mathcal{J}_{h}\}$. The set of interior edges in $\mathcal{J}_{h}$ is denoted by $\Gamma_{I}$ and the set of boundary edges is denoted by $\Gamma_{B}$. Set $\Gamma = \Gamma_{I} \cup \Gamma_{B}$. Denote the length of any edge $e\in \Gamma$ by $|e|$. For any $e\in\Gamma_{I}$ and $e = \partial \kappa^+\bigcap \partial\kappa^-$, the jump of the derivative of $\psi\in V_h$ on $e$ is defined as

where $\nu$ denotes a unit normal vector on $e$, which is directed outward from $\kappa^+$; for $e\in \Gamma_{B} = \partial\kappa \bigcap \partial\Omega$,

where $\nu$ denotes a unit normal vector directed outward from the boundary $\partial\Omega$.

Define the finite element spaces as

where $P_m(\kappa)$ is the space of polynomials of degree $\leq m\; (m\geqslant 2)$.

Define the broken Sobolev space

with the mesh-dependent norm

Define the following norm on product space $W = L^2(\Omega)\times H^2(\Omega, \mathcal{J}_{h})$ as

Based on the mixed formulation (2.3) and (2.4), we can get the C-R mixed finite element approximation: find $(\lambda_h, u_h, v_h)\in \mathbb{R}\times V^0_h \times V_h$, $\|u_{h}\|_{0} = 1$, such that

Consider the following fourth-order problem:

We assume the following regularity assumption is valid:

For given $g\in L^{2}(\Omega)$, there is a unique solution $(\omega, \varphi)\in H^{2}_{0}(\Omega)\times H^{1}(\Omega)$ to the problem (2.9) satisfying the following elliptic regularity estimate:

When $\Omega$ is a smooth domain, (2.10) is valid. However, when $\Omega \subset {\mathbb{R}}^{2}$ is a bounded convex domain, Grisvard ^[22] only stated that $\Delta^2:H^3(\Omega)\rightarrow H^{-1}(\Omega)$ is isomorphic, and Blum et al. ^[23] stated that (2.10) is true if the maximum interior angle of $\Omega$ is less than $126.283696\cdots$. This assumption is made only to reduce the technical complexity of the error analysis.

Let $\lambda$ and $\lambda_{h}$ be the $k$th eigenvalue of (2.3), (2.4) and (2.7), (2.8), respectively. The algebraic multiplicity of $\lambda$ is $q$, $\lambda = \lambda_{k} = \lambda_{k+1} = ... = \lambda_{k+q-1}$. Let $V_{\lambda}$ denote the space spanned by all eigenfunctions corresponding to $\lambda$, and let $V_{\lambda}(h)$ denote the space spanned by all eigenfunctions corresponding to the eigenvalues $\lambda_{j, h}$ that converge to $\lambda$.

Lemma 2.1. Let $\lambda$ be the $k$th eigenvalue of (2.3) and (2.4), $V_{\lambda} \subset H^{m+1}(\Omega)$, and $\left(\lambda_{h}, v_{h}, u_{h}\right)$ be the $k$th eigenpair of (2.7) and (2.8) with $\left\|u_{h}\right\|_{0} = 1$, then there exists an eigenfunction $(v, u)$ corresponding to $\lambda$, such that $\|u\|_0 = 1$ and

where $\varepsilon = 0$ when $m = 2$ and $\varepsilon = 1$ when $m\geqslant3$. Let $u \in V_{\lambda}$ and $\|u\|_{0} = 1$, then there exists $u_{h} \in V_{\lambda}(h)$ such that $\left\|u-u_{h}\right\|_{1} \lesssim h^{m}$.

Proof. We know that (2.11), (2.12) and (2.14) are valid from Theorem 11.4 in ^[8]. We obtain the conclusion (2.13) from ^[4].

Lemma 2.2. Suppose $(\lambda, u, v)$ and $(\lambda_{h}, u_{h}, v_{h})$ are the eigenpairs of (2.3), (2.4) and (2.7), (2.8), respectively. Then

Proof. By (2.3) and (2.4) we deduce that

By (2.7) and (2.8) we have

Then, dividing by $-(u_{h}, u_{h})$ on both sides of (2.16), we obtain (2.15).

To discuss the error estimates, we state some results on the approximation properties of interpolation in ^[24] without proof, which will play a crucial role in our analysis.

Lemma 2.3. For any $\phi\in H_{0}^{2}(\Omega)$, let $\phi_{h}\in V_{h}$ be the Lagrange interpolant of $\phi$. Then, for any $\kappa\in\mathcal{J}_{h}$, there exists a positive constant $C$ which is independent of $h$ such that

Denote the piecewise (element-wise) Laplacian of $v\in V_{h}$ by $\Delta_{h}v$.

Lemma 2.4. For all $q_{h} \in V_{h}$ there exists a positive constant $C$ independent of $h$ such that

where $E_{h}: V_{h}\rightarrow \widetilde{V_{h}}\subset H_{0}^2(\Omega)$ is a recovery operator defined as in ^[21], $\widetilde{V_{h}}$ is a Hsieh-Clough-Tocher (HCT) finite element space associated with $\mathcal{J}_{h}$.

Proof. Charbonneau et al. ^[20] and Gudi ^[21] proved the above conclusion for $m = 2$ and $3$. From Lemma 1 in ^[25], we know the above conclusions are valid for $m\geq 2$.

3.   A posteriori error estimates

Based on the a posteriori error analysis of the source problem corresponding to the biharmonic eigenvalue problem (2.1) in ^[21], the local estimator can be defined as follows:

For $\kappa \in \mathcal{J}_{h}$,

for $e \in \Gamma_{I}$,

and for $e \in \Gamma$

Let

and

We can get the following theorem.

Theorem 3.1. Let $(\lambda, u, v)$ and $(\lambda_{h}, u_{h}, v_{h})$ be the $k$th eigenpairs of (2.3), (2.4) and (2.7), (2.8), respectively. Then it holds that

Proof. From the definitions of the norm $\|\cdot\|_W$ and $|||\cdot|||$, we know that

Now we estimate $\||u-u_{h}\||$. Since $[\frac{\partial u}{\partial \nu}] = 0$ on $e$, we have

Using the triangle inequality and Lemma 2.4 we obtain

Note that by the dual argument we have

Let $\phi \in H_{0}^{2}(\Omega)$. Then

Let $\phi_{h} \in V^0_{h}$ be the Lagrange interpolant of $\phi$, then we can deduce that

Using the Cauchy-Schwarz inequality and Lemma 2.3, we know

and

Substituting (3.9)–(3.11) into (3.8), we obtain

Using the triangle inequality and Lemma 2.4, we obtain

Substituting (3.12) and (3.13) into (3.7), and using (3.6), we deduce

Then, from (3.3)–(3.5) and (3.14), we can get

Using the triangle inequality (3.5) and (3.14), we obtain

The proof is complete.

The following theorem gives the error bounds for the approximate eigenvalue.

Theorem 3.2. Let $(\lambda, u, v)$ and $(\lambda_{h}, u_{h}, v_{h})$ be the $k$th eigenpairs of (2.3), (2.4) and (2.7), (2.8), respectively. Then it holds that

where $I_{h}v\in V_{h}$ is the Lagrange interpolant of $v$.

Proof. From (2.3) and (2.7), we get

Thus, using (2.15) and integrating by parts, we deduce that

Using the definition of norm $\|\cdot\|_{W}$ and (3.1), we can get (3.15). The proof is complete.

Now, based on ^[16,21] we study the efficiency of the error estimator.

Let $e$ represent a common edge shared by the two elements $\kappa^{+}$ and $\kappa^{-}$, and denote $\omega_{e} = \kappa^{+}\cup\kappa^{-}$.

Theorem 3.3. Let $(\lambda, u, v)$ and $(\lambda_{h}, u_{h}, v_{h})$ be the $k$th eigenpairs of (2.3), (2.4) and (2.7), (2.8), respectively. Then it holds that

Proof. Using bubble function techniques (see ^[16,21]), we first estimate (3.17).

Let $b_{\kappa}\in H_{0}^{2}(\kappa)$ be a bubble polynomial defined on $\kappa$. Then

Let $\phi = b_{\kappa}(\lambda_{h} u_{h}-\Delta v_{h})$. Then

Integrating by parts twice and using the inverse inequality, we get

Combining the above three estimates, we get (3.17).

In the proof of Lemma 3.3 in ^[21], let $f = \lambda_{h} u_{h}$, then we can get (3.18).

It is clear that

and using (3.17), (3.18) and the definition of norm $\|\cdot\|_{W}$, we can get (3.19). The proof is complete.

Remark 3.1. From Lemma 2.1, we know that $\|u_h-u\|_0$ is a higher-order term than $\|\Delta u-v_h\|_0$. And, interpolation theory shows that the estimate of the error $\sum\limits_{\kappa}\sum\limits_{j = 0}^{2}h_{\kappa}^{2j}\|I_{h}v-v\|_{j, \kappa}^{2}$ is optimal with respect to $h$, so we can expect to get

So, substituting (3.21) into (3.15), we obtain

Therefore, the estimator $\eta_{h}^{2}(\Omega)$ of the eigenvalue error $|\lambda_{h}-\lambda|$ is reliable up to the higher-order term $\lambda \|u_h-u\|^2_{0}$.

4.   Numerical experiments

In this section, we will present some numerical results to validate our theoretical analysis. We calculate the smallest eigenvalue of the biharmonic eigenvalue problem on adaptive meshes in three domains: the unit square $\Omega_S = (0, 1)^2$, the regular hexagon $\Omega_H$ with side length of 1, and the L-shaped domain $\Omega_L = (-\frac{1}{2}, \frac{1}{2})^2/[0, \frac{1}{2})\times(-\frac{1}{2}, 0]$. For $\Omega_S$, we choose the reference value $\lambda_{1}\approx1294.93397959171$ (see ^[26]), and take the reference value $\lambda_{1}\approx 163.59756815825$ in $\Omega_{H}$ and $\lambda_{1}\approx 6703.6047044786$ in $\Omega_{L}$ (see ^[19]).

The computations are implemented according to the following algorithm, and for $\Omega_{S}$ our calculations refer to Algorithm 2 in ^[18] when the P4 element is used. All computations are easily realized under the packages of the FEM ^[27,28].

The adaptive algorithm of the mixed conforming finite element method:

Choose the parameter $0 < \theta < 1$.

Step 1. Pick any initial mesh $\mathcal{J}_{h_{0}}$ with initial mesh size $h_{0}$.

Step 2. Solve (2.7)-(2.8) on $\mathcal{J}_{h_{0}}$ for discrete solution $(\lambda_{h_{0}}, u_{h_{0}}, v_{h_{0}})$.

Step 3. Let iterations $l = 0$.

Step 4. Compute the local estimator $\eta_{h_{l}}(\kappa)$.

Step 5. Construct $\widehat{\mathcal{J}_{h_l}}\subset \mathcal{J}_{h_l}$ by Marking Strategy E and parameter $\theta$.

Step 6. Refine $\mathcal{J}_{h_l}$ to get a new mesh $\mathcal{J}_{h_{l+1}}$ by procedure REFINE.

Step 7. Solve (2.7)-(2.8) on $\mathcal{J}_{h_{l+1}}$ for discrete solution $(\lambda_{h_{l+1}}, u_{h_{l+1}}, v_{h_{l+1}})$.

Step 8. Let $l\Leftarrow l+1$ and go to Step 4.

Marking Strategy E:

Step 1. Construct a minimal $\widehat{\mathcal{J}_{h_l}}\subset \mathcal{J}_{h_l}$ by selecting some elements in $\mathcal{J}_{h_l}$ such that

Step 2. Mark all elements in $\widehat{\mathcal{J}_{h_l}}$.

The value of $\theta$ is set to 0.5. The results computed by the adaptive algorithm with P2, P3 and P4 elements in $\Omega_{S}$, $\Omega_{H}$ and $\Omega_{L}$ are listed in Tables 1–3, respectively. We also depict the curves of absolute error $|\lambda_{h}-\lambda_{1}|$ in the three domains in Figures 1–3 and show the adaptive meshes obtained by P2, P3 and P4 elements in Figures 4–6.

For $\Omega_{S}$, from Table 1 we can obverse that the approximate eigenvalues of high accuracy can be obtained when using higher degree polynomials. From Table 4, compared with the results obtained by the DG mixed method in ^[19], we can conclude that with the same degree of freedom, using the mixed conforming finite element method can achieve higher accuracy. And, compared with the results calculated in ^[11], we can conclude that with the same degree of freedom, the approximations obtained by the adaptive algorithm with P3 element have higher precision than those computed by the C-R mixed method with P3 element on uniform meshes.

Figure 1 shows that the error curves are approximately parallel to the line with slope $-2$, $-3$ and $-4$, and the algorithm can achieve the optimal convergence order $O(dof^{-2})$, $O(dof^{-3})$ and $O(dof^{-4})$ when P2, P3 and P4 elements are used, respectively. This means that the results obtained in numerical experiments have higher order convergence than theoretical analysis, and we think the reason is that $\Delta u \in H^2(\Omega)$ when $u\in H^4(\Omega)$, thus the regularity of $v = \Delta u$ is underestimated in the theoretical analysis of the C-R mixed method.

For $\Omega_{H}$ and $\Omega_{L}$, we can observe similar conclusions. Although we only analyze the C-R mixed method for convex or smooth domains, we also implement adaptive calculations in the L-shaped domain, and the results in Table 3 and Figure 3 indicate that our method is still convergent.

Remark 4.1. There are usually two ways to determine when to terminate the iteration. One is by the error estimator. The adaptive procedure will continue until the error estimator is less than a prefixed tolerance. The other is by the difference between adjacent two or several iterations. When the difference is less than a prefixed tolerance, the iteration will be terminated. However, in this paper, since our error estimator is not asymptotically accurate and the error curves fluctuate, we judge whether the calculation result is accurate by observing the changing trend of the error.

5.   Conclusions

In this paper, we study the a posteriori error estimates and adaptive calculation of the C-R mixed method for the biharmonic eigenvalue problem on convex polygon domains. We propose a posteriori error estimators, prove the reliability and efficiency of the error estimator of the approximate eigenfunction, and analyze the reliability of the error estimator of the approximate eigenvalues. Numerical experiments confirm our theoretical analysis and indicate that our adaptive algorithm is efficient. Meanwhile, the results in Table 3 and Figure 3 show that the C-R mixed method in adaptive fashion is convergent and efficient on nonconvex domains. It is a challenging and valuable work to prove the convergence of C-R mixed method on nonconvex domains.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors cordially thank the editor and the referees for their valuable comments and suggestions that lead to the improvement of this paper.

This work was supported by the National Natural Science Foundation of China (Nos. 11561014 and 11761022).

Conflict of interest

The authors declare that this work does not have any conflicts of interest.