; Singularity-Free Boundary Methods for Electrostatics and Wave Scattering
Documents
Resources
Learning Center
Upload
Plans & pricing Sign in
Sign Out
Your Federal Quarterly Tax Payments are due April 15th Get Help Now >>

Singularity-Free Boundary Methods for Electrostatics and Wave Scattering

VIEWS: 11 PAGES: 96

  • pg 1
									SINGULARITY-FREE BOUNDARY METHODS FOR ELECTROSTATICS AND

                       WAVE SCATTERING




                              A Thesis

                            Presented to

           The Graduate Faculty of The University of Akron




                        In Partial Fulfillment

                  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 suffer from the singularity of the Green function

kernels. The boundary difference method (BDM) implemented and applied in this

thesis alleviates this drawback by replacing the singular kernels with Green’s functions

defined 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 finite difference and

finite element analysis, can be handled in a computationally efficient way using Fast

Multipole acceleration or other existing techniques, with relatively minor modifica-

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 effectiveness of BDM for elec-

trostatic and wave scattering problems, including field 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 Difference 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 Difference 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     Definitions 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 field 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 field 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 field. . . . . . .      52

6.1      An elliptic coordinate system. Reprinted from [3]. . . . . . . . . . . . .        62




                                             ix
6.2   Dielectric cylinder in the presence of a static uniform field: the
      electrostatic potential with grid size h = 0.01. . . . . . . . . . . . . . .     66
6.3   Dielectric cylinder in the presence of a static uniform field: 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 field with grid size h = 0.01. . . . . . . . . . . . . . . . . . . . .   67
6.5   Scattering by a circular dielectric cylinder: relative error of the
      electric field vs. grid size h in BDM. . . . . . . . . . . . . . . . . . . .      68
6.6   Scattering by a circular PEC cylinder: the real part of the electric
      field with grid size h = 0.01. The special behavior of the field 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
      field vs. grid size h in BDM and MoM. . . . . . . . . . . . . . . . . . .         69
6.8   Scattering by an elliptic dielectric cylinder: the real part of the
      electric field with grid size h = 0.01. . . . . . . . . . . . . . . . . . . . .   71
6.9   Scattering by an elliptic dielectric cylinder: relative error of the
      electric field vs. grid size h in BDM. . . . . . . . . . . . . . . . . . . .      71
6.10 Scattering by an elliptic PEC cylinder: the real part of the electric
     field with grid size h = 0.01. Note the special behavior of the electric
     field 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
     field 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 field 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 field 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 field 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 field 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 differential equations in the volume

or, alternatively for linear piecewise-homogeneous media, via surface integral equa-

tions. The first approach, when implemented numerically, gives rise to finite element

(FE) or finite difference (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 first 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 finite

difference method (FD) and the finite element method (FEM) – that are used more

extensively. The unknowns are typically defined in the volume for FEM and FD

and on surfaces for BEM. This difference 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 infinite 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.

        Confining the computations to the boundaries is a desirable feature, but the

singular kernels are often problematic. Therefore, efforts have been made to develop

new methods that enjoy the advantages of BEM, but do not include singular terms.

                                          2
Boundary difference equations (BDM) are an interesting alternative from this per-

spective. Two different methodologies exist for BDM. The first one is known as the

method of difference 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 define the difference 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 artificial 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 difference method extending the ap-
                                         3
proach of [1]. For EM scattering, different problems are considered: scattering by

an infinitely 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 field singularities at the corners. An electrostatic example is also

considered: a dielectric cylinder is immersed in a static field. In all cases, field plots

are examined and convergence rates analyzed.



1.3 Thesis Organization


After the introduction and overview in the first 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 field, 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 fields, waves and the Helmholtz equation.


2.1.1 Maxwell’s Equations

Electromagnetic fields are mathematically represented by four field vectors: electric

and magnetic field intensities E, H, electric displacement and magnetic field densities




                                         6
D, B. The differential 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 fields. In a linear isotropic homogeneous medium

the relations between the field 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 fields. At a fixed 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 field (e.g. E(t)) can be found from its phasor as


                                     E(t) = Re{Eejωt }                          (2.11)


From now on, all electromagnetic fields 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 fields

can be algebraically eliminated, leading to a second-order equation for the remaining

field. 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 defined 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 find an expression that involves H only.
                                               8
2.1.3 The Poynting Vector

The Poynting vector specifies power flow 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 fields.


2.1.4 Interface of Materials

EM scattering occurs when an electromagnetic wave propagates in a medium and

hits an obstacle with different material parameters. In the presence of two different

dielectric media with permittivities    1   and   2,   the electric fields 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 field at the interface, respectively. Note that (2.17)

predicts a discontinuity in the electric field. Similar expressions are available to

describe the magnetic field 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 field, respectively.


2.1.5 TE and TM Polarizations

The term polarization specifies the direction of the electric field in certain regions

of space. At material interfaces, the plane of incidence is defined 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 field vector is parallel to the plane of incidence,

then the wave is called transverse-magnetic (TM). Otherwise, it is transverse-electric

(TE). Other definitions for polarizations are not considered, as they are not used in

this work.



2.2 Electrostatic Formulations


The electric field 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 field is zero in this case, and the electric and magnetic fields 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 field 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 fields 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 first one is the Dirichlet boundary condition, in which the

value of the function to be solved for is defined on the boundary. The second one

is the Neumann boundary condition, where the normal derivative of the function is

defined on the boundary. Finally, the Robin boundary condition involves a linear

combination of the first two conditions.
                                           11
2.3.1 The Sommerfeld Radiation Condition

Wave problems require different types of boundary conditions, an example of which

is the Sommerfeld radiation condition. The Helmholtz equation defined in (2.14)

describes the behavior of the time-harmonic electric field in a source-free region. If

this field 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 find

that

                                   2
                                       E + k 2 E = jµωJ                         (2.25)

In general the field 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 field 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 fields (E, H, D, B) and sources (J, ρ). From
                                            12
Maxwell’s equations, we derived the Helmholtz equation for the electric field 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 field 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 finite difference

method and the flexible local approximation method. Additionally, it provides a very

brief summary of the finite element method. The boundary equation methods are

introduced in a separate chapter.



3.1 The Finite Difference method


In finite difference (FD) analysis a differential equation is discretized and solved on a

grid. Even though the finite 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 Difference 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) defined 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 finite amount of information can be processed (in FD, solution on a finite 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 difference 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 field, 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 difference scheme as a

5-point scheme. A similar approach is followed to derive a difference 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 finite difference 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. Refining 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 flexible approximation,

as discussed in the following section.



3.2 The Flexible Local Approximation Method


The flexible 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 difference 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 differential equation.

        In this thesis FLAME is used with the boundary difference method BDM in

last two chapters.




                                         17
3.2.1 Trefftz 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 Trefftz schemes for homogeneous

equations, consider, following Tsukerman [21], a homogeneous differential equation


                                     Lu = 0 in Ω(i)                                    (3.8)


with a differential 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 coefficients cα

                                              m(i)
                                    (i)               (i) (i)
                                   uh     =          cα ψα                           (3.10)
                                              α=1

                                                                (i)
For a M -point stencil the coefficient 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 difference 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 coefficient 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) find 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) find the scheme coefficients 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 Trefftz 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 differential 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 finite element method (FEM) is applied to a variational problem, typically ob-

tained from the differential 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, infinite 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 first one is the finite differ-

ence method. In particular, 5-point difference 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 flexible local approximation method (FLAME), reviewed in the


                                         21
chapter for both homogeneous and inhomogeneous (nonzero right hand side) equa-

tions. Local Trefftz 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 difference schemes in the last two chapters. A brief description of the finite

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 find 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 differential form into an

integral one, with the Green function as the kernel.


4.1.1 General Review

The Green function is mathematically defined as the solution of the equation


                              LG(r, r ) = − δ(r − r )                            (4.1)


where L represents a linear differential operator, G is the Green function, and δ is

the Dirac delta function. Equation (4.1) is considered in a finite or infinite 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
field 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 first 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 field theory in

general and in the analysis of EM wave scattering in particular. Green’s function for

the Helmholtz problem satisfies the following equation by definition:


                              2
                                  G(r, r ) + k 2 G(r, r ) = − δ(r − r )                (4.4)


where k is the wave number. G satisfies the Sommerfeld radiation condition (Sec.

2.3.1). We are not defining the boundary conditions yet, so G is not uniquely defined.

                                                 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) simplifies if G is translationally invariant, i.e.


                                    G(r, r ) = G(r − r )                                      (4.9)


                                                  26
Translational invariance holds if the differential operator has constant coefficients,

i.e. if the medium is homogeneous throughout. Solutions for Green’s functions are

provided in many textbooks for different 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

different, 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

difference 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 field approximation using two different methods; the

first method is due to Koster, while the second one uses an idea of Lifshitz. Finally,

[1] employs the finite difference (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 find the LGF by the finite difference method.

        In this thesis, the LGF g (lower case g is used to differentiate LGF from the

continuous Green function G) is found via the FD method. The LGF by definition

satisfies:

                                  Lg(n) = − δ(n)                                (4.12)

with a double index integer n ≡ (n1 , n2 ). In general L is any difference operator

approximating the Helmholtz operator. The boundary conditions can be approxi-

mated by the continuous Green function just as in [1]. In the far field g approaches

the continuous one, and therefore, the LGF can be accurately computed in a finite

size spatial window with G imposed as a boundary condition on that window. An


                                         29
alternative somewhat different 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 finite 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 find 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 difficulty 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 fields. 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 defined 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: Definitions 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 field

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

field 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 figures

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 finding the LGF associated with the

discrete Poisson equation in 2D: the Discrete-Time Fourier Transform (DTFT) anal-

ysis, and the Finite Difference 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 defining the boundary conditions just as we did in the Helmholtz

case.

          To apply the Fourier transform method consider equation (4.12), with the

five-point difference 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)

simplifies 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 defined 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 defined 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 specifically,

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 difference methods, as the one implemented in this thesis. There-

fore, in this chapter both the continuous-domain and the lattice Green functions

have been defined 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 different 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 confined 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 differential 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 defined on Γ:


                                   u =        αn ψn                             (5.1)
                                          n


where αn are some coefficients. Naturally, in numerical simulation the series (in

general infinite) is truncated to a finite 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 affected 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 infinitely long cylinder that is a perfect conductor with

an arbitrary cross section. We need to solve for the electric field outside the cylinder

(inside the cylinder the electric field is zero).

        Since the cylinder is infinitely 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 field Einc is
                                          z



                              Einc = E0 exp(−jkr cos φ)
                               z                                                  (5.8)
                                           42
in polar coordinates, with the field magnitude E0 , and the wavenumber k. The total

electric field on the boundary C admits


                            Ez = Einc + Esc = 0,
                                  z      z                     on C                      (5.9)


where Esc refers to the scattered field.
       z


          A filament current density J, defined on a contour boundary L, radiates the

electric field 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 field induces a current density on the contour of the

scatterer. Accordingly, a scattered field is produced due to the current. Based on

that, by using equations (5.9) and (5.10), the incident field 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 field 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 satisfied 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 filament 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 field at any point outside the cylinder is obtained by

plugging Jz along with (5.16) in (5.10).



5.2 The Boundary Difference Method


The boundary difference method (BDM) translates a boundary value problem into

boundary difference equations. As noted in Sec. 1.1, there are two methodologies

developed for boundary difference equations. The first one is the method of difference

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 satisfies the following integral by definition:

                                                                     ∂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 defined as

                                                           ∂u
                                  ξΓ = T r u ≡        u,        |Γ                    (5.21)
                                                           ∂n


                                                45
Accordingly, PΩ · ξ Γ (r) = u(r). The Calderon projection PΓ is defined 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

difference 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 infinitely long with an arbitrary cross-section, and zero free charges/currents

on its surface. The electric field 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 field (E) satisfies the

boundary condition

                                E out = E in     on C                           (5.24)

on the contour boundary C (the cylinder is infinitely long). Additionally, the scattered

field Esc satisfies the radiation condition (2.26). The incident electric field 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 field.




                                         47
        On a grid, the Helmholtz equation (5.23) is approximated by difference equa-

tions – for example a 5-point difference 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 defined based on the difference

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 field
                        
                         [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 field sources used in boundary integral methods, but they are fictitious

and need not have a physical interpretation. In contrast, the fields in (5.28) and

(5.29) are physical and represent the solution on the discrete boundary.


5.2.3 EM Scattering: Boundary Difference Equations

The next step in BDM is to construct boundary difference equations, from which the

sources f are expressed explicitly. To do so, one uses a difference scheme S valid on

the discrete boundary, and applies the boundary condition (5.24) on the fields (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 field 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 different schemes at different boundary points.

        The lattice-based field 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 field 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 field 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 difference schemes S (n) varies in sizes

at different 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 field 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 field 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 field.




         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 difference equations for electrostatics. To express the electric potential V on

a lattice via the discrete convolution (5.27), one assigns an appropriate difference

scheme S (used to approximate the Laplace equation), and defines 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 difference 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 defined 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 difference method are reviewed in

this chapter. In Sec. 5.1.1, a general procedure for solving differential 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 difference method.

        There are two advantages of using the boundary difference method. The first

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 fields,

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 field. 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 difference 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 field.

        This chapter summarizes the numerical results for these experiments, dis-

cusses the respective error and field 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 field. 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 difference 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 finest grid is considered in lieu of the exact one.



6.2 FLAME Basis Functions and Exact Solutions


A difference scheme S employed in BDM must be accurate at the boundary. FLAME

(Sec. 3.2) satisfies this consideration, as it utilizes accurate local Trefftz functions to

approximate the solution (the electric field/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 coefficients of these schemes, one should seek appropriate four local

Trefftz 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 infinitely long circular dielectric cylinder of radius

rcyl , immersed in a uniform static electric field. The mathematical formulations of

this problem are provided in Sec. 5.2.5. Local analytical functions are available

for FLAME via cylindrical harmonics. More specifically, 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 coefficients 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 field Es (the static uniform electric field 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 field Es , and is defined in

(5.34). We impose the boundary conditions (2.23) and (2.24) to find 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 field magnitude. Note that the exact solution cannot be used to assign

the FLAME basis functions, as it only affords 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 field 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 field Einc (defined 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

first kind, and the order l. In the vicinity of the cylinder, the electric field 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 coefficients 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

field inside the cylinder can be expanded into series (6.11), and the field outside into

(6.12). The coefficients 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 first derivative with respect to the radial variable r.

        For our second case study, the scatterer is a perfect conductor. Therefore,

the field 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 field summations. One should choose the most significant 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 defined 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)


Difficulties arise if one tries to find the exact solution for the electric field 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 field Einc (defined 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 first 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 field 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 coefficients 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 first 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 coefficient γ(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 coefficients

(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 defined

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
defined 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 significant terms from the series (6.11) and (6.12).

         For the second case study here, the scatterer is a perfect conductor. The

electric field inside is zero; therefore, only equation (6.24) holds. but the coefficients

am are now
                                              Jm (qout , u0 )
                              am = −                                              (6.32)
                                          Nm (qout )Hm (qout , u0 )

For these coefficients the series terms in (6.24) become linearly independent. There-

fore, they can be used to generate FLAME basis functions, by retaining the first four

terms of the field 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 five 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

field vector has a single component in the positive x-direction. The boundary scheme

S is a 5-point FLAME (defined 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 figure it is clear that the convergence rate is of order

1.95, with two sources of error: the boundary difference 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 first 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 field: the electro-

static potential with grid size h = 0.01.




Figure 6.3: Dielectric cylinder in the presence of a static uniform field: 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 field obtained by BDM at grid size h = 0.01.

        As a verification, the solutions extracted from BDM at different 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 figure, 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 field is represented in Figure 6.6 at grid size h = 0.01, where

the field 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

field with grid size h = 0.01.




                                          67
Figure 6.5: Scattering by a circular dielectric cylinder: relative error of the electric

field vs. grid size h in BDM.




strange at first glance that a few of the field contour lines “pierce” the boundary

of the conductor (e.g. at point P indicated in the figure). However, the picture is

correct. The behavior of the field 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 satisfies the Laplace equation and hence approximately

satisfies 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 (filament 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 field

with grid size h = 0.01. The special behavior of the field 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 field

vs. grid size h in BDM and MoM.




                                         69
move towards finer 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 first 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 coefficients 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

field 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 field 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

field with grid size h = 0.01.




Figure 6.9: Scattering by an elliptic dielectric cylinder: relative error of the electric

field 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 (defined in Sec.

6.2.3) in approximating the elliptic boundary; however, this noise does not affect the

convergence rate.

        For the same configuration of the example above, if the dielectric cylinder is

replaced by a PEC one, the real part of the electric field extracted from BDM at grid

size h = 0.01 can be viewed as in Figure 6.10. As in the dielectric case, if the field

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 difference 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 (defined 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 field 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 field

with grid size h = 0.01. Note the special behavior of the electric field 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 field

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 Trefftz basis for FLAME. However,

local approximations of the field at the corners can be used instead. More specifically,

Figure 6.13 represents a model problem for EM scattering by a PEC wedge. For an

arbitrary value of the angle α the total field 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 defined 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 defined 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 first 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 field distribution as a function of coordinates, where the

field was obtained by BDM at grid size h = 0.01.

        As a verification, Figure 6.15 shows the field 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

field 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 defined with respect to the corners as in

Figure 6.12. The field 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 difficulty can be avoided by decreasing

the wavenumber or the radian frequency of the incident field, 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 field 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 different 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 field 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 field 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 field 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 field. 2D plots for the electric potential and the real

part of the electric field 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 Trefftz 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 field singularities at the corners. In this case, an additional plot for the

electric field 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 difference 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 defined on the discrete boundaries.

        The validity and the accuracy of the boundary difference method are demon-

strated for six numerical experiments. The case studies considered two types of

boundaries. The first 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 field

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 difference 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 difference 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 first time (in addition to [1]) that boundary-

difference 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 Trefftz bases to generate FLAME boundary schemes. This can

        be accomplished by putting together “libraries” of FLAME basis functions for

        different geometric entities: at flat surfaces, near corners, for osculating circle

        approximations, etc.


  2. The boundary difference 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. Diffraction of Pulse by Cylindrical Obstacles of
     Arbitrary Cross Section. J. Appl. Mech. Trans. ASME, 29:40–6, 1962.

 [8] RP. Banaugh and W. Goldsmith. Diffraction 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 Differ-
     ence Potentials. Dokl. Nauk SSSR, 263:1318–1321, 1982.

[13] S. Tsynkov. On the Definition of Surface Potentials for Finite-Difference 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
     Difference 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-Difference Calculus.
     IEEE Trans. Magnetic, 41:2206–2225, 2005.

[24] I. Tsukerman. A Class of Difference 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. Trefftz Difference 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 Multifield Materials.
     Elsevier Ltd., 2007.

[29] G. Barton. Elements of Green’s Functions and Propagation: Potentials, Diffu-
     sion, and Waves (Oxford Science Publications). Oxford University Press, USA,
     July 1989.

[30] D. Duffy. 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 Difference Potentials and Its Applications. Springer,
     1st edition, 2001.

[35] M. Medvinsky, S. Tsynkov, and E. Turkel. The Method of Difference 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

								
To top