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

Suppose that $\Omega$ is a bounded domain in $\mathbb{R}^{3}$, $\partial\Omega$ is the boundary of $\Omega$ and $\hat{n}$ is the outward normal unit vector on $\partial\Omega$. Let $\epsilon_{0}$ and $\mu_{0}$ be the permeability and permittivity in vacuum, respectively. The relative permeability and permittivity tensor of an anisotropic medium are denoted by $\overline{\overline{\mu}}_{r}$ and $\overline{\overline{\epsilon}}_{r}$, respectively. The angular frequency of electromagnetic wave is denoted by $\omega$.

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 ${{\bf E}}$ in the source-free Maxwell’s equations, one can get a second-order partial differential equations (PDEs) for the magnetic field ${{\bf H}}$:

where $\Lambda=\omega^{2}\epsilon_{0}\mu_{0}$ is the square of the wavenumber in vacuum. We would like to seek $(\Lambda,{{\bf H}})$ with ${{\bf H}}\neq{\bf{0}}$ such that PDEs (1) holds. In electromagnetics, $(\Lambda,{{\bf H}})$ with ${{\bf H}}\neq{\bf{0}}$ is called a physical resonant eigenmode in the resonant cavity. In mathematics, $(\Lambda,{{\bf H}})$ with ${{\bf H}}\neq{\bf{0}}$ is called an eigenpair of PDEs (1), and $\Lambda$ and ${{\bf H}}$ 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 $\overline{\overline{\mu}}_{r}$ and $\overline{\overline{\epsilon}}_{r}$ are both Hermitian [10], that is $\overline{\overline{\mu}}_{r}^{{\dagger}}=\overline{\overline{\mu}}_{r}$ and $\overline{\overline{\epsilon}}_{r}^{{\dagger}}=\overline{\overline{\epsilon}}_{r}$, where the superscript ${{\dagger}}$ denotes the conjugate transposition. Moreover, when the anisotropic material is lossless, assuming that $\overline{\overline{\epsilon}}_{r}$ and $\overline{\overline{\mu}}_{r}$ 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 ${M}$ is denoted by $\mbox{{\boldmath${M}$}}^{{\dagger}}=\mbox{{\boldmath${M}$}}>0$. In terms of the lossy characteristics of the anisotropic and nonconductive material, it is usually divided into the following four categories:

1. 1.

   Case 1: $\overline{\overline{\epsilon}}_{r}^{{\dagger}}=\overline{\overline{\epsilon}}_{r}>0$ and $\overline{\overline{\mu}}_{r}^{{\dagger}}=\overline{\overline{\mu}}_{r}>0$. The medium is lossless.

2. 2.

   Case 2: $\overline{\overline{\epsilon}}_{r}^{{\dagger}}\neq\overline{\overline{\epsilon}}_{r}$ and $\overline{\overline{\mu}}_{r}^{{\dagger}}=\overline{\overline{\mu}}_{r}>0$. The medium is electric lossy, but is not magnetic lossy.

3. 3.

   Case 3: $\overline{\overline{\epsilon}}_{r}^{{\dagger}}=\overline{\overline{\epsilon}}_{r}>0$ and $\overline{\overline{\mu}}_{r}^{{\dagger}}\neq\overline{\overline{\mu}}_{r}$. The medium is magnetic lossy, but is not electric lossy.

4. 4.

   Case 4: $\overline{\overline{\epsilon}}_{r}^{{\dagger}}\neq\overline{\overline{\epsilon}}_{r}$ and $\overline{\overline{\mu}}_{r}^{{\dagger}}\neq\overline{\overline{\mu}}_{r}$. 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.

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 $\mathbb{C}$:

Define the continuous sesquilinear forms:

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:

Under Case 1, the eigenvalues $\Lambda$ are made up of countable nonnegative real numbers. Under Case 2, 3 and 4, the eigenvalues $\Lambda$ 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 $\mathcal{T}_{h}$ be a regular tetrahedral mesh of the cavity $\Omega$. Here $h$ is the length of the longest edge in tetrahedral mesh $\mathcal{T}_{h}$. As usual, the edge element space ${{\bf W}}^{h}$ of the lowest-order under the mesh $\mathcal{T}_{h}$ is used to approximate the Hilbert space ${{\bf H}}(\mbox{curl},\Omega)$ and standard linear element space $S^{h}$ under the mesh $\mathcal{T}_{h}$ is used to approximate the Hilbert space $H^{1}(\Omega)$. From [11], we know $S^{h}\subsetneq H^{1}(\Omega)$ and ${{\bf W}}^{h}\subsetneq{{\bf H}}(\mbox{curl},\Omega)$.

The linear element space $S^{h}$ can be written as:

where $L_{i}^{K}~{}(i=1,2,3,4)$ are four local nodal basis functions on 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 four local basis functions are defined on the four vertices of the tetrahedral element $K$.

The edge element space ${{\bf W}}^{h}$ of the lowest-order can be written as:

where $\mbox{{\boldmath${N}$}}_{i}^{K}~{}(i=1,2,\cdots,6)$ are six local edge basis functions on the tetrahedral element $K$ and $\mbox{{\boldmath${N}$}}_{i}^{K}$ is of the form $\vec{\alpha}_{i}+\vec{\beta}_{i}\times\vec{r}$, where $\vec{\alpha}_{i}$ and $\vec{\beta}_{i}$ are two constant vectors and $\vec{r}$ is the position vector. The concrete expressions of $\mbox{{\boldmath${N}$}}_{i}^{K}~{}(i=1,2,\cdots,6)$ are as follows:

These six local vector basis functions are defined on the six edges of the tetrahedral element $K$. In Fig. 1, we give a local nodal numbering in the tetrahedral element $K$, and specify the local reference direction for each edge in $K$.

Restricting (2) on ${{\bf W}}^{h}\times S^{h}$, we get the discrete variational formulation associated with (2):

Here $\Lambda_{h}$ and ${{\bf H}}_{h}$ are an approximation of the exact eigenvalue $\Lambda$ and exact eigenfunction ${{\bf H}}$ in (2), respectively.

Suppose that $S^{h}=\textrm{span}\big{\{}L_{i}\big{\}}_{i=1}^{m}$, where $L_{i}$ is the $i$-th global nodal basis function associated with the node $i$ and the integer $m$ is the number of the total nodes in $\mathcal{T}_{h}$. Assuming that ${{\bf W}}^{h}=\textrm{span}\big{\{}\mbox{{\boldmath${N}$}}_{i}\big{\}}_{i=1}^{n}$, where $\mbox{{\boldmath${N}$}}_{i}$ is the $i$-th global edge basis function associated with the edge $i$ and the integer $n$ is the number of the total edges in $\mathcal{T}_{h}$. Since ${{\bf H}}_{h}\in{{{\bf W}}^{h}}$, then it can be written as

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

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 $(\Lambda_{h},\xi)$ is found from (5), then the numerical eigenvalue in (3) is given by $\Lambda_{h}$ and the corresponding numerical eigenfunction ${{\bf H}}_{h}$ 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 ${A}$, ${M}$ and ${C}$, then the constrained generalized eigenvalue problem (5) may have no solution.

Let $\{A_{1},A_{2},\cdots,A_{m-1},A_{m}\}$ and $\{e_{1},e_{2},\cdots,e_{n-1},e_{n}\}$ be the sets consisting of all nodes and edges in $\mathcal{T}_{h}$ respectively, where $1,2,\cdots,m$ and $1,2,\cdots,n$ are the global labels of all nodes and edges in $\mathcal{T}_{h}$, respectively. Note that the global vector basis function $\mbox{{\boldmath${N}$}}_{i}$ has a local direction associated with the edge $e_{i}$. If two vertices of the edge $e_{i}$ are $A_{i_{1}}$ and $A_{i_{2}}$, then we state that the direction of $\mbox{{\boldmath${N}$}}_{i}$ is from the node $A_{i_{1}}$ to the node $A_{i_{2}}$, where $i_{1}<i_{2}$.

In accordance with deRham-complex [12], $\nabla{S^{h}}\subsetneq{{\bf W}}^{h}$ holds, where $\nabla{S^{h}}=\{\nabla f:\,\forall f\in{S^{h}}\}$. This implies that

The above formula (6) is also introduced in [13] and [14]. It is easy to know $\dim(\nabla{S^{h}})=m-1$ since $\sum_{i=1}^{m}\nabla{L_{i}}={\bf{0}}$. Set

where $\mbox{{\boldmath${y}$}}_{k}$ and $\mbox{{\boldmath${d}$}}_{k}$ are $k$-th row and column vector in the matrix ${Y}$, respectively.

In fact, the each entry $y_{ki}$ in ${Y}$ is $-1$, $1$ or $0$, and ${Y}$ is quite sparse. Let us consider the formula (6) under the case of an arbitrary tetrahedral element $K$ in $\mathcal{T}_{h}$ (see Fig. 1). It is easy to verify that the following formulas are valid:

Here we need to make use of the relation $\sum_{i=1}^{4}L_{i}^{K}=1$. Consider the $k$-th row $\mbox{{\boldmath${y}$}}_{k}$ in the matrix ${Y}$, which is related to 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 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 $\mbox{{\boldmath${y}$}}_{k}$ is equal to the degree of $A_{k}$. Assuming that the degree of $A_{k}$ is $v$ and $\{e_{i_{1}},e_{i_{2}},\cdots,e_{i_{v}}\}$ have the common vertex $A_{k}$. If the direction of $e_{i_{s}}$ points to $A_{k}$, then $y_{k,i_{s}}=1$, otherwise $y_{k,i_{s}}=-1$. Consider the $k$-th column vector $\mbox{{\boldmath${d}$}}_{k}$ in the matrix ${Y}$, which is related to the edge $e_{k}$. If $A_{i}$ is not a vertex of the edge $e_{k}$, then $y_{ik}=0$. Since each edge in a tetrahedral mesh only has two vertices, the column vector $\mbox{{\boldmath${d}$}}_{k}$ only has two nonzero entries. Assuming that the initial and terminal points of the edge $e_{k}$ are $A_{i_{1}}$ and $A_{i_{2}}$ respectively, then $y_{i_{1}k}=-1$ and $y_{i_{2}k}=1$, where $i_{1}<i_{2}$. Obviously, the sparse matrix ${Y}$ can be easily obtained by the mesh data in $\mathcal{T}_{h}$. The sparse matrix ${Y}$ is usually called the directional connectivity matrix associated with a tetrahedral mesh $\mathcal{T}_{h}$.

The sum of all entries in the column vector $\mbox{{\boldmath${d}$}}_{k}$ is equal to zero. This implies that the homogenous linear equation

has a special nonzero solution $\beta=[1,1,\cdots,1]^{T}\in{\mathbb{C}^{m}}$. An easy computation shows that $\textrm{rank}(\mbox{{\boldmath${Y}$}})=\textrm{rank}(\mbox{{\boldmath${Y}$}}^{{\dagger}})=m-1$. Therefore all solutions of (7) form a linear space of one dimension, that is

where Null stands for taking the nullspace of a matrix.

The matrix ${Y}$ builds a relation among ${A}$, ${M}$ and ${C}$. It can be proved that the following matrix identities are valid:

where ${O}$ is a null $m\times n$ matrix. In fact, for $\forall~{}1\leq i\leq m$ and $\forall~{}1\leq l\leq n$, we have

where $(\mbox{{\boldmath${S}$}})_{il}$ is the entry at $i$-th row and $l$-th column of the matrix ${S}$. 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 ${C}$ by using the sparse matrices ${Y}$ and ${M}$, instead of the calculation by using the continuous sesquilinear forms $\mathcal{C}$ directly.

If the eigenvalue $\Lambda_{h}$ is nonzero in (5a), then (5b) can be deduced from (5a). As a matter of fact, $\mbox{{\boldmath${M}$}}\xi=\Lambda_{h}^{-1}\mbox{{\boldmath${A}$}}\xi$ can be derived from (5a), and then we get $\mbox{{\boldmath${C}$}}\xi=\mbox{{\boldmath${Y}$}}\mbox{{\boldmath${M}$}}\xi=\Lambda_{h}^{-1}\mbox{{\boldmath${Y}$}}\mbox{{\boldmath${A}$}}\xi={\bf{0}}$ by (9), which is just (5b). It is worthwhile to point out the number of zero eigenvalues in (5a) is equal to $\dim(\nabla{S^{h}})$, and these zero eigenvalues are all spurious. The dimension of the linear space $\nabla{S^{h}}$ is equal to $m-1$, which shows that the larger the number of mesh nodes $m$ 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.

Consider the following generalized eigenvalue problem:

where $\alpha$ is a penalty constant, which is required the user to set and $\|\varphi\|_{2}$ is the Euclidean norm of a given vector $\varphi$. 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 $\alpha$ is usually set to a large positive real number. The reason is as follows: Obviously, from (10), it is easy to deduce that

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

Multiplying both sides of (12) on the left with ${\xi}^{{\dagger}}$, then we can get $\|\mbox{{\boldmath${C}$}}\xi\|_{2}=0$. As a result, the homogeneous linear equation (5b) is obtained again. Substituting (5b) into (10) gives $\mbox{{\boldmath${A}$}}\xi=\Lambda_{h}^{\prime}\mbox{{\boldmath${M}$}}\xi$. However, if $\alpha$ takes sufficiently large, then this will lead to $\mbox{{\boldmath${A}$}}+\alpha\mbox{{\boldmath${C}$}}^{{\dagger}}\mbox{{\boldmath${C}$}}\approx\alpha\mbox{{\boldmath${C}$}}^{{\dagger}}\mbox{{\boldmath${C}$}}$ because of the finite word length in a computer. The result is that the information of matrix ${A}$ 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 $\alpha$ in (10). It is worthwhile that choosing an appropriate penalty parameter $\alpha$ is important to the penalty method. Webb [17] has studied this problem, and support a method to select this suitable parameter $\alpha$.

It can be seen that the eigenpair $(\Lambda_{h},\xi)$ of (5) is always the eigenpair $(\Lambda_{h},\xi)$ of (10). However, the eigenpair $(\Lambda_{h}^{\prime},\xi)$ 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. 1.

   Assuming that $(\Lambda_{h}^{\prime},\xi)$ with $\|\xi\|_{2}=1$ is an eigenpair of (10). The quantity $\|\mbox{{\boldmath${C}$}}\xi\|_{2}$ can be used to identify the eigenvalues with physical significance. If $\|\mbox{{\boldmath${C}$}}\xi\|_{2}$ is very small, then $\Lambda_{h}^{\prime}$ is an eigenvalue corresponding to a physical mode. Otherwise, $\Lambda_{h}^{\prime}$ is an eigenvalue associated with the spurious mode.

2. 2.

   The eigenvalues corresponding to physical eigenmodes will not change as the penalty constant $\alpha$ changes, but the eigenvalues corresponding to spurious modes will change as the penalty constant $\alpha$ changes. Therefore, under the same mesh, set $\alpha=\alpha_{1},\alpha_{2},\alpha_{3},\cdots$ 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.

Consider the generalized eigenvalue problem:

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

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

Assuming that $(\Lambda_{h}^{\prime},[\xi;\zeta])$ is the eigenpair of (14), where $[\xi;\zeta]$ is the corresponding eigenvector associated with $\Lambda_{h}^{\prime}$ in (14). Multiplying both sides of (13a) on the left with ${Y}$, we obtain

Substituting (9) and into (16) gives

If (15) can be derived from (17), then the eigenpair $(\Lambda_{h}^{\prime},[\xi;\zeta])$ of (14) is also an eigenpair $(\Lambda_{h}^{\prime},\xi)$ of (5). Conversely, assuming that $(\Lambda_{h},\xi)$ is the eigenpair of (5). If we take a vector $\zeta$ such that (15) is valid (Obviously, this is achievable), then the eigenpair $(\Lambda_{h},\xi)$ of (5) is an eigenpair $(\Lambda_{h},[\xi;\zeta])$ 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

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

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

We now prove that if $\overline{\overline{\mu}}^{{\dagger}}=\overline{\overline{\mu}}>0$, then (20) holds. In fact, it follows that $\mbox{{\boldmath${M}$}}^{{\dagger}}=\mbox{{\boldmath${M}$}}>0$ from $\overline{\overline{\mu}}^{{\dagger}}=\overline{\overline{\mu}}>0$. Hence, there exists an invertible square matrix ${X}$ such that $\mbox{{\boldmath${M}$}}=\mbox{{\boldmath${M}$}}^{{\dagger}}=\mbox{{\boldmath${X}$}}^{{\dagger}}\mbox{{\boldmath${X}$}}$. In terms of the fundamental theory in linear algebra, it is known that

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 $\overline{\overline{\mu}}^{{\dagger}}=\overline{\overline{\mu}}>0$.

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 $\zeta$ in (14) is the same. Next, we briefly prove this conclusion. In fact, since the matrix ${M}$ is invertible, $\mbox{{\boldmath${Y}$}}^{{\dagger}}\zeta=\bf{0}$ can be derived from $\mbox{{\boldmath${C}$}}^{{\dagger}}\zeta=\bf{0}$ and $\mbox{{\boldmath${C}$}}=\mbox{{\boldmath${Y}$}}\mbox{{\boldmath${M}$}}$. By means of (8), $\zeta=c\beta$ is valid, which shows that each entry of eigenvector $\zeta$ 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.

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 $\mathcal{V}$ in $\mathbb{C}^{n}$. Set $\mathcal{V}=\textrm{span}\{q_{1},q_{2},\cdots,q_{r}\}$, where $r=\dim{\mathcal{V}}$ and $q_{i}\in{\mathbb{C}^{n}}$. If the matrix ${M}$ is invertible, then $\textrm{rank}(\mbox{{\boldmath${C}$}})=\textrm{rank}(\mbox{{\boldmath${Y}$}})=m-1$ is valid, which implies $r=\dim{\mathcal{V}}=n-m+1$. Set $\mbox{{\boldmath${Q}$}}=[q_{1},q_{2},\cdots,q_{r}]\in{\mathbb{C}^{n\times r}}$.

The basic idea of solving (5) is that choosing $\Lambda_{h}\in{\mathbb{C}}$ and $\xi\in{\mathcal{V}}$ such that

This is called the Galerkin condition [18]. Set $\xi=\mbox{{\boldmath${Q}$}}y$, where $y\in{\mathbb{C}^{r}}$. The Galerkin condition (22) can be equivalently expressed the following equation

In order to compute the eigenpair of (23), the matrix ${Q}$ 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 ${Q}$ for the null space of the above-mentioned large sparse matrix ${C}$. However, since Coleman and Pothen [19] prove that finding the sparsest basis for the null space of an underdetermined matrix is NP-complete hard, we cannot seek the sparsest matrix ${Q}$ for the null space of ${C}$.

Based on the importance of numerical stability, a set of normalized orthogonal bases $\{q_{1},q_{2},\cdots,q_{r}\}$ is used in this numerical calculation. In such a case, $\mbox{{\boldmath${Q}$}}^{{\dagger}}\mbox{{\boldmath${Q}$}}=\mbox{{\boldmath${I}$}}_{r}$ holds, where $\mbox{{\boldmath${I}$}}_{r}$ is the identity matrix of order $r$. Here, we employ singular value decomposition (SVD) technique to seek ${Q}$. Suppose that

is the SVD of the matrix ${C}$. We take $\mbox{{\boldmath${Q}$}}=\mbox{{\boldmath${V}$}}(:,n-r+1:n)$, where $\mbox{{\boldmath${V}$}}(:,n-r+1:n)$ is a submatrix consisting of the last $r$ columns of ${V}$, then $\mbox{{\boldmath${C}$}}\mbox{{\boldmath${Q}$}}=\mbox{{\boldmath${O}$}}$ is valid.

If the eigenpair $(\Lambda_{h},y)$ of (23) is obtained by using the implicitly restarted Arnoldi methods [15], then $(\Lambda_{h},\mbox{{\boldmath${Q}$}}y)$ 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 $\alpha$ is not an easy work.

The main advantage of augmented method can preserve the sparsity of the matrix. Moreover, when $\overline{\overline{\mu}}_{r}^{{\dagger}}=\overline{\overline{\mu}}_{r}>0$, 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 ${Q}$ is usually dense.

In a word, these three methods have their own advantages and disadvantages.

Let us consider an empty spherical resonator with the radius $r=1$ m. The exact eigenvalue associated with the dominant mode is $\Lambda=7.52793\,\textrm{m}^{-2}$ [20]. Furthermore, the algebraic multiplicity of the exact eigenvalue $\Lambda$ is 3. Suppose that the numerical eigenvalues $\Lambda_{h}^{(1)}$, $\Lambda_{h}^{(2)}$ and $\Lambda_{h}^{(3)}$ are the approximation of the exact eigenvalue $\Lambda$. Set $\Lambda_{h}=(\Lambda_{h}^{(1)}+\Lambda_{h}^{(2)}+\Lambda_{h}^{(3)})/3$. We employ the penalty method, augmented method and projection method to solve this spherical resonant cavity problem, and then list the numerical eigenvalues $\Lambda_{h}(\textrm{pe},\alpha)$, $\Lambda_{h}(\textrm{au})$ and $\Lambda_{h}(\textrm{pr})$ in Table I. In order to compare with the efficiency of these three numerical methods, the CPU time is also given in Table I.

In these three numerical methods, the time and memory consumed by the projection method is the largest since the dense matrix ${Q}$ obtained by SVD is used. This shows that the projection method is not efficient. In addition, the CPU time and memory consumed by the penalty method and the augmented method are roughly equivalent.

In Table I, we can also see that $\Lambda_{h}(\textrm{pe},800)=\Lambda_{h}(\textrm{au})\approx\Lambda_{h}(\textrm{pr})\approx\Lambda$ under the finest mesh. This shows that our numerical implementations are correct. In this example, we do find that there are many eigenvalues associated with the spurious modes. These numerical eigenvalues are less than $\Lambda$ provided that the penalty parameter $\alpha$ is less than $700$. However, the numerical eigenvalues obtained by the augmented method and projection method are all physical.

In this subsection we consider a cylindrical cavity with the radius $r=0.2$ m and the height $h=0.5$ m. Assuming that the relative permittivity and permeability tensor of the medium in the whole cylindrical cavity are

The first three numerical eigenvalues $\Lambda_{h}^{i}(\textrm{au})$ and $\Lambda_{h}^{i}(\textrm{pr})$ ($i=1,2,3$) are shown in Table II. In Table II, it can be observed that $\Lambda_{h}^{i}(\textrm{au})\approx\Lambda_{h}^{i}(\textrm{pr})$, $i=1,2,3$, and they coincide with the eigenvalues corresponding to physical modes from COMSOL Multiphysics 5.2a. In the COMSOL simulation, the eigenvalues associated with physical modes are obtained by the fourth mesh ($h=0.0428$ m). Notice that there are many spurious zero modes in the numerical results of COMSOL Multiphysics 5.2a. However, there are no any spurious modes in the numerical results of the augmented method and projection method.

Under the second mesh ($h=0.0714$ m), we employ the penalty method to solve this cylindrical cavity problem, where $\alpha=800$ is taken in the numerical calculation, and then list the first forty numerical eigenvalues in Fig. 2. In Fig. 2, one can see that there are only two eigenvalues with physical significance, and the rest are all the eigenvalues without physical significance, whose imaginary part is zero. Furthermore, these two eigenvalues with physical significance obtained by the penalty method are equal to the ones obtained by augmented method. In addition, in (14), we do find that the eigenvector $\zeta$ associated with each eigenvalue is a vector consisting of the same entry.

In this subsection we try to find the eigenmode of resonant cavity filled with an electric and magnetic lossy medium. The penalty method, augmented method and projection method are used to solve this problem, and then we list the numerical eigenvalues associated with the first two physical modes in Table III.

Suppose that the geometric shape of the cavity in this example is the same as the one in the example B. In the cylindrical cavity, the relative permittivity and permeability tensor of the medium are

Obviously, the above material is both electric and magnetic lossy. Since the exact solution to this problem is unknown, we employ COMSOL Multiphysics 5.2a to simulate this problem, and then obtain the approximate eigenvalues of certain accuracy. The eigenvalues with physical significance from COMSOL are obtained by the fourth mesh ($h=0.0428$ m). Notice that many spurious zero modes appear in the numerical results of COMSOL.

The numerical eigenvalues from the penalty method and projection method are listed in Table III. In Table III, we can see that $\Lambda_{h}^{i}(\textrm{pe},1000)\approx\Lambda_{h}^{i}(\textrm{au})\approx\Lambda_{h}^{i}(\textrm{pr})$, $i=1,2$, which coincide with the eigenvalues corresponding to physical significance from COMSOL. Here it is worthwhile to emphasize that the projection method can remove all the spurious modes.