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

This paper discusses highly ﬁnite element algorithms for the eigenvalue problem of electric ﬁeld. Combining the mixed ﬁnite 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 ﬁeld. Theoretical analysis and numerical results show that the computational schemes established in the paper have high e ﬃ ciency.


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 θ × |λ h l − λ h l−1 | θ ∈ 0, 1 as a posteriori error estimator of λ h l 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.

Preliminaries
Let V , W, and D be three real Hilbert spaces with inner products and norms ·, · V , · V , ·, · W , · W , ·, · D , and · D , respectively. Suppose that V → D continuously imbedded , a ·, · is a symmetric, continuous, and V -elliptic bilinear form on V × V , that is, It is obvious that a ·, · is an inner product on V and · a a ·, · is an equivalent norm to · V . Abstract and Applied Analysis 3 In scientific and engineering computations, many eigenvalue problems have the following first type of mixed variational formulation: find λ, u, σ ∈ R × V × W, u, σ / 0, 0 , such that In order to solve problem 2.
Consider the associated source and approximate source problems.

2.8
Note that the source term f 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.
then there exists a unique solution w h , p h to the problem 2.8 and the following error estimate is valid: where C e just depends on ν, ν 2 and M 1 , M 2 .
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: It is easy to verify that T : D → D, T h : D → D are self-adjoint operators in the sense of inner product ·, · D see 7 .

2.13
Abstract and Applied Analysis 5 Assume V c → D compactly embedded , then it's easy to prove that T : D → D is completely continuous, T : V → V is completely continuous, and T h is a finite rank operator. Combining 2.3 -2.4 , 2.5 -2.6 , and the V -ellipticity of a ·, · , we deduce 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 and the corresponding eigenfunctions are The eigenvalues of 2.5 -2.6 can be sorted as Let μ be the kth eigenvalue with algebraic multiplicity q, μ μ k μ k 1 · · · μ k q−1 .
M μ is the space spanned by all eigenfunctions {u j }

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.

2.27
Then Proof. From V c → D, we derive that T : D → D is a completely continuous operator. It is obvious that T h : D → D is a finite rank operator. From 2.12 , 2.26 , and 2.27 , we deduce It shows that T h : D → V pointwisely converges to T .
where u depends on h in general, and C 1 , C 2 , and C 3 are constants independent of h.
Thus from Theorem 2.2 in 7 , we see that the desired results are valid. The proof is completed.

2.35
Proof. The proof is completed by using the same argument as that of Lemma 9.1 see 7, 18 . 2.36 Taking u * , σ * u h , σ h in 2.35 and using 2.4 and 2.6 , we deduce the following.

8
Abstract and Applied Analysis

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 , − λ k 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.

3.2
Set u h i u / u a , σ h i σ / u a .
Step 4. Compute the Rayleigh quotient

3.3
Step 5. If i l, then output λ h l , u h l , σ h l , stop; else, i ⇐ i 1, and return to Step 3.
Let λ H , u H , σ H be the kth eigenpair, and we use λ h l , u h l , σ h l as the kth approximation eigenpair of 2.3 -2.4 .
Next, we will discuss the validity of Scheme 1.
Proof. See 8 Since the convergence rate of V h l−1 and W h l−1 approximating eigenfunctions is lower than that of V h l and W h l approximating eigenfunctions, respectively, the approximation order of λ h l−1 , u h l−1 is lower than that of λ h l , u h l . However, in general, the accuracy order of λ h l−1 , u h l−1 will not exceed that of λ h l−1 , u h l−1 ; therefore in the following Theorem 3.3 we assume that the accuracy order of λ h l−1 , u h l−1 is lower than that of λ h l , u h l .
, H is small properly, and Condition 1 holds. Let λ h l , u h l , σ h l be the approximate eigenpair obtained by Scheme 1, and let u h l−1 approximate u ∈ M λ , λ h l−1 approximate λ, and the accuracy order of λ h l−1 , u h l−1 be lower than that of λ h l , u h l . Then there exists u ∈ M λ such that where C 5 and C 6 are determined, respectively, by 3.10 and 3.14 in the following proof.

3.9
From Lemma 2.1, there exists a positive constant C 5 depending only on ν, ν 1 , ν 2 , M 1 , and M 2 such that 3.10 Then

3.12
Using the triangle inequality and 2.33 , we deduce According to the hypotheses of the theorem, we know that λ h l−1 → λ and λ h l−1 − λ are an infinitesimal of lower order comparing with λ j,h l − λ. Hence, there exists a positive constant C 6 independent of h l l 1, 2, . . . such that for j k, k 1, . . . , k q − 1 we have 3.14 Note that H is small enough and h l h l−1 ; from 3.13 and 3.12 , we obtain Noticing that λ λ k λ k 1 · · · λ k q−1 , we have which together with 3.14 , noting that λ j,h l − λ is an infinitesimal of higher order comparing with λ h l−1 − λ, yields Since ρ is the separation constant, H is small enough, and h l h l−1 , there holds For u in Step 3 of Scheme 1, from the definition of T h and S h taking i l , we have Hence, when i l, Step 3 of Scheme 1 is equivalent to the following: find u , σ ∈ V h l × W h l such that And set u h l u / u a , σ h l σ / u a . From 3.23 , we obtain By 3.26 and taking ψ u − λ h l−1 T h l u − T h l u h l−1 in 3.25 , we obtain From 3.28 we know that the first term on the left-hand side of 3.25 is equal to 0, thus then, using discrete inf-sup condition we get σ λ h l−1 S h l u S h l u h l−1 .

3.30
Thus Step 3 of Scheme 1 is equivalent to 3.28 , 3.30 , u h l u / u a , and σ h l σ / u a . Noting that 1/λ h l−1 T h l u h l−1 1/λ h l−1 T h l u h l−1 a u 0 differs from u 0 by only a constant and selecting u s λ h l−1 u / T h l u h l−1 a , we have

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  Combining 3.35 with the above inequality, we have

3.38
Substituting 3.12 into 3.38 , we get 3.7 . From Step 3 of Scheme 1, we know that b u h l , σ h l 0, thus

3.39
Select λ r λ h l , u * u h l , and σ * σ h l . From Lemma 2.4, we get

3.40
Noting that ∀v ∈ W h l , b u h l − u, v 0, and u h l , u h l D 1/λ h l a u h l , u h l 1/λ h l , we have

3.41
Note u h l − u D ≤ C 4 u h l − u a . From 3.41 we obtain 3.8 .
For l 1, Scheme 1 is actually the two-grid discretization scheme established in 7 . By Theorem 3.3 we get the following conclusion directly.
be an approximate eigenpair obtained by Scheme 1 (l 1). Then there exists u ∈ M λ such that

3.42
Proof. Consider Scheme 1. Here l 1, u h 0 u H . By Lemma 2.3, we know that u h 0 approximates u ∈ M λ , and the accuracy order of λ h 0 , u h 0 is lower than λ h 1 , u h 1 . Hence, for l 1, the conditions of Theorem 3.3 hold. Select u in the proof of Theorem 3.3 such that u h 0 − u satisfies Lemma 2.3. Then from 2.30 , we obtain 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 λ h l → λ l → ∞ and the convergence rate is high. Thus, we use θ × |λ h l − λ h l−1 | as a posteriori error indicator of λ h l − λ 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, H, t 1 , and h 1 H t 1 . Step

3.46
Let u h l u / u a , σ h l σ / u a .
Step 4. Compute the Rayleigh quotient λ h l a u h l , u h l u h l , u h l D .

Remark 3.5.
In Scheme 2, we use θ × |λ h l − λ h l−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.

The Eigenvalue Problem of Electric Field
Consider the eigenvalue problem of electric field: where Ω is a polyhedron in R 3 , γ is the unit outward normal to ∂Ω.
Physically u denotes the electric field, ω denotes the time frequency, and c is the speed of the light velocity. Usually, let λ ω 2 /c 2 named eigenvalue. Let

16 Abstract and Applied Analysis
When Ω is a convex polyhedron, we define the following function space:
When Ω is a nonconvex polyhedron, the problem is relatively complicated. Let E denote a set of edges of reentrant dihedral angles on ∂Ω, and d d x denote the distance to the set E: d x dist x, ∪ e∈E e . We introduce a weight function ω r which is a nonnegative smooth function of x. It can be represented by d r in reentrant edge and angular domain. We shall write ω r d r . Define the weighted functional spaces: q, ψ χ r curl q, curl ψ 0 div q, div ψ L 2 r , q χ r q, q 1/2 χ r .

4.6
Let σ D Δ be the following smallest singular exponent in the Laplace problem with homogenous Dirichlet boundary condition:

4.7
From the regularity estimate, we know σ D Δ ∈ 3/2 , 2 . Let r min 2 − σ D Δ . From 11, 25 , we see that for all r ∈ r min , 1 , the seminorm q χ r is a norm on χ r , and In the following discussion, we will use χ r , L 2 r Ω for both nonconvex and convex domains. We select r ∈ r min , 1 for non-convex domain and χ r χ, L 2 r Ω L 2 Ω for convex domain.
By introducing the Lagrange multiplier σ, 6 changed 4.1 into the mixed variational formulation: find λ, u, σ ∈ R × χ r × L 2 r Ω such that

4.8
Let K h be a regular simplex partition tetrahedral partition of Ω with the mesh diameter h. Define the finite element space V h × W h ⊂ χ r × L 2 r Ω . Restricting 4.8 on the above-mentioned finite space, we obtain the discrete mixed variational form:

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., u, ψ should be written in the forms of u, ψ . 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 a ·, · and b ·, · are continuous bilinear forms on V × V and It is true obviously when Ω is convex; when Ω is non-convex, see 25 .
Consider the source problem associated with 4.8 : find w, p ∈ χ r × L 2 r Ω such that

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 T, S, T h , and S h . 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 . 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
Let λ h , u h , σ h be the kth eigenpair of 4.9 with u h χ r 1, and let λ be the kth eigenvalue of 4.8 . Then λ h → λ h → 0 , and there exists an eigenfunction u corresponding to λ with u χ r 1, such that where C 2 , C 3 , and C 4 are constants independent of h.
Proof. From the preceding discussion and hypotheses of the theorem, we know that V c → D, a u, v is symmetric, and the conditions of Lemma 2.1 hold. Besides, since K h is a regular partition, when Ω is a convex polyhedron, χ r ⊂ H 1 Ω 3 χ r χ ; when Ω is a non-convex polyhedron, χ r ∩ H 1 Ω 3 is dense in χ r . Since C ∞ Ω 3 is dense in H 1 Ω 3 , thus, no matter Ω is convex or non-convex, C ∞ Ω 3 is dense in χ r . For any given f ∈ D, we have T f ∈ χ r . Thus for any ε > 0, according to the density, we know that there exists w ∈ C ∞ Ω 3 such that Selecting h 0 > 0 being small properly, when 0 < h ≤ h 0 , we have

4.17
Namely, inf q∈V h T f − q χ r → 0 h → 0 . Hence 2.26 is true. Analogously, using the density

4.19
Theorem 4.3. Assume that the discrete inf-sup condition (C3) holds, h 0 H is small properly, Condition 1 holds and sup i t i < 3. Let λ h l , u h l be an approximate eigenpair obtained by Scheme 1, then there exists u ∈ M λ such that Proof. We use induction to complete the proof. Note that the conditions of Theorem 2.2 hold. For l 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 l 1 and the above two inequalities, we know that there exists u ∈ M λ such that Abstract and Applied Analysis

4.24
The above two inequalities show that Theorem 4.2 is true for l 1.
Suppose that the theorem is true for l − 1, then by Theorem 3.3, we get That is, 4.20 is valid. From 4.20 and 3.8 , we obtain 4.21 . The proof is completed.
Assume that K h is a regular simplex partition tetrahedral partition of Ω with the mesh diameter h. Let V h and W h be the P k 1 -P k finite element spaces as follows:

4.26
Here we set E h ∪ κ∈K h ,∂κ∩E / φ κ. v| E h 0 means that v 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-DP 2 -iso-P 1 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 P k 1 -P k 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 P k 1 -P k element approximation of 4.11 satisfies the conditions of Theorem 4.2.

4.28
When Ω is a non-convex domain, there exists u ∈ M λ such that When Ω is a non-convex domain, for any u ∈ M λ , by 36 in 3 we know that there exists C independent of h l l 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 l → ∞, we have λ h l → λ. Suppose that the precision order of 4.28 is optimal which cannot be improved any more, then Abstract and Applied Analysis 2 Figure 1 where for convex domain r 1, while for non-convex domain r μ ∈ 0, τ but approximates τ arbitrarily. Therefore we have that λ h l − λ h l 1 ≤ λ h l−1 − λ h l . Then
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 K h l with mesh diameter h l .

5.1
We make use of Matlab to compute the first five approximate eigenvalues by using Scheme 1 with P 2 -P 1 element. The numerical results are listed in Tables 1, 2, and 3.