IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 47, NO.

1, JANUARY 1999 9

A Finite-Difference Time-Domain Method

for Solving Electromagnetic Problems
with Bandpass-Limited Sources
J. D. Pursel, Member, IEEE, and P. M. Goggans, Member, IEEE

Abstract— The complex-envelope representation of bandpass- where and are the in-phase and quadrature portions
limited signals is used to formulate a bandpass-limited vector of and is the center frequency of [1]. The
wave equation and a new finite-difference time-domain (FDTD) bandwidth and center frequency of can be expressed
scheme that solves the bandpass-limited vector wave equation is
presented. For narrow-band electromagnetic systems, this new in terms of , the maximum frequency in , and ,
method allows the time step to be several orders of magnitude the minimum frequency in . For , the
larger than current FDTD formulations while maintaining an relationships are
amplification factor equal to one. Example results obtained by
this method are presented and compared with analytic solutions.
(2)
Index Terms— Bandpass-limited, complex-envelope, FDTD,
wave equation.
and
(3)
I. INTRODUCTION
The time functions and are real and low-pass
T HE usual finite-difference time-domain (FDTD) solution
of the electromagnetic wave equation uses time and spa-
tial sampling according to the lowpass-limited (LPL) sampling
limited with bandwidth /2.
Haykin [1] introduces the complex envelope notation to
make the bandpass-limited signal and system problem isomor-
theorem. For time sampling this means that the time step
phic to the base-band (LPL) signal and system problem. The
must be set in accordance with maximum frequency in the
complex envelope of is denoted as and is defined by
source signal. However, for radar and communication systems
problems the source signal is usually bandpass-limited (BPL)
rather than lowpass-limited. The bandpass-limited sampling (4)
theorem states that using the proper techniques, the signal can
be sampled in accordance with the bandwidth of the signal where the function can be recovered from using the
source rather than its maximum frequency. This suggests that expression
the FDTD method can be modified so that the required time
step for BPL signals is significantly increased compared to (5)
the conventional FDTD method.
One way of modifying the FDTD method is to use the The Fourier transform of the complex envelope can be deter-
complex envelope representation of BPL signals [1], [2]. The mined from the Fourier transform of using
complex envelope representation of a real BPL signal yields a
complex LPL signal with a maximum frequency equal to one sgn (6)
half the bandwidth of the original signal.
where

II. COMPLEX-ENVELOPE REPRESENTATION for

sgn for (7)
A real bandpass-limited signal source can be expressed for .
in the form [1], [3]
Here, the standard engineering definition of the Fourier trans-
(1) form is used so that

Manuscript received August 25, 1997; revised June 8, 1998.

J. D. Pursel was with the Department of Electrical Engineering, University (8)
of Mississippi, University, MS 38677 USA. He is now with the Riverside
Research Institute, New York, NY 10036 USA. and
P. M. Goggans is with the Department of Electrical Engineering, University
of Mississippi, University, MS 38677 USA. (9)
Publisher Item Identifier S 0018-926X(99)02207-3.

0018–926X/99$10.00  1999 IEEE

Authorized licensed use limited to: Maharaja Agrasen Institute of Technology. Downloaded on April 19,2010 at 14:05:55 UTC from IEEE Xplore. Restrictions apply.
10 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 47, NO. 1, JANUARY 1999

space dependent electric field as

(10)

where the LPL complex-envelope representation of is

(11)

In the expressions above, the position vector identifies the

field point. The real vector functions and are
the in-phase and quadrature parts of . Vector quantities
are printed in boldface. If the system were also spatially
bandpass-limited, a similar transformation could be used to
remove the high spatial-frequency variation. However, the
spatial-frequency content is dependent on the geometry of the
electromagnetic system. Hence, electromagnetic systems are
Fig. 1. The frequency-domain spectrum of the bandpass-limited signal
V (f ). not generally spatially BPL.
The standard real-valued low-pass-limited (base-band) wave
equation for linear isotropic time-invariant medium is as
follows:

(12)

where , , and are the spatially dependent permeability,

permittivity, and conductivity of the medium and where ,
, and are the impressed electric current, the impressed
magnetic source current and the total electric charge. In (12),
a dot above a time function denotes the first partial derivative
Fig. 2. The frequency-domain spectrum of the lowpass-limited signal V~ (f ).
with respect to time of the function. Similarly, a double dot
above a time function denotes the second partial derivative
It is instructional to consider the relationship between with respect to time of the function.
and . Fig. 1 illustrates the Fourier transform of Using the complex envelope representation of time func-
a real bandpass-limited time function . Because is tions given in (10) and (11), a new complex-valued bandpass-
real, where the superscript denotes the limited form of the vector wave equation
complex conjugate. As a result of this relationship, can
be determined from a knowledge of the positive frequency
values of alone. The transformation makes use of this
property by setting the negative frequency portion of to
zero, multiplying by two and then shifting the result to the left
by to yield . If the transformation of to
is viewed in the graphical way described above, it is clear
that can be obtained from by reversing the steps
above. Fig. 2 illustrates the Fourier transform of . It
is apparent from Fig. 2 that is low-pass limited. Because
in general , is in general a complex time (13)
function.
can be derived from (12).
III. COMPLEX-ENVELOPE WAVE EQUATIONS In (13), all of the low-pass-limited time functions in (12)
In electromagnetics, if the sources of the electric and are replaced with their complex-envelope representations. We
magnetic fields are bandpass-limited, all of the fields resulting refer to (13) as the complex-envelope vector wave equation.
from the sources are also bandpass-limited (for linear time- Rewriting (13) for the case where the impressed magnetic
invariant medium). As a result we can write the time and current and the conductivity are zero and the permeability is

Authorized licensed use limited to: Maharaja Agrasen Institute of Technology. Downloaded on April 19,2010 at 14:05:55 UTC from IEEE Xplore. Restrictions apply.
PURSEL AND GOGGANS: FDTD METHOD FOR SOLVING EM PROBLEMS WITH BANDPASS-LIMITED SOURCES 11

constant yields the following equation: Crank–Nicolson scheme:

(14)

IV. FINITE-DIFFERENCE TIME-DOMAIN FORMULATION

In this paper, we consider a FDTD solution of the complex- (19)
envelope representation of the wave equation [4]. The initial
work to develop a numerical method for the solution of Crank–Nicolson schemes are often used in the numerical so-
the complex-envelope wave equation was done using one- lution of parabolic partial differential equations. A discussion
dimensional (1-D) geometries [5]. A 1-D version of the on the use of this type of differencing scheme in the numerical
standard wave equation (where , , and is solution of the heat equation can be found in [7] and [9].
constant) can be obtained from (12) by assuming a solution for For initial-value problems whose finite-difference schemes
the electric field that has the single Cartesian vector direction have constant coefficients, the von Neumann stability con-
and is a function of a single spatial coordinate . This results dition [7], [10], [11] provides a necessary test of numerical
in the following differential equation: stability. If the scheme is also Hermitian as it is here, the
von Neumann stability condition is also a sufficient test of
numerical stability [7]. To test (19) for numerical stability,
(15)
consider a homogeneous source-free region as the domain of
interest and assume an initial electric field distribution
Converting (15) to the complex-envelope representation yields
(20)
where is a real-valued wave number. Now substitute an
assumed solution

(16) (21)
into (19). Doing so yields
It is interesting to note that (16) is very similar to the telegraph
equation given in [6].
There are many differencing schemes that can be used to
approximate (15) and (16). In addition to possessing some
specified level of accuracy, the scheme must also be numer-
ically stable. Some schemes will be unconditionally stable, (22)
some will be unstable, and some will be conditionally stable.
An explicit second order differencing scheme that approx- If (19) is a numerically stable FD scheme, the von Neumann
imates (15) is [7] condition requires that the amplification factor, which is the
magnitude of the roots of (22), must always be less than
(implies that the scheme is dissipative), or equal to one.
(17) After rewriting (22) in a simpler form
(23)
where the Courant–Friedrichs–Lewy (CFL) number [8] is
denoted by and where

(18)

In the equations above, the superscript on indicates that and is the complex conjugate of , (22) becomes easily
is evaluated at an integer multiple of the time step (that is at recognizable as quadratic in . Therefore, the roots of (23)
where is the time step). Similarly, the argument are
of indicates that is evaluated at an integer multiple of the
spatial step (that is at where is the spatial step).
This difference scheme is conditionally stable ( ).
After considerable experimentation, we discovered that
(24)
(16) can be approximated by the following second-order

Authorized licensed use limited to: Maharaja Agrasen Institute of Technology. Downloaded on April 19,2010 at 14:05:55 UTC from IEEE Xplore. Restrictions apply.
12 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 47, NO. 1, JANUARY 1999

While has not been analytically proven for any wave equation. This and an amplification factor of one are
range of center frequency, frequency bandwidth, time step, and the principle advantages of the complex-envelope formulation.
spatial step; numerical computations have shown that Some of the time-step advantage is lost when the complex-
for a wide range of center frequency, frequency bandwidth, envelope formulation is compared to the standard explicit
time step, and spatial step. Furthermore, the amplification FDTD method because of the additional calculations that
factor for every parameter combination tested. must be performed at each time step in the complex-envelope
This also indicates that this method does not suffer from the formulation.
numerical dispersion that is inherent in any multidimensional The comparison of the relative time step required in the
or multigridded explicit finite-difference scheme. standard wave- and the complex-envelope wave-equation so-
Implicit (Crank–Nicolson) finite-difference solutions can be lutions assumes that choosing to be the same in both
written in the following matrix form: formulations results in similar accuracy. Our initial results
confirm this. A fair comparisons between implicit and explicit
(25) methods is difficult because for explicit methods, the spatial
sampling rate and the time sampling rate are not independent
where and are constant, symmetric, and banded (tridi- and requirements for spatial sampling often force to be
agonal in the 1-D case) matrices, and where is a column much larger than is needed for an implicit solution method.
vector made up of the values of at all the spatial nodes at Additionally, numerical dispersion decreases the accuracy of
the time . Similarly, is a column vector made up of multidimensional and multigridded explicit finite-difference
source-related terms at all the spatial nodes at the time . schemes. However, one can conservatively state that in all
Equation (25) is used to determine at all of the spatial nodes cases
at time using the previously determined values of
at all spatial nodes for times and . Because (30)
%
and are constant and banded (tridiagonal in the 1-D case)
the computational burden of solving (25) is relatively small V. ANALYTIC SOLUTION OF THE COMPLEX WAVE EQUATION
compared to obtaining a solution of (25) if and were The solution of the 1-D LPL wave equation for a -directed
dense matrices. surface current at in an unbounded homogeneous
Although the numerical solution of the standard wave equa- region is as follows [12]:
tion and complex-envelope formulation of the wave equation
can both be obtained using (25), the required time step for
the numerical solution of the two different wave equations are
very different. For both equations, the selection of the time
step is based on the sampling theorem requirements. For the (31)
standard wave equation (SWE)
where is the characteristic impedance of the medium and
(26) is the unit step function. It is convenient to use a
bandpass-limited time function with a Gaussian envelope for
where is the number of samples per period. The sampling the source current. The source current used here is
theorem states that , however, numerical solutions
usually require . For the complex-envelope wave (32)
equation (CEWE)
where the Gaussian envelope function is
(27)
(33)
If in the solution of both wave equations is the same then
In (32), is scaled so that the peak value of
(28)
% is one. Note that the Gaussian-shaped envelope is a low-pass
limited time function with bandwidth . For reasons
where the percent bandwidth
related to the numerical solution of the wave equation, the
time delay is set as so that is on the order
% (29)
of 10 . Substituting (32) into (31) gives
For problems with small percent bandwidth sources, the
time step for the complex-envelope wave equation can be
several orders of magnitude greater than the time step for
the standard wave equation. As a result, the total number of
time steps required to obtain a solution to some maximum
time using the complex-envelope wave equation can be several (34)
orders of magnitude less than that required using the standard

Authorized licensed use limited to: Maharaja Agrasen Institute of Technology. Downloaded on April 19,2010 at 14:05:55 UTC from IEEE Xplore. Restrictions apply.
PURSEL AND GOGGANS: FDTD METHOD FOR SOLVING EM PROBLEMS WITH BANDPASS-LIMITED SOURCES 13

Rewriting the cosine terms in (34) as the real part of complex

exponential terms and then manipulating the result to obtain
the form of (10) yields

(35)

Equation (35) is the solution of the 1-D complex-envelope

wave equation for a directed complex-envelope surface
current at in an unbounded homogeneous region
where

(36)
Fig. 3. The BPL FDTD solution of the 1-D cavity problem.

A 1-D cavity with perfect electric conducting (PEC) walls

and a plane wave source centered between the walls is used as
an example problem in the following section. The PEC walls
are at and and the source is located at
. Using image theory and (35), the analytic solution
of the complex envelope wave equation can be determined. For
the plane wave source current of (36), the solution for
between the cavity walls is

(37)

where

(38)
Fig. 4. Comparison of the BPL FDTD and analytic solutions at 4.5 s.
and

and temporally . Note that the complex-

(39) envelope formulation is not particularly computationally ad-
vantageous for this combination of parameters
are the fields to the left and right of the cavity due to the , however, these values were used so that the spatial
source and its images. oscillations of the solution can be displayed clearly.
Fig. 3 illustrates the FDTD solution of the 1-D scalar
VI. EXAMPLE RESULTS complex-envelope wave equation. Because the cavity is sym-
A 1-D free-space filled cavity with perfect electric conduct- metric, the field distribution is only given for one half of the
ing (PEC) walls and a plane wave source centered between the cavity. At 1.5 and 3.0 s the traveling waves have not yet
walls is used as an example problem. The PEC walls are at propagated to the walls of the cavity. The rapid variation in
and where . The plane wave the field strength at 4.5 s is due to the backward traveling
source is simulated numerically by applying the impressed (reflected) wave interfering with the forward traveling wave.
current at the single spatial node corresponding to the position Given that Fig. 3 is a plot of the magnitude of the envelope of
. In the standard wave-equation solution, the impressed a radio frequency (RF) signal, this rapid variation is expected.
electric current is given by (32). In the complex-envelope A comparison of the analytic and BPL FDTD solutions is
wave equation, the impressed current is the complex-envelope shown in Fig. 4. The field distributions for 4.5 s at the edge
representation of as given in (36). of the cavity are plotted. The largest error is observed in this
The following example results have been computed with region of the cavity. The BPL FDTD solution is in excellent
(19), a center frequency, MHz and a bandwidth, agreement with the analytic solution. The maximum observed
MHz. The continuous signals were sampled spatially error in the complex-valued solution is less than 2%.

Authorized licensed use limited to: Maharaja Agrasen Institute of Technology. Downloaded on April 19,2010 at 14:05:55 UTC from IEEE Xplore. Restrictions apply.
14 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 47, NO. 1, JANUARY 1999

Fig. 5. The BPL FDTD solution of the 1-D cavity problem.. Fig. 7. The BPL FDTD solution of a region with a dielectric interface.

VII. CONCLUSIONS
Using the complex-envelope representation of bandpass-
limited signals, a new bandpass-limited formulation of the
vector wave equation was formulated and presented here. A
new bandpass-limited FDTD scheme that solves the bandpass-
limited vector wave equation has been developed and pre-
sented. An amplification factor of one and the ability to set
the time-step size based solely on the bandwidth of the system
are the principle advantages of this new method. For narrow-
band (small percent bandwidth) systems, the time step can
be several orders of magnitude larger than the time step in
conventional FDTD methods.
Results that illustrate the accuracy of this new method
were presented. While these results were for a 1-D cavity,
the formulation is valid in two and three dimensions as
Fig. 6. Comparison of the BPL FDTD and analytic solutions at 14.5 s. well. Preliminary two-dimensional (2-D) results indicate that
the method is numerically stable and accurate. The authors
The field distributions at later times are illustrated in Figs. 5 are currently working on a comprehensive error analysis of
and 6. In Fig. 5, the field distribution at 14.5 s has propagated the 1-D formulation, application of an absorbing boundary
more than 125 m and is in the process of completing its third condition, and the verification of 2-D results.
reflection. A comparison with the analytic solution at 14.5 s
is shown in Fig. 6. While the error in the BPL FDTD solution REFERENCES
has increased, it is still less than 5%. This cumulative error
[1] S. Haykin, Communication Systems, 2nd ed. New York: Wiley, 1983.
is due to the finite-difference approximations of the temporal- [2] P. M. Goggans and J. D. Pursel, “Radar time- and frequency-domain
and spatial-derivatives. received signals for realistic antennas and scatterers,” in 12th Annu. Rev.
A dielectric interface is used as a second example problem. Progress Appl. Computat. Electromagn., Monterey, CA, 1996, vol. 12.
[3] A. Papoulis, The Fourier Integral and Its Applications. New York:
To the left of the interface is free space. To the McGraw-Hill, 1962.
right the relative permittivity is two. Note that the [4] J. D. Pursel, “A finite-difference time-domain solution of the second-
spatial dimension of the 1-D cavity was extended so that dur- order bandpass-limited electromagnetic wave equation,” Ph.D. disserta-
tion, Univ. Mississippi, 1998.
ing the times of interest, there would be no reflection from the [5] P. M. Goggans and J. D. Pursel, “Using time-domain complex-envelope
PEC walls. The frequency and sampling parameters remain the representations of band-pass-limited signals in the finite-difference solu-
tion of the wave equation,” in North Amer. Radio Sci. Meet., Montreal,
same as in the first example. Fig. 7 illustrates the BPL FDTD Canada, July 1997, p. 128.
computed field distributions for various times in the region of [6] R. Courant and D. Hilbert, Methods of Mathematical Physics— Vol. II.
the interface. At 1.5 s, the field has not yet propagated to New York: Intersci., 1962.
[7] E. Isaacson and H. B. Keller, Analysis of Numerical Methods. New
the interface. The field is in the process of being transmitted York: Wiley, 1966.
through and reflected from the interface at 3.0 s. At 4.5 s, [8] R. Courant, K. Friedrichs, and H. Lewy, “On the partial difference
the reflected field is propagating to the left, while the transmit- equations of mathematical physics,” IBM J. Res. Develop., vol. 11, no.
1, pp. 215–234, 1967; “Uber die partiellen differenzengleichungen der
ted field is propagating to the right. The reflection and trans- mathematischen physik,” Math. Ann., vol. 100, pp. 32–74, 1928 (Engl.
mission coefficients are within 0.1% of the analytic values. transl.).

Authorized licensed use limited to: Maharaja Agrasen Institute of Technology. Downloaded on April 19,2010 at 14:05:55 UTC from IEEE Xplore. Restrictions apply.
PURSEL AND GOGGANS: FDTD METHOD FOR SOLVING EM PROBLEMS WITH BANDPASS-LIMITED SOURCES 15

[9] B. Gustafsson, H. Kreiss, and J. Oliger, Time Dependent Problems and P. M. Goggans (S’78–M’89) was born in Opelika,
Difference Methods. New York: Wiley, 1995. AL, on September 13, 1954. He received the B.S.
[10] P. D. Lax, “Hyperbolic difference equations: A review of the and M.S. degrees in electrical engineering in 1976
Courant–Friedrichs–Lewy paper in the light of recent developments,” and 1978, respectively and the Ph.D. degree in 1990,
IBM J. Res. Development, vol. 11, no. 1, pp. 235–238, 1967. all from Auburn University, Auburn, AL.
[11] R. D. Richtmyer and K. W. Morton, Difference Methods for Initial-Value From 1979 to 1985, he was employed by Sandia
Problems, 2nd ed. New York: Wiley, 1967. National Laboratories, Albuquerque, NM, in the
[12] I. V. Lindell, Methods for Electromagnetic Field Analysis. Oxford, Radar Signal Analysis Division. From 1985 to 1990
U.K.: Oxford Univ. Press, 1992. he was with Auburn University as an Instructor
while working toward the Ph.D. degree. In 1990 he
was appointed Assistant Professor in the Department
of Electrical Engineering at the University of Mississippi, Oxford, MS.
In 1994 he was promoted to Associate Professor. His research interests
J. D. Pursel (S’90–M’99) was born in Bakersfield, CA, on July 3, 1969. He include computational electromagnetics methods based on integral equation
received the B.S.E.E. degree from the California State University at Fresno, and differential equation formulations, the application of signal processing
in 1993, and the M.S. and Ph.D. degrees in engineering science from the techniques to the development of computational solutions for electromagnetic
University of Mississippi, University, MS, in 1995 and 1998, respectively. applications, and the electromagnetic and signal processing aspects of radar
He served as a Graduate Instructor and Research Assistant in the De- systems.
partment of Electrical Engineering at the University of Mississippi under Dr. Goggans is a member of Eta Kappa Nu, Phi Kappa Phi, the Audio Engi-
a Graduate Assistance Areas of National Need Fellowship. In May 1998 neering Society, and the IEEE Antennas and Propagation, Microwave Theory
he joined the Electromagnetics Directorate of the Riverside Research In- and Techniques, Electromagnetic Compatibility, and Education Societies.
stitute, New York, NY. His current research interests are optimization,
finite-difference time-domain methods, and the use of high-performance
computing in the solution of electromagnetics problems.
Dr. Pursel is a member of the Eta Kappa Nu and Tau Beta Pi honor
Societies, the Sigma Xi Society of scientific research, the Antennas and
Propagation, Electromagnetic Compatibility, and Microwave Theory and
Techniques Societies of IEEE, and the Society of Industrial and Applied
Mathematics.

Authorized licensed use limited to: Maharaja Agrasen Institute of Technology. Downloaded on April 19,2010 at 14:05:55 UTC from IEEE Xplore. Restrictions apply.