A ROBUST STABILIZATION METHODOLOGY FOR TIME DOMAIN INTEGRAL
EQUATIONS IN ELECTROMAGNETICS
By
Andrew J. Pray

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Electrical Engineering - Doctor of Philosophy
2014

ABSTRACT
A ROBUST STABILIZATION METHODOLOGY FOR TIME DOMAIN INTEGRAL
EQUATIONS IN ELECTROMAGNETICS
By
Andrew J. Pray
Time domain integral equations (TDIEs) are an attractive framework from which to analyze
electromagnetic scattering problems. Casting problems in the time domain enables study of
systems with nonlinearities, characterization of transient behavior both at the early and late
time, and broadband analysis within a single simulation. Integral equation frameworks have
the advantages of restricting the computational domain to the scatterer surface (boundary
integral equations) or volume (volume integral equations), implicitly satisfying the radiation
boundary condition, and being free of numerical dispersion error. Despite these advantages,
TDIE solvers are not widely used by computational practitioners; principally because TDIE
solutions are susceptible to late-time instability. While a plethora of stabilization schemes
have been developed, particularly since the early 1980s, most of these schemes either do
not guarantee stability, are difficult to implement, or are impractical for certain problems.
The most promising methods seem to be the space-time Galerkin schemes. These are very
challenging to implement as they require the accurate evaluation of 4-dimensional spatial
integrals. The most successful recent approach to implementing these schemes has been to
approximate a subset of these integrals, and evaluate the remaining integrals analytically.
This approach describes the quasi-exact integration methods [1, 2]. The method of [1] approximates 2 of the 4 dimensions using numerical quadrature. The remaining integrals are
evaluated analytically by determining shadow boundaries on the domain of integration. In
[2], only 1 dimension is approximated, but the procedure also relies on analytical integration

between shadow boundaries. These two characteristics-the need to find shadow boundaries
and develop analytical integration rules-prevent these methods from being extended to higher
order tessellations of scattering surfaces. This is an important restriction as the use of curvilinear elements can greatly improve the accuracy of the geometric representation. The need
for a method to accurately evaluate the spatial integrals involved in these formulations on
higher order surface tessellations motivates this thesis.
The major novelty of this thesis is a space-time separated expansion of the convolution
with the retarded potential Green’s function. This separation leads to integrands that are
smooth over the entire domain of integration. This implies that integration can be accurately
carried out via numerical quadrature (not analytically) and shadow boundaries do not need
to be found, unlike in the quasi-exact integration methods. The numerical nature of the
method allows it to trivially be implemented on higher order surface descriptions.
In this thesis, we will detail the procedure of the separable expansion and investigate,
both numerically and analytically, the error incurred in truncating the expansion to a given
upper limit. We will validate the stability of the resulting scheme by (1) observing the late
time behavior of solutions to scattering from a variety of objects and (2) deriving and implementing an eigenvalue analysis to demonstrate the absence of growing terms. Additionally,
this thesis will detail the use of the separable expansion in tandem with the plane wave time
domain algorithm to accelerate the solution. Also, we will present extension of the spacetime separation to analysis of penetrable materials using the PMCHWT formulation. A
prescription for integrating singular and near-singular kernels over curved elements will also
be given. The final contribution of this thesis is the application of the space-time separation
to the generalized method of moments (GMM) with smooth surface parameterization.

Dedicated to my family and to the memory of Joseph Pray.

iv

ACKNOWLEDGMENTS

To begin, I must give my sincerest thanks to my primary advisor, Shanker Balasubramaniam.
Perhaps the most important thing I learned from him over these years is perseverance. Over
the past five years there have been innumerable instances where an approach to a project has
hit a wall, but his perseverance and optimism in trying new approaches has been infectious.
Outside of the research arena, his door has always been open and he has been a great mentor.
He has helped me in so many ways that it would be senseless to try and list them here.
Over the first several months and beyond, I had the great fortune of being mentored by
another graduate student, Naveen Nair. This early period was invaluable to me in adjusting
my work practices and mindset from undergraduate to graduate study. In addition to being
an excellent mentor, he was and has continued to be a first-rate collaborator on a variety
of projects. Much of the work presented in this thesis would not be possible without his
contributions.
Next, I would like to thank the entire Electromagnetics Research Group for helping make
an environment so conducive to quality work. The professors have always had an open door
and their insights have always been very helpful. I would particularly like to thank Professor
Rothwell for drawing my attention to this research group when I was an undergraduate in
one of his courses. I’ve also benefited from being surrounded by such bright, motivated, and
affable students. In particular, the many discussions I’ve had with Andrew Baczewski and
Dan Dault have always been lively, be they research-related or on other topics. I’m very
fortunate to have had the two of them as lab-mates.
Lastly, I must acknowledge the support I’ve received from my family over the years.
v

None of my achievements would be remotely possible without their unwavering support.
They have always had the most devout belief in my ability to succeed, even when I perhaps
wasn’t giving them reason to do so. In particular I cannot give enough praise to my mother,
Marcia McCay, for the energy she has expended and the sacrifices she has made on my
behalf. I’m not sure how she managed all of this as a single mother, but her strength is a
constant inspiration. It is to my family that I dedicate this thesis.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . .
LIST OF FIGURES

x

. . . . . . . . . . . . . . . . . . . . . . . . . . .

xi

KEY TO SYMBOLS AND ABBREVIATIONS . . . . . . . . . . . . . .

xv

Chapter 1 Introduction
1.1 Advantages of Time Domain Methods . . .
1.2 Advantages of Integral Equation Methods
1.3 Drawbacks and Challenges of TDIE solvers
1.4 Instability of MOT solutions to TDIEs . .
1.5 Stabilization Schemes for TDIEs . . . . . .

.
.
.
.
.

1
2
2
3
5
6

Chapter 2 Mathematical Preliminaries
2.1 Maxwell’s Equations and Wave Equations . . . . . . . . . . . . . . . . . . .
2.2 Mixed Potential Solutions for Integral Equations . . . . . . . . . . . . . . . .
2.3 Electric and Magnetic Field Integral Equations . . . . . . . . . . . . . . . . .

11
11
13
16

Chapter 3 Method of Moments Solution of TDIEs
3.1 Collocation in time . . . . . . . . . . . . . . . . . .
3.2 Higher order Space-time Galerkin Discretization . .
3.2.1 Higher order temporal basis . . . . . . . . .
3.2.2 Interpolation properties of higher order basis

21
21
24
25
28

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.

Chapter 4 Separable Expansion Method
4.1 Accurate evaluation of retarded potential integrals . . . . . . . . . . . . . . .
4.2 Quasi-Exact Integration Schemes . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Integration on Curvilinear Elements . . . . . . . . . . . . . . . . . . . . . . .
4.4 Separable Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Truncation Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6 Incorporation into Method of Moments Scheme . . . . . . . . . . . . . . . .
4.6.1 Efficient evaluation of scalar potential . . . . . . . . . . . . . . . . . .
4.7 Modification to Higher order space-time Galerkin scheme . . . . . . . . . . .
4.7.1 Interpolation properties . . . . . . . . . . . . . . . . . . . . . . . . .
4.7.2 Higher order space-time Galerkin discretization with separable expansion
vii

31
31
32
37
38
40
46
48
50
50
54

Chapter 5 Numerical Results
5.1 Eigenvalue stability analysis . . . . . . . . . . . .
5.2 Low Order Discretizations . . . . . . . . . . . . .
5.2.1 Sphere . . . . . . . . . . . . . . . . . . . .
5.2.2 Parallel plates . . . . . . . . . . . . . . . .
5.2.3 Thin box . . . . . . . . . . . . . . . . . . .
5.2.4 Rectangular Cavity . . . . . . . . . . . . .
5.3 Higher Order Discretizations . . . . . . . . . . . .
5.3.1 Convergence in far field scattering . . . . .
5.3.2 Stability analysis for thin box . . . . . . .
5.3.3 Sphere with various spatio-temporal bases
5.3.4 Cone-sphere . . . . . . . . . . . . . . . . .
5.3.5 Almond . . . . . . . . . . . . . . . . . . .
5.3.6 VFY218 . . . . . . . . . . . . . . . . . . .
5.4 Summary . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

55
56
62
63
64
65
67
69
70
72
74
75
76
78
79

Chapter 6 Acceleration using Multilevel Plane Wave Time Domain Algorithm
81
6.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.2 Multilevel PWTD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.2.2 PWTD Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.2.3 Undifferentiated TD-CFIE using Runge Kutta Integration . . . . . . 90
6.2.4 Formulation using separable expansion in higher order space-time Galerkin
framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.3.1 Stability and accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.3.2 Timing results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.4 More Challenging Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Chapter 7 Implementation for
7.1 Motivation . . . . . . . . .
7.2 Formulation . . . . . . . .
7.3 Discretization . . . . . . .
7.4 Separable Expansion . . .
7.5 Low Frequency Instability
7.6 Numerical Results . . . . .
7.7 Discussion . . . . . . . . .

Penetrable Scatterers
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

105
105
106
109
111
112
112
115

Chapter 8 Generalized Method of Moments with Separable Expansion Method117
8.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
8.1.1 Generalized Method of Moments (GMM) . . . . . . . . . . . . . . . . 118
8.1.2 Locally smooth surface approximations . . . . . . . . . . . . . . . . . 118
viii

8.2

8.3
8.4
8.5
8.6

8.1.3 Time Domain GMM with Separable Expansion
Methodology . . . . . . . . . . . . . . . . . . . . . . .
8.2.1 Locally smooth surface approximation . . . . .
8.2.2 Definition of basis functions . . . . . . . . . . .
GMM System . . . . . . . . . . . . . . . . . . . . . . .
Accurate Integration via Separable Expansion . . . . .
Numerical Results . . . . . . . . . . . . . . . . . . . . .
TD-EFIE instability . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

Chapter 9 Conclusions and Future Work
9.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.1.1 PWTD acceleration of TD-EFIE with separable expansion
9.1.2 GMM discretization of TD-EFIE . . . . . . . . . . . . . .
9.1.3 Nonlinear solvers . . . . . . . . . . . . . . . . . . . . . . .
9.1.4 Antenna problems . . . . . . . . . . . . . . . . . . . . . .
APPENDIX

.
.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

BIBLIOGRAPHY

. . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

.
.
.
.
.
.
.
.

119
120
120
122
124
125
126
127

.
.
.
.
.

129
131
131
132
132
133
135

146

LIST OF TABLES

Table 4.1 Lower bound of M needed to reproduce interpolation accuracies (undifferentiated/1st derivative/integral) . . . . . . . . . . . . . . . . . . . . . .

53

Table 5.1

Rectangular cavity resonant frequencies (all frequencies in MHz) . . .

70

Table 6.1

Input parameters for plates . . . . . . . . . . . . . . . . . . . . . . . .

97

Table 6.2

Input parameters for spheres . . . . . . . . . . . . . . . . . . . . . . .

98

Table A.1 Comparison of efficiency of different rules . . . . . . . . . . . . . . . . 143
Table A.2 Number of points required for convergence of integration rules on higher
order triangles for various test point locations . . . . . . . . . . . . . . . . . 144

x

LIST OF FIGURES

Figure 2.1

General description of a PEC scattering problem . . . . . . . . . . . .

13

Figure 2.2

General domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

Figure 2.3

Scattering from PEC object . . . . . . . . . . . . . . . . . . . . . . .

18

Figure 3.1

Interpolation errors for higher order basis . . . . . . . . . . . . . . . .

30

Figure 4.1 1st order shifted Lagrange polynomial: (a) undifferentiated, (b) integral, (c) 1st derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

Figure 4.2

Regions of smoothness on source element (graphic from [1]) . . . . . .

35

Figure 4.3 Various intersections between spheres and source triangles (graphic
from [1]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

Figure 4.4

Spatial support of expansion (4.8) . . . . . . . . . . . . . . . . . . . .

39

Figure 4.5

Domain of the spatial integrals

. . . . . . . . . . . . . . . . . . . . .

41

Figure 4.6

Convergence in M

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

Figure 4.7

Spectra of approximating signals

. . . . . . . . . . . . . . . . . . . .

47

Figure 4.8

Interpolation errors for various upper limits of (4.8) . . . . . . . . . .

52

Figure 4.9

Interpolation errors for separable expansion . . . . . . . . . . . . . .

53

Figure 5.1

Current on sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

Figure 5.2

Eigenvalue analysis of sphere . . . . . . . . . . . . . . . . . . . . . . .

65

Figure 5.3

RCS of sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

Figure 5.4

Current on parallel plates . . . . . . . . . . . . . . . . . . . . . . . .

66

xi

Figure 5.5

Current on box . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

Figure 5.6

Rectangular cavity . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

Figure 5.7

Current on rectangular cavity . . . . . . . . . . . . . . . . . . . . . .

68

Figure 5.8

Modes of rectangular cavity . . . . . . . . . . . . . . . . . . . . . . .

69

Figure 5.9

Far field convergence for plate . . . . . . . . . . . . . . . . . . . . . .

71

Figure 5.10 Far field convergence for sphere . . . . . . . . . . . . . . . . . . . . .

72

Figure 5.11 Eigenvalue analysis of thin box discretized with (a) 1st order basis/delta testing, (b) 0th order basis/Galerkin testing, (c) 1st order basis/Galerkin
testing, (d) 2nd order basis/Galerkin testing . . . . . . . . . . . . . . . . . . 73
Figure 5.12 Current on sphere using various spatio-temporal basis functions . . .

75

Figure 5.13 Current on cone-sphere . . . . . . . . . . . . . . . . . . . . . . . . . .

76

Figure 5.14 RCS of cone-sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

Figure 5.15 Current on almond . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

Figure 5.16 RCS of almond . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

Figure 5.17 Current on VFY218 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

Figure 5.18 RCS of VFY218 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

Figure 6.1

Source and observer plane wave expansions . . . . . . . . . . . . . . .

85

Figure 6.2

Decomposition of 2 dimensional scattering geometry . . . . . . . . . .

88

Figure 6.3

Interpolation errors for projections onto Prolates (δ = 0) . . . . . . .

93

Figure 6.4

Interpolation errors for projections onto Prolates (δ = 0.5) . . . . . .

94

Figure 6.5

Current on plate . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

xii

Figure 6.6

Current on sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

Figure 6.7

RCS of sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

Figure 6.8

Timings of PWTD-accelerated plate simulations . . . . . . . . . . . .

98

Figure 6.9

Timings of PWTD-accelerated sphere simulations . . . . . . . . . . .

99

Figure 6.10 Current on cone-sphere compared to direct solution (p = 2) . . . . . . 100
Figure 6.11 RCS of cone-sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Figure 6.12 Current on almond compared to direct solution (p = 2) . . . . . . . . 101
Figure 6.13 RCS of almond . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
Figure 6.14 Spectrum of scattered field . . . . . . . . . . . . . . . . . . . . . . . . 103
Figure 6.15 Grouping of dipoles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Figure 7.1 Pictorial representation of scatterer with multiple homogeneous regions
(graphic from [1]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
Figure 7.2

Current on sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

Figure 7.3

Current on composite sphere compared against exact integration . . . 114

Figure 7.4

Current on thin box compared against exact integration . . . . . . . . 115

Figure 7.5

Current on cone-sphere compared against exact integration . . . . . . 116

Figure 8.1

Overlapping GMM Patches (Ωk ) . . . . . . . . . . . . . . . . . . . . 121

Figure 8.2

Patch normal and projection planes Γk

Figure 8.3

Partition of unity functions: (a) 1D, (b) 2D . . . . . . . . . . . . . . 123

Figure 8.4

Current on nacelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

Figure 8.5

Current on cone-sphere . . . . . . . . . . . . . . . . . . . . . . . . . . 128
xiii

. . . . . . . . . . . . . . . . 121

Figure A.1 Projection onto a curved surface from a projection plane . . . . . . . 137
Figure A.2 Tangent triangle for curvilinear, mapped triangles. . . . . . . . . . . . 139
Figure A.3 Tangent triangle for curvilinear, projection surfaces. . . . . . . . . . . 140
Figure A.4 Subdivision of source domain for observation points projecting outside 141
Figure A.5 Integration limits for one sub-triangle

xiv

. . . . . . . . . . . . . . . . . 142

KEY TO SYMBOLS AND ABBREVIATIONS

• ACA = Adaptive Cross Approximation
• ACE = Accelerated Cartesian Expansions
• AIM = Adaptive Integral Method
• BIE = Boundary Integral Equation
• CFIE = Combined Field Integral Equation
• CFL = Courant-Friedrich-Levy
• DE = Differential Equation
• EFIE = Electric Field Integral Equation
• EM = Electromagnetics
• FDIE = Frequency Domain Integral Equation
• FDTD = Finite Difference Time Domain
• FEM = Finite Element Method
• FFT = Fast Fourier Transform
• FMM = Fast Multipole Method
xv

• GMM = Generalized Method of Moments
• GWP = Graglia Wilton Peterson
• IE = Integral Equation
• IFFT = Inverse Fast Fourier Transform
• MFIE = Magnetic Field Integral Equation
• MLFMA = Multilevel Fast Multipole Algorithm
• MoM = Method of Moments
• MOT = Marching on in Time
• PEC = Perfect Electric Conductor
• PMCHWT = Poggio Miller Chew Harrington Wu Tsai
• PU = Partition of Unity
• PWTD = Plane Wave Time Domain
• RCS = Radar Cross Section
• RK = Runge-Kutta
• RWG = Rao Wilton Glisson
• SIE = Surface Integral Equation
• SVD = Singular Value Decomposition
• TDIE = Time Domain Integral Equation
xvi

• TD-FEM = Time Domain Finite Element Method
• VIE = Volume Integral Equation

xvii

Chapter 1

Introduction

The reliable analysis of EM systems via solutions to TDIEs has remained an elusive goal
since the 1960s [3]. The major culprit for this unreliability has been the instability of MOT
solutions to TDIEs. The issue of instability has garnered much attention, which has led
to a plethora of stabilization schemes, to be discussed in more detail later in this Chapter.
Despite the prevalence of such schemes, the most reliable TDIE solvers are still, for the most
part, inferior in terms of efficiency and applicability, to their frequency domain- or DE-based
counterparts in one way or another. This statement will be justified in the ensuing sections.
In this thesis, we seek to further bridge this gap. The 3 major goals of this thesis are to
develop a method with (1) applicability to higher order geometric discretizations (2) higher
order accuracy in the temporal dimension and (3) robust stability properties. We also aim
to achieve a fourth goal, i.e., to accelerate the method to extend its applicability to practical
applications. To begin, we motivate this work by (1) discussing the desirable properties
of MOT based TDIE analysis and (2) presenting the major challenges to performing such
analysis reliably. To this end, we now examine the benefits of time domain methods and IE
methods, separately. We will then look at the specific drawbacks and challenges of TDIE
solvers.
1

1.1

Advantages of Time Domain Methods

Time domain analysis presents a number of benefits over frequency domain analysis.
First, nonlinear lumped circuit models are easily coupled to EM systems in the time domain.
Such analysis in the frequency domain is much less straightforward. A second advantage is
that, in the time domain, broadband analysis can be performed within a single simulation
run. In the frequency domain, simulations must be run at each frequency of interest to
replicate the data. Therefore, a TDIE solution can capture finely sampled data across a
broad frequency band. Lastly, if the transient response of a system is desired, then the
analysis must be performed in the time domain. Frequency domain analysis is restricted to
steady state phenomena.

1.2

Advantages of Integral Equation Methods

IE frameworks for EM analysis have a number of attractive qualities. First, for BIEs, only
the boundary surrounding a homogeneous domain needs to be discretized and for VIEs only
the volume of a scatterer is discretized. For DEs, not only the scatterer must be discretized,
but the medium surrounding the scatterer needs to be discretized up to the boundary where
the computational domain is truncated. Consequently, significantly fewer unknowns are
typically needed for a given problem using IEs. This leads us to a second advantage of
IE methods. While the radiation boundary condition at infinity is implicitly satisfied by IE
formulations, for DE formulations the computational domain must be truncated in some way.
Typically, the domain is extended multiple wavelengths past the scatterer and truncated with
an artificial absorbing boundary layer. These are local operators. Alternatively, conformal
2

surfaces may be used to create a DE-BIE formulation, where the BIE is used to satisfy the
radiation boundary condition [4]. The former increases the number of unknowns and leads
to higher numerical dispersion. The latter is a commonly used and powerful method, but
has the same complexity as a standard BIE. A third benefit of IEs is that they do not suffer
from numerical dispersion error. Each point in space directly interacts with every other
point through the Green’s function in contrast to the local, nearest neighbor interactions,
characteristic of DE formulations.

1.3

Drawbacks and Challenges of TDIE solvers

Although TDIE methods enjoy the aforementioned benefits, they also suffer from a number of challenges. First, a single TDIE simulation is computationally expensive. The global
interactions of IEs, while removing dispersion error, yield system matrices with O((Nsie )2 )
non-zero entries. DE solvers are characterized by local interactions, which leads to O(Nsde )
non-zero entries (keeping in mind that for many problems Nsie << Nsde ). Additionally, the
cost of one TDIE simulation is more expensive than that of an FDIE solution. The computational cost of one TDIE solution scales as O(Nt Ns2 ) compared to O(Ns2 ) for FDIEs. This
of course ignores the fact that TDIE solutions capture data across a broad frequency range
while frequency domain simulations compute data only at one frequency. Along these lines,
it should be noted that acceleration techniques are well developed for TDIEs [5, 6] as well as
for FDIEs [7, 8]. An implementation of one such accelerator with the method presented in
this thesis is covered in Chapter 6.
A second challenge to TDIE solvers is proper temporal approximation of the surface
current. Two ways in which higher order accuracy in time can be achieved are (1) the sur3

face current can either be represented using higher order temporal basis functions [1, 9, 10]
or (2) a higher order time stepping scheme can be introduced, via passage to the Laplace
and Z transform domain [11]. The former is a challenge because a higher order temporal
basis/testing scheme with each belonging to the correct function spaces can lead to complications (we will return to this later). The latter is challenging as it leads to infinite tailed
interactions due to the passage to the Laplace domain, which increases the computational
complexity unless this tail is truncated. One last difficulty with the discretization of the time
dimension is the presence of a temporal derivative and integral in the TD-EFIE. The temporal integral in scalar potential computations yields infinite tailed interactions, although this
can be remedied, as will be shown later in this thesis. Many formulations instead discretize
the time derivative of the TD-EFIE, removing the temporal integral from the scalar potential
term, but this results in two temporal derivatives on the vector potential. For insufficiently
smooth temporal basis sets, this will introduce delta functions at the interpolation points,
which are typically ignored in implementation. The effects of ignoring these delta functions
are not understood. The presence of the second derivative will also limit the accuracy of the
solution.

The most important disadvantage of TDIEs is that the stability of time domain DE
solvers is better understood and, while CFL conditions are well established for both FDTD
and TD-FEM, an analogous condition for TDIEs does not exist [12]. As stability is the
biggest challenge and the main focus of this thesis, we will devote the next section to the
issue in more detail.
4

1.4

Instability of MOT solutions to TDIEs

As mentioned in section 1.3, the major challenge in developing a robust TDIE solution
methodology is instability. It is therefore prudent to define what is meant by instability.
First, instability is mainly an issue for MOT solutions of TDIEs. In an MOT solution, rather
than inverting an Ns Nt × Ns Nt matrix which contains interactions between all elements on
the scatterer for the entire time history, current coefficients at a given time step are found
by solving an Ns × Ns system. These coefficients are used to update the right hand side of
the equation at subsequent time steps. Thus, an Ns × Ns system (same matrix, different
right hand side) is solved at Nt time steps to compute the entire time history of the current.
This procedure can cause errors to exponentially build as the marching system continues
and an unstable solution can result which grows without bound long after the excitation has
vanished.
There are two classes of instability that should be distinguished. The first category
is instability due to the scattered field operators themselves, pre-discretization. The time
domain EFIE has a null space corresponding to DC solenoidal currents. If this is improperly
handled, the solution will show a term which grows linearly in time [13]. Additionally, both
the time domain EFIE and MFIE can suffer from interior resonance problems for closed
scatterers. Spurious, non-radiating currents can appear in the solution which correspond
to the actual currents appearing in the interior problem. Analysis of these operators in the
Laplace domain shows that poles corresponding to these resonances lie on the imaginary axis
and should therefore not cause instability [14]. However, it has been shown that discretization
errors can push these poles to the right half plane [15], resulting in instability. Both of these
sources of instability have well developed remedies [13, 15–17].
5

The second category of instability is that due to errors in the discretization, e.g., insufficient accuracy in spatial quadrature or poor mesh quality. The greatest challenge in
addressing this type of instability is the type of temporal representation of the current and
the type of testing in time that is used. Conventional temporal basis functions are not
smooth and therefore difficult to integrate accurately. Smooth temporal basis functions have
been presented in the literature [9], but these have drawbacks, which will be discussed in the
next section. Inaccuracy in the spatial integration of non-smooth temporal basis functions
is a major focus of this thesis.

1.5

Stabilization Schemes for TDIEs

Many remedies to instability have been developed over the past 30 years or so. Among
the earliest of these efforts were the studies by Tijhuis, first for the one dimensional case [18]
and later for two dimensions [19]. These works investigated the accumulation of errors in
marching solutions and concluded that inaccuracies in the computation of source integrals
were responsible for instability. While stability criteria were developed, they were for specific
geometries and only ensured the instability appeared “late enough” in the solution, rather
than yielding a truly stable solution. Later work consisted of methods predicated upon the
removal of high frequency content from MOT solution. These included the work of Rynne
[20] in which he attempted to establish conditions which, when met, ensured stability. Lower
order temporal finite differencing was found to be more robust, severely limiting the accuracy
of the solution. Averaging schemes were later developed in order to suppress the high frequency oscillations blamed for instability [21, 22]. These approaches, while attractive in their
simplicity and often effective, were not robust and failed in certain cases. Another method
6

[23] incorporated low pass finite impulse response (FIR) filters. The authors discussed the
exponential growth term in unstable solutions and showed that when this term was below
unity, the error in the current decayed. While the authors were thorough and demonstrated
numerically that the FIR filters attenuated the growth term below unity, the analysis was a
posteriori, i.e., there was no way of knowing prior to simulation how many FIR terms were
required to ensure stability for a given problem.

Other methods attempted to remove high frequency oscillations (or rather not introduce
them) through the temporal basis set. One approach was (quasi-) band-limited basis functions [9]. This method used approximate prolate spheroidal wave functions for the temporal
basis. These functions were carefully tailored to be approximately band-limited as well as
approximately time-limited. The major drawback of this method was the relatively large,
symmetric support in time of these functions which created a discretely non-causal system.
This mandated an extrapolation procedure so that an MOT solution could be obtained.
Good results were demonstrated, but the method was found to be unstable when using the
EFIE for closed structures and the extrapolation procedure placed restrictions on the time
step size that could be used. Another approach was the use of quadratic B-spline basis
functions [10]. While these functions could be precisely integrated, they tended to shift the
energy in the marching system matrix off the diagonal, which could be another source of
instability. Additionally, the 2nd derivative which appears in the time-differentiated TDEFIE is discontinuous for orders of 2 or less, which makes precise integration much more
complicated. One method has been presented using splines, which can be shown to be
equivalent to a space-time Galerkin discretization [24], i.e., the marching systems in the two
approaches are identical. As a result, it inherits the robust stabilization characteristics of
7

the full space-time Galerkin framework. This is specifically for quadratic splines, and it has
not been demonstrated how this can be extended to arbitrary order.
While each of these methods produced stable results for a number of scattering geometries, none were developed within a framework from which stability should be expected for
arbitrary scatterers. In [25] Ha-Duong stated with regard to the engineering literature on
TDIE stabilization schemes that “...one notes the absence of mathematical analysis of the
’integral equations’ involved in these papers”. He pointed to the need for a more rigorous mathematical treatment of TDIEs in order to develop a provably stable scheme. He
reached the conclusion that 1st kind TDIEs and their Galerkin approximations should be
used in favor of 2nd kind TDIEs and collocation. In fact, it was shown in his paper that
the former naturally arose from energy identities with which coercivity measures could be
developed. Such Galerkin schemes, he claimed, had been developed as early as 1989 [26],
though much of the early work by mathematicians on variational formulations for TDIEs
was largely unknown to the engineering community.
While these methods were provably stable in their analytical forms, their discretizations
were still quite difficult to implement exactly. The major difficulty resided in the need to exactly evaluate a four dimensional spatial integral (two source and two testing). The difficulty
arises due to the nature of the spatial integrands. Typically, the temporal basis, being local
and of infinite bandwidth, introduces non-smoothness into the integrands, which degrades
the accuracy of using standard quadrature rules. It becomes necessary, therefore, to evaluate
these integrals in some other way. One possibility is to locate shadow regions and integrate
over each subregion analytically. As mentioned before, these are 4 dimensional integrals, so
this is not practically possible for all 4 dimensions. This raised the following question: if
8

these integrals were evaluated quasi-analytically (3 or less out of the 4 dimensions evaluated
analytically), would the scheme still be stable? This question was examined in [27] for acoustics through a combined collocation/exact evaluation procedure. It was shown numerically
that the scheme produced stable results. This technique was subsequently extended to EM
through the work of [1, 2]. Reference [1] evaluates the source integrals analytically, while the
test integral is computed using quadrature. In [2], the 4D integral is transformed from 2
surface integrals to a line integral and a volume integral. The volume integral is evaluated
analytically and the line integral is evaluated using quadrature. While these quasi-exact
schemes showed robust stability properties, they were still challenging to implement, despite
being simplified through the use of collocation. Furthermore, their extension to higher order
geometries has yet to be accomplished, if it is even possible. The major difficulties lie in (1)
identifying shadow regions on a higher order surface and (2) developing analytical methods
for evaluating these integrals on such a surface. This will be examined in more detail in
Chapter 4.

This motivates the main contribution of this thesis, namely the development of a scheme
which can accurately evaluate the four dimensional spatial integrals involved in Galerkin
discretizations, which can be extended to higher order discretizations. Rather than resorting
analytical integration, we employ a separable expansion in space and time of the retarded
potential Green’s function, which allows all integrals (both testing and source) to be handled
numerically to arbitrary precision. Therefore, shadow regions do not need to be determined
and analytical methods do not need to be derived for higher order surfaces. The main details
of this method are presented in Chapter 4. In Chapter 4, we will discuss the error incurred
in the introduction and truncation of this expansion into the Galerkin framework. We will
9

show that it can be practically implemented, without loss of accuracy, while stabilizing the
solution. We will also give details for its implementation in various discretization frameworks
in Chapters 7 and 8 as well as its use with the Plane Wave Time Domain algorithm in Chapter
6.
This thesis is organized as follows. In Chapter 2 we begin from Maxwell’s equations
and derive the TDIEs we wish to discretize and solve. Chapter 3 summarizes the numerical
discretization of the TD-EFIE, MFIE, and CFIE using the Method of Moments (MoM).
This Chapter will also introduce a higher order space-time Galerkin framework. The major
contribution of this thesis is presented in Chapter 4, in which we outline the separable
expansion and demonstrate its validity, both analytically and numerically. We will also
discuss how singular and near-singular integrals can be accurately evaluated on higher order
tessellations. Chapter 6 will show how acceleration is achieved using the Multilevel Plane
Wave Time Domain Algorithm. In Chapter 5 we will present scattering results on both lower
and higher order surfaces to demonstrate the stability and accuracy of the scheme. Chapter
7 covers the extension of the separable expansion method to penetrable scatterers. Chapter
8 details the implementation of the scheme within the Generalized Method of Moments
discretization framework. Lastly, Chapter 9 will draw conclusions.

10

Chapter 2
Mathematical Preliminaries
In this Chapter, we begin from Maxwell’s Equations and derive the relevant TDIEs we
aim to solve, i.e., the TD-EFIE and TD-MFIE for PEC structures. These will then be
discretized via the Method of Moments. Discussion of the discretization procedure will be
presented in Chapter 3

2.1

Maxwell’s Equations and Wave Equations

We begin from the point form of Maxwell’s Equations in the time domain in free space
in the absence of magnetic sources:

∂
B (r, t)
∂t

(2.1a)

∂
D (r, t) + J (r, t)
∂t

(2.1b)

∇ × E (r, t) = −

∇ × H (r, t) =

∇ · D (r, t) = ρ (r, t)

(2.1c)

∇ · B (r, t) = 0

(2.1d)

where E denotes the electric field intensity, H the magnetic field intensity, B the magnetic
flux density, D the electric flux density, J the current density, and ρ the charge density. In
11

addition to Maxwell’s equations, we have the continuity condition

∇ · J (r, t) = −

∂
ρ (r, t)
∂t

(2.1e)

and the constitutive relations
D (r, t) = 0 E (r, t)

(2.1f)

B (r, t) = µ0 H (r, t)

(2.1g)

where 0 and µ0 denote the permittivity and permeability of free space, respectively. For
perfectly conducting scatterers, the fields satisfy the following boundary conditions on the
surface of the scatterer:
ˆ (r) × E (r, t)|r∈Ω = 0
n

(2.2a)

ˆ (r) × H (r, t)|r∈Ω = J (r, t)
n

(2.2b)

These equations completely characterize the fields for all space and time. To derive the
wave equations for E and H, we take ∇×(2.1a) and ∇×(2.1b), respectively. Beginning with
the wave equation for E,

∇ × (∇ × E (r, t)) = −

∂
∇ × B (r, t)
∂t

(2.3a)

Substituting (2.1b), (2.1f), and (2.1g) into (2.3a) yields

∇ × ∇ × E (r, t) = −µ0

∂
0 ∂t E (r, t) + J (r, t)

∂
∂t

1 ∂2
∂
→ ∇ × ∇×E (r, t) +
E (r, t) = −µ0 J (r, t)
2
2
∂t
c ∂t

12

(2.3b)

{E (r, t) , H (r, t)}

ˆ (r)
n

J (r, t)

{E (r, t) , H (r, t)} = 0
Ω

Figure 2.1: General description of a PEC scattering problem

√
where c = 1/ 0 µ0 is the speed of light in free space. A similar derivation for H results in
1 ∂2
∇ × ∇ × H (r, t) +
H (r, t) = ∇ × J (r, t)
c2 ∂t2

2.2

(2.4)

Mixed Potential Solutions for Integral Equations

Derivations of EM relations are often facilitated through the use of auxiliary functions
called potentials. We begin by examining equation (2.1d), which implies that B can be
written as the curl of a vector function, i.e.

B (r, t) = ∇ × A (r, t)
13

(2.5)

The quantity A is the magnetic vector potential. Substituting (2.5) into (2.1a) leads to
∂
∇ × A (r, t)
∂t
∂
→ ∇× E (r, t) + A (r, t)
∂t

∇ × E (r, t) = −

(2.6a)
=0

which implies the function E (r, t) + ∂ A (r, t) can be written as the gradient of a scalar
∂t
function, i.e.
E (r, t) +

∂
A (r, t) = −∇Φ (r, t)
∂t

(2.6b)

where the quantity Φ is the electric scalar potential. Substituting (2.5) and (2.6b) into (2.1b)
yields

∇ × H (r, t) =
=

1
∇ × ∇ × A (r, t)
µ0
1
∇ (∇ · A (r, t)) − ∇2 A (r, t)
µ0

∂
= 0 E (r, t) + J (r, t)
∂t
∂
∂
=− 0
A (r, t) + ∇Φ (r, t) + J (r, t)
∂t ∂t
1 ∂2
1 ∂
Φ (r, t) + ∇ · A (r, t) − µ0 J (r, t)
→ ∇2 A (r, t) −
A (r, t) = ∇
c2 ∂t2
c2 ∂t

(2.7a)

Likewise, (2.1c) becomes

∇ · E (r, t) = −

∂
∇ · A (r, t) − ∇2 Φ (r, t)
∂t

= ρ/ 0

14

(2.7b)

We introduce the Lorentz gauge condition

1 ∂
Φ (r, t) + ∇ · A (r, t) = 0
c2 ∂t

(2.8)

which reduces (2.7a) and (2.7b) to

∇2 A (r, t) −

1 ∂2
A (r, t) = −µ0 J (r, t)
c2 ∂t2

(2.9a)

1
1 ∂2
Φ (r, t) = − ρ (r, t)
c2 ∂t2
0

(2.9b)

and
∇2 Φ (r, t) −

respectively. The use of this gauge leads to a symmetric system of two inhomogeneous wave
equations. The Green’s function for (2.9b) or for each component of (2.9a) is given by (see
e.g. [28])
δ (t − |r|/c)
4π|r|

g (|r|, t) =

(2.10)

Equivalently, the vector and scalar potentials are given as
∞
A (r, t) =

R3

dr

−∞

dt g R, t − t

∞
Φ (r, t) = −

R3

dr

−∞

t
dt g R, t − t

−∞

J r ,t

dtˆ∇ · J r , tˆ

where R = |r − r | and ∇ means the vector derivative is taken with respect to r .

15

(2.11a)

(2.11b)

2.3

Electric and Magnetic Field Integral Equations

In this Section, we give a brief derivation of the TD-EFIE and TD-MFIE for PEC surfaces.
For a more general derivation we refer the reader to [28, 29]. As in [29], we begin by deriving
the frequency domain equations and applying an inverse Fourier transform to the results.
Consider a region V residing between two boundaries, S and S1 , and let Σ = S + S1 . The
vector Green’s theorem states that, for two vector functions with continuous derivatives up
to 2nd order, for r ∈ Σ

V

(P × ∇ × Q − Q × ∇ × P) · ds

(Q · ∇ × ∇ × P − P · ∇ × ∇ × Q) dv =
V

(2.12a)

Σ

where Q is chosen to be
e−jk|r−r |
Q=a
ˆψ˜ r, r , ω = a
ˆ
|r − r |

(2.12b)

√
where a
ˆ is an arbitrary vector and k = ω µ0 0 . The function P is chosen to be E or H for
the EFIE or MFIE, respectively, where
∞
E (r, ω) = F {E (r, t)} =

−∞

dtE (r, t) e−iωt
(2.13a)

= Einc (r, ω) + Es (r, ω)

H (r, ω) = F {H (r, t)}
(2.13b)
= Hinc (r, ω) + Hs (r, ω)
where {Einc , Hinc } and {Es , Hs } are the incident and scattered fields, respectively. For
16

ˆ (r)
n

S
V
ˆ (r)
n

S1

Figure 2.2: General domain

PEC structures (2.12a) reduces to

ˆ (r) ׈
ˆ (r) × Es (r, ω)
n
n (r) × Einc (r, ω) = −ˆ
n (r) × n
=

S

1
4πjω 0

(2.14a)

ω 2 µ0 0 J r , ω ψ˜ r, r , ω − ∇ · J r , ω ∇ ψ˜ r, r , ω

J (r, ω) = 2ˆ
n (r) × Hinc (r, ω) +

· ds

1
ˆ (r) × − J r , ω × ∇ ψ˜ r, r , ω · ds
n
2π
S

(2.14b)

where
J (r, ω) = F {J (r, t)}
which are the frequency domain EFIE and MFIE, respectively.

(2.14c)
In order to derive the

TD-EFIE and MFIE, we take the inverse Fourier transform of (2.14a) and (2.14b), using the
17

E inc (r, t) , Hinc (r, t)

ˆ (r)
n
{E s (r, t) , Hs (r, t)}

J (r, t)

{E (r, t) , H (r, t)} = 0
S

Figure 2.3: Scattering from PEC object

following definitions and properties
V (r, t) = F −1 {V (r, ω)}
∞
1
=
dωV (r, ω) ejωt
2π −∞

(2.15a)

F −1 {V1 (r, ω) V2 (r, ω)} = v1 (r, t) ∗ v2 (r, t)
(2.15b)

∞
=

−∞

dt v1 r, t

v2 r, t − t

∗
F −1 {V1 (r, ω) × V2 (r, ω)} = V 1 (r, t) × V 2 (r, t)
(2.15c)

∞
=

−∞

dt V 1 r, t

18

× V 2 r, t − t

All integrals are assumed to be interchangeable meaning the integrands are square integrable.
Inverse Fourier transforming (2.14a) and (2.14b) results in

ˆ (r) × n
ˆ (r) × E inc (r, t) = F −1 n
ˆ (r) × n
ˆ (r) × Einc (r, ω)
n
=−

1
dr
4π S
+

µ0 F −1 jωJ r , ω ψ˜ r, r , ω
1
0

=−

dr
S

µ0

F −1

1
∇ · J r , ω ∇ ψ˜ r, r , ω
jω

t
∂
1
dt ∇ · J r , t
J r , t ∗ g (R, t) + ∇ g (R, t) ∗
∂t
−∞
0

.
= L˜e {J (r, t)}

(2.16a)

Likewise, the TD-MFIE can be shown to be
J (r, t)
ˆ (r) × Hinc (r, t) =
ˆ (r) × Hs (r, t)
n
−n
2
∗
J (r, t)
ˆ (r) × − dr J r , t × ∇g (R, t)
=
−n
2
S

(2.16b)

.
= Lh {J (r, t)}
(2.16a). We note that, for open structures, only (2.16a) is valid. The integral equation
(2.16b) is of the 2nd kind and typically more stable than (2.16a).

As discussed in the introduction, the discrete solution of (2.16a) and (2.16b) can be
corrupted by resonant modes of the structure. In order to alleviate this problem a linear
19

ˆ ×(2.16a) and (2.16b) can be used [15]
combination of n

ˆ (r) × n
ˆ (r) × E inc (r, t) + (1 − αc )ˆ
αc /η0 n
n (r) × Hinc (r, t)
= αc /η0 Le {J (r, t)} + (1 − αc )Lh {J (r, t)}

(2.16c)

.
= Lc {J (r, t)}
.
ˆ × L˜e {·}, and η0 =
where αc ∈ [0, 1], Le {·} = n

µ0 / 0 is the intrinsic impedance of free

ˆ ×(2.16a) rather
space. Equation (2.16c) is the TD-CFIE. The justification for discretizing n
than discretizing (2.16a) is given in the next Chapter.

20

Chapter 3
Method of Moments Solution of TDIEs
The procedure of discretizing (2.16a), (2.16b), or (2.16c) is equivalent to reducing one
of these from an integral equation to a matrix equation. One such way of achieving this
reduction is via the method of moments. The method of moments is an error minimization
method based on projections of the operators and excitations onto testing function spaces. In
the time domain, these projections are performed in both space and time. In this chapter we
outline two procedures for setting up a symmetric MOT system. The first uses a collocation
testing procedure in time while the second uses Galerkin testing in the temporal dimension.
The former is attractive for its simplicity, but achieving higher order accuracy in time is
challenging if stability is to be maintained. The latter is more challenging and expensive
but, as will be shown numerically in chapter 5, can achieve higher order accuracy and late
time stability.

3.1

Collocation in time

We begin with equations (2.16a), (2.16b), and (2.16c) and assume the scatterer boundary,
ˆ composed of Ne non-overlapping facets, i.e.,
S, has been approximated by a surface, S,
Ne
Sˆ =

Tn
n=1

Ne
Tn = ∅
n=1
21

(3.1)

In this work, the facets, Tn , are assumed to be either flat triangles or curvilinear triangles
(see e.g. [30]). The surface current is approximated as an expansion of Ns spatial and Nt
temporal basis functions.

ˆ (r, t)
J (r, t) ≈ J
Nt

Ns
=

Sn (r)
n=1

j
In T (t − j∆t )

(3.2)

j=1

Where the functions T (t) are typically shifted Lagrange polynomials (see e.g. [1]), ∆t denotes
j
the time step, and In are the unknown coefficients. The functions Sn (r) are RWG functions
ˆ Therefore, Ns is determined by the number of
[31] associated with each interior edge on S.
ˆ
interior edges on S.
ˆ into
The next step in the discretization is to substitute the approximate current, J,
ˆ ×(2.16a), (2.16b), or (2.16c) for J . Focusing on (2.16a), we have
n

ˆ (r) ׈
ˆ (r) × n
ˆ (r) ×
n
n (r) × E inc (r, t) ≈ n

Sˆ

dr

t
1
∂ˆ
ˆ r ,t
dt ∇ · J
J r , t ∗ g (R, t) + ∇ g (R, t) ∗
∂t
−∞
0

Ns Nt
Sn r
j
ˆ (r) ×
=ˆ
n (r) × n
dr µ0
In
T (t − R/c − j∆t )
R
S
n
n=1 j=1

t
1
+ ∇ · Sn r ∇ g (R, t) ∗
dt T t − j∆t 
−∞
0
µ0

(3.3)

where Sn is the spatial support of Sn (r). The prime on T (g(t)) denotes the derivative
with respect to the argument. Next, the system is tested by taking a space-time inner
product with spatial and temporal testing functions (see (3.4b)). Typically, Galerkin testing
22

is used in space, i.e., the spatial testing functions are chosen to be the same as those used to
approximate the current. In time, point testing is used with delta functions. As discussed in
the introduction, Galerkin testing is the preferred method in developing a stable framework.
While the use of the basis and testing functions presented here seems to go against this, it
can be easily shown that a delta testing/linear Lagrange basis function pair is equivalent (in
the impedance matrix elements) to Galerkin testing with rectangular pulse basis functions.
Applying this procedure to (3.3) results in

ˆ (r, t)
δ (t − i∆t ) Sm (r) , E inc (r, t) ≈ δ (t − i∆t ) Sm (r) , Le J

1 
j
=
In
drSm (r) ·
dr µ0 Sn r T ((i − j)∆t − R/c)
4π  Sm
S
n
n,j
1 ∇ · Sn r
j
+
dr∇ · Sm (r)
In
dr
R
Sm
Sn
0
n,j

(i−j)∆t −R/c
−∞

(3.4a)



dt T (t )


where the inner product is defined as
∞
f (t)A (r) , B (r, t) =

−∞

dtf (t)

R3

drA (r) · B (r, t)

(3.4b)

The same testing procedure is applied to (2.16b) with the approximated current to yield

ˆ (r, t)
ˆ (r) × Hinc (r, t) ≈ δ (t − i∆t ) Sm (r) , Lh J
δ (t − i∆t ) Sm (r) , n
1
1
j
ˆ (r) ×
In Sn (r) T ((i − j) ∆t ) −
− drSm (r) · n
= − drSm (r) ·
2 Sm
4π Sm
n,j
j
In − dr
Sn
n,j

T ((i − j)∆t − R/c) T ((i − j)∆t − R/c)
+
cR
R2

23

Sn r

ˆ
×R

(3.4c)

where
ˆ = r−r
R
|r − r |

(3.4d)

This testing procedure is performed for m = 1, 2, ..., Ns to yield the system
i−1
Z 0I i = V i −

Z i−j I j

(3.5a)

j=1
at the ith time step, where

Z i−j

mn

= δ (t − (i − j) ∆t ) Sm (r) , Lc {Sn (r) T (t)}
j j
j T
I j = I1 I2 ... IN
s

ˆ (r) × Hinc (r, t)
V i m = δ (t − i∆t ) Sm (r) , αc /η0 E inc (r, t) + ((1 − αc ) n

(3.5b)

(3.5c)
(3.5d)

We have assumed, without loss of generality, that (2.16c) is being discretized. The current
coefficients I i are obtained by solving (3.5a) and substituted into the right hand side of
(3.5a) in subsequent steps. This procedure is known as marching on in time. It is clear from
j
(3.5a) that any errors in In are recursively fed into the system and can lead to exponential
growth terms. The next chapter reviews a recent method for stabilizing (3.5a) and outlines
an alternative, which is the centerpiece of this thesis.

3.2

Higher order Space-time Galerkin Discretization

As discussed in the introduction, the stability properties of Galerkin methods, being
derived from energy identities, are much better understood than those of methods based
24

on collocation. The stability of collocation based methods can be analyzed using Fourier
analysis [32, 33], but this must be done for specific problems and, in general, stability of such
techniques can not be guaranteed. The previous section outlined a discretization framework
based on collocation with a linear temporal basis. This specific setup can be shown to be
equivalent to Galerkin testing with 0th order temporal basis functions and, therefore, enjoys
the stability properties of Galerkin discretizations. Extension to higher order temporal basis
functions, however, is more involved.

3.2.1

Higher order temporal basis

To motivate the difficulties alluded to above, consider a pth order temporal basis function
T (t). Let supp(T (t)) = [a, b], where b > a. Now consider the inner product,

T (t − i∆t ), T (t − j∆t ) = δ(t), T (t) ∗ T (−t − (j − i)∆t ) .

(3.6)

where i and j are integers (uniform temporal discretization). It can be shown that the shifted
autocorrelation function,

T (t)∗ T (−t − (j − i)∆t )|t=0 = 0

(3.7)

for j − i ∈ (a − b, b − a) .

Causality of the marching system requires that (3.6) is zero for j > i. This is true if and
only if the only integer j − i ∈ (a − b, b − a) = 0, which implies b − a < 1.
This result means that the temporal interpolants most commonly used in MOT schemes,
i.e., shifted Lagrange interpolants (see e.g. [1]) and splines [10]) will yield non-causal MOT
25

systems if used in a space-time Galerkin framework, because they span multiple time steps.
One possible route to overcoming this barrier was published in [34]. In that work, the
authors employed a Petrov-Galerkin scheme by expanding the temporal basis using multiple
functions defined on the same interval of support. These temporal basis functions effectively
had support over two time steps. This method also used multiple temporal testing functions
within the same interval of support, but with support restricted to 1 time step, which did
not produce a discretely non-causal system. The testing functions were designed to be 1
degree less smooth than the basis functions. As this discretization was applied to the time
differentiated TD-EFIE, these were the proper spaces, as prescribed by Terrasse in [35].
In this thesis, we are interested in discretizing the undifferentiated TD-EFIE. Therefore,
a symmetric Galerkin system is more appropriate. To retain discrete causality in such a
system, it is necessary to restrict the basis function support to 1 time step. However, we
will follow the procedure of [34] by using multiple functions over the same temporal support.
In chapter 5 we will show one result in which the separable expansion is used to generate a
stable result for the temporal discretization scheme of [34].

With the restriction of temporal support in mind, the surface current can be approximated by these new basis functions as
Nt

Ns
ˆ (r, t) =
J

Sn (r)
n=1

where

p

j,l
In T l (t − j∆t )

(3.8a)

j=1 l=0





 l (t)
l
T (t) =



0
26

t ∈ (−∆t , 0]
otherwise

(3.8b)

Here p denotes the order of the temporal basis. There is a great amount of flexibility in
choosing the functions l (t). The crucial property, as covered previously, is their support
over one time interval. For use in the TD-EFIE, it is essential that the correlation functions,

T kl (t) = T k (t) ∗ T l (−t)

(3.8c)

support at least one derivative. In this thesis, the results presented will use, for l (t), either
Legendre polynomials or Lagrange polynomials given by
p
l (t) =
i=0
i=l

t − ti
tl − ti

(3.8d)

where tl = (l/p − 1)∆t , l = 0, 1, ..., p.

Direct substitution of (3.8a) into (3.5a) would result in an Ns Nt × (p + 1)Ns Nt system.
Therefore, we need to increase the number of equations by a factor of p + 1 as well. To
create a square system we test each equation using Sm (r) T k (t − i∆t ) for k = 0, 1, ...p,
m = 1, ..., Ns , and i = 1, ..., Nt . The modified matrix elements are given by










Z i−j = 









11
Z i−j
21
Z i−j

12
Z i−j
12
Z i−j

.

...
...

1Ns
Z i−j
1Ns
Z i−j
.

.

.

.

.

Ns 1
Z i−j

Ns Ns
... Z i−j
27





















(3.9a)

V im

k

j,0 j,1
j,p j,0
j,p T
I j = I1 I1 ... I1 I2 ... IN
s

(3.9b)

T
V i = V i1 V i2 ... V iNs

(3.9c)

ˆ (r) × Hinc (r, t)
= T k (t − i∆t ) Sm (r) , αc /η0 E inc (r, t) + ((1 − αc ) n

(3.9d)

mn
where k, l = 0, 1, ..., p. Each sub-block, Z i−j , is of dimension (p + 1) × (p + 1), with entries
mn
Z i−j

3.2.2

kl

= T k (t − (i − j) ∆t ) Sm (r) , Lc Sn (r) T l (t)

(3.9e)

Interpolation properties of higher order basis

In this subsection we quantify the interpolation accuracy of the higher order basis. To
this end, we interpolate a modulated Gaussian pulse
2
2
g(t) = cos(2πf0 t)e−(t−tp ) /2σ

(3.10)

where, given some modulating frequency f0 and maximum frequency fmax , the values σ
and tp are defined as σ = 3/(2πB) and tp = 6σ, where B = fmax − f0 . This pulse is shifted
and approximated using the higher order basis as
Nt

p
g j,l T l (t − j∆t )

g(t) ≈ g˜(t) =

(3.11)

j=1 l=0
The time step is defined as ∆t = 1/(2ksamp fmax ), where ksamp > 0 is the oversampling
factor. As the expansion is used in a space-time Galerkin scheme, we obtain the coefficients
28

g j,l by solving, at each time step i = 1, 2, ..., Nt , the system

bi = Axi

(3.12a)

bi k = T k (t − i∆t ) , g(t) ,
A

kl

= T k (t) , T l (t) ,

k = 0, 1, ..., p

k, l = 0, 1, ..., p

xi = g i,0 , g i,1 , ..., g i,p
.
f1 (t), f2 (t) =

T

(3.12b)
(3.12c)
(3.12d)

∞
−∞

dtf1 (t)f2 (t)

(3.12e)

Given a random shift, ∆ ∈ (0, ∆t ), the error between the approximate and exact solution is
computed as
Error =

g approximate − g analytic

(3.13a)

g analytic
g approximate = [˜
g (∆t − ∆) g˜(2∆t − ∆) ...˜
g (Nt ∆t − ∆)]

(3.13b)

g analytic = [g(∆t − ∆) g(2∆t − ∆) ...g(Nt ∆t − ∆)]

(3.13c)

where the 2 -norm is used. The result for 1st, 2nd, and 3rd order Lagrange polynomials is
shown in figure 3.1 with a shift of ∆/∆t = 0.906.

As expected by theory, the error scales

as O(p + 1). This experiment will be repeated in the following chapters with this result used
as a reference.

29

Figure 3.1: Interpolation errors for higher order basis

30

Chapter 4
Separable Expansion Method
As discussed in Chapter 1, accurate computation of matrix elements is vital to the stability of MOT solutions to TDIEs. In this Chapter, we discuss why this is a challenge, detail
some existing remedies and present the remedy that forms the crux of this thesis. We will
review quasi-exact integration schemes, which have been used to great effect in stabilizing
TDIE solutions [1, 2]. As will be shown, these methods are restricted to flat tessellations due
to their reliance on analytic integration over shadow regions. This will motivate the main
focus of this Chapter and central contribution of this thesis, i.e., a separable expansion in
space and time of the retarded potential Green’s function. The key feature of this method
is that it enables accurate evaluation of retarded potential integrals, as do quasi-exact integration schemes, but the procedure here is entirely numerical, i.e., none of the integrals are
calculated analytically. This opens avenues for the method to be used on higher order surface
descriptions, which has not been accomplished and is perhaps impossible using quasi-exact
integration schemes. Validation of the stability properties of this method will be made via
numerical experiment in Chapter 5.

4.1

Accurate evaluation of retarded potential integrals

As discussed in Chapter 1, inaccuracies in the computation of matrix elements is a primary cause of instability in MOT solutions to TDIEs. To elucidate the challenges faced in
31

computing these matrix elements, consider the 4-dimensional integral

ij
Imn =

dr∇ · Sm (r)

dr ∇ · Sn r

Tn
ψ ((i − j)∆t − R/c)
dr
dr
=C
R
Tn
Tm
Tm

ψ ((i − j)∆t − R/c)
R

(4.1)

t dt T (t ) and T ⊂ supp {S (r)} and T ⊂ supp {S (r)} are assumed
where ψ(t) = −∞
m
m
n
n
to be triangles. A constant C has been removed from the integrand as we assume Sn (r)
to be an RWG basis function, whose divergence is constant on Tn . Such integrals appear
in the scalar potential due to one spatio-temporal basis function with Galerkin testing applied in space and point testing applied in time. For clarity we restrict our discussion to
the scalar potential term. The procedure for applying these methods to the vector potential is straightforward. Note that RWG functions have support over two triangles, but we
restrict our discussion to one test and one source triangle for simplicity. Typically, in MOT
discretizations, T (t) and

dtT (t) are only smooth to 0th (continuous) and 1st order, respec-

tively, and ∂t T (t) is discontinuous. This is shown for the oft-used shifted 1st order Lagrange
polynomials (hat or chapeau functions) in figure 4.1.

Therefore, integral (4.1) cannot be

evaluate accurately using numerical quadrature. In the next Section we review a class of
methods for accurately evaluation (4.1)

4.2

Quasi-Exact Integration Schemes

Although the exact evaluation of the four-dimensional integral (4.1) is impractical, schemes
have been developed [1, 2, 36] to evaluate the integral to high precision by approximating a
subset of the four dimensions and evaluating the remaining dimensions exactly on flat trian32

(a) : Undifferentiated

(b) : Integral

(c) : First derivative
Figure 4.1: 1st order shifted Lagrange polynomial: (a) undifferentiated, (b) integral, (c) 1st
derivative

33

gles with RWG functions for the spatial basis. These are the quasi-exact integration methods.
In [1, 36], analytical expressions are given for evaluating the 2-dimensional source integral
exactly. In [1], the testing integral is approximated using numerical quadrature to yield a set
of observations points. Then, for a given test point r, the procedure for accurately evaluating
the source integral is as follows. The first step is to determine the subregions of Tn over
which T (k∆t − R/c) is smooth. Let supp {T (t)} = [−1, P ∆t ], where P is an integer. The
domains over which T (k∆t − R/c) is smooth are given by
R
k∆t − ∈ (l, l + 1) ∆t
c

(4.2)

or R ∈ c∆t (k − l − 1, k − l)

l = −1, 0, ..., P − 1

which describes the region between two spheres of radii c∆t (k − l − 1) and c∆t (k − l) centered at r. Denoting this region as Dk,l , we must find

Cn,k,l = Dk,l

Tn

(4.3)

Cn,k,l is the region on Tn over which T (k∆t − R/c) is smooth. Determining Cn,k,l is
tantamount to finding the arcs of intersection between spheres and source elements, Tn
(see figure 4.2).

Analytical expressions exist for the integrals over these subregions. The

procedure for finding the subregions, as described in [36], is to determine all the intersection
points between these spheres and the edges of Tn . Next, these points need to be paired such
that each pair corresponds to the end points of an arc of intersection (boundary between
smooth subregions) on Tn . This is a cumbersome task as suggested by figure 4.3. There
are multitudinous ways in which a sphere can intersect a triangle and there can be up to 6

34

Figure 4.2: Regions of smoothness on source element (graphic from [1])

points of intersection between the sphere and triangle edges.

A second way of computing (4.1) is presented in [2], in which only one dimension of (4.1)
ij
is approximated. The integral Imn is rewritten as

ij
Imn = C
=C

dr
Tm

dr
Tn

ψ ((i − j)∆t − R/c)
R

drΠTm (r)
dr ΠTn r
R3
ψ ((i − j)∆t − |r|/c)
=C
dr
Ω(r)
|r|
R3
R3

(4.4)

ψ ((i − j)∆t − R/c)
R

The integrals have been transformed from surface to volume integrals by introducing the
35

Figure 4.3: Various intersections between spheres and source triangles (graphic from [1])

indicator functions ΠTm (r), defined as

ΠTm (r) =





1 r ∈ Tm



0 otherwise
36

(4.5)

Ω(r) is termed the correlation function and is given by

Ω(r) =

R3

dr ΠTm r

ΠTn r − r

(4.6)

Next, the last line of (4.4) is decomposed into a line and surface integral. The line integral
is discretized and evaluated numerically. At each quadrature point on the line, the resulting
surface integral is taken over a slice of the correlation function Ω(r). The support of a
slice of Ω(r), being the correlation between two line segments, is composed of polygonal
blocks. The form of Ω(r) on these surfaces is a summation of multinomials. The form of
the integrand and the domain of integration lend themselves well to analytical evaluation.
The line integral itself, however, is not smooth and it is therefore necessary to find critical
points on the line, between which the integrand is smooth. The resulting line integrals are
evaluated numerically over each smooth sub-domain.

4.3

Integration on Curvilinear Elements

One restriction is common to both quasi-exact integration methods [1] and [2], i.e., they
apply only to flat domains. For [1], extension to curved elements leads to two complications.
The first complication is in determining the subregions Cn,k,l . The intersection between a
sphere and a curved triangle is no longer a simple circle, and the boundaries between subdomains are not simply the arcs of a circle. Finding the intersection points on the edge of the
curved triangle is made much more difficult, if not impossible. The second complication arises
in deriving analytical expressions for the integrals over these sub-domains. For flat triangles,
the integrals transform nicely into a polar coordinate system. For curved triangles, even if
37

the boundaries of each sub-domain could be found, the integrals will not lend themselves to
such a convenient transformation in general. The complications to [2] are similar and arise
from the form of Ω(r) for indicator functions on curved elements.
This motivates the main contribution of this thesis. To accurately evaluate (4.1) on
curvilinear elements it would appear necessary to move away from schemes that rely on
analytical integration. In what follows, we develop a method that can compute source
integrals to arbitrary precision via standard quadrature rules. We will show that the spatial
integrands become entirely smooth. Therefore, there are no shadow regions to be determined
and no need to evaluate the integrals analytically. Instead, standard quadrature rules for
triangles can be utilized.

4.4

Separable Expansion

In this Section, we present a scheme which can accurately evaluate retarded potential integrals without resorting to analytical integration. This will enable it to be trivially extended
to higher order geometric discretizations. We begin by considering the two-dimensional
source integral
δ (t − R/c) ∗ T (t)Sn r
ψ (r, t) =

dr

R

Tn

(4.7)

This type of integral appears in the undifferentiated vector potential due to one spatiotemporal basis function. Next, we approximate this integral as

ψ (r, t) = δ(t − ζ/c) ∗

dr
Tn

∞ a P (k (R − ζ)/c + k ) P (k t + k ) ∗ T (t)S r
n
2 l 1
2
l=0 l l 1
R
(4.8)
38

M
≈ δ(t − ζ/c) ∗

al Pl (k1 t + k2 ) ∗ T (t)
l=0

dr

Pl (k1 (R − ζ)/c + k2 ) Sn r
R

Tn

where al = k1 (2l + 1) /2, k1 α + k2 = −1, k1 β + k2 = 1, β > α ≥ 0, and Pl (k1 x + k2 )
are Legendre polynomials orthogonal over x ∈ [α, β]. ζ represents the largest multiple of
c∆t between r and Tn , as shown in figure 4.4. This shift is used to minimize the spatial
support of the Legendre polynomials. It is clear from (4.8) that, given proper choice of α
and β, the spatial integrand is smooth over the entire domain (see figure 4.4). This allows

β
α
ζ

r

Tn
R−ζ

r
c∆t

Figure 4.4: Spatial support of expansion (4.8)

for the accurate numerical evaluation of (4.8) without having to determine shadow regions
as in [1, 2]. In this method Tn must lie completely outside the sphere centered at r of radius
ζ and must lie completely inside the sphere centered at r of radius ζ + (β − α)c. In other
words, determining the support of the Legendre polynomials in (4.8) is equivalent to finding
a spherical shell, denoted as B[cα,cβ] (R − ζ), that completely encloses Tn , where




1 |r| ∈ [cα, cβ]

B[cα,cβ] (|r|) =



0 otherwise
39

(4.9)

This is a much simpler task than the determination of shadow regions. Finding the minimum
and maximum distances on a flat triangle from r is a simple exercise. These same distances
(the minimum and maximum distance to the flat triangle formed by the curved triangle’s
vertices) can be used as an initial estimate when curvilinear elements are used. The inner
and outer radii of the spherical shell can be decreased and increased based on the curvature
of the element. The simplicity of this procedure enables the algorithm to be extended to
higher order tessellations whereas [1] and [2] are restricted to flat tessellations.

4.5

Truncation Error

While the integrand in (4.8) is clearly smooth over Tn , it remains to be seen whether the
truncation in (4.8) is valid and its upper limit practical. In this Section we look at the error
incurred through truncation of (4.8) to determine its practicality. To derive error bounds, we
use Fourier analysis. Such analysis is facilitated through the use of indicator functions and
pulse functions, which allow us to express the integrand in (4.8) in terms of entire domain
functions. To begin we define




1 t ∈ supp {T (t)}

ΠT (t) =



0 otherwise

(4.10a)

which is an indicator function for T (t), similar to ΠTn (r) from Section 4.2. Additionally,
we define





1 t ∈ [α, β]

H[α,β] (t) =



0 otherwise
40

(4.10b)

which is a pulse function defined over the regions of orthogonality of Pl (k1 t + k2 ), similar to
the spherical shell of (4.9), which is displayed in figure 4.5. In what follows, we will assume

ζ + cβ
r
ζ + cα

Figure 4.5: Domain of the spatial integrals

ζ = 0. Using these definitions, (4.8) can be written as a space-time convolution, which will
transform to multiplications in the space-time frequency domain. Rewriting (4.8), we have

∞

δ t−t − R
c
T (t )
4πR

dt ΠT t
dr ΠTn r Sn r
−∞
R3
M
Pν (k1 R/c + k2 )
≈
al
dr ΠTn r Sn r B[cα,cβ] (R)
4πR
R3
l=0
∞
dt ΠT t T (t )H[α,β] t − t Pl k1 t − t + k2
−∞

ψ (r, t) =

41

(4.11)

In order to determine the truncation error, we wish to evaluate
∞

.
M (r, t) =

R3

dr ΠTn r

Sn r

−∞

dt ΠT t

δ t−t − R
c
T (t )
4πR

M
−

al
l=0

dr ΠTn r
R3

Sn r

B[cα,cβ] (R)

Pν (k1 R/c + k2 )
4πR

(4.12)

∞
−∞

dt ΠT t

T (t )H[α,β] t − t

Pl k1 t − t

+ k2

This error will be evaluated in the frequency domain using the following definitions:

V˜ = Ft {V (t)} =

W = Fr {W(r)} =

∞

dtV (t)ejωt

(4.13a)

drW(r)ejλ·r

(4.13b)

−∞

R3

Fr,t {X (r, t)} = Fr {Ft {X (r, t)}}

(4.13c)

We assume ΠTn (r) Sn (r) ΠT (t) T (t) is effectively band-limited in spatial and temporal
frequency to |λmax | and ωmax , respectively. In practice, the upper limit M required for
convergence will be determined by the smaller of these bandwidths. In a method of moments
discretization, |λmax | will depend on the size of a given element and the order of the spatial
basis function. Analogously, ωmax will depend on the time step size and order of the temporal basis function. The goal is to show a bounded error within {|λ| < |λmax |, ω < ωmax },

42

where the error for a given truncation limit is


˜M (λ, ω) = 
Fr,t

δ (t − |r|/c)
4π|r|

−






|r|
Pl k1 c + k2

M

Fr,t B[cα,cβ] (r) H[α,β] (t)
al



l=0

4π|r|





Pl (k1 t + k2 ) 





(4.14)

A straightforward derivation shows that

Fr,t






M
B[cα,cβ] (r) H[α,β] (t)
al



l=0

|r|
Pl k1 c + k2
Pl (k1 t + k2 )

4π|r|

ωk
−j 2 M
c
k1
=
e
(2l + 1)jl
2π|λ|k1
l=0

ω
k1









jl

(4.15)

|λ|c
k1

where jl (z) is a spherical Bessel function, which satisfies the monotonicity relationship (see
[37] 10.14.5)

|jl (z)| ≤

π
2z

z
˜l2 − z 2

˜l +

˜l

˜2 2
e l −z

0<z ≤l+

1
2

(4.16)

where ˜l = l + 12 . Substituting (4.15) into (4.14) and using the bound (4.5) yields

1

k1

1
(2l + 1)e
˜M (λ, ω) ≤

3/2
1/2
4|λ|
ω
l=M +1 
∞

43

2
k1 ˜l −(|λ|c)2 + 1
k1

2
k1 ˜l −ω 2


˜l


|λ|c


2

2
˜
k1 l − (|λ|c)

(4.17)

˜l
2
1
 e ω|λ|c 
˜M (λ, ω) <
(2l + 1) 
2
3/2
1/2
4|λ|
ω
˜
4
k
l
l=M +1
1

˜l
∞
1
 e2 ω|λ|c 
(2l + 1) 
<
2
4|λ|3/2 ω 1/2 l=M +1
˜
4 k1 M

(4.18)






ω
k1 ˜l +

k1 ˜l

2

− ω2

˜l 
 
 
 

k1 ˜l +

which can be relaxed to the simplified form



∞

˜ = M + 3/2. This is a standard arithmetico-geometric series, which converges to
where M

1
4|λ|3/2 ω 1/2

2
e2 ω|λ|c
˜ 2 − e2 ω|λ|c
4k12 M

˜
M
e2 ω|λ|c
˜2
4k12 M

˜ +3
2M

˜2
4k12 M
˜ +1
− 2M
e2 ω|λ|c
(4.19)

˜
Setting zt = eω/ 2k1 M

˜M (λ, ω) <

˜
and zs = e|λ|c/ 2k1 M

1
(zt zs )M
3/2
1/2
4|λ|
ω

and using (4.19) it can be shown that

7/2 (2M + 6)
zt zs
− (2M + 4)
1 − zt zs
zt zs

(4.20)

for (zt zs ) < 1 this series is rapidly convergent in M .

Next, we numerically test this bound. For simplicity, we only focus on the temporal
portion of the expansion. In this test, the expansion

fˆ(ω) = Ft


M


l=0



al Pν (k1 R/c + k2 ) Pl (k1 t + k2 )

44

(4.21)

is computed for a given value of R and compared against

f˜(ω) = Ft {δ (t − R/c)}

(4.22)

= e−jkR

over a sampling of Nf frequencies in the interval [fm , fM ], fM > fm . The error
fˆ(ω) − f˜(ω)
Error =

=

f˜(ω)
∞
l=M +1 al Pν (k1 R/c + k2 ) Ft Pl (k1 t + k2 )
f˜(ω)

.
f (ω) =

(4.23)

Nf
|f (ωi )|2
i=1

(4.24)

is computed for various values of M with the width of the support of H set to 3∆t . The
result is shown in figure 4.6, while figure 4.7 compares the analytical and approximated
signals convolved with a band-limited function in the frequency domain for various M .

As

is evident from figure 4.6, the expansion (4.8) quickly converges. Figure 4.6 suggests that
the expansion (4.8) is practical.

A remark should also be made regarding the choice of M in practice. Given α (typically
set to 0) and β for a given triangle, Tn , the value of M should be determined by (4.20)
to achieve a given bound. Empirically, we have found that assuming a direct dependence
M = Cβ is sufficient to ensure convergence, with C typically set to either 2 or 3 for a
typical simulation setup. Here, “typical” means edge lengths on the order of λmin /10 and
∆t = 1/ (20fmax ), where fmax = f0 + B is the frequency at which the incident field power
45

Figure 4.6: Convergence in M

is 160 dB below its peak, λmin = c/fmax , and f0 is the center frequency and B is the
effective bandwidth of the incident field described by (5.11).

4.6

Incorporation into Method of Moments Scheme

In this Section we examine (3.5a) modified using the expansion in (4.8). We begin with
the right hand side of (3.4a) for one spatio temporal basis function. Substituting in expansion
46

Figure 4.7: Spectra of approximating signals

(4.8), this becomes
.
δ (t − i∆t ) Sm (r) , LM
e {Sn (r) T (t − j∆t )} =
M

m,n ν
m,n ν
aν αν T˜i,j
+ φν Tˆi,j

(4.25a)

ν=0
where
m,n
= µ0
αν

Sm

drSm (r) ·

Sn r
dr
Sn

Pν (k1 R/c + k2 )
4πR

∇ · Sn r Pν (k1 R/c + k2 )
1
m,n
φν
=
dr∇ · Sm (r)
dr
4πR
Sn
0 Sm
ν = P (k t + k ) ∗ T (t − j∆ )
T˜i,j
ν 1
t t=i∆
2
t
47

(4.25b)

(4.25c)
(4.25d)

ν = P (k t + k ) ∗
Tˆi,j
ν 1
2

t
−∞

dτ T (τ − j∆t )

(4.25e)
t=i∆t

Similarly, the right hand side of the TD-MFIE (3.4c) for one spatio-temporal basis function
becomes
.
δ (t − i∆t ) Sm (r) , LM
h {Sn (r) T (t − j∆t )} =
(4.26a)

M
1
m,n ν
drSm (r) · Sn (r) T ((i − j) ∆t ) −
aν κν Ti,j
2 Sm
ν=0
where
1
m,n
κν
=
drSm (r) ·
4π Sm
1 Pν (k1 R/c + k2 ) k1
− Pl (k1 (R − ζ)/c + k2 ) Sn r
n
ˆ × − dr
R
R
c
Sn
ν = P (k t + k ) ∗ T (t − j∆ )|
Ti,j
ν 1
t t=i∆t
2

(4.26b)
ˆ
×R

(4.26c)

As the modifications from using (4.8) are restricted to the Green’s function, the definitions
of V i and I i are identical to those given in Chapter 3. The spatial integrals (4.25b), (4.25c),
and (4.26b) can be evaluated using standard quadrature rules. The upper limit M is used
to determine the integration rule need for a given expansion. The temporal convolutions
(4.25d), (4.25e), and (4.26c) can be evaluated analytically or via one-dimensional integration
rules and are precomputed for l = 0, 1, ..., M .

4.6.1

Efficient evaluation of scalar potential

A remark is in order for the temporal integration of the basis function in (4.25e). A na¨ıve
discretization of (3.5a) leads to a scaling of O(Nt2 Ns2 ). Instead we make use of a recursive
48

procedure which we outline here. We begin by noting that

ν = P (k t + k ) ∗
Tˆi,j
ν 1
2
i∆t

=

−∞

t
−∞

dτ T (τ − j∆t )
t=i∆t

(4.27a)

t−t
dtPl k1 t + k2

−∞

dτ T (τ − j∆t )

it can be shown that



∞ dτ T (τ )

 −∞
t ≥ γ∆t
t−t
dτ T (τ − j∆t ) =

−∞


 t−t dτ T (τ − j∆t ) t < γ∆t
−∞

(4.27b)

where γ = β − α + i + 1. Therefore, when j ≥ γ we have

ν =
Tˆi,j

∞

∞
−∞

dtT (t)

−∞

dtPν (k1 t + k2 )

(4.27c)

We introduce the charge term
∞
C j = C j−1 + I j−1

−∞

dtT (t), C 1 = 0

(4.28a)

∞ dtP (k t + k ) = δ (β − α)∆ where δ is the Kronecker
Using the property that −∞
ν 1
t
2
ν0
ij
delta function, we can rewrite (3.5a) (modified with (4.8)) as

M
Z0 Ii = V i −

i−1

M
Z i−j I j −

j=1

i−1

M,c
Z i−j C j

(4.28b)

j=2

where
M
Z i−j

mn

= δ (t − (i − j) ∆t ) Sm (r) , LM
c {Sn (r) T (t)}
49

(4.28c)

.
M
M
LM
c {·} = αc /η0 Le {·} + (1 − αc )Lh {·}




(β − α)φm,n
M,c
0
Z i−j
=

mn 

0

(4.28d)

j≥γ
(4.28e)
otherwise

The recursive computation of the charge terms C j keeps the scaling of the solver at O(Nt Ns2 ).

4.7

Modification to Higher order space-time Galerkin scheme

It has been discussed earlier in this thesis that the stability properties of space-time
Galerkin methods is much more well understood than the properties of collocation based
methods. The discretization framework outlined in 3.2.1 should yield stable solutions. This
stability, however, is predicated upon an accurate discretization of the operators involved.
As discussed previously in this Chapter, different methods exist to enable the accurate
evaluation of retarded potential integrals. The quasi-exact integration methodology could
be used within this framework, but would require the development of analytical integration
rules for (p + 1) × (p + 2)/2 different integrals. Additionally, the restriction to low order
tessellations still applies. These considerations make the use of the separable expansion
method attractive. In this Section we outline detail the implementation of the separable
expansion method within the higher order space-time Galerkin framework of the previous
Chapter.

4.7.1

Interpolation properties

As shown in the previous Chapter, higher order convergence of interpolation error can be
achieved by expanding the unknown quantity using multiple basis functions within a single
50

time step, and applying Galerkin testing to discretize the system. In 3.2.2, a modulated
Gaussian pulse was shifted by convolving with a delta function and interpolated with higher
order Lagrange polynomials. In this Section, we will quantify the interpolation error when
this shift is achieved using the truncated expansion (4.8) rather than a delta function. The
signal, which is approximated using the higher order basis set and shifted using (4.8), is
defined for a given truncation limit, M , as

M (t) =
g˜∆

Nt



p

g j,l T l (t − j∆t ) ∗ 

j=1 l=0
Nt p
=
j=1 l=0



M

aν Pν (k1 ∆ + k2 ) Pν (k1 t + k2 )

(4.29)

ν=0
g j,l

M

aν Pν (k1 ∆ + k2 ) T l (t − j∆t ) ∗ Pν (k1 t + k2 )

ν=0

For the first interpolation result, the interpolation test of 3.2.2 is repeated for various upper
limits of (4.8). Let the interpolation error for a given p be defined as in (3.13a) and define
the error using the separable expansion with upper limit M as

ErrorM =

gM
approximate − g analytic

(4.30a)

g analytic
M
M
M (N ∆ )
gM
g∆
t t
approximate = g˜∆ (∆t ) g˜∆ (2∆t ) ...˜

(4.30b)

The error (4.30a) is shown in figure 4.8 as a function of the truncation limit for different p
and with ksamp = 10 (as defined in 3.2.2). This is compared against the same error when
the shift is achieved by convolution with a delta function. For each value of ksamp , M is
increased until the error converged below that achieved in Section 3.2.2.

It is expected

that the error using (4.8) will converge to that using the delta shift and figure 4.8 shows

51

Figure 4.8: Interpolation errors for various upper limits of (4.8)

that the error indeed converges to this error (or the error is slightly better when (4.8) is
used). The error rapidly improves as M is increased, to a value well below that of the deltashifted interpolation. This is likely due to the use of a very high bandwidth interpolant to
interpolate a band-limited signal. If band-limited interpolants were used, we suspect the
error would decrease monotonically until it converged.
Next, we show the error as a function of the oversampling factor, ksamp , in figure 4.9.
The support of the Legendre polynomials is chosen to reflect what would be expected in a
typical method of moments discretization with edge lengths of approximately λ/10, i.e., the
support is proportional to ksamp . The error behaves similarly to 3.1. It may be noticed that
the error for p = 3 does not converge at the same rate for large values of ksamp . This is due
to the very large support of the Legendre polynomials. It was not possible to reproduce the
error from 3.2.2 for these cases as it required more harmonics than was practical.

Table

4.1 shows the number of harmonics required in (4.8) to reproduce the error produced using
52

Figure 4.9: Interpolation errors for separable expansion

ksamp
p
1
2
3

5

10

20

40

1/1/1
2/1/1
4/3/4

1/1/2
4/2/2
4/4/4

2/1/2
4/3/3
6/4/5

3/1/3
4/3/3
7/5/5

Table 4.1: Lower bound of M needed to reproduce interpolation accuracies (undifferentiated/1st derivative/integral)

a delta shift. The upper limit needed to reproduce the interpolation error is defined as

M (Error) = inf M ≥ 0 : ErrorM ≤ Error

(4.31)

This is shown in table 4.1 for the interpolation of the function, as well as for its derivative
and integral, as required for the TD-EFIE. It can be seen that the number of harmonics
required is numerically tractable.
53

4.7.2

Higher order space-time Galerkin discretization with separable expansion

Now that we have quantified the errors of the higher order temporal basis used in conjunction with the separable expansion, we present the modified space time Galerkin scheme
with the separable expansion. It can be shown that for a temporal basis function T l (t) and
testing function T k (t − i∆t ), the resulting matrix elements are equivalent to those obtained
via collocation with temporal basis function

ξ kl (t − i∆t ) =

∞
−∞

dt T k t

T l t − i∆t + t

(4.32)

Therefore, the matrix elements for the higher order space-time Galerkin scheme can be
derived by substituting ξ kl (t) for the temporal basis in the collocation scheme detailed in
Section 4.6. The vectors V i and I i are as defined in (3.9c) and (3.9b).

54

Chapter 5
Numerical Results
In this Chapter we will demonstrate, via numerical experiment, the stability and accuracy
properties of the separable expansion method as applied to the TD-EFIE. It is well known
that the TD-EFIE is much more prone to instability than the TD-MFIE, which makes it
a better measure of the effectiveness of the separable expansion method. These numerical
experiments will be performed in two discretization settings: (1) Low order geometric discretizations, analyzed using collocation in time with 1st order temporal basis functions of 3.1,
and (2) higher order geometric discretizations in tandem with the higher order space-time
Galerkin 3.2. For type (1) discretizations we will present the following:

1. In Section 5.1, we will derive an eigenvalue stability analysis for the undifferentiated
TD-EFIE. This is similar to the analysis presented in [38], but the time integral in the
scalar potential changes the form of the companion matrix. This analysis will be a
useful tool in analyzing various problems. The relationship between these eigenvalues
and the stability properties of a given problem will be discussed.

2. In the ensuing Sections we will perform scattering simulations from a variety of geometries to demonstrate stability and accuracy. The former will be demonstrated either
by the aforementioned eigenvalue analysis or by observing the magnitude of a current
coefficient well after the excitation has vanished. The latter will be demonstrated for a
subset of the problems via RCS comparisons with a validated FDIE solver or analytical
results where available.
55

For type (2) discretizations we will present
1. Convergence studies of the higher order temporal basis via scattering simulations from
plates and spheres.
2. A result to show the flexibility of the temporal discretization framework by presenting
multiple temporal basis function types, while performing Galerkin testing and keeping
the support of the temporal basis functions within one time step.
3. A preliminary result showing that higher order spatial basis functions can also be used
with the separable expansion method without sacrificing stability.
4. A number of other results demonstrating stability and accuracy on more challenging
targets using the same metrics as were used for type (1) discretizations.

5.1

Eigenvalue stability analysis

The most robust method for post-hoc stability analysis of an MOT solution of a problem
is to examine the eigenvalue spectrum of the marching system. The magnitude of the largest
eigenvalues-whether they lie inside, outside, or on the unit circle in the complex plane-reveal
whether the magnitude of the solution can be expected to decay, grow, or remain bounded
by a constant value (marginal stability) as the solution time continues towards infinity. A
method to perform such analysis was presented in [38] for a variety of IE discretizations,
including the time-differentiated EFIE. In this Section, we look to validate the stability of
the separable expansion procedure as applied to both the temporal basis/testing scheme
of 3.1 with hat functions as the temporal basis and that of 3.2 with various basis function
orders. The major difference between the analysis performed here is that the undifferentiated
56

form of the TD-EFIE is being analyzed. The modification to the procedure is not trivial
because, as discussed in 4.6.1, the infinite tail of the scalar potential needs to be properly
handled to maintain the computational scaling of O(Nt Ns 2 ). This procedure modifies the
companion matrix of [38]. We now detail the modified eigenvalue analysis. This derivation
can also be found in [39].
Consider equation (4.28b), which we rewrite here (changing summation variables) as

M
Z0 Ii = V i −

i−1

M
Z j I i−j −

j=1

i−1

M,c
Z j C i−j

(5.1)

j=2

We first note that the marching matrix








M 

Z =









M
Z0
M
Z1


0

. . .

0




M
Z0
0 . .
0 




.




.




.

M
M
M 
Z N Z N −1 . . . Z 0
t
t

(5.2)

M
is banded from below, i.e. Z i = 0 for i > Nmax , where Nmax is, approximately, the
maximum delay in time between any source and observer point on the discretized scatterer.
More precisely, given the maximum spatial extent of the scatterer, diam(Ω),

Nmax ≤ diam(Ω)/(c∆t ) + nT

(5.3)

where nT << diam(Ω)/(c∆t ) depends on the temporal basis/testing functions being used.
57

Therefore, (5.1) can be rewritten as

M
Z0 Ii = V i −

min(Nmax ,i−1)

M
Z j I i−j −

min(Nmax ,i−1)

j=1

M,c
Z j C i−j

(5.4)

j=2

Combining (5.4) with (4.28a), which relates the current coefficients I j with the charge coefficients C j , the system can be rewritten as

AI i = V i − BI i−1

(5.5a)

where




 A11 0 
A=

A21 I



(5.5b)



 B11 B12 
B=

B21 B22

A11 =

(5.5c)





...
 − [Z0 ] 0



0
[I] 0
...




0
0 [I] 0 . . .




.


...
[I]








 ,







58

(5.5d)





 − Ti



0



 ...


...

A21 =

0
0



... 


... 

 ,

... 


0

ZP

 [Z1 ] [Z2 ] . . .


 [I]
0



 0
[I]



 .



 .


...
[I]

B11 =

...
...

0










 ,









(5.5f)




 Z˜1


 [I]



 [I]

B12 = 

 .



 .


...


B21 =

(5.5e)

Z˜2

Z˜P

...

0

...

0

...

[I]

0


0 ... ... 


0
... 



 ,






0 0

 Ti−1



0




.




.


...

59










 ,









(5.5g)

(5.5h)



















B22 =

.
Ti =

[I] 0
[I] 0
.
.
...

... 


... 



 ,






[I] 0

(5.5i)

dt Ti (t ) ,

(5.5j)

i∆t
(i−p−1)∆t









Ii = 







I1
I2
.
.
I min(N
max ,i)









Vi = 























(5.5k)


V1
V2
.
.
V min(N
max ,i)
















(5.5l)

and I is the identity matrix. If the matrix A is invertible, then the vector of current coefficients is given as

I1

=

−1
A V 1 − CI 0

(5.6a)

I2

=

−1
−1
A V 2 − C A V 1 − CI 0

(5.6b)

60

I3

=

−1
−1
−1
A V 3 − CA V 2 + C 2 A V 1 − C 3 I 0

(5.6c)

...
...
Ii

=

−1
A Vi +

i−1

−1
(−1)k C k A V k + (−1)j C j I 0 .
k=0

. −1
where C = A B Using the definition from [40], the eigenvalue decomposition of C is given
as
.
C=

νq † λq νq ,

(5.7a)

q
where † is the conjugate transpose. The eigenvalues λq are assumed to be sorted in monotonically decreasing order in magnitude. It then follows that

.
Ck =

νq † λk
q νq ,

(5.7b)

q
This leads to the following bound of a matrix vector product with C k

C k V ≤ λk
0

νq † νq V

(5.8)

q

which in turn leads to the bound on the current coefficients

−1
Ii ≤ A Vi +

i−1

λk
0

k=0
+ λi0

q

q

−1
νq † νq A V k

(5.9)

νq † νq I 0

This bound can be further reduced if it is assumed that the incident field vanishes at some
61

time I∆t . For i > I,
−1
I i ≤ I|λ0 |I maxk A V k

q

νq † νq + |λ0 |i

q

νq † νq I 0

(5.10)

There are three scenarios that result from (5.10).
1. |λ0 | < 1: The solution is stable and will approach 0 as i → ∞
2. |λ0 | > 1: The solution is unstable and will grow exponentially towards ∞ as i → ∞
3. |λ0 | = 1: The solution is marginally stable. The magnitude of the solution will be
bounded by a constant as i → ∞.
The analysis of this Section will be applied to a number of the scattering simulations to
follow later in this Chapter. While this is the ideal method of post-hoc stability analysis,
it is limited in the size of problem it can analyze. As the method relies on a singular
value decomposition of the companion matrix, it scales as O(Nc3 ), where Nc × Nc are the
dimensions of the companion matrix. Therefore, only a subset of the ensuing simulations
will be analyzed by this method. We will demonstrate the stability characteristics of the
other results by plotting a current coefficient well beyond the time at which the excitation
has vanished.

5.2

Low Order Discretizations

In this Section we examine stability for low order discretizations and solutions. By
low order we refer to the geometric discretization (flat triangles) and the spatio-temporal
basis/testing scheme employed. For each simulation, the EFIE is discretized with the ba62

sis/testing procedure of 3.1 and the separable expansion method. The excitation is an
incident plane wave of the form
2
2
ˆ
E inc (r, t) = uˆ cos(2πf0 t)e−(t−r·k/c−tp ) /2σ

(5.11)

where uˆ is the electric field polarization, kˆ the direction of incidence, and f0 the center
frequency. The values σ and tp are calculated as σ = 3/(2πB) and tp = 6σ, where B
is the bandwidth of E inc (r, t). fmax = f0 + B is defined as the frequency at which the
incident field power is approximately 40 dB below its peak. The time step is given as
∆t = χ/(20(f0 + B)) where χ > 0 is some real number. Each simulation is run for many
transits across the scatterer to show the absence of late time instability. Very broadband
excitations are used in each example.

5.2.1

Sphere

Our first result is scattering from a 1 m radius sphere of 864 spatial unknowns. The
incident field has parameters kˆ = zˆ, uˆ = xˆ, f0 = 90.9 MHz, and B = 90.89 MHz. The
simulation is run using two different time steps corresponding to χ = {4/3, 4}. The current
coefficient of one basis function is shown in figure 5.1 and is stable throughout the simulation.
As can be seen in the inset, the response due to the incident field decays very early on in the
simulation. To demonstrate the absence of instability we have run the simulation for over
2,500 transits across the geometry. In figure 5.2, we show the eigenvalues of the companion
matrix. As discussed in 5.1, eigenvalues within the unit circle correspond to decaying modes,
while those on the circle correspond to marginally stable modes. It is seen that all eigenvalues
have a magnitude below or equal to unity, confirming the stability suggested by figure 5.1.
63

Figure 5.1: Current on sphere
The only eigenvalues with unity magnitude lie on the positive real line, corresponding to the
DC solenoidal currents that reside in the null space of the Le operator To show that the
result is valid as well as stable, the RCS was computed in the x − z plane and compared
with the Mie series solution in figure 5.3. Some disagreement is seen, which we attribute
to the relative coarseness of the geometric discretization and the use of low order temporal
basis functions.

5.2.2

Parallel plates

To test the method on an open geometry, we simulate scattering from two parallel plates.
The plates are each 1 × 1 m in dimension and are separated by 0.1 m. The overall geometry
has 1120 spatial unknowns. The incident field has a center frequency f0 = 132.1 MHz,
√
bandwidth B = 132 MHz, incident direction kˆ = 0.5(ˆ
x + yˆ + 2ˆ
z ), and polarization kˆ =
√
0.5(ˆ
x + yˆ − 2ˆ
z ). The simulation is run for 4,000 transits across the geometry using time
64

Figure 5.2: Eigenvalue analysis of sphere
steps of χ = {1, 2}. The solution can be seen in figure 5.4 to remain stable throughout the
duration and there is no evidence of any late time instability. In fact, the current coefficient
decays below 10−25 in magnitude, which is similar to what has been observed in simulations
of other open geometries using this method.

5.2.3

Thin box

Our final low order result is a thin box of dimensions 0.5 × 1 × 0.1 m, which is discretized
with 390 unknowns. The geometry is quite thin and has sharp corners, which makes it a
65

Figure 5.3: RCS of sphere

Figure 5.4: Current on parallel plates

more challenging problem to stabilize than the sphere or parallel plates. The parameters of
the incident field are f0 = 106 MHz, B = 105 MHz, kˆ = −ˆ
y , and uˆ = zˆ. Simulating using
time steps χ = {0.5, 1}, the current is found for 1,589 transits across the geometry. The
66

solution can be seen in figure 5.5 to be stable. Also, an eigenvalue analysis was performed
and the results are shown in figure 5.11a. Aside from the eigenvalue λi = 1 + 0i, |λi | < 1 for
all i.

Figure 5.5: Current on box

5.2.4

Rectangular Cavity

For our final low order result we examine the stability of a very challenging problem, i.e.
a rectangular PEC cavity. The cavity is of dimensions 2 × 1 × 0.5m and is discretized to
9,450 unknowns. The excitation is a current dipole placed in the center of the cavity. The
temporal signature of the excitation is a modulated Gaussian pulse, as with the plane wave
excitations. For this simulation the excitation frequencies are set to f0 =300 MHz, and B =
150 MHz. The simulation is performed for 27,000 time steps (378 transits). Figure 5.7 shows
the current for the duration of the simulation.

The plot has been magnified so as to show

behavior of the solution, which is remarkably stable. There also appears to be no loss in the
67

2m

0.5 m
1m

Figure 5.6: Rectangular cavity

Figure 5.7: Current on rectangular cavity

energy of the solution. A longer simulation or an eigenvalue analysis should be performed
to confirm this.
To validate this result, the current has been transformed to frequency domain (sampled
every 3 kHz) to check against the analytical modes. The spectral content is shown in figure
5.8 with the modes labeled (ignoring degeneracies). Denoting the supported wave numbers
in the cavity as kx = mc π/2, ky = nc π, and kz = 2pc π, the analytical modes (mc , nc , pc )
and their associated frequencies are listed in table 5.1.

Aside from the small (1,1,1) mode,

every resonant frequency predicted by the TDIE solution is within 1 MHz of the analytical
68

Figure 5.8: Modes of rectangular cavity

frequency.

5.3

Higher Order Discretizations

In this Section we present results for higher order discretizations. For most of these
results, this means higher order in geometric discretization and/or temporal basis only. One
notable exception is a result for a higher order spatial basis. In general, however, the results
will use RWG basis functions. Unless otherwise stated, the geometric discretizations will use
2nd order curvilinear triangular elements. To achieve higher order accuracy in time, we use
the new temporal basis and testing scheme, which was presented in subsection 3.2.1. For
each result, we attempt to stabilize solution to the more challenging TD-EFIE. Additionally,
the accuracy of the results later in this Section will be quantified by comparing the error in
69

Mode (mc , nc , pc )
(1,1,0)
(2,1,0)
(3,1,0)
(0,0,1)
(1,0,1)
(0,1,1)
(1,1,1)
(2,1,1)
(3,0,1)
(3,1,1)
(4,0,1)
(1,2,1)
(4,1,1)

fexact
167.59
211.99
270.23
299.80
309.02
335.18
343.46
367.17
374.75
403.61
423.98
430.55
449.70

ftdie
167.50
211.85
270.07
300.00
308.93
334.93
339.49
366.93
374.34
403.55
423.44
430.00
449.33

∆f
0.10
0.14
0.16
0.19
0.09
0.25
3.98
0.25
0.41
0.06
0.54
0.56
0.38

Table 5.1: Rectangular cavity resonant frequencies (all frequencies in MHz)

RCS to that of an FD-EFIE code. This error will be defined as

Error =

ζtdie − ζfdie
ζfdie

(5.12)

where RCStdie,fdie = 10log10 (ζtdie,fdie ) is the radar cross section obtained using the TDEFIE or FD-EFIE. Except where otherwise noted, time steps corresponding to χ = 1 and
temporal basis orders of p = 2 are used.

5.3.1

Convergence in far field scattering

Our first results look at convergence in far field scattering with respect to p and ksamp .
The first geometry is a 1×1 m plate, meshed with 200 triangular elements. The incident field
is a plane wave of f0 = 150 MHz, fmax = 225 MHz, kˆ = zˆ, and uˆ = xˆ. This convergence
study was performed by simulating for p = 1, 2, and 3 and ksamp = 5, 10, 20, and 40 (as
70

defined in Section 3.2.2. The convergence shown in figure 5.9 is defined as.
|E θ k+1 (θ, φ) − E θ k (θ, φ)|
Convergencek+1 =
|E θ k+1 (θ, φ)|

(5.13)

where the scattered field is approximated for the kth value of ksamp (ksamp (k + 1) >
ksamp (k), k = 1, 2, 3) as

ˆ φ (θ, φ) + θE
ˆ θ (θ, φ)
Es (θ, φ) ≈ φE
k
k

(5.14)

The far field approximations for E φ and E θ are used. For this result, θ = −130◦ and

Figure 5.9: Far field convergence for plate
φ = 0◦ . The convergence rates were nonuniform across the various θ values in the x − z
p+1
p+2
plane. The range of these rates were roughly between O(∆t ) and O(∆t ).
Our second result performs the same experiment, but for a 1 m radius sphere. This geometry is meshed with 384, 2nd order triangular elements. The incident wave has parameters
71

f0 = 60 MHz, fmax = 90 MHz, uˆ = xˆ, and kˆ = zˆ with the far field computed at θ = −180◦
and φ = 0◦ . For this result, we solved the TD-CFIE with the combined field parameter
αc = 0.5.

The result is shown in figure 5.10. As in the case of the plate, the rates of

Figure 5.10: Far field convergence for sphere
p+1
convergence were nonuniform in the x − z plane, but were in the range between O(∆t )
p+2
and O(∆t ).

5.3.2

Stability analysis for thin box

To characterize the stability properties of the higher order space-time Galerkin scheme,
the eigenvalue analysis of 5.1 is performed for the thin box of 5.2.3. Lagrange polynomials
of p =0, 1, and 2 are used for the temporal basis. A subset of the eigenvalues (those close to
unity magnitude) are shown in figure 5.11 and compared with those of the result from 5.2.3.
The eigenvalues for p = 0 are roughly equal to those for the collocation/hat function scheme.
This is expected as the matrix elements for these two schemes can be shown to be equivalent.
72

(a)

(b)

(c)

(d)

Figure 5.11: Eigenvalue analysis of thin box discretized with (a) 1st order basis/delta testing,
(b) 0th order basis/Galerkin testing, (c) 1st order basis/Galerkin testing, (d) 2nd order
basis/Galerkin testing

For the higher values of p, the eigenvalues corresponding to resonant modes approach unity
in magnitude, although they remain within the unit circle. We have observed in practice,
that if the expansion (4.8) has not converged or integrals are not accurately evaluated, these
73

eigenvalues will migrate outside of the unit circle and the solution will be unstable. The
magnitude of the eigenvalue at 1 + 0i was seen to be unaffected by the order of the basis
functions used in higher order space-time Galerkin scheme.

5.3.3

Sphere with various spatio-temporal bases

Our next result shows scattering from a sphere, but for a variety of spatio-temporal
discretizations of the surface current. For this result, we look to accomplish the following:
1. Show the flexibility in the choice of temporal basis when using the scheme of 3.2.1.
2. Demonstrate that higher order spatial basis functions can also be used with the spacetime separation to generate stable solutions.
3. Show that the separable expansion method can be used in the temporal discretization
framework of [34] to produce stable results.
We use the same excitation and geometric discretization as those for the sphere in 5.3.1. To
accomplish the first task, 2nd order Legendre polynomials are used as the temporal basis
functions and compared with the result using 2nd order Lagrange polynomials. The results
lie on top of each other, decaying to the same DC value. For the 2nd task, 2nd order divergence conforming GWP functions [41] are used as the spatial basis with 2nd order Lagrange
polynomials as the temporal basis. Again, the result is stable throughout the simulation.
The inclusion of higher order spatial basis functions presents no added complications to the
separable expansion method. The only caveat is that the upper limit of (4.8) may be affected
by the type of spatial basis function used. The value of M required for must be studied on
a case by case basis. For task 3, the basis/testing scheme of [34] is stabilized using the
74

separable expansion. A linear growth term is observed, which corresponds to the solenoidal
currents with linear time dependence, which lie in the null space of the Le operator. It
is often preferred to discretize the time-differentiated TD-EFIE to avoid the time integral
operator in the scalar potential. In this setting, the basis functions of [34] can be used to
achieve higher order accuracy in time and, as demonstrated here, the separable expansion
can be used to preclude any high frequency instability from appearing.

Figure 5.12: Current on sphere using various spatio-temporal basis functions

5.3.4

Cone-sphere

Our next result is scattering from a cone-sphere of 7,965 unknowns. The spherical portion
of the geometry has a radius of 1 m and the height of the cone is 4 m. The excitation has
f0 = 80 MHz, B = 70 MHz, kˆ = xˆ, and uˆ = yˆ. Both the TD-EFIE and TD-CFIE (αc = 0.5)
are solved for 80 transits across the geometry. This is a challenging problem due to the sharp
tip on the cone and the small edges required to discretize it, yet the current can be seen in
75

figure 5.13 to remain stable throughout the simulation. To demonstrate the accuracy of this

Figure 5.13: Current on cone-sphere

problem, the RCS for the TD-EFIE solution is computed in the x − y plane and compared
against a validated, higher order FDIE code. The results are shown in figure 5.14 and match
very well. The error in far field is computed at 11, 80, and 149 MHz to be 1.8%, 1.5%, and
0.09%, respectively.

5.3.5

Almond

Our next result is a thin almond of 6,642 unknowns. The almond fits in a box of dimension
5×4.3×0.87 m. The incident field has f0 = 35 MHz, B = 25 MHz, kˆ = zˆ, and uˆ = xˆ (towards
the tip of the almond). The challenge of this problem is the thinness of the geometry and
the sharp tip. The current is computed for 2,000 time steps (100 transits) using both the
TD-EFIE and TD-CFIE (αc = 0.5) and the magnitude of one current coefficient is shown in
figure 5.15. The current is stable throughout the simulation. The RCS is computed for the
76

Figure 5.14: RCS of cone-sphere

Figure 5.15: Current on almond
TD-EFIE solution in the x − z plane and compared against an FDIE solver. The result can
be seen in figure 5.16 to match very well.

The errors in far field are computed at f =11,
77

Figure 5.16: RCS of almond
40, and 69 MHz and found to be 0.14%, 0.035%, and 0.051%.

5.3.6

VFY218

Our final result is scattering from a VFY218 discretized using 6,498 flat elements. The
scatterer fits in a box of dimensions 15 × 9 × 4 m. The axis from the tail to the nose of the
aircraft is in the x-direction, whereas that from wing to wing is in the y − direction. The
VFY218 is illuminated using a plane wave with parameters f0 = 60 MHz, B = 40 MHz,
kˆ = −ˆ
y (side on) and uˆ = xˆ. The solution to the TD-EFIE is computed for 100 transits
across the geometry (7,000 time steps), one spatial unknown of which is plotted in figure
5.17. This simulation is performed for the p = 0 and p = 1 temporal basis. This is a very
challenging problem to stabilize due to the sharp edges and corners of the geometry. It
is particularly challenging when the TD-EFIE is used, yet stability is seen throughout the
78

simulation.

As expected, the 0th order basis decays much more rapidly than does the 1st

Figure 5.17: Current on VFY218

order basis. To validate the solution, the RCS for the p = 1 solution is computed in the x − y
plane and compared to a frequency domain EFIE solver. As seen in 5.18, the two solutions
agree very well. The error between the two solutions at 21, 60, and 80 MHz was found to be
1.9%, 0.67%, and 2.1%, respectively, for p = 0 and 1.9%, 0.46%, and 0.57% for p = 1.

5.4

Summary

In this Chapter we have presented a variety of scattering results to show that the separable expansion method can be applied to the TD-EFIE to yield stable results, without
sacrificing accuracy. The claims made earlier in this thesis that the separable expansion
method can be used to stabilize solutions on higher order geometric discretizations has been
justified by examining a number of complex and challenging targets. Additionally, the accu79

Figure 5.18: RCS of VFY218
racy and stability properties of the higher order space-time Galerkin framework of Section
3.2 have been quantified. The rest of this thesis will be devoted to the acceleration of systems discretized using the separable expansion method, and the extension to more complex
problems and discretization frameworks.

80

Chapter 6
Acceleration using Multilevel Plane Wave Time
Domain Algorithm
6.1

Background

Discretized IEs can be quite challenging to solve efficiently due to the dense makeup of
their system matrices. This characteristic is a result of the global nature of the Green’s
function. In the frequency domain, the system matrices are completely filled. In the time
domain, the entries are spread across a number of matrices due to the non-zero propagation
time from a given source to a given observer. The number of matrices depends on the
time step size, dimensions of the scatterer geometry, and the material properties. These
differences aside, the number of nonzero entries for both FDIEs and TDIEs scales as the
number of spatial unknowns squared. This necessitates algorithms to reduce this scaling if
many problems of practical interest are to be analyzed.
Fortunately, such algorithms are well developed for both FDIEs and TDIEs. Some of
the most widely used acceleration algorithms for FDIEs are the adaptive cross approximation algorithm (ACA) [42]; FFT-based algorithms, such as the adaptive integral method
(AIM) [8]; the fast multipole method (FMM) [43]; and the method of accelerated cartesian
expansions (ACE) [44]. ACA relies on the low rank properties of asymptotically smooth
functions. The algorithm works by hierarchically partitioning the computational domain
into clusters of sources and observers. Well separated clusters can then be represented by a
81

low rank approximation. This results in a system matrix with many low-rank blocks. The
acceleration is achieved through multiplication with these low rank representations, which
4/3
leads to a scaling of O(Ns log(Ns )) [45]. As an algebraic method, this algorithm is kernel
independent. In AIM, a uniformly spaced, auxiliary grid is introduced, onto which source
data is projected and from which observer data is projected. This auxiliary grid creates a
three-level Toeplitz interaction matrix (between grid points), which can be accelerated via
FFT. For surface integral equations, the number of auxiliary grid points, Ng , can be much
greater than Ns , which limits its effectiveness for such problems. The algorithm scales as
3/2
O(Ns log(Ns )) for SIEs. On the other hand, the method is very well suited for volume
integral equation solutions and scales as O(Ns log(Ns )). FMM relies on the analytical properties of the Green’s function to achieve separation between sources and observers. Again,
interactions are separated into near and far interactions. Far interactions are computed using
a truncated, approximated representation of the Green’s function. Although FMM is inexact, the error is controllable, so this does not limit its utility. However, as it is an analytical
rather than an algebraic method, it is kernel dependent. Lastly, ACE is similar to FMM,
but uses Taylor expansions rather than the analytical properties of the Green’s function to
approximate interactions. This has the benefit of making the method kernel independent,
but the method is best suited for non-oscillatory kernels. The attractive characteristics at
lower frequencies makes its hybridization with FMM a very powerful methodology.

Each of these methods have analogues in the time domain [5, 6, 46, 47]. That of ACA
[46] is based on a marching on in degree solution procedure. In this work, the authors
observed that the matrices in this discretization scheme are also low rank, allowing ACA
compression to be used. Again, the method is purely algebraic, making it kernel independent.
82

Time domain AIM also shares a number of features with its frequency domain counterpart
[6]. As with frequency domain AIM, a uniform auxiliary grid is used to create a Toeplitz
matrix structure. The main difference being that the Green’s function is now translationally
invariant in time, yielding matrices with a four level Toeplitz structure. These are accelerated
with FFTs, resulting in scalings of O(Nt Ns 3/2 log2 (Ns )) and O(Nt Ns log2 (Ns )) for SIEs
(for surfaces that are not quasi-planar) and VIEs, respectively. The plane wave time domain
algorithm (PWTD) is the time domain analogue of FMM [5]. FMM is based on expansions
in multipoles, whereas PWTD relies on plane wave expansions. As in FMM, interactions
in PWTD are characterized as being either near or far field, but in the time domain an
additional time gating procedure is needed to maintain causality. The algorithm reduces
the complexity to O(Nt Ns log2 Ns ). Details of the algorithm will be presented later in
this Chapter. Lastly, the time domain version of ACE was developed in [47]. As in the
frequency domain, it is based on Taylor expansions and kernel independent. Analogous to
its effectiveness in hybridizing FMM in the frequency domain, the time domain variant can
be effectively used with PWTD in a multiple time step marching scheme.
In this thesis we will utilize PWTD as our acceleration methodology. We will implement this method with the separable expansion, by only using the expansion for near field
interactions. Far field interactions will be implemented as prescribed in [48], but for the
undifferentiated TD-CFIE. This will be justified later in this chapter. We will use the higher
order spatial basis described in 3.2.1, in conjunction with PWTD as well.
This chapter is organized as follows. In section 6.2, we will detail the multilevel PWTD
algorithm. In section 6.2.4, we will elaborate on the use of PWTD with the separable
expansion, and also with the higher order temporal basis set introduced in 3.2.1. Lastly, in
83

section 6.3 we will present results to verify the stability and computational scaling of the
algorithm.

6.2

Multilevel PWTD

In this section we detail the multilevel PWTD algorithm. First we will look at a toy
problem to elucidate the mechanism upon which the PWTD algorithm depends. Next,
we will look at the TD-CFIE itself. Conventional PWTD implementations accelerate the
solution of the time-differentiated TD-EFIE and -MFIE, while we seek to accelerate their
undifferentiated forms. To this end, we will discuss a numerical integration scheme based
on Runge Kutta methods applied to the PWTD algorithm. Then we will consider how this
algorithm is used in tandem with the space-time separation of chapter 4. Lastly, we will
discuss the use of the higher order temporal basis functions of section 3.2.1 within the PWTD
algorithm.

6.2.1

Motivation

.
To begin, consider a set of source and observation points (assumed coincident) R =
{r1 , r2 , ..., rN }, ri ∈ S, where S is the surface of some scatterer. Assume there exist point
sources that are approximately band-limited and of limited time duration given as




f i f (t) r = ri , 0 ≤ t ≤ Ts
0
fi (r, t) =



0
otherwise
84

i = 1, 2, ..., N

(6.1)

The potential due to one source is given as

δ(t − |r|/c)
∗st fj (t, r)
4π|r|
r=ri
∞
δ(t − t − |ri − r |/c)
1
dt
=
dr
fj (t , r )
4π S
|ri − r |)
−∞
fj (t − |ri − rj |/c, rj )
=
4π|ri − rj |

φj (t, ri ) =

(6.2)

We wish to compute this for i, j = 1, 2, ..., N , and for some set of time samples t =
∆t , 2∆t , ..., Nt ∆t . This is an O(Nt Ns 2 ) process. Now consider a source and observer
sphere, which enclose source point rj and observation point ri , respectively (see Figure 6.1).
Assuming the source fj is no longer radiating, i.e. t > Ts , (6.2) can be rewritten as

Source Sphere

Observer Sphere

rsc

ri

roc

rj

Figure 6.1: Source and observer plane wave expansions

2π
π
1
kˆ · R
φj (t, ri ) =
dφ
sin θdθ fj (t, r) ∗st δ t −
2πc 0
c
0
=

1
2πc

d2 Ωδ t −

kˆ · (r − roc )
c
85

r=ri
kˆ · (roc − rsc )
∗t δ t −
c

(6.3)

∗t

=

1
2πc

dr δ t −
S

d2 Ωδ t −

kˆ · (rsc − r )
c

kˆ · (r − roc )
c

∗t fj (t, r )
r=ri

∗t T (t, roc − rsc )

where R = r−r . As can be seen in the last line of (6.3), there are no longer any dependencies
on R, but only on source and observer locations with respect to the center of their bounding
spheres. As mentioned earlier, the sources fj (t, r) should be band-limited and have finite
support in time. This is not possible to realize exactly, but can be approximately satisfied
through proper choice of basis functions. As is discussed in greater depth in [48], approximate
prolate spheroidal wave functions can be used in this capacity to great effect. In the MOT
solution, the surface current is expressed as a linear combination of Nseg currents of finite
temporal support, Ts .
Nseg

Ns
ˆ (r, t) =
J

Sn (r)
n=1

j
Where supp fn (t)

j
fn (t)

(6.4)

j=1

= [(j − 1)Ts , jTs ]. With the current written in this way, computing

the interaction between two well-separated clusters using the two-level PWTD algorithm
can be summarized as a three step process: (1) aggregating data from the source cluster
onto outgoing plane wave directions, (2) translating plane waves from the source sphere onto
the observer sphere, and (3) disaggregating from plane waves back onto the points in the
observer cluster. This is performed for all sufficiently separated clusters. The remainder
of the interactions are computed directly. The two level algorithm achieves a scaling of
O(Nt Ns 1.5 log(Ns )).

86

6.2.2

PWTD Algorithm

Before proceeding to describe the multilevel PWTD algorithm, we must introduce an
expression analogous to (6.3), but for the method of moments discretization of the time
differentiated TD-CFIE (the undifferentiated version will be discussed in 6.2.3). This is
given as

δ (t − (i − j) ∆t ) Sm (r) ,
1
8π 2 c2

M

M

∂
Lc {Sn (r) T (t)}
∂t

≈

(6.5)

o (kˆ , t, kˆ ) + (1 − α )S o (kˆ , t, n
akl αc /η0 Sm
c m kl ˆ )
kl
kl

l=0 k=−M

∗t TM (kˆkl , t, roc − rsc ) ∗t Sns (kˆkl , t, kˆkl ) ∗t T (t)

where

o,s ˆ
Sm (k,
t, n
ˆ) =
ˆ t, R) =
TM (k,

S

dr n
ˆ × Sm r

c∂t3 M
(2p + 1)Pp
2R
p=0

o,s
r − rc
ˆ
δ t±k·
c
ct
R

Pp

kˆ · R
R

(6.6)

(6.7)

The (+) and (−) in (6.6) are taken for source and observer, respectively, and Pp (·) are pth
order Legendre polynomials.
We will now summarize the steps of the multilevel PWTD algorithm. First, we require
that the computational domain be recursively subdivided into cubes. The largest cube (root
box) is subdivided into 8 smaller cubes until the side length of the smallest cube (leaf box)
is on the order of one wavelength at the highest frequency. Given a cube, the 8 cubes into
which it is divided are referred to as the children, whereas the original cube is the parent. In
87

source box

far field

near field
S

Figure 6.2: Decomposition of 2 dimensional scattering geometry

what follows, we assume that the root box has been recursively subdivided down to the leaf
level leading to Nl total levels. In addition to this partition of the domain, we require that
a criterion is specified which determines whether two boxes are in the near or far field of
one another. Let α (l) and α(l) represent a source and observer box at level l, respectively,
where l = 1 represents the leaf level and l = Nl represents the root level. Two boxes, with
α (l)
α(l)
centers rc
and rc , and bounded by spheres of radii Rα (l) = Rα(l) = R(l), are said
to be in each other’s far field at level l if (1) their parents are in the near field of each
α(l)
α (l)
other and (2) R{α(l), α (l)} = |rc
− rc
| > γR(l), where 4 < γ < 6. This criterion is
slightly different from that of MLFMA, in which γ = 2. The larger γ for PWTD is explained
as follows: Given the three step process of aggregation of source data, translation of plane
waves, and disaggregation onto observers, it can be shown that the field seen at the observers
consists of a superposition of two fields. The first is the true observer field and the second is
what is termed the ”ghost signal”. This ghost signal can be effectively time-gated out of the
total observed field if the following condition is met. Given that the ghost signal vanishes at
observation point r at time tghost and the true signal appears at r at time ttrue , the true
88

field can be extracted from the total field if ttrue > tghost . This is a consequence of the
temporal locality of the source signals. This condition is equivalent to

cTs < R{α(l), α (l)} − 2R(l)

(6.8)

(6.8) shows that the distance between source and observer spheres must satisfy an additional
distance criterion of cTs due in order to time gate the ghost signal. For details of this
derivation, the reader is referred to [49].
With these definitions in hand, the Multilevel PWTD algorithm can now be summarized.
1. For all source leaf boxes, α (1), compute the response at all observer leaf boxes, α(1),
in its near field directly using (3.5b).
2. For all leaf boxes, compute plane wave expansions from the sources using (6.6)
3. At each level l = 2, ..., Nl − 2, compute the plane wave expansion of α (l) using the
expansions of its children. As discussed in [48], this procedure is equivalent to the
evaluation of the far field scattering pattern of the parent box from the far field patterns
of its children. This requires a finer sampling of plane wave directions. The sampling
procedure used in this thesis is based on a scalar spherical filter [50]. The temporal
width of each ray also doubles for each higher level.
4. Compute the translations from the plane wave expansions of source boxes to the local
expansions of boxes in their far fields. This must be done for all boxes on all levels.
5. Compute the local expansions at each level l = 1, ..., Nl − 3 from those at the higher
level. Arriving at the local expansion of a child box from its parent is the dual of
89

step (3). The process of coarsening the sampling of plane wave directions is called
”anterpolation”.

6. For all leaf boxes, compute the fields at observers from the plane wave expansions using
equation (6.6).

For step 4, the convolutions between the outgoing plane wave expansions and the translation operator are performed in the frequency domain by taking an FFT of the plane wave
expansion and multiplying by the Fourier transform of the translation operator. The form
of the translation operator in the frequency domain is known. In fact, it is the same as the
regular part of the MLFMA translation operator. The result is then transformed back to
the time domain using an IFFT. The scaling of this entire process is O(Nt Ns log2 Ns ).

6.2.3

Undifferentiated TD-CFIE using Runge Kutta Integration

To accelerate the undifferentiated TD-CFIE, it is necessary to temporally integrate equation (6.5). This can be done numerically using Runge Kutta methods. We will now give
a brief overview of these schemes. Runge-Kutta methods are used to solve the ordinary
differential equation
y(t)
˙ = f (t, y(t))

(6.9)

Where the dot denotes a derivative with respect to t. A general ν-stage RK scheme is
characterized by its Butcher Tableaux

c

A
b

T

90

(6.10a)

Aij = aij

i, j = 1, 2, ..., ν

(6.10b)

b = [b1 b2 ... bν ]T

(6.10c)

c = [c1 c2 ... cν ]T

(6.10d)

For explicit RK solvers, c1 = 0, and aij = 0 for j ≥ i. Given a Butcher Tableaux for an
explicit RK method, the value of yn+1 = y(tn+1 ) is computed as
ν
yn+1 = yn + ∆t

bj kj

(6.11a)

j=1
k1 = f (tn , yn )

(6.11b)

k2 = f (tn + c2 ∆t , yn + a21 ∆t k1 )
.
.
.
ν−1
aνj kj )

kν = f (tn + cν ∆t , yn + ∆t
j=1

Within the specific setting of PWTD, y represents the scattered fields and f is the Lc
operator.

6.2.4

Formulation using separable expansion in higher order spacetime Galerkin framework

The PWTD algorithm is unchanged when using the separable expansion of chapter 4. For
near field interactions, as defined in the previous section, the procedure outlined in chapter
91

4 is used, whereas the expansion is not needed to accurately evaluate (6.5). This is true
because the temporal basis functions appearing in (6.5) are prolates, which are (approximately) smooth, band-limited functions. As such, their appearance in the spatial source
integrand does not prevent accurate integration via numerical quadrature, which is the motivation for the separable expansion. The part of the solution procedure associated with
near field interactions has a cost which is unaffected by the separable expansion. Therefore,
the computational scaling remains as O(Nt Ns log2 Ns ).

The only modification to the algorithm comes not from the separable expansion but from
the use of higher order temporal discretization framework from section 3.2.1. We will refer to
this representation as the polynomial representation as opposed to the prolate representation.
Given a temporal basis function order of p, there are a factor of p + 1 more unknowns for the
polynomial representation than for the prolate representation. It is therefore necessary to
add two more steps to the PWTD algorithm. The first is prior to the construction of outward
travelling rays and the second is during the projection of incoming rays onto observers.

The (approximate) equivalence between the two representations is given as
Nt

Ns
ˆ (r, t) =
J

Sn (r)

p

j,l
In T l (t − j∆t )

(6.12)

n=1

j=1 l=0
Nseg
Ns
j
≈
Sn (r)
fn (t)
n=1
j=1
Nt
Ns
j,prol
=
Sn (r)
In
T (t − j∆t )
n=1
j=1
Assuming the prolates, T (t), are interpolatory at t = −δ ∈ [0, ∆t ), delta testing both sides
92

at the interpolation points results in

j,prol
In
=

p

j,l
In T l (δ)

(6.13)

l=0
j,l
for j = 1, 2, ..., Nt , where In are known. Through interpolation tests, we have found that
choosing the interpolation point to lie in the center of a time step results in the highest
accuracy. To demonstrate this, the experiment of section 3.2.2 was repeated for ksamp = 10
and ∆ varied between 0 and 1. As shown in figure 6.4, the interpolation error is notably
better when δ = 0.5. The error even converges below the original error when the basis of
(3.8a) is used.

Figure 6.3: Interpolation errors for projections onto Prolates (δ = 0)

Lastly, because Galerkin testing is used in time, equation (6.5) must be modified to reflect
this. The temporal testing integral can be computed using one dimensional quadrature. The
step of projecting incoming rays onto observation points must now loop over these temporal
93

Figure 6.4: Interpolation errors for projections onto Prolates (δ = 0.5)
integration points.

6.3

Results

In this section we present PEC scattering results to validate the computational scaling
and soundness of the PWTD-accelerated algorithm with the separable expansion and spacetime Galerkin framework. To achieve this, we examine scattering from plates using the TDEFIE and spheres using the TD-CFIE. For both geometries, we look at various discretization
densities, incident wave frequencies, and simulation times. Later in the chapter we will look
at scattering from more challenging objects.

6.3.1

Stability and accuracy

We begin by characterizing the stability and accuracy of the results. The first results
are for scattering from plates of dimension 1 × 1 m. For each result, the excitation is a
94

√
√
uˆ = − 2/2(ˆ
x + zˆ)-polarized plane wave, incident in the kˆ = − 2/2(ˆ
x − zˆ) direction.
The frequencies of the excitation are set such that, for a given discretization, the largest
edge lengths are approximately λmin /10 where λmin = c/fmax . Each simulation lasts
for roughly 50 transits across the geometry. The precise input parameters are shown in
table 6.1. p = 0 temporal basis functions are used for all simulations. Figure 6.5 shows
the magnitude of one current coefficient through the simulation for each discretization and
excitation. A slowly growing, low frequency term can be seen in most of the results. It
does appear after the initial response decays, but this does suggest that the algorithm is
not robust for the TD-EFIE in its current implementation. Later in this chapter we will
discuss some possible remedies for stabilizing the EFIE. We do note here that this is a low
frequency, rather than a high frequency instability, so the problem does not appear to be
inaccurate spatial integration. For the rest of the results we will present TD-CFIE results
when PWTD is used and leave stabilization of the PWTD-accelerated TD-EFIE as future
work. At the end of this chapter we will discuss potential sources of instability.

Our next result is scattering from a 1 m sphere using the TD-CFIE (αc = 0.5). For
this result, the plane wave excitation is incident in the kˆ = zˆ direction and polarized with
uˆ = xˆ. As with the plate simulations, maximum frequencies are set to yield 10 edges per
wavelength. Again, each simulation lasts long enough such that a wave can travel 50 times
across the largest dimension of the geometry. The current in each result is discretized with
p = 2 temporal basis functions. Figure 6.6 shows the magnitude of a current coefficient on
the sphere and stability can be seen throughout the simulation. This result is validated via
comparison to the Mie series. The result is seen in figure 6.7 to agree well with the analytic
result.
95

Figure 6.5: Current on plate

Figure 6.6: Current on sphere

6.3.2

Timing results

For these same results, we now study the scaling of the algorithm. The total times of
the solution procedure are presented against the number of unknowns and compared to the
96

Figure 6.7: RCS of sphere
Ns
645
1,160
3,605
7,400
11,408
16,725

Nt
1,462
1,958
3,491
5,000
5,944
6,200

f0 (MHz)
185
215
385
600
700
900

B (MHz)
125
200
355
460
560
760

Nl
3
4
5
5
6
6

Table 6.1: Input parameters for plates

expected scaling of the algorithm.
We begin with the plate. Figure 6.8 shows the timing results for these simulations. The
total solution time is compared against the expected scaling of O(Nt Ns log2 Ns ). It can be
seen that the algorithm indeed scales as expected. When the higher order temporal basis
is used, the number of unknowns increases by a factor of p + 1. Assuming the order of the
basis set is fixed, this should not affect the scaling of the algorithm, which is consistent with
our observations.
For the second result, scattering from a sphere of radius 1 m is simulated. Again, 6
97

Figure 6.8: Timings of PWTD-accelerated plate simulations
Ns
576
1296
2304
3600
5184

Nt
1270
1800
2335
2870
3270

f0 (MHz)
100
140
180
220
250

B (MHz
90
130
170
210
240

Nl
1
2
2
2
3

Table 6.2: Input parameters for spheres
different discretizations are used, with the input parameters given in table 6.2

Figure 6.9

shows the scaling of the algorithm is exactly as expected.

6.4

More Challenging Structures

In this section we present scattering results for more complex structures than those of
the prior section. For each result we simulate using the TD-CFIE with αc = 0.5. The
geometries being illuminated have been studied in chapter 5 and we repeat the analysis
using the PWTD-accelerated solver. For each result, the incident wave parameters are
unchanged. RWG spatial basis functions and the temporal basis functions of 3.2.1 are used
for all simulations.
Our first result is scattering from the same cone-sphere as was shown for the direct
98

Figure 6.9: Timings of PWTD-accelerated sphere simulations

solution in 5.3.4. This simulation was performed with p = 2 temporal basis functions. The
current is shown in figure 6.10 for 2,400 time steps in comparison to the direct solution to
both the TD-EFIE and TD-CFIE Again, the results are seen to be stable throughout. It
can be seen that the early time of the current agrees well with the direct solution. The
PWTD solution to the TD-CFIE hits a floor several orders of magnitude of above that of
the direct solution. One possible explanation is that we are approximately integrating the
time-differentiated Lh operator. This form of the operator has a null space that includes
DC currents, which could perhaps be the source of the DC floor we observe. For the p = 2
solution, the RCS was computed in the x − y plane and compared with an FD-CFIE (αc =
0.5) solution. The result is shown in figure 6.11 and the results agree very well.
For our last PWTD-accelerated result, we look at the higher order almond of 5.3.5.
Again, the incident wave parameters are unchanged from the direct solution. The current
99

Figure 6.10: Current on cone-sphere compared to direct solution (p = 2)

Figure 6.11: RCS of cone-sphere
is discretized using p = 2 temporal basis functions. Figure 6.12 shows a current coefficient
for 1,800 time steps and the result is seen to be stable throughout. For this result, the DC
floor for the PWTD accelerated solution is much lower than for the previous result for the
100

Figure 6.12: Current on almond compared to direct solution (p = 2)
cone-sphere. For this result, we noted that there were no interactions above the leaf level,
whereas for the cone-sphere, there were interactions 1 level above the leaf level. This points
to either (1) some inaccuracies in the interpolation and anterpolation process or (2) not
properly handling the higher spectral content contained in the larger, level 2 boxes. Again,
to validate the result we compare the RCS, computed in the x − z plane, to an FD-CFIE
result. This is shown in figure 6.13 and the results agree quite well.

6.5

Discussion

As shown in the plate results of section 6.3.1, the results for the TD-EFIE have slowly
growing instabilities. There are a number of possible reasons why this may be so. One
reason is that the formulation relies on an RK scheme to recover the scattered field from its
derivative in time. Therefore, the field is only approximately recovered and the currents in
the null space of the ∂t Le operator may appear in the solution. We have also implemented
101

Figure 6.13: RCS of almond

PWTD with the separable expansion applied to the time differentiated TD-EFIE and have
seen a similar growing mode, which supports this explanation.
For more complicated scatterers, such as a NASA almond and a cone-sphere, we have
also seen high frequency instabilities for the TD-EFIE. When these scatterers are analyzed
without PWTD acceleration, the result is stable, so the instability must be introduced by
PWTD. To shed light on this, the spectral content of the scattered field due to an electric
dipole is shown in figure 6.14.

The current time signature is projected onto approximate

prolate spheroidal wave functions of various widths in time. It can be seen that if the time
width is not large enough, higher frequency content is present, which is outside of the band
of the excitation. This is because the prolate function is not exactly, but only approximately
time limited, and the windowing in time corresponds to a convolution with a sinc function
in frequency. It appears as though the time window may need to be widened so as not to
introduce spurious high frequency error, which can cause high frequency instability. More
102

Figure 6.14: Spectrum of scattered field
investigation is needed to verify that this approach can stabilize the solution.
One other possible source of instability is the way in which source dipoles are grouped
into bounding spheres. For a given basis function, every associated quadrature point will be
grouped together in a given box. A case may arise in which a triangle overlaps between two
boxes, in which case a dipole will be assigned to a sphere that it is actually external to. This
is shown pictorially in figure 6.15. The effect of such a case on the resulting scattered field
merits further study.

103

×
×

×

×

×

Figure 6.15: Grouping of dipoles

104

Chapter 7
Implementation for Penetrable Scatterers
7.1

Motivation

The previous chapters have been devoted to the development of a stabilization methodology for analyzing scattering from PEC structures. While this is a significant challenge for
reasons highlighted earlier in this thesis, in our experience the stable analysis of penetrable
scattering bodies is more challenging than for PEC structures. In [1], stable solutions were
presented for scattering from such bodies using quasi-exact integration, which is restricted
to flat tessellations. As with PEC scatterers, the accuracy of geometric discretizations can
be greatly improved by using curvilinear elements. This motivates the extension of the separable expansion method to penetrable scattering problems. In this Chapter, we outline
the procedure for applying the separable expansion method to analyzing bodies composed
of multiple penetrable regions. Although the results presented in this Chapter will be for
low order geometric discretizations, the extension of the separable expansion method to
higher order elements has been presented in Chapter 5. Such discretizations for penetrable
scatterers presents no added difficulty to applying the separable expansion. The Chapter is
organized as follows. We first will formulate the problem of scattering from a body composed
of penetrable and PEC regions. Next, we will look at its discretization using the method
of moments with the separable expansion to enable accurate integration. Lastly, we will
show scattering from various geometries. Stability comparisons will be made between the
separable expansion method and the quasi-exact integration method of [1]. Additionally, the
105

agreement of the two methods will be shown by comparing the early time signature of the
current against each other.

7.2

Formulation

Consider a body Ω composed of Q homogeneous regions, either penetrable ( r,q , µr,q ) or
PEC, illuminated by a field {Ei (r, t), Hi (r, t)} that is assumed band-limited and zero for
t ≤ 0. Let the region exterior to the scatterer (assumed to be free space) be denoted by
Ω0 . Within each region, Ωq , assume there exist an impressed source(s), giving rise to fields

Figure 7.1: Pictorial representation of scatterer with multiple homogeneous regions (graphic
from [1])

inc
{Einc
q , Hq }. The total fields in region in Ωq are then given by
inc
scat
Etot
q (r, t) = Eq (r, t) + Eq (r, t) ◦ Jq , Mq
inc
scat
Htot
q (r, t) = Hq (r, t) + Hq (r, t) ◦ Jq , Mq
106

(7.1)

scat
where {Escat
q , Hq } are the fields scattered by equivalent sources {Jq , Mq }, residing on
the boundary of Ωq , denoted by ∂Ωq . In region q, these scattered fields are given as
Escat
◦ Jq , Mq = −µq Le,q (r, t) ◦ Jq − Lh,q (r, t) ◦ Mq
q

(7.2)

Hscat
◦ Jq , Mq = − q Le,q (r, t) ◦ Mq + Lh,q (r, t) ◦ Jq
q
with the operators
.
Le,q ◦ {X} =
.
Lh,q ◦ {X} =

t
dt
0

X(r , τq )
dS ∂ 2 ¯I − c2q ∇∇ ·
t
4πR
∂Ωq

X(r , τq )
dS∇ ×
4πR
∂Ωq

(7.3)

where cq denotes the speed of light in medium q, q = r,q 0 , µq = µr,q µ0 , and τq = t−R/cq
is the retarded time.

To begin developing the PMCHWT formulation, the problem is recast as an equivalent
interior and exterior problem for each region. To simplify the discussion we focus on one
penetrable, homogeneous region with boundary ∂Ω. We denote the interior and exterior re−
+ +
gions as Ω− and Ω+ , with parameters ( −
r , µr ) and ( r , µr ), respectively. For the exterior
problem, all fields in Ω− are assumed to be zero, and all of R3 is assumed homogeneous
+
with parameters ( +
r , µr ). For the interior problem, all fields in Ω+ are zero, with all of
−
R3 homogeneous with parameters ( −
r , µr ). The equivalent currents for each problem are
given by
J± (r, t) = ±ˆ
n × Htot,±

(7.4a)

M± = Etot,± × ±ˆ
n

(7.4b)

107

where {Etot,± , Htot,± } are the total fields for the equivalent interior and exterior problems
and n
ˆ is the normal to surface ∂Ω pointing outward from region Ω− . The boundary conditions on the tangential fields can then relate the interior and exterior currents. At the
interface, ∂Ω, between homogeneous regions, we have continuity in the tangential components of the fields
n
ˆ × Etot,+
n
ˆ × Htot,+

r∈∂Ω
r∈∂Ω

=n
ˆ × Etot,−
=n
ˆ × Htot,−

r∈∂Ω
r∈∂Ω

(7.5a)
(7.5b)

which lead to the conditions

.
J(r, t) = J+ (r, t) = −J− (r, t)

(7.6a)

.
M(r, t) = M+ (r, t) = −M− (r, t)

(7.6b)

As a side note, for PEC boundaries the interior fields are zero and the boundary conditions
for the exterior fields are
n
ˆ q × Etot,+
n
ˆ q × Htot,+

=0

(7.7a)

= J(r, t)

(7.7b)

r∈∂Ω

r∈∂Ω

This means that there are no equivalent magnetic currents, and only exterior electric currents
as defined in (7.7b). This is consistent with what was shown in section 2.1.
Lastly, the tangential components of both Etot,± and Htot,± are equated on the boundary ∂Ω to yield the PMCHWT formulation for penetrable scatterers [51]

n
ˆ × Einc,+ (r, t) − Einc,− (r, t) =
108

(7.8a)

−n
ˆ × Escat,+ (r, t) ◦ {J, M} − Escat,− (r, t) ◦ {J, M}

n
ˆ × Hinc,+ (r, t) − Hinc,− (r, t) =

(7.8b)

−n
ˆ × Hscat,+ (r, t) ◦ {J, M} − Hscat,− (r, t) ◦ {J, M}

where
±
Escat,± ◦ {J, M} = −µ± L±
e (r, t) ◦ {J} − Lh (r, t) ◦ {M}

(7.9)

±
Hscat,± ◦ {J, M} = − ± L±
e (r, t) ◦ {M} + Lh (r, t) ◦ {J}
are the scattered electric and magnetic fields for the exterior (+) and interior (-) problems.
±
The operators {L±
e , Lh } are the same as {Le,q , Lh,q }, but for the homogeneous regions of
the interior and exterior problems.

7.3

Discretization

To discretize the system, the method of moments is used in a spatio-temporal basis/testing scheme similar to 3.1. First, the equivalent magnetic and electric currents on
∂Ωq are approximated as
Nq



Nt



Jq (r, t) =

Sn,q (r) 
Jq,n,i Ti (t)
n=1
i=1


Nq
Nt
Mq (r, t) =
Sn,q (r) 
Mq,n,i Ti (t)
n=1
i=1
109

(7.10)

where the spatial basis functions, Sn,q (r), are RWG basis functions and the temporal basis
functions, Ti (t) = T (t − i∆t ), are first order shifted Lagrange polynomials (hat functions).
Next, Galerkin testing in space and point testing in time is used to yield a marching system
with the same form as (3.5a), where




0
...
0
 Z i−j,0



0
Z i−j,1




.
.

Z i−j = 


.
.




.
.


Z i−j,Q
0




















(7.11a)




ee
eh
 Z i−j,q Z i−j,q 
Z i−j,q = 

he
hh
Z i−j,q Z i−j,q
ee
Z i−j,q
eh
Z i−j,q

mn

(7.11b)

−
= µq δ (t − (i − j) ∆t ) Sm,q (r) , {L+
e,q + Le,q } ◦ Sn,q (r) T (t)

= δ (t − (i − j) ∆t ) Sm,q (r) , {L+ + L− } ◦ Sn,q (r) T (t)
h,q
h,q
mn

(7.11c)
(7.11d)

ee
Z i−j,q

mn

= − δ (t − (i − j) ∆t ) Sm,q (r) , {L+ + L− } ◦ Sn,q (r) T (t)
h,q
h,q

(7.11e)

ee
Z i−j,q

mn

−
= q δ (t − (i − j) ∆t ) Sm,q (r) , {L+
e,q + Le,q } ◦ Sn,q (r) T (t)

(7.11f)

I j = I 0,j I 1,j ... I Q,j

T

I q,j = Jq,1,j Jq,2,j ... Jq,Nq ,j Mq,1,j Mq,2,j ... Mq,Nq ,j
V i = V 0,i V 1,i ... V Q,i
110

T

(7.11g)
T

(7.11h)
(7.11i)

V q,i =

inc,+
inc,−
δ (t − (i − j) ∆t ) S1,q (r) , Eq
(r, t) − Eq
(r, t)
inc,+
inc,−
δ (t − (i − j) ∆t ) S2,q (r) , Eq
(r, t) − Eq
(r, t)

(7.11j)
...

inc,+
inc,−
δ (t − (i − j) ∆t ) SNq ,q (r) , Eq
(r, t) − Eq
(r, t)
inc,+
inc,−
δ (t − (i − j) ∆t ) S1,q (r) , Hq
(r, t) − Hq
(r, t)

...

inc,+
inc,−
δ (t − (i − j) ∆t ) S2,q (r) , Hq
(r, t) − Hq
(r, t)
inc,+
inc,−
δ (t − (i − j) ∆t ) SNq ,q (r) , Hq
(r, t) − Hq
(r, t)

7.4

T

Separable Expansion

As was the case for PEC structures, it is not prudent to directly compute the matrix entries (7.11c)-(7.11f) using numerical quadrature, as this will not be accurate and the resulting
solutions will be prone to instability. To facilitate the accurate evaluation of (7.11c)-(7.11f),
we use the separable expansion method introduced in Chapter 4. The modification to each
of the operators through the expansion (4.8) can essentially be found in Chapter 4, so we
will not explicitly present them here. The only important modifications come in substituting
the proper material parameters q and µq into the relevant expressions. As with the PEC
case, the support of the expansion depends on (1) the number of edges per wavelength and
(2) the time step size with respect to the excitation frequency. Assuming the same level of
discretization, the required number of harmonics should be no different than in the PEC
case.
111

7.5

Low Frequency Instability

As discussed briefly in the introduction, the scattered field operators are prone to low
frequency breakdown, which manifests as low frequency instability in the time domain. This
is due to the fact that, at DC, the electric current produces no electric field (equivalently
the magnetic current produces no magnetic field). There exist a number of remedies to
low frequency instability [13, 17, 52]. The method in this thesis is not meant to address low
frequency instability, but rather the high frequency instability resulting from inaccuracy in
filling the system matrix. We have not used any method to address low frequency breakdown
in the implementation for penetrable scatterers, so low frequency breakdown has been seen
in some results. To address this, one of the aforementioned techniques could be used to
eliminate the low frequency instability in tandem with the separable expansion method
to eliminate the high frequency instability. We note that when low frequency instability
has been observed for the separable expansion, it has also been observed for quasi-exact
integration when everything has been discretized to the same level of precision.

7.6

Numerical Results

In this section we present results to validate the separable expansion method as applied to
penetrable scatterers. The objective is to show the absence of any high frequency instability
in the results and to compare the stability to solutions obtained via quasi-exact integration
method of [1]. Additionally, we aim to show the accuracy of the method in comparison the
quasi-exact integration by visually comparing the currents. For each simulation, the incident
field is a plane wave, as defined in Chapter 5.
112

Figure 7.2: Current on sphere

The first result is scattering from a homogeneous dielectric ( r = 2.0) sphere of radius
1 m. The sphere is discretized using 576 spatial unknowns and the simulation is run for
4,800 time steps (424 transits across the scatterer). The incident field has f0 = 40 MHz and
B = 20 MHz, kˆ = zˆ, and uˆ = xˆ. The result is shown in figure 7.2 and can be seen the be
stable throughout the simulation. Also plotted is the solution using quasi-exact integration.
The two plots agree extremely well. Both appear to be absent any type of instability.
Our next result is scattering from a PEC sphere, which is enclosed within a dielectric
sphere of r = 2.0. The PEC sphere has radius r = 0.5 m and the dielectric sphere has a
radius of 1 m. The PEC and dielectric boundaries are each discretized with 576 unknowns.
The parameters of the incident field are f0 = 60 MHz, B = 30 MHz, kˆ = zˆ, and uˆ = xˆ.
Figure 7.3 shows the surface current for 4,500 time steps. Both solutions appear stable and
the agreement between the two is excellent.
The next scatterer is the same thin box from section 5.2.3, but instead of PEC, is a
113

Figure 7.3: Current on composite sphere compared against exact integration
homogeneous dielectric region ( r = 1.5). This geometry is illuminated by an incident wave
of f0 =90 MHz, B =80 MHz, kˆ = zˆ, and uˆ = xˆ. This simulation was performed for 2,890
time steps (185 transits) and the result is shown in figure 7.4. This is a challenging problem
and some evidence of instability is seen. It, however, is of the low frequency variety, which
neither the separable expansion or quasi-exact integration methods are designed to remedy.
The performance of the two methods are very similar for this problem and the two results
agree very well.
The last result is a dielectric cone-sphere ( r = 1.5) of 1,008 spatial unknowns, with
the cone aligned in the x direction. The radius of the sphere is 0.25 m, and the height of
the cone is 1 m. The incident wave is described by f0 =250 MHz, B =150 MHz, kˆ = yˆ,
and uˆ = −ˆ
z . The resulting current shown in figure 7.5. As with the thin box, this is a
more challenging problem than the spherical scatterers due to the sharp geometric features.
Again, a low frequency instability is observed, but there is no evidence of high frequency
114

Figure 7.4: Current on thin box compared against exact integration

instability. This is not to say none exists as a small high frequency unstable mode may be
present, but dominated by the low frequency mode. An eigenvalue analysis would confirm
the absence of a high frequency instability. Alternatively, recourse to a method such as
[13, 17, 52] in tandem with either the separable expansion or quasi-exact integration method
could be used.

7.7

Discussion

In this Chapter we have presented the extension of the separable expansion method
to scattering from penetrable bodies. Validation of the procedure has been performed by
examining scattering from a variety of targets and comparing the stability and accuracy of
the result with quasi-exact integration. The stability of the two methods for these problems
is very similar, with no evidence of high frequency instability. A low frequency instability
115

Figure 7.5: Current on cone-sphere compared against exact integration
has been observed for both in two of the results. This is not a concern as remedies exist
for this type of instability [13, 17, 52] and can be implemented in tandem with the separable
expansion method (and the quasi-exact integration method). Eigenvalue analysis of these
problems should be performed to confirm the absence of any high frequency instability.

116

Chapter 8
Generalized Method of Moments with Separable Expansion Method

8.1

Motivation

In this chapter we look to apply the separable expansion method to a time domain implementation of the Generalized Method of Moments (GMM) [53, 54] with locally smooth
surface parameterization. The chapter is organized as follows. First, we will discuss the
advantages of using the GMM methodology (as opposed to the conventional method of moments approach) to discretize the TD-CFIE for PEC surfaces. Second, we will motivate
the locally smooth surface parameterizations from [54] by discussing some of the restrictions
of standard geometric discretizations and how these smooth approximations can overcome
these. Third, we will motivate the novelty of this chapter, i.e., the application of the separable expansion method to time domain-GMM and local surface parameterizations. At
that point we will give some details on the overall algorithm, examining the local surface
approximations, GMM, and the separable expansion separately. Finally, we will present
some preliminary scattering results for the TD-MFIE and discuss future work to extend the
method to the TD-EFIE.
117

8.1.1

Generalized Method of Moments (GMM)

Conventional method of moments discretizations for electromagnetic surface integral
equations are built on RWG functions [31] or their higher order [41] alternatives. These
basis functions are carefully designed to be normally continuous across the edges between
cells, which implies their definitions are closely linked to the tessellation of the geometry.
Generally, it is not easy to mix the types and orders of basis functions used in a single
discretization of a problem. Often, foreknowledge of local physics in a problem can suggest
one basis be used one one area of the tessellation that is more optimal than those used elsewhere. For example, it may be desirable to use singular functions near edges and corners of
a geometry, while using lower order or plane-wave type functions to approximate the current
on smooth regions. Furthermore, the conventional methodology of basis function construction does not lend itself well to p-adaptivity (adaptivity in basis function order). These
were two of the main motivations behind the development of GMM [53, 54]. GMM moves
away from the conventional methodology by using partition of unity functions to enforce
the appropriate continuity conditions, which allows great freedom in the choice of functions
to locally approximate the surface current. Additionally, the GMM framework opens up
avenues towards p-adaptivity. This will be discussed in more detail later in this chapter.

8.1.2

Locally smooth surface approximations

When considering conventional method of moments solutions, it is also important to
look at the quality of geometric discretization. Typical surface meshes are often satisfactory,
but they can suffer from a number of problems in many scenarios. First, these discretizations are not well suited to local h- and g-adaptivity (mesh size and order, respectively).
118

Second, the mesh quality can suffer when the geometry is complex which can lead to poor
conditioning. Lastly, meshing software can often perform local mesh refinement, but this
process can result in non-conformal meshes. The two motivations of (1) h- and g-adaptivity
and (2) independence from poorly generated meshes spurred the work of [54]. In [54] a
surface parameterization technique was introduced that was locally smooth with simple expressions for derivatives defined on these surfaces. These locally smooth domains could also
be merged to form larger domains or subdivided into smaller domains, which would allow for
h- and g-adaptivity. Moreover, the procedure could begin from primitives as simple as point
clouds, which are much simpler to generate than surface meshes. In this chapter we will
briefly outline the process of locally fitting a smooth approximation to the scatterer surface.
The coupling of GMM and locally smooth surface approximations, as reported in [54] (for
acoustics), paves the way for hp-adaptivity.

8.1.3

Time Domain GMM with Separable Expansion

In the time domain, the same issues of instability are confronted as were discussed earlier.
Due to the temporal basis function, the spatial integrands in the TD-EFIE and TD-MFIE are
not smooth, and if they are not evaluated to high precision, then the solution will be prone
to high frequency instability. Within the framework of this chapter, the spatial domains of
integration are the smooth higher order parameterizations discussed in the previous section.
As highlighted numerous times in this thesis, the quasi-exact methods cannot be extended
to higher order tessellations for reasons not restated here. The separable expansion method
is, therefore, a suitable choice for evaluating these integrals. In this chapter we will briefly
discuss how the expansion (4.8) is incorporated into the time domain GMM scheme on the
119

locally smooth parameterizations. It will be seen that this is not difficult as the expansion
has no effect on either the basis functions or the geometric discretization, but only on the
Green’s function.

8.2

Methodology

To begin, we look at the problem of scattering from a PEC scatterer in free space (for
details and definitions see section 2.3). We will first detail how the scatterer surface is
discretized using locally smooth surface parameterizations. Then, we will discuss the GMM
framework for representing the surface currents. Lastly, we will present the overall discretized
system and discuss how to accurately evaluate matrix entries using the separable expansion.

8.2.1

Locally smooth surface approximation

We now discuss how the locally smooth surface parameterizations are used to approximate
the scattering surface. First, the following quantities are assumed to be known
1. A collection of nodes, rk ∈ ∂Ω, k = 1, 2, ..., Nnode .
2. Nearest neighbor mappings between these nodes.
.
ˆ (rk )
3. The outward surface normals at these points n
ˆk = n
This information is already contained in surface triangulations, but the algorithm can also
start from something as simple as a point cloud. Given the scatterer boundary ∂Ω, the surface
is first subdivided into a Np overlapping patches Ωk (see figure 8.1 such that

Np
Ω = ∂Ω.
k=1 k

Each patch Ωk is associated with a node rk and formed by the neighborhood of rk defined
by its nearest neighbors. To construct a GMM patch local to node rk , the nearest neighbors
120

Figure 8.1: Overlapping GMM Patches (Ωk )

r
n
ˆk

r

Figure 8.2: Patch normal and projection planes Γk

of rk are first projected into the projection plane of rk , Γk (figure 8.2). Γk is the plane
passing through rk with normal n
ˆ k,avg . The normal vector, n
ˆ k,avg , is the patch normal,
which is defined as the average of the normals associated with all points belonging to the
patch. Next, a local, orthogonal coordinate system (u, v, w) is defined on Γk , where the
coordinate w is normal to the plane Γk . To construct a smooth, local approximation to Ωk ,
a least squares polynomial fit of order pk is performed to yield the locally smooth surface
Λk . For definitions of surface Jacobians, normal vectors, and differential operators on an
approximation surface, Λk , the reader is referred to [54].
121

8.2.2

Definition of basis functions

We will now discuss how basis functions are defined with the GMM framework on these
approximation surfaces. As discussed earlier, GMM is a discretization framework based on
partition of unity functions. In other words, the basis functions used in GMM are composed
of (1) a partition of unity function and (2) an approximation function. The former is used
to enforce the appropriate order of continuity between patches, while the latter is used
to accurately represent the local physics of the problem. More explicitly, the current is
approximated as
Np Bk
ˆ (r, t) =
J

Nt
jk,m (r)

k=1 m=1

j=1

j
I
T (t − j∆t )
k,m

(8.1a)

where Bk is the number of basis functions associated with patch k. GMM does not affect
the form of the temporal basis functions T (t), so the same temporal functions as those used
in typical method of moments discretizations, such as shifted Lagrange polynomials, splines,
or prolates can be used. For the results in this chapter we use hat functions as the temporal
basis. The spatial basis functions are defined as

.
jk,m (r) = ψk (r)ν k,m (r)

(8.1b)

where ψk is the partition of unity function and ν k,m is the approximation function. As
discussed in section 8.1.1, ν k,m can be chosen locally to optimally represent the physics
on a given patch. Some examples include polynomial functions, singular basis functions,
radial functions, or plane wave type functions. In this work, we use functions inspired
by the Helmholtz decomposition. Given two scalar functions φk,m (r) and ψk,m (r), the
122

(a)

(b)

Figure 8.3: Partition of unity functions: (a) 1D, (b) 2D

approximation function is designed to satisfy

νk,m ∈ span ∇s φk,m (r), ∇s × [ˆ
n(r)ψk,m ]

(8.2)

The explicit form of νk,m is determined by following the procedure given in [53].
The partition of unity is defined to have 2 characteristics. First, the partition of unity
function associated with a given node rk on patch Λk , must be equal to 1 at rk and equal
to zero at the patch boundaries. Second, consider a point r ∈ Λk and let the set of partition
of unity functions {ψk , ψk , ..., ψk } be the only partition of unity functions for which
1
2
N
r ∈ supp ψk . Then, the following condition must hold
i
N
i=1

ψk (r) = 1
i

(8.3)

1D and 2D examples of this type of function are shown in figure 8.3. The PU function can
be defined by starting with function λk , which is equal to 1 at node rk and zero at the
j
j
123

edge of patch Ωk . The PU is then defined as
j
λk (r)
j
N λ (r)
i=1 ki

ψk =
j

(8.4)

Different orders of ψk can be defined by using different orders of functions for λk .
j
j

8.3

GMM System

To fully discretize the system, Galerkin testing is employed in space and point testing is
used in time. The resulting marching system has the same form as (3.5a), with the matrices
Z i−j written in block form as










Z i−j = 









11
Z i−j
21
Z i−j


1Np
Z i−j 







.




.




.

Np Np 
. . . Z i−j

12
Z i−j . . .
22
Z i−j

.
.
.
Np 1
Z i−j

(8.5a)

where the entries of each block are given by

kl
Z i−j

mn

= δ (t − (i − j) ∆t ) jk,m (r), Lc jl,n (r)T (t)
m = 1, 2, ..., Bk

n = 1, 2, ..., Bl
124

(8.5b)

j
j
j
j
j
I j = I1,1 I1,2 ... I1,B I2,1 ... IN ,B
p Np
1

T

T
T
T
Vi = VT
i,1 V i,2 ...V i,Np

V i,k

m

(8.5c)

(8.5d)

ˆ (r) × Hinc (r, t)
= δ (t − (i − j) ∆t ) jk,m (r), αc /η0 E inc (r, t) + ((1 − αc ) n
m = 1, 2, ..., Bk

(8.5e)

It should be noted that although equation (8.5a) has been shown as a full matrix, it is
kl
actually sparse, i.e., Z i = 0 for all but a few values of i due to the temporal locality of T (t).
Another important observation is that the matrices in (8.5a) are very low rank. Therefore, it
is necessary to use a low rank solver to solve the system. As we only analyze small scatterers
in this chapter, the results use an SVD of the entire system. For larger problems, it is
necessary to orthogonalize sub-blocks of the system.

8.4

Accurate Integration via Separable Expansion

Again, we employ the separable expansion method to accurately compute the entries
(8.5b) and (8.5b). As both GMM and the local approximations only affect the basis function
expressions and the domains of integration, respectively, and have no affect on the Green’s
function, the precise modifications to these entries through the introduction of expansion
(4.8) should be no different than in previous sections. However, we note one important
implementation detail. The upper limit of (4.8) required for convergence is closely related to
the size of the support of the Legendre polynomials, which is determined by the spatial size of
a given source element, i.e., as the size of supp {P l(k1 R/c + k2 )} goes down, so does M . In
125

general when integrating over the projection planes Γk , it is first necessary to triangulate Γk
and then integrate over each sub-triangle. Therefore, the size of the support of the Legendre
polynomials is determined by the size of these sub-triangles. Due to the nature of the spatial
integrands, it may be desirable to make these sub-domains quite small. These integrands are
described by the following: (1) higher order approximation functions, (2) partition of unity
functions, potentially of high order, (3) a higher order surface over which to integrate, and
(3) higher order Legendre polynomials of (4.8). This can lead to integrands that, though
smooth, require a very high order of integration to evaluate accurately, which can severely
slow down the matrix fill. More study needs to be performed on the number of harmonics
required to converge for a given domain size, order of surface approximation, and type of
spatial basis function.

8.5

Numerical Results

In this section we present preliminary results for the stabilization of the time domain
GMM discretization of the TD-MFIE using the separable expansion method. The true test
of the algorithm is in its stability in solving the TD-EFIE. For reasons to be discussed in
the next section, we have not yet obtained stable TD-EFIE results here and leave this for
future work. Each result is for scattering from a plane wave with the same parameters of
chapter 5. The discretizations all begin from standard triangular meshes and the smooth
surface approximations are to 2nd order.
Our first result is for scattering from a nacelle, which is a cavity similar in shape to inlets
commonly found on aircraft. The nacelle has a radius of 10 cm, depth of 30 cm, and is
discretized using 840 triangles. Figure 8.4 shows a current coefficient for 800 time steps, as
126

Figure 8.4: Current on nacelle
the result continues to decay below the axis limits for the entire simulation. Our next result
is scattering from a cone-sphere, discretized with 1104 triangular elements. The spherical
section of the geometry has a radius of 29 cm, whereas the height of the cone is 76 cm.
The incident wave has a center frequencyf0 = 20 MHz, bandwidth B = 10 MHz, uˆ = −ˆ
z
polarization, and is incident in the kˆ = xˆ direction. Figure 8.5 shows a current coefficient
for 800 time steps and the result is seen to decay for the entire simulation.

8.6

TD-EFIE instability

Our simulations of PEC scattering using the TD-EFIE, discretized with GMM have not
yielded stable results, even when the space-time separation has been employed. Through
numerical experiments, we have concluded that the problem lies with the discretization of
the scalar potential on curved surfaces. This conclusion relies on a number of observations.
First, stable results have been obtained using the MFIE and by artificially setting the scalar
127

Figure 8.5: Current on cone-sphere
potential contributions to zero. In both the Lh and vector potential operators, no surface
gradients are required. The particular problem with the scalar potential could be due to
a number of things. First, the gradient of the scalar potential on the smooth surface parameterization may have been evaluated improperly. We have seen stable results for the
TD-EFIE on flat plates, which points to this possibility. Second, the integration of these
functions may not be accurate. Convergence studies should be performed to determine if
this is the case. Lastly, the recursive charge term that is moved to the right hand side of the
marching system could be causing a problem. The orthogonalization procedure should not
prohibit this exercise, but perhaps an implementation of the differentiated TD-EFIE should
be used to rule this out. Of course, in implementing this formulation, the space of temporal
basis and testing functions need to be changed.

128

Chapter 9
Conclusions and Future Work
In this thesis, we have presented a methodology to accurately compute matrix entries
involved in discretized TDIEs for electromagnetics. Accuracy in this procedure is vital if the
final solution is to be free of high frequency instabilities. The method is based on a separable expansion of convolutions with the retarded potential Green’s function. The expansion
results in smooth spatial integrands, which can be evaluated to arbitrary precision using
standard quadrature rules. This is in contrast to the quasi-exact integration methods, which
rely on analytical integration over shadow regions. By moving away from analytical integration and the identification of shadow regions, the method is trivially extended to higher
order tessellations, which greatly improve the accuracy in the geometric representation.
In this thesis we have contributed the following:
• Details for implementing the separable expansion within the TD-EFIE and TD-MFIE
for PEC scatterers.
• Studies of the error introduced in truncating the separable expansion as a function of
the truncation limit. An analytical expression for the bound on this error has been
presented and numerical experiments have been performed to show the relationship
between error and truncation limit.
• A novel space-time Galerkin basis/testing scheme which can achieve higher order accuracy in time. This is accomplished by restricting the support of the temporal functions
to one time step and using multiple functions within the same domain of support. Con129

vergence studies have been performed via interpolation tests and solutions to scattering
from simple objects. The implementation of the separable expansion method within
this framework has been detailed and results for scattering from various challenging
objects have been presented.

• Acceleration of the TD-EFIE and TD-CFIE using the plane wave time domain algorithm to compute far field interactions, in combination with the separable expansion
for near field interactions. Timing results have been shown to verify proper scaling and
scattering results from two challenging structures have been presented. Stability was
seen for the TD-CFIE results. Stabilization of the PWTD accelerated TD-EFIE with
the separable expansion is left for future work.

• Extension of the separable expansion method to the PMCHWT formulation for penetrable scatterers. Scattering results have been presented in comparison to the quasiexact integration scheme [1]. We have shown that the separable expansion method’s
stability properties are very similar for these problems and the solutions obtained by
the two methods agree very well. Although low frequency instability has been seen in
some of the results, the separable expansion method is not designed to address this
type of instability. Well developed methods exist for addressing this type of instability
[13, 17, 52].

• Implementation of the separable expansion within the GMM discretization methodology on surfaces approximated by locally smooth parameterizations. We have presented
preliminary scattering results for the TD-MFIE. Stable results have not yet been obtained for the TD-EFIE within this framework.
130

• A methodology for evaluating singular or near-singular integrals on higher order elements. Detailed prescriptions have been given for two classes of higher order surfaces:
mapped and projection surfaces. Comparisons between the presented methodology
and the sinh−1 rule [55] have been numerically investigated.

9.1

Future Work

The work presented in this thesis can be expanded into a number of different research
directions. We will now list (1) remaining areas in which the method can be strengthened
and (2) other settings in which the method has not been applied, but would be of great
utility.

9.1.1

PWTD acceleration of TD-EFIE with separable expansion

In this thesis we have presented results for both the TD-EFIE and the TD-CFIE. For
the TD-EFIE as applied to flat plates, we observed a growing low frequency term in most of
the solutions. One possible remedy is to accelerate the augmented TD-EFIE [16] to remove
the static solenoidal currents from the solution. This method does not affect the Green’s
function, so the extension of the separable expansion to this framework is trivial.
More troublesome is that we have observed high frequency instabilities in the accelerated
TD-EFIE for more complex targets. These same targets have been analyzed stably by
solving the system directly and using the separable expansion method. This suggests that
the instability is due to inaccuracies in the PWTD implementation. As shown in section
6.5, the approximate prolate spheroidal wave functions we are using are only approximately
time-limited, and can introduce erroneous high frequency content into the scattered fields if
131

they are multiplied by an insufficiently wide time window. We have observed that this is a
particular problem when interactions take place higher up the tree or when leaf box sizes are
made arbitrarily large. The two obvious paths to addressing this challenge are (1) to widen
the time window on the prolates or (2) devise a filtering scheme to remove the extraband
content. Extension to the TD-EFIE is attractive due to the superior accuracy properties of
the TD-EFIE and its applicability to open scatterers.

9.1.2

GMM discretization of TD-EFIE

The implementation of GMM presented in this thesis has been restricted to the TDMFIE due to stability considerations. Convergence studies of the scalar potential term need
to be performed as it appears to be the source of the TD-EFIE and TD-CFIE instability.
Extension of the method to the time-differentiated TD-EFIE would simplify the procedure
as this removes the time integral from the scalar potential. If this is done, the spaces of
temporal basis and testing functions must be carefully modified.

9.1.3

Nonlinear solvers

One of the main reasons for analyzing systems in the time domain is that systems with
nonlinearities can easily be analyzed. Solvers than can couple full wave electromagnetic
systems to circuit simulators have a wide range of applications. Work on extending the
separable expansion method to these problems is in the early stages.
132

9.1.4

Antenna problems

Time domain analysis of antenna problems is very useful. Often broadband data is
required, which time domain solutions are well suited for. The transient properties of the
system can be of interest as well. All of the problems presented in this work have been for
scattering problems, but the extension to radiation problems would be straightforward.

133

APPENDIX

134

Appendix A

A.1

Accurate computation of Singular and Near Singular Source
integrals

In this appendix we present a method for accurately evaluating source integrals when
R → 0, i.e. when a source and test triangle are close or coincide. This issue has received
much attention from the IE community [55–58]. In this appendix we present a singularity
cancellation method which is general for both mapped and projection surfaces. The latter
is of particular interest in chapter 8 of this thesis, when the separable expansion method is
applied to the Generalized Method of Moments with smooth surface approximations [53].
One commonly used method is that of [55]. This method relies on an arcsinh transformation
to cancel the 1/R singularity. One characteristic of this method is prohibitive to use on higher
order tessellations. When the projection of an observation point lies outside of the source
element, the domain of integration is extended outside of the element. The contributions
from these sub-domains do not corrupt the final answer, however, because their normals are
oriented in such a way as to perfectly cancel out. This is only the case when flat tessellations
are being used. It cannot be applied to curvilinear discretizations as the mappings from
the parent triangle are invalid for external points. All graphics and results presented in this
appendix are from [59].
Given a surface Ω and an observation point, r, the problem of interest is to evaluate the
135

integral
I=

dr
Ω

fn (r )
|r − r |

(A.1)

where fn (r ) is some known function. In the method of moments discretizations presented in
chapter 3, fn (r ) is either one component of a spatial basis function, or the surface divergence
of the spatial basis function. The surface Ω can be written as a transformation operator T
acting on some 2-dimensional plane, Γ, described by variables u1 and u2 . Assuming T to
be a polynomial transformation, two such classes should be considered separately.

A.1.1

Mapped Surfaces

The first type of surface is the mapped surface. Given a set of control points rn , a point
r ∈ Ω is given by
r (u) =
n

g
Pn (u1 , u2 )rn

(A.2)

g
where Pn (u1 , u2 ) are polynomial functions and g is the order of the mapping. The domain
of u1 and u2 is the unit right triangle, i.e. 0 ≤ u1 + u2 ≤ 1. One popular example of a
mapped transformation is that of the GWP type curvilinear elements as defined in [41].

A.1.2

Projection Surfaces

The second class of surface is the projection surface. Given a local, orthogonal coordinate
system, described by vectors {ˆ
u1 , uˆ2 , uˆ3 }, a point r ∈ Ω is given by

r (u) = u1 uˆ1 + u2 uˆ2 + P g (u1 , u2 )ˆ
u3

(A.3)

where P g (u1 , u2 ) is a polynomial in u1 and u2 , which is complete to order g. Here, the coor136

ˆ
n
ˆ
u
ˆ
v

Figure A.1: Projection onto a curved surface from a projection plane

dinates u1 and u2 are typically, though not necessarily, local to one facet of the tessellation.
An example in which these coordinates are confined to a triangle is shown in figure A.1. As
(A.3) shows, r is the projection along the uˆ3 direction from a plane Γ, which is spanned by
vectors uˆ1 and uˆ2 . Implementations of these higher order projection surfaces are found in
[53, 54, 60]. Before describing how to evaluate (A.1) on a higher order surface of each type,
we present the following definitions, which hold for both mapped and projection surfaces.
1. Given u ∈ Γ corresponding to a point r ∈ Ω according to the specific transformation
T being used, the tangent vectors at this point are given by

. dr (u)
,
i (u) = du
i−1

(A.4)

where we are taking i − 1 modulo 3.

2. At u, the surface transformation Jacobian is computed as

g(u) =

∂r (u) ∂r (u)
×
∂ui
∂ui+1
137

(A.5)

3. Lastly, the unit normal to the surface at this point is
∂r (u) ∂r (u)
×
∂ui
∂ui+1
.
ˆ (u) =
n
∂r (u) ∂r (u)
×
∂ui
∂ui+1

(A.6)

4. Let u0 = (0, 0, 0) denote the origin of the coordinate system (u1 , u2 , u3 ). Note that for
projection surfaces, Γ is tangent to the curved surface Ω at r (u0 ) (assuming P g (0, 0) =
0).

A.1.3

Integration on Mapped Surfaces

Given the preceding definitions, we will now outline the procedure for integration on a
mapped surface.
1. Find the tangent triangle: Given an observation point, r, the first step is to find the
plane that is tangent to the curved surface Ω at the projection of r onto Ω. Denoting
the projection as r0 , this is equivalent to solving the problem

[r − r0 (u0 )] · i±1 (u0 ) = 0, i ∈ {1, 2}.

(A.7)

which is a nonlinear equation for u0 . A starting guess for the nonlinear solution can
be taken to be the vertices of the original, curved triangle. The curved surface and
tangent triangle are pictured in figure A.2.

2. Transform the integral onto tangent triangle: Once the vertices of the planar, tangent
138

2
•r
l1
• r0
1

l3

3

l2

Figure A.2: Tangent triangle for curvilinear, mapped triangles.
triangle have been determined, the integral (A.1) can be rewritten as

dr
Ω

fn (r )
fn (r (u, v))
=
.
du g(u, v)
|r − r |
|r − r (u, v)|
Γ

(A.8)

3. Numerically compute integral using appropriate rule: The procedure for developing
singularity cancellation rules on this surface will be covered in Section A.

A.1.4

Integration on Projection Surfaces

Before outlining the procedure for integration on projection surfaces, we note that these
higher order surfaces are typically not formed from a set of triangles. Therefore, the preliminary step of generating a triangulation must be performed to yield the projection surface.
Assuming this has been done, integration on a particular element of such a surface is outlined
as follows.
1. Find the tangent triangle: Again, given an observation point, r, the first step is to
tangent triangle. However, we have observed that the effect on accuracy between using
the tangent triangle and using the original projection surface is minimal. Therefore, the
solution to the nonlinear system (A.7) is deemed unnecessary for projection surfaces.
139

•r
•r

•u

Ω0

∆0

• r0
Figure A.3: Tangent triangle for curvilinear, projection surfaces.
The projection onto the projection plane is given by

ˆ (ˆ
r0 = r − n
n · [r − r 0 ]).

(A.9)

2. Transform the integral onto tangent triangle: Again, the integral is transformed as in
(A.8). We note that both the domain of integration and g(u, v) are not the same as
for mapped surfaces.

3. Numerically compute integral using appropriate rule: Again, singularity cancellation
rules will be covered in Section A.

A.1.5

Singularity cancellation rules for curvilinear surfaces

As discussed in Section A, the implementation of singularity cancellation rules on curvilinear surfaces can be challenging, particularly when the observation point and its projection
do not lie on the source triangle. For points lying in the source domain, the procedure for
evaluating such integrals is covered in depth in [55, 61]. We now outline the procedure for
the more difficult case of an observation point not belonging to the source triangle.
140

• r0
• r0
• r0
(a) case 1: r0 on an edge (b) case 2: r0 on an
altitude

(c) case 3

Figure A.4: Subdivision of source domain for observation points projecting outside

The case for r0 lying entirely outside of the source triangle is shown in figure A.4. Typically, for singularity cancellation rules, the source domain is subdivided into three triangles.
Each of these triangles is formed by the 2 vertices from the source triangle and r0 as the
third vertex. When r0 lies outside of the source triangle, the entirety of one of these subtriangles lies outside of the source domain. This is not an issue for flat surface descriptions,
but higher order transformations do not properly map these types of points back onto the
original surface. Therefore, it is necessary to restrict the domain of integration entirely on
the source triangle. The radial-angular rule [61] will be used to achieve this. In the radial
angular rule, the integration on a sub-triangle is transformed as

du
Γ
=

φ2
φ1

fn (r (u , v ))
|r − r (u , v )|
R2 (φ)
fn (r )
,
dR J (R, φ)g(u , v )
|r − r |
R1 (φ)

g(u, v)
dφ

(A.10)

We now outline how these sub-triangles are formed for r0 lying outside the source triangle.

Figure A.4 presents the three possible cases for r0 outside the source triangle. These are
when the projection lies along the extension of an edge (A.4a), an altitude (A.4b), or neither
141

•r
φ2 − φ1

• r0
ρ1 (φ)
ρ2 (φ)

Figure A.5: Integration limits for one sub-triangle

(A.4c). For the first case, the triangle is not subdivided, whereas for the other 2 cases,
the triangle is subdivided once into two sub-triangles. For the second case, the sub-triangle
vertices are (1) r0 , (2) one of the two vertices furthest from r0 , and (3) the normal projection
of r0 onto the opposite edge of the triangle. For the third case, the vertices are (1) r0 , (2)
one of the two nearest vertices to r0 , and (3) the vertex furthest from r0 .
The next step is to identify the limits in (A.10) (see figure A.5). To begin, the limits in
φ are found for each sub-triangle. A standard Gaussian quadrature rule is used to discretize
the integral in φ. Next, for each value of φ, the limits of R(φ) are determined and a Gaussian
quadrature rule is used.

A.1.6

Numerical Results for Radial-Angular Rule

In this Section we investigate the accuracy and convergence of the radial angular rule
applied to curvilinear surfaces. First, we validate the performance of this implementation of
the radial-angular and compare to the sinh−1 rule. To this end, we look at the number of
quadrature points required to converge to a given error tolerance for integrating 1/|r − r | on
a flat triangle. The error is computed in comparison to an adaptive integration rule. Three
different observation point locations are considered, i.e. when r0 lies (1) inside the triangle,
142

Analytical Tol.
1.865212 10−4
10−8
Analytical Tol.
1.394142 10−4
10−8
Analytical Tol.
1.099330 10−4
10−8

r = (1/3, 1/3, 0.1)
Adaptive sinh−1
304
108
3328
432
r = (0, 1/2, 0.1)
Adaptive sinh−1
108
147
3340
675
r = (0, 0, 0.1)
Adaptive sinh−1
102
147
3330
588

Radial-Angular
15
30
Radial-Angular
8
16
Radial-Angular
3
7

Table A.1: Comparison of efficiency of different rules
(2) on an edge of the triangle, and (3) on a vertex of the triangle. The results are shown in
table A.1.
It can be seen that the radial-angular outperforms the sinh−1 rule in each case. This is
expected as, for flat triangles, there is no φ dependence and the radial-angular rule transforms
the integrand into a natural coordinate system. The performance is more competitive for
curvilinear elements.
The next result quantifies the error for the radial-angular rule applied to both projection
and mapped surfaces for different observation points r. 3 triangles are considered, which are
a 2nd order GWP triangle as well as 2nd and 4th order projection triangles. Three different
cases of observation points are also examined, which are r0 in the interior of the triangle
(case 1), r0 on an edge of the triangle (case 2), and r0 on a vertex of the triangle (case 3).
The distance above the triangle |r − r0 | is also varied. The number of quadrature points
required to converge to a given tolerance for each case is shown in figure A.2.
The performance of the algorithm for mapped and projection surfaces are very similar.
For the 4th order triangle, there were 3 cases where the result failed to converge to the given
143

Error
10−4
10−8
10−12

Mapped triangle (g = 2)
Case 1 Case 2 Case 3
15
20
12
24
24
12
228
210
168

10−4
10−8
10−12

24
60
432

26
48
432

16
72
276

10−4
10−8
10−12

24
42
294

26
20
242

28
66
286

d = 0.1
Projection triangle (g = 2)
Case 1 Case 2 Case 3
15
20
8
47
24
12
230
240
144
d = 0.01
24
16
16
72
42
96
462
396
264
d = 0.001
24
12
28
42
20
60
282
224
312

Projection triangle (g = 4)
Case 1 Case 2 Case 3
24
44
84
42
68
120
560
1160
1440
24
42
462

60
120
1332

84
148
-

48
60
480

60
108
1872

120
-

Table A.2: Number of points required for convergence of integration rules on higher order
triangles for various test point locations
tolerance, but otherwise the algorithm converged within a reasonable amount of points.

144

BIBLIOGRAPHY

145

BIBLIOGRAPHY
[1] B. Shanker, M. Lu, J. Yuan, and E. Michielssen, “Time domain integral equation analysis of scattering from composite bodies via exact evaluation of radiation fields,” Antennas and Propagation, IEEE Transactions on, vol. 57, no. 5, pp. 1506–1520, 2009.
[2] Y. Shi, M.-Y. Xia, R. shan Chen, E. Michielssen, and M. Lu, “Stable electric field tdie
solvers via quasi-exact evaluation of mot matrix elements,” Antennas and Propagation,
IEEE Transactions on, vol. 59, no. 2, pp. 574–585, 2011.
[3] M. Friedman and R. Shaw, “Diffraction of pulses by cylindrical obstacles of arbitrary
cross section,” Journal of Applied Mechanics, vol. 29, no. 1, pp. 40–46, 1962.
[4] J. Jin, The finite element method in electromagnetics. John Wiley & Sons, 2014.
[5] A. A. Ergin, B. Shanker, and E. Michielssen, “The plane-wave time-domain algorithm
for the fast analysis of transient wave phenomena,” Antennas and Propagation Magazine, IEEE, vol. 41, no. 4, pp. 39–52, 1999.
[6] A. E. Yilmaz, J.-M. Jin, and E. Michielssen, “Time domain adaptive integral method for
surface integral equations,” Antennas and Propagation, IEEE Transactions on, vol. 52,
no. 10, pp. 2692–2708, 2004.
[7] J. Song, C.-C. Lu, and W. C. Chew, “Multilevel fast multipole algorithm for electromagnetic scattering by large complex objects,” Antennas and Propagation, IEEE
Transactions on, vol. 45, no. 10, pp. 1488–1493, 1997.
[8] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, “Aim: Adaptive integral method for
solving large-scale electromagnetic scattering and radiation problems,” Radio Science,
vol. 31, no. 5, pp. 1225–1251, 1996.
[9] D. S. Weile, G. Pisharody, N.-W. Chen, B. Shanker, and E. Michielssen, “A novel scheme
for the solution of the time-domain integral equations of electromagnetics,” Antennas
and Propagation, IEEE Transactions on, vol. 52, no. 1, pp. 283–295, 2004.
[10] M. Xia, G. Zhang, G. Dai, and C. Chan, “Stable solution of time domain integral
equation methods using quadratic b-spline temporal basis functions,” JOURNAL OF
COMPUTATIONAL MATHEMATICS-INTERNATIONAL EDITION-, vol. 25, no. 3,
p. 374, 2007.
146

[11] X. Wang, R. A. Wildman, D. S. Weile, and P. Monk, “A finite difference delay modeling
approach to the discretization of the time domain integral equations of electromagnetics,” Antennas and Propagation, IEEE Transactions on, vol. 56, no. 8, pp. 2442–2452,
2008.
[12] M. Costabel, “Time-dependent problems with the boundary integral equation method,”
in Encyclopedia of Computational Mechanics, John Wiley and Sons, 2004.
[13] F. P. Andriulli, K. Cools, F. Olyslager, and E. Michielssen, “Time domain calder´on
identities and their application to the integral equation analysis of scattering by pec
objects part ii: Stability,” IEEE transactions on antennas and propagation, vol. 57,
no. 8, pp. 2365–2375, 2009.
[14] D. S. Jones, “Methods in electromagnetic wave propagation,” Chemical Physics, vol. 1,
1979.
[15] B. Shanker, A. A. Ergin, K. Aygun, and E. Michielssen, “Analysis of transient electromagnetic scattering from closed surfaces using a combined field integral equation,”
Antennas and Propagation, IEEE Transactions on, vol. 48, no. 7, pp. 1064–1074, 2000.
[16] G. Pisharody and D. S. Weile, “Electromagnetic scattering from perfect electric conductors using an augmented time-domain integral-equation technique,” Microwave and
optical technology letters, vol. 45, no. 1, pp. 26–31, 2005.
[17] G. Pisharody and D. S. Weile, “Robust solution of time-domain integral equations using
loop-tree decomposition and bandlimited extrapolation,” Antennas and Propagation,
IEEE Transactions on, vol. 53, no. 6, pp. 2089–2098, 2005.
[18] A. Tijhuis, “Iterative determination of permittivity and conductivity profiles of a dielectric slab in the time domain,” Antennas and Propagation, IEEE Transactions on,
vol. 29, no. 2, pp. 239–245, Mar.
[19] A. Tijhuis, “Toward a stable marching-on-in-time method for two-dimensional transient
electromagnetic scattering problems,” Radio Science, vol. 19, no. 5, pp. 1311–1317, 1984.
[20] B. Rynne, “Stability and convergence of time marching methods in scattering problems,”
IMA journal of applied mathematics, vol. 35, no. 3, pp. 297–310, 1985.
[21] B. Rynne, “Time domain scattering from arbitrary surfaces using the electric field integral equation,” Journal of Electromagnetic Waves and Applications, vol. 5, no. 1,
pp. 93–112, 1991.
147

[22] D. A. Vechinski and S. M. Rao, “A stable procedure to calculate the transient scattering by conducting surfaces of arbitrary shape,” Antennas and Propagation, IEEE
Transactions on, vol. 40, no. 6, pp. 661–665, 1992.
[23] A. Sadigh and E. Arvas, “Treating the instabilities in marching-on-in-time method from
a different perspective [electromagnetic scattering],” Antennas and Propagation, IEEE
Transactions on, vol. 41, no. 12, pp. 1695–1702, 1993.
[24] E. van’t Wout, H. van der Ven, D. R. van der Heul, and C. Vuik, “A provably stable
mot scheme based on quadratic spline basis functions,” in Antennas and Propagation
Society International Symposium (APSURSI), 2012 IEEE, pp. 1–2, IEEE, 2012.
[25] T. Ha-Duong, “On retarded potential boundary integral equations and their discretisation,” in Topics in computational wave propagation, pp. 301–336, Springer, 2003.
[26] Y. Ding, A. Forestier, and T. H. Duong, “A galerkin scheme for the time domain integral
equation of acoustic scattering from a hard surface,” The Journal of the Acoustical
Society of America, vol. 86, p. 1566, 1989.
[27] T. Ha-Duong, B. Ludwig, and I. Terrasse, “A galerkin bem for transient acoustic scattering by an absorbing obstacle,” International journal for numerical methods in engineering, vol. 57, no. 13, pp. 1845–1882, 2003.
[28] J. A. Stratton, Electromagnetic Theory. McGraw-Hill College, 1941.
[29] R. Mittra, Computer techniques for electromagnetics, vol. 7. Pergamon, 1973.
[30] A. F. Peterson, “Mapped vector basis functions for electromagnetic integral equations,”
Synthesis Lectures on Computational Electromagnetics, vol. 1, no. 1, pp. 1–124, 2005.
[31] S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering by surfaces of arbitrary
shape,” Antennas and Propagation, IEEE Transactions on, vol. 30, no. 3, pp. 409–418,
1982.
[32] P. J. Davies, “Numerical stability and convergence of approximations of retarded potential integral equations,” SIAM Journal on Numerical Analysis, vol. 31, pp. 856–20,
06 1994. Copyright - Copyright] 1994 Society for Industrial and Applied Mathematics;
Last updated - 2012-02-23.
148

[33] P. J. Davies, “A stability analysis of a time marching scheme for the general surface
electric field integral equation,” Applied Numerical Mathematics, vol. 27, no. 1, pp. 33–
57, 1998.
[34] Y. Beghein, K. Cools, H. Bagci, and D. De Zutter, “A space-time mixed galerkin
marching-on-in-time scheme for the time-domain combined field integral equation,” Antennas and Propagation, IEEE Transactions on, vol. 61, no. 3, pp. 1228–1238, 2013.
[35] I. Terrasse, R´esolution math´ematique et num’erique des ´equations de Maxwell instationnaires par une m´ethode de potentiels retard´es. PhD thesis, Ecole Polytechnique,
1993.
[36] A. Yucel and A. Ergin, “Exact evaluation of retarded-time potential integrals for the
rwg bases,” Antennas and Propagation, IEEE Transactions on, vol. 54, pp. 1496–1502,
May 2006.
[37] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, eds., NIST Handbook
of Mathematical Functions. New York, NY: Cambridge University Press, 2010. Print
companion to [62].
[38] S. P. Walker, M. J. Bluck, and I. Chatzis, “The stability of integral equation timedomain computations for three-dimensional scattering ; similarities and differences between electrodynamic and elastodynamic computations,” International Journal for Numerical Modelling: Electronic Networks, Device and Fields, vol. 474, no. 15, pp. 459–474,
2002.
[39] A. Pray, N. Nair, and B. Shanker, “Stability properties of the time domain electric
field integral equation using a separable approximation for the convolution with the
retarded potential,” Antennas and Propagation, IEEE Transactions on, vol. 60, no. 8,
pp. 3772–3781, 2012.
[40] G. H. Golub and C. F. van Van Loan, Matrix Computations (Johns Hopkins Studies
in Mathematical Sciences)(3rd Edition). The Johns Hopkins University Press, 3rd ed.,
Oct. 1996.
[41] R. D. Graglia, D. R. Wilton, and A. F. Peterson, “Higher order interpolatory vector
bases for computational electromagnetics,” Antennas and Propagation, IEEE Transactions on, vol. 45, no. 3, pp. 329–342, 1997.
[42] M. Bebendorf, “Approximation of boundary element matrices,” Numerische Mathematik, vol. 86, no. 4, pp. 565–589, 2000.
149

[43] R. Coifman, V. Rokhlin, and S. Wandzura, “The fast multipole method for the wave
equation: A pedestrian prescription,” Antennas and Propagation Magazine, IEEE,
vol. 35, no. 3, pp. 7–12, 1993.
[44] B. Shanker and H. Huang, “Accelerated cartesian expansions a fast method for computing of potentials of the form r for all real ,” Journal of Computational Physics, vol. 226,
no. 1, pp. 732 – 753, 2007.
[45] K. Zhao, M. N. Vouvakis, and J.-F. Lee, “The adaptive cross approximation algorithm
for accelerated method of moments computations of emc problems,” Electromagnetic
Compatibility, IEEE Transactions on, vol. 47, no. 4, pp. 763–773, 2005.
[46] Y. Shi and J.-M. Jin, “Openmp parallelized mod solution of the time-domain efie accelerated by the aca algorithm,” Microwave and Optical Technology Letters, vol. 54, no. 5,
pp. 1206–1212, 2012.
[47] M. Vikram and B. Shanker, “Fast evaluation of time domain fields in sub-wavelength
source/observer distributions using accelerated cartesian expansions (ace),” Journal of
Computational Physics, vol. 227, no. 2, pp. 1007–1023, 2007.
[48] B. Shanker, A. A. Ergin, M. Lu, and E. Michielssen, “Fast analysis of transient electromagnetic scattering phenomena using the multilevel plane wave time domain algorithm,” Antennas and Propagation, IEEE Transactions on, vol. 51, no. 3, pp. 628–641,
2003.
[49] A. A. Ergin, B. Shanker, and E. Michielssen, “Fast evaluation of three-dimensional
transient wave fields using diagonal translation operators,” Journal of Computational
Physics, vol. 146, no. 1, pp. 157–180, 1998.
[50] R. Jakob-Chien and B. K. Alpert, “A fast spherical filter with uniform resolution,”
Journal of Computational Physics, vol. 136, no. 2, pp. 580–584, 1997.
[51] L. N. Medgyesi-Mitchang, J. M. Putnam, and M. B. Gedera, “Generalized method
of moments for three-dimensional penetrable scatteres,” J. Opt. Soc. Am. A, vol. 11,
pp. 1383–1398, 1994.
[52] G. Pisharody and D. S. Weile, “Electromagnetic scattering from homogeneous dielectric bodies using time-domain integral equations,” Antennas and Propagation, IEEE
Transactions on, vol. 54, no. 2, pp. 687–697, 2006.
150

[53] N. V. Nair and B. Shanker, “Generalized method of moments: a framework for analyzing
scattering from homogeneous dielectric bodies,” JOSA A, vol. 28, no. 3, pp. 328–340,
2011.
[54] N. Nair, B. Shanker, and L. Kempel, “Generalized method of moments: A boundary
integral framework for adaptive analysis of acoustic scattering,” The Journal of the
Acoustical Society of America, vol. 132, no. 3, pp. 1261–1270, 2012.
[55] M. A. Khayat and D. R. Wilton, “Numerical evaluation of singular and near-singular
potential integrals,” Antennas and Propagation, IEEE Transactions on, vol. 53, no. 10,
pp. 3180–3190, 2005.
[56] M. G. Duffy, “Quadrature over a pyramid or cube of integrands with a singularity at a
vertex,” SIAM J. Numer. Anal., vol. 19, no. 6, pp. 1260–1262, 1982.
[57] W. Brown and D. Wilton, “Singular basis functions and curvilinear triangles in the solution of the electric field integral equation,” Antennas and Propagation, IEEE Transactions on, vol. 47, pp. 347–353, Feb 1999.
[58] P. W. Fink, D. R. Wilton, and M. a. Khayat, “Simple and Efficient Numerical Evaluation
of Near-Hypersingular Integrals,” IEEE Antennas and Wireless Propagation Letters,
vol. 7, pp. 469–472, 2008.
[59] N. Nair, A. Pray, J. Villa-Giron, B. Shanker, and D. Wilton, “A singularity cancellation
technique for weakly singular integrals on higher order surface descriptions,” Antennas
and Propagation, IEEE Transactions on, vol. 61, no. 4, pp. 2347–2352, 2013.
[60] N. Nair, B. Shanker, and L. Kempel, “A discretization framework for scalar integral
equations using the generalized method of moments and locally smooth surface approximations,” in Antennas and Propagation (APSURSI), 2011 IEEE International
Symposium on, pp. 3193–3196, IEEE, 2011.
[61] M. A. Khayat, D. R. Wilton, and P. W. Fink, “An improved transformation and optimized sampling scheme for the numerical evaluation of singular and near-singular potentials,” in Antennas and Propagation Society International Symposium, 2007 IEEE,
pp. 4845–4848, IEEE, 2007.
[62] “NIST Digital Library of Mathematical Functions.” http://dlmf.nist.gov/, Release 1.0.5
of 2012-10-01. Online companion to [37].

151