Solving Eigenvalue Problems in a Discontinuous Approximation Space by Patch Reconstruction

TL;DR

This paper introduces a high-order discontinuous Galerkin method using patch reconstruction for solving elliptic eigenvalue problems, achieving high accuracy with fewer degrees of freedom and outperforming traditional methods in certain cases.

Contribution

It extends a patch reconstructed DG method to eigenvalue problems, demonstrating high accuracy and efficiency with minimal degrees of freedom compared to classical methods.

Findings

Higher order methods yield more reliable eigenvalues.

The method attains similar accuracy with fewer degrees of freedom.

Performance surpasses conforming finite element methods at higher polynomial orders.

Abstract

We adapt a symmetric interior penalty discontinuous Galerkin method using a patch reconstructed approximation space to solve elliptic eigenvalue problems, including both second and fourth order problems in 2D and 3D. It is a direct extension of the method recently proposed to solve corresponding boundary value problems, and the optimal error estimates of the approximation to eigenfunctions and eigenvalues are instant consequences from existing results. The method enjoys the advantage that it uses only one degree of freedom on each element to achieve very high order accuracy, which is highly preferred for eigenvalue problems as implied by Zhang's recent study [J. Sci. Comput. 65(2), 2015]. By numerical results, we illustrate that higher order methods can provide much more reliable eigenvalues. To justify that our method is the right one for eigenvalue problems, we show that the patch reconstructed approximation space attains the same accuracy with fewer degrees of freedom than classical discontinuous Galerkin methods. With the increasing of the polynomial order, our method can even achieve a better performance than conforming finite element methods, such methods are traditionally the methods of choice to solve problems with high regularities.

Tables

Table 1. Table 4.8. The first eigenvalues of the second order and biharmonic equation in L-shaped domain.

Order	DOFs	N=2.00E+2	N=8.00E+2	N=3.20E+3	N=1.28E+4	N=5.12E+4
$m=1$	Laplace	10.786	9.9562	9.7396	9.6867	9.6733
$m=2$		10.403	9.7422	9.6780	9.6707	9.6692
$m=3$		9.8128	9.6811	9.6724	9.6700	9.6691
$m=2$	Biharmonic	179.98	171.55	167.78	166.75	165.43
$m=3$		168.82	166.78	165.77	165.12	164.68

Table 2. Table 4.9. The first eigenvalues of the Laplace problem in 3D, λ 1 = 3 ​ π 2 ​ ( 29.61 ) subscript 𝜆 1 3 superscript 𝜋 2 29.61 \lambda_{1}=3\pi^{2}(29.61) .

Order	Mesh Size	$h=$2.500E-1	$h=$1.250E-1	$h=$6.250E-2	$h=$3.125E-3
$m=1$	Value	45.43	34.99	31.03	29.97
	Error	5.33E-1	1.81E-1	4.82E-2	1.24E-2
	Order	-	1.55	1.91	1.96
$m=2$	Value	35.57	29.96	29.63	29.61
	Error	2.01E-1	1.21E-2	7.70E-4	4.98E-5
	Order	-	4.10	3.96	3.95
$m=3$	Value	31.34	29.63	29.61	29.61
	Error	5.85E-2	6.79E-4	1.07E-5	1.64E-7
	Order	-	6.42	5.99	6.02
$m=4$	Value	30.23	29.61	29.61	29.61
	Error	2.12E-2	8.24E-5	3.23E-7	1.23E-9
	Order	-	8.03	7.99	7.93

Table 3. Table 4.10. The first eigenvalues of the biharmonic problem in 3D, λ 1 = 9 ​ π 4 ​ ( 876.68 ) subscript 𝜆 1 9 superscript 𝜋 4 876.68 \lambda_{1}=9\pi^{4}(876.68) .

Order	Mesh Size	h=2.500E-1	h=1.250E-1	h=6.250E-2	h=3.125E-3
m=2	Value	1000.54	906.41	883.20	878.20
	Error	1.41E-1	3.39E-2	7.44E-3	1.73E-3
	Order	-	2.05	2.18	2.09
m=3	Value	942.74	879.49	876.84	876.69
	Error	7.54E-2	3.21E-3	1.88E-4	1.07E-5
	Order	-	4.55	4.09	4.12
m=4	Value	897.55	876.85	876.68	876.68
	Error	2.38E-2	2.00E-4	2.91E-6	4.33E-8
	Order	-	6.89	6.10	6.07

Table 4. Table 4.13. The number j N subscript 𝑗 𝑁 j_{N} of linear converged eigenvalues.

Order	$N$(#DOF)	242	1,046	4,278
$m=1$	Laplace	8 (3.3%)	17 (1.6%)	39 (0.9%)
$m=2$		32 (13.2%)	92 (8.8%)	270 (6.3%)
$m=3$		38 (15.7%)	147 (14.0%)	553 (12.9%)
$m=4$		96 (39.6%)	355 (33.9%)	1417 (33.1%)
$m=2$	Biharmonic	24(9.9%)	53(5.0%)	94(2.2%)
$m=3$		45(18.6%)	204(19.5%)	691(16.1%)
$m=4$		170 (70.2%)	705 (67.3%)	2798 (65.4%)

Equations

$$U_h:=\{v\in L^2(\Omega )\mid v{|}_K\in {\mathbb{P}}^0(K),\forall K\in {\mathcal{T}}_h\}.$$

$${\mathcal{R}}_{K}v=\arg \underset{p\in {\mathbb{P}}^m(S(K))}{\min }\sum _{x\in {\mathcal{I}}_K}{\left|v(x)-p(x)\right|}^2.$$

$$(\mathcal{R}v){|}_K:=({\mathcal{R}}_{K}v){|}_K,\forall K\in {\mathcal{T}}_h.$$

$$V_h:=\mathcal{R}(U_h).$$

$${\psi }_K:=\mathcal{R}{e}_K,$$

$$e_K(x)=\{\begin{array}{l}1,x\in K,\\ 0,x\notin K.\end{array}$$

$$\mathcal{R}g=\sum _{K\in {\mathcal{T}}_h}g(x_K){\psi }_K(x),\forall g\in U_h.$$

$$S(K_0)=\left\{K_0,K_1,K_2,K_3,K_4\right\},$$

$${\mathcal{I}}_{K_0}=\left\{(x_{K_i},y_{K_i},z_{K_i}),i=0,1,2,3,4\right\}.$$

$${\mathcal{R}}_{K_0}g={\mathrm{argmin}}_{(a,b,c,d)\in \mathbb{R}}\sum _{i=0}^4|g(x_{K_i},y_{K_i},z_{K_i})-(a+b{x}_{K_i}+c{y}_{K_i}+d{z}_{K_i}){|}^2.$$

$$[a,b,c,d{]}^T=(A^{T}A{)}^{-1}{A}^{T}q,$$

$$A=\left[\begin{array}{cccc}1& x_{K_0}& y_{K_0}& z_{K_0}\\ 1& x_{K_1}& y_{K_1}& z_{K_1}\\ 1& x_{K_2}& y_{K_2}& z_{K_2}\\ 1& x_{K_3}& y_{K_3}& z_{K_3}\\ 1& x_{K_4}& y_{K_4}& z_{K_4}\end{array}\right],q=\left[\begin{array}{c}g(x_{K_0},y_{K_0},z_{K_0})\\ g(x_{K_1},y_{K_1},z_{K_1})\\ g(x_{K_2},y_{K_2},z_{K_2})\\ g(x_{K_3},y_{K_3},z_{K_3})\\ g(x_{K_4},y_{K_4},z_{K_4})\end{array}\right].$$

$$g{|}_{\mathcal{I}(K)}=0\text{implies}g{|}_{S(K)}\equiv 0.$$

$$\Lambda (m,{\mathcal{I}}_K):=\underset{p\in {\mathbb{P}}^m(S(K))}{\max }\frac{\Vert p{\Vert }_{L^{\infty }(S(K))}}{\Vert p{|}_{{\mathcal{I}}_K}{\Vert }_{{\ell }_{\infty }}}.$$

$${\mathcal{R}}_{K}g=g\text{for all}g\in {\mathbb{P}}^m(S(K)).$$

$$\Vert {\mathcal{R}}_{K}g{\Vert }_{L^{\infty }(K)}\le \Lambda (m,{\mathcal{I}}_K)\sqrt{\#{\mathcal{I}}_K}\Vert g{|}_{\mathcal{I}(K)}{\Vert }_{{\ell }_{\infty }},$$

$$\Vert g-{\mathcal{R}}_{K}g{\Vert }_{L^{\infty }(K)}\le {\Lambda }_m\underset{p\in {\mathbb{P}}^m(S(K))}{\inf }\Vert g-p{\Vert }_{L^{\infty }(S(K))},\forall K\in {\mathcal{T}}_h,$$

$$\Vert g-\mathcal{R}g{\Vert }_{L^2(K)}\le C{\Lambda }_m{h}_K{d}_K^m{\left|g\right|}_{H^{m+1}(K)},$$

$$\Vert \nabla (g-\mathcal{R}g){\Vert }_{L^2(K)}\le C\left(h_K^m+{\Lambda }_m{d}_K^m\right){\left|g\right|}_{H^{m+1}(K)}.$$

$$\left\{\begin{array}[]{ll}-\Delta u=\lambda u,&\text{in }\Omega,\\[6.45831pt] u=0,&\text{on }\partial\Omega,\end{array}\right.$$

$$a(u,v)=\lambda (u,v),\forall v\in V,$$

$$\left\{\begin{array}[]{ll}\Delta^{2}u=\lambda u,&\text{in }\ \Omega,\\[6.45831pt] u=\dfrac{\partial u}{\partial\boldsymbol{\mathrm{n}}}=0,&\text{on }\partial\Omega,\end{array}\right.$$

$$a(u,v)=\lambda (u,v),\forall v\in V,$$

$$a_h(\mathcal{R}{u}_h,\mathcal{R}{v}_h)={\lambda }_h(\mathcal{R}{u}_h,\mathcal{R}{v}_h),\forall v_h\in U_h.$$

$$a_h(v,w):=$$

$$-\sum _{e\in {\mathcal{E}}_h}{\int }_e\left([[\nabla v]]\left\{w\right\}+[[\nabla w]]\left\{v\right\}\right)\mathrm{d}s$$

$$+\sum _{e\in {\mathcal{E}}_h}{\int }_e{\eta }_e{h}_e^{-1}[[v]]\cdot [[w]]\mathrm{d}s,$$

$$a_h(v,w):=$$

$$+\sum _{e\in {\mathcal{E}}_h}{\int }_e\left([[v]]\left\{\nabla \Delta w\right\}+[[w]]\left\{\nabla \Delta v\right\}\right)\mathrm{d}s$$

$$-\sum _{e\in {\mathcal{E}}_h}{\int }_e\left(\left\{\Delta w\right\}[[\nabla v]]+\left\{\Delta v\right\}[[\nabla w]]\right)\mathrm{d}s$$

$$+\sum _{e\in {\mathcal{E}}_h}{\int }_e\left({\alpha }_e{h}_e^{-3}[[v]]\cdot [[w]]+{\beta }_e{h}_e^{-1}[[\nabla u]][[\nabla v]]\right)\mathrm{d}s,$$

$$\{q\}=\frac{1}{2}(q^{+}+q^{-}),\{\mathbf{v}\}=\frac{1}{2}({\mathbf{v}}^{+}+{\mathbf{v}}^{-}),$$

$$[[q]]={\mathrm{n}}^{+}{q}^{+}+{\mathrm{n}}^{-}{q}^{-},[[\mathbf{v}]]={\mathrm{n}}^{-}\cdot {\mathbf{v}}^{+}+{\mathrm{n}}^{-}\cdot {\mathbf{v}}^{-}.$$

Full text

keyword: elliptic eigenvalue problem, discontinuous Galerkin method, patch reconstruction

MSC2010: 49N45; 65N21

1. Introduction

In this paper, we consider the numerical method for solving eigenvalue problems of $2p$-th order elliptic operator for $p=1$ and $2$. Those problems arise in many important applications. The Laplace eigenvalue problem occurs naturally in vibrating elastic membranes, electromagnetic waveguides and acoustic theory, and the biharmonic eigenvalue problem appears in mechanics and inverse scatting theory.

The conforming finite element method (FEM) for eigenvalue problems has been well investigated. We refer to the review papers of Kuttler and Sigillito [27] and Boffi [8] for the details. For the biharmonic operator, we have the commonly used $C^{1}$ Argyris element [2] and the $C^{0}$ interior penalty Galerkin method ($C^{0}$ IPG) [17, 9, 11]. An old but hot topic for eigenvalue problems is the upper and lower bounds since [18]. It is well known that the conforming FEM can easily achieve the upper bound of the eigenvalues. In [3] and [24], the lower bound was achieved by mass lumping, see also other methods in [33, 8, 4]. Hu et al. [25, 22, 23] proposed a systematic method to produce lower bounds by nonconforming approximation spaces. The discontinuous Galerkin (DG) method, see for example [15, 5, 10], has been applied to the Laplace eigenvalue problem [1] and the Maxwell eigenvalue problem [21, 36]. As a nonconforming approximation, the DG method admits the totally discontinuous polynomial space which leads to a great flexibility though it is challenged [26] on its efficiency in number of degrees of freedom (DOF).

In a recent work [37], Zhang studied an interesting issue on the number of ”trusted” eigenvalues by finite element approximation for the elliptic eigenvalue problems. It was pointed out therein that only eigenvalues lower in the spectrum can achieve optimal convergence rate. Furthermore, the percentage of reliable eigenvalues will decrease on a finer mesh even if we relax the convergence rate to linear. Typically, the optimal convergence rate of the elliptic eigenvalue problem is $h^{2(m+1-p)}$, where $m$ is the polynomial degree. It is implied that high order methods are more likely to provide a greater number of reliable eigenvalues, measured relatively to the DOFs used, than a lower order method.

Motivated by Zhang’s result, in this paper we aim to apply a symmetric interior penalty discontinuous Galerkin method to elliptic eigenvalue problems. The method adopts a discontinuous approximation space proposed in [28], where it was applied to solve elliptic boundary value problems. The core of the method is to construct an approximation space by the patch reconstruction technique in a way that one DOF is used in each element. The reconstructed space is a piecewise polynomial space and is discontinuous across the element face, thus it is a subspace of the traditional DG space. The idea has been applied smoothly to the biharmonic equation [29] and the Stokes equation [31, 32]. For elliptic eigenvalue problems, it is a direct extension of the method for boundary value problems. Consequently, the optimal error estimates of the approximation to eigenfunctions and eigenvalues can be obtained instantly from existing results for arbitrary order accuracy.

We present all details on the numerical results to verify that higher order methods can provide much more reliable eigenvalues, which perfectly agrees with the theoretical prediction in [37]. In comparison to the classical DG method, one may see that the patch reconstructed approximation space attains the same accuracy with much less degrees of freedom. In case of using higher order polynomials, the numerical results show that a better efficiency in number of DOFs can be achieved by our method even than conforming finite element methods. We note that for problems with high regularities, the conforming finite element methods traditionally outperform the other methods in number of DOFs. The new observation here in efficiency gives us an enthusiastic encouragement to apply our method with high order polynomials to elliptic eigenvalue problems.

The rest of this paper is organized as follows. To be self-contained, we describe in section 2 the detailed process to construct the approximation space and the approximation properties of the corresponding space. The symmetric interior penalty method for elliptic operators is presented in section 3, and the optimal error estimates are then given for the eigenvalues and eigenfunctions. In section 4, we present the numerical results to illustrate that the proposed method is efficient for elliptic eigenvalue problems.

2. Approximation Space

Let us consider a convex polygonal domain $\Omega$ in $\mathbb{R}^{D}$, $D=2,3$. $\mathcal{T}_{h}$ is a polygonal partition of the domain $\Omega$. For each polygon $K$, $h_{K}$ and $\left\lvert K\right\rvert$ denote its diameter and area, respectively. Besides, let $h{:}=\max_{K\in\mathcal{T}_{h}}h_{K}$. For the optimal convergence analysis, the partition $\mathcal{T}_{h}$ is assumed to satisfy some shape regularity conditions. Those regularity conditions are commonly used in mimetic finite difference schemes [12, 7, 13] and discontinuous Galerkin method [34], which are stated as follows:

1. A1 

Any element $K\in\mathcal{T}_{h}$ admits a sub-decomposition $\widetilde{\mathcal{T}_{h}}|_{K}$ that consists of at most $N_{s}$ triangles, where $N_{s}$ is an integer independent of $h$; 2. A2 

If all the triangles $T\in\widetilde{\mathcal{T}_{h}}$ are shape-regular in the sense of Ciarlet-Raviart [14]: there exists a real positive number $\sigma$ independent of $h$ such that $h_{T}/\rho_{T}\leq\sigma$, where $\rho_{T}$ is the radius of the largest ball inscribed in $T$. Then the $\widetilde{\mathcal{T}_{h}}$ is a compatible sub-decomposition.

The above regularity assumptions lead to some useful estimates, such as Agmon inequality, approximation property and inverse inequality. Those inequalities are the foundations to derive the approximation error estimates for the finite element method. We refer to [28] for the detailed discussion.

The reconstruction operator $\mathcal{R}$ can be constructed with the given partition $\mathcal{T}_{h}$. The degrees of freedom of $\mathcal{R}$ are located at one point $x_{K}\in K$ on each element which are called the sampling nodes or collocation points. We usually assign the barycenter of $K$ as the sampling node $x_{K}$. Furthermore, the reconstruction operator $\mathcal{R}$ is defined element-wise. An element patch denoted as $S(K)$ is constructed for each element $K$. $S(K)$ is an agglomeration of elements including $K$ itself and other elements nearby $K$. Let $\mathcal{I}_{K}$ denote the set of sampling nodes belonging to $S(K)$, $\#S(K)$ and $\#\mathcal{I}_{K}$ denote the number of elements belonging to $S(K)$ and the number of sampling nodes belonging to $\mathcal{I}_{K}$, respectively. Obviously, these two numbers are equal to each other. We define $d_{K}{:}=\text{diam}\;S(K)$ and $d{:}=\max_{K\in\mathcal{T}_{h}}d_{K}$.

Here we specify the way to construct the element patch while it can be quite flexible, see [30, 28] for the alternative approaches. First, a constant number $t$ is assigned to $\#S(K)$ which is determined by the degree of polynomials. Then we initialize $S(K)$ as $\{K\}$, and fill $S(K)$ by adding the nearest Von Neumann neighbor (adjacent edge-neighboring elements) of the current geometry $S(K)$. We terminate the recursive process until the number $\#S(K)$ reaches the number $t$. With such an approach, the element patches are obtained with a constant number, which is convenient for the implementation. Meanwhile, the shape regularity of the geometry of $S(K)$ preserves. All the sampling nodes $x_{K}$ are located in element $K$ and all element patches are connected set, that the stability of reconstruction is fair promising. The reconstruction process can be conducted element-wise after the sampling nodes $\mathcal{I}_{K}$ and element patch $S(K)$ are specified.

Let $U_{h}$ be the piecewise constant space associated with $\mathcal{T}_{h}$, i.e.,

[TABLE]

For a piecewise constant function $v\in U_{h}$ and an element $K$, a high-order approximation polynomial $\mathcal{R}_{K}v$ of degree $m$ can be obtained by solving the following discrete local least-squares:

[TABLE]

We assume the problem (2.1) has a unique solution [28]. Now, we concentrate on the reconstruction operator and the corresponding finite element space. Although $\mathcal{R}_{K}v$ gives an approximation polynomial on element patch $S(K)$, we only use it on element $K$. The global reconstruction operator $\mathcal{R}$ is defined as:

[TABLE]

The reconstruction operator $\mathcal{R}$ actually defines a linear operator which maps $U_{h}$ into a piecewise polynomial space, denoted by

[TABLE]

Here, $V_{h}$ is the reconstructed finite element space which is spanned by the basis functions $\{\psi_{K}\}$. Here the basis functions are defined by the reconstruction operator,

[TABLE]

where $e_{K}\in U_{h}$ is the characteristic function corresponding to $K$,

[TABLE]

Thereafter the reconstruction operator can be explicitly expressed

[TABLE]

We present a 3D example below to illustrate the implementation of reconstruction process, while the details for 1D implementation and 2D implementation can be found in [29] and [31], respectively. We consider a linear reconstruction on a cubic domain $[0,1]^{3}$. The domain is partitioned into quasi-uniform tetrahedron elements using Gmsh [19], which is shown in Figure 2.1. We take element $K_{0}$ as an instance (see Figure 2.1). The number of degrees of freedom demanded by linear reconstruction is $4$. Therefore, the $\#S(K_{0})$ could be taken as $5$. In this case, the element patch is containing the element itself and 4 Von Neumann neighbors coincidentally. Figure 2.2 shows the geometry of the element patch and the corresponding sampling nodes. The element patch $S(K_{0})$ is chosen as

[TABLE]

and the sampling nodes are as follows,

[TABLE]

For any continuous function $g$, we consider the linear approximation for an illustration. For the polynomial degree $m=1$, the least squares problem (2.1) is specified as

[TABLE]

The solution of the problem is given by the generalized inverse of matrix,

[TABLE]

where $A$ and $q$ are

[TABLE]

A direct observation is that matrix $(A^{T}A)^{-1}A^{T}$ is not relevant to the interpolation function $g$. Moreover, the matrix $(A^{T}A)^{-1}A^{T}$ actually stores the polynomial basis function coefficients corresponding to $\psi_{K_{i}}$, $i=0,\cdots,4$. All the basis functions and the finite element space $V_{h}$ are determined after the reconstruction process on each element $\forall K\in\mathcal{T}_{h}$. Clearly, the basis functions are discontinuous across the interface.

Next, for completeness, we report the results on the properties of the reconstruction operator. Following [30], we first make the following assumption.

Assumption A  For any $K\in\mathcal{T}_{h}$ and $g\in\mathbb{P}^{m}(S(K))$,

[TABLE]

This assumption implies the uniqueness for least squares problem (2.1). A necessary condition for Assumption A is that the number $\#\mathcal{I}_{K}$ needs to be greater than $\text{dim}(\mathbb{P}^{m})$, whose quantities are $m+1$, $(m+1)(m+2)/2$ and $(3m^{2}+3m+2)/2$ corresponding to 1D,2D and 3D, respectively. A constant $\Lambda(m,\mathcal{I}_{K})$ is defined as [30]:

[TABLE]

Then, the uniform upper bound can be obtained by adding some constrains on element patches and the partition, see also [28] for the details. We have the following properties of the reconstruction operator $\mathcal{R}_{K}$.

Lemma 2.1**.**

[30, Theorem 3.3]** If Assumption A holds, then there exists a unique solution to (2.1). Moreover $\mathcal{R}_{K}$ satisfies

[TABLE]

The stability property holds true for any $K\in\mathcal{T}_{h}$ and $g\in C^{0}(S(K))$ as

[TABLE]

and the quasi-optimal approximation property is valid in the sense

[TABLE]

where $\Lambda_{m}{:}=\max_{K\in\mathcal{T}_{h}}\{1+\Lambda(m,\mathcal{I}_{K})\sqrt{\#\mathcal{I}_{K}}\}$.

With Lemma 2.1 and the interpolation result in [16], the local estimates on element $K$ can be obtained.

Lemma 2.2**.**

[28, Lemma 2.4]** Let $u\in C^{0}\left(\Omega\right)\cap H^{m+1}(\Omega)$, then there exists a constant $C$ that depends on $N_{s}$ and $\sigma$, but independent of $h$, such that

[TABLE]

and

[TABLE]

3. Elliptic Eigenvalue Problems

Let us consider the $2p$-th $(p=1,2)$ order elliptic eigenvalue problems, for $p=1$, the second order elliptic eigenvalue problem reads:

[TABLE]

and the corresponding weak form is: find $\lambda\in\mathbb{R}$ and $u\in V=H^{1}_{0}(\Omega)$, with $u\neq 0$, such that

[TABLE]

where $a(u,v)=\int_{\Omega}\nabla u\cdot\nabla v\mathrm{d}x$ and $(u,v):=\int_{\Omega}uv\mathrm{d}x$.

For $p=2$, the biharmonic eigenvalue problem reads:

[TABLE]

and the corresponding weak form is: find $\lambda\in\mathbb{R}$ and $u\in V=H^{2}_{0}(\Omega)$, with $u\neq 0$, such that

[TABLE]

where $a(u,v)=\int_{\Omega}\Delta u\Delta v\mathrm{d}x$.

The discretized variational problem for equations (3.1) and (3.2) reads: find $\lambda_{h}\in\mathbb{R}$ and $u_{h}\in U_{h}$, with $u_{h}\neq 0$, such that

[TABLE]

Here we use the notations $a$, $a_{h}$ for unification. In the rest of the paper, we will specify the sense of the notation when a particular equation is considered.

The symmetric interior penalty method is employed to discretize the elliptic operators. For the second order elliptic operator, $a_{h}(\cdot,\cdot)$ is

[TABLE]

and for the biharmonic operator, $a_{h}(\cdot,\cdot)$ is

[TABLE]

where $\eta_{e},\alpha_{e},\beta_{e}$ are positive constants. Here we let $\mathcal{E}_{h}$ denote the collection of all the faces of $\mathcal{T}_{h}$, $\mathcal{E}_{h}^{i}$ denote the collection of the interior faces. The set of boundary faces is denoted as $\mathcal{E}_{h}^{b}$, and then $\mathcal{E}_{h}=\mathcal{E}_{h}^{i}\cup\mathcal{E}_{h}^{b}$. Let $e$ be an interior face shared by two neighbouring elements $K^{+},K^{-}$, and $\boldsymbol{\mathrm{n}}^{+}$ and $\boldsymbol{\mathrm{n}}^{-}$ denote the corresponding outward unit normal. For the scalar-valued function $q$ and the vector-valued function $\boldsymbol{v}$, the average operator $\left\{\,\cdot\,\right\}$ and the jump operator $[\![\cdot]\!]$ are defined as

[TABLE]

and

[TABLE]

Here $q^{+}=q|_{K^{+}}$, $\boldsymbol{v}^{+}=\boldsymbol{v}|_{K^{+}}$ and $q^{-}=q|_{K^{-}}$, $\boldsymbol{v}^{-}=\boldsymbol{v}|_{K^{-}}$. For $e\in\mathcal{E}_{h}^{b}$, we set

[TABLE]

and

[TABLE]

We note that the problem (3.3) is equivalent to the following problem: find $\lambda_{h}\in\mathbb{R}$ and $\varphi_{h}\in V_{h}$, with $\varphi_{h}\neq 0$, such that

[TABLE]

This is a more standard formulation for finite element methods. By the formulation (3.3), it is emphasized that the number of DOFs of the approximation space is always $\text{dim}(U_{h})$.

We define the energy norms $\|\cdot\|_{h}$ and $\|\!|\,\cdot\,\|\!|_{h}$ for any $v\in V_{h}=\mathcal{R}(U_{h})$ as:

[TABLE]

From the Lemma 2.2 and Agmon inequality, the following interpolation estimates are straightforward results for the reconstruction operator in the energy norm.

Lemma 3.1**.**

[28, Equation 3.4]** [29, Theorem 2.1] Let $u\in H^{m+1}(\Omega)$, and $\mathcal{R}u\in V_{h}$ be the interpolation polynomial of $u$, there exists a constant $C$ that depends on $N_{s}$, $\sigma$ and $m$, but independent of $h$, such that

[TABLE]

Next, the boundedness and coercivity of the bilinear operator $a_{h}(\cdot,\cdot)$ in (3.4) and (3.5) are as below.

Lemma 3.2**.**

[5, Equations 4.4,4.10]** If the penalty constant $\eta_{e}$ is sufficiently large, then the bilinear operator (3.4) is bounded and coercive, indeed there exist constants $C_{b}$ and $C_{s}$, such that

[TABLE]

[29, Lemmata 3.1, 3.2]** If the penalty constants $\alpha_{e},\beta_{e}$ are sufficiently large , then there exist constants $C_{b}$ and $C_{s}$, such that the bilinear operator (3.5) satisfies

[TABLE]

We refer to [5, 29] for the proof.

To derive the error estimates, we introduce the sum space $V(h)=V+\mathcal{R}(U_{h})$, and endows it with the energy norm (3.6), denoted as $\|\cdot\|_{V(h)}$ for unification,

[TABLE]

Let $\lambda^{(i)}$, $i\in\mathbb{N}$, denote the sequence of eigenvalues of (3.1) and (3.2) with the natural numbering

[TABLE]

and the corresponding eigenfunctions with the standard normalization $\|u^{(i)}\|=1$

[TABLE]

which are orthogonal to each other

[TABLE]

Let $N=\text{dim}(V_{h})$, thus the discrete eigenvalues of (3.3) can be ordered as follows:

[TABLE]

and the discrete eigenfunctions with the normalization $\|\mathcal{R}u^{(i)}_{h}\|=1$,

[TABLE]

which satisfy the same orthogonalities

[TABLE]

The convergence analysis for the eigenvalue problem (3.3) can be obtained by the Babuška-Osborn theory [6]. We define the following continuous and discrete solution operators:

[TABLE]

Obviously the operator $T$ and $T_{h}$ are self-adjoint and from the elliptic regularity, there exists $\epsilon>0$ such that

[TABLE]

And the operators have the gradual approximation property,

[TABLE]

Let $\sigma(T),\sigma(T_{h})$ and $\rho(T),\rho(T_{h})$ denote the spectrum and the resolvent set of the solution operator $T$ and $T_{h}$, respectively. Define the resolvent operators as follows

[TABLE]

Then the first result of convergence is that there is no pollution of the spectrum.

Theorem 3.3**.**

[8, Theorem 9.1]** Assume the convergence in norm (3.11) is satisfied, for any compact set $K\subset\rho(T)$, there exists $h_{0}>0$, such that, for all $h<h_{0}$, we have

[TABLE]

If $\mu\in\sigma(T)$ is a non-zero eigenvalue with algebraic multiplicity $k$, then exactly $k$ discrete eigenvalues of $T_{h}$, convergence to $\mu$ as $h$ tend to zero.

Let $\Gamma$ be an arbitrary closed smooth curve $\Gamma\in\rho(T)$ which encloses $\mu\in\sigma(T)$, and no other elements of $\sigma(T)$, we define the Riesz spectral projection operators $E$, $E_{h}$ by:

[TABLE]

When $h$ is sufficiently small, we have $\Gamma\in\rho(T_{h})$ and $\Gamma$ encloses exactly $k$ eigenvalues of $T_{h}$. More precisely, the dimension of $E(\mu)V$ and $E_{h}(\mu)\mathcal{R}(U_{h})$ is equal to $k$. Further we have

[TABLE]

The convergence of the generalized eigenvectors has been claimed.

The gap between the eigenspaces is defined as follows,

[TABLE]

Lemma 3.4**.**

[8, Theorem 9.3]** Let $\mu$ be a non-zero eigenvalue of $T$, let $E=E(\mu)V$ be its generalized eigenspace, and let $E_{h}=E_{h}(\mu)\mathcal{R}(U_{h})$. Then

[TABLE]

Lemma 3.5**.**

[8, Corollary 9.4]** Let $\lambda$ be a non-zero eigenvalue of (3.1) and (3.2), respectively, and $E=E(\lambda^{-1})V$ be its generalized eigenspace, and let $E_{h}=E_{h}(\lambda^{-1})\mathcal{R}(U_{h})$. Then

[TABLE]

We now claim the approximation estimate for the solution operator, and we refer to [28, 29] for more details.

Lemma 3.6**.**

Let $\lambda$ be a non-zero eigenvalue of (3.1) and (3.2), respectively, let $E$ be the eigenspace associated with $\lambda$, and its regularity satisfy $E\subset H^{m+1}(\Omega)$, $m\geq 2p-1$, then

[TABLE]

where $\tau=m+1-p$.

Proof.

The source problem corresponding to (3.1) takes the form

[TABLE]

and the source problem corresponding to (3.2) takes the form

[TABLE]

The discrete variational problem for the source problem reads: find $u_{h}\in U_{h}$ such that

[TABLE]

From [28, Theorem 3.1] [29, Theorem 3.1], we conclude that there exists a unique solution to (3.13). Furthermore, if $u_{s}\in H^{m+1}(\Omega)$, we have the following estimate:

[TABLE]

where $\tau=m+1-p$. This estimate directly implies

[TABLE]

which completes the proof. ∎

Then, the error estimates for the eigenfunctions can be directly derived.

Theorem 3.7**.**

Let $u^{(i)}$ be a unit eigenfunction associated with an eigenvalue $\lambda^{(i)}$ of multiplicity $k$, such that $\lambda^{(i)}=\cdots=\lambda^{(i+k-1)}$, and $\mathcal{R}u^{(i)}_{h},\cdots,\mathcal{R}u^{(i+k-1)}_{h}$ denote the discrete eigenfunctions associated with the $k$ discrete eigenvalues converging to $\lambda^{(i)}$. Then there exists

[TABLE]

such that

[TABLE]

Moreover, if the regularity of eigenspace is $E\subset H^{m+1}(\Omega)$, $m\geq 2p-1$, then

[TABLE]

where $\tau=m+1-p$.

Proof.

The results (3.14) and (3.15) are direct extensions of Lemma 3.5 and the estimate  (3.16) is directly derived from Lemma 3.6. ∎

Finally, the error estimates for the eigenvalues of (3.1) and (3.2) are the following.

Theorem 3.8**.**

Let $\lambda^{(i)}$ denote the eigenvalue of (3.1) and (3.2) with multiplicity $k$, $\lambda^{(i)}_{h}$ be the discrete eigenvalues and $E$ denote the eigenspace associated with $\lambda^{(i)}$, then we have

[TABLE]

Moreover, if eigenspace $E\subset H^{m+1}(\Omega)$, $m\geq 2p-1$, then the following optimal double order of convergence holds

[TABLE]

where $\tau=m+1-p$.

Proof.

Since the operator $a_{h}(\cdot,\cdot)$ is symmetric, i.e. $a_{h}(T_{h}f,\mathcal{R}v)=a_{h}(\mathcal{R}v,T_{h}f)$, the estimate (3.17) is a direct application of  [8, Theorem 9.13]. The estimate (3.18) is the combination of the inequality (3.17) and Lemma 3.6. ∎

4. NumericalResults

In this section, we present some numerical results to show that our method is efficient for eigenvalue problems if we use higher order approximation. We would like to emphasize two points:

• •

Less DOFs are used by our method for high order approximation comparing to the classical DG method and conforming finite element methods;

• •

More reliable eigenvalues can be obtained increasing the order of approximation.

Besides, we will compute the numerical order of convergence to verify the theoretical error estimates and give results on different domains and different meshes to demonstrate the flexibility of the implementation using our method.

4.1. Examples setup

First, let us list the setup of the examples to be investigated.

Example 1

We consider the two-dimensional square domain $\Omega=[0,\pi]^{2}$, the eigenpairs of problem (3.1) are given by

[TABLE]

and for the problem (3.2) with the boundary condition $u|_{\partial\Omega}=\Delta u|_{\partial\Omega}=0$, which is related to the bending of a simply supported plate [11], the eigenpairs are given by

[TABLE]

In this example, the computation involves a series of regular unstructured triangular meshes which are generated by Gmsh [19]. For the second order elliptic problem we take $(m=1,2,3,4,5)$ and for the biharmonic problem $m$ is taken as $(2,3,4,5)$.

Example 2

We consider the L-shaped domain $[-1,1]^{2}\backslash(0,1]\times(0,-1]$. The domain is partitioned into polygonal meshes by PolyMesher [35]. Figure 4.1 shows the initial mesh and the refined mesh. The meshes contain the elements with various geometries such as quadrilaterals, pentagons, hexagons, and so on. The first eigenfunction in L-sharped domain has a singularity at the reentrant corner and has no analytical expression. We note that the third eigenpair is smooth for L-shaped domain. For the second order elliptic equation, the third eigenvalue is $2\pi^{2}$ and the corresponding eigenfunction is $\sin(\pi x)\sin(\pi y)$, and we take $(m=1,2,3)$ to solve the eigenvalue problem. For the biharmonic equation, the third eigenpair is $4\pi^{4}$ and $\sin(\pi x)\sin(\pi y)$, and we choose $(m=2,3)$ to solve it.

Example 3

We solve the eigenvalue problem in three dimensions in this example. The computational domain is the unit cubic $\Omega=[0,1]^{3}$ which is partitioned into tetrahedral meshes by Gmsh. The eigenpairs of problem (3.1) are as follows:

[TABLE]

and for problem (3.2) with the simply supported plate boundary condition, the eigenpairs are given by

[TABLE]

4.2. Convergence order study

At first, we show that the numerical results verify the optimal convergence order as predicted by the theory.

For Example 1, Figure 4.2 shows the convergence rates of the eigenvalue and eigenfunction. The exact $20$-th eigenvalue is $32$ and the corresponding eigenfunction is $\sin(4x)\sin(4y)$. The eigenvalue converges to the exact one with $h^{2m}$ rate and for the eigenfunction the convergence rate is $h^{m}$. Figure 4.3 shows the convergence rates of the eigenvalue and eigenfunction of the biharmonic equation. The exact $20$-th eigenvalue is $1024$ and the corresponding eigenfunction is $\sin(4x)\sin(4y)$. The eigenvalue convergences to the exact one with $h^{2(m-1)}$ rate, and for the eigenfunction the convergence rate is $h^{m-1}$. The numerical results agree with Theorem 3.7 and 3.8 perfectly.

The Example 2 shows that the proposed method can handle these polygonal elements easily. First, we calculate the third smooth eigenpair to verify the analysis of the proposed method. Figures 4.4 and 4.6 show the numerical results that agree with the theoretical prediction. The values of the first eigenvalue of the Laplace/biharmonic equation are shown in Table 4.8. It is clear that the eigenvalues converge to the real eigenvalue as $h$ approaches [math]. The eigenfunctions corresponding to the first eigenvalue and third eigenvalue are presented in Figures 4.5 and 4.7.

For Example 3, the numerical results are presented in Tables  4.9 and  4.10 for the second order and biharmonic equation, respectively. The convergence order of the second order equation is $h^{2m}$, and of the biharmonic equation is $h^{2(m-1)}$. Obviously, the computational results agree with the error estimates.

Remark 4.1**.**

We note that all the eigenvalues obtained by the proposed method are greater than the exact eigenvalues. This behavior appears if conforming finite element method is used to solve the eigenvalue problem. However, the approximate space $V_{h}$ is not a subspace of the space $V=H^{1}_{0}$ or $H^{2}_{0}$. In DG framework, this phenomenon is related to the penalty parameter. Warburton and Embree studied the role of penalty in the LDG method for Maxwell’s eigenvalue problem in [36]. Giani et al. [20] used the asymptotic perturbation theory to analyze the dependence of eigenvalues and eigenspaces on the penalty parameter. We hope the reason why this happened in our method can be clarified in future study.

4.3. Efficiency in terms of number of DOFs

Next, we make a comparison in terms of number of DOFs among different methods. For the second order elliptic problem, we consider the conforming FEM, standard SIPDG method  [1] and our method. For the biharmonic problem, we consider the $C^{0}$ IPG, standard SIPDG and our method. Here we will study the numerical behavior for higher order approximation. We restrict to Example 1, since in this case the solution has enough regularity.

We calculate the first eigenvalue and eigenfunction on successively refined meshes. The errors of eigenvalue are measured in the relative error, and the errors of the eigenfunction are measured in $|\cdot|_{1,h}$ and $|\cdot|_{2,h}$ semi-norms, respectively.

For the Laplace problem, Figure 4.11 shows the performance of the conforming FEM, SIPDG method and our method. The approximation order $m$ is taken from $1$ to $4$. The convergence rate for the eigenvalue is $h^{2m}$ and for the eigenfunction the rate is $h^{m}$ which meet the theoretical predictions. The horizontal ordinate is the number of DOFs. The number of DOFs employed by our method is fixed while the approximation order increases. In all cases, the SIPDG method uses the maximum number of DOFs. As one’s expectation, the figure shows that the efficiency of FEM is higher than others for the low order approximation. Increasing of the approximation order, our method becomes the most efficient method among these three methods.

For the biharmonic problem, Figure 4.12 shows the error in terms of number of DOFs of the $C^{0}$ IPG method, standard SIPDG method and our method. The approximation order $m$ is taken as $2$, $3$, and $4$. The convergence rate for eigenvalue is $h^{2(m-1)}$ and the convergence rate is $h^{m-1}$ for eigenfunction which perfectly agree with the error estimates. The experiments show that our method performs better than the other methods in all cases. The advantage of our method in efficiency is more remarkable for higher order approximation.

4.4. Number of reliable eigenvalues

Zhang studied the number of reliable eigenvalues of the finite element method in [37], and the main result he gave is as below:

Theorem 4.1**.**

Suppose that we solve a $2p$-order elliptic equation on a domain $\Omega\in\mathbb{R}^{D}$ by the finite element method (conforming or non-conforming) of polynomial degree $m$ under a shape regular and quasi-uniform mesh with mesh-parameter $h$. Assume that the exact eigenvalue grows as $\lambda_{j}=O(j^{\frac{2p}{D}})$ and the relative error can be estimated by $\frac{\lambda_{i}^{h}-\lambda_{i}}{\lambda_{i}}=h^{m+1-p}\lambda_{i}^{\frac{m+1}{p-1}}$. Then there are about

[TABLE]

reliable numerical eigenvalues with the relative error of $\lambda_{jN}$, converging at rate $h^{\alpha}$ for $\alpha\in(0,2(m+1-p)]$. Here $N$ is the total degrees of freedom.

Theorem 4.1 implies that the quantity of the reliable numerical eigenvalues who have the optimal convergence rate $\alpha=2(m+1-p)$ is $O(1)$, which means only eigenvalues lower in the spectrum can achieve the optimal convergence rate. Therefore, for the eigenvalue problem, the number of eigenvalues that have the optimal convergence rate is very small. We here relax the convergence rate to linear, saying taking $\alpha=1$, to identify if a numerical eigenvalue is reliable. For the lowest order approximation of the eigenvalue problem, linear element for Laplace operator and quadratic element for biharmonic operator shall be involved. The predicted number of the reliable numerical eigenvalues from Theorem 4.1 is $O(N^{1/2})$, which implies that the percentage of the reliable numerical eigenvalues reduce rapidly as the number of DOFs of the system increases. For the higher order approximation, the percentage of the reliable numerical eigenvalues reduces much slower than the low order approximation.

To identify numerically if an eigenvalue is reliable, we define the relative error by $\frac{|\lambda-\lambda_{h}|}{|\lambda|}$, and the convergence rate by $\log_{2}\left(\frac{|\lambda-\lambda_{2h}|}{|\lambda-\lambda_{h}|}\right)$. If the convergence rate is not less than $1$, the eigenvalue is identified as reliable. We carry out a series of numerical experiments with various $m$, while the results are quite robust with almost the same efficiency.

Again we are limited to study the setup in Example 1 since we need reference solutions. We calculate $j_{N}$ eigenvalues whose relative errors are of order $O(h)$. Precisely, we enumerate the number of the eigenvalues that are at least linearly convergent, with the result given in Table 4.13. For the Laplace problem, there are $O(N^{1/2})$ reliable numerical eigenvalues. In this table, the percentage decreases rapidly as the computational scale $N$ increases. The number of the eigenvalues that are at least linearly convergent increases a lot if the higher order approximation is applied, which is as implied by Zhang’s result that the higher order method could produce more reliable numerical eigenvalues with the same $N$. Moreover, for the higher order method, the percentage of the reliable numerical eigenvalues reduces much slower than the lower order method.

The behavior of the number of reliable eigenvalues is similar for the biharmonic equation, as shown in Table 4.13. The numerical results confirm the prediction of Theorem 4.1 and emphasize that the higher order approximations are more robust and preferred for the eigenvalue problem.

5. Conclusion

We applied the symmetric interior penalty discontinuous Galerkin method based on a patch reconstructed approximation space for solving elliptic eigenvalue problem. The proposed method, when compared to other existing approximation methods, can be implemented in a more flexible way and its approximation properties are easier to analyse. Numerical results confirm the optimal convergence rates and emphasize the great efficiency of our method in number of DOFs. The great efficiency and convenient implementation is even remarkable in the case of higher order approximation. Since high order approximation is preferred for the elliptic eigenvalue problems, our method is a quite appropriate method to solve the elliptic eigenvalue problems.

Acknowledgment

The authors would like to thank the anonymous referees. They have very constructively helped to improve the original version of this paper.

The research is supported by the National Natural Science Foundation of China (Grant No. 91630310, 11421110001 and 11421101) and Science Challenge Project, No.TZ2016002.