Multigrid Discretization and Iterative Algorithm for Mixed Variational Formulation of the Eigenvalue Problem of Electric Field

Abstract

This paper discusses highly finite element algorithms for the eigenvalue problem of electric field. Combining the mixed finite element method with the Rayleigh quotient iteration method, a new multi-grid discretization scheme and an adaptive algorithm are proposed and applied to the eigenvalue problem of electric field. Theoretical analysis and numerical results show that the computational schemes established in the paper have high efficiency.

1. Introduction

The finite element method for eigenvalue problem of electric field has become a hot topic in the field of mathematics and physics (see, e.g., [1–7]). This paper discusses high efficient mixed finite element calculation schemes for the eigenvalue problem of electric field.

Kikuchi [6] introduced the first type of mixed variational formulation for the eigenvalue problem of electric field. Based on this formulation, in [3] Buffa et al.analyzed the approximation of nodal finite element. Boffi et al. [1] discussed the second type of mixed variational formulation for the eigenvalue problem of electric field and analyzed approximations of edge element and nodal element. Yang et al. [7] studied a two-grid discretization scheme of finite element for the first type of mixed variational formulation.

Based on the work mentioned above, in this paper a new multi-grid discretization scheme and an adaptive algorithm are proposed for the first type of mixed variational formulation of eigenvalue problem and applied to the eigenvalue problem of electric field. The main features of the research in this paper are as follows.(1)Our multi-grid discretization scheme and adaptive algorithm, which are the extension of conforming finite element multi-grid discretization scheme (see scheme  3 in [8] and scheme  1 in [9]), are a combination of the mixed finite element method and the Rayleigh quotient iteration method (see the algorithm  27.3 in [10]). With our algorithm one solves an eigenvalue problem on a coarse grid just at the first step and always solves a linear algebraic system on finer and finer grids at each following step. We derive the error estimates for the algorithm and prove that the constants appeared in the error estimates are independent of the iteration degrees. Thus we prove the convergence of iterations.(2)The eigenvalue problem of electric field is so complicated that it is very difficult to obtain local a posteriori error estimates for the eigenfunctions of mixed finite element. As yet, there is no achievement reported in this field. Our adaptive algorithm substitutes the weight method established by Costabel and Dauge (see [3, 11]) for local refinement, which uses 𝜃×|𝜆^ℎ𝑙−𝜆^ℎ𝑙−1|(𝜃∈(0,1]) as a posteriori error estimator of 𝜆^ℎ𝑙 instead of estimating local a posteriori error for the eigenfunction. And the results are satisfying. (3)We analyze the mixed finite element error for the eigenvalue problem of electric field (see Theorem 2.2 and Theorem 4.2). We refer to [12] to propose a new proof method which differs from the usual one in [13].

The rest of this paper is organized as follows. Some preliminaries of finite element approximations for eigenvalue problems which are needed in this paper are provided in the next section. In Section 3, for the first type of mixed variational formulation of eigenvalue problem, the finite element multi-grid discretization scheme and the adaptive algorithm are established and the validity of these schemes is proved theoretically. In Section 4, the multi-grid discretization scheme is applied to the eigenvalue problem of electric field. Finally, numerical experiments are presented in Section 5.

2. Preliminaries

Let 𝑉, 𝑊, and 𝐷 be three real Hilbert spaces with inner products and norms (⋅,⋅)_𝑉, ‖⋅‖_𝑉, (⋅,⋅)_𝑊, ‖⋅‖_𝑊, (⋅,⋅)_𝐷, and ‖⋅‖_𝐷, respectively. Suppose that 𝑉↪𝐷 (continuously imbedded), 𝑎(⋅,⋅) is a symmetric, continuous, and 𝑉-elliptic bilinear form on 𝑉×𝑉, that is, 𝑎𝑞,𝜓≤𝑀₁𝑞_𝑉𝜓_𝑉𝑎𝑞,∀𝑞,𝜓∈𝑉,𝑞,𝑞≥𝜈²_𝑉,∀0≠𝑞∈𝑉. (2.1)𝑏(⋅,⋅) is a continuous bilinear form on 𝑉×𝑊, that is, 𝑏𝜓,𝑣≤𝑀₂𝜓_𝑉‖𝑣‖_𝑊,∀𝜓∈𝑉,𝑣∈𝑊. (2.2) It is obvious that 𝑎(⋅,⋅) is an inner product on 𝑉 and ‖⋅‖_𝑎=𝑎(⋅,⋅) is an equivalent norm to ‖⋅‖_𝑉.

In scientific and engineering computations, many eigenvalue problems have the following first type of mixed variational formulation: find (𝜆,𝑢,𝜎)∈𝑅×𝑉×𝑊, (𝑢,𝜎)≠(0,0), such that 𝑎𝑢,𝜓+𝑏𝜓,𝜎=𝜆𝑢,𝜓_𝐷𝑏(),∀𝜓∈𝑉,𝑢,𝑣=0,∀𝑣∈𝑊. (2.3)(2.4)

In order to solve problem (2.3)-(2.4), one should construct finite element spaces 𝑉_ℎ⊂𝑉 and 𝑊_ℎ⊂𝑊. Restricting (2.3)-(2.4) on 𝑉_ℎ×𝑊_ℎ, we get the conforming mixed finite element approximation as follows: find (𝜆_ℎ,𝑢_ℎ,𝜎_ℎ)∈𝑅×𝑉_ℎ×𝑊_ℎ,(𝑢_ℎ,𝜎_ℎ)≠(0,0), such that 𝑎𝑢_ℎ,𝜓+𝑏𝜓,𝜎_ℎ=𝜆_ℎ𝑢_ℎ,𝜓_𝐷,∀𝜓∈𝑉_ℎ,𝑏𝑢_ℎ,𝑣=0,∀𝑣∈𝑊_ℎ. (2.5)(2.6)

Consider the associated source and approximate source problems. Given 𝑓∈𝐷, find (𝑤,𝑝)∈𝑉×𝑊 satisfying 𝑎𝑤,𝜓+𝑏𝜓,𝑝=𝑓,𝜓_𝐷𝑏(),∀𝜓∈𝑉,𝑤,𝑣=0,∀𝑣∈𝑊. (2.7) Given 𝑓∈𝐷, find (𝑤_ℎ,𝑝_ℎ)∈𝑉_ℎ×𝑊_ℎ satisfying 𝑎𝑤_ℎ,𝜓+𝑏𝜓,𝑝_ℎ=𝑓,𝜓_𝐷,∀𝜓∈𝑉_ℎ,𝑏𝑤_ℎ,𝑣=0,∀𝑣∈𝑊_ℎ. (2.8) Note that the source term 𝑓 is independent of the solution.

As for the mixed finite element method for boundary value problems, Brezzi and Fortinand so forth have established a systematic theory (see [14]). By Brezzi-Babuska Theorem, we have the following.

Lemma 2.1 (Brezzi-Babuska). Suppose that(C1) (2.1)-(2.2) hold;(C2) inf-sup condition holds, namely, there exists a constant 𝜈₁>0, such that sup_{𝜓∈𝑉,𝜓≠0}𝑏𝜓,𝑣𝜓_𝑉≥𝜈₁‖𝑣‖_𝑊,∀𝑣∈𝑊, (2.9) then there exists a unique solution (𝑤,𝑝) to the problem (2.7) and ‖𝑤‖_𝑎+𝑝_𝑊≤𝐶_𝑟𝑓_𝐷, (2.10) where 𝐶_𝑟 just depends on 𝜈,𝜈₁, and 𝑀₁. Moreover, suppose; (C3) discrete inf-sup condition holds, namely, there exists a constant 𝜈₂>0 independent of ℎ, such that sup_{𝜓∈𝑉ℎ,𝜓≠0}𝑏𝜓,𝑣𝜓_𝑉≥𝜈₂‖𝑣‖_𝑊,∀𝑣∈𝑊_ℎ, (2.11) then there exists a unique solution (𝑤_ℎ,𝑝_ℎ) to the problem (2.8) and the following error estimate is valid: 𝑤−𝑤_ℎ𝑎+𝑝−𝑝_ℎ𝑊≤𝐶_𝑒inf_{𝑞∈𝑉ℎ}𝑤−𝑞_𝑎+inf_{𝑣∈𝑊ℎ}𝑝−𝑣_𝑊, (2.12) where 𝐶_𝑒 just depends on 𝜈,𝜈₂ and 𝑀₁,𝑀₂.

Suppose conditions (C1)–(C3) hold in Lemma 2.1. Then there exist unique solutions to the problem (2.7) and (2.8), respectively. Thus, we can define linear bounded operators as follows: 𝑇∶𝐷→𝑉,𝑆∶𝐷→𝑊∶∀𝑓∈𝐷, 𝑎𝑇𝑓,𝜓+𝑏𝜓,𝑆𝑓=𝑓,𝜓_𝐷𝑏,∀𝜓∈𝑉,𝑇𝑓,𝑣=0,∀𝑣∈𝑊. (2.13)𝑇_ℎ∶𝐷→𝑉_ℎ⊂𝑉,𝑆_ℎ∶𝐷→𝑊_ℎ⊂𝑊∶∀𝑓∈𝐷, 𝑎𝑇_ℎ𝑓,𝜓+𝑏𝜓,𝑆_ℎ𝑓=𝑓,𝜓_𝐷,∀𝜓∈𝑉_ℎ;𝑏𝑇_ℎ𝑓,𝑣=0,∀𝑣∈𝑊_ℎ. (2.14) Obviously, (2.3)-(2.4) has the following equivalent operator form (),𝜆𝑇𝑢=𝑢,𝜎=𝑆𝜆𝑢 (2.15) and (2.5)-(2.6) has the following equivalent operator form 𝜆_ℎ𝑇_ℎ𝑢_ℎ=𝑢_ℎ,𝜎_ℎ=𝑆_ℎ𝜆_ℎ𝑢_ℎ. (2.16)(2.17)

It is easy to verify that 𝑇∶𝐷→𝐷,𝑇_ℎ∶𝐷→𝐷 are self-adjoint operators in the sense of inner product (⋅,⋅)_𝐷 (see [7]).

Assume 𝑉↪^𝑐𝐷 (compactly embedded), then it’s easy to prove that 𝑇∶𝐷→𝐷 is completely continuous, 𝑇∶𝑉→𝑉 is completely continuous, and 𝑇_ℎ is a finite rank operator. Combining (2.3)-(2.4), (2.5)-(2.6), and the 𝑉-ellipticity of 𝑎(⋅,⋅), we deduce 𝑎()𝜆=𝑢,𝑢(𝑢,𝑢)_𝐷>0,𝜆_ℎ=𝑎𝑢_ℎ,𝑢_ℎ𝑢_ℎ,𝑢_ℎ𝐷>0. (2.18) Then from the spectral theory of self-adjoint and completely continuous operator, we know that the eigenvalues of (2.3)-(2.4) can be sorted as 0<𝜆₁≤𝜆₂≤⋯≤𝜆_𝑘≤⋯↗+∞, (2.19) and the corresponding eigenfunctions are 𝑢₁,𝜎₁,𝑢₂,𝜎₂,…,𝑢_𝑘,𝜎_𝑘,…, (2.20) where (𝑢_𝑖,𝑢_𝑗)_𝐷=𝛿_𝑖𝑗.

The eigenvalues of (2.5)-(2.6) can be sorted as 0<𝜆_1,ℎ≤𝜆_2,ℎ≤⋯≤𝜆_𝑘,ℎ≤⋯≤𝜆_𝑁ℎ,ℎ, (2.21) and the corresponding eigenfunctions are 𝑢_1,ℎ,𝜎_1,ℎ,𝑢_2,ℎ,𝜎_2,ℎ,…,𝑢_𝑘,ℎ,𝜎_𝑘,ℎ,…,𝑢_𝑁ℎ,ℎ,𝜎_𝑁ℎ,ℎ, (2.22) where 𝑁_ℎ=dim𝑉_ℎ, (𝑢_𝑖,ℎ,𝑢_𝑗,ℎ)_𝐷=𝛿_𝑖𝑗.

From (2.5), by taking 𝑢_ℎ=𝑢_𝑖,ℎ, 𝜓=𝑢_𝑗,ℎ, we have 𝑎𝑢_𝑖,ℎ,𝑢_𝑗,ℎ+𝑏𝑢_𝑗,ℎ,𝜎_ℎ=𝜆_𝑖,ℎ𝑢_𝑖,ℎ,𝑢_𝑗,ℎ𝐷, (2.23) and from (2.6) we see that 𝑏(𝑢_𝑗,ℎ,𝜎_ℎ)=0, then a𝑢_𝑖,ℎ,𝑢_𝑗,ℎ=𝜆_𝑖,ℎ𝑢_𝑖,ℎ,𝑢_𝑗,ℎ𝐷=𝜆_𝑖,ℎ𝛿_𝑖𝑗, (2.24) thus {𝑢_𝑖,ℎ/‖𝑢_𝑖,ℎ‖_𝑎} is a completely normal eigenvector system on 𝑉_ℎ in the sense of the inner product 𝑎(⋅,⋅).

Denote 𝜆_𝑘=1/𝜇_𝑘, 𝜆_𝑘,ℎ=1/𝜇_𝑘,ℎ. In this paper, 𝜇_𝑘 and 𝜇_𝑘,ℎ, 𝜆_𝑘 and 𝜆_𝑘,ℎ are all called eigenvalues.

Let 𝜇 be the 𝑘th eigenvalue with algebraic multiplicity 𝑞, 𝜇=𝜇_𝑘=𝜇_𝑘+1=⋯=𝜇_{𝑘+𝑞−1}. 𝑀(𝜇) is the space spanned by all eigenfunctions {𝑢_𝑗}_𝑘^{𝑘+𝑞−1} corresponding to 𝜇 of 𝑇. 𝑀_ℎ(𝜇) is the space spanned by all eigenfunctions {𝑢_𝑗,ℎ}_𝑘^{𝑘+𝑞−1} corresponding to all eigenvalues of 𝑇_ℎ that converge to 𝜇. Let 𝑀(𝜇)={𝑣∶𝑣∈𝑀(𝜇),‖𝑣‖_𝑎=1}, 𝑀_ℎ(𝜇)={𝑣∶𝑣∈𝑀_ℎ(𝜇),‖𝑣‖_𝑎=1}. We call 𝜆=1/𝜇 as the 𝑘th eigenvalue, too. Denote 𝑀(𝜆)=𝑀(𝜇), 𝑀_ℎ(𝜆)=𝑀_ℎ(𝜇), and 𝑀(𝜆)=𝑀(𝜇). Define 𝑇−𝑇_ℎ|_{𝑀(𝜆)𝐷}=max_{𝑢∈𝑀(𝜆)}𝑇−𝑇_ℎ𝑢_𝐷‖𝑢‖_𝐷,𝑇−𝑇_ℎ|_{𝑀(𝜆)𝑎}=max_{𝑢∈𝑀(𝜆)}𝑇−𝑇_ℎ𝑢_𝑎‖𝑢‖_𝑎. (2.25)

The convergence and error estimates about the mixed element method of eigenvalue problem have been studied in [7, 12, 13, 15, 16]. From [12], we know that the following results are valid.

Theorem 2.2. Suppose 𝑉↪^𝑐𝐷, 𝑎(𝑢,𝑣) is symmetric, and the conditions of Lemma 2.1 hold; moreover, for any 𝑓∈𝐷, inf_{𝑞∈𝑉ℎ}𝑇𝑓−𝑞_𝑉→0(ℎ→0),inf_{𝑣∈𝑊ℎ}𝑆𝑓−𝑣_𝑊→0(ℎ→0). (2.26)(2.27) Then ‖𝑇_ℎ−𝑇‖_𝐷→0(ℎ→0).

Proof. From 𝑉↪^𝑐𝐷, we derive that 𝑇∶𝐷→𝐷 is a completely continuous operator. It is obvious that 𝑇_ℎ∶𝐷→𝐷 is a finite rank operator. From (2.12), (2.26), and (2.27), we deduce 𝑇𝑓−𝑇_ℎ𝑓_𝑎≤𝐶_𝑒inf_{𝑞∈𝑉ℎ}𝑇𝑓−𝑞_𝑎+inf_{𝑣∈𝑊ℎ}𝑆𝑓−𝑣_𝑊→0(ℎ→0). (2.28) It shows that 𝑇_ℎ∶𝐷→𝑉 pointwisely converges to 𝑇. From (2.10) and (2.12) we derive that both 𝑇∶𝐷→𝑉 and 𝑇_ℎ∶𝐷→𝑉 are linear bounded. Hence, from Banach-Steinhaus Theorem, we know that there exists a positive constant 𝑀 independent of ℎ, such that sup_ℎ𝑇_{ℎ𝐿(𝐷,𝑉)}≤𝑀. (2.29) Thus, ∪_ℎ>0(𝑇−𝑇_ℎ)𝐵 is a bounded set in 𝑉 with respect to the unit ball 𝐵 of 𝐷. From 𝑉↪^𝑐𝐷, we know that ∪_ℎ>0(𝑇−𝑇_ℎ)𝐵 is a relatively compact set in 𝐷, which proves that {𝑇_ℎ} is collectively compact. From (2.28), we know that 𝑇∶𝐷→𝐷,𝑇_ℎ∶𝐷→𝐷, and 𝑇_ℎ pointwisely converge to 𝑇. From [7], 𝑇∶𝐷→𝐷,𝑇_ℎ∶𝐷→𝐷 are self-adjoint operators in the sense of inner product (⋅,⋅)_𝐷. Then by Lemma 3.7 or Table  3.1 in [17], we get ‖𝑇_ℎ−𝑇‖_𝐷→0(ℎ→0). The proof is completed.

Lemma 2.3. Suppose that the conditions of Theorem 2.2 hold. Let (𝜆_ℎ,𝑢_ℎ,𝜎_ℎ) be the 𝑘th eigenpair of (2.5)-(2.6) with ‖𝑢_ℎ‖_𝑎=1, let 𝜆 be the 𝑘th eigenvalue of (2.3)-(2.4). Then 𝜆_ℎ→𝜆(ℎ→0), and there exists an eigenfunction (𝑢,𝜎) corresponding to 𝜆, such that 𝜆_ℎ−𝜆+𝑢_ℎ−𝑢_𝐷≤𝐶₁𝑇−𝑇_ℎ|_{𝑀(𝜆)𝐷},𝜎−𝜎_ℎ𝑊≤𝑆_ℎ(𝜆𝑢)−𝑆(𝜆𝑢)_𝑊+𝐶𝑇−𝑇_ℎ|_{𝑀(𝜆)𝐷},𝑢−𝑢_ℎ𝑎≤𝐶₂𝑇_ℎ−𝑇|_{𝑀(𝜆)𝑎}, (2.30)(2.31)(2.32) let 𝑢∈𝑀(𝜆), then there exists 𝑢_ℎ∈𝑀_ℎ(𝜆) such that 𝑢−𝑢_ℎ𝑎≤𝐶₃𝑇_ℎ−𝑇|_{𝑀(𝜆)𝑎}, (2.33) where 𝑢 depends on ℎ in general, and 𝐶₁, 𝐶₂, and 𝐶₃ are constants independent of ℎ.

Proof. By Theorem 2.2, we know ‖𝑇_ℎ−𝑇‖_𝐷→0(ℎ→0). Thus from Theorem  2.2 in [7], we see that the desired results are valid. The proof is completed.

For (𝑢^∗,𝜎^∗)∈𝑉×𝑊, 𝑢^∗≠0, define the Rayleigh quotient 𝜆^𝑟=𝑎𝑢^∗,𝑢^∗+2𝑏𝑢^∗,𝜎^∗(𝑢^∗,𝑢^∗)_𝐷. (2.34)

Lemma 2.4. Suppose (𝜆,𝑢,𝜎) is an eigenpair of (2.3)-(2.4), then for any (𝑢^∗,𝜎^∗)∈𝑉×𝑊, 𝑢^∗≠0, the Rayleigh quotient 𝜆^𝑟 satisfies 𝜆^𝑟𝑎𝑢−𝜆=^∗−𝑢,𝑢^∗−𝑢+2𝑏𝑢^∗−𝑢,𝜎^∗−𝜎(𝑢^∗,𝑢^∗)_𝐷𝑢−𝜆^∗−𝑢,𝑢^∗−𝑢_𝐷(𝑢^∗,𝑢^∗)_𝐷. (2.35)

Proof. The proof is completed by using the same argument as that of Lemma 9.1 (see [7, 18]).

Since 𝑉↪𝐷 (continuously imbedded), there exists a constant 𝐶₄ independent of ℎ such that ‖𝑣‖_𝐷≤𝐶₄‖𝑣‖_𝑎,∀𝑣∈𝑉. (2.36) Taking (𝑢^∗,𝜎^∗)=(𝑢_ℎ,𝜎_ℎ) in (2.35) and using (2.4) and (2.6), we deduce the following.

Lemma 2.5. Suppose (𝜆_ℎ,𝑢_ℎ,𝜎_ℎ) is an approximation of (𝜆,𝑢,𝜎)  and ‖𝑢_ℎ‖_𝑎=1, then 𝜆_ℎ𝑎𝑢−𝜆=^ℎ−𝑢,𝑢_ℎ−𝑢+2𝑏𝑢_ℎ−𝑢,𝑣−𝜎𝑢_ℎ,𝑢_ℎ𝐷𝑢−𝜆^ℎ−𝑢,𝑢_ℎ−𝑢_𝐷𝑢_ℎ,𝑢_ℎ𝐷,∀𝑣∈W_ℎ,𝜆_ℎ−𝜆≤𝜆_ℎ+𝜆𝜆_ℎ𝐶²₄𝑢_ℎ−𝑢²_𝑎+2𝜆_ℎ𝑀₂𝑢_ℎ−𝑢_𝑉‖𝜎−𝑣‖_𝑊,∀𝑣∈𝑊_ℎ. (2.37)(2.38)

3. Mixed Finite Element Multigrid Discretization Scheme Based on the Rayleigh Quotient Iteration

In this section, we develop the work in [8], noticing that in [8], when 𝑘>1,𝜆^ℎ1_𝑘−𝜆_𝑘,𝜆^ℎ𝑙−1_𝑘−𝜆_𝑘,ℎ𝑙 and 𝜆^ℎl−1_𝑘−𝜆_𝑘 sholud be modified to their absolute values, respectively, and establish the following mixed finite element multi-grid discretization scheme based on the Rayleigh quotient iteration, and give a rigorous theoretical analysis. Suppose the partition satisfies the following conditions.

Condition 1. {𝐾^ℎ𝑖} is a family of regular meshes (see [19]) with the mesh diameter {ℎ_𝑖} and ℎ_𝑖=ℎ^𝑡𝑖_𝑖−1, 𝑡_𝑖∈(1,3) is arbitrarily chosen, 𝑖=1,2,…, and inf_𝑖𝑡_𝑖>1.
Let {𝑉_ℎ𝑖}^𝑙₀ and {𝑊_ℎ𝑖}^𝑙₀ be finite element spaces on {𝐾^ℎ𝑖}^𝑙₀.

Scheme 1. Multigrid Discretization.

Step 1. Solve the eigenvalue problem (2.3)-(2.4) on 𝑉_𝐻×𝑊_𝐻: find (𝜆_𝐻,𝑢_𝐻,𝜎_𝐻)∈𝑅×𝑉_𝐻×𝑊_𝐻, ‖𝑢_𝐻‖_𝑎=1 such that 𝑎𝑢_𝐻,𝜓+𝑏𝜓,𝜎_𝐻=𝜆_𝐻𝑢_𝐻,𝜓_𝐷,∀𝜓∈𝑉_𝐻,𝑏𝑢_𝐻,𝑣=0,∀𝑣∈𝑊_𝐻. (3.1)

Step 2. 𝑢^ℎ0⇐𝑢_𝐻, 𝜆^ℎ0⇐𝜆_𝐻, 𝑖⇐1.

Step 3. Solve an equation on 𝑉_ℎ𝑖×𝑊_ℎ𝑖: find (𝑢^′,𝜎^′)∈𝑉_ℎ𝑖×𝑊_ℎ𝑖 such that 𝑎𝑢^′,𝜓+𝑏𝜓,𝜎^′−𝜆^ℎ𝑖−1𝑢^′,𝜓_𝐷=𝑢^ℎ𝑖−1,𝜓_𝐷,∀𝜓∈𝑉_ℎ𝑖,𝑏𝑢^′,𝑣=0,∀𝑣∈𝑊_ℎ𝑖. (3.2) Set 𝑢^ℎ𝑖=𝑢^′/‖𝑢^′‖_𝑎, 𝜎^ℎ𝑖=𝜎^′/‖𝑢^′‖_𝑎.

Step 4. Compute the Rayleigh quotient 𝜆^ℎ𝑖=𝑎𝑢^ℎ𝑖,𝑢^ℎ𝑖𝑢^ℎ𝑖,𝑢^ℎ𝑖_𝐷. (3.3)

Step 5. If 𝑖=𝑙, then output (𝜆^ℎ𝑙,𝑢^ℎ𝑙,𝜎^ℎ𝑙), stop; else, 𝑖⇐𝑖+1, and return to Step 3.
Let (𝜆_𝐻,𝑢_𝐻,𝜎_𝐻) be the 𝑘th eigenpair, and we use (𝜆^ℎ𝑙,𝑢^ℎ𝑙,𝜎^ℎ𝑙) as the 𝑘th approximation eigenpair of (2.3)-(2.4).
Next, we will discuss the validity of Scheme 1.

Lemma 3.1. For any nonzero 𝑢,𝑣∈𝑉, 𝑢‖𝑢‖_𝑎−𝑣‖𝑣‖_𝑎a‖‖≤2𝑢−𝑣_𝑎‖𝑢‖_𝑎,𝑢‖𝑢‖_𝑎−𝑣‖𝑣‖_𝑎a‖‖≤2𝑢−𝑣_𝑎‖𝑣‖_𝑎. (3.4)

Proof. See [8].

Denote dist(𝑢,𝑉)=inf_𝑣∈𝑉‖𝑢−𝑣‖_𝑎.

Consider the eigenvalue problem (2.16) on 𝑉_ℎ.

Lemma 3.2. Suppose that 𝜇 and 𝜇_ℎ are the 𝑘th eigenvalue of 𝑇 and 𝑇_ℎ, respectively, and (𝜇₀,𝑢₀) is an approximate eigenpair, where 𝜇₀ is not an eigenvalue of 𝑇_ℎ, 𝑢₀∈𝑉_ℎ,‖𝑢₀‖_𝑎=1, dist(𝑢₀,𝑀_ℎ(𝜇))≤1/2, max_{𝑘≤𝑗≤𝑘+𝑞−1}|(𝜇_𝑗,ℎ−𝜇_ℎ)/(𝜇₀−𝜇_𝑗,ℎ)|≤1/2, |𝜇₀−𝜇_𝑗,ℎ|≥(𝜌/2)(𝑗≠𝑘,𝑘+1,…,𝑘+𝑞−1), and 𝑢^𝑠∈𝑉_ℎ,𝑢^ℎ∈𝑉_ℎ satisfy 𝜇₀−𝑇_ℎ𝑢^𝑠=𝑢₀,𝑢^ℎ=𝑢^𝑠‖𝑢^𝑠‖_𝑎. (3.5) Then dist𝑢^ℎ,𝑀_ℎ𝜇≤16𝜌𝜇₀−𝜇_ℎdist𝑢₀,𝑀_ℎ𝜇, (3.6) where 𝜌=min_{𝜇𝑗≠𝜇}|𝜇_𝑗−𝜇| is the separation constant of the eigenvalue 𝜇.

Proof. See [8]

Since the convergence rate of 𝑉_ℎ𝑙−1 and 𝑊_ℎ𝑙−1 approximating eigenfunctions is lower than that of 𝑉_ℎ𝑙 and 𝑊_ℎ𝑙 approximating eigenfunctions, respectively, the approximation order of (𝜆_ℎ𝑙−1,𝑢_ℎ𝑙−1) is lower than that of (𝜆_ℎ𝑙,𝑢_ℎ𝑙). However, in general, the accuracy order of (𝜆^ℎ𝑙−1,𝑢^ℎ𝑙−1) will not exceed that of (𝜆_ℎ𝑙−1,𝑢_ℎ𝑙−1); therefore in the following Theorem 3.3 we assume that the accuracy order of (𝜆^ℎ𝑙−1,𝑢^ℎ𝑙−1) is lower than that of (𝜆_ℎ𝑙,𝑢_ℎ𝑙).

Theorem 3.3. Suppose that ‖𝑇_ℎ−𝑇‖_𝐷→0(ℎ→0), 𝐻 is small properly, and Condition 1 holds. Let (𝜆^ℎ𝑙,𝑢^ℎ𝑙,𝜎^ℎ𝑙) be the approximate eigenpair obtained by Scheme 1, and let 𝑢^ℎ𝑙−1 approximate 𝑢∈𝑀(𝜆), 𝜆^ℎ𝑙−1 approximate 𝜆, and the accuracy order of (𝜆^ℎ𝑙−1,𝑢^ℎ𝑙−1) be lower than that of (𝜆_ℎ𝑙,𝑢_ℎ𝑙). Then there exists 𝑢∈𝑀(𝜆) such that 𝑢^ℎ𝑙−𝑢_𝑎≤32𝜌𝐶₅𝐶₆𝜆^ℎ𝑙−1−𝜆²+𝜆^ℎ𝑙−1−𝜆𝑢^ℎ𝑙−1−𝑢_𝐷+𝐶₂×𝑞+3𝑇−𝑇_ℎ𝑙|_{𝑀(𝜆)𝑎},𝜆^ℎ𝑙−𝜆≤𝜆^ℎ𝑙+𝜆𝜆^ℎ𝑙𝐶²₄𝑢^ℎ𝑙−𝑢²_𝑎+2𝜆^ℎ𝑙𝑀₂𝑢^ℎ𝑙−𝑢_𝑉inf_{𝑣∈𝑊ℎ𝑙}‖𝜎−𝑣‖_𝑊, (3.7)(3.8) where 𝐶₅ and 𝐶₆ are determined, respectively, by (3.10) and (3.14) in the following proof.

Proof. Let 𝜇₀=1/𝜆^ℎ𝑙−1, and 𝑢₀=𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^ℎ𝑙−1/‖𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^ℎ𝑙−1‖_𝑎. Since 𝑢∈𝑀(𝜆), by calculation we deduce 𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^ℎ𝑙−1−𝑢_𝑎=𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^ℎ𝑙−1−𝜆𝑇𝑢_𝑎=𝜆^ℎ𝑙−1−𝜆𝑇_ℎ𝑙𝑢^ℎ𝑙−1_𝑎+𝜆𝑇_ℎ𝑙𝑢^ℎ𝑙−1−𝑢_𝑎+𝜆𝑇_ℎ𝑙−𝑇𝑢_𝑎. (3.9) From Lemma 2.1, there exists a positive constant 𝐶₅ depending only on 𝜈,𝜈₁,𝜈₂,𝑀₁, and 𝑀₂ such that 𝑇_ℎ𝑙𝑣_𝑎≤𝐶₅‖𝑣‖_𝐷,∀𝑣∈𝐷. (3.10) Then 𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^ℎ𝑙−1−𝑢_𝑎≤𝐶₅𝐶₄𝜆^ℎ𝑙−1−𝜆+𝜆𝑢^ℎ𝑙−1−𝑢_𝐷+𝜆𝑇_ℎ𝑙−𝑇|_{𝑀(𝜆)𝑎}. (3.11) By Lemma 3.1, we derive dist𝑢₀,𝑀(𝜆)≤𝑢₀−𝑢_𝑎≤2𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^ℎ𝑙−1−𝑢_𝑎≤2𝐶₅𝐶₄𝜆^ℎ𝑙−1−𝜆+𝜆𝑢^ℎ𝑙−1−𝑢_𝐷+2𝜆𝑇_ℎ𝑙−𝑇|_{𝑀(𝜆)𝑎}. (3.12) Using the triangle inequality and (2.33), we deduce dist𝑢₀,𝑀_ℎ𝑙(𝜆)≤dist𝑢₀,𝑀(𝜆)+𝐶₃𝑇_ℎ𝑙−𝑇|_{𝑀(𝜆)𝑎}. (3.13) According to the hypotheses of the theorem, we know that 𝜆^ℎ𝑙−1→𝜆 and 𝜆^ℎ𝑙−1−𝜆 are an infinitesimal of lower order comparing with 𝜆_𝑗,ℎ𝑙−𝜆. Hence, there exists a positive constant 𝐶₆ independent of ℎ_𝑙(𝑙=1,2,…) such that for 𝑗=𝑘,𝑘+1,…,𝑘+𝑞−1 we have 𝜇₀−𝜇_𝑗,ℎ𝑙=𝜆^ℎ𝑙−1−𝜆+𝜆−𝜆_𝑗,ℎ𝑙𝜆_𝑗,ℎ𝑙𝜆^ℎ𝑙−1≤𝐶₆𝜆^ℎ𝑙−1−𝜆. (3.14) Note that 𝐻 is small enough and ℎ_𝑙≪ℎ_𝑙−1; from (3.13) and (3.12), we obtain dist𝑢₀,𝑀_ℎ𝑙(𝜆≤1)2. (3.15) Noticing that 𝜆=𝜆_𝑘=𝜆_𝑘+1=⋯=𝜆_{𝑘+𝑞−1}, we have 𝜇_𝑗,ℎ𝑙−𝜇_ℎ𝑙=𝜆_ℎ𝑙−𝜆_𝑗,ℎ𝑙𝜆_ℎ𝑙𝜆_𝑗,ℎ𝑙=𝜆_ℎ𝑙−𝜆+𝜆_𝑗−𝜆_𝑗,ℎ𝑙𝜆_ℎ𝑙𝜆_𝑗,ℎ𝑙, (3.16) which together with (3.14), noting that 𝜆_𝑗,ℎ𝑙−𝜆 is an infinitesimal of higher order comparing with 𝜆^ℎ𝑙−1−𝜆, yields max_{𝑘≤𝑗≤𝑘+𝑞−1}𝜇_𝑗,ℎ𝑙−𝜇_ℎ𝑙𝜇₀−𝜇_𝑗,ℎ𝑙≤12. (3.17) Since 𝜌 is the separation constant, 𝐻 is small enough, and ℎ_𝑙≪ℎ_𝑙−1, there holds 𝜇₀−𝜇_𝑗,ℎ𝑙≥𝜌2,𝑗≠𝑘,𝑘+1,…,𝑘+𝑞−1. (3.18) For 𝑢^′ in Step 3 of Scheme 1, from the definition of 𝑇_ℎ and 𝑆_ℎ (taking 𝑖=𝑙), we have 𝑎𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^′,𝜓+𝑏𝜓,𝜆^ℎ𝑙−1𝑆_ℎ𝑙𝑢^′=𝜆^ℎ𝑙−1𝑢^′,𝜓_𝐷,∀𝜓∈𝑉_ℎ𝑙,𝑏𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^′,𝑣=0,∀𝑣∈𝑊_ℎ𝑙a𝑇_ℎ𝑙,𝑢^ℎ𝑙−1,𝜓+𝑏𝜓,𝑆_ℎ𝑙𝑢^ℎ𝑙−1=𝑢^ℎ𝑙−1,𝜓_𝐷,∀𝜓∈𝑉_ℎ𝑙,𝑏𝑇_ℎ𝑙𝑢^ℎ𝑙−1,𝑣=0,∀𝑣∈𝑊_ℎ𝑙. (3.19)(3.20)(3.21)(3.22) Hence, when 𝑖=𝑙, Step 3 of Scheme 1 is equivalent to the following: find (𝑢^′,𝜎^′)∈𝑉_ℎ𝑙×𝑊_ℎ𝑙 such that 𝑎𝑢^′,𝜓+𝑏𝜓,𝜎^′−𝜆^ℎ𝑙−1𝑎𝑇_ℎ𝑙𝑢^′,𝜓−𝜆^ℎ𝑙−1𝑏𝜓,𝑆_ℎ𝑙𝑢^′=𝑎𝑇_ℎ𝑙𝑢^ℎ𝑙−1,𝜓+𝑏𝜓,𝑆_ℎ𝑙𝑢^ℎ𝑙−1,∀𝜓∈𝑉_ℎ𝑙,𝑏𝑢^′,𝑣=0,∀𝑣∈𝑊_ℎ𝑙. (3.23)(3.24) And set 𝑢^ℎ𝑙=𝑢^′/‖𝑢^′‖_𝑎, 𝜎^ℎ𝑙=𝜎^′/‖𝑢^′‖_𝑎.
From (3.23), we obtain 𝑎𝑢^′−𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^′−𝑇_ℎ𝑙𝑢^ℎ𝑙−1,𝜓+𝑏𝜓,𝜎^′−𝜆^ℎ𝑙−1𝑆_ℎ𝑙𝑢^′−𝑆_ℎ𝑙𝑢^ℎ𝑙−1=0,∀𝜓∈𝑉_ℎ𝑙. (3.25) Combining (3.24), (3.20) and (3.22), we get 𝑏𝑢^′−𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^′−𝑇_ℎ𝑙𝑢^ℎ𝑙−1,𝑣=0,∀𝑣∈𝑊_ℎ𝑙. (3.26) By (3.26) and taking 𝜓=𝑢^′−𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^′−𝑇_ℎ𝑙𝑢^ℎ𝑙−1 in (3.25), we obtain 𝑎𝑢^′−𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^′−𝑇_ℎ𝑙𝑢^ℎ𝑙−1,𝑢^′−𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^′−𝑇_ℎ𝑙𝑢^ℎ𝑙−1=0. (3.27) Thus 1𝜆^ℎ𝑙−1−𝑇_ℎ𝑙𝑢^′=1𝜆^ℎ𝑙−1𝑇_ℎ𝑙𝑢^ℎ𝑙−1,𝑢^ℎ𝑙=𝑢^′𝑢^′_𝑎. (3.28) From (3.28) we know that the first term on the left-hand side of (3.25) is equal to 0, thus 𝑏𝜓,𝜎^′−𝜆^ℎ𝑙−1𝑆_ℎ𝑙𝑢^′−𝑆_ℎ𝑙𝑢^ℎ𝑙−1=0,∀𝜓∈𝑉_ℎ𝑙, (3.29) then, using discrete inf-sup condition we get 𝜎^′=𝜆^ℎ𝑙−1𝑆_ℎ𝑙𝑢^′+𝑆_ℎ𝑙𝑢^ℎ𝑙−1. (3.30) Thus Step 3 of Scheme 1 is equivalent to (3.28), (3.30), 𝑢^ℎ𝑙=𝑢^′/‖𝑢^′‖_𝑎, and 𝜎^ℎ𝑙=𝜎^′/‖𝑢^′‖_𝑎.
Noting that (1/𝜆^ℎ𝑙−1)𝑇_ℎ𝑙𝑢^ℎ𝑙−1=‖(1/𝜆^ℎ𝑙−1)𝑇_ℎ𝑙𝑢^ℎ𝑙−1‖_𝑎𝑢₀ differs from 𝑢₀ by only a constant and selecting 𝑢^𝑠=(𝜆^ℎ𝑙−1𝑢^′)/‖𝑇_ℎ𝑙𝑢^ℎ𝑙−1‖_𝑎, we have 1𝜆^ℎ𝑙−1−𝑇_ℎ𝑙𝑢^𝑠=𝑢₀,𝑢^ℎ𝑙=𝑢^𝑠‖𝑢^𝑠‖_𝑎. (3.31) By (3.15), (3.17), (3.18), and (3.31), we see that the conditions of Lemma 3.2 hold. Thus, substituting (3.13) and (3.14) into (3.6), we obtain dist𝑢^ℎ𝑙,𝑀_ℎ𝑙(𝜆)≤16𝜌𝐶₆𝜆^ℎ𝑙−1−𝜆dist𝑢₀,𝑀(𝜆)+𝐶₃𝑇−𝑇_ℎ𝑙|_{𝑀(𝜆)𝑎}. (3.32) Let the eigenfunctions {𝑢_𝑗,ℎ𝑙}_𝑘^{𝑘+𝑞−1} be an orthonormal system of 𝑀_ℎ𝑙(𝜆) in the sense of 𝑎(⋅,⋅), then dist𝑢^ℎ𝑙,𝑀_ℎ𝑙(𝜆)=𝑢^ℎ𝑙−^{𝑘+𝑞−1}_𝑗=𝑘𝑎𝑢^ℎ𝑙,𝑢_𝑗,ℎ𝑙𝑢_{𝑗,ℎ𝑙𝑎}. (3.33) Let 𝑢^∗=^{𝑘+𝑞−1}_𝑗=𝑘𝑎𝑢^ℎ𝑙,𝑢_𝑗,ℎ𝑙𝑢_𝑗,ℎ𝑙, (3.34) noting ‖𝑢^ℎ𝑙−𝑢^∗‖_𝑎≤dist(𝑢^ℎ𝑙,𝑀_ℎ𝑙(𝜆)), from (3.32) we deduce 𝑢^ℎ𝑙−𝑢^∗_𝑎≤16𝜌𝐶₆𝜆^ℎ𝑙−1−𝜆dist𝑢₀,𝑀(𝜆)+𝐶₃𝑇−𝑇_ℎ𝑙|_{𝑀(𝜆)𝑎}. (3.35) By Lemma 2.3, there exist {𝑢⁰_𝑗}_𝑘^{𝑘+𝑞−1}⊂𝑀(𝜆) such that 𝑢_𝑗,ℎ𝑙−𝑢⁰_𝑗 satisfy (2.32). Let 𝑢=^{𝑘+𝑞−1}_𝑗=𝑘𝑎𝑢^ℎ𝑙,𝑢_𝑗,ℎ𝑙𝑢⁰_𝑗, (3.36) then 𝑢∈𝑀(𝜆). Using (2.32), we deduce 𝑢^∗−𝑢_𝑎=^{𝑘+𝑞−1}_𝑗=𝑘𝑎𝑢^ℎ𝑙,𝑢_𝑗,ℎ𝑙𝑢_𝑗,ℎ𝑙−𝑢⁰_𝑗𝑎≤𝐶₂𝑞𝑇_ℎ𝑙−𝑇|_{𝑀(𝜆)𝑎}. (3.37) Combining (3.35) with the above inequality, we have 𝑢^ℎ𝑙−𝑢_𝑎≤16𝜌𝐶₆𝜆^ℎ𝑙−1−𝜆dist𝑢₀,𝑀(𝜆)+𝐶₂𝑞+1𝑇−𝑇_ℎ𝑙|_{𝑀(𝜆)𝑎}. (3.38) Substituting (3.12) into (3.38), we get (3.7).
From Step 3 of Scheme 1, we know that 𝑏(𝑢^ℎ𝑙,𝜎^ℎ𝑙)=0, thus 𝜆^ℎ𝑙=𝑎𝑢^ℎ𝑙,𝑢^ℎ𝑙𝑢^ℎ𝑙,𝑢^ℎ𝑙_𝐷=𝑎𝑢^ℎ𝑙,𝑢^ℎ𝑙+2𝑏𝑢^ℎ𝑙,𝜎^ℎ𝑙𝑢^ℎ𝑙,𝑢^ℎ𝑙_𝐷. (3.39) Select 𝜆^𝑟=𝜆^ℎ𝑙, 𝑢^∗=𝑢^ℎ𝑙, and 𝜎^∗=𝜎^ℎ𝑙. From Lemma 2.4, we get 𝜆^ℎ𝑙𝑎𝑢−𝜆=^ℎ𝑙−𝑢,𝑢^ℎ𝑙−𝑢+2𝑏𝑢^ℎ𝑙−𝑢,𝜎^ℎ𝑙−𝜎𝑢^ℎ𝑙,𝑢^ℎ𝑙_𝐷𝑢−𝜆^ℎ𝑙−𝑢,𝑢^ℎ𝑙−𝑢_𝐷𝑢^ℎ𝑙,𝑢^ℎ𝑙_𝐷. (3.40) Noting that ∀𝑣∈𝑊_ℎ𝑙, 𝑏(𝑢^ℎ𝑙−𝑢,𝑣)=0, and (𝑢^ℎ𝑙,𝑢^ℎ𝑙)_𝐷=(1/𝜆^ℎ𝑙)𝑎(𝑢^ℎ𝑙,𝑢^ℎ𝑙)=1/𝜆^ℎ𝑙, we have 𝜆^ℎ𝑙−𝜆=𝜆^ℎ𝑙𝑎𝑢^ℎ𝑙−𝑢,𝑢^ℎ𝑙−𝑢+2𝜆^ℎ𝑙𝑏𝑢^ℎ𝑙−𝑢,𝑣−𝜎−𝜆𝜆^ℎ𝑙𝑢^ℎ𝑙−𝑢,𝑢^ℎ𝑙−𝑢_𝐷,∀𝑣∈𝑊_ℎ𝑙. (3.41) Note ‖𝑢^ℎ𝑙−𝑢‖_𝐷≤𝐶₄‖𝑢^ℎ𝑙−𝑢‖_𝑎. From (3.41) we obtain (3.8).

For 𝑙=1, Scheme 1 is actually the two-grid discretization scheme established in [7]. By Theorem 3.3 we get the following conclusion directly.

Theorem 3.4. Suppose ‖𝑇_ℎ−𝑇‖_𝐷→0(ℎ→0). Let ℎ₀=𝐻 be small properly. Let (𝜆^ℎ1,𝑢^ℎ1,𝜎^ℎ1) be an approximate eigenpair obtained by Scheme 1 (𝑙=1). Then there exists 𝑢∈𝑀(𝜆) such that 𝑢^ℎ1−𝑢_𝑎≤32𝜌𝐶₅𝐶₆𝜆^ℎ0−𝜆²+𝜆^ℎ0−𝜆𝐶₁𝑇−𝑇_ℎ0|_{𝑀(𝜆)𝐷}+𝐶₂×𝑞+3𝑇−𝑇_ℎ1|_{𝑀(𝜆)𝑎},𝜆^ℎ1−𝜆≤𝜆^ℎ1+𝜆𝜆^ℎ1𝐶²₄𝑢^ℎ1−𝑢²_𝑎+2𝜆^ℎ1𝑀₂𝑢^ℎ1−𝑢_𝑉inf_{𝑣∈𝑊ℎ1}‖𝜎−𝑣‖_𝑊. (3.42)(3.43)

Proof. Consider Scheme 1. Here 𝑙=1, 𝑢^ℎ0=𝑢_𝐻. By Lemma 2.3, we know that 𝑢^ℎ0 approximates 𝑢∈𝑀(𝜆), and the accuracy order of (𝜆^ℎ0,𝑢^ℎ0) is lower than (𝜆_ℎ1,𝑢_ℎ1). Hence, for 𝑙=1, the conditions of Theorem 3.3 hold. Select 𝑢 in the proof of Theorem 3.3 such that 𝑢^ℎ0−𝑢 satisfies Lemma 2.3. Then from (2.30), we obtain 𝑢^ℎ0−𝑢_𝐷≤𝐶₁𝑇−𝑇_ℎ0|_{𝑀(𝜆)𝐷}, (3.44) substituting it into (3.7), we obtain (3.42). From (3.8), we deduce (3.43).

Theorem 3.4 is actually Theorem  3.3 in [7], but we analyze in detail the constants appeared in the error estimates.

From Theorem 3.4 and Theorem 3.3, we know that 𝜆^ℎ𝑙→𝜆 (𝑙→∞) and the convergence rate is high. Thus, we use 𝜃×|𝜆^ℎ𝑙−𝜆^ℎ𝑙−1| as a posteriori error indicator of 𝜆^ℎ𝑙−𝜆 (details can be seen in Remark 4.5). Then we establish the following adaptive algorithm.

Scheme 2 (Adaptive Algorithm). Give an error tolerance 𝜀 and choose the parameter 0<𝜃≤1, 𝐻, 𝑡₁, and ℎ₁=𝐻^𝑡1.

Step 1. Solve (2.3)-(2.4) on 𝑉_𝐻×𝑊_𝐻: find (𝜆_𝐻,𝑢_𝐻,𝜎_𝐻)∈𝑅×𝑉_𝐻×𝑊_𝐻, ‖𝑢_𝐻‖_𝑎=1 such that 𝑎𝑢_𝐻,𝜓+𝑏𝜓,𝜎_𝐻=𝜆_𝐻𝑢_𝐻,𝜓_𝐷,∀𝜓∈𝑉_𝐻,𝑏𝑢_𝐻,𝑣=0,∀𝑣∈𝑊_𝐻. (3.45)

Step 2. 𝑢^ℎ0⇐𝑢_𝐻, 𝜆^ℎ0⇐𝜆_𝐻, 𝑙⇐1.

Step 3. Solve an equation on 𝑉_ℎ𝑙×W_ℎ𝑙: find (𝑢^′,𝜎^′)∈𝑉_ℎ𝑙×𝑊_ℎ𝑙 such that 𝑎𝑢^′,𝜓+𝑏𝜓,𝜎^′−𝜆^ℎ𝑙−1𝑢^′,𝜓_𝐷=𝑢^ℎ𝑙−1,𝜓_𝐷,∀𝜓∈𝑉_ℎ𝑙,𝑏𝑢^′,𝑣=0,∀𝑣∈𝑊_ℎ𝑙. (3.46) Let 𝑢^ℎ𝑙=𝑢^′/‖𝑢^′‖_𝑎, 𝜎^ℎ𝑙=𝜎^′/‖𝑢^′‖_𝑎.

Step 4. Compute the Rayleigh quotient 𝜆^ℎ𝑙=𝑎𝑢^ℎ𝑙,𝑢^ℎ𝑙𝑢^ℎ𝑙,𝑢^ℎ𝑙_𝐷. (3.47)

Step 5. If 𝜃×|𝜆^ℎ𝑙−𝜆^ℎ𝑙−1|>𝜀, then select 𝑡_𝑙+1, ℎ_𝑙+1=ℎ^𝑡𝑙+1_𝑙, 𝑙⇐𝑙+1 and return to Step 3; else output (𝜆^ℎ𝑙,𝑢^ℎ𝑙,𝜎^ℎ𝑙), stop.

Remark 3.5. In Scheme 2, we use 𝜃×|𝜆^ℎ𝑙−𝜆^ℎ𝑙−1| as a posteriori error indicator which is global. In order to cope with difficulties caused by local singularity of a complicated problem in calculation, so far, most algorithms designed a local a posteriori error indicator to establish adaptive algorithm with local mesh refinement (e.g., see [20–22]). However, because the eigenvalue problem of electric field is so complicated, that it is very difficult to obtain a local a posteriori error indicator of eigenfunction. Fortunately, the influence of local singularity can be avoided by using the weight method which is established by Costabel and Dauge to discrete the eigenvalue problem of electric field. And the performance of the weigh method is very good (see [3, 11]). Hence, without local mesh refinement, by using the weight method mentioned above our algorithm can also guarantee its high efficiency.

4. The Eigenvalue Problem of Electric Field

Consider the eigenvalue problem of electric field: 𝑐²𝐜𝐮𝐫𝐥𝐜𝐮𝐫𝐥^→𝑢=𝜔^2→𝑢,inΩ,div^→𝑢=0,inΩ,_→𝑢×^→𝛾=0,on∈𝜕Ω, (4.1) where Ω is a polyhedron in 𝐑^𝟑, ^→𝛾 is the unit outward normal to 𝜕Ω.

Physically ^→𝑢 denotes the electric field, 𝜔 denotes the time frequency, and 𝑐 is the speed of the light velocity. Usually, let 𝜆=𝜔²/𝑐² named eigenvalue.

Let 𝐻(𝐜𝐮𝐫𝐥,Ω)=^→𝑞∈𝐿₂(Ω)³∶𝐜𝐮𝐫𝐥^→𝑞∈𝐿₂(Ω)³,𝐻₀(𝐜𝐮𝐫𝐥,Ω)=^→𝑞∈𝐻(𝐜𝐮𝐫𝐥,Ω)∶^→𝑞×^→𝛾|_𝜕Ω=0. (4.2)

When Ω is a convex polyhedron, we define the following function space: 𝜒=^→𝑞∈𝐻₀(𝐜𝐮𝐫𝐥,Ω)∶div^→𝑞∈𝐿²(Ω). (4.3) Denote ^→𝑞,^→𝜓₀=_Ω^→𝑞⋅^→𝜓d𝑥,^→𝑞₀=^→𝑞,^→𝑞₀^1/2._→𝑞,^→𝜓_𝜒=𝐜𝐮𝐫𝐥^→𝑞,𝐜𝐮𝐫𝐥^→𝜓₀+div^→𝑞,div^→𝜓₀,^→𝑞_𝜒=^→𝑞,^→𝑞_𝜒^1/2, (4.4) From [23, 24], we know that 𝜒⊂𝐻¹(Ω)³; (^→𝑞,^→𝜓)_𝜒 is a coercive bilinear form on 𝜒, and ‖^→𝑞‖_𝜒 is a norm.

When Ω is a nonconvex polyhedron, the problem is relatively complicated. Let 𝐸 denote a set of edges of reentrant dihedral angles on 𝜕Ω, and 𝑑=𝑑(𝑥) denote the distance to the set 𝐸: 𝑑(𝑥)=dist(𝑥,∪_𝑒∈𝐸𝑒). We introduce a weight function 𝜔_𝑟 which is a nonnegative smooth function of 𝑥. It can be represented by 𝑑^𝑟 in reentrant edge and angular domain. We shall write 𝜔_𝑟≃𝑑^𝑟. Define the weighted functional spaces: 𝐿²_𝑟(Ω)=𝑣∈𝐿²_loc(Ω)∶𝜔_𝑟𝑣∈𝐿₂(Ω,𝜒_𝑟)=^→𝑞∈𝐻₀(𝐜𝐮𝐫𝐥,Ω)∶div^→𝑞∈𝐿²_𝑟(Ω). (4.5) Denote ^→𝑞,^→𝜓_L2𝑟=_Ω𝜔²_𝑟^→𝑞⋅^→𝜓d𝑥,^→𝑞_𝐿2𝑟=^→𝑞,^→𝑞_𝐿^1/2_2𝑟._→𝑞,^→𝜓_𝜒𝑟=𝐜𝐮𝐫𝐥^→𝑞,𝐜𝐮𝐫𝐥^→𝜓₀+div^→𝑞,div^→𝜓_𝐿2𝑟,^→𝑞_𝜒𝑟=^→𝑞,^→𝑞_𝜒^1/2_𝑟. (4.6) Let 𝜎^𝐷_Δ be the following smallest singular exponent in the Laplace problem with homogenous Dirichlet boundary condition: 𝜙∈𝐻¹₀(Ω)∶Δ𝜙∈𝐿₂(Ω)⊂∩_𝑠<𝜎^𝐷_Δ𝐻^𝑠(Ω),𝜙∈𝐻¹₀(Ω)∶Δ𝜙∈𝐿₂(Ω)⊂^¨¨𝐻^𝜎𝐷Δ(Ω). (4.7) From the regularity estimate, we know 𝜎^𝐷_Δ∈((3/2),2). Let 𝑟_min=2−𝜎^𝐷_Δ.

From [11, 25], we see that for all 𝑟∈(𝑟_min,1), the seminorm ‖^→𝑞‖_𝜒𝑟 is a norm on 𝜒_𝑟, and 𝜒_𝑟∩𝐻¹(Ω)³ is dense in 𝜒_𝑟.

In the following discussion, we will use 𝜒_𝑟,𝐿²_𝑟(Ω) for both nonconvex and convex domains. We select 𝑟∈(𝑟_min,1) for non-convex domain and 𝜒_𝑟=𝜒,𝐿²_𝑟(Ω)=𝐿²(Ω) for convex domain.

By introducing the Lagrange multiplier 𝜎, [6] changed (4.1) into the mixed variational formulation: find (𝜆,^→𝑢,𝜎)∈𝐑⁺×𝜒_𝑟×𝐿²_𝑟(Ω) such that ^→𝑢,^→𝜓_𝜒𝑟+div^→𝜓,𝜎_𝐿2𝑟=𝜆^→𝑢,^→𝜓₀,∀^→𝜓∈𝜒_𝑟,div^→𝑢,𝑣_𝐿2𝑟=0,∀𝑣∈𝐿²_𝑟(Ω). (4.8)

Let 𝐾^ℎ be a regular simplex partition (tetrahedral partition) of Ω with the mesh diameter ℎ. Define the finite element space 𝑉_ℎ×𝑊_ℎ⊂𝜒_𝑟×𝐿²_𝑟(Ω).

Restricting (4.8) on the above-mentioned finite space, we obtain the discrete mixed variational form: find (𝜆_ℎ,^→𝑢_ℎ,𝜎_ℎ)∈𝐑⁺×𝑉_ℎ×𝑊_ℎ such that ^→𝑢_ℎ,^→𝜓_𝜒𝑟+div^→𝜓,𝜎_ℎ𝐿2𝑟=𝜆_ℎ^→𝑢_ℎ,^→𝜓₀,∀^→𝜓∈𝑉_ℎ,div^→𝑢_ℎ,𝑣_𝐿2𝑟=0,∀𝑣∈𝑊_ℎ. (4.9)

Set 𝑉=𝜒_𝑟,‖⋅‖_𝑉=‖⋅‖_𝑎=‖⋅‖_𝜒𝑟,𝑊=𝐿²_𝑟(Ω),‖⋅‖_𝑊=‖⋅‖_𝐿²_𝑟,𝐷=𝐿₂(Ω)³,‖⋅‖_𝐷=‖⋅‖₀,𝑎^→𝑞,^→𝜓=^→𝑞,^→𝜓_𝜒𝑟𝑏^→𝜓,𝑣=div^→𝜓,𝑣_𝐿2𝑟. (4.10) Then (4.8) and (4.9) can be written as (2.3)-(2.4) and (2.5)-(2.6), respectively (it is needed to add  →  for the vector function, e.g., 𝑢,𝜓 should be written in the forms of ^→𝑢,^→𝜓).

We apply Schemes 1 and 2 to the eigenvalue problem of electric field (4.8). Adding the symbol  →  for the vector function, we get a multi-grid discretization scheme and adaptive algorithm for mixed finite element of the eigenvalue problem of electric field which are still called Schemes 1 and 2.

It is easy to know that 𝑎(⋅,⋅) and 𝑏(⋅,⋅) are continuous bilinear forms on 𝑉×𝑉 and 𝑉×𝑊, respectively. 𝑉↪^𝑐𝐷. (It is true obviously when Ω is convex; when Ω is non-convex, see [25].)

Consider the source problem associated with (4.8): find (𝑤^→,𝑝)∈𝜒_𝑟×𝐿²_𝑟(Ω) such that 𝑤^→,^→𝜓_𝜒𝑟+div^→𝜓,𝑝_𝐿2𝑟=^→𝑓,^→𝜓₀,∀^→𝜓∈𝜒_𝑟,div𝑤^→,𝑣_𝐿2𝑟=0,∀𝑣∈𝐿²_𝑟(Ω). (4.11)

For (4.11) and its conforming finite element approximation, condition (C1) of Lemma 2.1 holds obviously; [3] has proved that condition (C2) holds; assume that the discrete inf-sup condition (C3) holds, then conditions of Lemma 2.1 hold. Thus we can define operators 𝑇,𝑆,𝑇_ℎ, and 𝑆_ℎ. (4.8) and (4.9) can be written as (2.15) and (2.16)-(2.17), respectively.

The following Lemma 4.1 is cited from [3, 25].

Lemma 4.1. Equation (4.1) is equivalent to (4.8), and the solutions of (4.8) 𝜎=𝑆(𝜆^→𝑢)=0 and ^→𝑢∈𝜒_𝑟, div^→𝑢=0.

Note that Lemma 4.1 shows that for the eigenvalue problem of electric field, the second term on the right-hand side in (3.8) is equal to 0.

Theorem 4.2. Suppose that the discrete inf-sup condition (C3) holds. Then ‖𝑇−𝑇_ℎ‖_𝐷→0(ℎ→0). Let (𝜆_ℎ,^→𝑢_ℎ,𝜎_ℎ) be the 𝑘th eigenpair of (4.9) with ‖^→𝑢_ℎ‖_𝜒𝑟=1, and let 𝜆 be the 𝑘th eigenvalue of (4.8). Then 𝜆_ℎ→𝜆(ℎ→0), and there exists an eigenfunction ^→𝑢 corresponding to 𝜆 with ‖^→𝑢‖_𝜒𝑟=1, such that 𝜆_ℎ−𝜆≤𝜆_ℎ+𝜆𝜆_ℎ𝐶²₄𝐶²₂𝑇_ℎ−𝑇|_𝑀(𝜆)²_𝜒𝑟,_→𝑢−^→𝑢_ℎ𝜒𝑟≤𝐶₂𝑇_ℎ−𝑇|_{𝑀(𝜆)𝜒𝑟}, (4.12)(4.13) let ^→𝑢∈𝑀(𝜆), then there exists ^→𝑢_ℎ∈𝑀_ℎ(𝜆) such that ^→𝑢−^→𝑢_ℎ𝜒𝑟≤𝐶₃𝑇_ℎ−𝑇|_{𝑀(𝜆)𝜒𝑟}, (4.14) where 𝐶₂, 𝐶₃, and 𝐶₄ are constants independent of ℎ.

Proof. From the preceding discussion and hypotheses of the theorem, we know that 𝑉↪^𝑐𝐷, 𝑎(𝑢,𝑣) is symmetric, and the conditions of Lemma 2.1 hold. Besides, since 𝐾^ℎ is a regular partition, when Ω is a convex polyhedron, 𝜒_𝑟⊂𝐻¹(Ω)³(𝜒_𝑟=𝜒); when Ω is a non-convex polyhedron, 𝜒_𝑟∩𝐻¹(Ω)³ is dense in 𝜒_𝑟. Since 𝐶^∞(Ω)³ is dense in 𝐻¹(Ω)³, thus, no matter Ω is convex or non-convex, 𝐶^∞(Ω)³ is dense in 𝜒_𝑟. For any given ^→𝑓∈𝐷, we have 𝑇^→𝑓∈𝜒_𝑟. Thus for any 𝜀>0, according to the density, we know that there exists 𝑤^→∈𝐶^∞(Ω)³ such that 𝑇^→𝑓−𝑤^→_𝜒𝑟≤𝜀2. (4.15) Selecting ℎ₀>0 being small properly, when 0<ℎ≤ℎ₀, we have 𝑤^→−𝐼𝑤_ℎ^→_𝜒𝑟≤𝐶ℎ𝑤^→₂≤𝜀2, (4.16) where 𝐼_ℎ∶𝐶^∞(Ω)³→𝑉_ℎ is an interpolation operator. Thus inf_{→𝑞∈𝑉ℎ}𝑇^→𝑓−^→𝑞_𝜒𝑟≤𝑇^→𝑓−𝐼𝑤_ℎ^→_𝜒𝑟≤𝑇^→𝑓−𝑤^→_𝜒𝑟+𝑤^→−𝐼𝑤_ℎ^→_𝜒𝑟≤𝜀2+𝜀2=𝜀. (4.17) Namely, inf^→_{𝑞∈𝑉ℎ}‖𝑇^→𝑓−^→𝑞‖_𝜒𝑟→0(ℎ→0). Hence (2.26) is true. Analogously, using the density of 𝐶^∞(Ω) in 𝐿²_𝑟(Ω), we deduce that inf_{𝑣∈𝑊ℎ}‖𝑆^→𝑓−𝑣‖_𝐿²_𝑟→0(ℎ→0), namely, (2.27), is true. Hence, from Theorem 2.2 and Lemma 2.3, we know that ‖𝑇−𝑇_ℎ‖_𝐷→0(ℎ→0), 𝜆_ℎ→𝜆(ℎ→0), (4.13), and (4.14) hold. From (2.38), (4.13), and Lemma 4.1, we get (4.12). The proof is completed.

Denote 𝜀_𝜆(ℎ)=sup_→𝑢∈𝑀_(𝜆)inf_{→𝜓∈𝑉ℎ}^→𝑢−^→𝜓_𝜒𝑟. (4.18) From Lemma 2.1, noting 𝜎=0, we deduce 𝑇_ℎ−𝑇|_{𝑀(𝜆)𝜒𝑟}=sup_→𝑢∈𝑀₍^𝜆₎𝑇_ℎ−𝑇𝑢_𝜒𝑟≤sup_→𝑢∈𝑀₍^𝜆₎𝐶_𝑒inf_{→𝜓∈𝑉ℎ}𝑇^→𝑢−^→𝜓_𝜒𝑟≤𝜆⁻¹𝐶_𝑒𝜀_𝜆(ℎ)≡𝐶₇𝜀_𝜆(ℎ). (4.19)

Theorem 4.3. Assume that the discrete inf-sup condition (C3) holds, ℎ₀=𝐻 is small properly, Condition 1 holds and sup_𝑖𝑡_𝑖<3. Let (𝜆^ℎ𝑙,𝑢^ℎ𝑙) be an approximate eigenpair obtained by Scheme 1, then there exists ^→𝑢∈𝑀(𝜆) such that ^→𝑢^ℎ𝑙−^→𝑢_𝜒𝑟≤2𝐶₂×𝑞+3𝐶₇𝜀_𝜆ℎ_𝑙,𝜆^ℎ𝑙−𝜆≤4𝜆^ℎ𝑙+𝜆𝜆^ℎ𝑙𝐶²₄𝐶₂×𝑞+3²𝐶²₇𝜀_𝜆ℎ_𝑙². (4.20)(4.21)

Proof. We use induction to complete the proof. Note that the conditions of Theorem 2.2 hold. For 𝑙=1, Scheme 1 is actually two-grid discretization scheme. Substituting (4.19) into (4.12) and (4.13), we derive ^→𝑢_𝐻−^→𝑢_𝜒𝑟≤𝐶₂𝐶₇𝜀_𝜆(𝐻),𝜆−𝜆_𝐻≤𝜆_𝐻+𝜆𝜆_𝐻𝐶²₄𝐶²₂𝐶²₇𝜀_𝜆(𝐻)². (4.22) Combining (3.7) with 𝑙=1 and the above two inequalities, we know that there exists ^→𝑢∈𝑀(𝜆) such that 𝑢^ℎ1−^→𝑢_𝜒𝑟≤32𝜌𝐶₅𝐶₆𝜆^ℎ0−𝜆²+𝜆^ℎ0−𝜆𝑢^ℎ0−𝑢_𝐷+𝐶₂×𝑞+3𝑇−𝑇_ℎ1|_{𝑀(𝜆)𝜒𝑟}≤32𝜌𝐶₅𝐶₆𝜆_𝐻+𝜆𝜆_𝐻𝐶²₄²𝐶⁴₂𝐶⁴₇𝜀_𝜆(𝐻)⁴+𝜆_𝐻+𝜆𝜆_𝐻𝐶²₄𝐶²₂𝐶²₇𝜀_𝜆(𝐻)²𝐶₂𝐶₄𝐶₇𝜀_𝜆(𝐻)+𝐶₂×𝑞+3𝐶₇𝜀_𝜆ℎ₁≤2𝐶₂×𝑞+3𝐶₇𝜀_𝜆ℎ₁. (4.23) Since 𝜎=0, substituting (4.23) into (3.43), we deduce 𝜆^ℎ1−𝜆≤𝜆^ℎ1+𝜆𝜆^ℎ1𝐶²₄^→𝑢^ℎ1−^→𝑢²_𝜒𝑟≤𝜆^ℎ1+𝜆𝜆^ℎ1𝐶²₄4𝐶₂×𝑞+3²𝐶²₇𝜀_𝜆ℎ₁². (4.24) The above two inequalities show that Theorem 4.2 is true for 𝑙=1.
Suppose that the theorem is true for 𝑙−1, then by Theorem 3.3, we get ^→𝑢^ℎ𝑙−^→𝑢_𝜒𝑟≤32𝜌𝐶₅𝐶₆𝜆^ℎ𝑙−1+𝜆𝜆^ℎ𝑙−1𝐶²₄²4²𝐶₂×𝑞+3⁴𝐶⁴₂𝐶⁴₇𝜀_𝜆ℎ_𝑙−1⁴+𝜆^ℎ𝑙−1+𝜆𝜆^ℎ𝑙−1𝐶²₄4𝐶₂×𝑞+3²𝐶²₇𝜀_𝜆ℎ_𝑙−1²2𝐶₂×𝑞+3𝐶₇𝜀_𝜆ℎ_𝑙−1+𝐶₂×𝑞+3𝐶₇𝜀_𝜆ℎ_𝑙≤2𝐶₂×𝑞+3𝐶₇𝜀_𝜆ℎ_𝑙, (4.25) That is, (4.20) is valid.
From (4.20) and (3.8), we obtain (4.21). The proof is completed.

Assume that 𝐾^ℎ is a regular simplex partition (tetrahedral partition) of Ω with the mesh diameter ℎ. Let 𝑉_ℎ and 𝑊_ℎ be the 𝑃_𝑘+1-𝑃_𝑘 finite element spaces as follows: 𝑉_ℎ=^→𝑞∈𝐶⁰Ω³∶^→𝑞×^→𝛾|_𝜕Ω=0,^→𝑞|_𝜅∈𝑃_𝑘+1(𝜅)³⋅∀𝜅∈𝐾^ℎ,𝑊_ℎ=𝑣∈𝐶⁰Ω∶𝑣|_𝜅∈𝑃_𝑘(𝜅),∀𝜅∈K^ℎ,𝑣|_𝐸ℎ=0. (4.26) Here we set 𝐸_ℎ=∪_𝜅∈𝐾^ℎ_{,𝜕𝜅∩𝐸≠𝜙}𝜅. 𝑣|_𝐸ℎ=0 means that 𝑣 is equal to 0 on the tetrahedron where reentrant edge and angular point are adjacent. Considering finite element approximation of (4.11), for the 3-𝐷𝑃₂-iso-𝑃₁ Taylor-Hood finite element, Ciarlet and Girault [26] have discussed that the discrete inf-sup condition (C3) holds when Ω is a convex domain; for the 𝑃_𝑘+1-𝑃_𝑘 element, Ciarlet and Hechme have proved that the discrete inf-sup condition (C3) holds when Ω is a polyhedron (see Section  2.2 in [4] and pp. 509 in [3]).

From the above mentioned, we know that the 𝑃_𝑘+1-𝑃_𝑘 element approximation of (4.11) satisfies the conditions of Theorem 4.2.

Let 𝜎^𝑁_Δ be the smallest singular exponent in the Laplace problem with homogenous Neumman boundary condition, then 𝜎^𝑁_Δ∈(3/2,2). Denote 𝜏=min(𝑟−𝑟_min,𝜎^𝑁_Δ−1).

Corollary 4.4. Assume that ℎ₀=𝐻 is small properly, 𝑡_𝑖∈(1,3](𝑖=1,2,.,𝑙), and ℎ_𝑖=ℎ^𝑡𝑖_𝑖−1(𝑖=1,2,…,𝑙). Let (𝜆^ℎ𝑙,𝑢^ℎ𝑙) be an approximate eigenpair of the 𝑃_𝑘+1-𝑃_𝑘 element obtained by Scheme 1. Then when Ω is a convex domain, there exists ^→𝑢∈𝑀(𝜆) such that ^→𝑢^ℎ𝑙−^→𝑢_𝜒𝑟≤2𝐶₂×𝑞+3𝐶₇𝐶^′ℎ_𝑙,𝜆^ℎ𝑙−𝜆≤4𝜆^ℎ𝑙+𝜆𝜆^ℎ𝑙𝐶²₄𝐶₂×𝑞+3²𝐶²₇𝐶^′2ℎ²_𝑙. (4.27)(4.28) When Ω is a non-convex domain, there exists ^→𝑢∈𝑀(𝜆) such that ^→𝑢^ℎ𝑙−^→𝑢_𝜒𝑟≤2𝐶₂×𝑞+3𝐶₇𝐶^′′ℎ^𝜇_𝑙(),𝜆,∀𝜇∈0,𝜏_ℎ𝑙−𝜆≤4𝜆^ℎ𝑙+𝜆𝜆^ℎ𝑙𝐶²₄𝐶₂×𝑞+3²𝐶²₇𝐶^′′2ℎ_𝑙^2𝜇,∀𝜇∈(0,𝜏), (4.29)(4.30) where 𝐶^′ and 𝐶^′′ are determined by (4.31) and (4.32) in the proof, respectively.

Proof. The hypotheses of Corollary 4.4 imply that the conditions of Theorem 4.3 hold. When Ω is convex, for any ^→𝑢∈𝑀(𝜆) we have ^→𝑢=𝑇(𝜆^→𝑢)∈𝐻²(Ω) (see (44) in [1]). Thus there exists 𝐶^′ independent of ℎ_𝑙(𝑙=1,2,…) such that 𝜀_𝜆ℎ_𝑙≤𝐶^′ℎ_𝑙. (4.31) Substituting the above inequality into (4.20) and (4.21), we get (4.27) and (4.28), respectively.
When Ω is a non-convex domain, for any ^→𝑢∈𝑀(𝜆), by (36) in [3] we know that there exists 𝐶^′′ independent of ℎ_𝑙(𝑙=1,2,…), such that 𝜀_𝜆ℎ_𝑙≤𝐶^′′ℎ^𝜇_𝑙. (4.32) Substituting the above inequality into (4.20) and (4.21), we derive (4.29) and (4.30), respectively.

Remark 4.5. From Corollary 4.4, we see that the constants in the error estimates are not only independent of the mesh diameter but also independent of the iterative degrees. Thus, when 𝑙→∞, we have 𝜆^ℎ𝑙→𝜆. Suppose that the precision order of (4.28) is optimal which cannot be improved any more, then 𝜆^ℎ𝑙−𝜆^ℎ𝑙+1=𝜆^ℎ𝑙−𝜆+𝜆−𝜆^ℎ𝑙+1=𝑂ℎ_𝑙^2𝑟,𝜆^ℎ𝑙−1−𝜆^ℎ𝑙=𝜆^ℎ𝑙−1−𝜆+𝜆−𝜆^ℎ𝑙=𝑂ℎ^2𝑟_𝑙−1, (4.33) where for convex domain 𝑟=1, while for non-convex domain 𝑟=𝜇∈(0,𝜏) but approximates 𝜏 arbitrarily. Therefore we have that 𝜆^ℎ𝑙−𝜆^ℎ𝑙+1≤𝜆^ℎ𝑙−1−𝜆^ℎ𝑙. Then 𝜆^ℎ𝑙−𝜆≤𝜆^ℎ𝑙−𝜆^ℎ𝑙+1+𝜆^ℎ𝑙+1≤𝜆−𝜆_ℎ𝑙−1−𝜆^ℎ𝑙+𝜗ℎ^2𝑟_𝑙+1, (4.34) thus we can use 𝜂(𝜆^ℎ𝑙)=𝜃×|𝜆^ℎ𝑙−1−𝜆^ℎ𝑙| as a posteriori error indicator of |𝜆^ℎ𝑙−𝜆|, where 𝜃∈(0,1].

5. Numerical Experiments

Consider the eigenvalue problem (4.1) of electric field, where Ω=[0,𝜋]×[0,𝜋] is a square domain or Ω=[−1,0]×[−1,0]∪[−1,1]×[0,1] is an L-shaped domain. For the square domain, the first five exact eigenvalues are 𝜆₁=𝜆₂=1, 𝜆₃=2, and 𝜆₄=𝜆₅=4; for the L-shaped domain, the first five eigenvalues are 𝜆₁≈1.47562182408, 𝜆₂≈3.53403136678, 𝜆₃=𝜆₄=𝜋²≈9.86960440109, and 𝜆₅≈11.3894793979, and the first eigenfunction has a strong singularity (see, e.g., [3]).

We adopt a uniform isosceles right triangulation for Ω (the edge in each element is along three fixed directions, see Figure 1) to produce the meshes 𝐾^ℎ𝑙 with mesh diameter ℎ_𝑙.

Figure 1 

The definition of 𝑃₂-𝑃₁ mixed finite element spaces is given by 𝑉_ℎ𝑙=^→𝑞∈𝐶⁰Ω²∶^→𝑞×^→𝛾|_𝜕Ω=0,^→𝑞|_𝜅∈𝑃₂(𝜅)²,∀𝜅∈𝐾^ℎ𝑙,𝑊_ℎ𝑙=𝑣∈𝐶⁰Ω∶𝑣|_𝜅∈𝑃₁(𝜅),∀𝜅∈𝐾^ℎ𝑙,𝑣|_𝐸ℎ𝑙=0. (5.1) We make use of Matlab to compute the first five approximate eigenvalues by using Scheme 1 with 𝑃₂-𝑃₁ element. The numerical results are listed in Tables 1, 2, and 3.

Table 1 

The results on the square by Scheme 1 (𝑃₂-𝑃₁ element) for the eigenvalue problem of electric field (𝑟=0).

Table 2 

The results on the L-shaped domain by Scheme 1 (𝑃₂-𝑃₁ element) for the eigenvalue problem of electric field (𝑟=0.5).

Table 3 

The results on the L-shaped domain by Scheme 1 (𝑃₂-𝑃₁ element) for the eigenvalue problem of electric field (𝑟=0.95).

In Tables 1–3, 𝜆_1,ℎ𝑙,𝜆_2,ℎ𝑙,…,𝜆_5,ℎ𝑙 denote the first five eigenvalues obtained by using the mixed element on 𝐾^ℎ𝑙 directly; 𝜆^ℎ𝑙₁,𝜆^ℎ𝑙₂,…,𝜆^ℎ𝑙₅ denote the first five eigenvalues obtained by Scheme 1: 𝑒𝜆_𝑘,ℎ𝑙=𝜆_𝑘,ℎ𝑙−𝜆_𝑘,𝑒𝜆^ℎ𝑙_𝑘=𝜆^ℎ𝑙_𝑘−𝜆_𝑘,𝜂𝜆_𝑘,ℎ𝑙=𝜃𝜆_𝑘,ℎ𝑙−𝜆_{𝑘,ℎ𝑙−1},𝜂𝜆^ℎ𝑙_𝑘=𝜃𝜆^ℎ𝑙_𝑘−𝜆^ℎ𝑙−1_𝑘,(𝜃=1). (5.2)

From Tables 1–3, we can see that (1)  𝜆^ℎ𝑙_𝑘 and 𝜆_𝑘,ℎ𝑙 have the same accuracy. (2) It fails to find 𝜆_𝑘,ℎ𝑙 by direct computation by using the mixed element on 𝐾^ℎ𝑙 with ℎ_𝑙=√2/128 in the case of the square domain and ℎ_𝑙=√2/160 in the case of the L-shaped domain (here Matlab shows that Sparse lu with 4 outputs (UMFPACK) failed), but Scheme 1 still works. (3)𝜂(𝜆^ℎ𝑙_𝑘)=|𝜆^ℎ𝑙_𝑘−𝜆^ℎ𝑙−1_𝑘| is an efficient and reliable error indicator of 𝜆^ℎ𝑙_𝑘.

It can be seen from Tables 1–3 that, in the calculation of error indicators, 𝜃 should be selected as follows: in the case of the square domain, 𝜃=1/15; in the case of the L-shaped domain, 𝜃 is equal to 1,3/5,1/15, and 3/5, respectively when 𝑟=0.5, and 𝜃 is equal to 1,1/4,1/15 and 1/4 respectively, when 𝑟=0.95.

Remark 5.1. Taking Table 1, for example, we illustrate how to select 𝑡_𝑖 next. In Table 1, ℎ₀=√2/8, ℎ₁=√2/32, ℎ₂=√2/64, and ℎ₃=√2/128. According to ℎ_𝑖=ℎ^𝑡𝑖_𝑖−1, by calculation, we get that 𝑡₁≈1.80, 𝑡₂≈1.22, and 𝑡₃≈1.18.

Acknowledgments

The authors cordially thank the editor and the referees for their valuable comments which lead to the great improvement of this paper. This work is supported by National Natural Science Foundation of China (Grant no. 11161012).