Improving condition number and convergence of the surface integral-equation method of moments for penetrable bodies

Luis Landesa; Marta Gómez Araújo; José Manuel Taboada; Luis Bote; Fernando Obelleiro

doi:10.1364/OE.20.017237

1. Introduction

Most practical electromagnetic scattering and radiation problems can be defined as a combination of conducting and penetrable objects. The surface integral equation (SIE) method [1] is perhaps the most powerful numerical method in the electromagnetic analysis of this kind of problems. It has demonstrated to be very accurate and efficient for the analysis of composite objects with real conductors and dielectrics in radio frequency and microwave regimes [2–10]. Recently, this approach was successfully extended to the analysis of homogenized metamaterials and plasmonic problems [11–18] in the context of nanoscience applications. However, traditional SIE methods are still restricted to electrically small problems (small in terms of the wavelength) because the implicit discretization process of the SIE results in a dense matrix equation which is very expensive to solve and to store. In order to solve electrically large problems, advanced methods like the fast multipole method (FMM) [19] and its multilevel version (MLFMA) [20,21] are required in combination with iterative solutions of the matrix equation [22–24]. Nevertheless, many of the traditional SIE formulations lead to ill-conditioned matrix equations. A high condition number may cause problems in the numerical solution, such as slow convergence of iterative solvers (also depending on the clustering of eigenvalues) and even inaccuracy of the solution when the condition number is really high.

The most typical way to improve the condition number of a matrix equation, and so its convergence in the context of an iterative algorithm, is preconditioning [25]. Some popular preconditioners are for example: the use of the diagonal elements of the matrix; the near-field, block-diagonal [26]; the incomplete LU factorization [27,28]; the multiplicative preconditioner using Calderon identities [29,30]; or inner/outer GMRES preconditioners [31,32]. Effective preconditioners are, however, computationally costly.

An alternative way to improve the condition of the matrix is a proper choice of the SIE to be solved. The main reason for the poor conditioning of the SIE matrix is that the unknowns (the electric and magnetic surface current densities), and the integral equations (the electric and magnetic field integral equations), are of the very different scales. Usually the balance between the unknowns and between the equations is improved, for example, by multiplying the magnetic field integral equation (MFIE) with the wave impedance [33]; or scaling the magnetic, or electric, surface current density by the wave impedance [34–36] or both [37,38]. In these cases, the conditioning of the matrix is improved by using proper weighting coefficients and combination of the equations [8]. In fact, the widespread combined tangential formulation (CTF) can be seen as a weighted version of the original Poggio–Miller–Chang–Harrington–Wu–Tsai (PMCHWT) formulation [34], in which the balance between the diagonal blocks of the matrix has been improved. And the same goes for the combined normal formulation (CNF), which can be seen as a variation of the original Müller formulation [2,39]. Anyway, the resultant matrices in all these SIE formulations are still unbalanced because, for example, when balancing the diagonal blocks of the matrices, the off-diagonal blocks become unbalanced. Otherwise, modifying the formulation through the appropriate selection of combination coefficients to improve block balance may affect other desirable properties. This is the case of CTF and PMCHWT when dealing with plasmonic problems. While CTF has a much better convergence than PMCHWT, this is at the expense of a lower accuracy, as has been demonstrated in [18].

In a recent paper [40], an efficient method to improve the balance between the unknowns and the equations and, hence, also between the matrix elements has been presented. The method is based on the use of normalized field quantities and unknowns together with carefully chosen scaling factors, achieving a reduction in the condition number of the SIE matrices of several orders of magnitude.

In this paper we propose a similar solution, although in this case the problem is addressed from novel and a strictly numerical point of view. The SIE matrix formulation is not modified anymore; namely, the balance between matrix elements is straightforwardly addressed during the iterative solution process by means of a left and right (L&R) diagonal matrix preconditioner. This allows improving the matrix balance independently of the formulation, which can now be selected solely in terms of accuracy for a given application. The proposed preconditioning formulation is developed in a general way, so that particular expressions can be derived for any particular SIE formulation.

Numerical results show that, owing to the good balance between the matrix elements, the proposed preconditioner provides a significant reduction in the condition number. In this paper, the proposed preconditioner has been tested for the conventional SIE formulations previously mentioned: PMCHWT formulation, CTF, the normal Müller formulation and CNF, and for the electric and magnetic current combined field integral equation (JMCFIE) [9], although it can be straight forwardly extended to any other formulation.

2. Surface integral-equation formulations for composite penetrable bodies

Let us consider a homogeneous penetrable object in a homogeneous (unbounded) medium. The exterior region is denoted by $R_{1}$ and the interior region is denoted by $R_{2}$ . Let us define with $ε_{i} = ε_{r i} ε_{0}$ the complex permittivity and with $μ_{i} = μ_{r i} μ_{0}$ the complex permeability of region $R_{i}$ with i = 1,2. $ε_{r i}$ and $μ_{r i}$ are the complex valued relative permittivity and permeability constants of region $R_{i}$ and $ε_{0}$ , $μ_{0}$ are the constitutive parameters of vacuum. A time harmonic variation exp(jωt) is assumed and suppressed from the formulation. Let us denote with S the interface surface between regions $R_{1}$ and $R_{2}$ . ${\hat{n}}_{i}$ denotes the unit vector normal to S and pointing toward $R_{i}$ , hence ${\hat{n}}_{2} = - {\hat{n}}_{1}$ .

According to the equivalence principle [41] the electric field integral equation (EFIE) and the magnetic field integral equation (MFIE) can be formulated separately in each region $R_{i}$ . Two alternative formulations can be derived depending on how the fields are projected onto the surfaces surrounding this region. The tangential (T) equations in $R_{i}$ are given by the following expressions [1,8]:

{T-EFIE}_{i} : {(η_{i} ℒ_{i} (J_{i}) - K_{i} (M_{i}))}_{t a n} + \frac{1}{2} {\hat{n}}_{i} \times M_{i} = {(E_{i}^{i n c})}_{t a n}

{T-MFIE}_{i} : {(K_{i} (J_{i}) + η_{i}^{- 1} ℒ_{i} (M_{i}))}_{\tan} - \frac{1}{2} {\hat{n}}_{i} \times J_{i} = {(H_{i}^{i n c})}_{t a n}

In a similar way, the normal (N) equations in $R_{i}$ can be written as

{N-EFIE}_{i} : {\hat{n}}_{i} \times (η_{i} ℒ_{i} (J_{i}) - K_{i} (M_{i})) - \frac{1}{2} M_{i} = {\hat{n}}_{i} \times E_{i}^{i n c}

{N-MFIE}_{i} : {\hat{n}}_{i} \times (K_{i} (J_{i}) + η_{i}^{- 1} ℒ_{i} (M_{i})) + \frac{1}{2} J_{i} = {\hat{n}}_{i} \times H_{i}^{i n c}

Here $η_{i}$ is the intrinsic impedance in medium $R_{i}$ . The integro-differential $ℒ_{i}$ and $K_{i}$ operators are defined as

ℒ_{i} (X_{i}) = j k_{i} [\int_{S}^{} X_{i} (r') G_{i} (r, r') d S' + \frac{1}{k_{i}^{2}} \nabla \int_{S}^{} \nabla' \cdot X_{i} G_{i} (r, r') d S']

K_{i} (X_{i}) = \int_{S, P V}^{} X_{i} (r') \times \nabla G_{i} (r, r') d S'

where r denotes the observation points approaching to

S

from the inner of region

R_{i}

and _{$r' \in S$} denotes the source points.

\nabla' \cdot

is the divergence in the primed (source) coordinates, PV denotes the principal value of the integral in Eq. (6),

k_{i}

is the wavenumber in

R_{i}

, and

G_{i} (r, r') = \frac{\exp (- j k_{i} | r - r' |)}{4 π | r - r' |}

is the homogeneous Green’s function in

R_{i}

.

E_{i}^{i n c} (r')

and

H_{i}^{i n c} (r')

are the incident fields due to the sources located inside region

R_{i}

, and

J_{i} (r')

and

M_{i} (r')

are the equivalent electric and magnetic currents placed over the boundary surface S, as given by the equivalence principle for each region

R_{i}

. According to the boundary conditions (continuity of tangential field components across the boundary surface S), the surface current densities in both sides of S should satisfy

J_{1} = - J_{2}

and

M_{1} = - M_{2}

, thus the problem can be described using only one set of equivalent electric and magnetic currents on S, namely

J = J_{1} = - J_{2}

and

M = M_{1} = - M_{2}

.

Among the various possibilities of combination of the integral equations Eqs. (1)–(4), we have chosen the procedure of [7] to derive two stable and well tested SIEs. The $1 / η T-EFIE$ is combined with the N-MFIE, and the N-EFIE is combined with the $η T-MFIE$ leading to the electric current combined field integral equation (JCFIE) and the magnetic current combined field integral equation (MCFIE) in region $R_{i}$ :

{JCFIE}_{i} = a_{i} η_{i}^{- 1} {T-EFIE}_{i} + b_{i} {N-MFIE}_{i}

{MCFIE}_{i} = - c_{i} {N-EFIE}_{i} + d_{i} η_{i} {T-MFIE}_{i}

The JCFIE and MCFIE equations for regions $R_{1}$ and $R_{2}$ are then combined respectively to obtain two independent integral equations for the two unknown functions J and M:

\begin{array}{l} JCFIE: a_{1} {(ℒ_{1} (J) - η_{1}^{- 1} K_{1} (M))}_{t a n} + a_{2} {(ℒ_{2} (J) - η_{2}^{- 1} K_{2} (M))}_{t a n} \\ + b_{1} \hat{n} \times (K_{1} (J) + η_{1}^{- 1} ℒ_{1} (M)) - b_{2} \hat{n} \times (K_{2} (J) + η_{2}^{- 1} ℒ_{2} (M)) \\ + \frac{1}{2} (a_{1} η_{1}^{- 1} - a_{2} η_{2}^{- 1}) \hat{n} \times M + \frac{1}{2} (b_{1} + b_{2}) J \\ = a_{1} η_{1}^{- 1} {(E_{1}^{i n c})}_{t a n} - a_{2} η_{2}^{- 1} {(E_{2}^{i n c})}_{t a n} + b_{1} \hat{n} \times H_{1}^{i n c} + b_{2} \hat{n} \times H_{2}^{i n c}, r \in S \end{array}

\begin{array}{l} MCFIE : - c_{1} \hat{n} \times (η_{1} ℒ_{1} (J) - K_{1} (M)) + c_{2} \hat{n} \times (η_{2} ℒ_{2} (J) - K_{2} (M)) \\ + d_{1} {(η_{1} K_{1} (J) + ℒ_{1} (M))}_{t a n} + d_{2} {(η_{2} K_{2} (J) + ℒ_{2} (M))}_{t a n} \\ + \frac{1}{2} (c_{1} + c_{2}) M - \frac{1}{2} (d_{1} η_{1} - d_{2} η_{2}) \hat{n} \times J \\ = - c_{1} \hat{n} \times E_{1}^{i n c} - c_{2} \hat{n} \times E_{2}^{i n c} + d_{1} η_{1} {(H_{1}^{i n c})}_{t a n} - d_{2} η_{2} {(H_{2}^{i n c})}_{t a n}, r \in S \end{array}

3. Discretization of the integral equations

To obtain the equivalent currents the standard method of moments (MoM) procedure is applied to the two previous equations. The equivalent currents are expanded into a sum of known vector basis functions $f_{n}$ in the form

J = \sum_{n} J_{n} f_{n}; M = \sum_{n} M_{n} f_{n}; r \in S_{i j}

where

J_{n}

and

M_{n}

are the unknown expansion complex coefficients. Substituting Eq. (12) into Eq. (10) and Eq. (11) and applying the Galerkin testing procedure, a system of linear equations is derived from the integral equations, as follows:

\begin{array}{l} a_{1} \sum_{n \in S} (A_{m n}^{1} J_{n} - η_{1}^{- 1} B_{m n}^{1} M_{n}) + a_{2} \sum_{n \in S} (A_{m n}^{2} J_{n} - η_{2}^{- 1} B_{m n}^{2} M_{n}) \\ + b_{1} \sum_{n \in S} (B'_{m n}^{1} J_{n} + η_{1}^{- 1} A'_{m n}^{1} M_{n}) - b_{2} \sum_{n \in S} (B'_{m n}^{2} J_{n} + η_{2}^{- 1} A'_{m n}^{2} M_{n}) \\ + \frac{1}{2} \sum_{n \in S} [(a_{1} η_{1}^{- 1} - a_{2} η_{2}^{- 1}) I'_{m n}^{} M_{n}] + \frac{1}{2} \sum_{n \in S} [(b_{1} + b_{2}) I_{m n}^{} J_{n}] \\ = E_{m}^{} + H'_{m}^{}, m \in S \end{array}

\begin{array}{l} - c_{1} \sum_{n \in S} (η_{1} A'_{m n}^{1} J_{n} - B'_{m n}^{1} M_{n}) + c_{2} \sum_{n \in S} (η_{2} A'_{m n}^{2} J_{n} - B'_{m n}^{2} M_{n}) \\ + d_{1} \sum_{n \in S} (η_{1} B_{m n}^{1} J_{n} + A_{m n}^{1} M_{n}) + d_{2} \sum_{n \in S} (η_{2} B_{m n}^{2} J_{n} + A_{m n}^{2} M_{n}) \\ + \frac{1}{2} \sum_{n \in S} [(c_{1} + c_{2}) I_{m n} M_{n}] - \frac{1}{2} \sum_{n \in S} [(d_{1} η_{1} - d_{2} η_{2}) I'_{m n}^{} J_{n}] \\ = E'_{m}^{} + H_{m}^{}, m \in S \end{array}

where the following quantities have been defined:

A_{m n}^{i} = \int_{Δ_{m}} f_{m} \cdot ℒ_{i} (f_{n}) d S

B_{m n}^{i} = \int_{Δ_{m}} f_{m} \cdot K_{i} (f_{n}) d S

A'_{m n}^{i} = \int_{Δ_{m}} f_{m} \cdot {\hat{n}}_{m} \times ℒ_{i} (f_{n}) d S

B'_{m n}^{i} = \int_{Δ_{m}} f_{m} \cdot {\hat{n}}_{m} \times K_{i} (f_{n}) d S

I_{m n}^{} = \int_{Δ_{m}} f_{m} \cdot f_{n} d S

I'_{m n}^{} = \int_{Δ_{m}} f_{m} \cdot {\hat{n}}_{m} \times f_{n} d S .

E_{m}^{} = \int_{Δ_{m}} f_{m} \cdot [a_{1} {(η_{i}^{- 1} E_{i}^{i n c})}_{t a n} - a_{2} {(η_{2}^{- 1} E_{j}^{i n c})}_{t a n}] d S

H_{m}^{} = \int_{Δ_{m}} f_{m} \cdot [d_{1} {(η_{1} H_{1}^{i n c})}_{t a n} - d_{2} {(η_{2} H_{2}^{i n c})}_{t a n}] d S

E'_{m}^{} = - \int_{Δ_{m}} f_{m} \cdot [c_{1} {\hat{n}}_{m} \times E_{1}^{i n c} + c_{2} {\hat{n}}_{m} \times E_{2}^{i n c}] d S

H'_{m}^{} = \int_{Δ_{m}} f_{m} \cdot [b_{1} {\hat{n}}_{m} \times H_{1}^{i n c} + b_{2} {\hat{n}}_{m} \times H_{2}^{i n c}] d S

Regrouping terms in Eq. (13) and Eq. (14), the system of linear equations can be expressed in the usual dense matrix system form:

\bar{Z} I = V

where the impedance matrix

\bar{Z}

can be written as

{\bar{Z}}_{}^{} = [\begin{matrix} {\bar{Z}}_{}^{1 J} & {\bar{Z}}_{}^{1 M} \\ {\bar{Z}}_{}^{2 J} & {\bar{Z}}_{}^{2 M} \end{matrix}]

with the elements given by

{\bar{Z}}_{m n}^{1 J} = a_{1} A_{m n}^{1} + a_{2} A_{m n}^{2} + b_{1} B'_{m n}^{1} - b_{2} B'_{m n}^{2} + \frac{1}{2} (b_{1} + b_{2}) I_{m n}

{\bar{Z}}_{m n}^{1 M} = - a_{1} η_{1}^{- 1} B_{m n}^{1} - a_{2} η_{2}^{- 1} B_{m n}^{2} + b_{1} η_{1}^{- 1} A'_{m n}^{1} - b_{2} η_{2}^{- 1} A'_{m n}^{2} + \frac{1}{2} (a_{1} η_{1}^{- 1} - a_{2} η_{2}^{- 1}) I'_{m n}

{\bar{Z}}_{m n}^{2 J} = d_{1} η_{1} B_{m n}^{1} + d_{2} η_{2} B_{m n}^{2} - c_{1} η_{1} A'_{m n}^{1} + c_{2} η_{2} A'_{m n}^{2} - \frac{1}{2} (d_{1} η_{1} - d_{2} η_{2}) I'_{m n}

{\bar{Z}}_{m n}^{2 M} = d_{1} A_{m n}^{1} + d_{2} A_{m n}^{2} + c_{1} B'_{m n}^{1} - c_{2} B'_{m n}^{2} + \frac{1}{2} (c_{1} + c_{2}) I_{m n}

Similarly we define $I = {[\begin{matrix} J & M \end{matrix}]}^{T}$ and $V = {[\begin{matrix} V^{1} & V^{2} \end{matrix}]}^{T}$ , with superscript T indicating vector transposition, and where

V_{m}^{1} = E_{m}^{} + H'_{m}^{}

V_{m}^{2} = E'_{m}^{} + H_{m}^{} .

Different known formulations can be obtained depending on the selection of the complex combination parameters $a_{i}, b_{i}, c_{i}, d_{i}$ for $i = 1, 2$ (see Table 1 ).

Table 1. Combination parameters $a_{i}, b_{i}, c_{i}, d_{i}$ and $O_{n o r m} (\bar{Z})$ for the five formulations considered

View Table | View all tables in this article

4. Numerical balance of the impedance matrix

The impedance matrix in Eq. (26) has a four block structure. At this point, we will address the balance between the elements inside this matrix. Looking at the definition of the elements given in Eqs. (27)–(30), and considering that the quantities given in Eqs. (15)–(20) are the same order of magnitude, namely, $O (A_{m n}^{i}) \approx O (A_{m n}^{' i}) \approx O (B_{m n}^{i}) \approx O (B_{m n}^{' i}) \approx O (I_{m n}^{i}) \approx O (I_{m n}^{' i}) \approx Θ$ for $i = 1, 2$ , the order of the four blocks of $\bar{Z}$ becomes

O (\bar{Z}) = O ([\begin{matrix} {\bar{Z}}_{}^{1 J} & {\bar{Z}}_{}^{1 M} \\ {\bar{Z}}_{}^{2 J} & {\bar{Z}}_{}^{2 M} \end{matrix}]) \approx (\begin{matrix} Θ f (a_{1}, a_{2}, b_{1}, b_{2}) & Θ g (a_{1}, a_{2}, b_{1}, b_{2}, η_{1}, η_{2}) \\ Θ h (c_{1}, c_{2}, d_{1}, d_{2}, η_{1}, η_{2}) & Θ f (c_{1}, c_{2}, d_{1}, d_{2}) \end{matrix})

where f, g and h are functions which note linear combinations of the combination parameters (_{$a_{i}, b_{i}, c_{i}, d_{i}$}) and the intrinsic impedances (_{$η_{i}$}) of the media on both sides of surface _$S$ (i = 1,2) according to the expressions of Eqs. (27)–(30).

Thereafter, since the combination parameters $a_{i}, b_{i}, c_{i}, d_{i}$ are defined in the same way inside each region (although they might actually have different values when they depend on the constitutive parameters of the media), we can assume that they are the same order of magnitude. So, let us assume that $O (a_{1}) \approx O (a_{2}) \approx a$ , $O (b_{1}) \approx O (b_{2}) \approx b$ , $O (c_{1}) \approx O (c_{2}) \approx c$ , and $O (d_{1}) \approx O (d_{2}) \approx d$ . Similarly, we can assume that $O (η_{1}) \approx O (η_{2}) \approx η$ .

With the previous considerations we can rewrite the expression of Eq. (33) as follows:

O ([\begin{matrix} {\bar{Z}}_{}^{1 J} & {\bar{Z}}_{}^{1 M} \\ {\bar{Z}}_{}^{2 J} & {\bar{Z}}_{}^{2 M} \end{matrix}]) \approx (\begin{matrix} Θ (a + b) & \frac{Θ (a + b)}{η} \\ Θ (c + d) η & Θ (c + d) \end{matrix})

In order to appreciate/study the imbalance in the order of magnitude between the different blocks of matrix $\bar{Z}$ , let us consider a normalized version of the above expression, in which the order of magnitude has been normalized with respect to the order of magnitude of ${\bar{Z}}^{1 J}$ :

O_{n o r m} (\bar{Z}) = \frac{O (\bar{Z})}{O ({\bar{Z}}_{}^{1 J})} \approx (\begin{matrix} 1 & \frac{1}{η} \\ \frac{(c + d)}{(a + b)} η & \frac{(c + d)}{(a + b)} \end{matrix})

In the next table we present $O_{n o r m} (\bar{Z})$ for the most usual SIE formulations.

Looking at Table 1, it becomes clear that in the PMCHWT formulation the diagonal blocks differ by a factor of $η^{2}$ , while the off-diagonal blocks are the same order and differ by a factor of η with the diagonal blocks. In the case of CTF, the diagonal blocks are the same order (which is the main objective of this formulation, thus greatly improving the convergence with respect to PMCHWT), but consequently the off-diagonal blocks differ by a factor of $η^{2}$ , and by a factor of η with respect to the diagonal blocks. Otherwise, CNF behaves similar to CTF, with equal diagonal blocks but with off-diagonal blocks differing by a factor of $η^{2}$ and by a factor of η with regard to the diagonal blocks. Finally, Müller formulation behaves in the same way that the PMCHWT formulation. It can be concluded that in all formulations the strong imbalance between the impedance matrix blocks leads to high condition numbers, which causes bad convergence and/or lose of accuracy.

To circumvent the previous problems, it would be desirable to simultaneously balance the four blocks of each impedance submatrix, thus improving the condition number and convergence of the overall matrix system. Based on Eq. (35), such well-balanced submatrices could be straightforwardly obtained by applying the following block-wise Hadamard product:

\tilde{Z} = [\begin{matrix} {\bar{Z}}_{}^{1 J} & {\bar{Z}}_{}^{1 M} \\ {\bar{Z}}_{}^{2 J} & {\bar{Z}}_{}^{2 M} \end{matrix}] ⊙ (\begin{matrix} 1 & η \\ \frac{(a + b)}{(c + d)} \frac{1}{η} & \frac{(a + b)}{(c + d)} \end{matrix})

Different designations can be made for a, b, c, d and η. In this work, we have selected for these parameters the outer medium values: a = a₁, b = b₁, c = c₁, d = d₁, and η = η₁, although other possibilities also work properly, such as choosing the arithmetic or geometric means or the vacuum equivalent values, among others.

5. Left and right preconditioner

The previous balancing scheme of the impedance matrix can be straightforwardly included in the matrix equation of Eq. (25) for any given formulation by applying the following change of basis:

\tilde{I} = M_{R}^{- 1} I

\tilde{V} = {\bar{M}}_{L} V

where

{\bar{M}}_{L}

and

{\bar{M}}_{R}

are diagonal matrices (2N × 2N) as follows:

{\bar{M}}_{L}^{} = diag {{[α_{11}^{} ... α_{11}^{}]}_{1 \times N}, {[α_{22}^{} ... α_{22}^{}]}_{1 \times N}}

{\bar{M}}_{R}^{} = diag {{[β_{11}^{} ... β_{11}^{}]}_{1 \times N}, {[β_{22}^{} ... β_{22}^{}]}_{1 \times N}}

With these changes Eq. (25) can be written as

\tilde{Z} \tilde{I} = \tilde{V}

where we have defined

\tilde{Z} = {\bar{M}}_{L} \bar{Z} {\bar{M}}_{R}

The diagonal values of the diagonal matrices in Eq. (39) and Eq. (40), $α_{11}^{}, α_{22}^{}, β_{11}^{}$ and $β_{22}^{}$ , can be obtained by comparing Eq. (42) to Eq. (36). For simplicity in the derivation of the wanted elements, we will take $α_{11}^{} = β_{11}^{} = 1$ . So, the coefficients for the left and right diagonal matrices that balance the impedance matrix blocks are given by

\begin{array}{l} α_{11}^{} = β_{11}^{} = 1 \\ β_{22}^{} = η \\ α_{22}^{} = \frac{a + b}{c + d} \frac{1}{η} \end{array}

In Table 2 we summarize the combination and preconditioning coefficients for the usual formulations, where μ = μ₁ and ε = ε₁ according to the previous choice made for a, b, c, d and η.

Table 2. Combination ( $a, b, c, d$ ) and preconditioning ( $α_{11}, α_{22}, β_{11}, β_{22}$ ) coefficients for the five formulations considered

View Table | View all tables in this article

Using the above definitions, Eq. (41) is a well-balanced system with low condition number which can be easily solved using any iterative scheme. Once Eq. (41) is solved, the solution to the original electromagnetic problem is finally obtained from Eq. (37) as $I = {\bar{M}}_{R} \tilde{I}$ .

Regarding the computational cost, the matrices involved in the preconditioner are diagonal, so the required number of multiplications become O(N), N being the number of unknowns.

6. Numerical examples

The condition number of the matrix corresponding to each conventional SIE formulation under consideration is studied next. Spheres of different size and different permittivity values have been considered. The Rao-Wilton-Glisson (RWG) [42] basis functions have been used for the discretization of the problems. The analyses have been carried out at an exterior (vacuum) wavelength of λ₀ = 546 nm and a mean mesh size of λ₀/15 has been considered in all cases except for the smallest sphere, where a mesh refinement was required in order to guarantee the geometrical sphericity.

Figure 1 shows the matrix condition number with and without preconditioner for several dielectric spheres characterized by ε_r = 2.1 that have been analyzed using the five formulations. The condition number is represented as a function of k₀r, where k₀ is the wave number in vacuum and r the radius of the sphere. The results of Fig. 2 and Fig. 3 correspond to a gold plasmonic sphere with ε_r = −5.84−j2.11 [43] and a high-contrast dielectric sphere of ε_r = 25, respectively.

Fig. 1 Condition number without/with preconditioner vs. sphere size for dielectric spheres of ε_r = 2.1 analyzed with the five formulations.

Download Full Size | PDF

Fig. 2 Condition number without/with preconditioner vs. sphere size for Au spheres (ε_r = −5.84−j2.11) analyzed with the five formulations.

Download Full Size | PDF

Fig. 3 Condition number without/with preconditioner vs. sphere size for dielectric spheres of ε_r = 25 analyzed with the five formulations.

Download Full Size | PDF

Despite the electrical properties of the sphere material, the general conclusion that can be extracted for all formulations is that a strong reduction in the condition number is obtained when the proposed preconditioning scheme is applied. Of course, other considerations regarding the different behavior of each formulation depending on the material properties could be pointed out. However, this issue has been the subject of previous studies [8,17,18,23,39,40,44] and it is out of the scope of this work, whose main objective is to show that the proposed L&R preconditioner provides a reduction of several orders of magnitude in the condition number when considering any of the usual SIE formulations and different permittivity values as is illustrated in Figs. 1, 2 and 3.

Due to the reduction in the condition number, the proposed preconditioner makes even possible to reach fast convergence with typically slowly converging formulations such as Müller and especially PMCHWT. To illustrate this fact, the analysis of a gold sphere involving 24000 unknowns (larger sample of Fig. 2) is presented next as an example of an advantageous use of the preconditioner. PMCHWT would be the preferable choice to solve such a problem, because it has demonstrated to be the most accurate formulation among the considered approaches to deal with plasmonic problems [18]. However, when the electrical size increases and the direct solution is not possible, iterative schemes are required and the analysis through PMCHWT becomes unfeasible due to its well-known convergence problems. In the case of conventional materials, CTF is usually proposed as an alternative to PMCHWT since its weighting coefficients lead to better convergence [8,23]. This behavior is clearly illustrated in Fig. 4 , where the plasmonic problem has been solved iteratively by means of GMRES (restart 30). The convergence problems of PMCHWT (without preconditioner) are overcome when CTF is applied. Nevertheless, additionally to convergence, the accuracy of the obtained results must be also considered when it comes to selecting a formulation. With this aim, the following definition for the normalized root mean square (RMS) error of the bistatic radar cross section (RCS, σ) with respect to the analytical Mie’s series result [45] has been considered:

e_{r m s} = \frac{\sqrt{\sum {(σ_{M i e} - σ)}^{2} / N}}{\max (σ_{M i e})}

with N the number of RCS samples. An error of 2.01⋅10⁻⁴ has been obtained for PMCHWT with both direct and precorrected iterative solutions, while an error of 0.013 has been obtained for CTF. These values show that, unlike in conventional problems, CTF leads to an important lack of accuracy in the context of plasmonic problems despite its good iterative convergence. This is in accordance with the behavior reported in the thorough comparison of [18].

Fig. 4 Iterative convergence of PMCHWT formulation and CTF solved with GMRES(30) with and without preconditioner.

Download Full Size | PDF

Under the above conditions, the use of a low-cost L&R preconditioner that improves the convergence without affecting the nature of the SIE formulation constitutes the best option. Figure 4 shows that the pursued improvement of the PMCHWT convergence is attained with the preconditioner proposed in this work, enabling the accurate analysis of the problem. In fact, even taking the result after only 8 outer iterations (relative residue below 10⁻³) an error level of 2.06⋅10⁻⁴ is already obtained. As expected, improvement on the iteration convergence is not observed with the preconditioned CTF since convergence essentially depends on the good balance of the diagonal blocks, which is already satisfied by CTF without applying the preconditioner.

7. Conclusion

An efficient left and right (L&R) diagonal preconditioner has been presented to balance the different-scale blocks of the impedance matrix in the context of SIE-MoM formulations. The good balance between matrix elements and between electric and magnetic unknowns is straightforwardly obtained during the iterative solution, without altering the formulation at all. This allows fixing the condition number and convergence issues regardless the applied formulation, which can now be selected solely in terms of accuracy for a given application, without worrying at all about its convergence properties. The results demonstrate that the preconditioner leads to significantly better conditioned systems, making even possible to reach fast convergence with typically very slowly converging formulations. Indeed, this is the case with PMCHWT. Although this formulation has demonstrated to be very accurate in plasmonics, its lack of convergence prevented its use in the case of large problems where the direct solution is not feasible. The proposed preconditioner overcomes this difficulty achieving the desirable combination of fast convergence and high accuracy with the PMCHWT formulation.

Finally, the proposed preconditioning scheme is completely general. Proper expressions can be easily derived for any particular SIE formulation. Besides, due to its simplicity and low cost of O(N), it can be easily and efficiently integrated into any existing MoM implementation.

Acknowledgments

This work was supported by the Spanish Government and FEDER (projects TEC2011-28784-C02-01, TEC2011-28784-C02-02, CONSOLIDER-INGENIO2010 CSD2008-00068).

References and links

1. B. M. Kolundzija and A. R. Djordjevic, Electromagnetic Modeling of Composite Metallic and Dielectric Structures (Artech House, 2002).

2. C. Müller, Foundations of the Mathematical Theory of Electromagnetic Waves (Springer, 1969).

3. A. J. Poggio and E. K. Miller, Computer Techniques for Electromagnetics (Permagon, 1973), Chap. 4.

4. Y. Chang and R. F. Harrington, “A surface formulation for characteristic modes of material bodies,” IEEE Trans. Antenn. Propag. AP-25(6), 789–795 (1977). [CrossRef]

5. T. K. Wu and L. L. Tsai, “Scattering from arbitrarily-shaped lossy dielectric bodies of revolution,” Radio Sci. 12(5), 709–718 (1977). [CrossRef]

6. S. M. Rao and D. R. Wilton, “E-field, H-field, and combined field solution for arbitrarily shaped three-dimensional dielectric bodies,” Electromagetics 10(4), 407–421 (1990). [CrossRef]

7. M. S. Yeung, “Single integral equation for electromagnetic scattering by three-dimensional homogeneous dielectric objects,” IEEE Trans. Antenn. Propag. 47(10), 1615–1622 (1999). [CrossRef]

8. P. Ylä-Oijala, M. Taskinen, and S. Järvenpää, “Surface integral equation formulations for solving electromagnetic scattering problems with iterative methods,” Radio Sci. 40(6), RS6002 (2005). [CrossRef]

9. P. Ylä-Oijala and M. Taskinen, “Application of combined field integral equation for electromagnetic scattering by dielectric and composite objects,” IEEE Trans. Antenn. Propag. 53(3), 1168–1173 (2005). [CrossRef]

10. P. Ylä-Oijala, M. Taskinen, and J. Sarvas, “Surface integral equation method for general integral equation method for general composite metallic and dielectric structures with junctions,” Prog. Electromagn. Res. 52, 81–108 (2005). [CrossRef]

11. D. L. Smith, L. N. Medgyesi-Mitschang, and D. W. Forester, “Surface integral equation formulations for left-handed materials,” Prog. Electromagn. Res. 51, 27–48 (2005). [CrossRef]

12. Y. A. Liu and W. C. Chew, “Stability of surface integral equation for left-handed materials,” IET Microwaves Antenn. Propag. 1(1), 84–89 (2007). [CrossRef]

13. B. Gallinet, A. M. Kern, and O. J. F. Martin, “Accurate and versatile modeling of electromagnetic scattering on periodic nanostructures with a surface integral approach,” J. Opt. Soc. Am. A 27(10), 2261–2271 (2010). [CrossRef] [PubMed]

14. Ö. Ergül and L. Gürel, “Efficient solutions of metamaterial problems using a low-frequency multilevel fast multipole algorithm,” Prog. Electromagn. Res. 108, 81–99 (2010). [CrossRef]

15. J. Rivero, J. M. Taboada, L. Landesa, F. Obelleiro, and I. García-Tuñón, “Surface integral equation formulation for the analysis of left-handed metamaterials,” Opt. Express 18(15), 15876–15886 (2010).

16. J. M. Taboada, J. Rivero, F. Obelleiro, M. G. Araújo, and L. Landesa, “Method-of-moments formulation for the analysis of plasmonic nano-optical antennas,” J. Opt. Soc. Am. A 28(7), 1341–1348 (2011). [CrossRef] [PubMed]

17. M. G. Araújo, J. M. Taboada, J. Rivero, and F. Obelleiro, “Comparison of surface integral equations for left-handed materials,” Prog. Electromagn. Res. 118, 425–440 (2011). [CrossRef]

18. M. G. Araújo, J. M. Taboada, D. M. Solís, J. Rivero, L. Landesa, and F. Obelleiro, “Comparison of surface integral equation formulations for electromagnetic analysis of plasmonic nanoscatterers,” Opt. Express 20(8), 9161–9171 (2012). [CrossRef] [PubMed]

19. R. Coifman, V. Rokhlin, and S. Wandzura, “The fast multipole method for the wave equation: a pedestrian prescription,” IEEE Antenn. Propag. Mag. 35(3), 7–12 (1993). [CrossRef]

20. J. M. Song and W. C. Chew, “Multilevel fast multipole algorithm for solving combined field integral equations of electromagnetic scattering,” Microw. Opt. Technol. Lett. 10(1), 14–19 (1995). [CrossRef]

21. J. M. Song, C. C. Lu, W. C. Chew, and S. Lee, “Fast Illinois solver code (FISC),” IEEE Antenn. Propag. Mag. 40(3), 27–34 (1998). [CrossRef]

22. K. C. Donepudi, J.-M. Jin, and W. C. Chew, “A higher order multilevel fast multipole algorithm for scattering from mixed conducting/dielectric bodies,” IEEE Trans. Antenn. Propag. 51(10), 2814–2821 (2003). [CrossRef]

23. Ö. Ergül and L. Gürel, “Comparison of integral-equation formulations for the fast and accurate solution of scattering problems involving dielectric objects with the multilevel fast multipole algorithm,” IEEE Trans. Antenn. Propag. 57(1), 176–187 (2009). [CrossRef]

24. M. G. Araújo, J. M. Taboada, J. Rivero, D. M. Solís, and F. Obelleiro, “Solution of large-scale plasmonic problems with the multilevel fast multipole algorithm,” Opt. Lett. 37(3), 416–418 (2012). [CrossRef] [PubMed]

25. Y. Saad, Iterative Methods for Sparse Linear Systems (PWS Publishing, 1996).

26. J. M. Song, C. C. Lu, and W. C. Chew, “Multilevel fast multipole algorithm for electromagnetic scattering by large complex objects,” IEEE Trans. Antenn. Propag. 45(10), 1488–1493 (1997). [CrossRef]

27. K. Sertel and J. L. Volakis, “Incomplete LU preconditioner for FMM implementation,” Microw. Opt. Technol. Lett. 26(4), 265–267 (2000). [CrossRef]

28. J. Lee, J. Zhang, and C.-C. Lu, “Incomplete LU preconditioner for large scale dense complex linear systems from electromagnetic wave scattering problems,” J. Comput. Phys. 185(1), 158–175 (2003). [CrossRef]

29. R. J. Adams, “Physical and analytical properties of a stabilized electric field integral equation,” IEEE Trans. Antenn. Propag. 52(2), 362–372 (2004). [CrossRef]

30. F. P. Andriulli, K. Cools, H. Bagci, F. Olyslager, A. Buffa, S. Christiansen, and E. Michielssen, “A multiplicative calderon preconditioner for the electric field integral equation,” IEEE Trans. Antenn. Propag. 56(8), 2398–2412 (2008). [CrossRef]

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

32. J. M. Bértolo, M. G. Araújo, J. M. Taboada, L. Landesa, F. Obelleiro, and J. L. Rodríguez, “Extended near field preconditioner for the analysis of large problems using the Nested-FMM-FFT algorithm,” Microw. Opt. Technol. Lett. 53(2), 430–433 (2011). [CrossRef]

33. X.-Q. Sheng, J.-M. Jin, J. Song, W. C. Chew, and C.-C. Lu, “Solution of combined-field integral equation using multilevel fast multipole algorithm for scattering by homogeneous bodies,” IEEE Trans. Antenn. Propag. 46(11), 1718–1726 (1998). [CrossRef]

34. J. R. Mautz and R. F. Harrington, “Electromagnetic scattering from a homogeneous material body of revolution,” Arch. Elektron. Ubertragungstechn. (Electron. Commun.) 33, 71–80 (1979).

35. L. N. Medgyesi-Mitschang, J. M. Putnam, and M. B. Gedera, “Generalized method of moments for three-dimensional penetrable scatterers,” J. Opt. Soc. Am. A 11(4), 1383–1398 (1994). [CrossRef]

36. X.-Q. Sheng, J.-M. Jin, J. Song, C.-C. Lu, and W. C. Chew, “On the formulation of hybrid finite-element and boundary-integral methods for 3-D scattering,” IEEE Trans. Antenn. Propag. 46(3), 303–311 (1998). [CrossRef]

37. A. Zhu, S. D. Gedney, and J. L. Visher, “A study of combined field formulations for material scattering for a locally corrected Nyström discretization,” IEEE Trans. Antenn. Propag. 53(12), 4111–4120 (2005). [CrossRef]

38. S. Chen, J.-S. Zhao, and W. C. Chew, “Analyzing low-frequency electromagnetic scattering from a composite object,” IEEE Trans. Geosci. Rem. Sens. 40(2), 426–433 (2002). [CrossRef]

39. P. Ylä-Oijala and M. Taskinen, “Well-conditioned Müller formulation for electromagnetic scattering by dielectric objects,” IEEE Trans. Antenn. Propag. 53(10), 3316–3323 (2005). [CrossRef]

40. P. Ylä-Oijala and M. Taskinen, “Improving conditioning of electromagnetic surface integral equations using normalized field quantities,” IEEE Trans. Antenn. Propag. 55(1), 178–185 (2007). [CrossRef]

41. R. F. Harrington, Time-Harmonic Electromagnetic Fields (McGraw-Hill, 1961).

42. S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antenn. Propag. 30(3), 409–418 (1982). [CrossRef]

43. P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B 6(12), 4370–4379 (1972). [CrossRef]

44. T. W. Lloyd, J. M. Song, and M. Yang, “Numerical study of surface integral formulations for low-contrast objects,” IEEE Antennas Wirel. Propag. Lett. 4(1), 482–485 (2005). [CrossRef]

45. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, 1983).

Formulation	a_i for i = 1,2	b_i for i = 1,2	c_i for i = 1,2	d_i for i = 1,2	$O_{n o r m} (\bar{Z})$
PMCHWT	η_i	0	0	1/η_i	$(\begin{matrix} 1 & 1 / η \\ 1 / η & 1 / η^{2} \end{matrix})$
CTF	1	0	0	1	$(\begin{matrix} 1 & 1 / η \\ η & 1 \end{matrix})$
CNF	0	1	1	0	$(\begin{matrix} 1 & 1 / η \\ η & 1 \end{matrix})$
JMCFIE	1	1	1	1	$(\begin{matrix} 1 & 1 / η \\ η & 1 \end{matrix})$
Müller	0	μ_i	ε_i	0	$(\begin{matrix} 1 & 1 / η \\ 1 / η & 1 / η^{2} \end{matrix})$

Formulation	a_i for i = 1,2	b_i for i = 1,2	c_i for i = 1,2	d_i for i = 1,2	$O_{n o r m} (\bar{Z})$
PMCHWT	η_i	0	0	1/η_i	$(\begin{matrix} 1 & 1 / η \\ 1 / η & 1 / η^{2} \end{matrix})$
CTF	1	0	0	1	$(\begin{matrix} 1 & 1 / η \\ η & 1 \end{matrix})$
CNF	0	1	1	0	$(\begin{matrix} 1 & 1 / η \\ η & 1 \end{matrix})$
JMCFIE	1	1	1	1	$(\begin{matrix} 1 & 1 / η \\ η & 1 \end{matrix})$
Müller	0	μ_i	ε_i	0	$(\begin{matrix} 1 & 1 / η \\ 1 / η & 1 / η^{2} \end{matrix})$

Improving condition number and convergence of the surface integral-equation method of moments for penetrable bodies

Abstract

1. Introduction

2. Surface integral-equation formulations for composite penetrable bodies

3. Discretization of the integral equations

4. Numerical balance of the impedance matrix

5. Left and right preconditioner

6. Numerical examples

7. Conclusion

Acknowledgments

References and links

Cited By

Figures (4)

Tables (2)

Equations (44)

Optics Express

Formulation	a	b	c	d	$α_{11}$	$α_{22}$	$β_{11}$	$β_{22}$
PMCHWT	η	0	0	1/η	1	η	1	η
CTF	1	0	0	1	1	1/η	1	η
CNF	0	1	1	0	1	1/η	1	η
JMCFIE	1	1	1	1	1	1/η	1	η
Müller	0	μ	ε	0	1	η	1	η