Three Numerical Eigensolvers for 3-D Cavity Resonators Filled With Anisotropic and Nonconductive Media
11 Three Numerical Eigensolvers for 3-D CavityResonators Filled With Anisotropic andNonconductive Media
Wei Jiang and Jie Liu
Abstract —This paper mainly investigates the classic resonantcavity problem with anisotropic and nonconductive media, whichis a linear vector Maxwell’s eigenvalue problem. The finiteelement method based on edge element of the lowest-order andstandard linear element is used to solve this type of 3-D closedcavity problem. In order to eliminate spurious zero modes in thenumerical simulation, the divergence-free condition supportedby Gauss’ law is enforced in a weak sense. After the finiteelement discretization, the generalized eigenvalue problem witha linear constraint condition needs to be solved. Penalty method,augmented method and projection method are applied to solvethis difficult problem in numerical linear algebra. The advantagesand disadvantages of these three computational methods are alsogiven in this paper. Furthermore, we prove that the augmentedmethod is free of spurious modes as long as the anisotropicmaterial is not magnetic lossy. The projection method basedon singular value decomposition technique can be used to solvethe resonant cavity problem. Moreover, the projection methodcannot introduce any spurious modes. At last, several numericalexperiments are carried out to verify our theoretical results.
Index Terms —Augmented method, penalty method, projectionmethod, resonant cavity, spurious mode.
I. I
NTRODUCTION M ICROWAVE resonant cavity is an important passivedevice in microwave engineering. It has many applica-tions in many projects, such as particle accelerator, microwaveoven and microwave filter. The microwave resonant cavityproblem usually needs to solve the eigenmodes of source-freeMaxwell’s equations. If the cavity is filled with inhomoge-neous media and/or has a complex geometry, then finding itsresonant modes by an analytical method is impossible. In orderto get the resonant mode, numerical methods must be appliedto solve the problem. The main numerical methods includefinite element method, finite difference method and boundaryelement method.Solving 3-D closed cavity problem will introduce manyspurious modes if the numerical method cannot preserve thephysical property of electromagnetic field. These spurious
This work was supported by the National Natural Science Foundation ofChina under Grant 61901131, the Natural Science Foundation of GuizhouMinzu University under Grant GZMU[2019]YB07 and the China PostdoctoralScience Foundation under Grant 2019M662244 (Corresponding author: WeiJiang).
W. Jiang is with the School of Mechatronics Engineering, Guizhou MinzuUniversity, Guiyang 550025, China (e-mail: [email protected]).J. Liu is with the Postdoctoral Mobile Station of Information and Com-munication Engineering, School of Informatics, and Institute of Electromag-netics and Acoustics, Xiamen University, Xiamen 361005, China (e-mail:[email protected]). modes are usually divided into two kinds. One is spuriousnonzero mode and the other one is spurious zero mode.In general, the spurious zero mode in numerical results iscaused by the neglection of divergence-free condition, andthe introduction of the spurious nonzero mode in numericalresults is the result of improper discretization method. It is ourgoal to design a numerical method that can eliminate thesetwo kinds of spurious modes together. We now know edgeelement method can remove all spurious nonzero modes insolving electromagnetic resonant cavity problems. There aremany references on this subject, such as Lee and Mittra [1],Wang and Ida [2], Pichon and Razek [3], and so on. However,the edge element method cannot eliminate spurious zeromodes because it cannot guarantee the solenoidal propertiesof electric and magnetic flux densities in the whole cavity.The main difficulty of solving resonant cavity problems ishow to enforce the divergence-free condition given by Gauss’slaw in electromagnetics. For the first time, F. Kikuchi [4]introduces a Lagrange multiplier to enforce the divergence-free condition in a weak sense, and propose a mixed finiteelement method (MFEM) to solve 3-D empty cavity problem.As a consequence, Kikuchi’s method is free of all the spuriousmodes, including spurious zero modes. Based on Kikuchi’sidea, Jiang et al. [5], [6] make use of MFEM to solve 3-D reso-nant cavity problems with anisotropic media, and successfullyremove spurious zero and nonzero modes together. However,the MFEM supported in [5], [6] cannot deal with 3-D closedcavity problems filled with the anisotropic media, which isboth electric and magnetic lossy. On the basis of the reference[5], Liu et al. [7] give a two-grid vector discrete scheme for3-D cavity problems with lossless media. The scheme givenin [7] is free of all spurious modes and is very efficient if wejust need to know the first few physical modes. Using edgeelement and linear element, Jiang et al. [8] successfully solve3-D closed cavity problem filled with fully conductive mediain the whole cavity. The numerical method given in [8] canalso remove the spurious zero and nonzero modes together.This paper continues to study the microwave cavity prob-lem filled with anisotropic and nonconductive media. Aftereliminating the electric field in source-free Maxwell’s equa-tions, one can get a vector Maxwell’s eigenvalue problemfor magnetic field. This problem can be transformed intoa corresponding variational formulation by using Green’sformulaes. The edge basis functions of the lowest-order andstandard nodal basis functions of linear element are used todiscretize the variational formulation. Finally, we need to solve a r X i v : . [ m a t h . NA ] J u l the generalized eigenvalue problem with a linear constraintcondition, which is a very difficult problem in numerical linearalgebra. In order to overcome this difficult problem, penaltymethod, augmented method and projection method reduce it tothe generalized eigenvalue problem without any constraint. Inaddition, the advantages and disadvantages among these threecomputational methods are also given in the paper.The outline of the paper is as follows. The governingequations and finite element discretization of 3-D resonantcavity problem are given in Section II. In Section III, weprovide the penalty method, augmented method and projectionmethod to solve the constrained generalized eigenvalue prob-lem and discuss the advantages and disadvantages among thesethree numerical computational methods. In Section IV, threenumerical experiments are carried out to verify our theoreticalresults.II. F INITE E LEMENT D ISCRETIZATION OF
ESONANT C AVITY P ROBLEM
A. Governing Equations for 3-D Resonant Cavity Problem
Suppose that Ω is a bounded domain in R , ∂ Ω is theboundary of Ω and ˆ n is the outward normal unit vector on ∂ Ω . Let (cid:15) and µ be the permeability and permittivity invacuum, respectively. The relative permeability and permittiv-ity tensor of an anisotropic medium are denoted by µ r and (cid:15) r ,respectively. The angular frequency of electromagnetic waveis denoted by ω .The governing equations of 3-D closed cavity problem arethe source-free Maxwell’s equations of the first-order. Aftereliminating the electric field E in the source-free Maxwell’sequations, one can get a second-order partial differentialequations (PDEs) for the magnetic field H : ∇ × (cid:16) (cid:15) − r ∇ × H (cid:17) = Λ µ r H in Ω (1a) ∇ · (cid:0) µ r H (cid:1) = 0 in Ω (1b) ˆ n × ( (cid:15) − r ∇ × H ) = on ∂ Ω (1c) ˆ n · ( µ r H ) = 0 on ∂ Ω , (1d)where Λ = ω (cid:15) µ is the square of the wavenumber invacuum. We would like to seek (Λ , H ) with H (cid:54) = such thatPDEs (1) holds. In electromagnetics, (Λ , H ) with H (cid:54) = iscalled a physical resonant eigenmode in the resonant cavity.In mathematics, (Λ , H ) with H (cid:54) = is called an eigenpairof PDEs (1), and Λ and H are called the eigenvalue andeigenfunction of PDEs (1), respectively. In addition, we knowPDEs (1) only has a discrete point spectrum. Note that theremay be several zero eigenmodes in PDEs (1), for details,please see [9].If the anisotropic material is lossless, then µ r and (cid:15) r areboth Hermitian [10], that is µ † r = µ r and (cid:15) † r = (cid:15) r , where thesuperscript † denotes the conjugate transposition. Moreover,when the anisotropic material is lossless, assuming that (cid:15) r and µ r are Hermitian positive definite since the anisotropicand lossless material in nature usually has this property. Forthe sake of simplicity, a Hermitian positive definite matrix M is denoted by M † = M > . In terms of the lossy characteristics of the anisotropic and nonconductive material,it is usually divided into the following four categories:1) Case 1: (cid:15) † r = (cid:15) r > and µ † r = µ r > . The medium islossless.2) Case 2: (cid:15) † r (cid:54) = (cid:15) r and µ † r = µ r > . The medium iselectric lossy, but is not magnetic lossy.3) Case 3: (cid:15) † r = (cid:15) r > and µ † r (cid:54) = µ r . The medium ismagnetic lossy, but is not electric lossy.4) Case 4: (cid:15) † r (cid:54) = (cid:15) r and µ † r (cid:54) = µ r . The medium is bothelectric and magnetic lossy.Under Case 1, 2 and 3, the MFEM can deal with thesetypes of the resonant cavity problems very well, and it isfree of all the spurious modes in numerical results [5], [6].However, the MFEM is not suitable for 3-D closed cavityproblem under Case 4, because it is difficult to propose anappropriate mixed variational formulation. For the 3-D closedcavity problem under Case 4, the projection method introducedin this paper can deal with this problem very well, and theprojection method can remove all the spurious modes. B. Finite Element Discretization
It is well-known that the finite element method is a varia-tional method, and only operates on the weak form of PDE,instead of the strong form of PDE. Hence, the correspondingweak form associated with PDEs (1) is given in the subsection.To get this weak form, we introduce the following Hilbertspaces over C : L (Ω) = (cid:8) f : (cid:90) Ω | f ( x, y, z ) | d x d y d z < + ∞ (cid:9) H (Ω) = (cid:8) f ∈ L (Ω) : ∇ f ∈ ( L (Ω)) (cid:9) H ( curl , Ω) = (cid:8) F ∈ ( L (Ω)) : ∇ × F ∈ ( L (Ω)) (cid:9) . Define the continuous sesquilinear forms: A : H ( curl , Ω) × H ( curl , Ω) → C ( H , F ) → (cid:90) Ω (cid:15) − r ∇ × H · ∇ × F ∗ d Ω M : ( L (Ω)) × ( L (Ω)) → C ( H , F ) → (cid:90) Ω µ r H · F ∗ d Ω C : H ( curl , Ω) × H (Ω) → C ( H , q ) → (cid:90) Ω µ r H · ∇ q ∗ d Ω where the symbol ∗ stands for the complex conjugation of agiven complex-valued function.Using the Green’s formulas, the weak form of PDEs (1)reads as: Seek Λ ∈ C , H ∈ H ( curl , Ω) , H (cid:54) = such that A ( H , F ) = Λ M ( H , F ) , ∀ F ∈ H ( curl , Ω) (2a) C ( H , q ) = 0 , ∀ q ∈ H (Ω) (2b)Under Case 1, the eigenvalues Λ are made up of countablenonnegative real numbers. Under Case 2, 3 and 4, the eigen-values Λ are made up of countable complex numbers. Thephysical interpretation is such a physical fact that there is K Fig. 1. Local nodal numbering for the element K and the local referencedirection for the edge are chosen by means of local nodal numbering. no electromagnetic energy loss in the resonant cavity if thematerial is lossless and there exists electromagnetic energyloss in the resonant cavity provided that the material has adielectric loss.We now consider the conforming finite element discretiza-tion of (2). Let T h be a regular tetrahedral mesh of the cavity Ω . Here h is the length of the longest edge in tetrahedral mesh T h . As usual, the edge element space W h of the lowest-orderunder the mesh T h is used to approximate the Hilbert space H ( curl , Ω) and standard linear element space S h under themesh T h is used to approximate the Hilbert space H (Ω) .From [11], we know S h (cid:40) H (Ω) and W h (cid:40) H ( curl , Ω) .The linear element space S h can be written as: S h = (cid:8) φ : φ | K ∈ span { L K , L K , L K , L K } (cid:9) where L Ki ( i = 1 , , , are four local nodal basis functionson the tetrahedral element K and of the form a i + b i x + c i y + d i z , where a i , b i , c i , d i are four constants. These fourlocal basis functions are defined on the four vertices of thetetrahedral element K .The edge element space W h of the lowest-order can bewritten as: W h = (cid:8) F : F | K ∈ span { N K , N K , · · · , N K } (cid:9) where N Ki ( i = 1 , , · · · , are six local edge basis functionson the tetrahedral element K and N Ki is of the form (cid:126)α i + (cid:126)β i × (cid:126)r , where (cid:126)α i and (cid:126)β i are two constant vectors and (cid:126)r isthe position vector. The concrete expressions of N Ki ( i =1 , , · · · , are as follows: N K = L K ∇ L K − L K ∇ L K , N K = L K ∇ L K − L K ∇ L K N K = L K ∇ L K − L K ∇ L K , N K = L K ∇ L K − L K ∇ L K N K = L K ∇ L K − L K ∇ L K , N K = L K ∇ L K − L K ∇ L K These six local vector basis functions are defined on the sixedges of the tetrahedral element K . In Fig. 1, we give a localnodal numbering in the tetrahedral element K , and specify thelocal reference direction for each edge in K .Restricting (2) on W h × S h , we get the discrete variationalformulation associated with (2): Seek Λ h ∈ C , H h ∈ W h , H h (cid:54) = such that A ( H h , F ) = Λ h M ( H h , F ) , ∀ F ∈ W h (3a) C ( H h , q ) = 0 , ∀ q ∈ S h (3b) Here Λ h and H h are an approximation of the exact eigenvalue Λ and exact eigenfunction H in (2), respectively.Suppose that S h = span (cid:8) L i (cid:9) mi =1 , where L i is the i -th globalnodal basis function associated with the node i and the integer m is the number of the total nodes in T h . Assuming that W h = span (cid:8) N i (cid:9) ni =1 , where N i is the i -th global edge basis functionassociated with the edge i and the integer n is the number ofthe total edges in T h . Since H h ∈ W h , then it can be writtenas H h = n (cid:88) i =1 ξ i N i . (4)Finally, the discrete variational formulation (3) can bereduced to the following generalized eigenvalue problem witha linear constraint: (cid:26) A ξ = Λ h M ξ (5a) C ξ = (5b)where A = ( a ik ) ∈ C n × n , M = ( m ik ) ∈ C n × n , C = ( c ik ) ∈ C m × n , ξ = [ ξ , ξ , · · · , ξ n ] T ∈ C n ,a ik = A ( N k , N i ) , m ik = M ( N k , N i ) , c ik = C ( N k , L i ) . Solving (5a) directly and ignoring (5b) will introduce a lot ofspurious zero modes. In order to eliminate these spurious zeromodes, we must enforce the constraint (5b) in solving (5a).How to enforce the constraint (5b) in solving (5a) is a keyproblem. In next section, we shall deal with this troublesomeproblem. Once the eigenpair (Λ h , ξ ) is found from (5), thenthe numerical eigenvalue in (3) is given by Λ h and thecorresponding numerical eigenfunction H h in (3) is given by(4).In the community of numerical linear algebra, the problem(5) is a constrained generalized eigenvalue problem. Obvi-ously, its numerical computation is much more difficult thanthat of the generalized eigenvalue problem without any con-straint. In Section III, the problem (5) is reduced to three typesof generalized eigenvalue problem without any constraint. Itis important to point out if there is no relation among thesethree matrices A , M and C , then the constrained generalizedeigenvalue problem (5) may have no solution.III. T HREE N UMERICAL S OLVERS OF THE C ONSTRAINED G ENERALIZED E IGENVALUE P ROBLEM (5)In this section, we first try to give a relation among thematrices A , M and C , and then support three numerical com-putational methods of solving (5). They are penalty method,augmented method and projection method, respectively. A. The Relation Among the Matrices A , M and C Let { A , A , · · · , A m − , A m } and { e , e , · · · , e n − , e n } be the sets consisting of all nodes and edges in T h respectively,where , , · · · , m and , , · · · , n are the global labels of allnodes and edges in T h , respectively. Note that the global vectorbasis function N i has a local direction associated with theedge e i . If two vertices of the edge e i are A i and A i , then we state that the direction of N i is from the node A i to thenode A i , where i < i .In accordance with deRham-complex [12], ∇ S h (cid:40) W h holds, where ∇ S h = {∇ f : ∀ f ∈ S h } . This implies that ∇ L i = n (cid:88) k =1 y ik N k , i = 1 , , · · · , m. (6)The above formula (6) is also introduced in [13] and [14]. Itis easy to know dim( ∇ S h ) = m − since (cid:80) mi =1 ∇ L i = .Set Y = y y · · · y n y y · · · y n ... ... ... ... y m y m · · · y mn = y y ... y m = [ d , d , · · · , d n ] , where y k and d k are k -th row and column vector in the matrix Y , respectively.In fact, the each entry y ki in Y is − , or , and Y isquite sparse. Let us consider the formula (6) under the caseof an arbitrary tetrahedral element K in T h (see Fig. 1). It iseasy to verify that the following formulas are valid: ∇ L K = − N K − N K − N K , ∇ L K = N K − N K − N K ∇ L K = N K + N K − N K , ∇ L K = N K + N K + N K Here we need to make use of the relation (cid:80) i =1 L Ki = 1 .Consider the k -th row y k in the matrix Y , which is relatedto the node A k . If A k is not a vertex of the i -th edge e i ,then y ki = 0 . Let us recall the basic concept of degree of avertex in graph theory. The degree of a vertex is defined bythe number of edges connecting it. The number of the nonzeroentries of y k is equal to the degree of A k . Assuming that thedegree of A k is v and { e i , e i , · · · , e i v } have the commonvertex A k . If the direction of e i s points to A k , then y k,i s = 1 ,otherwise y k,i s = − . Consider the k -th column vector d k inthe matrix Y , which is related to the edge e k . If A i is nota vertex of the edge e k , then y ik = 0 . Since each edge in atetrahedral mesh only has two vertices, the column vector d k only has two nonzero entries. Assuming that the initial andterminal points of the edge e k are A i and A i respectively,then y i k = − and y i k = 1 , where i < i . Obviously,the sparse matrix Y can be easily obtained by the mesh datain T h . The sparse matrix Y is usually called the directionalconnectivity matrix associated with a tetrahedral mesh T h .The sum of all entries in the column vector d k is equal tozero. This implies that the homogenous linear equation Y † ζ = (7)has a special nonzero solution β = [1 , , · · · , T ∈ C m . Aneasy computation shows that rank ( Y ) = rank ( Y † ) = m − . Therefore all solutions of (7) form a linear space of onedimension, that is N ull ( Y † ) = span { β } , (8)where Null stands for taking the nullspace of a matrix.The matrix Y builds a relation among A , M and C . It canbe proved that the following matrix identities are valid: Y A = O , C = Y M , (9) where O is a null m × n matrix. In fact, for ∀ ≤ i ≤ m and ∀ ≤ l ≤ n , we have ( Y A ) il = n (cid:88) k =1 y ik a kl = n (cid:88) k =1 y ik A ( N l , N k )= A ( N l , n (cid:88) k =1 y ik N k ) = A ( N l , ∇ L i ) = 0 , ( Y M ) il = n (cid:88) k =1 y ik m kl = n (cid:88) k =1 y ik M ( N l , N k )= M ( N l , n (cid:88) k =1 y ik N k ) = M ( N l , ∇ L i )= C ( N l , L i ) = c il , where ( S ) il is the entry at i -th row and l -th column of thematrix S . It is important to emphasize that the matrix identities(9) are always valid for the medium under Case 1, 2, 3 and4. Hence, we can obtain the matrix C by using the sparsematrices Y and M , instead of the calculation by using thecontinuous sesquilinear forms C directly.If the eigenvalue Λ h is nonzero in (5a), then (5b) can bededuced from (5a). As a matter of fact, M ξ = Λ − h A ξ canbe derived from (5a), and then we get C ξ = Y M ξ =Λ − h Y A ξ = by (9), which is just (5b). It is worthwhileto point out the number of zero eigenvalues in (5a) is equal to dim( ∇ S h ) , and these zero eigenvalues are all spurious. Thedimension of the linear space ∇ S h is equal to m − , whichshows that the larger the number of mesh nodes m is, themore the number of spurious zero modes in (5a) is. Based onthis reason, we need to remove these spurious zero modes byintroducing an appropriate numerical method. B. Penalty Method
Consider the following generalized eigenvalue problem: ( A + α C † C ) ξ = Λ (cid:48) h M ξ, (cid:107) ξ (cid:107) = 1 (10)where α is a penalty constant, which is required the user toset and (cid:107) ϕ (cid:107) is the Euclidean norm of a given vector ϕ . Theproblem (10) is a generalized eigenvalue problem without anyconstraint, and it can be solved by the numerical softwarepackage ARPACK [15].The parameter α is usually set to a large positive realnumber. The reason is as follows: Obviously, from (10), itis easy to deduce that α (cid:107) A ξ − Λ h M ξ (cid:107) = (cid:107) C † C ξ (cid:107) . (11)When α approaches positive infinity, we can obtain the fol-lowing homogeneous linear equation from (11): C † C ξ = . (12)Multiplying both sides of (12) on the left with ξ † , then wecan get (cid:107) C ξ (cid:107) = 0 . As a result, the homogeneous linearequation (5b) is obtained again. Substituting (5b) into (10)gives A ξ = Λ (cid:48) h M ξ . However, if α takes sufficiently large,then this will lead to A + α C † C ≈ α C † C because of the finiteword length in a computer. The result is that the information of matrix A is completely submerged. Consequently, thenumerical accuracy of the eigenvalues associated with thephysical modes in (10) will become very poor [16]. Therefore,we cannot take a sufficiently large penalty parameter α in(10). It is worthwhile that choosing an appropriate penaltyparameter α is important to the penalty method. Webb [17]has studied this problem, and support a method to select thissuitable parameter α .It can be seen that the eigenpair (Λ h , ξ ) of (5) is alwaysthe eigenpair (Λ h , ξ ) of (10). However, the eigenpair (Λ (cid:48) h , ξ ) of (10) is not always the eigenpair of (5). That is to say thatpenalty method will introduce many spurious modes in solving(10). Therefore, this method is not perfect. The penalty methodis applicable to 3-D closed cavity problem under Case 1, 2,3 and 4, except that it cannot remove all the spurious modes.If the eigenvalues in (10) are known, then how to choosethe eigenvalues with physical significance is an importantproblem. Here, we introduce two methods to identify theseeigenvalues associated with physical modes:1) Assuming that (Λ (cid:48) h , ξ ) with (cid:107) ξ (cid:107) = 1 is an eigenpairof (10). The quantity (cid:107) C ξ (cid:107) can be used to identifythe eigenvalues with physical significance. If (cid:107) C ξ (cid:107) is very small, then Λ (cid:48) h is an eigenvalue correspondingto a physical mode. Otherwise, Λ (cid:48) h is an eigenvalueassociated with the spurious mode.2) The eigenvalues corresponding to physical eigenmodeswill not change as the penalty constant α changes, butthe eigenvalues corresponding to spurious modes willchange as the penalty constant α changes. Therefore,under the same mesh, set α = α , α , α , · · · and thensolve (10) repeatedly, finally select the eigenvalues thatremain unchanged in this processing. These unchangedeigenvalues are physical, whereas the changed eigenval-ues are spurious. C. Augmented Method
Consider the generalized eigenvalue problem: (cid:26) A ξ + C † ζ = Λ (cid:48) h M ξ (13a) C ξ = (13b)Obviously, the problem (13) can be rewritten as the followingmatrix form: (cid:20) A C † C O (cid:21) (cid:20) ξζ (cid:21) = Λ (cid:48) h (cid:20) M OO O (cid:21) (cid:20) ξζ (cid:21) . (14)In the numerical linear algebra, the software packageARPACK [15] can be used to solve (14).It is clear that the necessary and sufficient condition for theequivalence of the eigenpair between (5) and (14) is C † ζ = . (15)Assuming that (Λ (cid:48) h , [ ξ ; ζ ]) is the eigenpair of (14), where [ ξ ; ζ ] is the corresponding eigenvector associated with Λ (cid:48) h in (14).Multiplying both sides of (13a) on the left with Y , we obtain ( Y A ) ξ + ( Y C † ) ζ = Λ (cid:48) h ( Y M ) ξ (16) Substituting (9) and into (16) gives Y C † ζ = Λ (cid:48) h C ξ = . (17)If (15) can be derived from (17), then the eigenpair (Λ (cid:48) h , [ ξ ; ζ ]) of (14) is also an eigenpair (Λ (cid:48) h , ξ ) of (5). Conversely, assum-ing that (Λ h , ξ ) is the eigenpair of (5). If we take a vector ζ such that (15) is valid (Obviously, this is achievable), thenthe eigenpair (Λ h , ξ ) of (5) is an eigenpair (Λ h , [ ξ ; ζ ]) of(14). This is to say that each eigenvalue of (5) is always theeigenvalue of (14).According to the above discussion, it can be concluded thatthe necessary and sufficient condition for the equivalence ofthe eigenpair between (5) and (14) is C † ζ = ⇐⇒ Y C † ζ = . (18)The fundamental theory in linear algebra tells us that (18)holds if and only ifrank ( C † ) = rank ( Y C † ) . (19)By substituting (9) into (19), we obtain the necessary and suf-ficient condition for the equivalence of the eigenpair between(5) and (14) isrank ( M † Y † ) = rank ( Y M † Y † ) . (20)We now prove that if µ † = µ > , then (20) holds. In fact,it follows that M † = M > from µ † = µ > . Hence, thereexists an invertible square matrix X such that M = M † = X † X . In terms of the fundamental theory in linear algebra,it is known thatrank ( M † Y † ) = rank ( Y † ) = rank ( XY † )= rank ( Y X † XY † ) = rank ( Y M † Y † ) . (21)The rank equality (21) shows that the rank equality (20) holds.As a result, we have already proved that the eigenpair between(5) and (14) is equivalent provided that µ † = µ > .According to the above discussion, when the material isnot magnetic lossy, i.e., under Case 1 and 2, each eigenpairbetween (5) and (14) is equivalent. In this case, it is easy toshow that each entry of eigenvector ζ in (14) is the same. Next,we briefly prove this conclusion. In fact, since the matrix M is invertible, Y † ζ = can be derived from C † ζ = and C = Y M . By means of (8), ζ = cβ is valid, which showsthat each entry of eigenvector ζ is the same.For all the magnetic lossy materials, we cannot make sureeach eigenvalue of (14) is physical, because we cannot provethat (20) is valid in this case. However, after we carry out manynumerical experiments, these numerical results show that theaugmented method is still free of all the spurious modes forcertain magnetic lossy media. D. Projection Method
In this subsection, we apply the projection method to solveresonant cavity problems under Case 1, 2, 3 and 4. Sincethe divergence-free condition is enforced in this numericalmethod, as a consequence, the projection method can removeall the spurious modes, including spurious zero modes.
It is known that all the solutions to (5b) form a linearsubspace V in C n . Set V = span { q , q , · · · , q r } , where r = dim V and q i ∈ C n . If the matrix M is invertible,then rank ( C ) = rank ( Y ) = m − is valid, which implies r = dim V = n − m + 1 . Set Q = [ q , q , · · · , q r ] ∈ C n × r .The basic idea of solving (5) is that choosing Λ h ∈ C and ξ ∈ V such that ( A ξ − Λ h M ξ ) ⊥ V . (22)This is called the Galerkin condition [18]. Set ξ = Q y , where y ∈ C r . The Galerkin condition (22) can be equivalentlyexpressed the following equation ( Q † AQ ) y = Λ h ( Q † M Q ) y. (23)In order to compute the eigenpair of (23), the matrix Q mustbe given in (23).It is well-known that the popular numerical method offinding several eigenpairs of large-scale sparse matrices is aniterative method, since a fundamental operation in the iterativemethod is the matrix-vector multiplication and this operationis very efficient to the sparse matrix and vector. The efficientprojection method needs to find a sparse matrix Q for the nullspace of the above-mentioned large sparse matrix C . However,since Coleman and Pothen [19] prove that finding the sparsestbasis for the null space of an underdetermined matrix is NP-complete hard, we cannot seek the sparsest matrix Q for thenull space of C .Based on the importance of numerical stability, a set ofnormalized orthogonal bases { q , q , · · · , q r } is used in thisnumerical calculation. In such a case, Q † Q = I r holds, where I r is the identity matrix of order r . Here, we employ singularvalue decomposition (SVD) technique to seek Q . Suppose that C = U DV ∗ (24)is the SVD of the matrix C . We take Q = V (: , n − r + 1 : n ) ,where V (: , n − r + 1 : n ) is a submatrix consisting of the last r columns of V , then CQ = O is valid.If the eigenpair (Λ h , y ) of (23) is obtained by using theimplicitly restarted Arnoldi methods [15], then (Λ h , Q y ) isjust an eigenpair of (5), which corresponds to a physicalnumerical mode of (1). E. The Advantage and Disadvantage among the Above ThreeMethods
The main advantage of penalty method can preserve matrixsize comparing with the augmented method. In addition,penalty method does not destroy the sparsity of the matricescomparing with the projection method. However, the maindisadvantage is that penalty method will introduce spuriousmodes in solving 3-D closed cavity problem. In addition,selecting an appropriate penalty parameter α is not an easywork.The main advantage of augmented method can preserve thesparsity of the matrix. Moreover, when µ † r = µ r > , theaugmented method is free of all the spurious modes. However,the main disadvantage of the augmented method is that the sizeof the matrix has increased. The main advantage of projection method based on SVDcan eliminate all the spurious modes even if the material isboth electric and magnetic lossy. However, the main disadvan-tage of projection method based on SVD is that this methodis not efficient since the matrix Q is usually dense.In a word, these three methods have their own advantagesand disadvantages.IV. N UMERICAL E XPERIMENTS
In this section, we simulate three cavity problems by theabove penalty method, augmented method and projectionmethod. In order to distinguish the numerical eigenvalues asso-ciated with penalty method, augmented method and projectionmethod, the numerical eigenvalues obtained by the penaltymethod, augmented method and projection method are denotedby Λ h ( pe , α ) , Λ h ( au ) and Λ h ( pr ) , respectively. Here α is theparameter in the penalty method and it is usually a positivereal number. It is worthwhile to point out that our adoptedcomputational strategy is serial, instead of parallel. TABLE IT HE N UMERICAL E IGENVALUES ( Λ h , m − ) A SSOCIATED W ITH THE D OMINANT M ODE FROM AN E MPTY S PHERICAL R ESONANT C AVITY AND
CPU T
IME (U NDER C ASE h ( m ) Λ h ( pe , t ( s ) Λ h ( au ) t ( s ) Λ h ( pr ) t ( s ) A. Empty Spherical Resonant Cavity
Let us consider an empty spherical resonator with the radius r = 1 m. The exact eigenvalue associated with the dominantmode is Λ = 7 . m − [20]. Furthermore, the algebraicmultiplicity of the exact eigenvalue Λ is 3. Suppose that the nu-merical eigenvalues Λ (1) h , Λ (2) h and Λ (3) h are the approximationof the exact eigenvalue Λ . Set Λ h = (Λ (1) h +Λ (2) h +Λ (3) h ) / . Weemploy the penalty method, augmented method and projectionmethod to solve this spherical resonant cavity problem, andthen list the numerical eigenvalues Λ h ( pe , α ) , Λ h ( au ) and Λ h ( pr ) in Table I. In order to compare with the efficiencyof these three numerical methods, the CPU time is also givenin Table I.In these three numerical methods, the time and memoryconsumed by the projection method is the largest since thedense matrix Q obtained by SVD is used. This shows thatthe projection method is not efficient. In addition, the CPUtime and memory consumed by the penalty method and theaugmented method are roughly equivalent.In Table I, we can also see that Λ h ( pe , h ( au ) ≈ Λ h ( pr ) ≈ Λ under the finest mesh. This shows that our numer-ical implementations are correct. In this example, we do findthat there are many eigenvalues associated with the spuriousmodes. These numerical eigenvalues are less than Λ provided TABLE IIT HE E IGENVALUES Λ h (au) AND Λ h (pr) ( m − ) WITH P HYSICAL S IGNIFICANCE FROM C YLINDRICAL C AVITY (U NDER C ASE h ( m ) 0 . . . . COMSOL Λ h ( au ) 24 . . j . . j . . j . . j . . j Λ h ( pr ) 24 . . j . . j . . j . . j . . j Λ h ( au ) 26 . . j . . j . . j . . j . . j Λ h ( pr ) 26 . . j . . j . . j . . j . . j Λ h ( au ) 38 . . j . . j . . j . . j . . j Λ h ( pr ) 38 . . j . . j . . j . . j . . jTABLE IIIT HE E IGENVALUES Λ h (pe , , Λ h (au) AND Λ h (pr) ( m − ) WITH P HYSICAL S IGNIFICANCE FROM C YLINDRICAL C AVITY (U NDER C ASE h ( m ) Λ h ( pe , . − . j . − . j . − . j . − . j . − . j Λ h ( au ) 24 . − . j . − . j . − . j . − . j . − . j Λ h ( pr ) 24 . − . j . − . j . − . j . − . j . − . j Λ h ( pe , . − . j . − . j . − . j . − . j . − . j Λ h ( au ) 25 . − . j . − . j . − . j . − . j . − . j Λ h ( pr ) 25 . − . j . − . j . − . j . − . j . − . j I m ag i ne pa r t o f e i gen v a l ue s (23.8807,11.9245) (26.4780,13.2215) Eigenvalues with Physical SignificanceEigenvalues without Physical Significance
Fig. 2. Under the second mesh ( h = 0 . m), the eigenvalues associatedwith physical modes and spurious modes obtained by the penalty method with α = 800 . that the penalty parameter α is less than . However, thenumerical eigenvalues obtained by the augmented method andprojection method are all physical. B. The Resonant Cavity Filled With Electric Lossy Media
In this subsection we consider a cylindrical cavity with theradius r = 0 . m and the height h = 0 . m. Assuming that therelative permittivity and permeability tensor of the medium inthe whole cylindrical cavity are (cid:15) r = − j − j
00 0 2 , µ r = − . j . j . The first three numerical eigenvalues Λ ih ( au ) and Λ ih ( pr ) ( i = 1 , , ) are shown in Table II. In Table II, it can beobserved that Λ ih ( au ) ≈ Λ ih ( pr ) , i = 1 , , , and they coincidewith the eigenvalues corresponding to physical modes fromCOMSOL Multiphysics 5.2a. In the COMSOL simulation, theeigenvalues associated with physical modes are obtained bythe fourth mesh ( h = 0 . m). Notice that there are many spurious zero modes in the numerical results of COMSOLMultiphysics 5.2a. However, there are no any spurious modesin the numerical results of the augmented method and projec-tion method.Under the second mesh ( h = 0 . m), we employ thepenalty method to solve this cylindrical cavity problem, where α = 800 is taken in the numerical calculation, and thenlist the first forty numerical eigenvalues in Fig. 2. In Fig.2, one can see that there are only two eigenvalues withphysical significance, and the rest are all the eigenvalueswithout physical significance, whose imaginary part is zero.Furthermore, these two eigenvalues with physical significanceobtained by the penalty method are equal to the ones obtainedby augmented method. In addition, in (14), we do find thatthe eigenvector ζ associated with each eigenvalue is a vectorconsisting of the same entry. C. The Resonant Cavity Filled With both Electric and Mag-netic Lossy Media
In this subsection we try to find the eigenmode of resonantcavity filled with an electric and magnetic lossy medium. Thepenalty method, augmented method and projection method areused to solve this problem, and then we list the numericaleigenvalues associated with the first two physical modes inTable III.Suppose that the geometric shape of the cavity in thisexample is the same as the one in the example B . In thecylindrical cavity, the relative permittivity and permeabilitytensor of the medium are (cid:15) r = j j
00 0 2 , µ r = − j . j . j − j
00 0 2 . Obviously, the above material is both electric and magneticlossy. Since the exact solution to this problem is unknown, we employ COMSOL Multiphysics 5.2a to simulate this prob-lem, and then obtain the approximate eigenvalues of certainaccuracy. The eigenvalues with physical significance fromCOMSOL are obtained by the fourth mesh ( h = 0 . m).Notice that many spurious zero modes appear in the numericalresults of COMSOL.The numerical eigenvalues from the penalty method andprojection method are listed in Table III. In Table III, we cansee that Λ ih ( pe , ≈ Λ ih ( au ) ≈ Λ ih ( pr ) , i = 1 , , whichcoincide with the eigenvalues corresponding to physical sig-nificance from COMSOL. Here it is worthwhile to emphasizethat the projection method can remove all the spurious modes.V. C ONCLUSION
The finite element method can be applied to solve 3-Dclosed cavity problem filled with anisotropic and nonconduc-tive media. The matrix system resulting from the finite elementmethod is a constrained generalized eigenvalue problem. Thisdifficult problem can be solved by the penalty method, aug-mented method and projection method. The penalty methodcannot remove all the spurious modes. We prove that theaugmented method is free of all the spurious modes if themedium is not magnetic lossy. When the medium is bothelectric and magnetic lossy, the projection method based onSVD technique can deal with this type of resonant cavityproblem very well. However, the projection method basedon SVD technique is not efficient. In future, we would liketo give an efficient iterative method to solve the constrainedgeneralized eigenvalue problem.A
CKNOWLEDGEMENT
We gratefully acknowledge the help of Prof. Qing Huo Liu,who has offered us valuable suggestions in the revision andDr. Yuanguo Zhou for improving our English writing.R
EFERENCES[1] J. F. Lee and R. Mittra, “A note on the application of edge-elementsfor modeling three-dimensional inhomogeneously-filled cavities,”
IEEETrans. Microw. Theory Techn. , vol. 40, no. 9, pp. 1767–1773, Sep. 1992.[2] J. Wang and N. Ida, “Eigenvalue analysis in anisotropically loadedelectromagnetic cavities using ’edge’ finite elements,”
IEEE Trans.Magn. , vol. 28, no. 2, pp. 1438–1441, Mar. 1992.[3] L. Pichon and A. Razek, “Three dimensional resonant mode analysisusing edge elements,”
IEEE Trans. Magn. , vol. 28, no. 2, pp. 1493–1496, Mar. 1992.[4] F. Kikuchi, “Mixed and penalty formulations for finite element analysisof an eigenvalue problem in electromagnetism,”
Comput. Methods Appl.Mech. Eng. , vol. 64, pp. 509–521, 1987.[5] W. Jiang, N. Liu, Y. Yue, and Q. H. Liu, “Mixed finite-element methodfor resonant cavity problem with complex geometric topology andanisotropic lossless media,”
IEEE Trans. Magn. , vol. 52, no. 2, pp. 1–8,Feb. 2016.[6] W. Jiang, X. Wang, and L. Wu, “Mixed finite-element method for 3-Dclosed cavity problem with anisotropic and lossy media,”
IEEE Trans.Magn. , vol. 55, no. 5, pp. 1–6, May 2019.[7] J. Liu, W. Jiang, F. Lin, N. Liu, and Q. H. Liu, “A two-grid vectordiscretization scheme for the resonant cavity problem with anisotropicmedia,”
IEEE Trans. Microw. Theory Techn. , vol. 65, no. 8, pp. 2719–2725, Aug. 2017.[8] W. Jiang, J. Liu, X. Xiong, and Q. H. Liu, “Finite element methodfor resonant cavity problem with complex geometrical structure andanisotropic fully conducting media,”
IEEE Trans. Microw. TheoryTechn. , vol. 65, no. 7, pp. 2240–2248, Jul. 2017. [9] W. Jiang, “Physical DC modes in the microwave resonator with complexgeometric topology,”
IEEE Trans. Magn. , vol. 55, no. 11, pp. 1–7, Nov.2019.[10] W. C. Chew,
Waves and Fields in Inhomogeneous Media . New York:Van Nostrand Reinhold, 1990.[11] R. Hiptmair, “Finite elements in computational electromagnetism,”
ActaNumerica , vol. 11, pp. 237–339, 2002.[12] A. Bossavit, “Mixed finite elements and the complex of whitney forms,” in The Mathematics of Finite Elements and Applications VI. London,U.K.: Academic , pp. 137–144, 1988.[13] R. Geus, “The Jacobi-Davidson algorithm for solving large sparse sym-metric eigenvalue problems with application to the design of acceleratorcavities,” Ph.D. dissertation, ETH Zurich, 2002.[14] D. A. White and J. M. Koning, “Computing solenoidal eigenmodes ofthe vector helmholtz equation: a novel approach,”
IEEE Trans. Magn. ,vol. 38, no. 5, pp. 3420–3425, Sep. 2002.[15] R. Lehoucq, D. Sorensen, and C. Yang,
ARPACK User’s Guide: Solutionof Large-Scale Eigenvalue Problems With Implicitly Restarted ArnoldiMethods . Philadelphia, PA: SIAM, 1998.[16] K. Hayata, M. Eguchi, M. Koshiba, and M. Suzuki, “Vectorial waveanalysis of side-tunnel type polarization-maintaining optical fibers byvariational finite elements,”
J. Lightw. Technol. , vol. 4, no. 8, pp. 1090–1096, Aug. 1986.[17] J. P. Webb, “Efficient generation of divergence-free fields for the finiteelement analysis of 3D cavity resonances,”
IEEE Trans. Magn. , vol. 24,no. 1, pp. 162–165, Jan. 1988.[18] Y. Saad,
Numerical Methods for Large Eigenvalue Problems . Philadel-phia: PA: SIAM, 2011.[19] T. Coleman and A. Pothen, “The null space problem I. Complexity,”
SIAM J. Algebraic Discrete Methods , vol. 7, no. 4, pp. 527–537, 1986.[20] J. M. Jin,
Theory and Computation of Electromagnetic Fields . JohnWiley & Sons, 2011.
Wei Jiang was born in Wuhan, Hubei province,China, in 1984. He received the B.S. degree inmathematics and applied mathematics from HubeiNormal University, Huangshi, China, in 2008, theM.S. degree in computational mathematics fromGuizhou Normal University, Guiyang, China, in2011, and the Ph.D. degree in radio physics fromXiamen University, Xiamen, China, in 2016.From September 2016 to July 2018, he has en-gaged in postdoctoral research at the Institute ofGeophysics and Geomatics, China University ofGeosciences, Wuhan, China. He is currently a lecturer at the School ofMechatronics Engineering, Guizhou Minzu University, Guiyang, China. Untilnow, he has published over 10 papers in refereed journals. His current researchinterests include the finite-element method for partial differential equations,applied numerical algebra, and computational electromagnetics, especially foreigenvalue problems in electromagnetics.