A posteriori error analysis for the mixed Laplace eigenvalue problem

TL;DR

This paper develops a posteriori error estimates for mixed finite element methods solving the Laplace eigenvalue problem, providing reliable and efficient error bounds with local reconstructions to improve numerical accuracy.

Contribution

It introduces a new local reconstruction technique for the primal variable, enabling effective a posteriori error estimation in mixed Laplace eigenvalue approximations.

Findings

Error estimator provides an upper bound for the approximation error.

The error estimator is shown to be locally efficient.

Reconstruction is performed on vertex patches for improved accuracy.

Abstract

This paper derives a posteriori error estimates for the mixed numerical approximation of the Laplace eigenvalue problem with homogeneous Dirichlet boundary conditions. In particular, the resulting error estimator constitutes an upper bound for the error and is shown to be local efficient. Therefore, we present a reconstruction in the standard $H^1_0$-conforming space for the primal variable of the mixed Laplace eigenvalue problem. This reconstruction is performed locally on a set of vertex patches.

Figures

Equations

$$\{\begin{array}{rlrl}& (\mathbf{\sigma },\mathbf{\tau })+(\mathrm{div}\mathbf{\tau },u)=0& & \forall \mathbf{\tau }\in \mathbf{H}(\mathrm{div};\Omega )\\ & (\mathrm{div}\mathbf{\sigma },v)=-\lambda (u,v)& & \forall v\in L^2(\Omega ).\end{array}$$

$${\Sigma }_h=\{\mathbf{\tau }\in \mathbf{H}(\mathrm{div};\Omega ):\mathbf{\tau }{|}_K\in [P_k(K){]}^d\oplus \mathbf{x}{\tilde{P}}_k(K)\forall K\in {\mathcal{T}}_h\},$$

$$U_h=\{v\in L^2(\Omega ):v{|}_K\in P_k(K)\forall K\in {\mathcal{T}}_h\}.$$

$$\{\begin{array}{rlrl}& ({\mathbf{\sigma }}_h,\mathbf{\tau })+(\mathrm{div}\mathbf{\tau },u_h)=0& & \forall \mathbf{\tau }\in {\Sigma }_h\\ & (\mathrm{div}{\mathbf{\sigma }}_h,v)=-{\lambda }_h(u_h,v)& & \forall v\in U_h.\end{array}$$

$$|\lambda -{\lambda }_h|\le C{h}^{2r}|u{|}_{H^{r+1}(\Omega )}$$

$$\Vert u-u_h{\Vert }_{L^2(\Omega )}+\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}\le C{h}^r|u{|}_{H^{r+1}(\Omega )}$$

$$\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{\mathbf{H}(\mathrm{div};\Omega )}\le C{h}^r|u{|}_{H^{r+1}(\Omega )}$$

$$\lambda -{\lambda }_h=\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}^2-{\lambda }_h\Vert u-u_h{\Vert }_{L^2(\Omega )}^2.$$

$$\Vert {\mathcal{P}}_{h}u-u_h{\Vert }_{L^2(\Omega )}\le \rho (h)\left(\Vert u-u_h{\Vert }_{L^2(\Omega )}+\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}\right).$$

$$U_h^{\ast }=\{v\in L^2(\Omega ):v{|}_K\in P_{k+1}(K)\forall K\in {\mathcal{T}}_h\}.$$

$$\{\begin{array}{rl}& P_h{u}_h^{\ast }=u_h\\ & (\nabla u_h^{\ast },\nabla v{)}_K=({\mathbf{\sigma }}_h,\nabla v{)}_K\forall v\in (I-P_h)U_h^{\ast }{|}_K\forall K\in {\mathcal{T}}_h.\end{array}$$

$${\left(\sum _{K\in {\mathcal{T}}_h}\Vert \nabla u-\nabla u_h^{\ast }{\Vert }_{L^2(K)}^2\right)}^{1/2}$$

$$\Vert u-u_h^{\ast }{\Vert }_{L^2(\Omega )}$$

$$U_h^{\ast \ast }=\{v\in H_0^1(\Omega ):v{|}_K\in P_{k+1}(K)\forall K\in {\mathcal{T}}_h\}$$

$$\Vert \nabla u-\nabla u_h^{\ast \ast }{\Vert }_{{\mathbf{L}}^2(\Omega )}\le C{h}^r|u{|}_{H^{r+1}(\Omega )}$$

$$\Vert u-u_h^{\ast \ast }{\Vert }_{L^2(\Omega )}\le C\rho (h)h^r|u{|}_{H^{r+1}(\Omega )}.$$

$$\eta (K)=\Vert \nabla u_h^{\ast \ast }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(K)}$$

$$\eta ={\left(\sum _{K\in {\mathcal{T}}_h}\eta (K{)}^2\right)}^{1/2}=\Vert \nabla u_h^{\ast \ast }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}.$$

$$\eta (K)\le \sqrt{2}{\left(\Vert \nabla u-\nabla u_h^{\ast \ast }{\Vert }_{{\mathbf{L}}^2(K)}+\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(K)}\right)}^{\frac{1}{2}}.$$

$$\eta (K)=\Vert \nabla u_h^{\ast \ast }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(K)}\le \Vert \nabla u-\nabla u_h^{\ast \ast }{\Vert }_{{\mathbf{L}}^2(K)}+\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(K)}.$$

$${\left(\Vert \nabla u-\nabla u_h^{\ast \ast }{\Vert }_{{\mathbf{L}}^2(\Omega )}^2+\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}^2\right)}^{\frac{1}{2}}\le \eta +hot(h),$$

$$\Vert \nabla u_h^{\ast \ast }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}^2$$

$$=\Vert \nabla u_h^{\ast \ast }-\nabla u{\Vert }_{{\mathbf{L}}^2(\Omega )}^2+\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}^2$$

$$+2(\nabla u_h^{\ast \ast }-\nabla u,\mathbf{\sigma }-{\mathbf{\sigma }}_h{)}_{{\mathbf{L}}^2(\Omega )}$$

$$=\Vert \nabla u_h^{\ast \ast }-\nabla u{\Vert }_{{\mathbf{L}}^2(\Omega )}^2+\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}^2$$

$$-2(u_h^{\ast \ast }-u,\mathrm{div}(\mathbf{\sigma }-{\mathbf{\sigma }}_h){)}_{L^2(\Omega )}.$$

$$\Vert \nabla u-\nabla u_h^{\ast \ast }{\Vert }_{{\mathbf{L}}^2(\Omega )}^2+\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}^2\le {\eta }^2+2\Vert u-u_h^{\ast \ast }{\Vert }_{L^2(\Omega )}\Vert \lambda u-{\lambda }_h{u}_h{\Vert }_{L^2(\Omega )}.$$

$$\Vert \lambda u-{\lambda }_h{u}_h{\Vert }_{L^2(\Omega )}$$

$$\le \lambda \Vert u-u_h{\Vert }_{L^2(\Omega )}+|\lambda -{\lambda }_h|\Vert u_h{\Vert }_{L^2(\Omega )}$$

$$\le C(h^r+h^{2r})|u{|}_{H^{r+1}(\Omega )},$$

$$\Vert \nabla u-\nabla u_h^{\ast \ast }{\Vert }_{{\mathbf{L}}^2(\Omega )}^2+\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}^2\le {\eta }^2+O(\rho (h)h^{2r}).$$

Full text

1 Introduction

Eigenvalue problems arise in countless applications and the precise numerical approximation of eigenvalues and eigenvectors of elliptic operators is therefor crucial. Finite element methods for these problems have been widely used and analyzed under a general framework in many works, e.g. in BabOsb:91 and the references therein. The Laplace problem was shown to be an interesting model problem worth to start with. Convergence and optimal a priori estimates for both eigenvalues and eigenfunctions are given e.g. in acta ; bbf .

It is well-known that on general domains, the eigenfunction can not excepted to be sufficiently regular to use these a priori estimates. Therefore, it is important to use adaptive procedures based on a posteriori error estimators in order to retain the optimal convergence order and several approaches have been considered to construct estimators based on the residual equation presented in AinOde:93 and MR3059294 . In particular DurPadRod:03 , obtained an a posteriori error estimator for the linear finite element approximation of second order elliptic eigenvalue problems using a residual type error estimator.

For the study of the optimal convergence rates of these procedures, important progress has been made during the last decades. The crucial marking provided in MR1393904 allows the error measured in the energy norm to decrease at a constant rate in each step until a prescribed error bound is reached. MR1770058 showed that the strict error reduction used therfor cannot be expected in general and proved convergence by introducing the concept of data oscillation and the interior node property. MR2324418 provided a new overall theoretical understanding of adaptive FEM in order to realize optimal computational complexity. MR2421046 proved the optimal cardinality of the AFEM using a decay between consecutive loops. A survey and an axiomatic presentation of the proof of optimal convergence rates for adaptive finite element methods can be found in MR2648380 ; MR3076038 ; CarFeiPagPra:14 . For non-conforming elements, see MR2595052 .

An alternative approach for a posteriori error estimation is based on the hypercirle identity dating back at least as far as LadLeg:83 and PraSyn:47 . In fact, using a flux reconstruction of the primal variable of the source problem usually leads to guaranteed, easily, fully, and locally computable, upper bound on the error measured in the energy norm, see e.g. BraSch:08 ; ErnVoh:15 . A unified framework for a posteriori error estimation based on stress reconstruction for the Stokes system is carried out in HanSteVoh:12 , and this was extended to the linear elasticity problem in BerMolSta:17 .

In addition to standard conforming Galerkin approximations, mixed finite elements are very popular since they provide a good approximation of the eigenvalues and eigenfunctions, see acta ; bbf . A study of optimal rates of convergence was considered in MR2772091 using a residual based error estimator. For the mixed scheme however, the hypercircle identity can be used as an alternativ approach in order to obtain a considerably better approximation for the scalar variable of the source problem, see StM2an91 . This paper aims to extend this approach to the eigenvalue.

Mixed finite element methods are an alternative to standard conforming methods. The methods yield an accurate approximation to the flux, whereas that of the scalar variables is less accurate. However, the discrete solution can be used in order to obtain a considerable better approximation to the scalar variable by local post processing schemes, as first shown in MR813687 . The postprocessing can also be used to obtain an efficient a posteriori error estimator for the source problem, cf. MR2240629 . In this paper, we will use the post processing combined with the idea behind the hypercircle method in order to design an a posteriori estimator for the eigenvalue problem. We use the Raviart-Thomas-Nedelec mixed methods since for other main family, i.e. Brezzi-Douglas-Marini, the approximation properties needed are not valid, cf. Remark LABEL:BDMrem below.

2 Setting of the problem and preliminary results

Let $\Omega$ be a polygon in $\mathbb{R}^{2}$ or a polyhedron in $\mathbb{R}^{3}$. We consider the standard mixed formulation of the Laplace eigenvalue problem with homogeneous Dirichlet boundary conditions: find $\lambda\in\mathbb{R}$ and a non vanishing $u\in L^{2}(\Omega)$, and $\boldsymbol{\sigma}\in\mathbf{H}(\operatorname*{\textrm{div}};\Omega)$, such that it holds

[TABLE]

Given a triangulation $\mathcal{T}_{h}$ of $\Omega$ and an integer $k\geq 0$, we define $\Sigma_{h}$ as the standard Raviart–Thomas space of order $k$ (sometimes called Raviart–Thomas-Nédélec space)

[TABLE]

where $\tilde{P}_{k}(K)$ is the space of homogeneous polynomials of degree $k$, and $U_{h}$ as the space of discontinuous piecewise polynomials of $k$

[TABLE]

It is well known (see for instance acta ; bbf for a review) that the following scheme provides a good approximation of the eigenvalues and eigenfunctions of (2.1): find $\lambda_{h}\in\mathbb{R}$ and a non vanishing $u_{h}\in U_{h}$, and $\boldsymbol{\sigma}_{h}\in\Sigma_{h}$, such that it holds

[TABLE]

Let $(\lambda,u,\boldsymbol{\sigma})$ be a solution of (2.1) corresponding to a simple eigenvalue. In general there exists $s>1/2$ such that $u\in H^{1+s}(\Omega)$ and $\boldsymbol{\sigma}\in\mathbf{H}^{s}(\Omega)\cap\mathbf{H}(\operatorname*{\textrm{div}};\Omega)$. We have that there exist a discrete solution $(\lambda_{h},u_{h},\boldsymbol{\sigma}_{h})$ of (2.4) such that

[TABLE]

with $r=\min(s,k+1)$.

Remark 1

Estimate (2.5) requires a careful choice of $u_{h}$ and $\boldsymbol{\sigma}_{h}$. It is well known that this result can be achieved with an appropriate normalization and choice of the sign. The case of eigenvalues with multiplicity higher than one needs more attention $($see, for instance, the discussion in (acta, , Sec. 3.1)$)$ and it has been recently observed $($see gallistl $)$ that an effective a posteriori analysis should consider clusters of eigenvalues. In bggg a complete analysis of an adaptive scheme for the problem under consideration has been performed in the framework of clusters approximation; for the sake of readability, in this paper we discuss the case of a simple eigenvalue, being understood that our result can be extended to the more realistic situation using the techniques of gallistl ; bggg .

Remark 2

In (MR1722056, , Lemma 3.2) the following equality is proved for mixed methods in general

[TABLE]

From this it follows that for the BDM method with the same space for the scalar unknown, contrary to what is the case for the source problem, the accuracy is the same as for the RT. Therefore we in this paper do not consider the BDM alternative.

Let $\mathcal{P}_{h}:L^{2}(\Omega)\to U_{h}$ denote the $L^{2}(\Omega)$ projection. It has been shown in roberts that for the source problems associated with (2.1) and (2.4) the term $\|\mathcal{P}_{h}u-u_{h}\|_{L^{2}(\Omega)}$ achieves a higher rate of convergence than the actual error $\|u-u_{h}\|_{L^{2}(\Omega)}$ (see also (bbf, , Sec. 7.4)). Not surprisingly, the same result holds true also for the eigenfunctions of  (2.1) and (2.4), even if the proof in the case of the eigenvalue problem is not a trivial extension of the original one. The following estimate has been first found in gardini (lowest order case, using the equivalence with a non-conforming approximation) and then proved in (bggg, , Sec. 6.1) (general case, notation suited for clusters of eigenvalue) and in (rodolfo1, , Sec. 3) and (rodolfo2, , Lemma 10) (these last two results are for a mixed formulation associated with the Maxwell eigenvalue problem, but the same proof applies to our situation with the natural modifications): there exists $\rho(h)$, tending to zero as $h$ goes to zero, such that

[TABLE]

The value of $\rho(h)$ depends on the regularity of our problem. In our case it is at least $\rho(h)=O(h^{1/2})$ and, when the domain is convex, we have $\rho(h)=O(h)$.

We now recall the postprocessing procedure which has been introduced in StM2an91 in order to improve the approximation of $u$ in the case of the source problem associated with (2.1). Let $U_{h}^{*}$ be the space of polynomials of degree up to $k+1$

[TABLE]

The postprocessed solution $u_{h}^{*}\in U_{h}^{*}$ is defined such that

[TABLE]

A proof analogous to the one presented in StM2an91 for the source problem, together with the estimate (2.7), can be applied to the eigenvalue problem in order to obtain the following bounds for all $k\geq 0$:

[TABLE]

where $\mathcal{E}_{h}$ is the set of edges in the triangulation $\mathcal{T}_{h}$, $h_{\ell}$ is the length of $\ell$, $[[\cdot]]$ denotes the jump, $\rho(h)$ comes from (2.7), and $r$ has been defined after Equation (2.5). In particular, if the domain is convex, the right hand side of the last estimate is $O(h^{r+1})$.

Remark 3

In StM2an91 the last estimate has been proved for $k\geq 1$. However, it was observed that if the source term belongs to the discrete space $U_{h}$, then the optimal bound holds true also in the lowest order case $($see Theorem 2.1, Theorem 2.2, and Remark 2.1$)$. In our case, the right-hand side $($corresponding to the source term in the framework of StM2an91 $)$ is equal to $\lambda_{h}u_{h}\in U_{h}$, so that the optimal bound is valid also when $k=0$.

The next step consists in building a second postprocessed solution which is globally continuous. This will allow to consider its gradient over the entire domain $\Omega$. The standard procedure consists in averaging the degrees of freedom and it is sometimes called Oswald interpolation, being a particular case of the more general case considered in oswald . For more details the reader is referred to (dipietroern, , Sec. 5.5.2), where the construction is explained together with the relevant estimates. Using the Oswald interpolation, starting from $u_{h}^{*}$ we can build an element $u_{h}^{**}$ in

[TABLE]

such that

[TABLE]

Moreover, the degrees of freedom on the boundary can be set equal to zero so that $u_{h}^{**}$ belongs to $H^{1}_{0}(\Omega)$.

3 A posteriori analysis

We define the following local error estimator

[TABLE]

and the corresponding global error estimator

[TABLE]

It is immediate to see that the above error estimator is locally efficient .

Theorem 1** (Local efficiency)**

Let $(\lambda,u,\boldsymbol{\sigma})$ be a solution of (2.1) corresponding to a simple eigenvalue $($see Remark 1$)$ and $(\lambda_{h},u_{h},\boldsymbol{\sigma}_{h})$ the associated discrete solution in the spirit of estimate (2.5). Then the following estimate holds true for each element $K$ in $\mathcal{T}_{h}$

[TABLE]

Proof

The proof is immediate from $\operatorname*{\boldsymbol{\nabla}}u=\boldsymbol{\sigma}$ and the triangular inequality

[TABLE]

The following theorem provides a reliability estimate with constant equal to one.

Theorem 2** (Global reliability)**

Let $(\lambda,u,\boldsymbol{\sigma})$ be a solution of (2.1) corresponding to a simple eigenvalue $($see Remark 1$)$ and $(\lambda_{h},u_{h},\boldsymbol{\sigma}_{h})$ the associated discrete solution in the spirit of estimate (2.5). Then the following estimate holds true

[TABLE]

where $hot(h)$ is a higher order term with respect to $h^{r}$ as $h$ goes to zero.

Proof

From $\operatorname*{\boldsymbol{\nabla}}u=\boldsymbol{\sigma}$ we have

[TABLE]

Since $\operatorname*{\textrm{div}}\boldsymbol{\sigma}=-\lambda u$ and $\operatorname*{\textrm{div}}\boldsymbol{\sigma}_{h}=-\lambda_{h}u_{h}$, we obtain

[TABLE]

We now conclude the proof by showing that the last term of the previous equation is a higher order term.

We have already seen in (2.9) that $\|u-u_{h}^{**}\|_{L^{2}(\Omega)}$ is of order $O(\rho(h)h^{r})$. The other term can be estimated as follows

[TABLE]

so that we finally get

[TABLE]

4 Numerical illustration

This section is devoted to numerical results using uniform and adaptive $h$-refinements of triangular meshes for two-dimensional domains.

The aim of our tests is twofold: on one hand we want to confirm our theory developed in the previous section; on the other hand we want to investigate the convergence of the adaptive scheme based on the error indicator that we have introduced.

We compute the solution of problem (2.4) with Raviart–Thomas elements with $k=0$ and $k=1$.

The first computations are carried on the square domain $[0,2\pi]^{2}$ with homogeneous Dirichlet boundary conditions imposed on the entire boundary, similarly to acta . In this case, the smallest eigenvalue of (2.1) is known to be $\lambda=2$. A corresponding eigenfunction is given by $u(x,y)=\sin x\sin y$. Since the eigenfunction is smooth, we can take $s=2$ in the error estimates (2.5). The approximation of the eigenvalue is illustrated in Figure 1. Since on a uniform mesh the number of degrees of freedom $N$ is proportional to $h^{-2}$, we can see that the eigenvalues converge with optimal order $O(h^{2})$ and $O(h^{4})$ for $k=0$ and $k=1$, respectively. In Figure 2, we see that the error estimator converges with the same rate. The corresponding convergence rates are summarized in Figure 3.

In order to investigate the efficiency of the error estimator with respect to the natural norms, the efficiency index is illustrated in Figure 4.

In Figure 5 we check the validity of Theorem 2: in particular, we compute the rate of the higher order terms involved with the reliability estimate. More precisely, the lines in Figure 5 refer to the difference between the left and right hand sides of equation (3.1).

In order to illustrate the adaptive scheme, a second test case is concerned with the L-shaped domain $[-1,1]^{2}\backslash[0,1]^{2}$. Since the Sobolev regularity is restricted by the biggest interior angle of the domain (see classical regularity theory e.g. KozMazRos:97 ), the lowest eigenfunction can only be shown to be in $H^{s}$ for $s<4/3$ for our L-shaped case. The convergence rate for the approximation of $\lambda$ for a uniform refinement will be at best ${N}^{-2/3}$, as illustrated in Figure 6. There, we also see that the adaptive scheme is able to retain the optimal convergence order. Similarly, the convergence rate for the approximation of the eigenfunction for a uniform refinement will be at best ${N}^{-1/3}$. In Figure 7, we see that the error estimator converges with the optimal convergence rate. Figure 8 shows that the adaptive strategy refines the mesh where the singularity is expected, i.e. in the reentrant corner.