Multigrid methods based on shifted inverse iteration for the Maxwell eigenvalue problem
Jiayu Han

 Discovery of optical flickering from the symbiotic star EF Aquilae


    R. K. Zamanov1Corresponding authors:
        rkz@astro.bas.bg,  kstoyanov@astro.bas.bg


S. Boeva 1
Y. M. Nikolov1
B. Petrov1
R. Bachev1
G. Y. Latev1
V. A. Popov2
K. A. Stoyanov1
M. F. Bode3
J. Martí4
T. Tomov5  
A. Antonova6

    December 30, 2023
======================================================================================================================================================================================================================================

In this paper  two types of multgrid methods, i.e., the Rayleigh quotient iteration
and the inverse iteration with fixed shift,  are developed
for solving the Maxwell eigenvalue problem with  discontinuous relative magnetic permeability
and electric permittivity. With the aid of the mixed form of source problem associated with the eigenvalue problem, we prove
 the uniform convergence of
 the discrete solution operator  to the solution operator in ł2 using discrete compactness of edge element space.
Then we prove the asymptotically optimal error estimates for both multigrid methods. Numerical experiments
confirm our theoretical analysis.


Maxwell
eigenvalue problem, multigrid method, edge element, error analysis


65N25, 65N30

myheadings plain JIAYU HAN MULTIGRID METHODS FOR MAXWELL EIGENVALUES


§ INTRODUCTION

The Maxwell eigenvalue problem is of basic importance in designing resonant structures for advanced waveguide.
Up to now, the communities from  numerical mathematics and  computational electromagnetism have developed plenty of
numerical methods  for solving this problem (see, e.g., <cit.>).

The difficulty of numerically solving the  eigenvalue problem lies in imposing the divergence-free constraint.
For this purpose, the nodal finite element methods utilize the filter, parameterized and mixed approaches to
find the true eigenvalues <cit.>.
The researchers in  electromagnetic field usually   adopt edge finite elements due to   the property of
tangential continuity of electric field <cit.>.
 Using edge finite element methods, when one only considers to compute the nonzero eigenvalue, the divergence-free
 constraint can be dropped from the weak form and satisfied naturally (see <cit.>).
 But this will introduce spurious zero eigenvalues.
Since the eigenspace  corresponding to zero  is infinite-dimensional,
usually, the finer the mesh is, the more the spurious eigenvalues there are.
However, using this form there is no difficulty in computing   eigenvalues on a very coarse mesh.
So the work of <cit.> subtly applies this weak form to two grid method for the Maxwell eigenvalue problem.
That is, they first solve a Maxwell eigenvalue problem on a coarse mesh and then solve a linear Maxwell
equation on a fine mesh.
  Another approach is the mixed form of saddle point type in which a Lagrange multiplier is introduced to impose
  the divergence-free constraint
  (see <cit.>).
  A remarkable feature of the mixed form is
   its equivalence to the weak form in <cit.> for nonzero eigenvalues.
   The mixed form is well known as having a good property of no spurious eigenvalues being introduced. However,
   it is not an easy task to   solve it on a fine mesh (see <cit.>).

The multigrid methods for solving eigenvalue problems originated from the idea of two grid method proposed by <cit.>.
Afterwards, this work was further developed by <cit.>. Among them the more recent work   <cit.>
makes a relatively systematical research on  multigrid methods based on shifted inverse iteration especially on
its  adaptive fashion.

Inspired by the above works, this paper is devoted to developing  multigrid methods for solving the Maxwell eigenvalue problem. We first use the mixed form to solve the eigenvalue problem on a coarser mesh and then
solve a series of Maxwell equations on finer and finer meshes without using the mixed form. Roughly speaking,
we develop the two grid method in <cit.> into multigrid method where only  nonzero eigenvalues are focused on.
We prefer to use the mixed form instead of the one in <cit.> on a coarse mesh to capture the information of the true
 eigenvalues.  One reason is that  the mixed form can include the physical zero eigenvalues and rule out the spurious zero ones simultaneously.
 Using it, the physical zero eigenvalues can be captured on a very coarse mesh,   
which is necessary when  the resonant cavity has   disconnected boundaries(see, e.g., <cit.>).
Another reason lies in that the mixed discretization of saddle point type is  not difficult to solve on a coarse mesh.

In this paper, we study two types of multigrid methods based on shifted inverse iteration:
Rayleigh quotient iteration and inverse iteration with fixed shift. The former  is a well-known method
 for solving matrix eigenvalues but the corresponding coefficient matrix is nearly singular and
 difficult to solve to some extend. To overcome this difficulty the latter first performs the
 Rayleigh quotient iteration at previous few steps and then fixes the shift at the following steps as the estimated eigenvalue
 obtained by the former.
 Referring to the error analysis framework in <cit.> and using compactness property of edge element space,
  we first prove the  uniform convergence of
 the discrete solution operator  to the solution operator in ł2 and then the error estimates of eigenvalues and eigenfunctions
 for the mixed discretization;
 then we adopt  the analysis tool in <cit.> that is different from the one in <cit.> and  prove the  asymptotically optimal error estimates for both
 multigrid methods.
     In addition, this paper is concerned about the  theoretical analysis   for the case of the
     discontinuous electric permittivity  μ and magnetic permeability ϵ in complex matrix form, which has important applications
for the resonant cavity being filled with different dielectric materials invariably.
It is noticed that our multigrid methods and  theoretical results are not only  valid for the lowest order edge elements
but also for high order ones.

More importantly, based on the work of this paper, once given an a posteriori error indicator
of eigenpair one can further develop the adaptive algorithms
of shifted inverse iteration type for the problem.
In the last
section of this paper, we present several numerical examples to validate the efficiency of our methods in different
cases.

Throughout this paper,   we use
the symbol a ≲ b   to mean that a ≤ Cb, where C denotes a positive constant independent of   mesh parameters and iterative times
and may not be the same   in different places.



§ PRELIMINARIES

Consider the Maxwell eigenvalue problem in electric field


    curl(μ^-1curl𝐮)=ω^2ϵ𝐮   in Ω,
     div(ϵ𝐮)=0    in Ω,
    𝐮_t=0    in ∂Ω,

where Ω is a bounded Lipschitz polyhedron domain in ℝ^d(d=2,3), 𝐮_t is the tangential
trace of 𝐮.
The coefficient μ is the electric permittivity, and ϵ is the magnetic permeability and  piecewise smooth.


    In this paper, λ=ω^2 with ω being the angular frequency, is defined as the eigenvalue of
     this problem. We assume that μ,ϵ are two positive definite
      Hermite matrices such that μ^-1,ϵ∈ (L^∞(Ω))^d× d
   and there exist  two
positive numbers γ,β satifying
 
    ξ·μ^-1ξ≥γξ·ξ, ξ·ϵξ≥βξ·ξ, ∀ 0≠ξ∈ℂ^d.






 §.§ Some weak forms

Let

    𝐇_0(curl,Ω)={𝐮∈𝐋^2(Ω): curl(𝐮)∈𝐋^2(Ω),𝐮_t|_∂Ω=0},

equipped with the norm 𝐮_curl:=curl𝐮_0+𝐮_0.
Throughout this paper,
·_0 and ·_0,ϵ denote the norms in ł2 induced by the inner products (·,·)
and (ϵ·,·)    respectively.
Define the  divergence-free space:

    𝐗:={𝐮∈𝐇_0(curl,Ω): div(ϵ𝐮)=0}.






The standard weak form of the Maxwell eigenvalue problem (<ref>)-(<ref>) is as follows:
Find (λ,𝐮)∈ℝ×𝐗 and 𝐮≠0 such that

    a(𝐮,𝐯)=λ  (ϵ𝐮,𝐯),   ∀𝐯∈𝐗,

where a(𝐮,𝐯)=(μ^-1curl𝐮,curl𝐯).
Denote v̱_a:=√(a(v̱,v̱)),∀v̱∈0̧.

As the divergence-free space  in (<ref>) is difficult to  discretize, alternatively, we would like to solve
the eigenvalue problem (<ref>)-(<ref>) in the larger space 𝐇_0(curl;Ω), that is: 
 Find (λ,𝐮)∈ℝ×𝐇_0(curl;Ω) and 𝐮≠0 such that

    a(𝐮,𝐯)=λ  (ϵ𝐮,𝐯),   ∀𝐯∈𝐇_0(curl;Ω).

Note that when λ≠0 (<ref>) and (<ref>) are equivalent, since (<ref>) implies the divergence-free
condition holds for λ≠0 (e.g., see <cit.>).

According to (<ref>), we have
  
    √(γ)curlu̱_0≤u̱_a

In order to study the  eigenvalue problem in 0̧ we need the auxiliary bilinear form

    A(w̱,v̱)=a(w̱,v̱)+ γ/β(ϵw̱,v̱),

which defines an equivalent norm ·_A=√(A(·,·)) in 0̧.














 By Lax-Milgram Theorem we can define the solution operator T:ł2→𝐗 as

    A(T𝐟,𝐯)=   (ϵ𝐟,𝐯), ∀𝐯∈𝐗.



Then  the eigenvalue problem (<ref>)   has the operator form

    Tu̱=λ^-1u̱ with λ=λ+γ/β.








The following mixed weak form of saddle point type   can be found in <cit.>:
Find (λ,𝐮,σ)∈ℝ×0̧× H^1_0(Ω) with
𝐮0  such that

    a(𝐮,𝐯)+b(𝐯,σ)   =   λ(ϵ𝐮,𝐯),  ∀𝐯∈0̧,
    
    b(𝐮,p)   =   0,  ∀ p∈ H^1_0(Ω),

where b(𝐯, p):=(ϵv̱,∇ p) for any v̱∈0̧, p∈ H^1_0(Ω).


We introduce the corresponding mixed equation: Find Tf̱∈0̧ and S f̱∈  H^1_0(Ω) for f̱∈ł2 such that

    A(Tf̱,𝐯)+b(𝐯, S f̱)   =    (ϵ𝐟,𝐯),  ∀𝐯∈0̧,
    
    b(T𝐟,p)   =   0,  ∀ p∈  H^1_0(Ω);



the following LBB condition
can be verified by taking w̱=∇v̱

    sup_w̱∈0̧|b(w̱,v̱)|/w̱_curl≥β |v̱|_1, ∀v̱∈ H^1_0(Ω).

This yields the existence and uniqueness of linear bounded operators T and S (see <cit.>).
Due to Helmholtz decomposition 𝐇_0(curl;Ω)=∇H_0^1(Ω)⊕𝐗, it is easy to
see R(T)⊂ and S f̱=0, Tf̱=Tf̱  for any f̱∈. Hence  T and T
share the same eigenpairs. More importantly,  the operator T: ł2ł2 is self-adjoint. In fact,
∀𝐯,𝐰∈ł2,

    (ϵw̱, Tv̱)=A( T𝐰, T𝐯)
    =A( T𝐯, T𝐰)
    =(ϵ𝐯, T𝐰)= (ϵ T𝐰,𝐯).

Note that T is compact as a operator from ł2 to ł2 and from
𝐗 to 𝐗 since 𝐗↪ł2 compactly (see Corollary 4.3 in <cit.>).






 §.§ Edge element discretizations and error estimates

 We will consider the edge element approximations based on the  weak forms (<ref>),
 (<ref>) and (<ref>)-(<ref>). Let π_h be a shaped-regular triangulation of Ω composed of the elements κ.
Here we restrict our attention to  edge   elements on tetrahedra because the argument for edge   elements on hexahedra
is the same.
The k(k≥0) order edge element of the first family <cit.> generates the space

    𝐕_h={𝐮_h∈𝐇_0(curl,Ω):𝐮_h|_κ∈ [P_κ(k)]^d⊕𝐱× [P_κ(k)]^d},

where P_κ(k) is the polynomial space of  degree less than or equal to k on κ, P_κ(k) is the
homogeneous polynomial space of  degree k  on κ, and x̱ =(x_1,⋯,x_d)^T.

We also introduce the discrete divergence-free space

    𝐗_h={𝐮_h∈𝐕_h: (ϵ𝐮_h,∇ p)=0,  ∀ p ∈ U_h},

where U_h is the standard  Lagrangian finite element space vanishing on ∂Ω of total degree less than or
equal to k+1 and ∇ U_h⊂𝐕_h. 

The standard finite element discretization of (<ref>) is stated as: Find (λ_h,𝐮_h)∈ R×_h and 𝐮_h≠0 such that

    a(𝐮_h,𝐯_h)=λ_h  (ϵ𝐮_h,𝐯_h),  ∀𝐯_h∈_h.

It is also equivalent to the following form for  nonzero λ_h (see <cit.>):  
Find (λ_h,𝐮_h)∈ℝ×V̱_h and 𝐮_h≠0 such that

    a(𝐮_h,𝐯_h)=λ_h  (ϵ𝐮_h,𝐯_h),  ∀𝐯_h∈𝐕_h.


In order to investigate the convergence of edge element discretization   (<ref>), we have to study the convergence of edge element discretization for
the associated Maxwell source problem. 







 Then by Lax-Milgram Theorem we can define the solution operator T_h:ł2
→𝐗_h as

    A(T_h𝐟,𝐯)=   (ϵ𝐟,𝐯), ∀𝐯∈𝐗_h.

Then  the eigenvalue problem (<ref>)  has the operator form

    T_hu̱_h=λ_h^-1u̱_h with λ_h=λ_h+γ/β.

Introduce the  discrete form of (<ref>)-(<ref>):
Find (λ_h,𝐮_h,σ_h)∈ℝ×V̱_h× U_h,
𝐮_h0 such that

    a(𝐮_h,𝐯)+b(𝐯,σ_h)   =   λ_h(ϵ𝐮_h,𝐯),
      ∀𝐯∈V̱_h,
    
    b(𝐮_h,p)   =   0,  ∀ p∈ U_h.

Introduce the corresponding operators: Find T_hf̱∈V̱_h and S_hf̱∈   U_h for any f̱∈ł2

    A(T_hf̱,𝐯)+b(𝐯,S_hf̱)   =    (ϵ𝐟,𝐯),  ∀𝐯∈V̱_h,
    
    b(T_hf̱,p)   =   0,  ∀ p∈ U_h.

Due to discrete Helmholtz decomposition 𝐕_h=∇ U_h⊕𝐗_h, it is easy to know R(T_h)⊂_h and S_hf̱=0, T_hf̱=T_hf̱ for any f̱∈+_h.
Hence T_h and T_h share the same eigenpairs.






 Similar to (<ref>)-(<ref>), one can  verify the corresponding LBB condition for the discrete mixed  form  (<ref>)-(<ref>).
 According to the theory of mixed finite elements (see <cit.>), we get for all f̱∈ł2,

    T_h 𝐟_A+ T 𝐟_A+|S_hf̱|_1+|Sf̱|_1≤𝐟_0,ϵ,


    T 𝐟- T_h𝐟_A+|S f̱-S_h f̱|_1≤(inf_𝐯_h∈𝐕_hT𝐟-𝐯_h_curl+inf_v_h∈  U_h|S f̱-v_h|_1).





Similar to (<ref>) we can prove T_h:ł2ł2 is self-adjoint in the sense of (ϵ·,·)_0. In fact, ∀w̱,v̱∈ł2,

    (ϵw̱, T_hv̱)=A( T_h w̱, T_hv̱)=A( T_h v̱, T_hw̱)=
    (ϵv̱, T_hw̱)= (ϵ T_h w̱,v̱).









 The discrete compactness is a very interesting and important property in edge elements because it is intimately related to the property of the collective compactness. Kikuchi <cit.> first successfully applied this property to numerical analysis of electromagnetic problems,
 and more recently it was further developed by <cit.> and so on.
  The following lemma, which states the discrete compactness of _h into ł2, is a direct citation of Theorem 4.9 in <cit.>.
(Discrete compactness property) Any sequence {v̱_h}_h>0 with v̱_h∈_h that is uniformly bounded in H̱(curl,Ω) contains a subsequence that converges strongly in Ḻ^2(Ω).


In the remainder of this subsection, we will prove the error estimates for the discrete forms (<ref>)-(<ref>),
(<ref>) or (<ref>) with λ_h0.
The authors in <cit.> have built a general analysis framework for
the a priori error estimates of mixed form (see Theorem 2.2 and Lemma 2.3 therein). Although we cannot directly apply their theoretical results to
the mixed discretization (<ref>)-(<ref>), we can use its proof idea
to derive the following Lemma 2.2 and Theorem 2.3.
The following uniform convergence provides us with the possibility
 to use the spectral approximation theory in <cit.>.
 There holds the uniform convergence

    T-T_h_ł2→0, h→0.






Since ∪_h>0V̱_h and ∪_h>0  U_h are dense in 0̧ and H^1_0(Ω), respectively, we deduce from (<ref>) for any f̱∈ł2

    T 𝐟- T_h𝐟_A≤(inf_𝐯_h∈𝐕_hT𝐟-𝐯_h_curl+inf_v_h∈  U_h|Sf̱-v_h|_1) 0.

 That is,  T_h converges to T pointwisely.
 Since T, T_h : ł20̧ are linear bounded uniformly with respect to h, ∪_h>0( T- T_h)B is a bounded set in 0̧ where B is the unit ball in ł2. From
↪ł2 compactly and the discrete compactness property of _h in Lemma 2.1, we
know that ∪_h>0( T- T_h)B is a relatively compact set in ł2, which implies  collectively
compact convergence T_h T.  Noting T,  T_h: ł2ł2 are self-adjoint,   due to Proposition 3.7 or Table 3.1 in <cit.> we get  T-T_h_ł2→0, h→0.
 This ends the proof.

































 



























































































































Prior to proving the error estimates for   edge element discretizations, we define some notations as follows.
Let λ be the kth eigenvalue of (<ref>) or (<ref>)-(<ref>) of multiplicity q.
Let λ_j,h (j=k,k+1,⋯,k+q-1) be eigenvalues of T_h that converge to the eigenvalue λ=λ_k=⋯=λ_k+q-1.
Here and hereafter we use M(λ) to denote the space spanned by all eigenfunctions  corresponding to the eigenvalue λ,
and M_h(λ) to denote the direct sum of all eigenfunctions corresponding to the eigenvalues λ_j,h (j=k,k+1,⋯,k+q-1).
For argument convenience,
hereafter we denote λ_j=λ_j+γ/β and
 λ_j,h=λ_j,h+γ/β.
 Now we introduce the following small quantity:

    δ_h(λ)=sup_𝐮∈M(λ),u̱_A=1inf_𝐯∈V̱_h𝐮-𝐯_A.


Thanks to (<ref>) and (<ref>) we have

    λ(   T_h-   T)|_M(λ)_A≤λsup_u̱∈M(λ),u̱_A=1inf_v̱_h∈V̱_h  Tu̱-  v̱_h_curl≤√(γ)δ_h(λ).












The error estimates of edge elements for  the Maxwell eigenvalue problem have been obtained in, e.g., <cit.>. Here we would like to use the quantity δ_h(λ) to characterize the error
for eigenpairs.
From the  spectral approximation, we actually derive the a priori error estimates for the discrete eigenvalue
problem  (<ref>) with λ_h0, (<ref>) or (<ref>)-(<ref>).

Let λ  be the  eigenvalue of (<ref>) or (<ref>)-(<ref>) and let  λ_h be the discrete eigenvalue of
(<ref>) or (<ref>)-(<ref>) converging to
λ.
There exist   h_0 > 0 such that if h ≤ h_0 then
for any eigenfunction u̱_h corresponding to λ_h with u̱_h_A = 1 there exists u̱∈
M(λ) such that

    u̱-u̱_h_A≤δ_h(λ)

and for any u̱∈ M(λ) with  u̱_A=1 there exists u̱_h∈ M_h(λ) such that

    u̱-u̱_h_A≤δ_h(λ),

where the positive constant  is independent of   mesh parameters.

 We take λ=λ_k.
Suppose  u̱_h is an eigenfunction of (<ref>)-(<ref>) corresponding to λ_h satisfying
u̱_h_A=√(λ_h)u̱_h_0,ϵ = 1.
Then according to Theorems 7.1 and 7.3 in <cit.> and Lemma 2.2 there exists u̱∈ M(λ) satisfying

    u̱_h-u̱_0,ϵ≲( T- T_h)|_M(λ)_0,ϵ,
          |λ_j,h-λ|≲( T- T_h)|_M(λ)_0,ϵ for j=k,⋯,k+q-1.

By a simple calculation, we deduce

    |u̱_h-u̱_A- λ( T- T_h)u̱_A|   =   |λ_h  T_hu̱_h-λ Tu̱_A-
    λ( T- T_h)u̱_A|
       ≤    T_h(λ_hu̱_h-λu̱)_A
       ≲   λ_hu̱_h-λu̱_0,ϵ
       ≲   |λ_h-λ|u̱_h_0,ϵ+λu̱_h- u̱_0,ϵ.

Since the equality  (<ref>) implies ( T- T_h)|_M(λ)_0,ϵ≲δ_h(λ),
this together with (<ref>)-(<ref>)  yields (<ref>).
Conversely, suppose u̱ is an eigenfunction of (<ref>)-(<ref>) corresponding to λ
 satisfying u̱_A =√(λ)u̱_0,ϵ =  1. Then according to    Theorems 7.1 in <cit.> and Lemma 2.1 there exists u̱_h∈ M_h(λ) satisfying
 
    u̱_h-u̱_0,ϵ≲( T- T_h)|_M(λ)_0,ϵ.

 Let u̱_h=∑_j=k^k+q-1u̱_j,h where
 u̱_j,h  is the eigenfunction corresponding to λ_j,h such that
{u̱_j,h}_j=k^k+q-1 constitutes an orthogonal basis of M_h(λ) in
(ϵ·,·). Then

    |u̱_h-u̱_A- λ( T- T_h)u̱_A|≤u̱_h- T_h(λu̱)_A
             ≲ T_h(∑_j=k^k+q-1λ_j,hu̱_j,h-λu̱)_A
             ≲∑_j=k^k+q-1 (λ_j,h-λ_h)u̱_j,h+λ_hu̱_h-λu̱_0,ϵ.

Since  ( T- T_h)|_M(λ)_0,ϵ≲δ_h(λ),
this together with (<ref>)-(<ref>) yields (<ref>).


Remark 2.1. Based on the estimate (<ref>), one can naturally obtain the optimal convergence order
𝒪(δ_h^2(λ)) for  |λ_h-λ| using the Rayleigh quotient relation (<ref>)
in the following section.
In addition, note that when λ0 in Theorem 2.2 the estimate (<ref>) implies u̱_h_a
converges to  u̱_a=√(λ)u̱_0,ϵ>0.
Here we introduce u̱_h= u̱_h/u̱_h_a then u̱_h= u̱_h/u̱_h_A and
(<ref>) gives

    u̱-u̱_h_A≤δ_h(λ).

For simplicity of notation, we still use the same C_3 and u̱ in the above estimate as in (<ref>)-(<ref>).

Remark 2.2.



When Ω is a Lipschitz polyhedron and ϵ,μ are properly smooth, it is known that
⊂(̱H^σ(Ω))^3 (σ∈(1/2,1]) (see <cit.>) and
 δ_h(λ)≲ h^σ. In particular, if ⊂{v̱∈H̱^s(Ω):curlv̱∈H̱^s(Ω)} (1≤ s≤ k+1) then
δ_h(λ)≲ h^s (see  Theorem 5.41 in <cit.>).


§ MULTIGRID SCHEMES BASED ON SHIFTED INVERSE ITERATION






 §.§ Multigrid Schemes

In practical computation, the information on the physical zero eigenvalues can be easily captured on a coarse mesh H
 using  the mixed discretization (<ref>)-(<ref>).
In this section we shall present our multigrid methods for solving nonzero Maxwell eigenvalue.
The following schemes are proposed by  <cit.>. Note that we assume in the following schemes the numerical eigenvalue
λ_H approximates the nonzero eigenvalue λ.

Scheme 3.1.   Rayleigh quotient iteration.

Given the maximum number of iterative times l.

Step 1. Solve the eigenvalue problem (<ref>)-(<ref>) on coarse finite element space V̱_H× U_H: find (λ_H,𝐮_H,σ_H)∈ R×V̱_H× U_H,
𝐮_H_a=1 such that

    a(𝐮_H,𝐯)+b(𝐯,σ_H)   =   λ_H(ϵ𝐮_H,𝐯),  ∀𝐯∈V̱_H,
    
    b(𝐮_H,p)   =   0,  ∀ p∈ U_H.

Step 2. 𝐮^h_0⇐𝐮_H, λ^h_0⇐λ_H, i⇐ 1.

Step 3. Solve an equation on V̱_h_i: find (𝐮',σ')∈V̱_h_i such that

    a(𝐮',𝐯)-λ^h_i-1(ϵ𝐮',𝐯)   =   (ϵ𝐮^h_i-1,𝐯),  ∀𝐯∈V̱_h_i.

Set 
𝐮^h_i=𝐮'/𝐮'_a .
    


Step 4. Compute the Rayleigh quotient

    λ^h_i=a(𝐮^h_i,𝐮^h_i)/(ϵ𝐮^h_i,𝐮^h_i).

Step 5. If i=l, then output (λ^h_i,𝐮^h_i), stop; else, i⇐ i+1, and return to step 3.


In Step 3 of the above Scheme, when the shift λ^h_l-1 is close to the exact eigenvalue enough,
the coefficient matrix of linear equation is nearly singular. Hence the following algorithm gives a natural
way of handling this problem.

Scheme 3.2.   Inverse iteration with fixed shift.

Given the maximum number of iterative times l and i0.

Step 1-Step 4. The same as Step 1-Step 4 of Scheme 3.1.


Step 5. If i>i0 then i⇐ i+1 and return to Step 6; else i⇐ i+1 and return to Step 3.

Step 6. Solve an equation on V̱_h_i: find (𝐮',σ')∈V̱_h_i such that

    a(𝐮',𝐯)-λ^h_i0(ϵ𝐮',𝐯)   =   (ϵ𝐮^h_i-1,𝐯),  ∀𝐯∈V̱_h_i.

Set 
𝐮^h_i=𝐮'/𝐮'_a.
    


Step 7. Compute the Rayleigh quotient

    λ^h_i=a(𝐮^h_i,𝐮^h_i)/(ϵ𝐮^h_i,𝐮^h_i).

Step 8. If i=l, then output (λ^h_i,𝐮^h_i), stop; else, i⇐ i+1, and return to step 6.


Remark 3.1. The mixed discretization (<ref>)-(<ref>) was adopted by the literatures
 <cit.>. As is proved in Theorem 2.2, using this discretization  we can compute the Maxwell eigenvalues
 without introducing    spurious eigenvalues. However, it is also a saddle point problem that is difficult
 to solve on a fine mesh (see <cit.>). Therefore, the multigrid schemes can properly overcome the difficulty since
 we only solve (<ref>)-(<ref>)
  on a coarse mesh, as shown in step 1 of Schemes 3.1 and 3.2. Moreover, in order to further improve the
   efficiency of solving the equation in Steps 3 and 6 in Schemes 3.1 and 3.2 the HX preconditioner in
   <cit.> is a good choice (see <cit.>).










 §.§ Error Analysis


In  this subsection, we aim to prove the error estimates for Schemes 3.1 and 3.2.
We shall analyze the constants in the error estimates are independent of mesh parameters and iterative times l.
First of all, we give
two useful lemmas.

For any nonzero 𝐮,𝐯∈0̧,   there hold

    𝐮/𝐮_A-𝐯/𝐯_A_A≤2𝐮-𝐯_A/𝐮_A,    𝐮/𝐮_A-𝐯/𝐯_A_A≤2𝐮-𝐯_A/𝐯_A.


See <cit.>.



Let (λ,u̱) be an eigenpair of
(<ref>) or of (<ref>) with λ0, then for any v ∈0̧\{0}, the
Rayleigh quotient R(v̱)=a(v̱,v̱)/v̱_0,ϵ^2 satisfies

    R(v̱)-λ=v̱-u̱_a^2/v̱_0,ϵ^2-λv̱-u̱_0,ϵ^2/v̱_0,ϵ^2.



See pp.699 of <cit.>.




 The basic relation in Lemma 3.2 cannot be directly applied to our theoretical analysis, so in the following  we shall
further simplify   the estimate (<ref>).
 Let C= (γ/β)^1/2 then according to the definition of A(·,·),

    v̱_0,ϵ≤C^-1v̱_A,   ∀v̱∈0̧.

If u̱∈ M(λ), v̱_h∈V̱_h, v̱_h_A=1 and
v̱_h-u̱_A≤C(4√(λ))^-1, then by Lemma 3.1 we deduce

    v̱_h-u̱/u̱_A_A≤ 2v̱_h-u̱_A≤C(2√(λ))^-1,
          v̱_h-u̱/u̱_A_0,ϵ≤C^-1v̱_h-u̱/u̱_A_A≤ (2√(λ))^-1,

which together with
u̱_A=√(λ)u̱_0,ϵ  yields

    v̱_h_0,ϵ≥u̱_0,ϵ/u̱_A-v̱_h-u̱/u̱_A_0,ϵ≥
    (2√(λ))^-1.

Hence, from Lemma <ref> we get the following estimate





    |R(v̱_h)-λ|≤ C_4v̱_h-u̱_A^2,

where 
C_4=4λ(1+λC^-2).

Define the    operators T: ł20̧ and T_h:ł2
→𝐕_h  as

    A(Tf̱,v̱)=  (ϵf̱,v̱), ∀v̱∈0̧,
          A(T_h𝐟,𝐯_h)=   (ϵ𝐟,𝐯_h), ∀𝐯_h∈𝐕_h.


The following lemma   turns our attention from the spectrum  of T and T_h
to that of T  and T_h.
  T, T and T share the    eigenvalues   greater than γ/β
 and the associated eigenfunctions.
 The same conclusion is valid for   T_h, T_h and T_h.
 Moreover, T|_= T|_ =T|_ and
 T_h|_= T_h|_  =T_h|_.
  
 The assertions  regarding the relations among T, T, T_h and T_h have been described in section 2.
Next we shall prove the relations between T and  T and between T_h and T_h. By the definition of T and
T, the eigenpair (λ,u̱) of T satisfies
A(𝐮,𝐯)=   λ(ϵ𝐮,𝐯) for all 𝐯∈𝐗
and the eigenpair (λ,u̱) of T satisfies
A(𝐮,𝐯)=   λ(ϵ𝐮,𝐯) for all 𝐯∈0̧.
Note that the above two weak forms are equivalent when λ>γ/β (since this implies
the eigenfunction u̱ of the latter satisfies   divergence-free constraint).
Hence T and T share the    eigenvalues  λ>γ/β
 and the associated eigenfunctions. Similarly one can check T_h and T_h share the    eigenvalues
 λ_h>γ/β
 and the associated eigenfunctions.
 Thanks to Helmholtz decomposition 0̧=∇H_0^1(Ω)⊕𝐗 and (<ref>), we also have for any f̱∈

    A(T𝐟,𝐯)=(ϵ𝐟,𝐯), ∀𝐯∈0̧.

This together with (<ref>) yields T|_= T|_.
 Thanks to discrete Helmholtz decomposition 𝐕_h=∇ U_h⊕𝐗_h and (<ref>), we also have for any f̱∈_h+

    A(T_h𝐟,𝐯)=   (ϵ𝐟,𝐯), ∀𝐯∈𝐕_h.

This together with (<ref>) yields T_h|_= T_h|_.



Denote dist(𝐰,W)=inf_𝐯∈ W𝐰-𝐯_A.
 For  better understanding of notations, hereafter
 we write ν_k=λ^-1,ν_j,h=λ_j,h^-1, and
 M_h(ν_k)=M_h(λ_k).

The following lemma  (see <cit.>) is valid since T_h and T_h share the same eigenpairs.
It provides  a crucial   tool for  analyzing the error of multigrid Schemes 3.1 and 3.2.




Let (ν_0,𝐮_0) is an approximate eigenpair of (ν_k,𝐮_k), where ν_0 is not an eigenvalue of T_h and  𝐮_0∈V̱_h with 𝐮_0_a=1. Suppose that

    dist(𝐮_0,M_h(ν_k))≤1/2, 
          |ν_0-ν_k|≤ρ/4 ,   |ν_j,h-ν_j|≤ρ/4(j= k-1,k,k+q,j≠0),

where ρ=min_ν_j≠ν_k|ν_j-ν_k|.
Let 𝐮^s∈V̱_h,𝐮_k^h∈V̱_h satisfy

    (ν_0- T_h)𝐮^s=𝐮_0,   𝐮_k^h=𝐮^s/𝐮^s_a.

Then

    dist(𝐮_k^h, M_h(ν_k))≤4/ρmax_k≤ j≤ k+q-1|ν_0-ν_j,h|dist(𝐮_0,M_h(ν_k)).











Let δ_0 and δ'_0 be two positive constants such that






    δ_0≤min{C(4√(λ_k))^-1,1/2},  δ'_0≤λ_k/2,  δ'_0/(λ_k- δ'_0 )λ_k≤ρ/4,
          δ_0^2<λ_j, ^2δ_0^2   /(λ_j-δ_0^2 )λ_j≤ρ/4, 
    j=k-1,k,k+q,j≠0,
           (3+) δ_0+
    3C^-2δ_0^2 +C^-2(3λ_k+2C^2) δ_0≤ 1/2.



In the coming theoretical analysis, in the step 3 of Scheme 3.1 and step 6 of Scheme  3.2 we introduce a new
auxiliary  variable 𝐮^h_i satisfying

    𝐮^h_i=𝐮'/𝐮'_A.

Then it is clear that 𝐮^h_i=𝐮^h_i/𝐮^h_i_a
and λ^h_i=a( 𝐮^h_i, 𝐮^h_i)/(ϵ 𝐮^h_i, 𝐮^h_i).


C̱o̱ṉḏi̱ṯi̱o̱ṉ ̱3̱.̱1̱.̱ There exists u̱_k   ∈ M(λ_k) 

such that for some   i∈{1,2,⋯,l}

    u̱_k^h_l-1-u̱_̱ḵ_A≤δ_0,  δ_h_l (λ_j ) ≤δ_0 (j = k - 1, k,k+q, j ≠ 0),
          |λ^h_i_k-λ_k| ≤δ'_0,

where λ^h_i_k and (λ^h_l-1_k,u̱_k^h_l-1)
 are   approximate eigenpairs corresponding to the eigenvalue λ_k
obtained by Scheme 3.1 or Scheme 3.2.

We are in a position to prove a critical theorem which establishes the error relation for approximate eigenpairs between
two adjacent iterations.
Our proof
shall   sufficiently make use of the relationship among the operators  T, T_h_l, T, T_h_l, T and
 T_h_l, as
shown in Lemma 3.3, and   the proof method is an extension of that   in <cit.>.


Let (λ_k^h_l,u̱_k^h_l) be an approximate eigenpair obtained by Scheme 3.1
or Scheme 3.2. Suppose Theorem 2.2 holds with  λ=λ_k-1,λ_k,λ_k+q,  and  Condition 3.1   holds with i=l-1 for Scheme 3.1 or with i=i0,l-1 for
Scheme 3.2. Let λ_0=λ_k^h_l-1 for
Scheme 3.1 or λ_0=λ_k^h_i0 for Scheme 3.2.   


Then there  exists 𝐮_k∈ M(λ_k) such that

    𝐮^h_l_k-𝐮_k_A≤/2(|λ_0-λ_k|(|λ^h_l-1_k-
      λ_k|+𝐮^h_l-1_k-𝐮_k_0,ϵ)
      +δ_h_l(λ_k)),

where   is independent of the mesh parameters  and the iterative times l.

 
 






 Step 3 of Scheme 3.1 with i = l  is equivalent to: find (𝐮_h_l',σ')∈ U_h_l× V_h_l such that
 
    A(𝐮',𝐯)-(λ_0+C^2)A( T_h_l𝐮',𝐯)=A( T_h_l𝐮^h_l-1_k,𝐯),   ∀𝐯∈V̱_h_l,

and  𝐮^h_l_k=𝐮'/𝐮'_a,  𝐮^h_l_k=𝐮'/𝐮'_A.














That is

    ((λ_0+C^2)^-1- T_h_l)𝐮'=(λ_0+C^2)^-1 T_h_l𝐮^h_l-1_k,   𝐮̱^h_l_k=𝐮'/𝐮'_A.














Denote

    ν_0=(λ_0+C^2)^-1,   𝐮_0=(λ^h_l-1_k+C^2) T_h_l𝐮^h_l-1_k
    /(λ^h_l-1_k+C^2) T_h_l𝐮^h_l-1_k_A,
    
          𝐮^s=(λ_0+C^2)𝐮'/(λ^h_l-1_k+C^2) T_h_l𝐮^h_l-1_k_A,    ν_h_l=1/λ^h_l _k .

Noting 𝐮_k^h_l-1=𝐮^h_l-1_k/𝐮^h_l-1_k_a, then Step 3 of Scheme 3.1   is equivalent to:

    (ν_0- T_h_l)𝐮^s=𝐮_0,   𝐮^h_l_k=𝐮^s/𝐮^s_A.










Noting 𝐮_k_A≤1+δ_0≤ 3/2, using Lemma 3.3 we derive from (<ref>) and (<ref>)

    (λ^h_l-1_k+C^2) T_h_l𝐮_k-𝐮_k_A=
      (λ_k+C^2) T𝐮_k-(λ^h_l-1_k+C^2) T_h_l𝐮_k_A 
             ≤(λ_k+C^2)( T- T_h_l)𝐮_k_A +(λ_k-λ^h_l-1_k) T_h_l𝐮_k_A 
             ≤δ_h_l(λ_k)𝐮_k_A +
      C^-1|λ^h_l-1_k-λ_k|𝐮_k_0,ϵ
             ≤3/2δ_h_l(λ_k) +
      3/2C^-2|λ^h_l-1_k-λ_k|.







By (<ref>), (<ref>) and (<ref>), we have

    𝐮_0-𝐮_k/𝐮_k_A_A≤ 2
    (λ^h_l-1_k+C^2) T_h_l𝐮^h_l-1_k-𝐮_k_A
             ≤2((λ^h_l-1_k+C^2) T_h_l𝐮_k-𝐮_k_A+(λ^h_l-1_k+C^2) T_h_l(𝐮_k-𝐮^h_l-1_k)_A )
             ≤2(λ^h_l-1_k+C^2) T_h_l𝐮_k-𝐮_k_A+C^-1(3λ_k+2C^2)𝐮_k-𝐮^h_l-1_k_0,ϵ .










We shall verify the conditions of Lemma <ref>.
Recalling (<ref>), (<ref>) and (<ref>), the estimates (<ref>) and (<ref>) lead to

    dist(𝐮_0,M_h_l(λ_k))≤𝐮_0-𝐮_k/𝐮_k_A_A
    +dist(𝐮_k/𝐮_k_A, M_h_l(λ_k) )
             ≤ (3+) δ_h_l(λ_k)+
    3C^-2|λ_k^h_l-1-λ_k| +C^-1(3λ_k+2C^2)𝐮_k-𝐮^h_l-1_k_0,ϵ
             ≤ (3+) δ_0+
    3C^-2δ_0^2 +C^-2(3λ_k+2C^2) δ_0
             ≤ 1/2.

Due to Condition 3.1 we have from (<ref>)

    |ν_k-ν_0|=|λ_0-λ_k|/|(λ_0+C^2)λ_k|≤δ'_0/(λ_k- δ'_0 )λ_k≤ρ/4.

Since by (<ref>), (<ref>) and (<ref>) we get

    λ_j,h_l≥λ_j-|λ_j-λ_j,h_l|
    
    ≥λ_j-δ_h_l^2(λ_j)
    ≥λ_j-δ_0^2>0

and then for j=k-1,k,k+q,j≠0

    |ν_j-ν_j,h_l|=|λ_j-λ_j,h_l|/|λ_j λ_j,h_l|≤^2δ_0^2   /(λ_j-δ_0^2 )λ_j≤ρ/4.

Therefore the   conditions of Lemma <ref> hold, and we have

    dist(𝐮_k^h_l, M_h_l(λ_k))≤4/ρmax_k≤ j≤ k+q-1|ν_j,h_l-ν_0|dist(𝐮_0,M_h_l(λ_k)).


Applying (<ref>), (<ref>) and (<ref>) we have for j=k,k+1,⋯,k+q-1

    |ν_j,h_l-ν_0|=|λ_0-λ_j,h_l|/|(λ_0+C^2) λ_j,h_l|≤|λ_0-λ_k|+|λ_k-λ_j,h_l|/(λ_k- δ'_0 ) (λ_k-δ_0^2)
             ≤|λ_0-λ_k|+δ_h_l^2(λ_k)/(λ_k- δ'_0 ) (λ_k-δ_0^2).

 Substituting (<ref>) and (<ref>) into (<ref>), we have

    dist(𝐮^h_l_k, M_h_l(λ))≤4/ρ(
    |λ_0-λ_k|+δ_h_l^2(λ_k)/|(λ_k- δ'_0 ) (λ_k-δ_0^2)|) ×
             ((3+) δ_h_l(λ_k)+
    3C^-2|λ^h_l-1_k-λ_k| +C^-1(3λ_k+2C^2)𝐮_k-𝐮^h_l-1_k_0,ϵ).

Let 𝐮_j,h be the eigenfunction corresponding to λ_j,h such that
{𝐮_j,h}_j=k^k+q-1 constitutes an orthonormal basis of M_h(λ) in the sense of norm ·_A.
Let 𝐮^* =  ∑_j=k^k+q-1A(𝐮^h_l_k,𝐮_j,h_l) 𝐮_j,h_l then 𝐮^h_l_k-𝐮^*_A=  dist(𝐮^h_l_k, M_h_l(λ_k)).

 From Theorem 2.2, we know there exists {𝐮_j^0}_k^k+q-1⊂ M(λ_k) such that 𝐮_j,h_l-𝐮_j^0 satisfies (<ref>) and it holds
by taking 𝐮_k =  ∑_j=k^k+q-1A(𝐮^h_l_k,𝐮_j,h_l)𝐮_j^0

    𝐮_k-𝐮^*_A   =   ∑_j=k^k+q-1A(𝐮^h_l_k,𝐮_j,h_l)(𝐮_j^0-𝐮_j,h_l)_A
       ≤    (∑_j=k^k+q-1𝐮_j^0-𝐮_j,h_l^2_A)^1/2
       ≤   √(q)δ_h_l(λ_k).






Therefore,  summing up (<ref>) and (<ref>), we know 




there exists a positive constant ≥ that is independent of mesh parameters
and l such that (<ref>)  holds.

















































































Condition 3.2.  For any given ε∈ (0,2), there exist
 t_i∈ (1,3-ε] such that δ_h_i(λ_k)=δ^t_i_h_i-1(λ_k)
 and δ_h_i(λ_k)→ 0 (i→∞).


Condition 3.2 is easily satisfied. For example, for smooth solution, by using
the uniform mesh, let h_0 =√(2)/8, h_1 =√(2)/32, h_2 =√(2)/64 and h_3=√(2)/128, we have h_i = h_i-1^t_i,
i.e., δ_h_i = δ_h_i-1 ^t_i, where t_1≈ 1.80,
t_2 ≈ 1.22, t_3 ≈ 1.18. For non-smooth solution, the condition could be satisfied
when the local refinement is performed near   the singular points.


Let (λ_k^h_l, u̱_k^h_l) be the approximate eigenpairs
obtained by Scheme 3.1. Suppose  Condition
3.2 holds. Then there exist  u̱_k∈ M(λ_k) and H_0>0 such that if H≤ H_0 then

    u̱_k^h_l-u̱_k_A≤  C_0δ_h_l(λ_k),
    
    
          |λ_k^h_l-λ_k|≤ C_4 C_0^2δ_h_l^2(λ_k),    l≥1.



The proof is completed by using induction and Theorem 3.5 with
λ_0=λ_k^h_l-1.
 Noting that δ_H(λ_k) 0 as H 0,  there exists H_0>0 such that if H<H_0 then Theorem 2.2
 holds for λ=λ_k-1,λ_k,λ_k+q  and

    δ_H(λ_k)≤δ_0, ^2δ_H^2(λ_k)≤δ'_0, δ_H(λ_j)≤δ_0,(j=k-1,k,k+q,j0),
          C_4^2C_0^4δ_H^1+ε(λ_k)+C_4C_0^3C^-1δ_H^ε(λ_k)≤ 1.

 When l=1,
(λ_k^h_l-1, u̱_k^h_l-1)=(λ_k,H, u̱_k,H),
from (<ref>) and (<ref>) we know that
there exists u̱_k∈ M(λ_k) such that

    u̱_k,H-u̱_k_A≤ C_3δ_H(λ_k),
    
          |λ_k,H-λ_k|≤
    C_4C_3^2δ_H^2(λ_k).



Then u̱_k^h_0-u̱_k_A≤δ_H(λ_k)≤δ_0, |λ_k^h_0-λ_k|≤
C_4^2δ_H^2(λ_k)≤δ'_0 and
δ_h_1(λ_j)≤δ_0 (j=k-1,k,k+q,j0), i.e.,
Condition 3.1 holds for l=1. Thus, by Theorem 3.5 and C_3≤ C_0 we
get

    u̱_k^h_1-u̱_k_A≤/2{C_4^2C_0^4δ_H^4(λ_k)
    +C_4C_0^3C^-1δ_H^3(λ_k)+δ_h_1(λ_k)}
             ≤/2{C_4^2C_0^4δ_H^4-t_1(λ_k)+C_4C_0^3C^-1δ_H^3-t_1(λ_k)+1}δ_h_1(λ_k)
             ≤/2{C_4^2C_0^4δ_H^1+ε(λ_k)+C_4C_0^3C^-1δ_H^ε(λ_k)+1}δ_h_1(λ_k),































where we have used  the fact  3-t_1≥ε. This yields (<ref>) and (<ref>)
for l=1. 




Suppose that Theorem 3.6 holds for l-1, i.e.,   there exists
u̱_k∈ M(λ_k) such that

    u̱_k^h_l-1-u̱_k_A≤  C_0δ_h_l-1(λ_k),
    
    
          |λ_k^h_l-1-λ_k|≤ C_4 C_0^2δ_h_l-1^2(λ_k),

then u̱_k^h_l-1-u̱_A≤δ_0,
|λ_k^h_l-1-λ_k|≤δ'_0
 and δ_h_l(λ_j)≤δ_0 (j=k-1,k,k+q,j0),
and the conditions of Theorem 3.5 hold. Therefore, for l, by
(<ref>) we deduce

    u̱_k^h_l-u̱_k_A≤/2{C_4^2C_0^4δ_h_l-1^4(λ_k)
    +C_4C_0^3C^-1δ_h_l-1^3(λ_k)
    +δ_h_l(λ_k)}
             ≤/2{C_4^2C_0^4δ_h_l-1^4-t_i(λ_k)
    +C_4C_0^3C^-1δ_h_l-1^3-t_i(λ_k)
    +1}δ_h_l(λ_k)
             ≤/2{C_4^2C_0^4δ_H^1+ε(λ_k)
    +C_4C_0^3C^-1δ_H^ε(λ_k)
    +1}δ_h_l(λ_k)
             ≤ C_0δ_h_l(λ_k),

i.e., (<ref>)  are valid. And from
(<ref>) and (<ref>) we get (<ref>). This ends the proof.


 Condition 3.3. 
There exist β_0∈ (0,1) and
β_i∈ [β_0, 1) (i=1,2,⋯) such that
δ_h_i(λ_k)=β_iδ_h_i-1(λ_k)
 and δ_h_i(λ_k)→ 0 (i→∞).


Remark 3.2.
 Note that if Condition 3.3 is valid, Condition 3.2 holds for H properly small;
however, the inverse is not true. So in Theorem 3.6, (<ref>) and (<ref>) still hold if we replace
Condition 3.2 with  Condition 3.3.



Let (λ_k^h_l, u̱_k^h_l)  be an approximate
eigenpair  obtained by Scheme 3.2. Suppose that  Condition 3.2 holds for i≤ i0
and Condition 3.3 holds for i> i0. Then there exist  u̱_k∈ M(λ_k) and H_0>0 such that if H≤ H_0 then

    u̱_k^h_l-u̱_k_A≤
    C_0δ_h_l(λ_k),
          |λ_k^h_l-λ_k|≤
    C_4C_0^2δ_h_l^2(λ_k),   l> i0.



The proof is completed by using induction. Noting    δ_H(λ_k) 0 as H 0,
there exists H_0>0 such that if H<H_0  then  Theorems 2.2 holds for λ=λ_k-1,λ_k,λ_k+q,
 Theorem 3.6  holds  and

    δ_H(λ_k)≤δ_0, ^2δ_H^2(λ_k)≤δ'_0, δ_H(λ_j)≤δ_0,(j=k-1,k,k+q,j0),
          C_4^2C_0^4δ_H^3(λ_k)1/β_0
    +C_4C_0^3/Cδ_H^2(λ_k)1/β_0≤ 1.


 From Theorem 3.6 we know that  when l=i0,i0+1 there exists u̱_k∈ M(λ_k) such that

    u̱_k^h_l-u̱_k_A≤
    C_0δ_h_l(λ_k),
           |λ_k^h_l-λ_k|≤
    C_4C_0^2δ_h_l^2(λ_k).

Suppose  Theorem 3.7 holds for l-1, i.e., there exists
u̱_k∈ M(λ_k) such that

    u̱_k^h_l-1-u̱_k_A≤
    C_0δ_h_l-1(λ_k),
          |λ_k^h_l-1-λ|≤
    C_4C_0^2δ_h_l-1^2(λ_k).

Then the conditions of Theorem 3.5 hold, therefore, for l, observing that in (<ref>)
u̱_k^h_l-1-u̱_k_0,ϵ can be replaced by
C^-1u̱_k^h_l-1-u̱_k_A, we deduce

    u̱_k^h_l-u̱_k_a≤C_0/2{C_4^2C_0^4δ_h_i0^2(λ_k)(δ_h_l-1^2(λ_k)
    +C_0/Cδ_h_l-1(λ_k))
    +δ_h_l(λ_k)}
             ≤C_0/2{C_4^2C_0^4δ_h_i0^2(λ_k)δ_h_l-1(λ_k)1/β_0
    +C_4C_0^3/Cδ_h_i0^2(λ_k)1/β_0
    +1}δ_h_l(λ_k),

noting that
δ_h_l-1(λ_k)≤δ_h_i0(λ_k)≤δ_H(λ_k), we get
(<ref>)
immediately.
(<ref>) can be obtained from (<ref>) and (<ref>).
The proof is completed.


Remark 3.3.
The error estimates (<ref>) and (<ref>) for u̱^h_l_k can lead to the error estimates for u̱^h_l_k. In fact,
under the conditions of Theorem 3.6 or Theorem 3.7, we have

    u̱_k_A≥u̱^h_l_k_A-δ_h_l(λ_k)≥1-δ_0≥1/2,

then  u̱_k_a=√(λ_k)/√(λ_k)u̱_k_A≥√(λ_k)/2√(λ_k).
We further assume δ_0≤√(λ_k)/4√(λ_k) then


    u̱^h_l_k_a≥u̱_k_a-δ_h_l(λ_k)≥√(λ_k)/2√(λ_k)-δ_0
    ≥√(λ_k)/4√(λ_k).

Therefore we derive from (<ref>) or (<ref>)

    u̱^h_l_k -u̱_k/u̱_k_a_A≤u̱^h_l_k-u̱_k_Au̱_k_a+u̱^h_l_k-u̱_k_au̱_k_A/u̱^h_l_k_a𝐮_k_a
             ≤(√(λ_k)+√(λ_k))u̱^h_l_k-u̱_k_A/u̱^h_l_k_a√(λ_k)≤4(√(λ_kλ_k)+λ_k)δ_h_l(λ_k)/λ_k,

i.e., u̱^h_l_k has the same convergence order as u̱^h_l_k in the sense of ·_A.





















































































§ NUMERICAL EXPERIMENT


In this section, we will report  several numerical experiments for solving the Maxwell eigenvalue
problem by multigrid Scheme 3.2 using the lowest order edge element
to validate our theoretical results.  We use MATLAB 2012a to compile our program codes  and adopt the data structure of finite elements in  the package
of iFEM <cit.> to generate and refine the meshes.
 





























































We use the sparse solver eigs(A,B,k,'sm')
to solve (<ref>) for   k lowest eigenvalues.
 In our tables 4.1-4.6   we use the
notation λ_k^h_i to denote the numerical eigenvalue approximating    λ_k
   obtained by multigrid methods at ith iteration   on the mesh π_h_i (with number of
   degree of freedom N^h_i), and R to denote the convergence rate with respect to Dof^-1/3 where Dof is  the number of degrees of freedom.
 For comparative purpose, λ_k,h_j denotes the numerical eigenvalue approximating λ_k
   computed by the direct solver eigs on the mesh π_h_j.



Example 4.1. Consider the Maxwell eigenvalue problem with μ=ϵ=I on
the unit cube Ω=(-1/2,1/2)^3. 
We use Scheme 3.2 with i0=0 to compute the  eigenvalues λ_1=2π^2,
λ_4=3π^2,λ_6=5π^2 (of multiplicity 3, 2 and 6 respectively).
 The numerical results are shown in Table  4.1, which indicates that numerical eigenvalues obtained by multigrid methods achieve the optimal convergence rate R≈2.



















































































Example 4.2. Consider the  eigenvalue problem with ϵ=I and μ=I or

    μ=([       2    1-2j      -j; 1 + 2 j       4       j;       j      -j       5 ])

  on the thick L-shaped Ω=((-1,1)^2\ (-1,0]^2)× (0,1).
  When μ=ϵ=I λ_1≈9.6397, λ_2≈11.3452
 and λ_3≈13.4036 (see <cit.>).  We use Scheme 3.2
  with i0=1 to compute the lowest three eigenvalues for both cases. The numerical results are shown
in Tables 4.2-4.3. From Table 4.2 we see that the   eigenvalue errors  obtained by Scheme 3.2
after 2nd iteration are respectively
0.019, 0.012 and 7.0e-04, which indicates the accuracy of  the  lowest two eigenvalues is affected
by the singularity of the associated eigenfunctions  in the directions
perpendicular to the reentrant edge and the convergence rate R is usually less than 2. Alternately, we adopt the meshes  locally refined towards
the reentrant edge (see Figure 4.1) to perform the
iterative procedure. And the associated numerical results are listed in Table 4.4,  which implies
the errors of λ^h_2_1 and λ^h_2_2 are significantly
 decreased to 0.0033 and 0.0066 respectively  with less degrees of freedom.



Example 4.3. Consider the Maxwell eigenvalue problem with Ω=(-1/2,1/2)× (0,0.1)× (-1/2,1/2)
 where μ=I and if x_3>0  then ϵ=2I otherwise ϵ=I.
This is a practical problem in engineering computed in <cit.>.
We use Scheme 3.2
  with i0=1 to compute the three lowest  eigenvalues for both cases. The numerical results are shown
in Table 4.5.  The relatively accurate eigenvalues reported in <cit.> are respectively  3.538^2(12.5174), 5.445^2(29.6480)
and  5.935^2(35.2242). Using them as the reference values, the relative errors of numerical eigenvalues after 3rd iteration
are respectively 6.0877e-05, 3.7524e-05 and  0.0105. Obviously we get the good approximations of
the eigenvalues λ_1 and λ_2. Regarding the computation for λ_3, we refer to
 Table II in <cit.> whose relative error for computing λ_3 is 0.0107 using a higher order
 method. This is a computational result very close to ours.
Hence we think our method is also efficient for solving the problem.



Example 4.4. Consider the Maxwell eigenvalue problem with μ=ϵ=I and Ω=(-1,1)^3\ [-1/2,1/2]^3. In this example, we capture a physical zero eigenvalue on a
coarse mesh with number of degrees of freedom 6230, i.e., λ_1,H=1.9510e-12.
We use Scheme 3.2
  with i0=0 to compute the   eigenvalues λ_2(of multiplicity 3), λ_5(of multiplicity 2) and λ_7(of multiplicity 3). The numerical results are shown
in Table 4.6. Note that the coarse mesh
seems slightly “fine”. This is because we would like to capture all  information of the lowest
 eight eigenvalues (some of them would not be captured on a very coarse mesh). This is an example of handling the problem
in a cavity with two disconnected boundaries. For more numerical examples of the cavity with disconnected
boundaries,
we refer the readers to the work of <cit.>.



Acknowledgment The author wishes to thank Prof. Yidu Yang for many
valuable comments on this paper.



s30

arbenzP. Arbenz, R. Geus, and S. Adam, Solving Maxwell eigenvalue problems for accelerating
cavities, Phys. Rev. Accelerators and Beams, 4 (2001), 022001.

arbenz1P. Arbenz and R. Geus, Multilevel preconditioned iterative eigensolvers for Maxwell eigenvalue
problems, Appl. Numer. Math., 54 (2005), pp. 107-121.













ainsworthM. Ainsworth, J. Coyle, 
 

Computation of Maxwell eigenvalues on curvilinear domains using hp-version Nédélec elements,
Numer. Math. Adv. Appl.,  2003, pp. 219-231.


















boffi1D. Boffi, P. Fernandes, L. Gastaldi, and I. Perugia, Computational models of electromagnetic
resonantors: Analysis of edge element approximation, SIAM J. Numer. Anal.,
36 (1999), pp. 1264-1290.
boffi2D. Boffi, Fortin operator and discrete compactness for edge elements, Numer. Math., 87
(2000), pp. 229-246.
boffi3D. Boffi, Finite element approximation of eigenvalue problems, Acta. Numerica, 19
(2010), pp. 1-120.
buffa1 A. Buffa, P. Ciarlet, and E. Jamelot, Solving electromagnetic eigenvalue problems in
polyhedral domains with nodal finite elements, Numer. Math., 113 (2009), pp. 497-518.


buffa2 A. Buffa, P. Houston, and I. Perugia, Discontinuous Galerkin computation of the Maxwell
eigenvalues on simplicial meshes, J. Comput. Appl. Math., 204 (2007), pp. 317-333.
brezziF. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, vol. 15, Springer, New York, NY, USA, 1991.
brenner
S. C. Brenner, F. Li , L. Sung,
Nonconforming Maxwell Eigensolvers, J. Sci. Comput., 40(2009),   pp. 51-85.

babuskaI. Babuska, J. Osborn, Eigenvalue Problems, in: P .G. Ciarlet, J. L. Lions, (Ed.), Finite Element Methods (Part 1),
Handbook of Numerical Analysis, vol.2, Elsevier Science Publishers,
North-Holand, 1991,  pp. 640-787.



caorsiS. Caorsi, P. Fernandes, M. Raffetto, On the Convergence of Galerkin Finite Element Approximations of Electromagnetic Eigenproblems, SIAM J. Numer. Anal.,  38(2000), pp. 580-607.
chenL. Chen, iFEM: An Integrated Finite Element Methods Package in MATLAB, Technical
report, University of California at Irvine, 2009.
ciarlet1 P. Ciarlet, Jr., and G. Hechme, Computing electromagnetic eigenmodes with continuous
Galerkin approximations, Comput. Methods Appl. Mech. Engrg., 198 (2008),
pp. 358-365.
chen2J. Chen, Y. Xu, and J. Zou, An adaptive inverse iteration for Maxwell eigenvalue problem
based on edge elements, J. Comput. Phys., 229 (2010), pp. 2649-2658.
chatterjeeA. Chatterjee, J. M. Jin, and J. L Volakis, Computation of cavity resonances using edge based
finite elements, IEEE Trans. Microwave Theory Tech., 40 (1992), pp. 2106-2108.
chatelin F. Chatelin, Spectral Approximations of Linear Operators, Academic Press, New York, 1983.
dauge2M. Dauge, Benchmark computations for Maxwell equations for the approximation of
highly singular solutions (2004). See Monique Dauges personal web page at the location
http://perso.univ-rennes1.fr/monique.dauge/core/index.html
girault V. Girault and P. A. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory
and Algorithms, Springer-Verlag, Berlin, 1986.
hiptmair1 R. Hiptmair, Finite elements in computational electromagnetism, Acta Numer., 11 (2002),
pp. 237-339.
hiptmair2 R. Hiptmair and J. Xu, Nodal auxiliary space preconditioning in H(curl) and H(div) spaces,
SIAM J. Numer. Anal., 45 (2007), pp. 2483-2509.
huX. Hu and X. Cheng, Acceleration of a two-grid method for eigenvalue problems, Math. Comp.,
80 (2011), pp. 1287-1301.








milanM. M. Ilić, B. M. Notaroš,  Higher order hierarchical curved hexahedral vector
finite elements for electromagnetic modeling, IEEE Trans. Microwave Theory Tech., 51(2003), pp. 1026-1033.


jiangW. Jiang, N. Liu, Y. Yue,   Q. Liu, Mixed finite element method for resonant cavity problem with complex geometric
topology and anisotropic media, IEEE Trans.   Magnetics, 52(2016), 7400108.
kikuchiF. Kikuchi, Weak formulations for finite element analysis of an electromagnetic eigenvalue
problem, Scientific Papers of the College of Arts and Sciences, University of Tokyo, 38
(1988), pp. 43-67.
kikuchi1F. Kikuchi, On a discrete compactness property for the Nedelec finite elements, J.fac.sci.univ.tokyo Sect. 1A Math, 36(1989), pp. 479-490.
kirschA. Kirsch, P. Monk,    A finite element method for approximating electromagnetic
scattering from a conducting object, Numer. Math.,  92(2002), pp. 501-534.






monkP. Monk, Finite Element Methods for Maxwells Equations, Oxford University Press, Oxford,
UK, 2003.
monk2P. Monk, L. Demkowicz, Discrete compactness and the approximation of Maxwell's equations in ℝ^3, Math. Comp., 70(2000), pp. 507-523.












nedelec  J. C. Nédélec, Mixed finite elements in R^3. Numer. Math., 35(1980), pp. 315-341.

reddy C. J.  Reddy, M. D. Deshpande, C. R. Cockrell, F. B. Beck, Finite Element Method for Eigenvalue Problems in Electromagnetics, Nasa Sti/recon Technical Report N, 95(1995).

russoA. D. Russo, A. Alonso, Finite element approximation of Maxwell eigenproblems on curved Lipschitz polyhedral domains,
Appl. Numer. Math.,  59(2009), pp. 1796-1822
xu2 J. Xu, A. Zhou, A two-grid discretization scheme for eigenvalue problems,
Math. Comp. 70(1999), pp. 17-25.
yang2 Y. Yang, H. Bi, J. Han, Y. Yu, The Shifted-Inverse Iteration Based on the Multigrid Discretizations
for Eigenvalue Problems, SIAM J. Sci. Comput., (37)(2015), pp. A2583-A2606.
























































































































yang3Y. Yang, H. Bi, A two-grid discretization scheme based on shifted-inverse power method,
SIAM J. Numer. Anal.  49 (2011), pp. 1602-1624.

yang4Y. Yang, Y. Zhang, H Bi, Multigrid Discretization and Iterative Algorithm for Mixed Variational Formulation of the Eigenvalue Problem of Electric Field, Abstr. Appl. Anal., 2012(2012), pp. 1-25, doi:10.1155/2012/190768.

zhou J. Zhou, X. Hu, L. Zhong, S. Shu, L. Chen. Two-Grid Methods for Maxwell Eigenvalue Problem,
SIAM J.   Numer. Anal., 52(4)(2014), pp. 2027-2047.



