I Introduction
Microwave resonant cavity is an important passive device in microwave engineering. It has many applications in many projects, such as particle accelerator, microwave oven and microwave filter. The microwave resonant cavity problem usually needs to solve the eigenmodes of sourcefree Maxwell’s equations. If the cavity is filled with inhomogeneous media and/or has a complex geometry, then finding its resonant modes by an analytical method is impossible. In order to get the resonant mode, numerical methods must be applied to solve the problem. The main numerical methods include finite element method, finite difference method and boundary element method.
Solving 3D closed cavity problem will introduce many spurious modes if the numerical method cannot preserve the physical property of electromagnetic field. These spurious modes are usually divided into two kinds. One is spurious nonzero mode and the other one is spurious zero mode. In general, the spurious zero mode in numerical results is caused by the neglection of divergencefree condition, and the introduction of the spurious nonzero mode in numerical results is the result of improper discretization method. It is our goal to design a numerical method that can eliminate these two kinds of spurious modes together. We now know edge element method can remove all spurious nonzero modes in solving electromagnetic resonant cavity problems. There are many 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 zero modes because it cannot guarantee the solenoidal properties of electric and magnetic flux densities in the whole cavity.
The main difficulty of solving resonant cavity problems is how to enforce the divergencefree condition given by Gauss’s law in electromagnetics. For the first time, F. Kikuchi [4] introduces a Lagrange multiplier to enforce the divergencefree condition in a weak sense, and propose a mixed finite element method (MFEM) to solve 3D empty cavity problem. As a consequence, Kikuchi’s method is free of all the spurious modes, including spurious zero modes. Based on Kikuchi’s idea, Jiang et al. [5, 6] make use of MFEM to solve 3D resonant cavity problems with anisotropic media, and successfully remove spurious zero and nonzero modes together. However, the MFEM supported in [5, 6] cannot deal with 3D closed cavity problems filled with the anisotropic media, which is both electric and magnetic lossy. On the basis of the reference [5], Liu et al. [7] give a twogrid vector discrete scheme for 3D cavity problems with lossless media. The scheme given in [7] is free of all spurious modes and is very efficient if we just need to know the first few physical modes. Using edge element and linear element, Jiang et al. [8] successfully solve 3D closed cavity problem filled with fully conductive media in the whole cavity. The numerical method given in [8] can also remove the spurious zero and nonzero modes together.
This paper continues to study the microwave cavity problem filled with anisotropic and nonconductive media. After eliminating the electric field in sourcefree Maxwell’s equations, one can get a vector Maxwell’s eigenvalue problem for magnetic field. This problem can be transformed into a corresponding variational formulation by using Green’s formulaes. The edge basis functions of the lowestorder and standard nodal basis functions of linear element are used to discretize the variational formulation. Finally, we need to solve the generalized eigenvalue problem with a linear constraint condition, which is a very difficult problem in numerical linear algebra. In order to overcome this difficult problem, penalty method, augmented method and projection method reduce it to the generalized eigenvalue problem without any constraint. In addition, the advantages and disadvantages among these three computational methods are also given in the paper.
The outline of the paper is as follows. The governing equations and finite element discretization of 3D resonant cavity problem are given in Section II. In Section III, we provide the penalty method, augmented method and projection method to solve the constrained generalized eigenvalue problem and discuss the advantages and disadvantages among these three numerical computational methods. In Section IV, three numerical experiments are carried out to verify our theoretical results.
Ii Finite Element Discretization of 3D Resonant Cavity Problem
Iia Governing Equations for 3D Resonant Cavity Problem
Suppose that is a bounded domain in , is the boundary of and is the outward normal unit vector on . Let and
be the permeability and permittivity in vacuum, respectively. The relative permeability and permittivity tensor of an anisotropic medium are denoted by
and , respectively. The angular frequency of electromagnetic wave is denoted by .The governing equations of 3D closed cavity problem are the sourcefree Maxwell’s equations of the firstorder. After eliminating the electric field
in the sourcefree Maxwell’s equations, one can get a secondorder partial differential equations (PDEs) for the magnetic field
:(1a)  
(1b)  
(1c)  
(1d) 
where is the square of the wavenumber in vacuum. We would like to seek with such that PDEs (1) holds. In electromagnetics, with is called a physical resonant eigenmode in the resonant cavity. In mathematics, with is called an eigenpair of PDEs (1), and and
are called the eigenvalue and eigenfunction of PDEs (
1), respectively. In addition, we know PDEs (1) only has a discrete point spectrum. Note that there may be several zero eigenmodes in PDEs (1), for details, please see [9].If the anisotropic material is lossless, then and are both Hermitian [10], that is and , where the superscript denotes the conjugate transposition. Moreover, when the anisotropic material is lossless, assuming that and are Hermitian positive definite since the anisotropic and lossless material in nature usually has this property. For the sake of simplicity, a Hermitian positive definite matrix is denoted by . In terms of the lossy characteristics of the anisotropic and nonconductive material, it is usually divided into the following four categories:

Case 1: and . The medium is lossless.

Case 2: and . The medium is electric lossy, but is not magnetic lossy.

Case 3: and . The medium is magnetic lossy, but is not electric lossy.

Case 4: and . The medium is both electric and magnetic lossy.
Under Case 1, 2 and 3, the MFEM can deal with these types of the resonant cavity problems very well, and it is free of all the spurious modes in numerical results [5, 6]. However, the MFEM is not suitable for 3D closed cavity problem under Case 4, because it is difficult to propose an appropriate mixed variational formulation. For the 3D closed cavity problem under Case 4, the projection method introduced in this paper can deal with this problem very well, and the projection method can remove all the spurious modes.
IiB Finite Element Discretization
It is wellknown that the finite element method is a variational method, and only operates on the weak form of PDE, instead of the strong form of PDE. Hence, the corresponding weak form associated with PDEs (1) is given in the subsection. To get this weak form, we introduce the following Hilbert spaces over :
Define the continuous sesquilinear forms:
where the symbol stands for the complex conjugation of a given complexvalued function.
Using the Green’s formulas, the weak form of PDEs (1) reads as:
(2a)  
(2b) 
Under Case 1, the eigenvalues are made up of countable nonnegative real numbers. Under Case 2, 3 and 4, the eigenvalues are made up of countable complex numbers. The physical interpretation is such a physical fact that there is no electromagnetic energy loss in the resonant cavity if the material is lossless and there exists electromagnetic energy loss in the resonant cavity provided that the material has a dielectric loss.
We now consider the conforming finite element discretization of (2). Let be a regular tetrahedral mesh of the cavity . Here is the length of the longest edge in tetrahedral mesh . As usual, the edge element space of the lowestorder under the mesh is used to approximate the Hilbert space and standard linear element space under the mesh is used to approximate the Hilbert space . From [11], we know and .
The linear element space can be written as:
where are four local nodal basis functions on the tetrahedral element and of the form , where are four constants. These four local basis functions are defined on the four vertices of the tetrahedral element .
The edge element space of the lowestorder can be written as:
where are six local edge basis functions on the tetrahedral element and is of the form , where and are two constant vectors and is the position vector. The concrete expressions of are as follows:
These six local vector basis functions are defined on the six edges of the tetrahedral element . In Fig. 1, we give a local nodal numbering in the tetrahedral element , and specify the local reference direction for each edge in .
Restricting (2) on , we get the discrete variational formulation associated with (2):
(3a)  
(3b) 
Here and are an approximation of the exact eigenvalue and exact eigenfunction in (2), respectively.
Suppose that , where is the th global nodal basis function associated with the node and the integer is the number of the total nodes in . Assuming that , where is the th global edge basis function associated with the edge and the integer is the number of the total edges in . Since , then it can be written as
(4) 
Finally, the discrete variational formulation (3) can be reduced to the following generalized eigenvalue problem with a linear constraint:
(5a)  
(5b) 
where
Solving (5a) directly and ignoring (5b) will introduce a lot of spurious zero modes. In order to eliminate these spurious zero modes, we must enforce the constraint (5b) in solving (5a). How to enforce the constraint (5b) in solving (5a) is a key problem. In next section, we shall deal with this troublesome problem. Once the eigenpair is found from (5), then the numerical eigenvalue in (3) is given by and the corresponding numerical eigenfunction in (3) is given by (4).
In the community of numerical linear algebra, the problem (5) is a constrained generalized eigenvalue problem. Obviously, its numerical computation is much more difficult than that of the generalized eigenvalue problem without any constraint. In Section III, the problem (5) is reduced to three types of generalized eigenvalue problem without any constraint. It is important to point out if there is no relation among these three matrices , and , then the constrained generalized eigenvalue problem (5) may have no solution.
Iii Three Numerical Solvers of the Constrained Generalized Eigenvalue Problem (5)
In this section, we first try to give a relation among the matrices , and , and then support three numerical computational methods of solving (5). They are penalty method, augmented method and projection method, respectively.
Iiia The Relation Among the Matrices , and
Let and be the sets consisting of all nodes and edges in respectively, where and are the global labels of all nodes and edges in , respectively. Note that the global vector basis function has a local direction associated with the edge . If two vertices of the edge are and , then we state that the direction of is from the node to the node , where .
In accordance with deRhamcomplex [12], holds, where . This implies that
(6) 
The above formula (6) is also introduced in [13] and [14]. It is easy to know since . Set
where and are th row and column vector in the matrix , respectively.
In fact, the each entry in is , or , and is quite sparse. Let us consider the formula (6) under the case of an arbitrary tetrahedral element in (see Fig. 1). It is easy to verify that the following formulas are valid:
Here we need to make use of the relation . Consider the th row in the matrix , which is related to the node . If is not a vertex of the th edge , then . Let us recall the basic concept of degree of a vertex in graph theory. The degree of a vertex is defined by the number of edges connecting it. The number of the nonzero entries of is equal to the degree of . Assuming that the degree of is and have the common vertex . If the direction of points to , then , otherwise . Consider the th column vector in the matrix , which is related to the edge . If is not a vertex of the edge , then . Since each edge in a tetrahedral mesh only has two vertices, the column vector only has two nonzero entries. Assuming that the initial and terminal points of the edge are and respectively, then and , where . Obviously, the sparse matrix can be easily obtained by the mesh data in . The sparse matrix is usually called the directional connectivity matrix associated with a tetrahedral mesh .
The sum of all entries in the column vector is equal to zero. This implies that the homogenous linear equation
(7) 
has a special nonzero solution . An easy computation shows that . Therefore all solutions of (7) form a linear space of one dimension, that is
(8) 
where Null stands for taking the nullspace of a matrix.
The matrix builds a relation among , and . It can be proved that the following matrix identities are valid:
(9) 
where is a null matrix. In fact, for and , we have
where is the entry at th row and th column of the matrix . It is important to emphasize that the matrix identities (9) are always valid for the medium under Case 1, 2, 3 and 4. Hence, we can obtain the matrix by using the sparse matrices and , instead of the calculation by using the continuous sesquilinear forms directly.
If the eigenvalue is nonzero in (5a), then (5b) can be deduced from (5a). As a matter of fact, can be derived from (5a), and then we get by (9), which is just (5b). It is worthwhile to point out the number of zero eigenvalues in (5a) is equal to , and these zero eigenvalues are all spurious. The dimension of the linear space is equal to , which shows that the larger the number of mesh nodes is, the more the number of spurious zero modes in (5a) is. Based on this reason, we need to remove these spurious zero modes by introducing an appropriate numerical method.
IiiB Penalty Method
Consider the following generalized eigenvalue problem:
(10) 
where is a penalty constant, which is required the user to set and is the Euclidean norm of a given vector . The problem (10) is a generalized eigenvalue problem without any constraint, and it can be solved by the numerical software package ARPACK [15].
The parameter is usually set to a large positive real number. The reason is as follows: Obviously, from (10), it is easy to deduce that
(11) 
When approaches positive infinity, we can obtain the following homogeneous linear equation from (11):
(12) 
Multiplying both sides of (12) on the left with , then we can get . As a result, the homogeneous linear equation (5b) is obtained again. Substituting (5b) into (10) gives . However, if takes sufficiently large, then this will lead to because of the finite word length in a computer. The result is that the information of matrix is completely submerged. Consequently, the numerical accuracy of the eigenvalues associated with the physical 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 penalty parameter is important to the penalty method. Webb [17] has studied this problem, and support a method to select this suitable parameter .
It can be seen that the eigenpair of (5) is always the eigenpair of (10). However, the eigenpair of (10) is not always the eigenpair of (5). That is to say that penalty method will introduce many spurious modes in solving (10). Therefore, this method is not perfect. The penalty method is applicable to 3D 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 choose the eigenvalues with physical significance is an important problem. Here, we introduce two methods to identify these eigenvalues associated with physical modes:

Assuming that with is an eigenpair of (10). The quantity can be used to identify the eigenvalues with physical significance. If is very small, then is an eigenvalue corresponding to a physical mode. Otherwise, is an eigenvalue associated with the spurious mode.

The eigenvalues corresponding to physical eigenmodes will not change as the penalty constant changes, but the eigenvalues corresponding to spurious modes will change as the penalty constant changes. Therefore, under the same mesh, set and then solve (10) repeatedly, finally select the eigenvalues that remain unchanged in this processing. These unchanged eigenvalues are physical, whereas the changed eigenvalues are spurious.
IiiC Augmented Method
Consider the generalized eigenvalue problem:
(13a)  
(13b) 
Obviously, the problem (13) can be rewritten as the following matrix form:
(14) 
In the numerical linear algebra, the software package ARPACK [15] can be used to solve (14).
It is clear that the necessary and sufficient condition for the equivalence of the eigenpair between (5) and (14) is
(15) 
Assuming that is the eigenpair of (14), where
is the corresponding eigenvector associated with
in (14). Multiplying both sides of (13a) on the left with , we obtain(16) 
Substituting (9) and into (16) gives
(17) 
If (15) can be derived from (17), then the eigenpair of (14) is also an eigenpair of (5). Conversely, assuming that is the eigenpair of (5). If we take a vector such that (15) is valid (Obviously, this is achievable), then the eigenpair of (5) is an eigenpair of (14). This is to say that each eigenvalue of (5) is always the eigenvalue of (14).
According to the above discussion, it can be concluded that the necessary and sufficient condition for the equivalence of the eigenpair between (5) and (14) is
(18) 
The fundamental theory in linear algebra tells us that (18) holds if and only if
(19) 
By substituting (9) into (19), we obtain the necessary and sufficient condition for the equivalence of the eigenpair between (5) and (14) is
(20) 
We now prove that if , then (20) holds. In fact, it follows that from . Hence, there exists an invertible square matrix such that . In terms of the fundamental theory in linear algebra, it is known that
(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 is not magnetic lossy, i.e., under Case 1 and 2, each eigenpair between (5) and (14) is equivalent. In this case, it is easy to show that each entry of eigenvector in (14) is the same. Next, we briefly prove this conclusion. In fact, since the matrix is invertible, can be derived from and . By means of (8), is valid, which shows that each entry of eigenvector is the same.
For all the magnetic lossy materials, we cannot make sure each eigenvalue of (14) is physical, because we cannot prove that (20) is valid in this case. However, after we carry out many numerical experiments, these numerical results show that the augmented method is still free of all the spurious modes for certain magnetic lossy media.
IiiD Projection Method
In this subsection, we apply the projection method to solve resonant cavity problems under Case 1, 2, 3 and 4. Since the divergencefree condition is enforced in this numerical method, as a consequence, the projection method can remove all the spurious modes, including spurious zero modes.
It is known that all the solutions to (5b) form a linear subspace in . Set , where and . If the matrix is invertible, then is valid, which implies . Set .
The basic idea of solving (5) is that choosing and such that
(22) 
This is called the Galerkin condition [18]. Set , where . The Galerkin condition (22) can be equivalently expressed the following equation
(23) 
In order to compute the eigenpair of (23), the matrix must be given in (23).
It is wellknown that the popular numerical method of finding several eigenpairs of largescale sparse matrices is an iterative method, since a fundamental operation in the iterative method is the matrixvector multiplication and this operation is very efficient to the sparse matrix and vector. The efficient projection method needs to find a sparse matrix for the null space of the abovementioned large sparse matrix . However, since Coleman and Pothen [19] prove that finding the sparsest basis for the null space of an underdetermined matrix is NPcomplete hard, we cannot seek the sparsest matrix for the null space of .
Based on the importance of numerical stability, a set of normalized orthogonal bases is used in this numerical calculation. In such a case, holds, where
is the identity matrix of order
. Here, we employ singular value decomposition (SVD) technique to seek . Suppose that(24) 
is the SVD of the matrix . We take , where is a submatrix consisting of the last columns of , then is valid.
IiiE The Advantage and Disadvantage among the Above Three Methods
The main advantage of penalty method can preserve matrix size comparing with the augmented method. In addition, penalty method does not destroy the sparsity of the matrices comparing with the projection method. However, the main disadvantage is that penalty method will introduce spurious modes in solving 3D closed cavity problem. In addition, selecting an appropriate penalty parameter is not an easy work.
The main advantage of augmented method can preserve the sparsity of the matrix. Moreover, when , the augmented method is free of all the spurious modes. However, the main disadvantage of the augmented method is that the size of the matrix has increased.
The main advantage of projection method based on SVD can eliminate all the spurious modes even if the material is both electric and magnetic lossy. However, the main disadvantage of projection method based on SVD is that this method is not efficient since the matrix is usually dense.
In a word, these three methods have their own advantages and disadvantages.
Iv Numerical Experiments
In this section, we simulate three cavity problems by the above penalty method, augmented method and projection method. In order to distinguish the numerical eigenvalues associated with penalty method, augmented method and projection method, the numerical eigenvalues obtained by the penalty method, augmented method and projection method are denoted by , and , respectively. Here is the parameter in the penalty method and it is usually a positive real number. It is worthwhile to point out that our adopted computational strategy is serial, instead of parallel.
0.38493  0.27062  0.22416  0.16258  Exact  

7.71147  7.62386  7.59006  7.55655  7.52793  
10.3  18.5  28.7  40.6  –  
7.71147  7.62386  7.59006  7.55655  7.52793  
12.7  21.8  32.6  60.7  –  
7.71147  7.62386  7.59006  7.55654  7.52793  
60.5  100.7  180.9  360.8  – 
Iva Empty Spherical Resonant Cavity
Let us consider an empty spherical resonator with the radius m. The exact eigenvalue associated with the dominant mode is [20]. Furthermore, the algebraic multiplicity of the exact eigenvalue is 3. Suppose that the numerical eigenvalues , and are the approximation of the exact eigenvalue . Set . We employ the penalty method, augmented method and projection method to solve this spherical resonant cavity problem, and then list the numerical eigenvalues , and in Table I. In order to compare with the efficiency of these three numerical methods, the CPU time is also given in Table I.
COMSOL  

0.1043  0.0714  0.0580  0.0428  COMSOL  

Comments
There are no comments yet.