Numerical Method for Solving Electromagnetic Wave Scattering by One and Many Small Perfectly Conducting Bodies

Nhan T. Tran

Department of Mathematics, Kansas State University, Manhattan, USA.

DOI: 10.4236/ajcm.2017.74030                

Abstract

In this paper, we investigate the problem of electromagnetic (EM) wave scattering by one and many small perfectly conducting bodies and present a numerical method for solving it. For the case of one body, the problem is solved for a body of arbitrary shape, using the corresponding boundary integral equation. For the case of many bodies, the problem is solved asymptotically under the physical assumptions a << d << λ, where a is the characteristic size of the bodies, d is the minimal distance between neighboring bodies, λ = 2π/k is the wave length and k is the wave number. Numerical results for the cases of one and many small bodies are presented. Error analysis for the numerical method is also provided.

Keywords

Electromagnetic Scattering, Many Bodies, Perfectly Conducting Body, Integral Equation, EM Waves

1. Introduction

Many real-world electromagnetic (EM) problems like EM wave scattering, EM radiation, etc. [1] , cannot be solved analytically and exactly to get a solution in a closed form. Thus, numerical methods have been developed to tackle these problems approximately. Computational Electromagnetics (CEM) has evolved enormously in the past decades to a point that its methods can solve EM problems with extreme accuracy. These methods can be classified into two categories: Integral Equation (IE) method and Differential Equation (DE) method. Typical IE methods include: Method of Moment (MoM) developed by Roger F. Harrington (1968) [2] , Fast Multipole Method (FMM) first introduced by Greengard and Rokhlin (1987) [3] and then applied to EM by Engheta et al. (1992) [4] , Partial Element Equivalent Circuit (PEEC) method [5] , and Discrete Dipole Approximation [6] . Typical DE methods are: Finite Difference Time Domain (FDTD) developed by Kane Yee (1966) [7] , Finite Element Method (FEM) [8] , Finite Integration Technique (FIT) proposed by Thomas Weiland (1977) [9] , Pseudospectral Time Domain (PSTD) [10] , Pseudospectral Spatial Domain (PSSD) [11] , and Transmission Line Matrix (TLM) [12] . Among these methods, FDTD has emerged as one of the most popular techniques for solving EM problems due to its simplicity and ability to provide animated display of the EM field, see [12] . However, FDTD requires the entire computational domain be gridded [13] , which results in very long solution times. Furthermore, as a DE method, it does not take into account the radiation condition in exact sense [14] [15] , which leads to certain error in the solution. On the other hand, spurious solutions might exist in DE methods [16] [17] [18] . Most importantly, most of DE methods are not suitable if the number of bodies is very large.

In [19] - [24] , A. G. Ramm has developed a theory of EM wave scattering by many small perfectly conducting and impedance bodies. In this theory, the EM wave scattering problem is solved asymptotically under the physical assumptions: $a\ll d\ll \lambda$ , where a is the characteristic size of the bodies, d is the minimal distance between neighboring bodies, $\lambda =2\text{π}/k$ is the wave length and k is the wave number. In [25] , a numerical method is developed for solving EM wave scattering by many small impedance bodies. In this paper, the problem of EM wave scattering by one and many small perfectly conducting bodies is considered. A numerical method for solving this problem asymptotically based on the above theory is presented. For the case of one body, the problem is solved for a body of arbitrary shape, using the corresponding boundary integral equation. For the case of many small bodies, the problem is solved under the basic assumptions $a\ll d\ll \lambda$ and the assumption about the distribution of the small bodies

$\mathcal{N}\left(\Delta \right)=\frac{1}{a^3}{\displaystyle {\int }_{\Delta }}N\left(x\right)\text{d}x\left[1+o\left(1\right)\right]\mathrm{,}a\to \mathrm{0,}$ (1)

in which $\Delta$ is an arbitrary open subset of the domain $\Omega$ that contains all the small bodies, $\mathcal{N}\left(\Delta \right)$ is the number of the small bodies in $\Delta$ , and $N\left(x\right)$ is the distribution function of the bodies

$N\left(x\right)\ge \mathrm{0,}N\left(x\right)\in C\left(\Omega \right)\mathrm{.}$ (2)

In Sections 2 and 3, the theory of EM wave scattering by one and many small perfectly conducting bodies is presented. The numerical methods for solving these problems are also described in details. Furthermore, error analysis for the numerical methods of solving the EM scattering problem are also provided. In Section 4 these methods are tested and numerical results are discussed.

2. EM Wave Scattering by One Perfectly Conducting Body

Let D be a bounded perfectly conducting body, $a=\frac{1}{2}\text{diam}D$ , S be its C²-smooth

boundary, and $D^{\prime }\mathrm{<mo>\; :\; </mo>}={ℝ}^3\backslash D$ . Let $ϵ$ and $\mu$ be the dielectric permittivity and magnetic permeability constants of the medium in $D^{\prime }$ . Let E and H denote the electric and magnetic fields, respectively, $E_0$ be the incident field and $v_E$ be the scattered field. The problem of electromagnetic wave scattering by one perfectly conducting body can be stated as follows

$\nabla \times E=i\omega \mu H\mathrm{,}\text{in}{D}^{\prime }\mathrm{<mo>\; :\; </mo>}={ℝ}^3\backslash D\mathrm{,}$ (3)

$\nabla \times H=-i\omega ϵE\mathrm{,}\text{in}{D}^{\prime }\mathrm{,}$ (4)

$\left[N\mathrm{,}\left[E\mathrm{,}N\right]\right]=0\mathrm{,}\text{on}S\mathrm{<mo>\; :\; </mo>}=\partial D\mathrm{,}$ (5)

$E=E_0+v_E,$ (6)

$E_0=\mathcal{E}{\text{e}}^{ik\alpha \cdot x},\mathcal{E}\cdot \alpha =0,\alpha \in S^2,$ (7)

$\frac{\partial v_E}{\partial r}-ik{v}_E=o\left(\frac{1}{r}\right)\mathrm{,}r\mathrm{<mo>\; :\; </mo>}=\left|x\right|\to \infty \mathrm{,}$ (8)

where $\omega >0$ is the frequency, $k=2\text{π}/\lambda =\omega \sqrt{ϵ\mu }$ is the wave number, $ka\ll 1$ , $\lambda$ is the wave length, $\mathcal{E}$ is a constant vector, and $\alpha$ is a unit vector that indicates the direction of the incident wave $E_0$ . This incident wave satisfies the relation $\nabla \cdot E_0=0$ . The scattered field $v_E$ satisfies the radiation Condition (8). Here, N is the unit normal vector to the surface S, pointing out of D. By $\left[\cdot \mathrm{,}\cdot \right]$ the vector product is denoted and $\alpha \cdot x$ is the scalar product of two vectors.

The solution to problem (3)-(8) can be found in the form

$E\left(x\right)=E_0\left(x\right)+\nabla \times {\displaystyle {\int }_S}g\left(x,t\right)J\left(t\right)\text{d}t,g\left(x,t\right):=\frac{{\text{e}}^{ik\left|x-t\right|}}{4\text{π}\left|x-t\right|},$ (9)

see [19] . Here, E is a vector in ${ℝ}^3$ and $\nabla \times E$ is a pseudo-vector, that is a vector-like object which changes sign under reflection of its coordinate axes. $E_0$ is the incident plane wave defined in (7) and J is an unknown pseudo-vector that is to be found. J is assumed to be tangential to S and continuous. J can be found by applying the boundary Condition (5), or equivalently $\left[N\mathrm{,}E\right]=0$ , to (9) and solving the resulting boundary integral equation

$\frac{J}{2}+AJ\mathrm{<mo>\; :\; </mo>}=\frac{J\left(s\right)}{2}+{\displaystyle {\int }_S}\left[N_s\mathrm{,}\left[{\nabla }_{s}g\left(s\mathrm{,}t\right)\mathrm{,}J\left(t\right)\right]\right]\text{d}t=-\left[N_s\mathrm{,}{E}_0\right]\mathrm{,}$ (10)

or, equivalently

$\left(I+2A\right)J=F\mathrm{,}$ (11)

where $F\mathrm{<mo>\; :\; </mo>}=-2\left[N_s\mathrm{,}{E}_0\right]$ . Equation (11) is of Fredholm type since A is compact, see [21] .

Once we have J, E can be computed by formula (9) and H can be found by the formula

$H=\frac{\nabla \times E}{i\omega \mu }.$ (12)

If D is sufficiently small, then Equation (11) is uniquely solvable in $C\left(S\right)$ and its solution J is tangential to S, see [19] . The asymptotic formula for E when the radius a of the body D tends to zero can be derived as follows, see [19] . Rewrite Equation (9) as

$E\left(x\right)=E_0\left(x\right)+\left[\nabla g\left(x,x_1\right),Q\right]+\nabla \times {\displaystyle {\int }_S}\left[g\left(x,t\right)-g\left(x,x_1\right)\right]J\left(t\right)\text{d}t,$ (13)

where $x_1\in D$ , an arbitrary point inside the small body D, and

$Q:={\displaystyle {\int }_S}J\left(t\right)\text{d}t.$ (14)

Since

$\left|\nabla g\left(x,x_1\right)\right|=O\left(\frac{k}{d}+\frac{1}{d^2}\right),d=\left|x-x_1\right|,$ (15)

$\left|g\left(x,t\right)-g\left(x,x_1\right)\right|=O\left(\left(\frac{k}{d}+\frac{1}{d^2}\right)a\right),a=\left|t-x_1\right|,$ (16)

and

$\left|\nabla \left[g\left(x,t\right)-g\left(x,x_1\right)\right]\right|=O\left(\frac{a{k}^2}{d}+\frac{ak}{d^2}+\frac{a}{d^3}\right),$ (17)

the second term in (13) is much greater than the last term

$\left|\left[\nabla g\left(x\mathrm{,}{x}_1\right)\mathrm{,}Q\right]\right|\gg \left|\nabla \times {\displaystyle {\int }_S}\left[g\left(x\mathrm{,}t\right)-g\left(x\mathrm{,}{x}_1\right)\right]J\left(t\right)\text{d}t\right|\mathrm{,}a\to 0.$ (18)

Then, the asymptotic formula for E when a tends to zero is

$E\left(x\right)=E_0\left(x\right)+\left[{\nabla }_{x}g\left(x\mathrm{,}{x}_1\right)\mathrm{,}Q\right]\mathrm{,}a\to \mathrm{0,}$ (19)

where $\left|x-x_1\right|\gg a$ , $x_1\in D$ . Thus, when D is sufficiently small, instead of finding J, we can just find one pseudovector Q.

The analytical formula for Q is derived as follows, see [21] . By integrating both sides of (10) over S, one gets

${\displaystyle {\int }_S}\frac{J\left(s\right)}{2}\text{d}s+{\displaystyle {\int }_S}\text{d}s{\displaystyle {\int }_S}\text{d}t\left[N_s\mathrm{,}\left[{\nabla }_{s}g\left(s\mathrm{,}t\right)\mathrm{,}J\left(t\right)\right]\right]=-{\displaystyle {\int }_S}\left[N_s\mathrm{,}{E}_0\right]\text{d}s\mathrm{.}$ (20)

This is equivalent to

$\frac{Q}{2}+{\displaystyle {\int }_S}\text{d}t{\displaystyle {\int }_S}\text{d}s{\nabla }_{s}g\left(s\mathrm{,}t\right)N_s\cdot J\left(t\right)-{\displaystyle {\int }_S}\text{d}tJ\left(t\right){\displaystyle {\int }_S}\text{d}s\frac{\partial g\left(s\mathrm{,}t\right)}{\partial N_s}=-{\displaystyle {\int }_D}\nabla \times E_0\text{d}x\mathrm{.}$ (21)

When $a\to 0$ , this equation becomes

$\frac{Q}{2}+e_p{\displaystyle {\int }_S}\text{d}t{\displaystyle {\int }_S}\text{d}s\frac{\partial g\left(s\mathrm{,}t\right)}{\partial s_p}{N}_q\left(s\right)J_q\left(t\right)+\frac{1}{2}{\displaystyle {\int }_S}\text{d}tJ\left(t\right)=-\left|D\right|\nabla \times E_0\mathrm{,}1\le p\mathrm{,}q\le \mathrm{3,}$ (22)

where in the second term, summations over the repeated indices are understood, $e_p,1\le p\le 3$ , are the orthogonal unit vectors in ${ℝ}^3$ , $\left|D\right|$ is the volume of D, $\left|D\right|=c_D{a}^3$ , and in the third term we use this estimate

${\displaystyle {\int }_S}\text{d}s\frac{\partial g\left(s\mathrm{,}t\right)}{\partial N_s}\simeq {\displaystyle {\int }_S}\text{d}s\frac{\partial g_0\left(s\mathrm{,}t\right)}{\partial N_s}=-\frac{1}{2}\mathrm{,}{g}_0\left(s\mathrm{,}t\right)\mathrm{<mo>\; :\; </mo>}=\frac{1}{4\text{π}\left|s-t\right|}\mathrm{.}$ (23)

Let

${\Gamma }_{pq}\left(t\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_S}\text{d}s\frac{\partial g\left(s\mathrm{,}t\right)}{\partial s_p}{N}_q\left(s\right)\mathrm{,}$ (24)

then Equation (22) can be rewritten as follows

$\frac{Q}{2}+e_p{\displaystyle {\int }_S}\text{d}t{\Gamma }_{pq}\left(t\right)J_q\left(t\right)+\frac{Q}{2}=-\left|D\right|\nabla \times E_0\mathrm{,}$ (25)

or

$Q+\Gamma Q=-\left|D\right|\nabla \times E_0,$ (26)

where $\Gamma$ is a $3\times 3$ constant matrix and it is defined by

$\Gamma Q=e_p{\displaystyle {\int }_S}\text{d}t{\Gamma }_{pq}\left(t\right)J_q\left(t\right),$ (27)

in which summations are understood over the repeated indices. Thus, Q can be written as

$Q=-\left|D\right|{\left(I+\Gamma \right)}^{-1}\nabla \times E_0,a\to 0,$ (28)

where $I:=I_3$ , the $3\times 3$ identity matrix. This formula is asymptotically exact as $a\to 0$ .

2.1. Numerical Method for Solving EM Wave Scattering by One Perfectly Conducting Spherical Body

In this section, we consider the EM wave scattering problem by a small perfectly conducting spherical body. Instead of solving the problem (3)-(8) directly, we will solve its corresponding boundary integral Equation (10) for the unknown vector J

$\frac{J\left(s\right)}{2}+{\displaystyle {\int }_S}\left[N_s\mathrm{,}\left[{\nabla }_{s}g\left(s\mathrm{,}t\right)\mathrm{,}J\left(t\right)\right]\right]\text{d}t=-\left[N_s\mathrm{,}{E}_0\right]\mathrm{.}$ (29)

Then the solution E to the EM wave scattering problem by one perfectly conducting body can be computed by either the exact formula (9) or the asymptotic formula (19).

Scattering by a sphere has been discussed in many papers, for example [26] in which Mie solves the EM wave scattering problem by separation of variables. The EM field, scattered by a small body, is proportional to $O\left(a^3\right)$ .

Suppose S is a smooth surface of a spherical body. Let S be partitioned into P non-intersecting subdomains $S_{ij}\mathrm{,1}\le i\le m_{\theta }\mathrm{,1}\le j\le m_{\varphi }$ , using spherical coordinates, where $m_{\theta }$ is the number of intervals of $\theta$ between 0 and 2π and $m_{\varphi }$ defines the number of intervals of $\varphi$ between 0 and π. Then $P=m_{\theta }{m}_{\varphi }+2$ , which includes the two poles of the sphere. $m_{\theta }$ is defined in this way:

$m_{\theta }\mathrm{<mo>\; =\; </mo>}{m}_{\varphi }+\left|\varphi -\frac{\text{π}}{2}\right|6m_{\varphi }$ . This means the closer it is to the poles of the sphere, the

more intervals for $\theta$ are used. Then the point $\left({\theta }_i\mathrm{,}{\varphi }_j\right)$ in $S_{ij}$ is chosen as follows

${\theta }_i=i\frac{2\text{π}}{m_{\theta }},1\le i\le m_{\theta },$ (30)

${\varphi }_j=j\frac{\text{π}}{m_{\varphi }+1},1\le j\le m_{\varphi }.$ (31)

Note that there are many different ways to distribute collocation points. However, the one that we describe here will guarantee convergence to the solution to (29) with fewer collocation points used from our experiment. Furthermore, one should be careful when choosing the distribution of

collocation points on a sphere. If one chooses ${\varphi }_j=j\frac{\text{π}}{m_{\varphi }},1\le j\le m_{\varphi }$ , then when

$j=m_{\varphi }$ , ${\varphi }_j=\text{π}$ and thus there is only one point for this $\varphi$ regardless of the value of $\theta$ as shown in (32). The position of a point in each $S_{ij}$ can be computed by

${\left(x\mathrm{,}y\mathrm{,}z\right)}_{ij}=a\left(cos{\theta }_{i}sin{\varphi }_j\mathrm{,}sin{\theta }_{i}sin{\varphi }_j\mathrm{,}cos{\varphi }_j\right)\mathrm{,}$ (32)

and the outward-pointing unit normal vector N to S at this point is

$N_{ij}=N\left({\theta }_i,{\varphi }_j\right)=\left(\cos {\theta }_i\sin {\varphi }_j,\sin {\theta }_i\sin {\varphi }_j,\cos {\varphi }_j\right).$ (33)

For a star-shaped body with a different shape, only the normal vector N needs to be recomputed. Rewrite the integral Equation (29) as

$\frac{J\left(s\right)}{2}+{\displaystyle {\int }_S}{\nabla }_{s}g\left(s\mathrm{,}t\right)N_s\cdot J\left(t\right)\text{d}t-{\displaystyle {\int }_S}\frac{\partial g\left(s\mathrm{,}t\right)}{\partial N_s}J\left(t\right)\text{d}t=-\left[N_s\mathrm{,}{E}_0\right]\mathrm{.}$ (34)

This integral equation can be discretized as follows

$J\left(i\right)+2{\displaystyle \underset{j\ne i}{\overset{P}{\sum }}}\left[{\nabla }_{s}g\left(i\mathrm{,}j\right)N_s\left(i\right)\cdot J\left(j\right)-J\left(j\right){\nabla }_{s}g\left(i\mathrm{,}j\right)\cdot N_s\left(i\right)\right]{\Delta }_j=F\left(i\right)\mathrm{,}1\le i\le P\mathrm{,}$

(35)

in which by i the point $\left(x_i\mathrm{,}{y}_i\mathrm{,}{z}_i\right)$ is denoted, $F\left(i\right)\mathrm{<mo>\; :\; </mo>}=-2\left[N_s\mathrm{,}{E}_0\right]\left(i\right)$ , and ${\Delta }_j$ is the surface area of the subdomain j. This is a linear system with unknowns $J\left(i\right):=\left(X_i,Y_i,Z_i\right),1\le i\le P$ . This linear system can be rewritten as follows

$X_i+{\displaystyle \underset{j\ne i}{\overset{P}{\sum }}}{a}_{ij}{X}_j+b_{ij}{Y}_j+c_{ij}{Z}_j=F_x\left(i\right),$ (36)

$Y_i+{\displaystyle \underset{j\ne i}{\overset{P}{\sum }}}{a^{\prime }}_{ij}{X}_j+{b^{\prime }}_{ij}{Y}_j+{c^{\prime }}_{ij}{Z}_j=F_y\left(i\right),$ (37)

$Z_i+{\displaystyle \underset{j\ne i}{\overset{P}{\sum }}}{a^{″}}_{ij}{X}_j+{b^{″}}_{ij}{Y}_j+{c^{″}}_{ij}{Z}_j=F_z\left(i\right),$ (38)

where by the subscripts $x\mathrm{,}y\mathrm{,}z$ the corresponding coordinates are denoted, e.g. $F\left(i\right)=\left(F_x\mathrm{,}{F}_y\mathrm{,}{F}_z\right)\left(i\right)$ , and

$a_{ij}:=2\left[\nabla g{\left(i,j\right)}_x{N}_x\left(i\right)-\nabla g\left(i,j\right)\cdot N\left(i\right)\right]{\Delta }_j,$ (39)

$b_{ij}:=2\nabla g{\left(i,j\right)}_x{N}_y\left(i\right){\Delta }_j,$ (40)

$c_{ij}:=2\nabla g{\left(i,j\right)}_x{N}_z\left(i\right){\Delta }_j,$ (41)

for $i\ne j$ ; when $i=j$ : $a_{ii}=1,b_{ii}=0$ , and $c_{ii}=0$ ,

${a^{\prime }}_{ij}:=2\nabla g{\left(i,j\right)}_y{N}_x\left(i\right){\Delta }_j,$ (42)

${b^{\prime }}_{ij}:=2\left[\nabla g{\left(i,j\right)}_y{N}_y\left(i\right)-\nabla g\left(i,j\right)\cdot N\left(i\right)\right]{\Delta }_j,$ (43)

${c^{\prime }}_{ij}:=2\nabla g{\left(i,j\right)}_y{N}_z\left(i\right){\Delta }_j,$ (44)

for $i\ne j$ ; when $i=j$ : ${a^{\prime }}_{ii}=0,{b^{\prime }}_{ii}=1$ , and ${c^{\prime }}_{ii}=0$ ,

${a^{″}}_{ij}:=2\nabla g{\left(i,j\right)}_z{N}_x\left(i\right){\Delta }_j,$ (45)

${b^{″}}_{ij}:=2\nabla g{\left(i,j\right)}_z{N}_y\left(i\right){\Delta }_j,$ (46)

${c^{″}}_{ij}:=2\left[\nabla g{\left(i,j\right)}_z{N}_z\left(i\right)-\nabla g\left(i,j\right)\cdot N\left(i\right)\right]{\Delta }_j,$ (47)

for $i\ne j$ ; when $i=j$ : ${a^{″}}_{ii}=0,{b^{″}}_{ii}=0$ , and ${c^{″}}_{ii}=1$ .

Solving system (36)-(38) yields $X_i\mathrm{,}{Y}_i\mathrm{,}{Z}_i$ , and $J\left(i\right)$ can be computed by $J\left(i\right)=\left(X_i,Y_i,Z_i\right),1\le i\le P$ .

2.2. Error Analysis

Recall the boundary integral Equation (10)

$\frac{J\left(s\right)}{2}+{\displaystyle {\int }_S}\left[N_s,\left[{\nabla }_{s}g\left(s,t\right),J\left(t\right)\right]\right]\text{d}t=-\left[N_s,E_0\right].$ (48)

Integrate both sides of this equation over S and get

$Q+\Gamma Q=-\left|D\right|\nabla \times E_0,$ (49)

see Section 2. Once J is found from solving (48), Q can be computed by $Q={\displaystyle {\int }_S}J\left(t\right)\text{d}t$ . Then one can validate the values of J and Q by checking the following things

・ Is J tangential to S as shown in Section 2? One needs to check $J\left(s\right)\cdot N_s$ .

・ Is $Q={\displaystyle {\int }_S}J\left(t\right)\text{d}t$ correct? The relative error of Q can be computed as follows

$\text{Error}=\frac{\left|Q+\Gamma Q-RHS\right|}{\left|RHS\right|}\mathrm{,}$ (50)

where $RHS\mathrm{<mo>\; :\; </mo>}=-\left|D\right|\nabla \times E_0$ . This will give the error of the numerical method for the case of one body.

Furthermore, one can also compare the value of the asymptotic $Q_a$ in formula (28) with the exact $Q_e$ defined in (14) by

$\text{Error}=\frac{\left|Q_e-Q_a\right|}{\left|Q_e\right|}\mathrm{,}$ (51)

and check the difference between the asymptotic $E_a$ in (19) and the exact $E_e$ defined in (9) by computing this relative error

$\text{Error}=\frac{\left|E_e-E_a\right|}{\left|E_e\right|}\mathrm{.}$ (52)

2.3. General Method for Solving EM Wave Scattering by One Perfectly Conducting Body

In this section, we present a general method for solving the EM wave scattering problem by one perfectly conducting body, whose surface is parametrized by $f\left(u,v\right)=\left(x\left(u,v\right),y\left(u,v\right),z\left(u,v\right)\right)$ .

・ Step 1: One needs to partition the surface of the body into P non-intersecting subdomains. In each subdomain, choose a collocation point. The position of the collocation points can be computed using $f\left(u,v\right)=\left(x\left(u,v\right),y\left(u,v\right),z\left(u,v\right)\right)$ , see for example (30)-(32).

・ Step 2: Find the unit normal vector N of the surface from the function f.

・ Step 3: Solve the linear system (36)-(38) for $X_i\mathrm{,}{Y}_i$ , and $Z_i\mathrm{,1}\le i\le P$ . Then vector J in the boundary integral Equation (10) is computed by $J\left(i\right):=\left(X_i,Y_i,Z_i\right)$ at the point i on the surface.

・ Step 4: Compute the electric field E using (9).

3. EM Wave Scattering by Many Small Perfectly Conducting Bodies

Consider a bounded domain $\Omega$ containing M small bodies $D_m$ , $1\le m\le M$ , and $S_m$ are their corresponding smooth boundaries. Let $D:={\displaystyle {\cup }_{m=1}^M}{D}_m\subset \Omega$ and $D^{\prime }$ be the complement of D in ${ℝ}^3$ . We assume that $S={\displaystyle {\cup }_{m=1}^M}{S}_m$ is C²-smooth. $ϵ$ is the dielectric permittivity constant and $\mu$ is the magnetic permeability constant of the medium. Let E and H denote the electric and magnetic fields, respectively. $E_0$ is the incident field and $v$ is the scattered field. The problem of electromagnetic wave scattering by many small perfectly conducting bodies involves solving the following system

$\nabla \times E=i\omega \mu H\mathrm{,}\text{in}{D}^{\prime }\mathrm{<mo>\; :\; </mo>}={ℝ}^3\backslash D\mathrm{,}D\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{m=1}{\overset{M}{\cup }}}{D}_m\mathrm{,}$ (53)

$\nabla \times H=-i\omega ϵE\mathrm{,}\text{in}{D}^{\prime }\mathrm{,}$ (54)

$\left[N\mathrm{,}\left[E\mathrm{,}N\right]\right]=0\mathrm{,}\text{on}S\mathrm{,}$ (55)

$E=E_0+v,$ (56)

$E_0=\mathcal{E}{\text{e}}^{ik\alpha \cdot x},\mathcal{E}\cdot \alpha =0,\alpha \in S^2.$ (57)

where $v$ satisfies the radiation Condition (8), $\omega >0$ is the frequency, $k=2\text{π}/\lambda$

is the wave number, $ka\ll 1$ , $a\mathrm{<mo>\; :\; </mo>}=\frac{1}{2}ma{x}_m\text{diam}{D}_m$ , and $\alpha$ is a unit vector that

indicates the direction of the incident wave $E_0$ . Furthermore,

$ϵ={ϵ}_0,\mu ={\mu }_0\text{in}{\Omega }^{\prime }:={ℝ}^3\backslash \Omega .$ (58)

Assume that the distribution of small bodies $D_m,1\le m\le M$ , in $\Omega$ satisfies the following formula

$\mathcal{N}\left(\Delta \right)=\frac{1}{a^3}{\displaystyle {\int }_{\Delta }}N\left(x\right)\text{d}x\left[1+o\left(1\right)\right]\mathrm{,}a\to \mathrm{0,}$ (59)

where $\mathcal{N}\left(\Delta \right)$ is the number of small bodies in $\Delta$ , $\Delta$ is an arbitrary open subset of $\Omega$ , and $N\left(x\right)$ is the distribution function

$N\left(x\right)\ge \mathrm{0,}N\left(x\right)\in C\left(\Omega \right)\mathrm{.}$ (60)

Note that E solves this equation

$\nabla \times \nabla \times E=k^{2}E,k^2={\omega }^2ϵ\mu ,$ (61)

if μ = const. Once we have E, then H can be found from this relation

$H=\frac{\nabla \times E}{i\omega \mu }.$ (62)

From (62) and (61), one can get (54). Thus, we need to find only E which satisfies the boundary condition (55). It was proved in [19] that under the radiation condition and the assumptions $a\ll d\ll \lambda$ , the problem (53)-(57) has a unique solution and its solution is of the form

$E\left(x\right)=E_0\left(x\right)+{\displaystyle \underset{m=1}{\overset{M}{\sum }}}\nabla \times {\displaystyle {\int }_{S_m}}g\left(x,t\right)J_m\left(t\right)\text{d}t,$ (63)

where $J_m$ are unknown continuous functions that can be found from the boundary condition. Let

$Q_m:={\displaystyle {\int }_{S_m}}{J}_m\left(t\right)\text{d}t.$ (64)

When $a\to 0$ , the asymptotic solution for the electric field is given by

$E\left(x\right)=E_0\left(x\right)+{\displaystyle \underset{m=1}{\overset{M}{\sum }}}\left[\nabla g\left(x,x_m\right),Q_m\right],a\to 0.$ (65)

Therefore, instead of finding $J_m\left(t\right)\mathrm{,}\forall t\in S\mathrm{,1}\le m\le M$ , to get the solution E, one can just find $Q_m$ . This allows one to solve the EM scattering problem with a very large number of small bodies which is impossible to do before. The analytic formula for $Q_m$ can be derived by using formula (28) and replacing $E_0$ in this formula by the effective field $E_e\left(x_m\right)$ acting on the m-th body

$Q_m=-\left|D_m\right|{\left(I+\Gamma \right)}^{-1}\nabla \times E_e\left(x_m\right)\mathrm{,}1\le m\le M\mathrm{,}{x}_m\in D_m\mathrm{,}$ (66)

where the effective field acting on the m-th body is defined as

$E_e\left(x_m\right)=E_0\left(x_m\right)+{\displaystyle \underset{j\ne m}{\overset{M}{\sum }}}\left[\nabla g\left(x_m,x_j\right),Q_j\right],1\le m\le M.$ (67)

When $a\to 0$ , the effective field $E_e\left(x\right)$ is asymptotically equal to the field $E\left(x\right)$ in (65) as proved in [19] and [21] .

Let $E_{em}:=E_e\left(x_m\right)$ , where $x_m$ is a point in $D_m$ . From (66), and (67), one gets

$E_{em}=E_{0m}-{\displaystyle \underset{j\ne m}{\overset{M}{\sum }}}\left[\nabla g\left(x_m,x_j\right),{\left(I+\Gamma \right)}^{-1}\nabla \times E_{ej}\right]\left|D_j\right|,1\le m\le M.$ (68)

3.1. Numerical Method for Solving EM Wave Scattering by Many Small Perfectly Conducting Bodies

For finding the solution to EM wave scattering in the case of many small perfectly conducting bodies, we need to find $E_{em}$ in (68). Apply the operator ${\left(I+\Gamma \right)}^{-1}\nabla \times$ to both sides of (68) and let $A_m\mathrm{<mo>\; :\; </mo>}={\left(I+\Gamma \right)}^{-1}\nabla \times E_{em}$ . Then

$A_m=A_{0m}-{\left(I+\Gamma \right)}^{-1}{\displaystyle \underset{j\ne m}{\overset{M}{\sum }}}\left|D_j\right|{\left({\nabla }_x\times \left[\nabla g\left(x,x_j\right),A_j\right]\right)|}_{x=x_m},1\le m\le M,$ (69)

Solving this system yields $A_m$ , for $1\le m\le M$ . Then E can be computed by

$E\left(x\right)=E_0\left(x\right)+{\displaystyle \underset{m=1}{\overset{M}{\sum }}}\left[\nabla g\left(x,x_m\right),Q_m\right],$ (70)

where

$Q_m=-\left|D_m\right|A_m,1\le m\le M.$ (71)

Equation (69) can be rewritten as follows

$A_m=A_{0m}-{\displaystyle \underset{j\ne m}{\overset{M}{\sum }}}\tau \left[k^{2}g\left(x_m,x_j\right)A_j+\left(A_j\cdot {\nabla }_x\right){\nabla g\left(x,x_j\right)|}_{x=x_m}\right]\left|D_j\right|,$ (72)

where $1\le m\le M$ , $\tau \mathrm{<mo>\; :\; </mo>}={\left(I+\Gamma \right)}^{-1}$ , and $A_m$ are vectors in ${ℝ}^3$ .

Let $A_i:=\left(X_i,Y_i,Z_i\right)$ then one can rewrite the system (72) as

$X_i+{\displaystyle \underset{j\ne i}{\overset{M}{\sum }}}{a}_{ij}{X}_j+b_{ij}{Y}_j+c_{ij}{Z}_j=F_x\left(i\right),$ (73)

$Y_i+{\displaystyle \underset{j\ne i}{\overset{M}{\sum }}}{a^{\prime }}_{ij}{X}_j+{b^{\prime }}_{ij}{Y}_j+{c^{\prime }}_{ij}{Z}_j=F_y\left(i\right),$ (74)

$Z_i+{\displaystyle \underset{j\ne i}{\overset{M}{\sum }}}{a^{″}}_{ij}{X}_j+{b^{″}}_{ij}{Y}_j+{c^{″}}_{ij}{Z}_j=F_z\left(i\right),$ (75)

in which by the subscripts $x\mathrm{,}y\mathrm{,}z$ the corresponding coordinates are denoted, e.g. $F\left(i\right)=\left(F_x\mathrm{,}{F}_y\mathrm{,}{F}_z\right)\left(i\right)$ , where $F\left(i\right):=A_{0i}$ and

$a_{ij}\mathrm{<mo>\; :\; </mo>}=\left[k^{2}g\left(i\mathrm{,}j\right)+{\partial }_x\nabla g{\left(i\mathrm{,}j\right)}_x\right]\left|D_j\right|\tau \left(\mathrm{1,1}\right)\mathrm{,}$ (76)

$b_{ij}\mathrm{<mo>\; :\; </mo>}={\partial }_y\nabla g{\left(i\mathrm{,}j\right)}_x\left|D_j\right|\tau \left(\mathrm{1,1}\right)\mathrm{,}$ (77)

$c_{ij}\mathrm{<mo>\; :\; </mo>}={\partial }_z\nabla g{\left(i\mathrm{,}j\right)}_x\left|D_j\right|\tau \left(\mathrm{1,1}\right)\mathrm{,}$ (78)

for $i\ne j$ , here $\tau \left(\mathrm{1,1}\right)$ is the entry (1,1) of matrix $\tau$ in (72); when $i=j$ : $a_{ii}=1,b_{ii}=0$ , and $c_{ii}=0$ ;

${a^{\prime }}_{ij}:={\partial }_x\nabla g{\left(i,j\right)}_y\left|D_j\right|\tau \left(2,2\right),$ (79)

${b^{\prime }}_{ij}:=\left[k^{2}g\left(i,j\right)+{\partial }_y\nabla g{\left(i,j\right)}_y\right]\left|D_j\right|\tau \left(2,2\right),$ (80)

${c^{\prime }}_{ij}:={\partial }_z\nabla g{\left(i,j\right)}_y\left|D_j\right|\tau \left(2,2\right),$ (81)

for $i\ne j$ ; when $i=j$ : ${a^{\prime }}_{ii}=0,{b^{\prime }}_{ii}=1$ , and ${c^{\prime }}_{ii}=0$ ;

${a^{″}}_{ij}:={\partial }_x\nabla g{\left(i,j\right)}_z\left|D_j\right|\tau \left(3,3\right),$ (82)

${b^{″}}_{ij}:={\partial }_y\nabla g{\left(i,j\right)}_z\left|D_j\right|\tau \left(3,3\right),$ (83)

${c^{″}}_{ij}:=\left[k^{2}g\left(i,j\right)+{\partial }_z\nabla g{\left(i,j\right)}_z\right]\left|D_j\right|\tau \left(3,3\right),$ (84)

for $i\ne j$ ; when $i=j$ : ${a^{\prime }}_{ii}=0,{b^{\prime }}_{ii}=0$ , and ${c^{\prime }}_{ii}=1$ .

3.2. Error Analysis

The error of the solution to the EM wave scattering problem by many small perfectly conducting bodies can be estimated as follows. From the solution E of the electromagnetic scattering problem by many small bodies given in (63)

$E\left(x\right)=E_0\left(x\right)+{\displaystyle \underset{m=1}{\overset{M}{\sum }}}\nabla \times {\displaystyle {\int }_{S_m}}g\left(x,t\right)J_m\left(t\right)\text{d}t,$ (85)

we can rewrite it as

$E\left(x\right)=E_0\left(x\right)+{\displaystyle \underset{m=1}{\overset{M}{\sum }}}\left[\nabla g\left(x,x_m\right),Q_m\right]+{\displaystyle \underset{m=1}{\overset{M}{\sum }}}\nabla \times {\displaystyle {\int }_{S_m}}\left[g\left(x,t\right)-g\left(x,x_m\right)\right]J_m\left(t\right)\text{d}t.$ (86)

Comparing this with the asymptotic formula for E when $a\to 0$ given in (65)

$E\left(x\right)=E_0\left(x\right)+{\displaystyle \underset{m=1}{\overset{M}{\sum }}}\left[\nabla g\left(x,x_m\right),Q_m\right],$ (87)

we have the error of this asymptotic formula is

$\text{Error}=\left|{\displaystyle \underset{m=1}{\overset{M}{\sum }}}\nabla \times {\displaystyle {\int }_{S_m}}\left[g\left(x\mathrm{,}t\right)-g\left(x\mathrm{,}{x}_m\right)\right]J_m\left(t\right)\text{d}t\right|$ (88)

$~\frac{1}{\text{4π}}\left(\frac{a{k}^2}{d}+\frac{ak}{d^2}+\frac{a}{d^3}\right){\displaystyle \underset{m=1}{\overset{M}{\sum }}}\left|Q_m\right|\mathrm{,}$ (89)

where $d={\min }_m\left|x-x_m\right|$ and

$Q_m=-\left|D_m\right|{\left(I+\Gamma \right)}^{-1}\nabla \times E_e\left(x_m\right),1\le m\le M,x_m\in D_m,a\to 0,$ (90)

because

$\left|\nabla \left[g\left(x\mathrm{,}t\right)-g\left(x\mathrm{,}{x}_m\right)\right]\right|=O\left(\frac{a{k}^2}{d}+\frac{ak}{d^2}+\frac{a}{d^3}\right)\mathrm{,}a=\underset{m}{max}\left|t-x_m\right|\mathrm{.}$ (91)

4. Experiments

4.1. EM Wave Scattering by One Perfectly Conducting Spherical Body

To illustrate the idea of the numerical method, we use the following physical parameters to solve the EM wave scattering problem by one small perfectly conducting sphere, i.e. solving the linear system (36)-(38)

・ Speed of wave, $c=\left(3.0E+10\right)\text{cm}/\text{sec}$ .

・ Frequency, $\omega =\left(5.0E+14\right)\text{Hz}$ .

・ Wave number, $k=\left(1.05E+05\right){\text{cm}}^{-1}$ .

・ Wave length, $\lambda =\left(6.00E-05\right)\text{cm}$ .

・ Direction of incident plane wave, $\alpha =\left(0,1,0\right)$ .

・ Magnetic permeability, $\mu =1$ .

・ Vector $\mathcal{E}=\left(\mathrm{1,0,0}\right)$ .

・ Incident field vector, $E_0$ : $E_0\left(x\right)=\mathcal{E}{\text{e}}^{ik\alpha \cdot x}$ .

・ The body is a sphere of radius a, centered at the origin.

We use GMRES iterative method, see [27] , to solve the linear system (36)-(38). For a spherical body, matrix $\Gamma$ in (26) can be computed analytically as follows. Recall that

${\Gamma }_{pq}\left(t\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_S}\frac{\partial g\left(s\mathrm{,}t\right)}{\partial s_p}{N}_q\left(s\right)\text{d}s\mathrm{,}1\le p\mathrm{,}q\le \mathrm{3,}$ (92)

where

$N=\left(cos\theta sin\varphi \mathrm{,}sin\theta sin\varphi \mathrm{,}cos\varphi \right)$ (93)

and

$\frac{\partial g\left(s\mathrm{,}t\right)}{\partial s_p}\simeq \frac{\partial g_0\left(s\mathrm{,}t\right)}{\partial s_p}=-\frac{s_p-t_p}{4\text{π}{\left|s-t\right|}^3}\mathrm{,}{g}_0\left(s\mathrm{,}t\right)\mathrm{<mo>\; :\; </mo>}=\frac{1}{4\text{π}\left|s-t\right|}\mathrm{.}$ (94)

We choose a coordinate system centered at the center of the sphere such that $t=\left(0,0,a\right)$ and $s=aN$ . Then

${\Gamma }_{pq}\left(t\right):=-\frac{a^2}{4\text{π}}{\displaystyle {\int }_0^{\text{2π}}}\text{d}\theta {\displaystyle {\int }_0^{\text{π}}}\text{d}\varphi \sin \varphi \frac{\left(s_p-t_p\right)N_q}{a^{3}8{\sin }^3\frac{\varphi }{2}},1\le p,q\le 3.$ (95)

When

$p=q=1:{\Gamma }_{11}\left(t\right)=-1/3,$ (96)

$p=q=2:{\Gamma }_{22}\left(t\right)=-1/3,$ (97)

$p=q=3:{\Gamma }_{33}\left(t\right)=1/6,$ (98)

$p\ne q:{\Gamma }_{pq}\left(t\right)=0.$ (99)

Therefore, matrix $\Gamma$ is

$\Gamma \simeq \left[\begin{array}{ccc}-1/3& 0& 0\\ 0& -1/3& 0\\ 0& 0& 1/6\end{array}\right]$ (100)

For example, Table 1 shows the exact and asymptotic vector Q when the radius of the body is $a=\left(1.0E-09\right)\text{cm}$ and the number of collocation points used to solve the integral Equation (29) is $P=766$ . Note that $a=\left(1.0E-09\right)\text{cm}$ satisfies $ka\ll 1$ . The point $x_1$ in (19) is taken at the center of the body, the origin. Table 2 and Table 3 show the exact and asymptotic vector $E=\left(E_x,E_y,E_z\right)$ , the electric field, at the point x outside of the body, respectively. The distance $\left|x-x_1\right|$ is measured in cm in these tables.

Table 1. Vector $Q_e$ and $Q_a$ when $P=766$ collocation points and $a=\left(1.0E-09\right)\text{cm}$ .

Table 2. Vector $E_e$ for one perfectly conducting body with $a=\left(1.0E-09\right)\text{cm}$ and $P=766$ collocation points.

Table 3. Vector $E_a$ for one perfectly conducting body with $a=\left(1.0E-09\right)\text{cm}$ and $P=766$ collocation points.

In this case, we also verify the following things:

1) Is J tangential to S?

In fact, this vector J is tangential to the surface S of the body, $J\cdot N_s=O\left({10}^{-14}\right)$ .

2) How accurate is the asymptotic formula (28) for Q?

We check the accuracy of the asymptotic formula for Q in (28) by comparing it with the exact formula (14), see Section 2.2, and the relative error is $4.21E-02$ . The more collocation points used, the little this relative error is.

3) How accurate is the asymptotic formula (19) for E?

The accuracy of the asymptotic formula for E in (19) can be checked by comparing it with the exact formula (9) at several points x outside of the body, $\left|x-x_1\right|\gg a$ where $x_1$ is the center of the body, see the error analysis in Section 2.2. The relative errors are given in Table 4.

Table 5 compares the asymptotic $Q_a$ versus exact $Q_e$ and asymptotic $E_a$ versus exact $E_e$ , when $P=1386$ collocation points, $\left|x-x_1\right|=1.73E-05$ cm, and with various a. The errors shown in this table are relative errors, see the error analysis in Section 2.2. As one can see from this table, the smaller the radius a is, compared to the distance from the point of interest to the center of the body, the more precise the asymptotic formulas of E and Q are.

Furthermore, the numerical results also depend on the number of collocation points used. The more collocation points used, the more accurate the results is.

4.2. EM Wave Scattering by One Perfectly Conducting Ellipsoid Body

In this section, we consider the EM wave scattering problem by a small perfectly conducting ellipsoid body. The method for solving the problem in this setting is the same as that of Section 2.1 except that one needs to recompute the unit normal vector N.

To get the solution of this problem, one can follow the steps in Section 2.1. In particular, one needs to solve the linear system (36)-(38).

To illustrate the idea, we use the same physical parameters as described in Section 4.1, except that the body now is an ellipsoid. Let S be its smooth surface. The way we partition S into many subdomains $S_{ij}$ is the same as the way we partition a spherical body as described in Section 2.1. Then, the position of the collocation point in each subdomain $S_{ij}$ is defined by

${\left(x\mathrm{,}y\mathrm{,}z\right)}_{ij}=\left(acos{\theta }_{i}sin{\varphi }_j\mathrm{,}bsin{\theta }_{i}sin{\varphi }_j\mathrm{,}ccos{\varphi }_j\right)\mathrm{,}$ (101)

Table 4. Relative errors between the asymptotic and exact formulas for E when $P=766$ collocation points and $a=\left(1.0E-09\right)\text{cm}$ .

Table 5. Relative errors of the asymptotic E and Q when $P=1386$ collocation points.

where a, b and c are the lengths of the semi-principal axes of the ellipsoid. The outward-pointing normal vector n to S at this point is

$n_{ij}=n\left({\theta }_i,{\varphi }_j\right)=2\left(\frac{\cos {\theta }_i\sin {\varphi }_j}{a},\frac{\sin {\theta }_i\sin {\varphi }_j}{b},\frac{\cos {\varphi }_j}{c}\right),$ (102)

and the corresponding unit normal vector N is

$N_{ij}=N\left({\theta }_i,{\varphi }_j\right)=n_{ij}/\left|n_{ij}\right|.$ (103)

For example, Table 6 and Table 7 show the exact and asymptotic vector $E=\left(E_x,E_y,E_z\right)$ , the electric field, got from solving this EM wave scattering problem with one perfectly conducting ellipsoid body, when the semi-principle axes of the body are $a=\left(1.0E-08\right)\text{cm}$ , $b=\left(1.0E-09\right)\text{cm}$ , $c=\left(1.0E-09\right)\text{cm}$ , and the number of collocation points is $P=1052$ . Note that a, b, and c satisfy $kmax\left(a\mathrm{,}b\mathrm{,}c\right)\ll 1$ . The point $x_1$ in (19) is taken at the center of the ellipsoid body. Each row in Table 6 and Table 7 shows the exact and asymptotic $E=\left(E_x\mathrm{,}{E}_y\mathrm{,}{E}_z\right)$ , respectively, at the point x outside of the body. The distance $\left|x-x_1\right|$ is measured in cm in these tables.

As for the case of one body, we need to verify the following things:

a) Is J tangential to S?

In fact, this vector J is tangential to the surface S of the body, $J\cdot N_s=O\left({10}^{-13}\right)$ .

b) Are Q and J correct? We check the relative error described in Section 2.2,

$\text{Error}=\frac{\left|Q+\Gamma Q-RHS\right|}{\left|RHS\right|}=\mathrm{14\; <mi>\; \%\; </mi>}$ . The more collocation points used, the smaller

this error is, for example, with $P=1762$ collocation points, this error is only 3.6%.

c) How accurate is the asymptotic formula (19) for E?

The accuracy of the asymptotic formula for E in (19) can be checked by comparing it with the exact formula (9) at several points x outside of the body, $\left|x-x_1\right|\gg max\left(a\mathrm{,}b\mathrm{,}c\right)$ where $x_1$ is the center of the body. The relative errors

Table 6. Vector $E_e$ for one perfectly conducting ellipsoid body with $a=\left(1.0E-08\right)\text{cm}$ , $b=\left(1.0E-09\right)\text{cm}$ , $c=\left(1.0E-09\right)\text{cm}$ , and $P=1052$ collocation points.

Table 7. Vector $E_a$ for one perfectly conducting ellipsoid body with $a=\left(1.0E-08\right)\text{cm}$ , $b=\left(1.0E-09\right)\text{cm}$ , $c=\left(1.0E-09\right)\text{cm}$ , and $P=1052$ collocation points.

are given in Table 8.

Table 9 shows the relative errors between the asymptotic E versus exact E, when $P=1052$ collocation points, $\left|x-x_1\right|=1.73E-07$ cm, and with various semi-principle axes a, b, and c. As one can see from this table, the smaller the semi-principle axes are, compared to the distance from the point of interest to the center of the body, the more accurate the asymptotic formulas of E is.

4.3. EM Wave Scattering by One Perfectly Conducting Cubic Body

In this section, we consider the EM wave scattering problem by a small perfectly conducting cubic body. Again, the method for solving the problem in this setting is the same as that of Section 2.1 except that one needs to recompute the unit normal vector N.

One can follow the steps outlined in Section 2.1 to solve this problem. That means, one needs to solve the linear system (36)-(38).

For illustration purpose, we use the same physical parameters as described in Section 4.1, except that the body now is a cube. Suppose the cube is placed in the first octant where the origin is one of its vertices, one can use the standard unit vectors in ${ℝ}^3$ as the unit normal vectors to the surfaces of the cube.

For example, Table 10 and Table 11 show the exact and asymptotic vector $E=\left(E_x,E_y,E_z\right)$ , the electric field, got from solving this EM wave scattering problem with one perfectly conducting cubic body, when the half side of the body is $a=\left(1.0E-07\right)\text{cm}$ and the number of collocation points is $M=600$ . Note that $a=1.0E-07\text{cm}$ satisfies $ka\ll 1$ . The point $x_1$ in (19) is taken at the center of the cubic body. Each row in Table 10 and Table 11 shows the exact and asymptotic $E=\left(E_x,E_y,E_z\right)$ , respectively, at the point $x$ outside of the body. The distance $\left|x-x_1\right|$ is measured in cm in these tables.

As before, for the case of one body, we need to verify the following things:

a) Is J tangential to S?

Table 8. Relative errors between the asymptotic and exact formulas for E when $P=1052$ collocation points, $a=\left(1.0E-08\right)\text{cm}$ , $b=\left(1.0E-09\right)\text{cm}$ , and $c=\left(1.0E-09\right)\text{cm}$ .

Table 9. Relative errors of the asymptotic E when $P=1052$ collocation points.

Table 10. Vector $E_e$ for one perfectly conducting body with $a=\left(1.0E-07\right)\text{cm}$ and $M=600$ collocation points.

Table 11. Vector $E_a$ for one perfectly conducting body with $a=\left(1.0E-07\right)\text{cm}$ and $M=600$ collocation points.

In fact, this vector J is tangential to the surface S of the body, $J\cdot N_s=O\left({10}^{-13}\right)$ .

b) How accurate is the asymptotic formula (28) for Q?

We check the relative error described in Section 2.2,

$\text{Error}=\frac{\left|Q+\Gamma Q-RHS\right|}{\left|RHS\right|}=\mathrm{1.13\; <mi>\; \%\; </mi>}$ .

c) How accurate is the asymptotic formula (19) for E?

The accuracy of the asymptotic formula for E in (19) can be checked by comparing it with the exact formula (9) at several points x outside of the body, $\left|x-x_1\right|\gg a$ where $x_1$ is the center of the body. The relative errors are given in Table 12.

Table 13 shows the relative errors between the asymptotic E versus exact E, when $M=600$ collocation points, $\left|x-x_1\right|=\left(1.73E-06\right)\text{cm}$ , and with various

Table 12. Relative errors between the asymptotic and exact formulas for E when $M=600$ collocation points and $a=\left(1.0E-07\right)\text{cm}$ .

Table 13. Relative errors of the asymptotic E when $M=600$ collocation points.

a. From this table, we can see that the smaller the side of the cube is, compared to the distance from the point of interest to the center of the body, the more accurate the asymptotic formula of E is.

4.4. EM Wave Scattering by Many Small Perfectly Conducting Bodies

To illustrate the idea, consider a domain $\Omega$ as a unit cube placed in the first octant such that the origin is one of its vertices. This domain $\Omega$ contains M small bodies. Suppose these small bodies are particles. We use GMRES iterative method, see [27] , to solve the linear system (73)-(75). The following physical parameters are used to solve the EM wave scattering problem

・ Speed of wave, $c=\left(3.0E+10\right)\text{cm}/\text{sec}$ .

・ Frequency, $\omega =\left(5.0E+14\right)\text{Hz}$ .

・ Wave number, $k=\left(1.05E+05\right){\text{cm}}^{-1}$ .

・ Wave length, $\lambda =\left(6.00E-05\right)\text{cm}$ .

・ Direction of incident plane wave, $\alpha =\left(\mathrm{0,1,0}\right)$ .

・ Magnetic permeability, $\mu =1$ .

・ Volume of the domain $\Omega$ that contains all the particles, $\left|\Omega \right|=1$ cm³.

・ The distance between two neighboring particles, $d=\left(1.00E-07\right)\text{cm}$ .

・ Vector $\mathcal{E}=\left(\mathrm{1,0,0}\right)$ .

・ Vector $A_0$ : $A_{0m}\mathrm{<mo>\; :\; </mo>}={\left(I+\Gamma \right)}^{-1}\nabla \times {E_0\left(x\right)|}_{x=x_m}={\left(I+\Gamma \right)}^{-1}\nabla \times {\mathcal{E}{\text{e}}^{ik\alpha \cdot x}|}_{x=x_m}$ .

Note that the distance d satisfies the assumption $d\ll \lambda$ . The radius a of the particles is chosen variously so that it satisfies the assumption $ka\ll 1$ . For illustration purpose, the problem of EM wave scattering by many small perfectly conducting bodies is solved with $M=27$ and 1000 particles.

For example, Table 14 show the result of solving the EM wave scattering problem with $M=27$ particles in the unit cube in which the distance between neighboring particles is $d=\left(1.0E-07\right)\text{cm}$ and the radius of the particles is $a=\left(1.0E-09\right)\text{cm}$ . Each row in Table 14 is a vector $E\left(i\right)=\left(E_x\mathrm{,}{E}_y\mathrm{,}{E}_z\right)\left(i\right)$ at

Table 14. Vector E when M = 27 particles, $d=\left(1.0E-07\right)\text{cm}$ and $a=\left(1.0E-09\right)\text{cm}$ .

the point i in the cube. The norm of this asymptotic solution E is $5.20E+00$ and the error of the solution is $8.16E-10$ . This error is computed using (89).

Table 15 shows the relative errors of E when there are $M=27$ particles in the cube, the distance between neighboring particles is $d=\left(1.0E-07\right)\text{cm}$ , and with various radius a. Figure 1 shows the relative error of the asymptotic E. From this figure, one can see that when the ratio a/d decreases from $1.0E-01$ to $1.0E-04$ , the error of the asymptotic solution decreases linearly and rapidly from $8.16E-06$ to about $8.16E-18$ . The smaller the ratio a/d is, the better the asymptotic formula (65) approximates E.

The slope of the curve in Figure 1 indicates the rate of decay of the error when the ratio a/d decreases.

Table 16 and Figure 2 show the results of solving the problem with

Table 15. Error of the asymptotic solution E when $M=27$ and $d=\left(1.0E-07\right)\text{cm}$ .

Table 16. Error of the asymptotic solution E when $M=1000$ and $d=\left(1.0E-07\right)\text{cm}$ .

Figure 1. Error of the asymptotic solution E when $M=27$ and $d=\left(1.0E-07\right)\text{cm}$ .

Figure 2. Error of the asymptotic solution E when $M=1000$ and $d=\left(1.0E-07\right)\text{cm}$ .

$M=1000$ particles, when the distance between neighboring particles is $d=\left(1.0E-07\right)\text{cm}$ , and with different radius a. From these table and figure, one can see that the relative error of the asymptotic solution in this case is also very small, less than $3.02E-04$ , when the ratio $a/d<1.0E-01$ . In this case, the error of the asymptotic E is greater than that of the previous case when $M=27$ . However, this time, the error is also decreasing quickly and linearly when the ratio a/d decreases from $1.0E-01$ to $1.0E-04$ . Therefore, the asymptotic formula (65) for the solution E is applicable when $a\ll d$ .

5. Conclusions

In this paper, we present a numerical method for solving the EM wave scattering by one and many small perfectly conducting bodies. One of the advantages of this method is that it is relatively easy to implement. Furthermore, one can get an asymptotically exact solution to the problem when the characteristic size of the bodies tends to zero. To illustrate the applicability and efficiency of the method, we use it to solve the EM wave scattering problem by one and many small perfectly conducting bodies. Numerical results of these experiments are presented and error analysis of the asymptotic solutions for the case of one and many bodies are also discussed. For the case of one small body, one can always find the exact solution using the described method. For the case of many small bodies, the accuracy of our method is high if $a\ll d\ll \lambda$ .

The problem of EM wave scattering is much harder to treat, compared to scalar wave scattering [28] [29] [30] [31] . For scalar wave scattering problem, a fast algorithm is developed in [31] to deal with billions of particles. This is still open to EM wave scattering. One might consider this as future research.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Tsang, L., Kong, J.A. and Ding, K.H. (2004) Scattering of Electromagnetic Waves, Theories and Applications. John Wiley & Sons, Vol. 27.
[2]	Harrington, R.F. and Harrington, J.L. (1996) Field Computation by Moment Methods. Oxford University Press.
[3]	Greengard, L. and Rokhlin, V. (1987) A Fast Algorithm for Particle Simulations. Journal of Computational Physics, 73, 325-348.
[4]	Engheta, N., Murphy, W.D., Rokhlin, V. and Vassiliou, M. (1985) The Fast Multipole Method for Electromagnetic Scattering Computation. IEEE Transactions on Antennas and Propagation, 40, 634-641.
[5]	Ruehli, A.E. (1974) Equivalent Circuit Models for Three-Dimensional Multiconductor Systems. IEEE Transactions on Microwave Theory and Techniques, 22, 216-221. https://doi.org/10.1109/TMTT.1974.1128204
[6]	DeVoe, H. (1964) Optical Properties of Molecular Aggregates. I. Classical Model of Electronic Absorption and Refraction. The Journal of Chemical Physics, 41, 393-400. https://doi.org/10.1063/1.1725879
[7]	Yee, K.S., et al. (1966) Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media. IEEE Transactions on Antennas and Propagation, 14, 302-307. https://doi.org/10.1109/TAP.1966.1138693
[8]	Zienkiewicz, O.C., Taylor, R.L., Zienkiewicz, O.C. and Taylor, R.L. (1977) The Finite Element Method. McGraw-Hill, London, 3.
[9]	Weiland, M.C.T. (2001) Discrete Electromagnetism with the Finite Integration Technique. Progress in Electromagnetics Research, 32, 65-87.
[10]	Liu, Q.H. (1997) The PSTD Algorithm: A Time-Domain Method Requiring Only Two Cells per Wavelength. Microwave and Optical Technology Letters, 15, 158-165. https://doi.org/10.1002/(SICI)1098-2760(19970620)15:3<158::AID-MOP11>3.0.CO;2-3
[11]	Tyrrell, J.C.A., Kinsler, P. and New, G.H.C. (2005) Pseudospectral Spatial-Domain: A New Method for Nonlinear Pulse Propagation in the Few-Cycle Regime with Arbitrary Dispersion. Journal of Modern Optics, 52, 973-986. https://doi.org/10.1080/09500340512331334086
[12]	Hoefer, W.J.R. (1985) The Transmission-Line Matrix Method-Theory and Applications. IEEE Transactions on Microwave Theory and Techniques, 33, 882-893.
[13]	Stavroulakis, P. (2013) Biological Effects of Electromagnetic Fields: Mechanisms, Modeling, Biological Effects, Therapeutic Effects, International Standards, Exposure Criteria. Springer Science & Business Media.
[14]	Kunz, K.S. and Luebbers, R.J. (1993) The Finite Difference Time Domain Method for Electromagnetics. CRC Press.
[15]	Buchanan, W.J., Gupta, N.K. and Arnold, J.M. (1993) Simulation of Radiation from a Microstrip Antenna using Three-Dimensional Finite-Difference Time-Domain (FDTD) Method, in Antennas and Propagation, 1993. 8th International Conference on IET, 639-642.
[16]	Mur, G. and Lager, I.E. (2002) On the Causes of Spurious Solutions in Electromagnetics. Electromagnetics, 22, 357-367. https://doi.org/10.1080/02726340290083950
[17]	Lu, M., Shanker, B. and Michielssen, E. (2004) Elimination of Spurious Solutions Associated with Exact Transparent Boundary Conditions in FDTD Solvers. IEEE Antennas and Wireless Propagation Letters, 3, 59-62.
[18]	Zhao, S. (2007) On the Spurious Solutions in the High-Order Finite Difference Methods for Eigenvalue Problems. Computer Methods in Applied Mechanics and Engineering, 196, 5031-5046.
[19]	Ramm, A.G. (2012) Electromagnetic Wave Scattering by Many Small Perfectly Conducting Particles of an Arbitrary Shape. Optics Communications, 285, 3679-3683.
[20]	Ramm, A. (2013) Scattering of Acoustic and Electromagnetic Waves by Small Bodies of Arbitrary Shapes. Applications to Creating New Engineered Materials, Momentum, New York. https://doi.org/10.5643/9781606506226
[21]	Ramm, A.G. (2015) Scattering of Electromagnetic Waves by Many Small Perfectly Conducting or Impedance Bodies. Journal of Mathematical Physics, 56, Article ID: 091901. https://doi.org/10.1063/1.4929965
[22]	Andriychuk, M.I., Indratno, S.W. and Ramm, A.G. (2012) Electromagnetic Wave Scattering by a Small Impedance Particle: Theory and Modeling. Optics Communications, 285, 1684-1691.
[23]	Ramm, A.G. and Andriychuk, M.I. (2014) Application of the Asymptotic Solution to EM Field Scattering Problem for Creation of Media with Prescribed Permeability. Journal of Applied Mathematics and Computing, 45, 461-485.
[24]	Ramm, A.G. and Andriychuk, M.I. (2014) Calculation of Electromagnetic Wave Scattering by a Small Impedance Particle of an Arbitrary Shape. Mathematical Modelling of Natural Phenomena, 9, 254-269. https://doi.org/10.1051/mmnp/20149517
[25]	Tran, N.T. (2015) Numerical Method for Solving Electromagnetic Scattering Problem by Many Small Impedance Bodies. Computational Physics.
[26]	Mie, G. (1908) Beiträge zur optik trüber medien, speziell kolloidaler metallösungen. Annalen der physik, 330, 377-445. https://doi.org/10.1002/andp.19083300302
[27]	Saad, Y. and Schultz, M.H. (1986) GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 7, 856-869. https://doi.org/10.1137/0907058
[28]	Nakayama, J., Ogura, H. and Sakata, M. (1981) Scattering of a Scalar Wave from a Slightly Random Surface. Journal of Mathematical Physics, 22, 471-477. https://doi.org/10.1063/1.524933
[29]	Ito, S. (1985) Analysis of Scalar Wave Scattering from Slightly Rough Random Surfaces: A Multiple Scattering Theory. Radio Science, 20, 1-12. https://doi.org/10.1029/RS020i001p00001
[30]	Tran, N.T. (2013) Numerical Solution of Many-Body Wave Scattering Problem and Creating Materials with a Desired Refraction Coefficient. The International Journal of Structural Changes in Solids, 5, 27-38.
[31]	Ramm, A.G. and Tran, N.T. (2015) A Fast Algorithm for Solving Scalar Wave Scattering Problem by Billions of Particles. Journal of Algorithms and Optimization, 3, 1-13.