    # Three Numerical Eigensolvers for 3-D Cavity Resonators Filled With Anisotropic and Nonconductive Media

This paper mainly investigates the classic resonant cavity problem with anisotropic and nonconductive media, which is a linear vector Maxwell's eigenvalue problem. The finite element method based on edge element of the lowest-order and standard linear element is used to solve this type of 3-D closed cavity problem. In order to eliminate spurious zero modes in the numerical simulation, the divergence-free condition supported by Gauss' law is enforced in a weak sense. After the finite element discretization, the generalized eigenvalue problem with a linear constraint condition needs to be solved. Penalty method, augmented method and projection method are applied to solve this difficult problem in numerical linear algebra. The advantages and disadvantages of these three computational methods are also given in this paper. Furthermore, we prove that the augmented method is free of spurious modes as long as the anisotropic material is not magnetic lossy. The projection method based on singular value decomposition technique can be used to solve the resonant cavity problem. Moreover, the projection method cannot introduce any spurious modes. At last, several numerical experiments are carried out to verify our theoretical results.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 source-free 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 3-D 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 divergence-free 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, Wang and Ida , Pichon and Razek , 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 divergence-free condition given by Gauss’s law in electromagnetics. For the first time, F. Kikuchi  introduces a Lagrange multiplier to enforce the divergence-free condition in a weak sense, and propose a mixed finite element method (MFEM) to solve 3-D 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 3-D 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 3-D closed cavity problems filled with the anisotropic media, which is both electric and magnetic lossy. On the basis of the reference , Liu et al.  give a two-grid vector discrete scheme for 3-D cavity problems with lossless media. The scheme given in  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.  successfully solve 3-D closed cavity problem filled with fully conductive media in the whole cavity. The numerical method given in  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 source-free 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 lowest-order 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 3-D 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 3-D Resonant Cavity Problem

### Ii-a Governing Equations for 3-D 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 3-D closed cavity problem are the source-free Maxwell’s equations of the first-order. After eliminating the electric field

in the source-free Maxwell’s equations, one can get a second-order partial differential equations (PDEs) for the magnetic field

:

 ∇×(¯¯¯¯ϵ−1r∇×H)=Λ¯¯¯¯¯¯μrH  in Ω (1a) ∇⋅(¯¯¯¯¯¯μrH)=0  in Ω (1b) ^n×(¯¯¯¯ϵ−1r∇×H)=0  on ∂Ω (1c) ^n⋅(¯¯¯¯¯¯μrH)=0  on ∂Ω, (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 .

If the anisotropic material is lossless, then and are both Hermitian , 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:

1. Case 1: and . The medium is lossless.

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

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

4. 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 3-D closed cavity problem under Case 4, because it is difficult to propose an appropriate mixed variational formulation. For the 3-D 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.

### Ii-B Finite Element Discretization

It is well-known 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 :

 L2(Ω)={f:∫Ω|f(x,y,z)|2dxdydz<+∞} H1(Ω)={f∈L2(Ω):∇f∈(L2(Ω))3} H(curl,Ω)={F∈(L2(Ω))3:∇×F∈(L2(Ω))3}.

Define the continuous sesquilinear forms:

 A: H(curl,Ω)×H(curl,Ω)→C (H,F)→∫Ω¯¯¯¯¯¯ϵ−1r∇×H⋅∇×F∗dΩ M: (L2(Ω))3×(L2(Ω))3→C (H,F)→∫Ω¯¯¯¯¯¯¯μrH⋅F∗dΩ C: H(curl,Ω)×H1(Ω)→C (H,q)→∫Ω¯¯¯¯¯¯μrH⋅∇q∗dΩ

where the symbol stands for the complex conjugation of a given complex-valued function.

Using the Green’s formulas, the weak form of PDEs (1) reads as:

 Seek~{}Λ∈C, H∈H(curl,Ω), H≠0~{}such that A(H,F)=ΛM(H,F), ∀ F∈H(curl,Ω) (2a) C(H,q)=0, ∀ q∈H1(Ω) (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 lowest-order 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 , we know and .

The linear element space can be written as:

 Sh={ϕ:ϕ|K∈span{LK1,LK2,LK3,LK4}}

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 lowest-order can be written as:

 Wh={F:F|K∈span{{% \boldmathN}K1,{\boldmathN}K2,⋯,{% \boldmathN}K6}}

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:

 {\boldmathN}K1=LK1∇LK2−LK2∇LK1, {\boldmathN}K2=LK2∇LK3−LK3∇LK2 {\boldmathN}K3=LK1∇LK3−LK3∇LK1, {\boldmathN}K4=LK3∇LK4−LK4∇LK3 {\boldmathN}K5=LK1∇LK4−LK4∇LK1, {\boldmathN}K6=LK2∇LK4−LK4∇LK2

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 . Fig. 1: Local nodal numbering for the element K and the local reference direction for the edge are chosen by means of local nodal numbering.

Restricting (2) on , we get the discrete variational formulation associated with (2):

 Seek~{}Λh∈C, Hh∈Wh, Hh≠0~{}such that A(Hh,F)=ΛhM(Hh,F), ∀ F∈Wh (3a) C(Hh,q)=0, ∀ q∈Sh (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

 Hh=n∑i=1ξi{\boldmathN}i. (4)

Finally, the discrete variational formulation (3) can be reduced to the following generalized eigenvalue problem with a linear constraint:

 {\boldmathA}ξ=Λh{\boldmathM}ξ (5a) {\boldmathC}ξ=0 (5b)

where

 {\boldmathA}=(aik)∈Cn×n, {\boldmathM}=(mik)∈Cn×n, {\boldmathC}=(cik)∈Cm×n, ξ=[ξ1,  ξ2,⋯,ξn]T∈Cn, aik=A(Nk,Ni), mik=M(Nk,Ni), cik=C(Nk,Li).

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.

### Iii-a The Relation Among the Matrices A, M and C

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 deRham-complex , holds, where . This implies that

 ∇Li=n∑k=1yikNk,  i=1,2,⋯,m. (6)

The above formula (6) is also introduced in  and . It is easy to know since . Set

 {\boldmathY}=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣y11y12⋯y1ny21y22⋯y2n⋮⋮⋮⋮ym1ym2⋯ymn⎤⎥ ⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣{\boldmathy}% 1{\boldmathy}2⋮{\boldmathy}m⎤⎥ ⎥ ⎥ ⎥ ⎥⎦=[{\boldmathd}1,{% \boldmathd}2,⋯,{\boldmathd}n],

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:

 ∇LK1=−{\boldmathN}K1−{% \boldmathN}K3−{\boldmathN}K5,  ∇LK2={\boldmathN}K1−{\boldmathN}K2−{% \boldmathN}K6 ∇LK3={\boldmathN}K2+{\boldmath% N}K3−{\boldmathN}K4,  ∇LK4={% \boldmathN}K4+{\boldmathN}K5+{\boldmathN% }K6

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

 {\boldmathY}†ζ=0 (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

 Null({\boldmathY}†)=span{β}, (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:

 {\boldmathY}{\boldmathA}={\boldmathO},  {\boldmathC}={\boldmathY}{\boldmathM}, (9)

where is a null matrix. In fact, for and , we have

 ({\boldmathY}{\boldmathA})il = n∑k=1yikakl=n∑k=1yikA(Nl,Nk) = A(Nl,n∑k=1yikNk)=A(Nl,∇Li)=0, ({\boldmathY}{\boldmathM})il = n∑k=1yikmkl=n∑k=1yikM(Nl,Nk) = M(Nl,n∑k=1yikNk)=M(Nl,∇Li) = C(Nl,Li)=cil,

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.

### Iii-B Penalty Method

Consider the following generalized eigenvalue problem:

 ({\boldmathA}+α{\boldmathC}†{% \boldmathC})ξ=Λ′h{\boldmathM}ξ,  ∥ξ∥2=1 (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 .

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

 1α∥{\boldmathA}ξ−Λh{\boldmathM}% ξ∥2=∥{\boldmathC}†{\boldmathC}ξ∥2. (11)

When approaches positive infinity, we can obtain the following homogeneous linear equation from (11):

 {\boldmathC}†{\boldmathC}ξ=0. (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 . 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  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 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 choose the eigenvalues with physical significance is an important problem. Here, we introduce two methods to identify these eigenvalues associated with physical modes:

1. 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.

2. 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.

### Iii-C Augmented Method

Consider the generalized eigenvalue problem:

 {\boldmathA}ξ+{\boldmathC}†ζ=Λ′h{\boldmathM}ξ (13a) {\boldmathC}ξ=0 (13b)

Obviously, the problem (13) can be rewritten as the following matrix form:

 [{\boldmathA}{\boldmathC}†{\boldmathC}{\boldmathO}][ξζ]=Λ′h[{% \boldmathM}{\boldmathO}{\boldmathO}{\boldmathO}][ξζ]. (14)

In the numerical linear algebra, the software package ARPACK  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

 {\boldmathC}†ζ=0. (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

 ({\boldmathY}{\boldmathA})ξ+({\boldmathY}{\boldmathC}†)ζ=Λ′h({% \boldmathY}{\boldmathM})ξ (16)

Substituting (9) and into (16) gives

 {\boldmathY}{\boldmathC}†ζ=Λ′h{\boldmathC}ξ=0. (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

 {\boldmathC}†ζ=0⟺{% \boldmathY}{\boldmathC}†ζ=0. (18)

The fundamental theory in linear algebra tells us that (18) holds if and only if

 rank({\boldmathC}†)=rank({% \boldmathY}{\boldmathC}†). (19)

By substituting (9) into (19), we obtain the necessary and sufficient condition for the equivalence of the eigenpair between (5) and (14) is

 rank({\boldmathM}†{\boldmathY}†)=rank({\boldmathY}{\boldmathM}†{\boldmathY}†). (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

 rank({\boldmathM}†{\boldmath% Y}†)=rank({\boldmathY}†)=rank({\boldmathX}{\boldmathY}†) =rank({\boldmathY}{\boldmathX}†{\boldmathX}{\boldmathY}†)=% rank({\boldmathY}{\boldmathM}†{% \boldmathY}†). (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.

### Iii-D Projection Method

In this subsection, we apply the projection method to solve resonant cavity problems under Case 1, 2, 3 and 4. Since the divergence-free 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

 ({\boldmathA}ξ−Λh{\boldmathM}ξ)⊥V. (22)

This is called the Galerkin condition . Set , where . The Galerkin condition (22) can be equivalently expressed the following equation

 ({\boldmathQ}†{\boldmathA}{\boldmathQ})y=Λh({\boldmathQ}†{\boldmathM}{\boldmathQ})y. (23)

In order to compute the eigenpair of (23), the matrix must be given in (23).

It is well-known that the popular numerical method of finding several eigenpairs of large-scale sparse matrices is an iterative method, since a fundamental operation in the iterative method is the matrix-vector 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 above-mentioned large sparse matrix . However, since Coleman and Pothen  prove that finding the sparsest basis for the null space of an underdetermined matrix is NP-complete 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

 {\boldmathC}={\boldmathU}{\boldmathD}{% \boldmathV}∗ (24)

is the SVD of the matrix . We take , where is a submatrix consisting of the last columns of , then is valid.

If the eigenpair of (23) is obtained by using the implicitly restarted Arnoldi methods , then is just an eigenpair of (5), which corresponds to a physical numerical mode of (1).

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 3-D 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.

### Iv-a Empty Spherical Resonant Cavity

Let us consider an empty spherical resonator with the radius  m. The exact eigenvalue associated with the dominant mode is . 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.