VIEWS: 11 PAGES: 96 CATEGORY: Consumer Electronics POSTED ON: 11/15/2013
SINGULARITY-FREE BOUNDARY METHODS FOR ELECTROSTATICS AND WAVE SCATTERING A Thesis Presented to The Graduate Faculty of The University of Akron In Partial Fulﬁllment of the Requirements for the Degree Master of Science Osama Alkhateeb May, 2012 SINGULARITY-FREE BOUNDARY METHODS FOR ELECTROSTATICS AND WAVE SCATTERING Osama Alkhateeb Thesis Approved: Accepted: Advisor Dean of the College Dr. Igor Tsukerman Dr. George K. Haritos Faculty Reader Dean of the Graduate School Dr. Nathan Ida Dr. George R. Newkome Faculty Reader Date Dr. Joan Carletta Department Chair Dr. Alex De Abreu Garcia ii ABSTRACT Traditional boundary methods are widely used in applied electromagnetics and other areas of physics and engineering but suﬀer from the singularity of the Green function kernels. The boundary diﬀerence method (BDM) implemented and applied in this thesis alleviates this drawback by replacing the singular kernels with Green’s functions deﬁned on a discrete lattice. Although its key ideas are available in the literature, BDM has not been extensively applied in electromagnetics. The main objective of this thesis is to demonstrate the accuracy of BDM for electrostatic and wave scattering problems, and, furthermore, to extend the method to problems with edges and corners that present serious computational challenges. Another novelty of the thesis is the treatment of perfect electric conductors in the BDM framework. BDM shares one key advantage of the traditional boundary methods: the dimensionality of the problem is reduced (from 3D to 2D or from 2D to 1D). The resulting dense matrices, as opposed to sparse ones arising in ﬁnite diﬀerence and ﬁnite element analysis, can be handled in a computationally eﬃcient way using Fast Multipole acceleration or other existing techniques, with relatively minor modiﬁca- tions. The focus of the thesis, however, is on the development and construction of BDM rather than numerical solvers. iii Six case studies in the thesis demonstrate the eﬀectiveness of BDM for elec- trostatic and wave scattering problems, including ﬁeld singularities at corners and edges. The method is applied to both dielectric and perfectly conducting scatter- ers. Convergence rates of the method established via numerical experiments are in agreement with the theoretical predictions. iv ACKNOWLEDGEMENTS I would like to express my gratitude to my advisor, Dr. Igor Tsukerman, for all of his assistance, guidance, and expertise that added considerably to my research and graduate experience. I thank him for motivating and encouraging me throughout this journey, and I thank him for inspiring me and for providing me with all the support I needed. It was under his tutelage that I became interested in electromagnetics. I also appreciate his persistence, understanding, and vast knowledge through which I was able to complete this work. An appreciation also goes to my research committee for all of their support and kindness. A very special thank you goes to my family, mom and dad, brother, and sisters, for all their support, encouragement, and love. There is no doubt that it is hard for me to convey my appreciation fully to them, but I owe them my eternal love and gratitude, and as they read this I want them to know how lucky I feel for having them in my life and for being able to share this with them. v TABLE OF CONTENTS Page LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix CHAPTER I. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Overview of Boundary Equation Methods . . . . . . . . . . . . . . . 1 1.2 Objectives of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 II. ELECTROMAGNETIC THEORY . . . . . . . . . . . . . . . . . . . . . . 6 2.1 Mathematical Formulations of Electromagnetic Fields . . . . . . . . 6 2.2 Electrostatic Formulations . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 12 III. VOLUME-DISCRETIZATION METHODS . . . . . . . . . . . . . . . . . 14 3.1 The Finite Diﬀerence method . . . . . . . . . . . . . . . . . . . . . . 14 3.2 The Flexible Local Approximation Method . . . . . . . . . . . . . . 17 3.3 The Finite Element Method . . . . . . . . . . . . . . . . . . . . . . . 21 3.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 21 vi IV. THE GREEN FUNCTION . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.1 Green’s Function for Continuous Problems . . . . . . . . . . . . . . . 24 4.2 The Lattice Green Function . . . . . . . . . . . . . . . . . . . . . . . 28 4.3 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 38 V. BOUNDARY EQUATION METHODS . . . . . . . . . . . . . . . . . . . 40 5.1 The Boundary Element Method . . . . . . . . . . . . . . . . . . . . . 40 5.2 The Boundary Diﬀerence Method . . . . . . . . . . . . . . . . . . . . 45 5.3 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 53 VI. THE BOUNDARY DIFFERENCE METHOD: NUMERICAL EX- PERIMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.1 Introduction and General Assumptions . . . . . . . . . . . . . . . . . 55 6.2 FLAME Basis Functions and Exact Solutions . . . . . . . . . . . . . 56 6.3 Numerical Results for Smooth Boundaries . . . . . . . . . . . . . . . 64 6.4 Implementation of BDM with Non-Smooth Boundaries . . . . . . . . 72 6.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 VII. CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 7.1 Conclusions and Summary . . . . . . . . . . . . . . . . . . . . . . . . 80 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 vii LIST OF TABLES Table Page 4.1 Green’s function for Helmholtz’s equation in n dimensions, for n = 1, 2, 3. 28 4.2 Green’s function for Poisson’s equation in n dimensions, for n=1,2,3. . . 28 4.3 Deﬁnitions of variables in (4.16) . . . . . . . . . . . . . . . . . . . . . . 33 4.4 Terms in (4.22). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 viii LIST OF FIGURES Figure Page 3.1 5-point and 9-point stencils. . . . . . . . . . . . . . . . . . . . . . . . . 16 4.1 The real part of the lattice Green function for the 2D Helmholtz equation with k = 1, h = 1/7, and M = 50 [1]. . . . . . . . . . . . . . . 31 4.2 The real part of the lattice Green function for the 2D Helmholtz equation with k = 2, h = 1/7, and M = 50 [1]. . . . . . . . . . . . . . . 32 4.3 The real part of the lattice Green function for the 2D Helmholtz case with kh = 0.1334 and M = 50. The boundary conditions are imposed as the far ﬁeld approximation derived by Martin [2], and as the continuous Green function. . . . . . . . . . . . . . . . . . . . . . 33 4.4 The imaginary part of the lattice Green function for the 2D Helmholtz case with kh = 0.1334 and M = 50. The boundary conditions are imposed as the far ﬁeld approximation derived by Martin [2], and as the continuous Green function. . . . . . . . . . . . . 34 4.5 The lattice Green function for the 2D Poisson case with h = 1, M = 20. 37 4.6 The lattice Green function by the FD method when the boundaries are approximated by asymptotic expansions and continuous Green’s function for the 2D Poisson case with h = 1, M = 20. . . . . . . . . . . 38 5.1 Setup of an EM wave scattering problem with TM polarization. . . . . 47 5.2 Discrete double-layer boundary for simple shapes: a) Circle, b) El- lipse, and c) Square. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.3 3-point and 4-point stencils with a PEC cylinder. . . . . . . . . . . . . 51 5.4 Dielectric cylinder in presence of static uniform electric ﬁeld. . . . . . . 52 6.1 An elliptic coordinate system. Reprinted from [3]. . . . . . . . . . . . . 62 ix 6.2 Dielectric cylinder in the presence of a static uniform ﬁeld: the electrostatic potential with grid size h = 0.01. . . . . . . . . . . . . . . 66 6.3 Dielectric cylinder in the presence of a static uniform ﬁeld: relative error of the electric potential vs. grid size h in BDM. . . . . . . . . . . 66 6.4 Scattering by a circular dielectric cylinder: the real part of the electric ﬁeld with grid size h = 0.01. . . . . . . . . . . . . . . . . . . . . 67 6.5 Scattering by a circular dielectric cylinder: relative error of the electric ﬁeld vs. grid size h in BDM. . . . . . . . . . . . . . . . . . . . 68 6.6 Scattering by a circular PEC cylinder: the real part of the electric ﬁeld with grid size h = 0.01. The special behavior of the ﬁeld in the vicinity of point P is noted in the text. . . . . . . . . . . . . . . . . . . 69 6.7 Scattering by a circular PEC cylinder: relative error of the electric ﬁeld vs. grid size h in BDM and MoM. . . . . . . . . . . . . . . . . . . 69 6.8 Scattering by an elliptic dielectric cylinder: the real part of the electric ﬁeld with grid size h = 0.01. . . . . . . . . . . . . . . . . . . . . 71 6.9 Scattering by an elliptic dielectric cylinder: relative error of the electric ﬁeld vs. grid size h in BDM. . . . . . . . . . . . . . . . . . . . 71 6.10 Scattering by an elliptic PEC cylinder: the real part of the electric ﬁeld with grid size h = 0.01. Note the special behavior of the electric ﬁeld in the vicinity of point P; it is explained in Sec. 6.3.2. . . . . . . . 73 6.11 Scattering by an elliptic PEC cylinder: relative error of the electric ﬁeld vs. grid size h in BDM and MoM. . . . . . . . . . . . . . . . . . . 73 6.12 Scattering by a rectangular PEC cylinder: problem setup. . . . . . . . . 74 6.13 Scattering by a rectangular PEC cylinder: wedge problem. . . . . . . . 75 6.14 Scattering by a rectangular PEC cylinder: the real part of the elec- tric ﬁeld with grid size h = 0.01. The special behavior in the vicinity of point P outside the cylinder is discussed in Sec. 6.3.2. . . . . . . . . 76 6.15 Scattering by a rectangular PEC cylinder: the electric ﬁeld at a square located three grid layers away from the cylinder boundary by BDM and MoM, with grid size h = 0.01. . . . . . . . . . . . . . . . 77 6.16 Scattering by a rectangular PEC cylinder: convergence rate of the electric ﬁeld vs. grid size h in BDM, with wavenumber kout = 0.1. . . . 78 x 6.17 Scattering by a rectangular PEC cylinder: convergence rate of the electric ﬁeld vs. grid size h in BDM, with wavenumber kout = 0.5. . . . 78 xi CHAPTER I INTRODUCTION The Boundary Equation Methods are widely used in applications of electromagnetics such as antenna design, spectroscopic techniques of nanostructures, radar analysis, and integrated circuit design [4, 5, 6]. This introduction gives an overview of these methods and compares them with the volume-discretization methods commonly used in computational practice. The objectives of the work are stated. The last section outlines the organization of the thesis. 1.1 Overview of Boundary Equation Methods Electromagnetic problems can be formulated via diﬀerential equations in the volume or, alternatively for linear piecewise-homogeneous media, via surface integral equa- tions. The ﬁrst approach, when implemented numerically, gives rise to ﬁnite element (FE) or ﬁnite diﬀerence (FD) analysis, while the second one leads to the boundary integral / boundary element method (BIM/BEM) developed in the 1960-1970s. BEM transforms a linear boundary value problem in a piecewise-homogeneous domain into equivalent integral equations on the interface and domain boundaries. A study by Fredman [7] and Shaw in 1962 is recognized to be one of the ﬁrst implementations of BEM. They solved a scalar wave equation in the time domain. The same problem was 1 tackled by Banaugh and Goldsmith [8] one year later, but in the frequency domain. Later the method became more popular and widely used, especially after substantial contributions made by Harrington [9], [10]. According to Cheng and Cheng [11], BEM ranks as the third method amongst all the known numerical methods in terms of the number of publications for solving boundary value problems. BEM should be compared with the volume-discretization methods – the ﬁnite diﬀerence method (FD) and the ﬁnite element method (FEM) – that are used more extensively. The unknowns are typically deﬁned in the volume for FEM and FD and on surfaces for BEM. This diﬀerence can be especially dramatic for unbounded domains that must be truncated in FEM and FD using special – usually approximate – boundary conditions. In contrast, BEM handles inﬁnite domains in a natural way. In addition to this advantage, BEM reduces the dimensionality of the problem: 3D problems are reduced to 2D surfaces, etc. This is especially convenient for problems with moving boundaries. However, BEM has two major disadvantages. The boundary formulation results in fully populated matrices, in contrast with sparse matrices in the other two methods. This disadvantage can be alleviated by using the Fast Multipole Methods (FMM) [1]. Secondly, BEM works with the Green functions that appear as integral kernels and are singular. Conﬁning the computations to the boundaries is a desirable feature, but the singular kernels are often problematic. Therefore, eﬀorts have been made to develop new methods that enjoy the advantages of BEM, but do not include singular terms. 2 Boundary diﬀerence equations (BDM) are an interesting alternative from this per- spective. Two diﬀerent methodologies exist for BDM. The ﬁrst one is known as the method of diﬀerence potentials, developed by the Ryaben’kii school in [12], [13], and other publications. The Ryaben’kii theory leads to a system of algebraic equations on a discrete boundary. The second methodology is the method of boundary algebraic equations by Martinsson and Rodin [14]. It can be traced back to 1958 when Saltzer used a discrete version of Fredholm boundary integral equations to solve boundary value problems for the Poisson and Laplace equations [15]. In [14], [16], and [17] Martinsson and Rodin revisited the approach, with the main focus on mathematical analysis of the static boundary value problems. Finally, in [1] Tsukerman used a similar approach to solve an electromagnetic scattering problem. He used a numer- ical method called FLAME (see Sec. 3.2) to deﬁne the diﬀerence equations on the boundary. In the thesis we implement this latter approach and develop it further. The solution of the boundary value problem can be represented as a convolution of the lattice Green function and some artiﬁcial sources on the boundary. This approach has the advantages of BEM but does not involve any singularities because the lattice Green functions are nonsingular. 1.2 Objectives of the Thesis In this thesis electromagnetic (EM) scattering and electrostatic potential boundary value problems are solved using the boundary diﬀerence method extending the ap- 3 proach of [1]. For EM scattering, diﬀerent problems are considered: scattering by an inﬁnitely long cylinder with elliptic or circular cross sections, where the cylinder is either a dielectric or a perfect conductor. In addition to this, EM scattering by a rectangular perfectly conducting cylinder is analyzed numerically; this problem is challenging due to ﬁeld singularities at the corners. An electrostatic example is also considered: a dielectric cylinder is immersed in a static ﬁeld. In all cases, ﬁeld plots are examined and convergence rates analyzed. 1.3 Thesis Organization After the introduction and overview in the ﬁrst chapter, the second chapter reviews electromagnetic theory and Maxwell’s equations. In addition, a brief summary of the concept of boundary value problems and the most common boundary conditions are reviewed. Chapter 3 is an overview of the numerical methods used in the thesis: FD method, FLAME, and FEM. In Chapter 4 the Green function and its discrete counterpart are reviewed, since they are critical for the boundary equation method. Plots for the lattice Green function are presented, to illustrate their nonsingular behavior, in contrast with the continuous Green functions. Chapter 5 reviews BEM and BDM; they are explained in detail and model problems for EM scattering and electrostatics are considered. 4 Chapter 6 analyzes solutions of the model problems considered in Chapter 5 via ﬁeld, error rate, and comparison plots of the results obtained from BDM with those of BEM. The last chapter includes a summary of the work done in this thesis and a discussion of future work. 5 CHAPTER II ELECTROMAGNETIC THEORY The boundary equation methods can be used for solving boundary value problems of electrostatics and EM wave scattering. This chapter introduces the mathematical formulations of electrodynamics and electrostatics that are used with these problems, and it explains the boundary conditions that are commonly used with them. Our focus is on part of the theory immediately relevant to this thesis. For more details see [18, 19, 20]. 2.1 Mathematical Formulations of Electromagnetic Fields This section introduces major notation for ﬁelds, waves and the Helmholtz equation. 2.1.1 Maxwell’s Equations Electromagnetic ﬁelds are mathematically represented by four ﬁeld vectors: electric and magnetic ﬁeld intensities E, H, electric displacement and magnetic ﬁeld densities 6 D, B. The diﬀerential form of Maxwell’s equations in the time domain is ·D = ρ (2.1) ·B = 0 (2.2) ∂B ×E = − (2.3) ∂t ∂D ×H = J+ (2.4) ∂t J is a vector current density, while ρ is an electric charge density. Both J and ρ are the sources for the electromagnetic ﬁelds. In a linear isotropic homogeneous medium the relations between the ﬁeld pairs are D = r oE (2.5) B = µr µo H (2.6) Parameters o , µo , r , µr are called the permittivity and permeability in free space, and the relative permittivity and permeability of the medium, respectively. Equations (2.1) – (2.4) deal with time varying ﬁelds. At a ﬁxed frequency, one can switch to complex phasors, replace all time derivatives with jω, and obtain time-harmonic Maxwell equations: ·D = ρ (2.7) ·B = 0 (2.8) × E = − jωB (2.9) × H = J + jωD (2.10) 7 A time varying ﬁeld (e.g. E(t)) can be found from its phasor as E(t) = Re{Eejωt } (2.11) From now on, all electromagnetic ﬁelds are treated as phasor vectors. 2.1.2 The Helmholtz Equation in a Source-Free Region Fields in the Maxwell equations (2.7) – (2.10) are coupled, but either one of the ﬁelds can be algebraically eliminated, leading to a second-order equation for the remaining ﬁeld. In a source free region, if we apply the curl operation on (2.9) we get × × E = − jµω ×H (2.12) where µ = µo µr . The left hand side can be changed to another more useful form by using a vector identity, while the right hand side can be replaced using (2.10), the result is 2 ( · E) − E = µ ω2E (2.13) where = r o. By using (2.7) equation (2.13) reduces to 2 E + k2E = 0 (2.14) k is the wave number deﬁned as √ k = ω µ (2.15) Equation (2.14) is known as the Helmholtz equation; its solution represents electro- magnetic waves propagating in space, such as those in scattering problems. A similar procedure can be followed to ﬁnd an expression that involves H only. 8 2.1.3 The Poynting Vector The Poynting vector speciﬁes power ﬂow in electromagnetic waves. The instantaneous power density Poynting vector is expressed as P = E×H (2.16) Thus the direction of propagation in an EM wave is orthogonal to both the electric and the magnetic ﬁelds. 2.1.4 Interface of Materials EM scattering occurs when an electromagnetic wave propagates in a medium and hits an obstacle with diﬀerent material parameters. In the presence of two diﬀerent dielectric media with permittivities 1 and 2, the electric ﬁelds E1 and E2 for the two media, respectively, must satisfy the interface conditions n· 1 E1 = n· 2 E2 (2.17) n × E1 = n × E2 (2.18) where the surface free sources (charges and currents) are equal to zero, and n is a unit normal vector from medium 1 to medium 2. (2.17) and (2.18) refer to normal and tangential components of the ﬁeld at the interface, respectively. Note that (2.17) predicts a discontinuity in the electric ﬁeld. Similar expressions are available to describe the magnetic ﬁeld at the interface, under the assumptions made (J = 0 and 9 ρ = 0 at the interface). n · µ1 H1 = n · µ2 H2 (2.19) n × H1 = n × H2 (2.20) Similarly, equations (2.19) and (2.19) involve the normal and the tangential compo- nents of the magnetic ﬁeld, respectively. 2.1.5 TE and TM Polarizations The term polarization speciﬁes the direction of the electric ﬁeld in certain regions of space. At material interfaces, the plane of incidence is deﬁned as the plane that includes both a vector that represents the direction of wave propagation, and a vector normal to the interface. If the electric ﬁeld vector is parallel to the plane of incidence, then the wave is called transverse-magnetic (TM). Otherwise, it is transverse-electric (TE). Other deﬁnitions for polarizations are not considered, as they are not used in this work. 2.2 Electrostatic Formulations The electric ﬁeld that does not depend on time, such as the one considered in the model problem in Sec. 5.2.5, is called electrostatic. From equation (2.3), the curl of the electric ﬁeld is zero in this case, and the electric and magnetic ﬁelds are no longer coupled. According to vector calculus, a vector with a zero curl in a simply connected domain is a gradient of a scalar function. This function is called the 10 electrostatic potential: E = − V (2.21) Likewise, if the electric ﬁeld is known, the potential can also be evaluated by the line integral V = − E · dl (2.22) The boundary conditions in equations (2.17) and (2.18) are still valid in electrostatics. However, sometimes it is more useful to formulate those conditions in terms of the potential. V1 = V2 (2.23) ∂V1 ∂V2 1 = 2 (2.24) ∂n ∂n where the V1,2 refer to the potential along the interface in media 1 and 2, respectively. Expressions for time-independent magnetic ﬁelds are available in the literature but will not be included here, as they are not going to be used later. 2.3 Boundary Conditions There are several types of boundary conditions, three of them being most common in static problems. The ﬁrst one is the Dirichlet boundary condition, in which the value of the function to be solved for is deﬁned on the boundary. The second one is the Neumann boundary condition, where the normal derivative of the function is deﬁned on the boundary. Finally, the Robin boundary condition involves a linear combination of the ﬁrst two conditions. 11 2.3.1 The Sommerfeld Radiation Condition Wave problems require diﬀerent types of boundary conditions, an example of which is the Sommerfeld radiation condition. The Helmholtz equation deﬁned in (2.14) describes the behavior of the time-harmonic electric ﬁeld in a source-free region. If this ﬁeld is radiated by a current density J, a new term should be added to the equation to represent it. To do so, we use the derivation applied in Sec. 2.1.2, to ﬁnd that 2 E + k 2 E = jµωJ (2.25) In general the ﬁeld is a superposition of incoming and outgoing waves. The Sommer- feld radiation condition implies that the energy is radiated away from the sources (J in our example). In two dimensions, assuming the electric ﬁeld has only a component in the z direction, this condition has the form √ ∂ lim r + jk Ez = 0 (2.26) r→∞ ∂r with r = x2 + y 2 . This condition is imposed on the Green functions and in the scattering problems that are discussed in Chapters 4 and 5, respectively. 2.4 Summary and Conclusions Electromagnetic theory reviewed in this chapter forms a basis for the electrostatic and EM wave scattering problems considered in this thesis. Maxwell’s equations, in the time domain and in the time-harmonic form, are reviewed in Sec. 2.1. These equations relate the electromagnetic ﬁelds (E, H, D, B) and sources (J, ρ). From 12 Maxwell’s equations, we derived the Helmholtz equation for the electric ﬁeld in a source-free region and in a region with an electrical current density. This equation can be treated numerically via FD, FEM, BEM, or BDM analysis. Equations of electrostatics, as a particular case of electrodynamics, are given in Sec. 2.2. These equations and the related interface conditions can be written in terms of either the electric ﬁeld or, alternatively and often more conveniently, the electrostatic potential. Sec. 2.3 reviews common boundary conditions, including the Sommerfeld radiation condition for outgoing EM waves. These subjects form a basis for all subsequent chapters. 13 CHAPTER III VOLUME-DISCRETIZATION METHODS This chapter summarizes two methods that are used in the thesis: the ﬁnite diﬀerence method and the ﬂexible local approximation method. Additionally, it provides a very brief summary of the ﬁnite element method. The boundary equation methods are introduced in a separate chapter. 3.1 The Finite Diﬀerence method In ﬁnite diﬀerence (FD) analysis a diﬀerential equation is discretized and solved on a grid. Even though the ﬁnite element method (FEM) is widely recognized to be more powerful, FD has less overhead and is especially valuable for problems with simple geometries [21]. In this thesis FD schemes, with their respective lattice Green functions, are used in Chapter 4. Details about the subject are available in the literature, e.g. [21]. 3.1.1 Introduction to Diﬀerence Schemes On a 2D grid, the distance between two adjacent nodes is called the grid size h, with hx and hy for the two directions. Consider a function f (x, y) deﬁned in a domain D. 14 For some integers (n, m) ∈ Z2 the discretized version of f admits fn,m = f (xn = x0 + nhx , ym = y0 + mhy ) (n, m) ∈ Z2 (3.1) where (x0 , y0 ) ∈ D. Discretization is of course critical in computer simulations, as only a ﬁnite amount of information can be processed (in FD, solution on a ﬁnite grid). It is of course important to approximate the derivatives of functions on the 2 ∂2 ∂2 grid. For example, consider the 2D Laplace operator ≡ ∂x2 + ∂y 2 applied to a scalar function u(x, y). A diﬀerence operator can be derived on a grid using the Taylor series ∂u h2 ∂ 2 u h3 ∂ 3 u h4 ∂ 4 u un+1,m = un,m + hx |xn ,ym + x 2 |xn ,ym + x 3 |xn ,ym + x 4 |xn ,ym + . . . (3.2) ∂x 2! ∂x 3! ∂x 4! ∂x and ∂u h2 ∂ 2 u x h3 ∂ 3 u x h4 ∂ 4 u un−1,m = un,m − hx |xn ,ym + | 2 xn ,ym − | 3 xn ,ym + x 4 |xn ,ym + . . . (3.3) ∂x 2! ∂x 3! ∂x 4! ∂x ∂2u By performing straightforward algebra | ∂x2 xn ,ym admits ∂ 2u un−1,m − 2un,m + un+1,m | 2 xn ,ym = + O(h2 ) x (3.4) ∂x h2 x The O(h2 ) terms signify that the error due to the truncation of the series is of order x ∂2u 2. A similar expression for | ∂y 2 xn ,ym can be derived ∂ 2u un,m−1 − 2un,m + un,m+1 | 2 xn ,ym = 2 + O(h2 ) y (3.5) ∂y hy Hence, the Laplace operator is the sum of equations (3.4) and (3.5), and the overall error is of order 2. 15 Similarly, if a uniform grid h = hx = hy is considered, the 2D Helmholtz 2 equation ( E + k 2 E = 0) converts to En,m−1 + En−1,m − 4En,m + En+1,m + En,m+1 + h2 k 2 En,m = 0 (3.6) where E is a time-harmonic electric ﬁeld, and k is the wave number. From the scheme it is obvious that each node (n, m) is connected to four neighbor nodes. These nodes together are known as a 5-point stencil, and the diﬀerence scheme as a 5-point scheme. A similar approach is followed to derive a diﬀerence scheme for the 2 ∂2u 2D Poisson equation ( ∂ u + ∂x2 ∂y 2 = f ); this scheme reads un,m−1 + un−1,m − 4un,m + un+1,m + un,m+1 = h2 f (3.7) Both (3.6) and (3.7) are of order 2. To increase the accuracy of the schemes, the stencil has to be expanded and the derivation should involve more terms in the Taylor expansion. Figure 3.1 shows how a 5-point and 9-point stencils are formed at a center grid node (n, m). Figure 3.1: 5-point and 9-point stencils. 16 3.1.2 Assessment of FD The ﬁnite diﬀerence method can be easily applied to problems with simple geometries. However, performance of the method degrades for problems with curved or slanted boundaries, where the grid does not conform with the boundaries because the Taylor expansion breaks down at the interface. Reﬁning the grid to increase accuracy leads to higher computational cost. One alternative is to use grids conforming to complex boundaries, as done in FEM. Another possibility is to apply ﬂexible approximation, as discussed in the following section. 3.2 The Flexible Local Approximation Method The ﬂexible local approximation method (FLAME) is a generalized FD calculus where accurate local approximations of the solution, rather than generic Taylor expansions, are used to generate diﬀerence schemes. FLAME was implemented and used with many applications in [21, 22, 23, 24, 25, 26]. FLAME uses simple grids that may not conform with the interface boundaries between materials, but still provides accurate solutions. This can be done by replacing the Taylor expansion (used with FD) by some local functions, that satisfy the underlying diﬀerential equation. In this thesis FLAME is used with the boundary diﬀerence method BDM in last two chapters. 17 3.2.1 Treﬀtz Schemes for Homogeneous Equations In FLAME the computational domain Ω is covered by a set of overlapping patches Ω = ∪Ω(i) , i = 1, 2, . . . , n, with n being the number of patches. These patches are needed to approximate the solution locally. To derive FLAME Treﬀtz schemes for homogeneous equations, consider, following Tsukerman [21], a homogeneous diﬀerential equation Lu = 0 in Ω(i) (3.8) with a diﬀerential operator L. A local approximation space of the solution for each patch admits (i) Ψ(i) = span(ψα , α = 1, 2, . . . , m(i)) (3.9) (i) (i) where ψα is a basis function. The patch-wise approximation of the solution uh can (i) be written as a summation function including m(i) basis functions ψα accompanied (i) with some arbitrary coeﬃcients cα m(i) (i) (i) (i) uh = cα ψα (3.10) α=1 (i) For a M -point stencil the coeﬃcient vector c(i) ≡ cα ∈ Rm and the vector u(i) ≡ (i) uh ∈ RM (M and m need not to be the same) are related through a transformation matrix N (i) u(i) = N (i) c(i) (3.11) 18 N (i) contains the values of the basis functions for each node on the stencil; given the position vector rk located at node k, N (i) is represented explicitly (i) (i) (i) ψ1 (r1 ) ψ2 (r1 ) ... ψm (r1 ) ψ (i) (r ) ψ (i) (r ) (i) . . . ψm (r2 ) 1 2 2 2 N (i) = (3.12) . . . . .. . . . . . . (i) (i) (i) ψ1 (rM ) ψ2 (rM ) . . . ψm (rM ) The diﬀerence scheme vector s(i) shall be chosen such that: s(i)T u(i) = 0 (3.13) By replacing u(i) with its equivalent term in (3.11), equation (3.13) becomes s(i)T N (i) c(i) = 0 (3.14) For an arbitrary coeﬃcient vector c(i) , equation (3.14) holds if s(i) ∈ NullN (i)T (3.15) where Null refers to the null space operator. FLAME can be implemented in 3 steps: 1) ﬁnd an appropriate set of basis functions (m basis functions), 2) compute the transformation matrix N (i) (evaluate the basis functions on the M -point stencil), 3) ﬁnd the scheme coeﬃcients from the null space of N (i) . Usually, for a M -point stencil, M − 1 basis functions are needed (m = M − 1) so that the vector s(i) is one dimensional [21]. 19 3.2.2 Treﬀtz Schemes for Inhomogeneous Equations In the above section FLAME schemes are discussed for homogeneous problems (equa- tion (3.8)). The method is used also to solve an inhomogeneous diﬀerential problem of the form Lu = f in Ω(i) (3.16) where f is a known function of coordinates. The total solution for u can be di- vided into two: the homogeneous solution u0 , and the particular solution uf that corresponds to the inhomogeneous term f . u = u0 + uf (3.17) We apply the FLAME scheme derived for the homogeneous problem on equation (3.17) (i) s(i)T u(i) = s(i)T uf (3.18) uf is usually easy to be constructed because it is local (it does not have to satisfy a certain boundary condition). 3.2.3 Assessment of FLAME When compared with the FD method, FLAME gives a better approximation of the solution on the interface boundaries. The FLAME grid can be uniform and does not have to conform to the boundary. Thus, it can be a powerful tool in some applications, such as electromagnetic scattering problems for simple shapes, as in cylinders with circular or elliptic cross-sections. Therefore, this method is utilized in the subsequent chapters. 20 3.3 The Finite Element Method The ﬁnite element method (FEM) is applied to a variational problem, typically ob- tained from the diﬀerential formulation of a boundary value problem via a weighted residual (in particular, Galerkin) method. FE meshes are usually more complex than regular FD grids and conform to the interface and exterior boundaries. FE meshes are composed of discrete sub-domains known as elements. The number and the size of elements control the accuracy of the method. FEM is recognized as a powerful and a versatile method; it can be used to solve problems with complex geometries and materials with variable or even nonlinear parameters. In addition to this, FEM leads to sparse matrices – a desirable feature that reduces the computational cost and memory. However FEM needs relatively complex data structures, especially with prob- lems of complex geometries. Additionally, inﬁnite domains need to be truncated. More details about the method are provided in [21]. 3.4 Summary and Conclusions Three numerical methods are reviewed in this chapter. The ﬁrst one is the ﬁnite diﬀer- ence method. In particular, 5-point diﬀerence schemes for the Possion and Helmholtz equations, with their respective lattice Green functions, are used in Chapter 4. Per- formance of classical FD methods degrades at curved or slanted boundaries. This can be remedied by the ﬂexible local approximation method (FLAME), reviewed in the 21 chapter for both homogeneous and inhomogeneous (nonzero right hand side) equa- tions. Local Treﬀtz bases, used to construct FLAME schemes, provide an accurate approximation of the solution at interface boundaries. Examples of FLAME, with the respective basis functions, are given in Chapter 6. FLAME is used to construct boundary diﬀerence schemes in the last two chapters. A brief description of the ﬁnite element method is included in the chapter for completeness and comparison. 22 CHAPTER IV THE GREEN FUNCTION Green’s function is the response of a system to a delta-source (point source); it is use- ful as it can be used to ﬁnd the response to arbitrary sources. For example, in signal processing, the output signal of a time invariant system is evaluated by convolving the input signal and the Green function (impulse response). The Green function was originally introduced by George Green in [27]; it plays an important role in mathe- matics, where many methods, both analytical and numerical, use this function as a basic tool for solving initial value and boundary value problems. Integral equations are an example of Green-function-based methods that can be solved numerically, e.g. using boundary elements [28]. The review in this chapter focuses on Green’s function for the Poisson and Helmholtz equations, because they are very important in electromagnetics in general, and for the applications considered in this thesis in particular. In addition to the continuous Green function, its discrete analogue on the grid is also considered and is the key for numerical methods developed in Chapter 6. 23 4.1 Green’s Function for Continuous Problems The Green function is a useful tool for solving linear initial or boundary value prob- lems with sources. It helps to transform the problem from diﬀerential form into an integral one, with the Green function as the kernel. 4.1.1 General Review The Green function is mathematically deﬁned as the solution of the equation LG(r, r ) = − δ(r − r ) (4.1) where L represents a linear diﬀerential operator, G is the Green function, and δ is the Dirac delta function. Equation (4.1) is considered in a ﬁnite or inﬁnite domain, and so is subject to conditions imposed on the domain boundary. The solution of (4.1) depends on L and on the boundary conditions. Equation (4.1) describes the 2 ﬁeld created by a unit source located at a point r = r . If L is replaced with ,G would represent the electrostatic potential of a point charge. As already noted, once the Green function has been found, one can calculate the response to an arbitrary source by representing this source as a superposition of delta sources. More precisely, a general solution can be expressed by the convolution of the Green function with the sources. For more details about the Green function see [28, 29, 30]. 24 4.1.2 Green’s Formulas Green’s formulas play a central role in the analysis of boundary value problems for the Poisson and related problems. The ﬁrst Green formula is 2 (ϕ χ) · n dS = ( ϕ· χ+ϕ χ) dV (4.2) S V where χ and ϕ are scalar functions, and n represents an outward normal unit vector on a surface S enclosing a volume V . The second formula, or Green’s theorem, is 2 2 (ϕ χ − χ ϕ) · n dS = (ϕ χ−χ ϕ) dV (4.3) S V It is clear that equations (4.2) and (4.3) are written in 3D, but can easily be used in 2D by replacing volume integrals with surface integrals, and closed surface integrals with contour integrals. Green’s formulas will be used in subsequent sections to convert boundary value problems to integral equation problems. Details and proofs of those formulas are available in many textbooks of mathematical physics, e.g. [29, 30]. 4.1.3 Green’s Function for Helmholtz equation The Helmholtz equation plays an important role in electromagnetic ﬁeld theory in general and in the analysis of EM wave scattering in particular. Green’s function for the Helmholtz problem satisﬁes the following equation by deﬁnition: 2 G(r, r ) + k 2 G(r, r ) = − δ(r − r ) (4.4) where k is the wave number. G satisﬁes the Sommerfeld radiation condition (Sec. 2.3.1). We are not deﬁning the boundary conditions yet, so G is not uniquely deﬁned. 25 To understand how the Green function can be useful in this case, consider the problem 2 u + k 2 u = f in V (4.5) where u is the unknown solution within the volume V , and f is a source function. From straightforward algebra and Green’s second formula (4.3), equation (4.5) can be converted to an integral form u(r) = [G(r, r ) u(r ) − u(r ) G(r, r )] · n dS + f (r )G(r, r ) dV (4.6) S V Green’s function is known to have the symmetry (reciprocity) property: G(r, r ) = G(r , r) (4.7) At this stage, it is important to discuss the boundary conditions that are most com- monly imposed with (4.4). In the Dirichlet problem G vanishes at the boundary, while in the Neumann problem G · n vanishes at the boundary. The solution of equation (4.6) depends on the source function f and on the values of u and u on the boundary S. Assuming a nonhomogeneous Helmholtz equation with a homogeneous Dirichlet boundary condition, the surface integrals in (4.6) vanishes, as both u and G have zero values on the boundary, and (4.6) becomes u(r) = f (r )G(r, r ) dV (4.8) V The convolution of f and G in (4.8) simpliﬁes if G is translationally invariant, i.e. G(r, r ) = G(r − r ) (4.9) 26 Translational invariance holds if the diﬀerential operator has constant coeﬃcients, i.e. if the medium is homogeneous throughout. Solutions for Green’s functions are provided in many textbooks for diﬀerent problems, e.g. [28, 29, 30, 31]. For the phasor convention ejωt , Table 4.1 provides expressions for G in n dimensions (n = 1, 2, 3). They represent outward waves from a point source. Note that in 2D and 3D G is singular at the origin – a critical point for further discussion in the subsequent chapters of the thesis. 4.1.4 Green’s Function for Poisson’s equation The Poisson equation can be thought of as a special case of the Helmholtz equation, with k = 0. Therefore, Equations (4.6) and (4.8) still hold. Green’s function is diﬀerent, because (4.4) is replaced with 2 G(r, r ) = − δ(r − r ) in V (4.10) The Green function for this problem is known to be 1/4πr, where r = |(x − x )2 + (y − y )2 + (z − z )2 |. Substituting the Green function in equation (4.8) and assuming a charge distribution inside V , the result would be the well known potential integral equation ∞ ∞ ∞ 1 ρ(x , y , z ) u(x, y, z) = dx dy dz (4.11) 4π −∞ −∞ −∞ |(x − x )2 + (y − y )2 + (z − z )2 | This is another example of expressing the solution as a convolution of given sources with Green’s function. Table 4.2 expresses G in n dimensions (n = 1, 2, 3). As in the case of the Helmholtz equation (Sec. 4.1.3), G is singular at the origin in 2D and 3D. 27 Table 4.1: Green’s function for Helmholtz’s equation in n dimensions, for n = 1, 2, 3. n Green’s function −i −ik|x−x | 1 2k e i (2) 2 − 4 Ho (kr) e−ikR 3 4πR r= (x − x )2 + (y − y )2 , R = (x − x )2 + (y − y )2 + (z − z )2 Table 4.2: Green’s function for Poisson’s equation in n dimensions, for n=1,2,3. n Green’s function |x−x | 1 2 2 − ln(r) 2π 1 3 4πR r= (x − x )2 + (y − y )2 , R = (x − x )2 + (y − y )2 + (z − z )2 4.2 The Lattice Green Function The continuous Green function discussed in Sec. 4.1 forms a basis of integral equation / moment methods [9, 10]. The numerical method implemented in this thesis uses the discrete version of the Green function or simply known as the lattice Green function (LGF). The solution of the lattice Green function depends on the governing diﬀerence operators for each problem, the characteristics of the medium, and in the 28 case of the Helmholtz equation on the mesh size. Equation (4.8) remains valid for LGF if the continuous convolution is replaced with a discrete one. Therefore, the central question is to evaluate the LGF. Further discussion will be limited to the 2D Laplace and Helmholtz operators as they are the case studies taken in this thesis. 4.2.1 Lattice Green’s Function for Helmholtz’s Equation In [32, 33] the LGF is given in a general hypergeometric integral that can be evaluated numerically. [2] derives a far ﬁeld approximation using two diﬀerent methods; the ﬁrst method is due to Koster, while the second one uses an idea of Lifshitz. Finally, [1] employs the ﬁnite diﬀerence (FD) method mentioned in Chapter 3 to solve for the LGF numerically, where the boundary conditions are approximated from the continuous Green function. The expressions given in [2] can also be used on the boundaries when we ﬁnd the LGF by the ﬁnite diﬀerence method. In this thesis, the LGF g (lower case g is used to diﬀerentiate LGF from the continuous Green function G) is found via the FD method. The LGF by deﬁnition satisﬁes: Lg(n) = − δ(n) (4.12) with a double index integer n ≡ (n1 , n2 ). In general L is any diﬀerence operator approximating the Helmholtz operator. The boundary conditions can be approxi- mated by the continuous Green function just as in [1]. In the far ﬁeld g approaches the continuous one, and therefore, the LGF can be accurately computed in a ﬁnite size spatial window with G imposed as a boundary condition on that window. An 29 alternative somewhat diﬀerent way of imposing those boundary conditions is due to Martin [2], who relies on the method previously developed by Koster. Below the numerical results for those two boundary conditions are compared. The simplest FD scheme valid for equation (4.12) is the 5-point scheme Lg(n) ≡ g(n1 −1, n2 )+g(n1 +1, n2 )+g(n1 , n2 −1)+g(n1 , n2 +1)+(h2 k 2 −4)g(n1 , n2 ) (4.13) where h is a spatial grid size (for simplicity assumed the same in both directions), and δ(n) is the discrete delta function equal 1 at the origin, and zero everywhere else in the lattice. As mentioned above the lattice Green function can be computed in a ﬁnite spatial window with a continuous Green’s function imposed as a boundary i (2) condition on that window. For the Helmholtz equation g = − 4 Ho (kr) (see Table 4.1), within the spatial window g can be found by solving FD system of equations (4.12), (4.13). Figures 4.1 and 4.2 show the real part of the FD solution for g on a large square [−M, M ]2 where M = 50, with k = 1 and k = 2 respectively, and h = 1/7. Note that g is bounded everywhere, unlike the continuous counterpart that has a singularity at the origin. This is an advantage of the LGF especially in cases where the singular kernel of the continuous Green function cannot be handled analytically [1]. The second approach that will be used here to ﬁnd g is based on the Discrete Time Fourier Transform (DTFT) for equation (4.12); the result of the transform is F[Lg] ≡ (exp(−jξ1 ) + exp(jξ1 ) + exp(−jξ2 ) + exp(jξ2 ) − 4 + k 2 h2 )F(g) = 1 (4.14) 30 Figure 4.1: The real part of the lattice Green function for the 2D Helmholtz equation with k = 1, h = 1/7, and M = 50 [1]. because the Fourier transform for the discrete delta function equals 1. Then g is given by the Inverse Discrete Time Fourier Transform (IDTFT) 1 ej(n1 ξ1 +n2 ξ2 ) g(n1 , n2 ) = dξ (4.15) (2π)2 (−π,π)2 k 2 h2 − 4 + e−jξ1 + ejξ1 + e−jξ2 + ejξ2 ∀ξ ∈ (−π, π)2 . The diﬃculty here is that the expression above has several poles. For example, if kh → 2, the poles will occur at four points as (ξ1 , ξ2 ) = (l1 π , l2 π ), where 2 2 l1,2 = ±1. However, Martin [2] was able to handle the singularities in his formulations for the far ﬁelds. In the method of Koster he used, the oscillatory integral in (4.15) is converted to a three dimensional integral; then, the method of stationary phases is 31 Figure 4.2: The real part of the lattice Green function for the 2D Helmholtz equation with k = 2, h = 1/7, and M = 50 [1]. applied. The resulting expression for g is given as iπ 1 ei(n1 υ+n2 ν) e 4 {4 − k 2 h2 (cos4 θo + sin4 θo )} 4 g(n1 , n2 ) ∼ − √ (4.16) 2πkhR (4 − k 2 h2 )(2 − k 2 h2 cos2 θo sin2 θo ) The variables appearing in (4.16) are deﬁned in Table 4.3 as functions of n1 , n2 , and k. √ Note that the expressions for υ and ν hold true only if kh ∈ [0, 2]. For kh ∈ (2, 2 2] iπ −iπ the e 4 factor in (4.16) would be replaced by its complex conjugate e 4 , and for √ kh > 2 2 no real solution exists. We applied (4.16) as the boundary condition for the FD problem in (4.12). Figures 4.3 and 4.4 compare this solution with the one obtained previously, with the continuous Green function (i.e. the Hankel function) imposed as the boundary condition. Figure 4.3 shows the real part of g as a function 32 Table 4.3: Deﬁnitions of variables in (4.16) R = n2 + n2 1 2 α = tan−1 ( n2 ) n1 υ = 2 sin−1 ( kh cos θo ), ν = 2 sin−1 ( kh sin θo ), 0 < kh < 2 2 2 √ 8−k2 h2 √ 8−k2 h2 √ υ = 2 cos−1 ( 2 cos θo ), ν = 2 cos−1 ( 2 sin θo ), 2 < kh < 2 2 √ 2(1−tan2 α) θo = tan−1 ( −λ + λ2 + tan2 α), λ = 4−kh2 Figure 4.3: The real part of the lattice Green function for the 2D Helmholtz case with kh = 0.1334 and M = 50. The boundary conditions are imposed as the far ﬁeld approximation derived by Martin [2], and as the continuous Green function. 33 Figure 4.4: The imaginary part of the lattice Green function for the 2D Helmholtz case with kh = 0.1334 and M = 50. The boundary conditions are imposed as the far ﬁeld approximation derived by Martin [2], and as the continuous Green function. of n1 for n2 = 0, while Figure 4.4 shows the imaginary part of g. From the ﬁgures one can see that the solutions obtained from the two boundary conditions are in a very good agreement because the spatial window m = 50 is large enough for the asymptotic approximation on the boundary to be accurate. 4.2.2 Lattice Green’s Function for The Poisson Equation There are at least two general approaches for ﬁnding the LGF associated with the discrete Poisson equation in 2D: the Discrete-Time Fourier Transform (DTFT) anal- ysis, and the Finite Diﬀerence method. The former one was used extensively by 34 Martinsson in [14, 16, 17]. In [17] a numerical integration was done by handling the singularity of the DTFT analytically. In [14] the boundary condition for the FD solu- tion is estimated by asymptotic expansions in DTFT. The continuous Green function can also be used in deﬁning the boundary conditions just as we did in the Helmholtz case. To apply the Fourier transform method consider equation (4.12), with the ﬁve-point diﬀerence operator L Lg(n) ≡ g(n1 −1, n2 )+g(n1 +1, n2 )+g(n1 , n2 −1)+g(n1 , n2 +1)−4g(n1 , n2 ) (4.17) Following the procedure of Sec. 4.2.1, the IDTFT associated with this problem can be expressed as 1 e−j(n1 ξ1 +n2 ξ2 ) g(n1 , n2 ) = dξ (4.18) (2π)2 (−π,π)2 4 − e−jξ1 − ejξ1 − e−jξ2 − ejξ2 The right hand side of the equation (4.18) does not have a closed form solution, and it has removable singularity (see Sec. 4.2.1). To get rid of the singularity, one integral can be solved analytically, then the second one can be evaluated numerically. After some algebraic manipulations and applying the residue theorem, equation (4.18) simpliﬁes to π 1 g(n1 , n2 ) = [cos(n2 ξ2 )fn1 (1 − cos ξ2 ) − f0 (1 − cos ξ2 )]dξ2 (4.19) (2π)2 0 where the new function f is deﬁned as √ (1 + α − α2 + 2α)k fk (α) = π √ (4.20) α2 + 2α 35 fk (1 − cos ξ2 ) blows up with ξ2 → 0; still this can be handled in (4.19) as lim [cos(n2 ξ2 )fn1 (1 − cos ξ2 ) − f0 (1 − cos ξ2 )] = − πn1 (4.21) ξ2 →1 Hence, integration (4.19) can be performed numerically. The integration was imple- mented in Mathcad, with the precision set to 10−10 . Figure 4.5 shows the solution for the lattice Green function over square window [−20, 20]2 , with h = 1. g is purely real, just as is its continuous counterpart. But g is bounded everywhere; actually it is zero at the origin. Numerical integration is time consuming, and few programming platforms are capable of computing improper integrals. Therefore a good alternative way to solve for g is by using the FD method. The boundary conditions for FD can be approximated in two ways: (1) using expressions derived by Martinsson in [14]; (2) using the continuous Green function from Table 4.2. Martinsson’s expression comes from an asymptotic analysis of equation (4.18), where g is represented as a sum gQ = ΣQ gl l=0 (4.22) Q represents the number of terms that one can choose, depending on the desired accuracy. gl are functions deﬁned in the polar coordinates (r, θ) and given in Table 4.4. Two conclusions can be drawn. From the table it is clear that all g’s blow up as r → 0, but this is not a problem since we are going to use them only at the boundaries. Secondly g0 is in fact the continuous Green function (Table 4.2), Therefore, other terms can be viewed as corrections. 36 Figure 4.6 compares two solutions for g with the boundary conditions are taken from (4.22) and Table 4.2. The agreement between the two results is almost ideal because the spatial window was very large [−20, 20]2 . Figure 4.5: The lattice Green function for the 2D Poisson case with h = 1, M = 20. Table 4.4: Terms in (4.22). l Term 1 0 − 2π ln(r) + c, c = 0.2573434264137 1 cos(4θ) 1 24π r2 1 25 cos(8θ)+18 cos(4θ) 2 480π r4 1 490 cos(12θ)+459 cos(8θ) 3 2016π r6 37 Figure 4.6: The lattice Green function by the FD method when the boundaries are approximated by asymptotic expansions and continuous Green’s function for the 2D Poisson case with h = 1, M = 20. 4.3 Summary and Conclusions The Green functions are central in boundary equation methods. More speciﬁcally, the continuous-domain Green functions are the kernels in the boundary integral equa- tions, and the lattice Green functions play a central role in the discrete convolutions in the boundary diﬀerence methods, as the one implemented in this thesis. There- fore, in this chapter both the continuous-domain and the lattice Green functions have been deﬁned and reviewed. It is shown how discrete convolutions with lat- tice Green functions are used to approximate boundary value problems numerically. Ways of computing the lattice Green function for the discrete Helmholtz and Poisson equations in 2D were discussed in Subsections 4.2.1 and 4.2.2. The discrete Green 38 functions were obtained using two diﬀerent methods, the results being in a very good agreement. This makes the LGF a reliable tool to be used in the subsequent chapters in this thesis. Additionally, it is shown how the LGF is bounded everywhere, unlike its continuous-domain counterpart, making it convenient for numerical implementation. 39 CHAPTER V BOUNDARY EQUATION METHODS In boundary equation methods, all the unknowns are conﬁned to boundaries. There- fore, these methods handle unbounded domains naturally, and they reduce the dimen- sionality of the problem. In this chapter we summarize these methods and formulate model problems of EM scattering and electrostatics that are used in Chapter 6. 5.1 The Boundary Element Method The boundary element method translates a boundary value problem into boundary integral equations. In electromagnetics (electrostatics and wave scattering), it has been popularized by Harrington as the method of moments (MoM) [9, 10]. The boundary integral equation has a singular kernel – the Green function. Therefore, an analytic expression for the integration is usually unavailable; but it can be approxi- mated. In this section we provide the mathematical formulation of the method, and construct a model problem that utilizes approximations done by Harrington. The model problem is used in Chapter 6. 40 5.1.1 Mathematical Formulation of MoM Consider, following Harrington [9, 10], the inhomogeneous diﬀerential equation (3.16) with a domain Ω and a boundary Γ. In MoM u is expressed as a combination of some basis functions ψn , linearly independent and deﬁned on Γ: u = αn ψn (5.1) n where αn are some coeﬃcients. Naturally, in numerical simulation the series (in general inﬁnite) is truncated to a ﬁnite number of terms. To discretize (3.16), one applies expansion (5.1) to the weak (weighted-residual) form of (3.16), which leads to αn < Wm , Lψn > = < Wm , f > (5.2) n where < Wm , Lψn > is a suitable linear functional. This expression can be written in the matrix form Iα = F (5.3) with the entries < W1 , Lψ1 > < W1 , Lψ2 > . . . I = < W2 , Lψ1 > < W2 , Lψ2 > . . . (5.4) . . . . .. . . . α1 α = 2 α (5.5) . . . 41 and < W1 , Lf > F = < W2 , Lf > (5.6) . . . Finally, in the same manner if Ψ represents a column matrix of the entries ψn , from (5.1) and (5.3) the solution for u has the form u = ΨI −1 F (5.7) where I should be nonsingular, in order to apply the inverse operation. In the Galerkin method, the weighting functions Wn are chosen the same as the expansion functions, Wn = ψn . The accuracy of the solution is aﬀected by the choices of Wn and ψn [9]. 5.1.2 Example: 2D Electromagnetic Scattering by Conducting Cylinders In this section we consider a model electromagnetic scattering problem, an important example to be utilized in Chapter 6. The problem is summarized as follows: a TM polarized plane wave hits an inﬁnitely long cylinder that is a perfect conductor with an arbitrary cross section. We need to solve for the electric ﬁeld outside the cylinder (inside the cylinder the electric ﬁeld is zero). Since the cylinder is inﬁnitely long, the dimension of the problem is reduced to 2, while the boundary of cylinder is reduced to a contour boundary C. For TM polarization, the incident electric ﬁeld Einc is z Einc = E0 exp(−jkr cos φ) z (5.8) 42 in polar coordinates, with the ﬁeld magnitude E0 , and the wavenumber k. The total electric ﬁeld on the boundary C admits Ez = Einc + Esc = 0, z z on C (5.9) where Esc refers to the scattered ﬁeld. z A ﬁlament current density J, deﬁned on a contour boundary L, radiates the electric ﬁeld E(r) −kη (2) E(r) = J(r )H0 (k|r − r |) dl (5.10) 4 L where H (2) is the Hankel function of the second kind, and η is the intrinsic impedance. In our example, the incident ﬁeld induces a current density on the contour of the scatterer. Accordingly, a scattered ﬁeld is produced due to the current. Based on that, by using equations (5.9) and (5.10), the incident ﬁeld Einc has the form z kη (2) Einc (r) = z Jz (r )H0 (k|r − r |) dl , r on C (5.11) 4 C on C. The only unknown in this expression is the current density Jz ; once it is determined the total electric ﬁeld can also be evaluated everywhere by using equation (5.10). If the contour boundary C is divided into N segments, the very basic way to solve (5.11) by MoM is by using a pulse function for the basis functions 1, on ∆Cn (5.12) ψn = 0, on all other ∆Cm (5.13) From equation (5.1) we let Jz = n αn ψn , then we substitute it in equation (5.11). The resultant equation is usually chosen to be satisﬁed at the midpoint (xm , ym ) on 43 the segments ∆Cm . It has the same form of the matrix equation in (5.3), with the elements of the vector F are fm = Einc (xm , ym ) z (5.14) while the elements of the transformation matrix I are kη (2) Im,n = H0 (k (x − xm )2 + (y − ym )2 ) dl (5.15) 4 ∆Cn As in equation (5.7) the solution for the current density Jz is Jz = ΨI −1 F . Integral (5.15) does not have an analytic expression, but various approxima- tions can be derived for it. The zero-order approximation of this integral is furnished by piece-wise constant functions; i.e., (5.15) can be approximated by a ﬁlament cur- rent Jz ∆Cn when n = m. Accordingly, the integration reduces to kη (2) Im,n = ∆Cn H0 (k (xn − xm )2 + (yn − ym )2 ) (5.16) 4 (2) The diagonal In,n needs special treatment, because the Hankel function (H0 ) (see Tabel 4.1) has a removable singularity. Hence, it must be handled analytically. To accomplish this, the segments ∆Cn are approximated by straight lines, and a small argument approximation is employed with the Hankel function (2) 2 γz H0 (z) = 1 − j log (5.17) π 2 with Euler’s constant γ = 1.781.... The elements of In,n have the form kη 2 kγ∆Cn In,n = ∆Cn 1 − j log (5.18) 4 π 4e 44 This expression does not have a singularity, because the terms inside the log do not vanish. The total electric ﬁeld at any point outside the cylinder is obtained by plugging Jz along with (5.16) in (5.10). 5.2 The Boundary Diﬀerence Method The boundary diﬀerence method (BDM) translates a boundary value problem into boundary diﬀerence equations. As noted in Sec. 1.1, there are two methodologies developed for boundary diﬀerence equations. The ﬁrst one is the method of diﬀerence potentials (DPM), put forward by Ryabenkii [34]. To introduce the central idea of the method, consider, following Tsynkov [35], a second-order homogeneous PDE Lu = 0, valid on a domain Ω ⊂ R2 and boundary Γ = ∂Ω. From (4.6) u(r) admits the integral ∂u ∂G u(r) = [G(r, r ) (r ) − u(r ) (r, r )] dS (5.19) S ∂n ∂n For a vector density ξ Γ ≡ (ξ0 , ξ1 )|Γ on the boundary Γ, a generalized potential of the Calderon type satisﬁes the following integral by deﬁnition: ∂G PΩ · ξ Γ (r) = [ξ1 (r )G(r, r ) − ξ0 (r ) (r, r )] dS (5.20) S ∂n PΩ is the Calderon vector potential on the domain Ω. If the Dirichlet and Neumann data in (5.20) coincide with ξ0 and ξ1 , respectively, ξ Γ is said to be the trace T r of u on the boundary Γ, where T r is deﬁned as ∂u ξΓ = T r u ≡ u, |Γ (5.21) ∂n 45 Accordingly, PΩ · ξ Γ (r) = u(r). The Calderon projection PΓ is deﬁned as the trace of the potential PΩ on the boundary Γ PΓ ≡ T rPΩ ξ Γ (5.22) DPM is a discrete analog of the Calderon projections. A detailed mathematical derivation of DPM is quite involved and will not be considered here. The second methodology is called boundary algebraic equations by Martins- son and Rodin [14]. In this methodology the boundary integral equation (4.8) is replaced with a discrete convolution, leading to a system of boundary algebraic equa- tions. Martinsson considered a simple model problem: the Laplace equation with a homogeneous Dirichlet boundary conditions. In [1] Tsukerman applied a similar approach to a wave scattering problem. He utilized FLAME in formulating boundary diﬀerence equations. This section summarizes the approach proposed for BDM in [1]. 2D model problems of EM wave scattering and electrostatics are considered. 5.2.1 EM Scattering: Model Setup Figure 5.1 displays a classical 2D EM scattering problem, in which a TM polarized plane wave propagates in free space and impinges on a dielectric cylinder. The cylin- der is inﬁnitely long with an arbitrary cross-section, and zero free charges/currents on its surface. The electric ﬁeld E with a single z-component is governed by the Helmholtz equation 2 E(r) + k 2 E(r) = 0 r ≡ (x, y) (5.23) 46 Figure 5.1: Setup of an EM wave scattering problem with TM polarization. everywhere. k is the wavenumber, equals to kout and kcyl outside and inside the cylinder, respectively. On the cylinder interface the electric ﬁeld (E) satisﬁes the boundary condition E out = E in on C (5.24) on the contour boundary C (the cylinder is inﬁnitely long). Additionally, the scattered ﬁeld Esc satisﬁes the radiation condition (2.26). The incident electric ﬁeld Einc is a plane wave that has the phasor form Einc = E0 exp (−jkout · r) (5.25) for the phasor convention ejωt , where E0 is the magnitude of the incident ﬁeld. 47 On a grid, the Helmholtz equation (5.23) is approximated by diﬀerence equa- tions – for example a 5-point diﬀerence scheme S5 of the form S5 E ≡ S(−2)E(mx , my − 1) + S(−1)E(mx − 1, my ) + S(0)E(mx , my ) S(1)E(mx + 1, my ) + S(2)E(mx , my + 1) = 0 (5.26) with a double index integer m ≡ (mx , my ). The discrete boundary of the cylinder is deﬁned based on the diﬀerence scheme (5.26). If a central-stencil node m lies inside the cylinder but at least one of its neighbors (the other four nodes) is outside, it belongs to the discrete inner boundary γin . Likewise, if m lies outside and one or more of the neighbor nodes are inside, it joins the discrete outer boundary γout . Therefore, the discrete boundary is composed of two layers: γ ≡ γin ∪ γout , while the total number of boundary nodes is nγ = nγ,in + nγ,out . Figure 5.2 shows discrete boundaries for simple boundary shapes. 5.2.2 EM Scattering: Boundary Sources The boundary integral equation (4.8) is replaced with a convolution of discrete sources f and the corresponding lattice Green function g – the LGF for Helmholtz’s equation [f ∗ g](m) ≡ f (m )g(m − m ) (5.27) m∈γ The discrete convolution allows us to express the lattice based electric ﬁeld [f ∗ g(·, ·; kcyl )](m), m ∈ γin (5.28) E(m) = Einc (m) + [f ∗ g(·, ·; kout )](m) m ∈ γout (5.29) 48 (a) Circle (b) Ellipse (c) Square Figure 5.2: Discrete double-layer boundary for simple shapes: a) Circle, b) Ellipse, and c) Square. on the discrete boundary γ. According to [1], the sources f are indirectly related to the surface ﬁeld sources used in boundary integral methods, but they are ﬁctitious and need not have a physical interpretation. In contrast, the ﬁelds in (5.28) and (5.29) are physical and represent the solution on the discrete boundary. 5.2.3 EM Scattering: Boundary Diﬀerence Equations The next step in BDM is to construct boundary diﬀerence equations, from which the sources f are expressed explicitly. To do so, one uses a diﬀerence scheme S valid on the discrete boundary, and applies the boundary condition (5.24) on the ﬁelds (5.28) and (5.29), which leads to S[f ∗ g(·, ·; k(n)](n) = − SEinc,h ; n ≡ (nx , ny ) ∈ γ (5.30) 49 Einc,h is a lattice-based ﬁeld equal to Einc outside the cylinder, and zero otherwise. Equation (5.30) can be more meaningful if expressed in a matrix form Af = q (5.31) The elements of matrix A are Anm = Sn (n)g(n − m; k) (5.32) n and the column matrix q is (n,n ) qn = − Sn (n)Einc,h (5.33) n with n and m being integers that run over the discrete boundary. n is used for indexing the neighbor nodes on a stencil centered at n. The wavenumber k equals to kcyl and kout with n being inside and outside the scatterer, respectively. The term S(n) allows the selection of diﬀerent schemes at diﬀerent boundary points. The lattice-based ﬁeld can be found from (5.28) and (5.29) at grid points inside and outside the scatterer, respectively. 5.2.4 EM Scattering: Treatment of Conductors Consider the model setup in Sec. 5.2.1, but with a perfect electrically conducting (PEC) cylinder. The electric ﬁeld E is equal to zero on the cylinder, and hence, discrete sources do not exist on the discrete inner boundary γin . Therefore, the discrete boundary γ reduces to the discrete outer boundary γout , i.e. γ = γout . Based on that, one still can use (5.29) to express the electric ﬁeld on γ. Expressions (5.31) – 50 (5.33) are also still valid, but the integer n here does not run over lattice nodes that are located inside the cylinder. As a result, the diﬀerence schemes S (n) varies in sizes at diﬀerent boundary nodes. For example, Figure 5.3 shows how 5-point stencils are reduced to 3 and 4. 5.2.5 An Electrostatic Model Problem The formulations of BDM can be extended to various boundary value problems. In particular, consider the electrostatic ﬁeld problem shown in Figure 5.4. A dielectric cylinder with radius rcyl , and permittivity cyl , is placed inside a dielectric region with permittivity out . The free charges are equal to zero (i.e. ρ = 0). The cylinder is immersed in an external static uniform electric ﬁeld Es = xE0 . Figure 5.3: 3-point and 4-point stencils with a PEC cylinder. 51 Figure 5.4: Dielectric cylinder in presence of static uniform electric ﬁeld. Outside the cylinder, the electric potential Vs that corresponds to Es can be expressed as Vs (r, φ) = − E0 r cos φ + c0 (5.34) in polar coordinates (r, φ), with constant c0 set to zero for simplicity c0 = 0. In contrast with the wave case, the electric potential V is governed by the Laplace equation 2 V (r, φ) = 0 (5.35) with the boundary conditions ∂Vout ∂Vcyl Vout = Vcyl , out = cyl (5.36) ∂r ∂r at r = rcyl . 52 The procedure, introduced in the sections above, is followed to derive bound- ary diﬀerence equations for electrostatics. To express the electric potential V on a lattice via the discrete convolution (5.27), one assigns an appropriate diﬀerence scheme S (used to approximate the Laplace equation), and deﬁnes the corresponding discrete boundary γ (as in Sec. 5.2.1), leading to [f ∗ g](m), m ∈ γin (5.37) V (m) = Vs (m) + [f ∗ g](m) m ∈ γout (5.38) for m ∈ Z2 . g is the lattice Green function for the Poisson and Laplace equations. In contrast with the wave case, g here is independent of the characteristics of the medium – the permittivity ( ) in electrostatics. The sources f are expressed via the matrix form of the boundary diﬀerence equation (5.31), but with A Anm = Sn (n)g(n − m) (5.39) n and q (n,n ) qn = − Sn (n)Vs,h (5.40) n where n, m, and n are deﬁned as in equations (5.32), (5.33). The lattice based potential Vs,h equals to Vs and zero outside and inside the cylinder, respectively. Finally, the potential V is evaluated using (5.37) and (5.38). 5.3 Summary and Conclusions The boundary element method and the boundary diﬀerence method are reviewed in this chapter. In Sec. 5.1.1, a general procedure for solving diﬀerential equations with 53 BEM is reviewed. Additionally, an important example of TM wave scattering from a PEC cylinder is derived in 2D, in which piece-wise constant approximations are provided for the singular integrals. This example is used in Chapter 6 to compare the results obtained from MoM with those of the boundary diﬀerence method. There are two advantages of using the boundary diﬀerence method. The ﬁrst one is that, as in the boundary element method, BDM reduces the computational analysis to the boundaries. Secondly, the lattice Green functions utilized in BDM are nonsingular, unlike their continuous-domain counterparts used in BEM. In Sec. 5.2 it is shown how BDM can be used to solve EM wave scattering and electrostatics problems in 2D. In the wave case, model problems of TM wave scattering by dielectric and PEC cylinders are considered, and expressions for the lattice-based electric ﬁelds, involving the lattice Green function and the discrete boundary sources, are derived. One model problem is considered for electrostatics in 2D: a dielectric cylinder with circular cross-section is immersed in a uniform electric ﬁeld. In this model, similar formulations are derived but for the electric potential. Both models form the basis of the case studies in Chapter 6. 54 CHAPTER VI THE BOUNDARY DIFFERENCE METHOD: NUMERICAL EXPERIMENTS To verify the accuracy and the convergence of the boundary diﬀerence method, it is applied in six numerical experiments. The selected case studies address the model problems proposed in Chapter 5. Five of these cases deal with cylinders of smooth boundaries (no edges or corners). In these cases exact solutions exist, allowing us to evaluate the accuracy of numerical methods. One case study involves a cylinder with a rectangular cross-section, which is substantially more challenging numerically due to the corner singularities in the ﬁeld. This chapter summarizes the numerical results for these experiments, dis- cusses the respective error and ﬁeld plots, and compares them with the results ob- tained with the method of moments. 6.1 Introduction and General Assumptions Five case studies are considered for the model wave scattering problem. Two of the cases involve dielectric cylinders with either a circular or elliptic cross-section. The remaining cases deal with perfectly electrically conducting (PEC) cylinders of circular, elliptic, or rectangular cross-sections. The electrostatic case study deals with 55 a circular dielectric cylinder immersed in a uniform external ﬁeld. For convenience, the cylinders are centered at the origin in all experiments. The BDM (5.32),(5.33),(5.39),(5.40) has two ingredients: the lattice Green function g, and the boundary diﬀerence scheme S. In the experiments, we used 5- point FD schemes and the corresponding lattice Green functions, and utilized the continuous-domain Green functions to estimate the boundary conditions. 5-point FLAME schemes were used at the interface boundary. Accuracy is evaluated by the relative error ANum − AExact Relative Error = (6.1) AExact where ANum and AExact are the numerical and the exact solutions on the grid for the quantity A, respectively. In the experiments that deal with PEC cylinders the numerical results are compared with those of the moment methods. The last case study presented in this chapter does not have an exact solution, and therefore, an ”overkill” solution on the ﬁnest grid is considered in lieu of the exact one. 6.2 FLAME Basis Functions and Exact Solutions A diﬀerence scheme S employed in BDM must be accurate at the boundary. FLAME (Sec. 3.2) satisﬁes this consideration, as it utilizes accurate local Treﬀtz functions to approximate the solution (the electric ﬁeld/potential in our case), especially at the boundaries, and yet the FLAME grid can be uniform and does not have to conform (n) to the boundary. We used 5-point FLAME schemes (S5 ) with our case studies. 56 To compute the coeﬃcients of these schemes, one should seek appropriate four local Treﬀtz functions on the grid node n, for each case study. In this section FLAME bases are provided for the case studies that deal with smooth boundaries. 6.2.1 FLAME Bases for Circular Cylinders – Electrostatics In this case study, we consider an inﬁnitely long circular dielectric cylinder of radius rcyl , immersed in a uniform static electric ﬁeld. The mathematical formulations of this problem are provided in Sec. 5.2.5. Local analytical functions are available for FLAME via cylindrical harmonics. More speciﬁcally, the solution of Laplace’s equation for the electric potential inside and outside the cylinder, is used to assign (n) the FLAME basis functions ψα rl cos lφ, (n) r ≤ rcyl ψα=2l = , l = 0, 1, 2, ... (6.2) (a rl + b r−l ) cos lφ, l r ≥ rcyl l and rl sin lφ, (n) r ≤ rcyl ψα=2l−1 = , l = 1, 2, 3, ... (6.3) (a rl + b r−l ) sin lφ, l r ≥ rcyl l where rcyl is radius of the cylinder. The coeﬃcients al and bl are determined via the boundary conditions (2.23) and (2.24), they admit out + cyl al = (6.4) 2 out out − cyl 2l bl = rcyl (6.5) 2 out where out and cyl are the permittivities outside and inside the cylinder, respectively. To obtain the four basis functions, we retain the harmonic l = 0, two harmonics of orders l = 1, and one of the two harmonics of the order l = 2. 57 The exact solution of the electric potential is available for this case study. It is the solution of the Laplace equation, as it was with the FLAME basis, but with taking the applied ﬁeld Es (the static uniform electric ﬁeld used with model setup in Figure 5.4) into account. By doing so, the electric potential V is of the form −1 c1 r cos φ + Vs (r, φ), r ≥ rcyl (6.6) V (r, φ) = c2 r cos φ, r ≤ rcyl (6.7) The electric potential Vs (r, φ) is induced by the applied ﬁeld Es , and is deﬁned in (5.34). We impose the boundary conditions (2.23) and (2.24) to ﬁnd that the con- stants c1 and c2 are cyl − out 2 c1 = E0 rcyl (6.8) cyl + out −2 out c2 = E0 (6.9) cyl + out where E0 is the ﬁeld magnitude. Note that the exact solution cannot be used to assign the FLAME basis functions, as it only aﬀords one harmonic (the dipole moment). Details about this problem are available in [36]. 6.2.2 FLAME Bases for Circular Cylinders – EM Scattering In this subsection we consider FLAME bases for two case studies: scattering of an E- polarized (TM) wave from a cylinder with a circular cross-section, where the cylinder is either a dielectric or a perfect conductor. A complete formulation of these problems is given in Sec. 5.2.1. For the dielectric cylinder, the exact solution of the electric ﬁeld is available, and can be used in assigning local analytic basis functions to generate the FLAME 58 scheme. To do so, we recall a wave transformation identity known as the Jacobi- Anger expansion, to express the incident ﬁeld Einc (deﬁned in (5.25)) via cylindrical z harmonics. This identity represents the solution of the 2D Helmholtz equation in polar coordinates. Accordingly, Einc admits z ∞ Einc = E0 z j −l Jl (kout r) exp(jlφ) (6.10) l=−∞ where kout is the wavenumber outside the cylinder. Jl is the Bessel function of the ﬁrst kind, and the order l. In the vicinity of the cylinder, the electric ﬁeld Ez admits ∞ E0 j l al Jl (kcyl r) exp(jlφ), r ≤ rcyl (6.11) l=−∞ Ez = ∞ (2) E0 j −l {Jl (kout r) + bl Hl (kout r)} exp(jlφ), r > rcyl (6.12) l=−∞ where al and bl are coeﬃcients that can be determined via the interface conditions (2) (2.17) and (2.18), kcyl is the wavenumber inside the cylinder, Hl is the Hankel function of the second kind and the order l, and rcyl is the cylinder radius. The ﬁeld inside the cylinder can be expanded into series (6.11), and the ﬁeld outside into (6.12). The coeﬃcients al and bl have the form √ √ out Jl (kcyl rcyl )Jl (kout rcyl ) − cyl Jl (kcyl rcyl )Jl (kout rcyl ) al = √ (2) √ (2) (6.13) cyl Jl (kcyl rcyl )Hl (kout rcyl ) − out Jl (kcyl rcyl )Hl (kout rcyl ) √ (2) √ (2) out Jl (kout rcyl )Hl (kout rcyl ) − out Jl (kout rcyl )Hl (kout rcyl ) bl = √ (2) √ (2) (6.14) cyl Jl (kcyl rcyl )Hl (kout rcyl ) − out Jl (kcyl rcyl )Hl (kout rcyl ) where out and cyl are the permittivities outside and inside the cylinder, respectively. The subscript ( ) refers to the ﬁrst derivative with respect to the radial variable r. For our second case study, the scatterer is a perfect conductor. Therefore, the ﬁeld inside the cylinder is zero. Outside the cylinder, the electric is still expressed 59 by equation (6.12), but with the constant bl Jl (kout rcyl ) bl = − (2) (6.15) Hl (kout rcyl ) Detailed derivations of these expressions are available in [37],[38]. For the two case studies, the four FLAME basis functions are assigned from the corresponding ﬁeld summations. One should choose the most signiﬁcant terms, and these are the ones of the orders l = 0, −1, 1 (i.e. monopole, and two dipoles), and one of the quadrupoles (l = −2, 2). 6.2.3 Bases for Elliptic Cylinders – EM Scattering As in the subsection above, consider the model problem in Sec. 5.2.1, but now with an elliptic cylinder that is either a dielectric or a perfect conductor. In the Cartesian coordinates (x, y), an ellipse centered at the origin with the semi-major axis a on the x-axis, and the semi-minor axis b on the y-axis, is deﬁned by the equation (x/a)2 + (y/b)2 = 1, a > b (6.16) For this ellipse, the semi-focal distance f has the form √ f= a2 − b 2 (6.17) Diﬃculties arise if one tries to ﬁnd the exact solution for the electric ﬁeld in the polar or the Cartesian coordinates. Instead, an elliptic coordinates system (u, v) can be used, so that u represents the unit normal vector on the ellipse curve, and v is the tangential one. The new coordinate system is found mapping the xy-plane on a 60 complex one C2 , with the mapping function z x − iy z= , − ∞ < x, y < ∞ (6.18) f z ∈ C2 . Accordingly, the elliptic coordinates (u,v) are functions of the complex variable z, and expressed by the relations √ u = |Im{−i log(z + z 2 − 1)}|, z ∈ C2 (6.19) √ Re[−i log(z + z 2 − 1)], Im(z) ≤ 0 (6.20) v= √ 2π − Re[−i log(z + z 2 − 1)], Im(z) > 0 (6.21) for u and v, respectively, where u ∈ [0, ∞], and v ∈ [0, 2π]. Re and Im refer to the real and the imaginary parts, respectively. Figure 6.1 shows the new coordinates system. Note that u is constant on each elliptic curve. Thus, it can be used to determine the two sides of the cylinder boundary. e.g. if the boundary is located at u0 , the inside and outside regions of the cylinder are u < u0 and u > u0 , respectively, and u0 is found by setting z = a (the major semi-axis) in equation (6.19). For the new coordinates system, the incident Electric ﬁeld Einc (deﬁned in z (5.25)) can be expanded into a series ∞ √ Einc z = 8πE0 im Sm (qout , v)Jm (qout , u)/Nm (qout ) e e e (6.22) m=0 e where the subscript (e) refers to even, and Nm (qout ) is a normalization factor, with 2 e e qout = kout f 2 /4. Sm (qout , v) and Jm (qout , u) are the m-th order even angular and radial Mathieu functions of the ﬁrst kind, respectively. Mathieu functions are elliptic harmonics; they are the solutions for the Helmholtz equation in elliptic coordinates 61 Figure 6.1: An elliptic coordinate system. Reprinted from [3]. (Details of the functions are provided in [39],[40],[41]). The normalization factor e Nm (qout ) can be found by applying the integral 2π e Nm (q) = e [Sm (qout , v)]2 dv (6.23) 0 For convenience, the subscript (e) will be neglected in the further discussion. In the vicinity of the dielectric cylinder, the total electric ﬁeld Ez is ∞ inc √ E + 8πE z 0 im am Sm (qout , v)Hm (qout , u), u > ucyl (6.24) m=0 Ez = ∞ √ 8πE0 im bm Sm (qcyl , v)Jm (qcyl , u), u ≤ ucyl (6.25) m=0 where am and bm are coeﬃcients to be determined, and Hm (qout , u) is the radial Mathieu function of the fourth kind. Elliptic and cylindrical harmonics are similar; Radial Mathieu functions of the ﬁrst and the fourth kinds are analogous to the Bessel 62 functions for the same kinds, while the angular Mathieu’s play the same role as the trigonometric functions in circular cylinders. If the interface conditions (2.17) and (2.18) are met at the boundary of the cylinder, the mathematical expressions for am and bm are Jm (qout , u0 ) Jm (qout , u0 ) − Hm (qout , u0 ) am = (6.26) γ(m)Jm (qcyl , u0 ) Jm (qcyl , u0 ) − Hm (qout , u0 ) Jm (qout , u0 ) Jm (qout , u0 ) − Jm (qcyl , u0 ) bm = (6.27) Nm (qout )Hm (qout , u0 ) Jm (qcyl , u0 ) − Hm (qout , u0 ) where J = J /J, and H = H /H. The derivatives denoted by the prime are taken with respect to u. The new coeﬃcient γ(m) appeared as a result of applying the interface conditions. It is found by evaluating the integration 2π γ(m) = Sm (qout , v)Sm (qcyl , v)dv (6.28) 0 Therefore, unlike the trigonometric harmonics exp(jlφ)(in (6.11) and (6.12)), the angular functions Sm (qout , v) and Sm (qcyl , v) do not form an orthogonal set. For this reason, the series terms in (6.25) and (6.24) cannot be used to assign the FLAME basis funtions in this case study, as they are linearly dependent for the coeﬃcients (6.26) and (6.26), but still they are important, as they handle the exact solution. A possible way to generate a FLAME scheme for this case study is by using an osculating circle that approximates the ellipse boundary, while still using expressions (6.11) and (6.12) as basis functions. The osculating circle in this approach is deﬁned by its radius R and center (xosc , yosc ) at each point on the ellipse boundary. Consider the parametrization of the ellipse x = a cos t and y = b sin t; the osculating circle is 63 deﬁned as ab R= (6.29) (b cos t)2 + (a sin t)2 for the radius, and cos t xosc = a cos t − {(a sin t)2 + (b sin t)2 } (6.30) a sin t yosc = b sin t − {(a sin t)2 + (b sin t)2 } (6.31) b the center. As in the circular cylinders to obtain the FLAME basis functions in this case, we retain the most signiﬁcant terms from the series (6.11) and (6.12). For the second case study here, the scatterer is a perfect conductor. The electric ﬁeld inside is zero; therefore, only equation (6.24) holds. but the coeﬃcients am are now Jm (qout , u0 ) am = − (6.32) Nm (qout )Hm (qout , u0 ) For these coeﬃcients the series terms in (6.24) become linearly independent. There- fore, they can be used to generate FLAME basis functions, by retaining the ﬁrst four terms of the ﬁeld summation. Details about these scattering problems are available in [39],[40]. 6.3 Numerical Results for Smooth Boundaries In this section, numerical results are provided for ﬁve problems in terms of the relative error, as well as plots of solutions obtained from the BDM. 64 6.3.1 Circular Cylinder in the Presence of Static Uniform Field As an electrostatic example, the problem described in Sec. 5.2.5 is solved by BDM with the physical parameters: in = 5, out = 1, E0 = 1, and rcyl = 0.855. The electric ﬁeld vector has a single component in the positive x-direction. The boundary scheme S is a 5-point FLAME (deﬁned in the Sec 6.2.1). Figure 6.2 shows the solution obtained by BDM of the electric potential inside and outside the cylinder at grid size h = 0.01. The relative error computed over all grid nodes is plotted in Figure 6.3 as a function of the grid size. From the ﬁgure it is clear that the convergence rate is of order 1.95, with two sources of error: the boundary diﬀerence scheme (5-point FLAME scheme), and the discrete Green function (5-point FD scheme). Both schemes are of order 2, and both lead to almost the same rate of convergence in BDM. 6.3.2 Electromagnetic Scattering by Circular Cylinders Two experiments are performed for scattering by circular cylinders; the ﬁrst one deals with a dielectric scatterer, while the second one addresses a perfectly electrically conducting (PEC) cylinder. The FLAME bases that are used with the boundary schemes for the two cases are mentioned in Sec. 6.2.2. For the dielectric case, the scatterer is of radius rcyl = 0.8075, with the wave numbers kcyl = 2.96 inside the cylinder, and kout = 2.1 everywhere outside. The incident wave is assumed to be TM polarized with magnitude E0 = 1; it propagates 65 Figure 6.2: Dielectric cylinder in the presence of a static uniform ﬁeld: the electro- static potential with grid size h = 0.01. Figure 6.3: Dielectric cylinder in the presence of a static uniform ﬁeld: relative error of the electric potential vs. grid size h in BDM. 66 in the positive x-direction. Figure 6.4 represents the solution for the real part of the electric ﬁeld obtained by BDM at grid size h = 0.01. As a veriﬁcation, the solutions extracted from BDM at diﬀerent grid sizes are compared with the quasi-exact one, where the expansion representation of the exact solution in equations (6.11) and (6.12) is computed up to order 50. Figure 6.5 shows the relative error of the method vs. grid size (h). From the ﬁgure, the convergence rate in this case is 2, for the same reasons as in the electrostatic example (FLAME and FD schemes are of order 2). For the same physical parameters as above, but with a PEC scatterer, the real part of the electric ﬁeld is represented in Figure 6.6 at grid size h = 0.01, where the ﬁeld inside the scatterer is zero due to the conducting material. It may look Figure 6.4: Scattering by a circular dielectric cylinder: the real part of the electric ﬁeld with grid size h = 0.01. 67 Figure 6.5: Scattering by a circular dielectric cylinder: relative error of the electric ﬁeld vs. grid size h in BDM. strange at ﬁrst glance that a few of the ﬁeld contour lines “pierce” the boundary of the conductor (e.g. at point P indicated in the ﬁgure). However, the picture is correct. The behavior of the ﬁeld in the vicinity of point P outside the conductor can be qualitatively described as E ∼ x y in the coordinate system (x , y ) centered at P. Indeed, the product x y satisﬁes the Laplace equation and hence approximately satisﬁes the Helmholtz equation for long wavelengths. Figure 6.7 shows two curves for relative error; one of them is obtained by BDM, and the second is for the method of moments MoM. In BDM the boundary sources were estimated on the grid as in Figure 5.2, but for MoM they were assigned on the contour of the scatterer. Additionally, a piecewise constant approximation for the boundary sources (ﬁlament currents) was used in MoM just as in the scattering model in Sec. 5.1.2. Here the asymptotic regime of the relative error tends to 68 Figure 6.6: Scattering by a circular PEC cylinder: the real part of the electric ﬁeld with grid size h = 0.01. The special behavior of the ﬁeld in the vicinity of point P is noted in the text. Figure 6.7: Scattering by a circular PEC cylinder: relative error of the electric ﬁeld vs. grid size h in BDM and MoM. 69 move towards ﬁner grid sizes, if compared with the dielectric case; i.e., the error plot includes additional points in the range [0.005 − 0.01]. The order of convergence in BDM is still around 2; it is 2.11. In comparison, BDM gives better results than MoM, as the convergence rate for the latter is of order 0.73. As in the above case the expansion of the quasi-exact solution was computed to order 50. 6.3.3 Electromagnetic Scattering by Elliptic Cylinders The numerical experiments addressing EM scattering by smooth boundaries are re- peated, but with an elliptic cylinder. Consider the same electrical characteristics as in the subsection above. The ellipse is assumed to have a major semi-axis a = 0.8075, and a minor semi-axis b = 0.6075. As above, the incident plane wave is TM polar- ized, and propagates in the positive x-direction. FLAME basis functions and the exact solution are already discussed and provided in Sec. 6.2.3. The ﬁrst 10 terms in the series representation in equations (6.25) and (6.24) are taken as the quasi-exact solution. Mathieu functions and their related coeﬃcients are computed in MATLAB with the help of a toolbox provided by Cojocaru [41]. For the dielectric scatterer, Figure 6.8 is a plot of the real part of the electric ﬁeld at grid size h = 0.01. By comparing this plot with the one obtained for the circular dielectric (Figure 6.4), it is obvious that the total ﬁeld throughout the region behaves in a similar manner; however, for the elliptic scatterer, it is stretched out horizontally and squeezed in vertically due to the geometry of the ellipse. In Figure 6.9 the relative error of BDM is plotted vs. the grid size h. The order of convergence 70 Figure 6.8: Scattering by an elliptic dielectric cylinder: the real part of the electric ﬁeld with grid size h = 0.01. Figure 6.9: Scattering by an elliptic dielectric cylinder: relative error of the electric ﬁeld vs. grid size h in BDM. 71 is 2.19, which is still around 2 as in the previous cases. The error plot shows some noise added to the curve, as a result of using the osculating circle (deﬁned in Sec. 6.2.3) in approximating the elliptic boundary; however, this noise does not aﬀect the convergence rate. For the same conﬁguration of the example above, if the dielectric cylinder is replaced by a PEC one, the real part of the electric ﬁeld extracted from BDM at grid size h = 0.01 can be viewed as in Figure 6.10. As in the dielectric case, if the ﬁeld around this scatterer in now compared with one for the PEC circular cylinder (in Figure 6.6), it is clear how it behaves in a similar way, but being stretched in (with respect to the y-axis) and out (to the x-axis) due to the diﬀerence in the geometry of the scatterers. From Figure 6.11, the relative errors for BDM and MoM are of orders 2.11 and 0.85, respectively. As in circular cylinders, BDM shows better rate of conver- gence than in MoM. Still, it is important to remember that a piecewise constant approximation (deﬁned in Sec. 5.1.2) is used with MoM. 6.4 Implementation of BDM with Non-Smooth Boundaries In this section one case study is considered, as an example of boundaries with corners. The results obtained with BDM are displayed in terms of the electric ﬁeld and the convergence rate. An excellent agreement of the results with the ones from MoM is demonstrated. 72 Figure 6.10: Scattering by an elliptic PEC cylinder: the real part of the electric ﬁeld with grid size h = 0.01. Note the special behavior of the electric ﬁeld in the vicinity of point P; it is explained in Sec. 6.3.2. Figure 6.11: Scattering by an elliptic PEC cylinder: relative error of the electric ﬁeld vs. grid size h in BDM and MoM. 73 6.4.1 FLAME Basis Functions Before proceeding with FLAME basis functions, consider the model problem in Figure 6.12. A rectangular cylinder is illuminated by a TM polarized plane wave. For Figure 6.12: Scattering by a rectangular PEC cylinder: problem setup. simplicity, the height and the width of the cylinder are the same (a square). No closed-form solution is available to generate a Treﬀtz basis for FLAME. However, local approximations of the ﬁeld at the corners can be used instead. More speciﬁcally, Figure 6.13 represents a model problem for EM scattering by a PEC wedge. For an arbitrary value of the angle α the total ﬁeld outside the wedge has the form ∞ E(r, φ) = 4νE0 im Jmν (kout r) sin(mνφinc ) sin(mνφ) (6.33) m=0 where φinc is the the angle of incidence measured with respect to the line of symmetry in the wedge. ν is constant value and deﬁned as ν = π/(2π −α). Naturally, expansion (6.33) can be written in a local coordinate system centered at any corner, to generate FLAME bases near that corner. For the PEC rectangular scatterer deﬁned in Figure 74 Figure 6.13: Scattering by a rectangular PEC cylinder: wedge problem. 6.12 only three basis functions are needed at the boundary nodes, since only 4-point stencils exist due to the conducting material of the cylinder (see Figures 5.2c and 5.3). The basis functions are found by retaining the ﬁrst three terms of the expansion in (6.33). 6.4.2 Numerical Results Consider the problem setup displayed in Figure 6.12. Again the incident TM polarized wave travels in the positive x-axis, with the wavenumber kout = 2.1. The cross-section of the PEC cylinder is assumed to be a square with side length L = 1.6. The FLAME basis functions are produced as discussed in the subsection above. Figure 6.14 shows the real part of the electric ﬁeld distribution as a function of coordinates, where the ﬁeld was obtained by BDM at grid size h = 0.01. As a veriﬁcation, Figure 6.15 shows the ﬁeld on a square contour whose sides are three grid layers away from the outer boundary of the scatterer. The orientation 75 Figure 6.14: Scattering by a rectangular PEC cylinder: the real part of the electric ﬁeld with grid size h = 0.01. The special behavior in the vicinity of point P outside the cylinder is discussed in Sec. 6.3.2. of the contour (a → b → c → d → a) is deﬁned with respect to the corners as in Figure 6.12. The ﬁeld is computed at grid size h = 0.01 by two methods: BDM, and MoM. It is clear that the results are in perfect agreement. Finally, the convergence rate is examined in this case. From the circular and elliptic case studies the asymptotic rate at wavenumber kout = 2.1 tends to manifest itself at small grid sizes (in the range h = [0.005 − 0.04]), but here the quasi-exact solution is assumed at h = 0.005 (the smallest grid size available in the experiments due to some constraints in the resources). This diﬃculty can be avoided by decreasing the wavenumber or the radian frequency of the incident ﬁeld, where the convergence rate is evident at much larger grid sizes. Figures 6.16 and 6.17 are plots of the 76 Figure 6.15: Scattering by a rectangular PEC cylinder: the electric ﬁeld at a square located three grid layers away from the cylinder boundary by BDM and MoM, with grid size h = 0.01. relative error vs. grid size at two diﬀerent wavenumbers kout = 0.1 and kout = 0.5, respectively. It is obvious that the order of convergence in BDM, which is 1.56, is smaller than the order of convergence for the case studies with smooth boundaries. The degradation occurred due to the singularities of the ﬁeld at the corners. 6.5 Conclusions The accuracy and convergence of BDM are demonstrated via six case studies con- sidered in this chapter. Five of the numerical experiments address problems of TM wave scattering from cylinders that are either a dielectric or a perfect conductor, 77 Figure 6.16: Scattering by a rectangular PEC cylinder: convergence rate of the elec- tric ﬁeld vs. grid size h in BDM, with wavenumber kout = 0.1. Figure 6.17: Scattering by a rectangular PEC cylinder: convergence rate of the elec- tric ﬁeld vs. grid size h in BDM, with wavenumber kout = 0.5. 78 with circular, elliptic, or rectangular cross-sections. The sixth case study considers an electrostatic problem in which a dielectric cylinder with a circular cross-section is immersed in a static uniform ﬁeld. 2D plots for the electric potential and the real part of the electric ﬁeld are provided for the electrostatic and the wave scattering case studies at grid size h = 0.01, respectively. For the case studies of smooth boundaries, the Treﬀtz bases of FLAME boundary schemes along with the exact solutions are provided in Sec. 6.2. The convergence rate in these cases is of order 2. In the PEC cylinders the convergence rates are compared with those of MoM. It is evident that BDM gives better results, as the simplest version of MoM is used. In the case study that deals with a non-smooth boundary, the FLAME basis functions of the boundary scheme are provided in Sec. 6.4.1, and the convergence rate is of order 1.56, as a result of the ﬁeld singularities at the corners. In this case, an additional plot for the electric ﬁeld shows the agreement of the results obtained from BDM and MoM. 79 CHAPTER VII CONCLUSIONS Conclusions and possible future work are summarized in this chapter. 7.1 Conclusions and Summary The boundary diﬀerence method is explained and implemented in three model prob- lems of EM scattering and electrostatics. It is shown how BDM avoids the singular Green functions, inherent in traditional boundary element methods, by replacing them with the lattice Green functions that are nonsingular. In addition to this, simi- larly to boundary element methods, it is shown how BDM treats unbounded domains naturally, as the unknowns are deﬁned on the discrete boundaries. The validity and the accuracy of the boundary diﬀerence method are demon- strated for six numerical experiments. The case studies considered two types of boundaries. The ﬁrst one is smooth boundaries in which all the results show that the convergence rate is of order ∼ 2. Secondly, a boundary with four corners (a square) is considered, where the rate of convergence is of order 1.56, as a result of the ﬁeld singularities at the corners. All experiments show that the convergence rates in BDM are commensurate with that of the schemes chosen for the lattice Green functions and for the boundary diﬀerence equations. To increase the order of convergence in 80 this method, both schemes are to be replaced with some higher order schemes, e.g. the 5-point FD and FLAME schemes can be replaced by the corresponding 9-point schemes. The formulations derived for the boundary diﬀerence method in Chapter 5 can be extended for 3D problems in a straightforward way. However, developing com- puter algorithms for vectorial 3D problems is not trivial. Additionally, fast multipole methods (FMM) should be applied to speed up the time of the 3D matrix operations. To our knowledge, this is the ﬁrst time (in addition to [1]) that boundary- diﬀerence methods have been applied to electromagnetic scattering 1 . 7.2 Future Work 1. The method can be extended to scatterers of arbitrary shapes, while still using local analytic Treﬀtz bases to generate FLAME boundary schemes. This can be accomplished by putting together “libraries” of FLAME basis functions for diﬀerent geometric entities: at ﬂat surfaces, near corners, for osculating circle approximations, etc. 2. The boundary diﬀerence method can be extended to electromagnetic problems in three dimensions, which would be of substantial interest. 3. The method can be generalized to solve other classes of linear problems, e.g. heat and elasticity. 1 See also [35] 81 REFERENCES [1] I. Tsukerman. A Singularity-Free Boundary Equation Method for Wave Scat- tering. IEEE Trans. Antennas and Propagation, 59(2):555–562, 2011. [2] P.A. Martin. Discrete Scattering Theory: Green’s Function for a Square Lattice. Wave Motion, 43(7):619–629, 2006. [3] E. Weisstein. Elliptic Cylindrical Coordinates. From MathWorld – A Wolfram Web Resource. http://mathworld.wolfram.com/EllipticCylindricalCoordinates.html. [4] J. Volakis and K. Sertel. Integral Equation Methods for Electromagnetics. SkiTech Publishing, Inc., 2012. [5] W. Chew, M. Tong, and B. Hu. Integral Equation Methods for Electromagnetic and Elastic Waves. Morgan and Claypool Publishers, 2008. [6] W. Gibson. The Method of Moments in Electromagnetics. Chapman and Hall/CRC, 2007. [7] MB. Friedman and R. Shaw. Diﬀraction of Pulse by Cylindrical Obstacles of Arbitrary Cross Section. J. Appl. Mech. Trans. ASME, 29:40–6, 1962. [8] RP. Banaugh and W. Goldsmith. Diﬀraction of Steady Acoustic Waves by Sur- faces of Arbitrary Shape. J. Acoust. Soc. Am., 35:1590–601, 1963. [9] R. Harrington. Field Computation by Moment Methods. Krieger Publishing Co., Inc., Melbourne, FL, USA, 1st edition, 1982. [10] R. Harrington. Matrix Methods for Field Problems. Proceedings of the IEEE, 55(2):136–149, 1967. [11] A. Cheng and D. Cheng. Heritage and Early History of the Boundary Element Method. Engineering Analysis with Boundary Elements, 29:268–302, 2005. 82 [12] A. Reznik. Approximation of Surface Potentials of Elliptic Operators by Diﬀer- ence Potentials. Dokl. Nauk SSSR, 263:1318–1321, 1982. [13] S. Tsynkov. On the Deﬁnition of Surface Potentials for Finite-Diﬀerence Oper- ators. J. Sci. Comput., 18:155–189, April 2003. [14] P.G. Martinsson and G. Rodin. Asymptotic Expansions of Lattice Greens Func- tions. 458:2609–2622, 2002. [15] C. Saltzer. Discrete Potential Theory for Two-Dimensional Laplace and Poisson Diﬀerence Equations. Technical Report 4086, National Advisory Committee on Aeronautics, 1958. [16] P.G. Martinsson and G. Rodin. Boundary Algebraic Equations for Lattice Prob- lems. 465(2108):2489–2503, 2009. [17] P.G. Martinsson. Fast Multiscale Methods for Lattice Equations. PhD thesis, The University of Texas at Austin, 2002. [18] J. Jackson. Classical Electrodynamics. Wiley, third edition, August 1998. [19] J. Stratton. Electromagnetic Theory. McGraw-Hill: New York, 1941. [20] N. Ida. Engineering Electromagnetics. Springer-Verlag, 2004. [21] I. Tsukerman. Computational Methods for Nanoscale Applications: Particles, Plasmons, and Waves. Springer, 2008. [22] I. Tsukerman. Flexible Local Approximation Method for Electro- and Magneto- statics. IEEE Trans. Magnetic, 40:941–944, 2004. [23] I. Tsukerman. Electromagnetic Applications of a New Finite-Diﬀerence Calculus. IEEE Trans. Magnetic, 41:2206–2225, 2005. [24] I. Tsukerman. A Class of Diﬀerence Schemes with Flexible Local Approximation. J. Comput. Phys., 211:659–699, January 2006. [25] I. Tsukerman. Quasi-Homogeneous Backward-Wave Plasmonic Structures: The- ory and Accurate Simulation. J. Opt. A: Pure Appl. Opt., 11:114025, 2009. [26] I. Tsukerman. Treﬀtz Diﬀerence Schemes on Irregular Stencils. J. Comput. Phys., 229:2948–2963, 2010. 83 [27] G. Green. An Essay on the Application of Mathematical Analysis to the Theories of Electricity and Magnetism. Nottingham, 1828. [28] Q.H. Qin. Green’s Function and Boundary Elements of Multiﬁeld Materials. Elsevier Ltd., 2007. [29] G. Barton. Elements of Green’s Functions and Propagation: Potentials, Diﬀu- sion, and Waves (Oxford Science Publications). Oxford University Press, USA, July 1989. [30] D. Duﬀy. Green’s Functions with Applications. Chapman and Hall/CRC, 2001. [31] E. Economou. Green’s Functions in Quantum Physics. Springer, Heidelberg, 2006. [32] S. Katsura, T. Morita, S. Inawashiro, T. Horiguchi, and Y. Abe. Lattice Green’s Function. Introduction. J. Math. Phys., 12(5):892–895, 1971. [33] S. Katsura and S. Inawashiro. Lattice Green’s Functions for the Rectangular and the Square Lattices at Arbitrary Points. J. Math. Phys., 12(8):1622–1630, 1971. [34] V. Ryaben”kii. Method of Diﬀerence Potentials and Its Applications. Springer, 1st edition, 2001. [35] M. Medvinsky, S. Tsynkov, and E. Turkel. The Method of Diﬀerence Potentials for the Helmholtz Equation Using Compact High Order Schemes. submitted for publication to the J. Sci. Comput. [36] K. Roy. Potential Theory in Applied Geophysics. Springer, 2008. [37] R. Harrington. Time-Harmonic Electromagnetic Fields. Wiley-IEEE Press, 2001. [38] J.M. Jin. Theory and Computation of Electromagnetic Fields. Wiley-IEEE Press, 2010. [39] J. Van Bladel. Electromagnetic Fields. Wiley-IEEE Press, 2007. [40] E. Cojocaru. Mathieu Functions Approach to Bidimensional Scattering by Di- electric Elliptical Cylinders. e-print arXiv0808.2123v1. 84 [41] E. Cojocaru. Matlab free available computer code Mathieu Functions Toolbox v. 1.0. 85