## Accelerated solution of the frequency-domain Maxwell’s equations by engineering the eigenvalue distribution of the operator |

Optics Express, Vol. 21, Issue 19, pp. 22578-22595 (2013)

http://dx.doi.org/10.1364/OE.21.022578

Acrobat PDF (1565 KB)

### Abstract

We introduce a simple method to accelerate the convergence of iterative solvers of the frequency-domain Maxwell’s equations for deep-subwavelength structures. Using the continuity equation, the method eliminates the high multiplicity of near-zero eigenvalues of the operator while leaving the operator nearly positive-definite. The impact of the modified eigenvalue distribution on the accelerated convergence is explained by visualizing residual vectors and residual polynomials.

© 2013 OSA

## 1. Introduction

*e*

^{+iωt}and nonmagnetic materials, Maxwell’s equations reduce to where

*ε*is the electric permittivity (which can be complex);

*μ*

_{0}is the magnetic permeability of vacuum;

*ω*is the angular frequency;

**E**and

**J**are the electric field and electric current source density, respectively.

1. N. J. Champagne II, J. G. Berryman, and H. M. Buettner, “FDFD: A 3D finite-difference frequency-domain code for electromagnetic induction tomography,” J. Comput. Phys. **170**, 830–848 (2001). [CrossRef]

3. U. S. Inan and R. A. Marshall, *Numerical Electromagnetics: The FDTD Method* (Cambridge University, 2011). Ch. 14. [CrossRef]

*A*is a matrix representing the operator [∇× (∇× ) −

*ω*

^{2}

*μ*

_{0}

*ε*];

*x*is an unknown column vector representing

**E**;

*b*is a column vector representing −

*iωμ*

_{0}

**J**. The matrix

*A*thus constructed is sparse (with only 13 nonzero elements per row when generated by the FDFD method) and typically very large (often with more than 10 million rows and columns for three-dimensional (3D) problems). To solve a system with such a large and sparse matrix, iterative methods are usually preferred to direct methods [5

5. V. Simoncini and D. B. Szyld, “Recent computational developments in Krylov subspace methods for linear systems,” Numer. Linear Algebra Appl. **14**, 1–59 (2007). [CrossRef]

6. W. Cai, W. Shin, S. Fan, and M. L. Brongersma, “Elements for plasmonic nanocircuits with three-dimensional slot waveguides,” Adv. Mater. **22**, 5120–5124 (2010). [CrossRef] [PubMed]

8. L. Verslegers, P. Catrysse, Z. Yu, W. Shin, Z. Ruan, and S. Fan, “Phase front design with metallic pillar arrays,” Opt. Lett. **35**, 844–846 (2010). [CrossRef] [PubMed]

9. J. T. Smith, “Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator,” Geophys. **61**, 1319–1324 (1996). [CrossRef]

11. C. J. Weiss and G. A. Newman, “Electromagnetic induction in a generalized 3D anisotropic earth, Part 2: The LIN preconditioner,” Geophys. **68**, 922–930 (2003). [CrossRef]

10. G. A. Newman and D. L. Alumbaugh, “Three-dimensional induction logging problems, Part 2: A finite-difference solution,” Geophys. **67**, 484–491 (2002). [CrossRef]

11. C. J. Weiss and G. A. Newman, “Electromagnetic induction in a generalized 3D anisotropic earth, Part 2: The LIN preconditioner,” Geophys. **68**, 922–930 (2003). [CrossRef]

*E*-field as

**E**= Ψ + ∇

*φ*, where Ψ is a divergence-free vector field and

*φ*is a scalar field [9

9. J. T. Smith, “Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator,” Geophys. **61**, 1319–1324 (1996). [CrossRef]

15. R. Hiptmair, F. Kramer, and J. Ostrowski, “A robust Maxwell formulation for all frequencies,” IEEE Trans. Magn. **44**, 682–685 (2008). [CrossRef]

*·*Ψ = 0, Eq. (1) is written as where the operator ∇× (∇× ), which has a huge null space, is replaced with the negative Laplacian −∇

^{2}, which is positive-definite for appropriate boundary conditions and thus has the smallest possible null space. However, these techniques either solve an extra equation for the extra unknown

*φ*at every iteration step [9

9. J. T. Smith, “Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator,” Geophys. **61**, 1319–1324 (1996). [CrossRef]

12. V. L. Druskin, L. A. Knizhnerman, and P. Lee, “New spectral Lanczos decomposition method for induction modeling in arbitrary 3-D geometry,” Geophys. **64**, 701–706 (1999). [CrossRef]

13. E. Haber, U. M. Ascher, D. A. Aruliah, and D. W. Oldenburg, “Fast simulation of 3D electromagnetic problems using potentials,” J. Comput. Phys. **163**, 150–171 (2000). [CrossRef]

15. R. Hiptmair, F. Kramer, and J. Ostrowski, “A robust Maxwell formulation for all frequencies,” IEEE Trans. Magn. **44**, 682–685 (2008). [CrossRef]

**J**= 0) position, where Eq. (1) can be modified to for an arbitrary constant

*s*; note that the right-hand side is 0 because

**J**= 0. In this class of techniques, Eqs. (1) and (5) are solved at positions with and without sources, respectively.

16. K. Beilenhoff, W. Heinrich, and H. Hartnagel, “Improved finite-difference formulation in frequency domain for three-dimensional scattering problems,” IEEE Trans. Microwave Theory Tech. **40**, 540–546 (1992). [CrossRef]

*s*= +1 to boundary value problems described in [17

17. A. Christ and H. L. Hartnagel, “Three-dimensional finite-difference method for the analysis of microwave-device embedding,” IEEE Trans. Microwave Theory Tech. **35**, 688–696 (1987). [CrossRef]

**J**= 0 everywhere, so Eq. (5) was solved throughout the entire simulation domain.

16. K. Beilenhoff, W. Heinrich, and H. Hartnagel, “Improved finite-difference formulation in frequency domain for three-dimensional scattering problems,” IEEE Trans. Microwave Theory Tech. **40**, 540–546 (1992). [CrossRef]

*s*. It also did not report whether its technique leads to accelerated convergence for problems with sources, even though many problems have nonzero electric current sources

**J**inside the simulation domain. Reference [1

1. N. J. Champagne II, J. G. Berryman, and H. M. Buettner, “FDFD: A 3D finite-difference frequency-domain code for electromagnetic induction tomography,” J. Comput. Phys. **170**, 830–848 (2001). [CrossRef]

*s*= +1 to problems with sources, but only in order to suppress spurious modes rather than to accelerate convergence.

**J**exist inside the simulation domain [18]. Unlike the previous technique that made the modification only at source-free positions, our technique modifies Eq. (1) everywhere including positions with sources. For the modification, we utilize the continuity equation which can be derived by taking the divergence of Eq. (1). When Eq. (6) is manipulated appropriately and then added to Eq. (1), we obtain for a constant

*s*. The modified equation (7) is the equation to solve in this paper.

*E*-field of Eq. (7) is the same as the solution of the original equation (1) regardless of the value of

*s*, because the solution of Eq. (1) always satisfies Eq. (6). However, the choice of

*s*affects the convergence speed of iterative methods significantly. In this paper, we demonstrate that

*s*= −1 induces faster convergence speed than other values of

*s*by comparing the convergence behavior of iterative methods for

*s*= −1, 0, +1; the latter two values of

*s*are of particular interest, because

*s*= 0 reduces Eq. (7) to the original equation (1) and

*s*= +1 is the value that Ref. [16

16. K. Beilenhoff, W. Heinrich, and H. Hartnagel, “Improved finite-difference formulation in frequency domain for three-dimensional scattering problems,” IEEE Trans. Microwave Theory Tech. **40**, 540–546 (1992). [CrossRef]

*s*. There are many general mathematical studies about the dependence of the convergence behavior on the eigenvalue distribution [19

19. A. Jennings, “Influence of the eigenvalue spectrum on the convergence rate of the conjugate gradient method,” IMA J. Appl. Math. **20**, 61–72 (1977). [CrossRef]

26. O. Axelsson, “Iteration number for the conjugate gradient method,” Math. Comput. Simulat. **61**, 421–435 (2003). [CrossRef]

5. V. Simoncini and D. B. Szyld, “Recent computational developments in Krylov subspace methods for linear systems,” Numer. Linear Algebra Appl. **14**, 1–59 (2007). [CrossRef]

*s*= 0 because the operator has eigenvalues clustered near zero, and for

*s*= +1 because the operator is strongly indefinite.

*s*= 0, −1, +1 for a simple homogeneous system. We also define the low-frequency regime rigorously in the section. In Sec. 3, we relate the eigenvalue distribution with the convergence behavior of an iterative method. In Sec. 4, we solve Eq. (7) for a wide range of realistic 3D problems to compare the convergence behavior of an iterative method for the three values of

*s*, and we conclude in Sec. 5.

## 2. Eigenvalue distribution of the operator for a homogeneous system

*s*. The impact of

*s*on the eigenvalue distribution has been studied in detail in the literature of the deflation method (also known as the penalty method) [27

27. J. P. Webb, “The finite-element method for finding modes of dielectric-loaded cavities,” IEEE Trans. Microwave Theory Tech. **33**, 635–639 (1985). [CrossRef]

29. D. A. White and J. M. Koning, “Computing solenoidal eigenmodes of the vector Helmholtz equation: a novel approach,” IEEE Trans. Magn. **38**, 3420–3425 (2002). [CrossRef]

*ε*is constant, Eq. (7) is simplified to where the operator is Hermitian for real

*ε*. Because

*ε*is constant in this section, the eigenvalue distribution of

*T*is shifted from the eigenvalue distribution of a Hermitian operator by a constant −

*ω*

^{2}

*μ*

_{0}

*ε*. In the low-frequency regime such shift is negligible, and thus the eigen-value distribution of

*T*

_{0}approximates that of

*T*very well. Hence, we examine the eigenvalue distribution of

*T*

_{0}below to investigate the eigenvalue distribution of

*T*.

**F**

_{k}*e*

^{−ik·r}with are the three eigenfunctions of both ∇× (∇× ) and ∇(∇· ) for each wavevector

**k**. We also show in the same appendix that the corresponding three eigenvalues are for ∇×(∇× ), and for ∇(∇· ). Therefore,

*T*

_{0}has as three eigenvalues for each wavevector

**k**.

*T*

_{0}is greatly affected by the value of

*s*. Specifically, the multiplicity of the eigenvalue 0 depends critically on whether

*s*is 0 or not: for

*s*= 0

*T*

_{0}has a very high multiplicity of the eigenvalue 0 because Eq. (14) has 0 as an eigenvalue for every

**k**, whereas for

*s*≠ 0

*T*

_{0}does not have such a high multiplicity of the eigenvalue 0. The definiteness of

*T*

_{0}also depends on the value of

*s*: for

*s*≤ 0

*T*

_{0}is positive-semidefinite because the three eigenvalues in Eq. (14) are always nonnegative, whereas for

*s*> 0

*T*

_{0}is indefinite because Eq. (14) has both positive and negative numbers as eigenvalues. The different properties of the eigenvalue distributions of

*T*

_{0}for

*s*= 0,

*s*< 0, and

*s*> 0 are summarized in Table 1.

*T*

_{0}should approximately hold for the eigenvalue distributions of

*T*as well in the low-frequency regime, as mentioned in the discussion of Eqs. (9) and (10). Moreover, even though

*T*is a differential operator defined in an infinite space, it turns out that the description also applies to the matrix

*A*discretized from

*T*that is defined in a spatially bounded simulation domain.

*A*for a two-dimensional (2D) system shown in Fig. 1, a square domain filled with vacuum. The domain is discretized on a finite-difference grid with

*N*×

_{x}*N*= 50 × 50 cells and cell size Δ = 2 nm. Therefore, the matrix

_{y}*A*for each

*s*has 3

*N*= 7500 rows and columns, where the extra factor 3 accounts for the three Cartesian components of the

_{x}N_{y}*E*-field. We choose

*ω*corresponding to the vacuum wavelength

*λ*

_{0}= 1550nm, which puts the system in the low-frequency regime as will be seen at the end of this section. The matrices

*A*are constructed for three values of

*s*: 0, −1, and +1, each of which represents each category of

*s*in Table 1.

*A*for

*s*= 0, −1, +1 are shown as three plots in Fig. 2. In each plot, the horizontal axis represents eigenvalues, and it is divided into 41 intervals

*t*

_{−20},...,

*t*

_{0},...,

*t*

_{20}where

*t*

_{0}∋ 0. The height of the column on each interval corresponds to the number of the eigenvalues in the interval.

*A*shown in Fig. 2 agree well with the description of the eigenvalues of

*T*

_{0}in Table 1: the very tall column on

*t*

_{0}in Fig. 2(a) indicates the very high multiplicity of

*λ*≃ 0 for

*s*= 0, and the eigenvalues distributed over

*t*

_{j<0}and

*t*

_{j>0}in Fig. 2(c) indicate a strongly indefinite operator for

*s*= +1. In addition, the height of the column on

*t*

_{0}in Fig. 2(a) is about 2500, or one third of the total number of eigenvalues, which agrees with Eq. (14) for

*s*= 0 where one of the three eigenvalues is 0 for each

**k**; the columns on

*t*

_{j>0}are about 1.5 times taller in Fig. 2(b) than in Fig. 2(a), which also agrees with Eq. (14) where the number of |

**k**|

^{2}increases from two for

*s*= 0 to three for

*s*= −1.

*A*

_{0}is the matrix discretized from

*T*

_{0}of Eq. (10). For

*s*= 0, the eigenvalues of

*A*

_{0}range from 0 to

_{min}is the minimum grid cell size [31]; note that the range agrees with Fig. 2(a). The eigenvalue distribution of

*A*is the shifted eigenvalue distribution of

*A*

_{0}by −

*ω*

^{2}

*μ*

_{0}

*ε*. The low-frequency regime is where the magnitude of the shift is so small that

*A*has an almost identical eigenvalue distribution as

*A*

_{0}. Therefore, the condition for the low-frequency regime is Equation (15) is consistent with the condition introduced in [10

10. G. A. Newman and D. L. Alumbaugh, “Three-dimensional induction logging problems, Part 2: A finite-difference solution,” Geophys. **67**, 484–491 (2002). [CrossRef]

*A*

_{0}. We can rewrite Eq. (15) in terms of the vacuum wavelength

*λ*

_{0}as where

*ε*=

_{r}*ε/ε*

_{0}is the relative electric permittivity. The system described in Fig. 1 satisfies Eq. (16), so it is in the low-frequency regime.

## 3. Impact of the eigenvalue distribution on the convergence behavior of GMRES

*s*examined in Sec. 2 influence the convergence behavior of an iterative method to solve Eq. (8).

*s*= −1, 0, +1, we discretize Eq. (8) using the FDFD method for the system illustrated in Fig. 1 with an

*x*-polarized electric dipole current source placed at the center of the simulation domain. We then solve the discretized equation by an iterative method to observe the convergence behavior. The iterative method to use in this section is the general minimal residual (GMRES) method [33

33. Y. Saad and M. H. Schultz, “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM J. Sci. Stat. Comp. **7**, 856–869 (1986). [CrossRef]

5. V. Simoncini and D. B. Szyld, “Recent computational developments in Krylov subspace methods for linear systems,” Numer. Linear Algebra Appl. **14**, 1–59 (2007). [CrossRef]

*x*of

_{m}*Ax*=

*b*at the

*m*th iteration step. As

*m*increases,

*x*eventually converges to the exact solution. We assume that convergence is achieved when the residual vector satisfies ||

_{m}*r*||/||

_{m}*b*|| <

*τ*, where ||·|| is the 2-norm and

*τ*is a user-defined small positive number. In practice,

*τ*= 10

^{−6}is sufficiently small for accurate solutions. We use

*x*

_{0}= 0 as an initial guess solution throughout the paper.

*r*||/||

_{m}*b*|| versus the number

*m*of iteration steps for the three values of

*s*. As can be seen in the figure, the convergence behavior of GMRES is quite different for different

*s*, with

*s*= −1 far more superior than the other two choices of

*s*.

*s*= 0, and according to [21

21. S. L. Campbell, I. C. F. Ipsen, C. T. Kelley, and C. D. Meyer, “GMRES and the minimal polynomial,” BIT Numer. Math. **36**, 664–675 (1996). [CrossRef]

*A*for

*s*= 0 (see Fig. 2(a)). Also, the convergence is very slow for

*s*= +1, and Ref. [23

23. M. Benzi, G. H. Golub, and J. Liesen, “Numerical solution of saddle point problems,” Acta Numer. **14**, 1–137 (2005). Sec. 9.2. [CrossRef]

*A*for

*s*= +1 (see Fig. 2(c)) than for definite matrices. In this section we provide a more intuitive explanation for the convergence behavior by using the residual polynomial.

*is the set of all polynomials*

_{m}*p̃*of degree at most

_{m}*m*such that For each

*p̃*∈ 𝒫

_{m}*, we can define a column vector At the*

_{m}*m*th iteration step of GMRES, the residual vector

*r*of Eq. (17) is the

_{m}*r̃*with the smallest 2-norm [5

_{m}**14**, 1–59 (2007). [CrossRef]

*p̃*for

_{m}*r̃*=

_{m}*r*as the residual polynomial

_{m}*p*. Therefore, from Eq. (19) we have

_{m}*A*influences

*p*at each iteration step and hence influences the convergence behavior of GMRES. The matrix

_{m}*A*∈ ℂ

^{n×n}for our homogeneous system described in Fig. 1 is Hermitian because it is discretized from the Hermitian operator

*T*of Eq. (8). Hence, the eigendecomposition of

*A*is where with real eigenvalues

*λ*and the corresponding normalized eigenvectors

_{i}*v*, and

_{i}*V*

^{†}is the conjugate transpose of

*V*; note that

*V*is unitary, i.e.,

*V*

^{†}

*V*=

*I*. Substituting Eq. (21) in Eq. (19), we obtain because (

*V*Λ

*V*

^{†})

*=*

^{k}*V*Λ

^{k}V^{†}for any nonnegative integer

*k*. We then define a column vector whose

*i*th element, which is referred to as

*z̃*below, is the projection of

_{mi}*r̃*/||

_{m}*b*|| onto the direction of the

*i*th eigenvector

*v*. Similarly, we also define From Eq. (25) for

_{i}*m*= 0 and Eqs. (23) and (24), we obtain which can be written element-by-element as Because ||

*z̃*|| = ||

_{m}*r̃*||/||

_{m}*b*||, GMRES minimizes ||

*z̃*|| to ||

_{m}*z*|| when it minimizes ||

_{m}*r̃*|| to ||

_{m}*r*|| at the

_{m}*m*th iteration step.

*z̃*| is minimized to 0 when

_{mi}*p̃*has

_{m}*λ*as a root. Thus, the most ideal

_{i}*p̃*has all the

_{m}*n*eigenvalues of

*A*as its roots, because it reduces ||

*z̃*|| to 0. However,

_{m}*p̃*has at most

_{m}*m*roots, and

*m*, which is the number of iteration steps, is typically far less than

*n*. Therefore,

*p̃*needs to have its roots optimally placed near the eigenvalues to minimize ||

_{m}*z̃*||. Hence, the eigenvalue distribution of

_{m}*A*greatly influences the convergence behavior of GMRES.

*s*. We begin with

*s*= −1. In Fig. 4 we plot

*r*/||

_{m}*b*|| for

*s*= −1 as bar graphs at the first few iteration steps. The horizontal axis in each plot represents eigenvalues. We divide the range of eigenvalues into the same 41 intervals

*t*

_{−20},...,

*t*

_{0},...,

*t*

_{20}used in Fig. 2; note that

*t*

_{0}∋ 0. The height of the column on each interval is the norm of the projection of

*r*/||

_{m}*b*|| onto the space spanned by the eigenvectors whose corresponding eigenvalues are contained in the interval. More specifically, the height of the column on

*t*after

_{j}*m*iteration steps is Note that

*r*/||

_{m}*b*|| for

*s*= −1 shown in Fig. 4 are readily predicted from the corresponding eigenvalue distribution of the matrix

*A*presented in Fig. 2(b). For instance,

*A*has no eigenvalues in

*t*

_{j<0}, and therefore

*r*/||

_{m}*b*|| has components only in

*t*

_{j}_{≥0}throughout the iteration process as demonstrated in Fig. 4. Also,

*A*has very few eigenvalues in

*t*

_{0}, and thus

*r*

_{0}/||

*b*|| has a very weak component in

*t*

_{0}as can be seen in the

*m*= 0 plot in Fig. 4.

*r*/||

_{m}*b*|| with the residual polynomial to explain the convergence behavior of GMRES for

*s*= −1. The residual polynomial

*p*(

_{m}*λ*), which is obtained by solving a least squares problem, is also plotted in Fig. 4 at each iteration step. As the iteration proceeds, the residual polynomial in Fig. 4 has more and more roots, but only in

*t*

_{j}_{≥0}, because the eigenvalues exist only in

*t*

_{j}_{≥0}and the roots of residual polynomials should stay close to the eigenvalues as mentioned in the discussion following Eq. (27). Also, as Eq. (27) predicts, the columns in each plot of Fig. 4 almost vanish at the roots of the residual polynomial. Therefore, all the columns quickly shrink as the number of the roots of the residual polynomial increases in the iteration process of GMRES. The fast reduction of the column heights provides visualization of the fast convergence of GMRES for

*s*= −1 shown in Fig. 3.

*s*= 0. Figure 5 shows

*r*/||

_{m}*b*|| for

*s*= 0 at the first few iteration steps. Note that

*r*

_{0}/||

*b*|| has a tall column on

*t*

_{0}because

*A*has many eigenvalues in

*t*

_{0}as shown in Fig. 2(a). Also, the tall column on

*t*

_{0}persists during the initial period of the iteration process.

*s*= 0, we show that for a nearly positive-definite matrix the column on

*t*

_{0}is persistent during the initial period of the iteration process of GMRES in general. For that purpose, we compare the three polynomials

*p̃*∈ 𝒫

_{m}*shown in Fig. 6. The three*

_{m}*p̃*are chosen as candidates for the residual polynomial

_{m}*p*for a nearly positive-definite matrix, and therefore the roots of the polynomials are placed in

_{m}*t*

_{j≥0}according to the discussion following Eq. (27). The three

*p̃*have the same roots except for their smallest roots:

_{m}*p̃*in Fig. 6(a) does not have its smallest root in

_{m}*t*

_{0}, whereas

*p̃*in Figs. 6(b) and 6(c) do. Note that the latter two

_{m}*p̃*can shrink the column on

_{m}*t*

_{0}more effectively than the first

*p̃*according to Eq. (27).

_{m}*p̃*are steeper than the slopes at the corresponding roots of the first

_{m}*p̃*as shown in Fig. 6. In Appendix B, we prove rigorously that the slopes of

_{m}*p̃*at all roots indeed increase as the smallest root decreases in magnitude. In general,

_{m}*p̃*with steeper slopes at the roots oscillates with larger amplitudes around the horizontal axis because it varies faster around the axis; compare the amplitudes of oscillation in Fig. 6(a) with those in Figs. 6(b) and 6(c). The increased amplitudes of oscillation amplify |

_{m}*z̃*| overall according to Eq. (27), and thus ||

_{mi}*z̃*|| as well.

_{m}*t*

_{0}(by placing the smallest root of

*p̃*in

_{m}*t*

_{0}) is achieved only at the penalty of amplifying the columns on

*t*

_{j>0}. This penalty is too heavy when the columns on

*t*

_{j>0}constitute a considerable portion of ||

*z̃*||. Therefore, roots of residual polynomials are not placed in

_{m}*t*

_{0}until the columns on

*t*

_{j>0}become quite small, which results in the persistence of the column on

*t*

_{0}during the initial period of the iteration process.

*t*

_{0}remains almost the same at the initial iteration steps of GMRES,

*h*

_{00}of Eq. (28), which is the initial height of this column, provides an approximate lower bound of ||

*z*|| = ||

_{m}*r*||/||

_{m}*b*|| during the initial period of the iteration process. A more accurate lower bound is calculated as the norm of

*r*

_{0}/||

*b*|| projected onto the eigenspace of the eigenvalue closest to 0. For our example system, for

*s*= 0 the calculated lower bound is 0.707. Note that ||

*r*||/||

_{m}*b*|| for

*s*= 0 indeed stagnates initially at this value in Fig. 3. For

*s*= −1 the calculated lower bound is 4.16 × 10

^{−7}, at which ||

*r*||/||

_{m}*b*|| also stagnates as shown in Fig. 3. However, this value is much smaller than the lower bound for

*s*= 0, because for

*s*= −1 the initial height of the column on

*t*

_{0}is almost negligible as shown in the

*m*= 0 plot in Fig. 4. In fact, the value is smaller than the conventional tolerance

*τ*= 10

^{−6}mentioned below Eq. (17), so the stagnation does not deteriorate the convergence speed for

*s*= −1.

*s*= +1. Figure 7 shows

*r*/||

_{m}*b*|| for

*s*= +1 at some first (

*m*= 0, 4, 7, 11) and later (

*m*= 120, 140) iteration steps. Because the matrix

*A*for

*s*= +1 has both positive and negative eigenvalues as indicated in Fig. 2(c),

*r*/||

_{m}*b*|| has components in both

*t*

_{j>0}and

*t*

_{j<0}, but in the present example the components of

*r*/||

_{m}*b*|| are concentrated in

*t*

_{j<0}initially (

*m*= 0 plot in Fig. 7). Thus GMRES begins with the roots of residual polynomials placed in

*t*

_{j<0}(

*m*= 4 plot in Fig. 7). However, such residual polynomials have large values in

*t*

_{j>0}, so they amplify the initially very small components of

*r*/||

_{m}*b*|| in

*t*

_{j>0}according to Eq. (27), and eventually we reach a point where the components of

*r*/||

_{m}*b*|| in

*t*

_{j>0}and

*t*

_{j<0}become comparable (

*m*= 7 plot in Fig. 7). Afterwards, GMRES places the roots of residual polynomials in both

*t*

_{j>0}and

*t*

_{j<0}so that the components of

*r*/||

_{m}*b*|| in both regions are reduced.

*s*= +1 is initially quite similar to that for

*s*= −1 because

*r*

_{0}/||

*b*|| for

*s*= +1 has components concentrated in

*t*

_{j<0}and only a very weak component in

*t*

_{0}. Therefore, ||

*r*||/||

_{m}*b*|| reduces quickly for

*s*= +1 without stagnation during the initial period of the iteration process as shown in Fig. 3.

*r*||/||

_{m}*b*|| for

*s*= +1 slows down significantly, and eventually

*s*= +1 produces the slowest convergence among the three values of

*s*as shown in Fig. 3. The slow reduction of ||

*r*||/||

_{m}*b*|| is due to the very persistent column on

*t*

_{0}: comparing the plots for

*m*= 120 and

*m*= 140 in Fig. 7 shows that the column barely reduces for 20 iteration steps.

*t*

_{0}is quite persistent for a nearly positive-definite matrix. The argument relied on the properties proved in Appendix B about a polynomial

*p̃*∈ 𝒫

_{m}*with only positive roots. We can easily extend the proof in the appendix to*

_{m}*p̃*with both positive and negative roots, and then show that the column on

_{m}*t*

_{0}is persistent also for a strongly indefinite matrix, which explains the slow convergence for

*s*= +1 described above. However, the explanation is insufficient to explain why the convergence is

*much*slower for

*s*= +1 than for

*s*= 0 as indicated in Fig. 3.

*t*

_{0}is in fact even more persistent for a strongly indefinite matrix than for a nearly positive-definite matrix. For that purpose, we compare the two polynomials

*p̃*∈ 𝒫

_{m}*shown in Fig. 8. As can be seen from the locations of their roots, they are candidates for the residual polynomials for different matrices:*

_{m}*p̃*shown in Fig. 8(a) is appropriate for a nearly positive-definite matrix (referred to as

_{m}*A*

_{def}below), and

*p̃*shown in Fig. 8(b) is appropriate for a strongly indefinite matrix (referred to as

_{m}*A*

_{ind}below). Moreover, we choose these two

*p̃*to have the same smallest-magnitude root

_{m}*ζ*

_{0}in

*t*

_{0}. Being elements of 𝒫

*, both*

_{m}*p̃*satisfy Eq. (18). Hence, we have |

_{m}*p̃′*(

_{m}*ζ*

_{0})| ≃ 1/|

*ζ*

_{0}| for both

*p̃*, where

_{m}*p̃′*is the first derivative of

_{m}*p̃*.

_{m}*p̃′*| evaluated at a root of

_{m}*p̃*tends to decrease as the root gets closer to the median of the roots; see Appendix C for a more rigorous explanation. Hence, |

_{m}*p̃′*| ≤ 1/|

_{m}*ζ*

_{0}| tends to hold at most roots of

*p̃*for

_{m}*A*

_{def}, because

*ζ*

_{0}is one of the farthest roots from the median of the roots. On the other hand, |

*p̃′*| ≥ 1/|

_{m}*ζ*

_{0}| tends to hold at most roots of

*p̃*for

_{m}*A*

_{ind}, because

*ζ*

_{0}is one of the closest roots to the median of the roots. Therefore,

*p̃*for

_{m}*A*

_{ind}has much steeper slopes at most roots than

*p̃*for

_{m}*A*

_{def}in general, and thus has larger amplitudes of oscillation around the horizontal axis, as demonstrated in Fig. 8.

*t*

_{0}(by placing the smallest-magnitude root of

*p̃*in

_{m}*t*

_{0}) increases ||

*z*|| much more for a strongly indefinite matrix than for a nearly positive-definite matrix. Therefore, the column on

_{m}*t*

_{0}should be much more persistent for a strongly indefinite matrix than for a nearly positive-definite matrix in general, which explains the much slower convergence for

*s*= +1 than for

*s*= 0 in Fig. 3.

*s*= −1 produces the most superior convergence behavior;

*s*= 0 induces stagnation during the initial period of the iteration process due to the high multiplicity of eigenvalues near zero;

*s*= +1 leads to the slowest convergence overall due to the strongly indefinite matrix. We have provided a graphical explanation of the difference in the convergence behavior of GMRES, for which a strong theoretical basis exists, using a simple system of a homogeneous dielectric medium.

*s*= −1 is indeed optimal among all values in general. Compared with the case of

*s*= −1, according to Eq. (14), for

*s*> 0

*A*is always more strongly indefinite and therefore the convergence should be slower; for −1 <

*s*< 0

*A*has more eigenvalues clustered near zero and thus the initial stagnation period should be longer; for

*s*< −1

*A*has a wider eigenvalue value range, so the condition number of

*A*should be larger and the convergence should be slower [23

23. M. Benzi, G. H. Golub, and J. Liesen, “Numerical solution of saddle point problems,” Acta Numer. **14**, 1–137 (2005). Sec. 9.2. [CrossRef]

*s*= −1 is the value that leaves

*A*nearly positive-definite while removing the eigenvalues clustered near zero sufficiently without increasing the condition number. With separate numerical experiments we have verified that

*s*= −1 is indeed superior to values other than

*s*= 0 and

*s*= +1 as well.

*s*is in fact quite general in practical situations.

## 4. Convergence behavior of QMR for 3D inhomogeneous systems

*s*= −1 still induces faster convergence than

*s*= 0 and

*s*= +1. We note that the systems examined in this section are inhomogeneous and have complex

*ε*in general. The analyses in Secs. 2 and 3, therefore, do not hold strictly here. Nevertheless, we will see that the analyses for the homogeneous system in the previous sections provide insight in understanding the convergence behavior for more general systems examined in this section.

32. W. Shin and S. Fan, “Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers,” J. Comput. Phys. **231**, 3406–3431 (2012). [CrossRef]

*Ax*=

*b*corresponding to

*s*= −1, 0, +1 by the FDFD method. The number of the grid cells

*N*,

_{x}*N*, and

_{y}*N*and the grid cell sizes Δ

_{z}*, Δ*

_{x}*, and Δ*

_{y}*of the finite-difference grid used to discretize each system are shown in Table 2. Considering the parameters summarized in Table 3 and the condition (16), all the three systems are in the low-frequency regime.*

_{z}37. R. Freund and N. Nachtigal, “QMR: a quasi-minimal residual method for non-hermitian linear systems,” Numer. Math. **60**, 315–339 (1991). [CrossRef]

**14**, 1–59 (2007). [CrossRef]

*s*= −1 results in the fastest convergence, and the choice of

*s*= +1 barely leads to convergence. The three systems shown in Fig. 9 are chosen deliberately to include different materials such as dielectrics and metals and geometries with different degrees of complexity. Therefore, Fig. 9 suggests that the superiority of

*s*= −1 over both

*s*= 0 and

*s*= +1 is quite general.

*ε*is complex, because it is approximated well by the Hermitian operator (10).

38. R. W. Freund, G. H. Golub, and N. M. Nachtigal, “Iterative solution of linear systems,” Acta Numer. **1**, 57–100 (1992). Secs. 2.4 and 3.3. [CrossRef]

*s*= −1 is no longer very different from that for

*s*= 0 in this case.

## 5. Conclusion and final remarks

*μ*=

*μ*

_{0}). For magnetic materials (

*μ*≠

*μ*

_{0}), we note that a similar equation which can also be formulated from Maxwell’s equations and the continuity equation, can be used instead of Eq. (7) to accelerate the convergence of iterative methods. We leave the discussion on the optimal value of

*s*in this equation for future work.

## Appendix A: Eigenvalues and eigenfunctions of ∇×(∇× ) and ∇(∇· )

**k**-space representations of the operators, in this section we derive the eigenvalues Eq. (12) of ∇×(∇× ) and Eq. (13) of ∇(∇· ) as well as their corresponding eigenfunctions.

**r**represents position,

**k**=

**x̂**

*k*+

_{x}**ŷ**

*k*+

_{y}**ẑ**

*k*is a real constant wavevector, and

_{z}**k**-dependent vector.

**F**) =

*λ*

**F**and ∇(∇

**·F**) =

*λ*

**F**by substituting Eq. (30) for

**F**. Then, the eigenvalue equation for ∇×(∇× ) is and the eigenvalue equation for ∇(∇· ) is

## Appendix B: Effect of the smallest root of *p̃*_{m} ∈ 𝒫_{m} on the slopes at the roots

_{m}

_{m}

*p̃*∈ 𝒫

_{m}*with all positive roots become steeper when the smallest root decreases in magnitude. This behavior is illustrated in Fig. 6.*

_{m}*p̃*∈ 𝒫

_{m}*satisfies the condition (18), it can be factored as where*

_{m}*d*≤

_{m}*m*is the degree of

*p̃*and

_{m}*ζ*’s are the roots of

_{i}*p̃*. Hence, the slope of

_{m}*p̃*at a root

_{m}*ζ*is

_{k}*ζ*

_{1}< ··· <

*ζ*

_{dm}. We can easily show that |

*p̃′*(

_{m}*ζ*)| increases for any

_{k}*k*when

*ζ*

_{1}decreases toward zero (while remaining positive) as follows. For

*k*= 1, we have which clearly increases as

*ζ*

_{1}decreases to 0. For

*k*> 1, we have where the parentheses enclose the only quantity that is dependent on

*ζ*

_{1}. We can therefore see that |

*p̃′*(

_{m}*ζ*)| increases as

_{k}*ζ*

_{1}decreases for

*k*> 1 as well. Therefore, for a given

*p̃*∈ 𝒫

_{m}*whose roots are all positive, the slopes of*

_{m}*p̃*at the roots become steeper if the smallest root decreases in magnitude while remaining positive. This situation is illustrated by the transition from Fig. 6(a) to Fig. 6(b).

_{m}*ζ*

_{1}is replaced by a negative value, as long as the negative value is smaller in magnitude than the original

*ζ*

_{1}. Replacing the originally positive

*ζ*

_{1}with a negative quantity that is smaller in magnitude is equivalent to first replacing

*ζ*

_{1}with a smaller positive value and then flipping its sign. Because we have already shown above that the slopes get steeper when the originally positive

*ζ*

_{1}is replaced by a smaller positive value, it is sufficient to show that flipping the sign of

*ζ*

_{1}makes the slopes even steeper. For a negative

*ζ*

_{1}, the slopes at the roots are and for

*k*> 1. These slopes are steeper than Eqs. (38) and (39), respectively, which are the slopes for a positive

*ζ*

_{1}with the same magnitude. Therefore, for a given

*p̃*∈ 𝒫

_{m}*whose roots are all positive, the slopes of*

_{m}*p̃*at the roots become steeper if the smallest root is replaced by the one that is smaller in magnitude but negative. This situation is illustrated by the transition from Fig. 6(a) to Fig. 6(c).

_{m}## Appendix C: Trend in the slopes of a polynomial at the roots

*p*with all real roots, and show that the slope of

*p*evaluated at a root closer to the median of the roots tends to be gentler than the slope evaluated at a root farther away from the median of the roots. This behavior is illustrated in Fig. 8.

*m*, with

*ζ*

_{1}< ··· <

*ζ*. The slope of

_{m}*p*at a root

*ζ*is

_{k}*p′*(

*ζ*

_{k}_{+1})|/|

*p′*(

*ζ*)|. We first consider the case where the roots are evenly spaced, i.e.,

_{k}*ζ*

_{i}_{+1}−

*ζ*= (const.), for which we have Equation (44) is an increasing function of

_{i}*k*for 1 ≤

*k*≤

*m*−1, and it is less than 1 for

*k*<

*m*/2 and greater than 1 for

*k*>

*m*/2. Therefore, |

*p′*(

*ζ*)| is largest for

_{k}*k*= 1 and

*k*=

*m*, and it decreases as

*k*becomes closer to

*k*= ⌈(

*m*+ 1)/2⌉ and

*k*= ⌊(

*m*+ 1)/2⌋, which are the medians of the indices. In other words, for

*p*with evenly spaced roots, the slopes of

*p*get gentler at the roots closer to the median of the roots.

*p*with unevenly spaced roots, unless the unevenness is too severe. To verify the expectation, we examine |

*p′*(

*ζ*

_{k}_{+1})|/|

*p′*(

*ζ*)| without assuming

_{k}*ζ*

_{i}_{+1}−

*ζ*= (const.):

_{i}*k*. Hence, |

*p′*(

*ζ*

_{k}_{+1})|/|

*p′*(

*ζ*)| tends to be less than 1 for smaller

_{k}*k*, and it tends to be greater than 1 for larger

*k*. This means that as

*k*increases |

*p′*(

*ζ*)| tends to decrease first and then tends to increase. Therefore, even if the roots of

_{k}*p*are unevenly spaced, the slopes of

*p*tend to get gentler at the roots closer to the median of the roots.

## Acknowledgment

## References and links

1. | N. J. Champagne II, J. G. Berryman, and H. M. Buettner, “FDFD: A 3D finite-difference frequency-domain code for electromagnetic induction tomography,” J. Comput. Phys. |

2. | G. Veronis and S. Fan, “Overview of simulation techniques for plasmonic devices,” in “ |

3. | U. S. Inan and R. A. Marshall, |

4. | J.-M. Jin, |

5. | V. Simoncini and D. B. Szyld, “Recent computational developments in Krylov subspace methods for linear systems,” Numer. Linear Algebra Appl. |

6. | W. Cai, W. Shin, S. Fan, and M. L. Brongersma, “Elements for plasmonic nanocircuits with three-dimensional slot waveguides,” Adv. Mater. |

7. | G. Veronis and S. Fan, “Bends and splitters in metal-dielectric-metal subwavelength plasmonic waveguides,” Appl. Phys. Lett. |

8. | L. Verslegers, P. Catrysse, Z. Yu, W. Shin, Z. Ruan, and S. Fan, “Phase front design with metallic pillar arrays,” Opt. Lett. |

9. | J. T. Smith, “Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator,” Geophys. |

10. | G. A. Newman and D. L. Alumbaugh, “Three-dimensional induction logging problems, Part 2: A finite-difference solution,” Geophys. |

11. | C. J. Weiss and G. A. Newman, “Electromagnetic induction in a generalized 3D anisotropic earth, Part 2: The LIN preconditioner,” Geophys. |

12. | V. L. Druskin, L. A. Knizhnerman, and P. Lee, “New spectral Lanczos decomposition method for induction modeling in arbitrary 3-D geometry,” Geophys. |

13. | E. Haber, U. M. Ascher, D. A. Aruliah, and D. W. Oldenburg, “Fast simulation of 3D electromagnetic problems using potentials,” J. Comput. Phys. |

14. | J. Hou, R. Mallan, and C. Torres-Verdin, “Finite-difference simulation of borehole EM measurements in 3D anisotropic media using coupled scalar-vector potentials,” Geophys. |

15. | R. Hiptmair, F. Kramer, and J. Ostrowski, “A robust Maxwell formulation for all frequencies,” IEEE Trans. Magn. |

16. | K. Beilenhoff, W. Heinrich, and H. Hartnagel, “Improved finite-difference formulation in frequency domain for three-dimensional scattering problems,” IEEE Trans. Microwave Theory Tech. |

17. | A. Christ and H. L. Hartnagel, “Three-dimensional finite-difference method for the analysis of microwave-device embedding,” IEEE Trans. Microwave Theory Tech. |

18. | At the final stage of our work, we were made aware of a related work by M. Kordy, E. Cherkaev, and P. Wannamaker, “Schelkunoff potential for electromagnetic field: proof of existence and uniqueness” (to be published), where an equation similar to our Eq. (7) with s= −1 was developed. |

19. | A. Jennings, “Influence of the eigenvalue spectrum on the convergence rate of the conjugate gradient method,” IMA J. Appl. Math. |

20. | A. van der Sluis and H. van der Vorst, “The rate of convergence of conjugate gradients,” Numer. Math. |

21. | S. L. Campbell, I. C. F. Ipsen, C. T. Kelley, and C. D. Meyer, “GMRES and the minimal polynomial,” BIT Numer. Math. |

22. | S. Goossens and D. Roose, “Ritz and harmonic Ritz values and the convergence of FOM and GMRES,” Numer. Linear Algebra Appl. |

23. | M. Benzi, G. H. Golub, and J. Liesen, “Numerical solution of saddle point problems,” Acta Numer. |

24. | B. Beckermann and A. B. J. Kuijlaars, “Superlinear convergence of conjugate gradients,” SIAM J. Numer. Anal. |

25. | B. Beckermann and A. B. J. Kuijlaars, “On the sharpness of an asymptotic error estimate for conjugate gradients,” BIT Numer. Math. |

26. | O. Axelsson, “Iteration number for the conjugate gradient method,” Math. Comput. Simulat. |

27. | J. P. Webb, “The finite-element method for finding modes of dielectric-loaded cavities,” IEEE Trans. Microwave Theory Tech. |

28. | F. Kikuchi, “Mixed and penalty formulations for finite element analysis of an eigenvalue problem in electromagnetism,” Comput. Method Appl. Mech. Eng. |

29. | D. A. White and J. M. Koning, “Computing solenoidal eigenmodes of the vector Helmholtz equation: a novel approach,” IEEE Trans. Magn. |

30. | N. W. Ashcroft and N. D. Mermin, |

31. | To obtain |

32. | W. Shin and S. Fan, “Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers,” J. Comput. Phys. |

33. | Y. Saad and M. H. Schultz, “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM J. Sci. Stat. Comp. |

34. | P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B |

35. | E. D. Palik, ed., |

36. | D. R. Lide, ed., |

37. | R. Freund and N. Nachtigal, “QMR: a quasi-minimal residual method for non-hermitian linear systems,” Numer. Math. |

38. | R. W. Freund, G. H. Golub, and N. M. Nachtigal, “Iterative solution of linear systems,” Acta Numer. |

39. | J. W. Goodman, |

**OCIS Codes**

(000.3860) General : Mathematical methods in physics

(000.3870) General : Mathematics

(000.4430) General : Numerical approximation and analysis

**ToC Category:**

Physical Optics

**History**

Original Manuscript: July 26, 2013

Revised Manuscript: September 8, 2013

Manuscript Accepted: September 9, 2013

Published: September 18, 2013

**Citation**

Wonseok Shin and Shanhui Fan, "Accelerated solution of the frequency-domain Maxwell’s equations by engineering the eigenvalue distribution of the operator," Opt. Express **21**, 22578-22595 (2013)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-21-19-22578

Sort: Year | Journal | Reset

### References

- N. J. Champagne, J. G. Berryman, and H. M. Buettner, “FDFD: A 3D finite-difference frequency-domain code for electromagnetic induction tomography,” J. Comput. Phys.170, 830–848 (2001). [CrossRef]
- G. Veronis and S. Fan, “Overview of simulation techniques for plasmonic devices,” in “Surface Plasmon Nanophotonics,” M. Brongersma and P. Kik, eds. (Springer, 2007), pp. 169–182. [CrossRef]
- U. S. Inan and R. A. Marshall, Numerical Electromagnetics: The FDTD Method (Cambridge University, 2011). Ch. 14. [CrossRef]
- J.-M. Jin, The Finite Element Method in Electromagnetics, 2nd ed. (Wiley, 2002).
- V. Simoncini and D. B. Szyld, “Recent computational developments in Krylov subspace methods for linear systems,” Numer. Linear Algebra Appl.14, 1–59 (2007). [CrossRef]
- W. Cai, W. Shin, S. Fan, and M. L. Brongersma, “Elements for plasmonic nanocircuits with three-dimensional slot waveguides,” Adv. Mater.22, 5120–5124 (2010). [CrossRef] [PubMed]
- G. Veronis and S. Fan, “Bends and splitters in metal-dielectric-metal subwavelength plasmonic waveguides,” Appl. Phys. Lett.87, 131102–3 (2005). [CrossRef]
- L. Verslegers, P. Catrysse, Z. Yu, W. Shin, Z. Ruan, and S. Fan, “Phase front design with metallic pillar arrays,” Opt. Lett.35, 844–846 (2010). [CrossRef] [PubMed]
- J. T. Smith, “Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator,” Geophys.61, 1319–1324 (1996). [CrossRef]
- G. A. Newman and D. L. Alumbaugh, “Three-dimensional induction logging problems, Part 2: A finite-difference solution,” Geophys.67, 484–491 (2002). [CrossRef]
- C. J. Weiss and G. A. Newman, “Electromagnetic induction in a generalized 3D anisotropic earth, Part 2: The LIN preconditioner,” Geophys.68, 922–930 (2003). [CrossRef]
- V. L. Druskin, L. A. Knizhnerman, and P. Lee, “New spectral Lanczos decomposition method for induction modeling in arbitrary 3-D geometry,” Geophys.64, 701–706 (1999). [CrossRef]
- E. Haber, U. M. Ascher, D. A. Aruliah, and D. W. Oldenburg, “Fast simulation of 3D electromagnetic problems using potentials,” J. Comput. Phys.163, 150–171 (2000). [CrossRef]
- J. Hou, R. Mallan, and C. Torres-Verdin, “Finite-difference simulation of borehole EM measurements in 3D anisotropic media using coupled scalar-vector potentials,” Geophys.71, G225–G233 (2006). [CrossRef]
- R. Hiptmair, F. Kramer, and J. Ostrowski, “A robust Maxwell formulation for all frequencies,” IEEE Trans. Magn.44, 682–685 (2008). [CrossRef]
- K. Beilenhoff, W. Heinrich, and H. Hartnagel, “Improved finite-difference formulation in frequency domain for three-dimensional scattering problems,” IEEE Trans. Microwave Theory Tech.40, 540–546 (1992). [CrossRef]
- A. Christ and H. L. Hartnagel, “Three-dimensional finite-difference method for the analysis of microwave-device embedding,” IEEE Trans. Microwave Theory Tech.35, 688–696 (1987). [CrossRef]
- At the final stage of our work, we were made aware of a related work by M. Kordy, E. Cherkaev, and P. Wannamaker, “Schelkunoff potential for electromagnetic field: proof of existence and uniqueness” (to be published), where an equation similar to our Eq. (7) with s= −1 was developed.
- A. Jennings, “Influence of the eigenvalue spectrum on the convergence rate of the conjugate gradient method,” IMA J. Appl. Math.20, 61–72 (1977). [CrossRef]
- A. van der Sluis and H. van der Vorst, “The rate of convergence of conjugate gradients,” Numer. Math.48, 543–560 (1986). [CrossRef]
- S. L. Campbell, I. C. F. Ipsen, C. T. Kelley, and C. D. Meyer, “GMRES and the minimal polynomial,” BIT Numer. Math.36, 664–675 (1996). [CrossRef]
- S. Goossens and D. Roose, “Ritz and harmonic Ritz values and the convergence of FOM and GMRES,” Numer. Linear Algebra Appl.6, 281–293 (1999). [CrossRef]
- M. Benzi, G. H. Golub, and J. Liesen, “Numerical solution of saddle point problems,” Acta Numer.14, 1–137 (2005). Sec. 9.2. [CrossRef]
- B. Beckermann and A. B. J. Kuijlaars, “Superlinear convergence of conjugate gradients,” SIAM J. Numer. Anal.39, 300–329 (2001). [CrossRef]
- B. Beckermann and A. B. J. Kuijlaars, “On the sharpness of an asymptotic error estimate for conjugate gradients,” BIT Numer. Math.41, 856–867 (2001). [CrossRef]
- O. Axelsson, “Iteration number for the conjugate gradient method,” Math. Comput. Simulat.61, 421–435 (2003). [CrossRef]
- J. P. Webb, “The finite-element method for finding modes of dielectric-loaded cavities,” IEEE Trans. Microwave Theory Tech.33, 635–639 (1985). [CrossRef]
- F. Kikuchi, “Mixed and penalty formulations for finite element analysis of an eigenvalue problem in electromagnetism,” Comput. Method Appl. Mech. Eng.64, 509–521 (1987). [CrossRef]
- D. A. White and J. M. Koning, “Computing solenoidal eigenmodes of the vector Helmholtz equation: a novel approach,” IEEE Trans. Magn.38, 3420–3425 (2002). [CrossRef]
- N. W. Ashcroft and N. D. Mermin, Solid State Physics, 1st ed. (SaundersCollege, 1976), Ch. 8.
- To obtain 8/Δmin2, take the first equation of Eq. (4.29) in [32] and then multiply the extra factor μ= μ0to account for the difference between Eq. (4.17a) in [32] and Eq. (1) in the present paper.
- W. Shin and S. Fan, “Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers,” J. Comput. Phys.231, 3406–3431 (2012). [CrossRef]
- Y. Saad and M. H. Schultz, “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM J. Sci. Stat. Comp.7, 856–869 (1986). [CrossRef]
- P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B6, 4370–4379 (1972). [CrossRef]
- E. D. Palik, ed., Handbook of Optical Constants of Solids (Academic Press, 1985).
- D. R. Lide, ed., CRC Handbook of Chemistry and Physics, 88th ed. (CRC, 2007).
- R. Freund and N. Nachtigal, “QMR: a quasi-minimal residual method for non-hermitian linear systems,” Numer. Math.60, 315–339 (1991). [CrossRef]
- R. W. Freund, G. H. Golub, and N. M. Nachtigal, “Iterative solution of linear systems,” Acta Numer.1, 57–100 (1992). Secs. 2.4 and 3.3. [CrossRef]
- J. W. Goodman, Introduction to Fourier Optics(Roberts & Company Publishers, 2005), 3rd ed. Sec. 2.3.2.

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.