The FDTD Grid and The Yee Algorithm PDF

4 The FDTD grid and the Yee algorithm
After a brief exposure to different finite difference algorithms and methods, we now
focus our attention on the so-called FDTD algorithm, or alternatively the Yee algorithm
[1], for time domain solutions of Maxwell’s equations. In this algorithm, the continuous
derivatives in space and time are approximated by second-order accurate, two-point
centered difference forms; a staggered spatial mesh is used for interleaved placement
of the electric and magnetic fields; and leapfrog integration in time is used to update
the fields. This yields an algorithm very similar to the interleaved leapfrog described in
Section 3.5.
The cell locations are defined so that the grid lines pass through the electric field
components and coincide with their vector directions, as shown in Figure 4.1. As a
practical note, the choice here of the electric field rather than the magnetic field is
somewhat arbitrary. However, in practice, boundary conditions imposed on the electric
field are more commonly encountered than those for the magnetic field, so that placing the
mesh boundaries so that they pass through the electric field vectors is more advantageous.
We will have more to say about the locations of field components later in this chapter.
Note that the Yee cell depicted in [2, Fig. 3.1] has the cell boundaries to be aligned
with the magnetic field components, rather than with the electric field components as
in Figure 4.1. This choice of the particular way of associating the spatial indices i, j,
and k with the field quantities is obviously arbitrary and should not matter in the final
analysis, as long as one is consistent.
The fundamental unit of our 3D grid, known as the Yee cell, is shown in Figure 4.1.
With the Yee cell so defined, the spatial derivatives of various quantities are evaluated
using the simple two-point centered difference method. For example, the z derivative
of any given field component evaluated at time nt and at the mesh point (i, j, k) is
given as

∂ n i,n j,k+1/2 − i,n j,k−1/2
= + O[(z)2 ]. (4.1)
∂z i, j,k z
This implies that the field is defined at integer grid points in x and y (i.e., i
and j) and at half grid points in z (i.e., z + 1/2). Now consider a real field, for
example, the z component of the magnetic field, z , which is defined at mesh points
(i ± 1/2, j ± 1/2, k) as shown in Figure 4.1. The derivative of z with respect to x at
Downloaded from https://www.cambridge.org/core. University of Sussex Library, on 28 Jun 2018 at 04:27:51, subject to the Cambridge Core terms of use, available at
https://www.cambridge.org/core/terms. https://doi.org/10.1017/CBO9780511921353.005
The FDTD grid and the Yee algorithm 73
y
x z
(i–1, j+1, k+1)
(i, j, k+1)
x
y
x z
z y
y
x
z (i–1, j+1, k)
(i, j, k) z x
y
(i, j+1, k)
x y
Figure 4.1 Placement of electric and magnetic field components in a 3D staggered mesh, known
as the Yee cell. The small vectors with thick arrows are placed at the point in the mesh at which
they are defined and stored. For example, y is defined/stored at mesh points (i + m 1 , j +
+ 1/2m 2 , k + m 3 ), where m 1,2,3 = 0, ±1, ±2, . . . , while y is defined/stored at mesh points
(i − 1/2 + m 1 , j + m 2 , k + 1/2 + m 3 ), where m 1,2,3 = 0, ±1, ±2, . . . .
the point (i, j + 1/2, k) (at which, by the way, y is defined) and at time step n is given
by:
n n
∂z n z i+1/2, j+1/2,k −z i−1/2, j+1/2,k
= + O[(x)2 ]. (4.2)
∂ x i, j+1/2,k x
The derivatives in time are also discretized with second-order centered difference
approximations such as that given in Equation (3.49), but the updating of and are
staggered in time by one half time step; in fact, the components of are defined at half
time steps, n + 1/2. In vector form we can write the discretized Maxwell’s equations
as:
n+1 n
∂ n+1/2 − 1
= [∇ × ]n+1/2
∂t t
n+1 n+3/2 n+1/2
∂ − 1
= − [∇ × ]n+1 .
∂t t μ
Leapfrog in time n+1
n+3/2
n t

n+1/2

Figure 4.2 Leapfrog time marching of electric and magnetic field vectors.
Rearranging terms, we have the time-marching procedure:
n+1 tn
= + [∇ × ]n+1/2 (4.3a)

n+3/2 n+1/2 t
= − [∇ × ]n+1 . (4.3b)
μ
This leapfrog updating in time is identical to that given in Figure 3.11 for voltage
instead of and current instead of . This is depicted in Figure 4.2.
The next step in the derivation of the FDTD algorithm is to discretize the spatial
components on the right-hand side of Equations (4.3). In Section 4.6, we will look
at cylindrical and spherical coordinate systems; for simplicity, however, we begin with
Cartesian coordinates. Maxwell’s equations in a source-free, simple medium in Cartesian
coordinates are given by:

∂x 1 ∂z ∂ y ∂x 1 ∂ y ∂z
= − = −
∂t ∂y ∂z ∂t μ ∂z ∂y

∂ y 1 ∂x ∂z ∂ y 1 ∂z ∂x
= − = − (4.4)
∂t ∂z ∂x ∂t μ ∂x ∂z

∂z 1 ∂ y ∂x ∂z 1 ∂x ∂ y
= − = − .
∂t ∂x ∂y ∂t μ ∂y ∂x
Next, we will look at the discretization of these equations, starting with 1D and
building up to the 2D and 3D cases.
4.1 Maxwell’s equations in one dimension
The simplest expression of Maxwell’s equations in the FDTD algorithm can be found
by limiting ourselves to the consideration of electromagnetic fields and systems with
no variations in two dimensions, namely, both z and y. We thus drop all of the y and
z derivatives in Equations (4.4) above. For now, we neglect magnetic or electric losses
4.1 Maxwell’s equations in one dimension 75
(i.e., let σ = 0 and σm = 0) and assume simple and source-free media, so that ρ̃ = 0,
= 0, and and μ are simple constants, independent of position (i.e., the medium
is homogeneous), direction (isotropic), or time. With these simplifications, Maxwell’s
Equations (2.1) and (2.3) simplify to:
∂x ∂x
=0 =0
∂t ∂t
∂ y 1 ∂z ∂ y 1 ∂z
=− = (4.5)
∂t ∂x ∂t μ ∂x
∂z 1 ∂ y ∂z 1 ∂ y
= =− .
∂t ∂x ∂t μ ∂x
Note that two of the equations have disappeared completely; the remaining equations
involve only derivatives with respect to x and with respect to t. We can now organize
the remaining four equations into pairs that are completely independent of each other:
∂ y 1 ∂z
1D TM = (4.6a)
∂t μ ∂x
∂z 1 ∂ y
= (4.6b)
∂t ∂x
and
∂ y 1 ∂z
1D TE =− (4.7a)
∂t ∂x
∂z 1 ∂ y
=− . (4.7b)
∂t μ ∂x
For example, note that Equations (4.6) involve x- and t-derivatives of z and y ,
but do not involve z or y ; the opposite is true of Equations (4.7). The first pair are
referred to as the TM mode, and the second comprise the TE mode.1 The placement
of the electric and magnetic field vectors on a one-dimensional Yee grid is shown in
Figure 4.3 for both TM and TE cases.
The corresponding finite difference equations are found by applying second-order
centered differences to both the time and space derivatives in Equations (4.6) and (4.7).
As shown in Figure 4.3, the components of are located at integer grid points in x, and
integer grid points in time t as well; the components of are located at half grid points.
1 Note that these TM and TE mode definitions are different from those in classic waveguide analysis. The
meaning of the nomenclature will become clear when we discuss the 2D FDTD algorithm.
z z
TE
TM
z
y z y
y y
z i
i
x x
i+1/2 i+1/2
z
i+1 y
z y i+1
y
(a) (b)
Figure 4.3 The placement of electric and magnetic field vectors for Maxwell’s equations in one
dimension. Both the TM and TE cases are shown.
After rearranging, we find:
1D TE mode
n+1 n n+1/2
t n+1/2
y = y − z − z
i i i x i+1/2 i−1/2
n+1/2 n−1/2 n (4.8)
t n
z = z − y − y
i+1/2 i+1/2 μi+1/2 x i+1 i
and
1D TM mode
n+1/2 n−1/2 t n n

y = y + z − z
i+1/2 i+1/2 μi+1/2 x i+1 i
n+1 n (4.9)
t n+1/2 n+1/2
z = z + y − y
i i i x i+1/2 i−1/2
Note that Equations (4.8) and (4.9) are identical to the interleaved leapfrog algorithm
given in Equation (3.86) in terms of voltage and current except for the substitution
of for C and μ for L.
A comment is necessary about starting the FDTD simulation. If we begin with
an initial electric field distribution 0z over the simulation space, then Equation (4.9)
suggests that we need −1/2 y in order to proceed with updating. Usually, one simply
estimates y instead, and proceeds to the update equation for nz . However, the accuracy
1/2
of the simulation is highly dependent on the first time step 1/2

y , so higher-order methods
are required to make that estimate.
4.1 Maxwell’s equations in one dimension 77
Sinusoidal Source at x = 0, Snapshot at t = 50 ns
1
z field (V/m)
0.5
−0.5
−1
−10 −5 0 5 10
Gaussian Pulse Source at x = 0
25.2 ns 10.8 ns t = 3.6 ns 18.0 ns t = 32.4 ns

1
z field (V/m)
0.5
0
−10 −5 0 5 10
x (meters)
Figure 4.4 1D FDTD simulation examples. The top panel shows a snapshot of z at t = 50 ns of
a simulation using a sinusoidal source at the center of the grid. The lower panel shows
successive snapshots of z where the source is a Gaussian pulse.
This problem goes away in most circumstances, when the simulation space is assumed
to have zero fields initially; then both 0z and −1/2
y are zero everywhere.
4.1.1 Example 1D simulations

Figure 4.4 shows example simulations in a simple, lossless 1D space for the TM mode.
In both simulations, t = 0.18 ns, x = v p t = 5.4 cm, and = 0 , μ = μ0 . The
grid is 401 cells in length, giving just under 22 meters in total length.
In the topmost panel, the space is excited by a sinusoidal source at the center point
given by z,in (n) = 0 sin(2π f 0 nt), with a center frequency f 0 = 278 MHz, corres-
ponding to exactly 20 grid cells per wavelength. In this way we resolve the wave
sufficiently to avoid numerical dispersion, which will be discussed in Chapter 6. The
amplitude of the input wave is 0 = 1 V/m. In this example, the wave has propagated
beyond the simulation space boundaries, but we have used an absorbing boundary
condition (ABC) to truncate the simulation space; these will be described in Chapters 8
and 9.
In the lower panel, the space is excited by a Gaussian pulse source at the center point
given by z,in (n) = 0 e−(n−8) /4 , i.e., the pulse has a 1/e half-width of four time steps,
2 2
and ramps up smoothly from zero. The amplitude is given as 0 = 1 V/m. The figure
shows that our Gaussian pulse propagates to the left and right at the speed of light:
specifically, at t = 32.4 ns the pulse has traveled exactly 9.72 meters, giving a numerical
phase velocity of 3 × 108 m/s.
4.2 Maxwell’s equations in two dimensions
We now proceed with discretization of Maxwell’s equations based on the staggered

mesh placement of the fields as given in Figure 4.1 for the 2D case. For this purpose,
we assume that there are no variations of either the fields or the excitation in one of
the directions, say the z direction. Thus, all derivatives with respect to z drop out from
the two curl equations. Again neglecting magnetic or electric losses (i.e., let σ = 0 and
σm = 0) and assuming simple and source-free media, so that ρ̃ = 0, = 0, and and μ
are simple constants, independent of position, direction, or time, Maxwell’s Equations
(2.1) and (2.3) can be written as:
∂x 1 ∂z
=− (4.10a)
∂t μ ∂y
∂ ∂ y 1 ∂z
−μ =∇ × → = (4.10b)
∂t ∂t μ ∂x

∂z 1 ∂x ∂ y
= − (4.10c)
∂t μ ∂y ∂x
and
∂x 1 ∂z
= (4.11a)
∂t ∂y
∂ ∂ y 1 ∂z
=∇ × → =− (4.11b)
∂t ∂t ∂x

∂z 1 ∂ y ∂x
= − . (4.11c)
∂t ∂x ∂y
As we did with the 1D case, we can group Equations (4.10) and (4.11) according to
field vector components, one set involving only x , y , and z , and another set involving
x , y , and z , referred to respectively as the TE and TM modes. The resulting sets of
equations are given as:
∂x 1 ∂z
=− (4.12a)
∂t μ ∂y
∂ y 1 ∂z
TM mode = (4.12b)
∂t μ ∂x

∂z 1 ∂ y ∂x
= − (4.12c)
∂t ∂x ∂y
4.2 Maxwell’s equations in two dimensions 79
and
∂x 1 ∂z
= (4.13a)
∂t ∂y
∂ y 1 ∂z
TE mode =− (4.13b)
∂t ∂x

∂z 1 ∂x ∂ y
= − . (4.13c)
∂t μ ∂y ∂x
It is important to note that the TE and TM modes above are completely uncoupled
from one another; in other words, they contain no common vector field components and
can therefore exist independently from one another.
This particular TE and TM designation of electromagnetic wave types is quite unre-
lated to the separation into TE and TM modes of classic guided electromagnetic waves.
For time-harmonic guided waves, the z direction is generally taken to be the direction
of propagation, with all field components thus exhibiting variations in z in the form
of e− jβz , where β is the phase constant for the particular waveguide, as determined by
the boundary conditions. Thus, all field components do in fact vary with z, and the TE
(or TM) modes are identified as those waves which have a magnetic (or electric) field
component in the axial (i.e., z) direction. For the 2D electromagnetic waves in the FDTD
grid, none of the field components varies with z so that there certainly is no propagation
in the z direction; instead, the TE or TM modes are identified on the basis of whether
the wave has a nonzero or component in the third dimension, which we decided
had no variation – in this case the z dimension. Hence, Equations (4.12) are identified
as the TM mode on account of the z component.
Another way of looking at the TM and TE modes is as follows: since we have decided
that there is no variation in the z direction (i.e., all derivatives with respect to z have been
set to zero), there can be no propagation in the z direction; but, in general, propagation
in either x- or y-directions (or both) is possible. In the TM mode, only x and y
components are nonzero, and are in the plane of propagation. The reverse is true in the
TE mode.
Although Equations (4.12) and (4.13) are entirely symmetric, the TE and TM modes
are physically quite different, due to orientation of electric and magnetic field lines with
respect to the infinitely long axis (i.e., the z axis) of the system. Note, for example, that
the electric fields for the TE mode lie in a plane perpendicular to the long axis, so that
when modeling infinitely long metallic structures (aligned along the long axis), large
electric fields can be supported near the metallic surface. For the TM mode, however,
the electric field is along the z direction and thus must be nearly zero at infinitely long
metallic surfaces (which must necessarily be aligned with the infinite dimension in the
system) being modeled. As such, when running 2D simulations, the choice of TE or TM
modes can depend on the types of structures in the simulation.
We now separately consider the discretization of Equations (4.13) and (4.12), respec-
tively, for TE and TM modes.
(i, j)
x y
x z y
i+1/2
j+1/2 (i+1, j+1)
Figure 4.5 An FDTD unit cell for transverse electric (TE) waves. The small vectors with thick
arrows are placed at the point in the mesh at which they are defined and stored. For example, y
is defined/stored at mesh point (i, j + 1/2), while z is defined/stored at mesh points
(i + 1/2, j + 1/2).
4.2.1 Transverse electric (TE) mode

A portion of the Yee cell constituting a unit cell for the TE case is depicted in Figure 4.5.
Note the spatial positions of the two electric field components and the single magnetic
field component. The x component is located at half x and integer y grid points, i.e.,
(i + 1/2, j), while the y component is located at integer x and half y grid points, i.e.,
(i, j + 1/2). The magnetic field component z is located at half x and half y grid points,
i.e., (i + 1/2, j + 1/2). These grid points are chosen to accommodate the interleaved
leapfrog algorithm that makes up the FDTD method; however, one should be able to
identify the trend that makes it easy to remember. This will be spelled out in the next
section, when we look at the 3D FDTD method.
The spatially discretized versions of the component Maxwell’s Equations (4.13) are:
⎡ ⎤

z − z
∂x 1⎢ i+1/2, j+1/2 i+1/2, j−1/2 ⎥
= ⎣ ⎦ (4.14a)
∂t i+1/2, j y
⎡ ⎤

z − z
∂ y 1⎢ i+1/2, j+1/2 i−1/2, j+1/2 ⎥
=− ⎣ ⎦ (4.14b)
∂t i, j+1/2 x
⎡ ⎤

x − x y − y
∂z 1⎢ i+1/2, j+1 i+1/2, j i+1, j+1/2 i, j+1/2 ⎥
= ⎣ − ⎦.
∂t i+1/2, j+1/2 μ y x
(4.14c)
z
y x
(i, j)
x y
(i, j+1)
(i+1, j) z
i+1/2
j+1/2 x
y
Figure 4.6 An FDTD unit cell for transverse magnetic (TM) waves. The small vectors with thick
arrows are placed at the point in the mesh at which they are defined and stored. For example, z
is defined/stored at mesh point (i, j), while y is defined/stored at mesh points (i + 1/2, j).
Using (4.3) for the time derivatives, we arrive at the FDTD algorithm for TE waves:
2D TE mode
n+1/2 n−1/2 t

z = z +
i+1/2, j+1/2 i+1/2, j+1/2 μi+1/2, j+1/2
⎡ n n n n ⎤

x −x y − y
⎢ i+1/2, j+1 i+1/2, j i+1, j+1/2 i, j+1/2 ⎥
×⎣ − ⎦
y x
n+1 n n+1/2 n+1/2

t
x = x + z − z
i+1/2, j i+1/2, j i+1/2, j y i+1/2, j+1/2 i+1/2, j−1/2
n+1 n n+1/2 n+1/2
t
y = y − z − z
i, j+1/2 i, j+1/2 i, j+1/2 x i+1/2, j+1/2 i−1/2, j+1/2
(4.15)
where we have explicitly acknowledged that the medium parameter values (, μ) to be
used are those at the spatial point at which the quantity is being updated. In this case, the
medium does not need to be homogeneous: the values of and μ can change throughout
the simulation space. However, this algorithm still requires that the medium is isotropic,
i.e., the parameters and μ do not vary with direction.
4.2.2 Transverse magnetic (TM) mode

A portion of the Yee cell constituting a unit cell for the TM case is depicted in Figure 4.6.
The spatially discretized equations for the TM mode can be written in an analogous
fashion to (4.14), by noting that the discrete electric field variables are located at the cell
edges, and including an appropriate half-cell shift. Or, one can use the pattern identified
for the TE mode (spelled out in the next section) to show that the z component is
located at integer grid points (i, j); the x component is located at integer x and half y
grid points, i.e., (i, j + 1/2); and the y component is located at half x and integer y
components, i.e., (i + 1/2, j).
The spatially discretized versions of Equations (4.12) are:
⎡ ⎤

− z
∂x 1 ⎢ z i, j+1 i, j ⎥
=− ⎣ ⎦ (4.16a)
∂t i, j+1/2 μ y
⎡ ⎤

z − z
∂ y 1⎢ i+1, j i, j ⎥
= ⎣ ⎦ (4.16b)
∂t i+1/2, j μ x
⎡ ⎤

y − y x − x
∂z 1⎢ i+1/2, j i−1/2, j i, j+1/2 i, j−1/2 ⎥
= ⎣ − ⎦. (4.16c)
∂t i, j x y
Again using analogies to Equations (4.3) for the time derivatives, we arrive at the
FDTD algorithm for TM waves:
2D TM mode
n+1/2 n−1/2 n
t n
x = x − z − z
i, j+1/2 i, j+1/2 μi, j+1/2 y i, j+1 i, j
n+1/2 n−1/2 n
t n
y = y + z − z
i+1/2, j μi+1/2, j x
i+1/2, j i+1, j i, j
⎡ n+1/2 n+1/2 n+1/2 n+1/2 ⎤

n+1 n y − y x − x
t ⎢ i+1/2, j i−1/2, j i, j+1/2 i, j−1/2 ⎥
z = z + ⎣ − ⎦
i, j i, j i, j x y
(4.17)
where we have again explicitly denoted the medium parameter values (, μ) as those
at the spatial point at which the quantity is being updated, removing the restriction of
homogeneity.
4.2.3 Example 2D simulations

Figure 4.7 shows example simulations
√ in a simple, lossless 2D space. In both simulations,
x = 5.4 cm, t = x/v p / 2 = 0.13 ns,2 and = 0 , μ = μ0 . The grid is 201 × 201
√
2 Note the 2 in this formula for t; this is a result of operating in 2D. We will derive this stability restriction
formally in Chapter 5.
4 Gaussian Source Sinusoidal Source

t = 50 Δt t = 50 Δ t
2
y (meters)
–2
–4
4 Gaussian Source Sinusoidal Source

t = 100 Δt t = 100 Δ t
2
y (meters)
–2
–4
–4 –2 0 2 4 –4 –2 0 2 4
x (meters) x (meters)
Snapshots at t = m × 20 Δ t Snapshot at t = 200 Δ t
0.2 r–1/2 0.8

trend line
z-field (V/m)
z-field (V/m)
0.4 r–1/2
0.1 trend line
0
0
−0.4
−4 −2 0 2 4 −4 −2 0 2 4
x (meters) x (meters)
Figure 4.7 2D FDTD simulation examples. In the left column, the space has been excited by the
derivative of Gaussian source at the center of the space; in the right column,
√ a sinusoidal source.
The lower panels show a slice through the center, demonstrating the 1/ r decay of the fields
with distance.
cells, giving ∼11 meters in each dimension; the figure zooms in to the range from −5 to
5 meters.
In the left column, the space is excited by the derivative of a Gaussian pulse at
the center point given by z,in = 0 e−(n−8) /4 , i.e., the pulse has a 1/e half-width of
2 2
four time steps, and ramps up smoothly from zero. Note that excitation with a Gaussian
Table 4.1 Field component locations in the Yee cell
Field x location y location z location
x i + 1/2 j k
y i j + 1/2 k
z i j k + 1/2
x i j + 1/2 k + 1/2
y i + 1/2 j k + 1/2
z i + 1/2 j + 1/2 k
leads to unphysical fields left in the space due to “grid capacitance” [2]. For this reason it
is prudent to use the first derivative of a Gaussian as an excitation source for broadband
simulations.
In the right column, the center point is excited by a sinusoidal source at 278 MHz,
ensuring exactly 20 grid cells per wavelength. The figures show the magnitude of the z
field (i.e., the absolute value), so each concentric ring is half a wavelength. The bottom
panels show slices through the center of the expanding concentric rings; in the case of
the Gaussian source, we show snapshots of the expanding field at successive times, while
for the sinusoidal source we show a single snapshot. In both cases, we can see that the
√
fields decay with radius at exactly 1/ r : recall that for cylindrical symmetry the power
√
decays as 1/r , so each of the fields and decay as 1/ r .
4.3 FDTD expressions in three dimensions
We have hinted up to this point that there is a pattern in the location of field values in the
Yee cell. Looking at Figure 4.1, and writing down the locations of each of the six field
components, we arrive at Table 4.1. The pattern is as follows:
r For the components, for a = x, y, z, the a-location is half-integer, while the other
a
two are integer.
r For the components, it is exactly the opposite: the a-location is integer, while the
a
other two are half-integer.
This is further illustrated in Figure 4.8, which expands Figure 4.1 and shows two
slices in the x-plane, at i and i − 1/2. These two plane slices are the same at i ± m and
i ± m − 1/2, for any integer m. Note that these two slices are, in fact, the TE and TM
mode patterns in the y-z plane.
For completeness, we present here the full FDTD expressions for Maxwell’s equa-
tions in simple media in three dimensions, derived by discretizing Maxwell’s Cartesian
4.3 FDTD expressions in three dimensions 85
(i, j, k+1) z
y
i plane y
x z (i–1, j+1, k+1)
(i, j, k+1)
z x z x
y
x z
y z y
(i, j, k) (i, j+1, k)
y
x
z (i–1, j+1, k)
(i–1/2, j, k+1)
x z x z
i–1/2 plane (i, j, k) x
y
(i, j+1, k)
y y z
x y
z x x y
x
(i–1/2, j, k) (i–1/2, j+1, k)
Figure 4.8 Yee cell in 3D. The cube at right is a repeat of Figure 4.1; at left, slices of the cube are
shown at i and i − 1/2, showing the interleaved field placement in each plane.
component Equations (4.4) with second-order centered differences:
n+1 n

x = x (4.18a)
i+1/2, j,k i+1/2, j,k
⎡ n+1/2 n+1/2

z −z
t ⎢ i+1/2, j+1/2,k i+1/2, j−1/2,k
+ ⎣
i+1/2, j,k y
n+1/2 n+1/2 ⎤

y − y
i+1/2, j,k+1/2 i+1/2, j,k−1/2⎥
− ⎦
z
n+1 n

y = y (4.18b)
i, j+1/2,k i, j+1/2,k
⎡ n+1/2 n+1/2

x −x
t ⎢ i, j+1/2,k+1/2 i, j+1/2,k−1/2
+ ⎣
i, j+1/2,k z
n+1/2 n+1/2 ⎤

z −z
i+1/2, j+1/2,k i−1/2, j+1/2,k ⎥
− ⎦
x
n+1 n

z = z (4.18c)
i, j,k+1/2 i, j,k+1/2
⎡ n+1/2 n+1/2

y − y
t ⎢ i+1/2, j,k+1/2 i−1/2, j,k+1/2
+ ⎣
i, j,k+1/2 x
n+1/2 n+1/2 ⎤

x −x
i, j+1/2,k+1/2 i, j−1/2,k+1/2 ⎥
− ⎦
y
n+1/2 n−1/2

x = x (4.19a)
i, j+1/2,k+1/2 i, j+1/2,k+1/2
⎡ n n

y − y
t ⎢ i, j+1/2,k+1 i, j+1/2,k
+ ⎣
μi, j+1/2,k+1/2 z
n n ⎤

z −z
i, j+1,k+1/2 i, j,k+1/2 ⎥
− ⎦
y
n+1/2 n−1/2

y = y (4.19b)
i+1/2, j,k+1/2 i+1/2, j,k+1/2
⎡ n n

z −z
t ⎢ i+1, j,k+1/2 i, j,k+1/2
+ ⎣
μi+1/2, j,k+1/2 x
n n ⎤

x −x
i+1/2, j,k+1 i+1/2, j,k ⎥
− ⎦
z
n+1/2 n−1/2

z = z (4.19c)
i+1/2, j+1/2,k i+1/2, j+1/2,k
⎡ n n

x −x
t ⎢ i+1/2, j+1,k i+1/2, j,k
+ ⎣
μi+1/2, j+1/2,k y
n n ⎤

y − y
i+1, j+1/2,k i, j+1/2,k ⎥
− ⎦.
x
Note that the above expressions are for the case where we have no electric or magnetic
current sources, i or i . Inclusion of the sources is straightforward and will be discussed
in later chapters. Together with the discretization pattern outlined in Table 4.1, the above
4.3 FDTD expressions in three dimensions 87
z Amplitude (dB)
−10
5.4
−20
2.7
−30
z (m) 0
−40
−2.7
−50
−5.4
5.4
2.7 5.4 −60
0 2.7
x (m) 0
−2.7 y (m)
−2.7 −70
−5.4 −5.4
Figure 4.9 3D FDTD simulation example. The sinusoidal source wave expands radially outward
from the center of the space, and the field amplitude decays as 1/r .
Equations (4.18) and (4.19) can be used to make the following observations about the
discretization process:
r The update equations for components are discretized in time at n + 1/2, so that
they use the time difference of the fields at times n + 1 and n.
r The update equations for components are discretized at n, so that they use the time
difference of the fields at times n + 1/2 and n − 1/2.
r In each equation, the values of μ and/or are used at the same location as the field
component being updated. This arises because the true Maxwell’s equations relate
to and μ to ; as such, and μ must be discretized along with the fields.
4.3.1 Example 3D simulation

Figure 4.9 shows an example simulation in 3D. This example uses a grid of 100 × 100 ×
100 grid cells, yielding 106 locations total. One can see from this relatively small
example that in 3D the spaces can quickly become prohibitively large. We excite a
sinusoidal z source at the center of the space with frequency 278 MHz just as in the
2D example earlier; however, in this case we have only 10 grid cells per wavelength, so
x = y = z = 10.8 cm.
The figure shows slices of the 3D z field amplitude through the center in each
dimension; this plot is created with the Matlab function slice. Because of the z
source at the center of the simulation space, the characteristic dipole radiation pattern
is evident, with very little field intensity radiated in the z direction above and below the
source. In the radial direction in the x–y plane, the field decays as 1/r as expected, since
the power decays as 1/r 2 .
4.4 FDTD algorithm for lossy media
We have so far derived the 1D, 2D, and 3D FDTD algorithms for simple media, by
neglecting the loss terms in the two curl equations. Incorporation of these losses into the
FDTD algorithm is straightforward. These losses are denoted by electric conductivity σ ,
and sometimes magnetic conductivity σm , which are parameters of the medium just like
and μ. The magnetic conductivity σm is, strictly speaking, not physically realistic, but
is often used for artificial absorption at simulation space boundaries, as will be discussed
in Chapter 9; as such, we include it here. Note that like and μ in Equations (4.18) and
(4.19) above, the conductivities σ and σm may vary throughout the space.
4.4.1 1D waves in lossy media: waves on lossy transmission lines

We start with the 1D case, which is posed in terms of uniform transmission line equations;
again, one should note that the algorithms we derive for and are directly applicable to
the 1D Maxwell’s equations by simple substitution of variables. In this case, the electric
and magnetic conductivities are replaced by electrical resistance R and conductance G.
Derivation begins with the lossy transmission line Equations (3.1), repeated here for
convenience:
∂ ∂
= −R − L (4.20a)
∂x ∂t
∂ ∂
= −G − C . (4.20b)
∂x ∂t
The discretization of these equations follows as before. The only difference is in
the first term on the right-hand sides of both Equations (4.20) above; these terms are
not derivatives, which means they can be discretized directly. The finite difference
approximations of these equations using the interleaved leapfrog algorithm are given by
modified versions of Equations (3.86):
Discretized at (i + 1/2, n):
n+1/2 n−1/2
i+1
n
− in i+1/2 − i+1/2
= −R i+1/2
n
−L . (4.21a)
x t
Discretized at (i, n + 1/2):
n+1/2 n+1/2
i+1/2 − i−1/2 n+1/2 in+1 − in
= −G i −C . (4.21b)
x t
4.4 FDTD algorithm for lossy media 89
Note that Equation (4.21a) is discretized at time step n and spatial step i + 1/2, while
Equation (4.21b) is discretized at time step n + 1/2 and spatial step i. Noting that the
current and voltage are advanced in time in a leapfrog manner, there is a problem in
using Equation (4.21a) directly, since the current is needed at time step n, but is only
available at half time steps n + 1/2. Similarly, Equation (4.21b) requires the voltage
at time step n + 1/2, but it is only available at integer time steps n. To estimate these
quantities with second-order accuracy, it is common to use two point averaging in time,
i.e.,
n+1/2 n−1/2
i+1/2 + i+1/2
i+1/2
n
(4.22a)
2
n+1/2 in+1 + in
i . (4.22b)
2
Substituting (4.22) in (4.21) and manipulating, we find the update equations for lossy
media:
Interleaved leapfrog for lossy line

n+1/2 2L − R t n−1/2 2t n
i+1/2 = i+1/2 − i+1 − in
2L + R t (2L + R t) x
(4.23)
2C − G t 2t n+1/2 n+1/2
in+1 = in − i+1/2 − i−1/2
2C + G t (2C + G t) x
The corresponding FDTD algorithms for 1D TE and 1D TM electromagnetic waves

in lossy media can be obtained from Equation (4.23) by making the appropriate substi-
tutions. In particular, one would substitute , μ, σ, σm for C, L , G, R, respectively; and
either y , z (TE Mode) or z , y (TM Mode) in place of and .
It should be apparent, and one can easily show by substitution, that when the resistance
and conductance R and G above are set to zero, the update equations in (4.23) collapse
to the same equations for lossless media, given in Equations (3.86).
Figure 4.10 shows an example of a lossy 1D simulation, using the same code as in
Figure 4.4. The results of Figure 4.4 are repeated in the top panel for comparison; below,
we have introduced a conductivity σ = 0.001 S/m into the FDTD algorithm.
Recall from electromagnetic theory that in general the decay rate of a wave is given
by [3]:
& & 1/2
μ σ 2
α=ω 1+ −1 nepers/m
2 ω
In our example in Figure 4.10, our Gaussian pulse has a width of about 0.6 m, for
which we can approximate a frequency of ∼500 MHz. Substituting this frequency into
the equation above, we find a decay rate of ∼ α = 0.189 np/m. The smallest pulse
at the far right of Figure 4.10 has traveled 9.2 m, so its theoretical amplitude should
be e−0.189·9.2 = 0.176 V/m — in fact, measurement shows that the amplitude in the
Lossless, σ = 0
25.2 ns 10.8 ns t = 3.6 ns 18.0 ns t = 32.4 ns

1
z field (V/m)
0.5
Lossy, σ = 0.001 S/m
1
t = 3.6 ns
z field (V/m)
10.8 ns
0.5 18.0 ns
25.2 ns
t = 32.4 ns
0
−10 −5 0 5 10
x (m)
Figure 4.10 Lossy 1D FDTD simulation example. The case with σ = 0.001 S/m and σm = 0 is
compared to the lossless case above it.
simulation shown is ∼ 0.178 V/m, so our lossy FDTD method has very accurately
simulated the decay of this pulse.
4.4.2 2D and 3D waves in lossy media

The procedure for incorporating the loss terms into the 2D FDTD algorithms follows
along very similar lines. For example, the corresponding versions of (4.15) (z compo-
nent) and (4.17) (z component) for lossy media are:
n+1/2
2μ − σm t n−1/2
z = z (4.24a)
i+1/2, j+1/2,k 2μ + σm t i+1/2, j+1/2,k
⎡ n n
x
−x
2t ⎢ i+1/2, j+1,k i+1/2, j,k
+ ⎣
2μ + σm t y
n n ⎤

y − y
i+1, j+1/2,k i, j+1/2,k ⎥
− ⎦.
x
4.4 FDTD algorithm for lossy media 91
Table 4.2 Coefficients for Lossy Media
Update Equation Lossless C1 Lossless C2 Lossy C1 Lossy C2

t 2 − σ t 2t
x,y,z (4.18) 1
2 + σ t 2 + σ t
t 2μ − σm t 2t
x,y,z (4.19) 1
μ 2μ + σm t 2μ + σm t
n+1
2 − σ t n
z = z (4.24b)
i, j,k+1/2 2 + σ t i, j,k+1/2
⎡ n+1/2 n+1/2
y
− y
2t ⎢ i+1/2, j,k+1/2 i−1/2, j,k+1/2
+ ⎣
2 + σ t x
n+1/2 n+1/2 ⎤

x − x
i, j+1/2,k+1/2 i, j−1/2,k+1/2 ⎥
− ⎦.
y
The similarity of the multiplying terms with the 1D version given in Equation (4.23)
is evident by making the substitutions ↔ C, μ ↔ L, σ ↔ G, and σm ↔ R.
From the 1D lossy Equations (4.23) and the two component equations above of the
3D algorithm, one should be able to distinguish a pattern in the coefficients of these
equations. If we write any particular component equation of the 3D algorithm in the
following form:
n+1 = C1 n + C2 ∇FD n+1/2
where we have used the notation ∇FD to denote the spatial differences on the right-hand
sides of Equations (4.18) and (4.19), noting that this includes the x, y, or z in the
denominator. The fields and refer to either or . We can then write the coefficients
C1 and C2 for each equation, as in Table 4.2. From these coefficients one can transform
from the lossless 3D update Equations (4.18) and (4.19) to their lossy versions; in the
same way one can find the lossy update equations for the 1D and 2D FDTD algorithms
for either TE or TM modes.
Figure 4.11 shows a 1D slice through the 3D simulation of Figure 4.9 along the x-axis.
The lossless case of Figure 4.9 is shown, as well as a lossy simulation in which we have
included σ = 10−4 S/m. One can easily observe that the field amplitudes decay with
distance faster than the lossless case, which approximately follows the 1/r trend.3
3 This trend line does not fit perfectly because we have not run the simulation long enough to smooth out any
start-up transients.
0.1
1/r trend line
0.05
Field (V/m)
−0.05
lossy σ = 0.0001 S/m
lossless σ = 0
−0.1
−50 −25 0 25 50
x (m)
Figure 4.11 Lossy 3D FDTD simulation example. The same simulation as Figure 4.9 is used, but
this time with σ = 10−4 S/m. This figure shows a 1D slice through the 3D simulation.
4.5 Divergence-free nature of the FDTD algorithm
Before proceeding further in our discussion, we must point out an important feature of
the FDTD algorithm as presented above. The integral and differential forms of Maxwell’s
curl equations (Faraday’s law and Ampère’s law, Equations 2.1 and 2.3) have been shown
to contain Gauss’s laws (Equations 2.2 and 2.4) implicitly; this implies that the fields are
divergence-free. However, after discretizing these equations and creating the difference
equations above, there is no guarantee that the divergence-free nature of Maxwell’s
equations is maintained. The conservation of the divergence-free nature of the fields is an
important property of the FDTD algorithm and the Yee grid, as we now demonstrate here.
Consider the integral form of Equation (2.2) in a charge-free region (i.e., ρ̃ = 0):

∇ ·=0 → · ds = 0. (4.25)
S
We apply this integral to the closed rectangular box surface (with sides x, y, and
z) shown in Figure 4.12. Noting that ds is by definition always outward from the
surface, we have

· ds = − x yz + x yz
S i+1/2, j i+3/2, j

− y xz + y xz . (4.26)
i+1, j−1/2 i+1, j+1/2
Taking the derivative of (4.26) with respect to time, and applying (4.15) to the ∂x /∂t
and ∂ y /∂t terms, we find:

∂ Bz − Cz Az − Dz
· ds = − yz + yz (4.27)
∂t S y y

Cz − Dz Bz − Az
− xz + xz = 0,
x x
where, for example, Bz simply denotes the z component

of the magnetic field at the
point marked B in Figure 4.12, as a shorthand for z i+1/2, j+1/2 . This result indicates
4.6 The FDTD method in other coordinate systems 93
y y
z x x
z
x
(i, j+1)
z B y
z z A
y B A z
y Δy
x x
x x
(i, j)
C y
z z z D
C D z
y Δx Δz
Figure 4.12 FDTD evaluation of the divergence of = . By integrating Gauss’s law over a
unit cell enclosed by ABCD, we demonstrate that the divergence-free nature of Maxwell’s
equations holds for the FDTD algorithm.
that if the initial conditions were such that the field was divergence free, then
' Gauss’s law
is automatically satisfied throughout any FDTD calculation, since ∂/∂t[ S · ds] = 0,
so that ∇ · = 0.
It is worth noting, however, that this divergence-free nature of the FDTD algorithm
is only valid in source-free media and for our second-order centered differencing. If we
were to use a higher-order method for spatial discretization and then apply the same
algorithm as above, in general we would find that the method is not divergence-free.
4.6 The FDTD method in other coordinate systems
It goes without saying that the FDTD method can easily be implemented in other
orthogonal curvilinear coordinate systems, such as cylindrical or spherical coordinates.
For this purpose, we can imagine a Yee cell drawn in such a coordinate system, although it
might well be a challenge to accurately draw one, and basically follow the same procedure
as before. This is particularly easy for the 2D systems with cylindrical symmetry, since
decoupled TE and TM modes can be identified in the same manner as was done above
for Cartesian coordinates. In this section we will derive the update equations for 2D
cylindrical coordinates, as well as for 3D cylindrical and spherical coordinates.
4.6.1 2D polar coordinates

In this section we will examine the FDTD method for 2D coordinate systems based
on the cylindrical coordinates r , φ, and z. In 2D, we can implement a polar coordinate
system (r − φ), in which the system is symmetric along the z-axis; or a cylindrical
system (r − z), in which the system is azimuthally symmetric. To illustrate these two
methods, we begin with Maxwell’s equations in 3D cylindrical coordinates, given by:

∂r 1 1 ∂z ∂φ ∂r 1 ∂φ 1 ∂z
= − = −
∂t r ∂φ ∂z ∂t μ ∂z r ∂φ

∂φ 1 ∂r ∂z ∂φ 1 ∂z ∂r
= − = − . (4.28)
∂t ∂z ∂r ∂t μ ∂r ∂z

∂z 11 ∂ ∂r ∂z 1 1 ∂r ∂
= (r φ ) − = − (r φ )
∂t r ∂r ∂φ ∂t μ r ∂φ ∂r
First, let us consider the case where there is no variation in the z direction, i.e., the
two-dimensional polar coordinate case. This r − φ coordinate system is useful when
simulating structures that lend themselves to circular cross-sections, such as circular
waveguides or coaxial cables [4, 5]. Setting all derivatives with respect to z to zero and
grouping the equations, we have, just as in the Cartesian system, two modes; the TE
mode:

∂r 1 1 ∂z
= (4.29a)
∂t r ∂φ

∂φ 1 ∂z
TE mode,∂/∂z = 0 =− (4.29b)
∂t ∂r

∂z 1 1 ∂r 1 ∂(r φ )
= − (4.29c)
∂t μ r ∂φ r ∂r
and the TM mode:

∂r 1 1 ∂z
=− (4.30a)
∂t μ r ∂φ

∂φ 1 ∂z
TM mode,∂/∂z = 0 = (4.30b)
∂t μ ∂r

∂z 1 1 ∂(r φ ) 1 ∂r
= − . (4.30c)
∂t r ∂r r ∂φ
As before, the TE designation implies that the singular component (z in this case)
in normal to the plane of the 2D coordinate system; or rather, that the components
are in the plane. Using the unit cells shown in Figure 4.13, we can discretize Equations
(4.29) and (4.30) in a straightforward manner, following a procedure exactly similar to
the Cartesian case, except for the (1/r ) factors. While these might seem at first glance
to be a potential difficulty, it turns out that the (1/r ) terms are quite simple: one merely
takes the value of r at the current location in the grid wherever r appears in the equations.
Similarly, the partial derivatives that include r terms, such as ∂(r φ )/∂r in Equation
y y
TM mode TE mode
z
ϕ ϕ r
r
z
(i, j+1) (i, j+1) z z
ϕ (i+1, j) (i+1, j)
r z
r ϕ
(i, j) (i, j)
z x z x
Figure 4.13 Two-dimensional FDTD cell in polar coordinates. As in Cartesian coordinates, fields
are located using the pattern in Table 4.1; for example, φ is located at an integer location in φ,
but at half-integer locations in r and z.
(4.30c), must be treated as if the “field” r φ is to be discretized; hence, our second-order

approximation of this partial derivative is:

∂(r φ ) ri+1 φ i+1, j − ri, j φ i, j
+ O(r 2 ). (4.31)
∂r r
i+1/2, j
Thus, including these partial derivatives, we find the update equations for the TE
mode by discretizing Equations (4.29):
⎡ n+1/2 n+1/2 ⎤

n+1 n z −z
t ⎢ i+1/2, j+1/2 i+1/2, j−1/2 ⎥
r = r + ⎣ ⎦ (4.32a)
i+1/2, j i+1/2, j i+1/2, j ri+1/2 φ
⎡ n+1/2 n+1/2 ⎤

n+1 n z −z
t ⎢ i+1/2, j+1/2 i−1/2, j+1/2 ⎥
φ = φ − ⎣ ⎦ (4.32b)
i, j+1/2 i, j+1/2 i, j+1/2 r
n+1/2 n−1/2

z = z (4.32c)
i+1/2, j+1/2 i+1/2, j+1/2
⎡ n n

r −r
t ⎢ i+1/2, j+1 i+1/2, j
+ ⎣
μi+1/2, j+1/2 ri+1/2 φ
n n ⎤

ri+1 φ −ri φ
i+1, j+1/2 i, j+1/2 ⎥
− ⎦
ri+1/2 r
0 3
Transformation to
60 polar coordinates
120
Angle (deg)
y (meters)
180 0
240
−1.5
300
360 −3
−3 −1.5 0 1.5 3
x (km)
0.2 r−1/2 trend line
z (V/m)
−0.2 z |z|
0 1 2 3
r (km)
Figure 4.14 2D FDTD simulation in polar coordinates, TM mode. The calculations yield fields in
a 2D r − φ grid, as shown in the top left; these are then transformed
√ to x − y to display in a
polar plot at right. The 1D slice shows how the fields decay as 1/ r .
where ri = ir and ri+1/2 = (i + 1/2)r , and i and j are respectively the indices in
the r and φ coordinates; later, in 3D cylindrical coordinates, we will use k as the index
in the z direction as usual. Similarly, for the TM mode we find:
⎡ n n ⎤

n+1/2 n−1/2 z −z
t ⎢ i, j+1 i, j ⎥
r = r − ⎣ ⎦ (4.33a)
i, j+1/2 i, j+1/2 μi, j+1/2 ri φ
⎡ n n ⎤

n+1/2 n−1/2 z −z
t ⎢ i, j i+1, j ⎥
φ = φ + ⎣ ⎦ (4.33b)
i+1/2, j i+1/2, j μi+1/2, j r
n+1 n

z = z (4.33c)
i, j i, j
⎡ n+1/2 n+1/2 n+1/2 n+1/2 ⎤

ri+1/2 φ −ri−1/2 φ r −r
t ⎢ i+1/2, j i−1/2, j i, j+1/2 i, j−1/2 ⎥
+ ⎣ − ⎦.
i, j ri r ri φ
Example polar simulation

Figure 4.14 shows an example simulation in 2D polar coordinates for the TM mode, in
which we have directly discretized Equations (4.33). For this example, we excite an z
source at the ‘center’ of a circular “grid” with center frequency f 0 = 1 MHz. We choose
r = 15 m, so that we have 20 cells per wavelength in r , and we break the azimuthal grid
into 200 cells, so that φ = 1.8 degrees or 0.0314 radians. The r , φ , and z fields
are stored in 2D arrays of 200 × 200 grid points, as they were for Cartesian coordinates;
what this means is that the spatial distance in the azimuthal direction grows with r . For
example, at the first grid point from the center where r = r , the steps in φ have length
r φ (in radians), which in our case is only 0.47 m. At the outer edge of the space
where r = 200r , the steps in φ are 94.2 m! Hence, while this grid is orthogonal, since
all of the coordinate vectors are orthogonal to each other, it is non-uniform.
Recall that for Cartesian coordinates, we found that stability was ensured in 1D when,
for a given spatial grid size x, t was chosen such that v p t/x ≤ 1. In 2D, √ this
requirement becomes more strict; if we choose x ≤ y, we require v p t/x ≤ 1/ 2;
i.e., the stability criterion is set by the smallest spatial grid size. Here, our smallest grid
cells are those closest to the center;√ with r = 15 m and r φ ≥ 0.47 m, we require
approximately v p t/r φ ≤ 1/ 2. The strict 2D stability criterion will be derived in
detail in Chapter 5.
The resulting z field in r − φ coordinates is shown in Figure 4.14, showing a snapshot
at time t = 8 μs. By transforming the coordinates into x − y pairs, we plot these fields
in their correct polar pattern on the right. As one can imagine by comparing the left
(uniform) grid with the right (nonuniform) grid, there is a much higher density of points
at the center of the grid compared to the edges! As such, this scheme is a simple example
of a nonuniform grid. The lower panel shows the 1D slice through a radius of the field
pattern, showing the r −1/2 falloff of the fields with distance as expected.
Finally, note the discontinuity between the first and last array indices in φ. In our polar
coordinate system in Figure 4.13, for the TM mode, the z component on the x-axis (i.e.,
at φ = 0) requires the r component below (i.e., at φ = −φ/2) in its update equation.
This minor correction can be implemented quite easily in the update equations; similar
corrections are required for the TE mode where the r component is on the φ = 0 axis.
4.6.2 2D cylindrical coordinates

A second situation that arises in many problems is that of cylindrical symmetry, in
which there is no variation in the azimuthal (i.e., φ) direction, but fields vary with r and
z. This scenario might occur, for example, with a vertically oriented antenna source,
whose radiation pattern is azimuthally symmetric; however, the medium and scattering
objects must also be azimuthally symmetric about the same axis. In this scenario, we set
derivatives with respect to φ to zero in Equations (4.33), and again find two separable
modes:

∂r 1 ∂φ
= (4.34a)
∂t μ ∂z

∂φ 1 ∂r ∂z
TM mode,∂/∂φ = 0 = − (4.34b)
∂t ∂z ∂r

∂z 11 ∂
=− (r φ ) (4.34c)
∂t μ r ∂r
z TM mode z TE mode
(i, k+1) (i, k+1)
ϕ
ϕ z r
r r z
z
ϕ
ϕ ϕ
(i, k) z (i+1, k) (i, k) r (i+1, k)
r r
Figure 4.15 Two-dimensional FDTD cell in 2D cylindrical coordinates, with ∂/∂φ = 0. The
same pattern of field placement holds as described in Figure 4.13.
and

∂r 1 ∂φ
=− (4.35a)
∂t ∂z

∂φ 1 ∂z ∂r
TE mode,∂/∂φ = 0 = − (4.35b)
∂t μ ∂r ∂z

∂z 11 ∂
= (r φ ) . (4.35c)
∂t r ∂r
Again, this latter set of equations is designated the TE mode because the components
of are in the r − z plane of the simulation. These equations can be discretized in the
same way as the previous case, using the unit cells shown in Figure 4.15. This procedure
yields the update equations, first for the TM mode:
⎡ n n ⎤

n+1/2 n−1/2 t φ −φ
⎢ i,k+1 i,k ⎥
r = r − ⎣ ⎦ (4.36a)
i,k+1/2 i,k+1/2 μi,k+1/2 z
⎡ n+1/2 n+1/2 n+1/2 n+1/2 ⎤

n+1 n t r −r z −z
⎢ i,k+1/2 i,k−1/2 i+1/2,k i−1/2,k ⎥
φ = φ + ⎣ − ⎦ (4.36b)
i,k i,k i,k z r
⎡ n n ⎤

n+1/2 n−1/2 t ri+1 φ −ri φ
⎢ i+1,k i,k ⎥
z = z − ⎣ ⎦, (4.36c)
i+1/2,k i+1/2,k μi+1/2,k ri+1/2 r
and for the TE mode:

⎡ n+1/2 n+1/2 ⎤

n+1 n φ −φ
t ⎢ i+1/2,k+1/2 i+1/2,k−1/2 ⎥
r = r + ⎣ ⎦ (4.37a)
i+1/2,k i+1/2,k i+1/2,k z
n+1/2 n−1/2

φ = φ (4.37b)
i+1/2,k+1/2 i+1/2,k+1/2
⎡ n n n n ⎤

z −z r −r
t ⎢ i+1,k+1/2 i,k+1/2 i+1/2,k+1 i+1/2,k ⎥
+ ⎣ − ⎦
μi+1/2,k+1/2 r z
⎡ n+1/2 n+1/2 ⎤

n+1 n ri+1/2 φ −ri−1/2 φ
t ⎢ i+1/2,k+1/2 i−1/2,k+1/2 ⎥
z = z + ⎣ ⎦.
i,k+1/2 i,k+1/2 i,k+1/2 ri r
(4.37c)
A problem with this algorithm arises when one considers the on-axis z (in the TE
mode) or z (in the TM mode) fields, i.e. at r = 0. Clearly, these update equations will
not work due to the 1/r terms, leading to division by zero. A correction to the on-axis
fields, however, is very straightforward. In the TE mode, for example, consider Ampère’s
law, Equation (2.3), applied in its integral formto a small circular region enclosing the
axis at a radius r = r/2, at the location of φ i=1/2,k+1/2 , as shown in the left panel of
Figure 4.16. In the absence of any source current, we have:

∂
· dl = · ds.
C S ∂t
For our purposes here, we assume that φ does not vary along the curve C, which
is a small circle in φ at r = r/2 (i.e., our pre-assumed azimuthal symmetry); and that
z does not vary inside the surface S, which is the area enclosed by C. Since z is
azimuthally symmetric, this latter assumption implies that z also does not vary with r
between r = 0 and r = r/2. Thus we have:
n+1/2 ∂z n+1/2

φ (2πri=1/2 ) = 2
(πri=1/2 ). (4.38)
1/2,k+1/2 ∂t 0,k+1/2
Noting that ri=1/2 = r/2, we find the update equation for z at r = 0:

n+1 n n+1/2
4t
z = z + φ . (4.39)
0,k+1/2 0,k+1/2 0 r 1/2,k+1/2
While this equation appears to suggest that the on-axis field depends only on one
value of φ , the application of the integral form of Ampère’s law shows that it depends
on the value of φ at the r = r/2 grid location, but at all locations in azimuth, over
On-axis cell Regular cell

z
zk+1 r
ϕ ϕ ϕj
r ϕ z
C = 2 Δr z r ϕj+1
k+3/2 r ϕ
z
ϕ z
zk
Δr/2 z ϕ r
ϕ r
ϕ z z
Δr
ϕ r ϕ
k−1/2
S = Δr2
Δr ri ri+1
Figure 4.16 Cylindrical coordinate system used in FDTD. The full 3D cylindrical coordinate
system is shown at right. At left, the method for determining the on-axis z component for the
2D TE mode, by directly applying Ampère’s law.
which the field does not vary. We will see a similar analysis in the 3D cylindrical case
next, in which φ is not azimuthally symmetric.
Example cylindrical simulation

Figure 4.17 shows an example simulation using this TE mode cylindrical coordi-
nate algorithm. This example uses grid cells of λ/20 with λ = 3 m (a frequency of
100 MHz). We excite z at r = 0, z = 0, and set z at the next five grid cells both
above and below the source to zero; this emulates a “half-wave” dipole antenna, 10 grid
cells long, with the source excited at the center as a voltage = z,in z. Note that
attempting to implement a small dipole in this scenario will not work as one would
expect from antenna theory; the “small dipole” may be only λ/20 or 15 cm in length,
but it is also 15 cm wide! In Chapter 7 we will introduce a method to circumvent this
problem known as the thin-wire approximation.
The right panel in Figure 4.17 plots the simulated dipole radiation pattern, measured
by sampling the z field amplitude along a circle of constant radius, shown by the dashed
line, and then normalizing. The ideal half-wave dipole radiation pattern is also shown
for comparison, given by [e.g., 6, p. 61]:
cos[(π/2) cos θ ]
F(θ ) = .
sin2 θ
It is thus evident that our half-wave dipole simulation approximately emulates the
behavior of a real half-wave dipole.
15 0
Measured
/2 1.0 radiation
/3
10 pattern
0.75 (normalized)
−1
Log10 (z) Amplitude

0.5 /6
5
z (meters)
0.25
0 −2
θ=0
−5
−3
− /6
−10
− /3
−15 −4
3 6 9 12 15 − /2 Ideal half-wave dipole
r (meters) (dashed)
(a) (b)
Figure 4.17 Example FDTD simulation in cylindrical coordinates. (a) 2D z field pattern in r − z
space at a snapshot in time. (b) The resulting radiation pattern as measured along the dashed
line in the left panel. The dashed line on the right is the ideal half-wave dipole radiation pattern.
Since we have chosen a dipole radiator at the center of what is effectively a 3D

space with azimuthal symmetry, our simulation is nearly identical to the 3D case in
Figure 4.9, except that we have reduced the simulation space from 100 × 100 × 100
grid cells to 100 × 200 grid cells, and the number of field components from six to three.
Hence, for problems which can be approximated with azimuthal symmetry, this 2D
formulation is invaluable.
4.6.3 3D cylindrical coordinates

Implementation of the 3D FDTD method in cylindrical coordinates simply requires
discretization of the full 3D cylindrical Maxwell’s equations, given in Equation (4.28)
above. As an example, we will show the discretization of the z component equation,
repeated below:

∂z 11 ∂ ∂r
= (r φ ) − .
∂t r ∂r ∂φ
The grid cell for the 3D cylindrical case is shown in Figure 4.16. As in Cartesian
coordinates, the locations of field components follow the same pattern as shown in
Table 4.1, where now (i, j, k) refer to indices in (r, φ, z). Applying the same discretization
rules as for the 2D cylindrical case, we find the update equation:

n+1 n

z = z (4.40)
i, j,k+1/2 i, j,k+1/2
⎡ n+1/2 n+1/2

ri+1/2 φ −ri−1/2 φ
t ⎢ i+1/2, j,k+1/2 i−1/2, j,k+1/2
+ ⎣
ri r
n+1/2 n+1/2 ⎤

r −r
i, j+1/2,k+1/2 i, j−1/2,k+1/2 ⎥
− ⎦.
ri φ
Similar discretization procedures lead to update equations for the other components
in Equation (4.28); the remaining five update equations are given by:
n+1 n

r = r (4.41a)
i+1/2, j,k i+1/2, j,k
⎡ n+1/2 n+1/2

z − z
t ⎢ i+1/2, j+1/2,k i+1/2, j−1/2,k
+ ⎣
ri+1/2 φ
n+1/2 n+1/2 ⎤

φ −φ
i+1/2, j,k+1/2 i+1/2, j,k−1/2 ⎥
− ⎦
z
n+1 n

φ = φ (4.41b)
i, j+1/2,k i, j+1/2,k
⎡ n+1/2 n+1/2

r − r
t ⎢ i, j+1/2,k+1/2 i, j+1/2,k−1/2
+ ⎣
z
n+1/2 n+1/2 ⎤

z − z
i+1/2, j+1/2,k i−1/2, j+1/2,k ⎥
− ⎦
r
n+1/2 n−1/2

r = r (4.41c)
i, j+1/2,k+1/2 i, j+1/2,k+1/2
⎡ n n

φ − φ
t ⎢ i, j+1/2,k+1 i, j+1/2,k
+ ⎣
μ z
n n ⎤

z − z
i, j+1,k+1/2 i, j,k+1/2 ⎥
− ⎦
ri φ
n+1/2 n−1/2

φ = φ (4.41d)
i+1/2, j,k+1/2 i+1/2, j,k+1/2
⎡ n n

z − z
t ⎢ i+1, j,k+1/2 i, j,k+1/2
+ ⎣
μ r
n n ⎤

r − r
i+1/2, j,k+1 i+1/2, j,k ⎥
− ⎦
z
n+1/2 n−1/2

z = z (4.41e)
i+1/2, j+1/2,k i+1/2, j+1/2,k
⎡ n n

r − r
t ⎢ i+1/2, j+1,k i+1/2, j,k
+ ⎣
μ ri+1/2 φ
n n ⎤

ri+1 φ − r i φ
i+1, j+1/2,k i, j+1/2,k ⎥
− ⎦.
ri+1/2 r
Note that we have assumed a constant here, but an inhomogeneous medium can
easily be implemented as we have shown earlier, by using the values of and μ at the
grid location at which the equation is discretized (i, j, k + 1/2 in Equation 4.40).
As we can see from Figure 4.16, however, there is a problem on the axis, as there was
in the 2D cylindrical case. In 3D, if we choose the axis to align with z , then we have
z , r , and φ components on the axis to deal with, as shown in Figure 4.16.
To update z on the axis, we apply Ampère’s law directly, similar to the 2D case.
In this case, however, we must integrate φ along the φ direction, at r = r/2, rather
than simply multiplying by 2π (r/2). Integrating in continous space is simply a sum in
discrete space, so we find [7]:
n+1 n 4t ( n+1/2
Nφ

z = z + φ . (4.42)
0, j,k+1/2 0, j,k+1/2 Nφ 0 r i=1/2, p,k+1/2
p=1
Note that the summation above, when divided by Nφ (the number of grid points in
φ), is actually the average value of φ around the contour; as such, Equation (4.42) is
nearly identical to Equation (4.39).
Now, despite the fact that φ and r have components on the axis, these axial
components need not be computed for the remainder of the space to be updated normally.
To see this, consider the on-axis cell in Figure 4.16. The fields that could require
knowledge of the on-axis fields are r , z , and φ at r = r/2. Now consider Equations
(4.41): in Equation (4.41a), r does not require any of the on-axis field components;
in Equation (4.41d), φ requires only z on the axis, which we have updated above in
Equation (4.42). Finally, in Equation (4.41e), z requires φ on the axis in the last term
of the equation, but this component is multiplied by ri which is equal to zero on the
axis, so it is not needed. In summary, the full set of update equations can be used for the
half-step away from the r = 0 axis, and the φ and r components on-axis are simply
ignored.
It is interesting to note that these axial r and φ components both point away from
the axis, and thus could equally be considered φ and r components. As such, the
computed φ and r fields at r = r/2 can be used to approximate these on-axis fields
if necessary.
Finally, note that in addition to the axial correction described above, the 3D cylindrical
coordinate problem has the same 2π discontinuity in φ which must be accounted for in
the update equations, as in the 2D polar coordinate case earlier.
4.6.4 3D spherical coordinates

Discretization to find the FDTD update equations in spherical coordinates is somewhat
more complicated than cylindrical coordinates, but follows the same principles. The
six component Maxwell’s equations in spherical coordinates are given below for simple
media [e.g., 3]:

∂r 1 ∂ ∂θ
+ σ r = (φ sin θ ) − (4.43a)
∂t r sin θ ∂θ ∂φ

∂θ 1 1 ∂r ∂
+ σ θ = − (r φ ) (4.43b)
∂t r sin θ ∂φ ∂r

∂φ 1 ∂ ∂r
+ σ φ = (r θ ) − (4.43c)
∂t r ∂r ∂θ

∂r 1 ∂ ∂θ
μ =− (φ sin θ ) − (4.43d)
∂t r sin θ ∂θ ∂φ

∂θ 1 1 ∂r ∂
μ =− − (r φ ) (4.43e)
∂t r sin θ ∂φ ∂r

∂φ 1 ∂ ∂r
μ =− (r θ ) − . (4.43f )
∂t r ∂r ∂θ
There certainly exist circumstances where a 2D spherical coordinate system is useful;
for example, the simulation shown in Figure 4.17 has azimuthal symmetry (∂/∂φ =
0), but would be more accurately described in r − θ space rather than r − z space.
The method for creating a 2D simulation follows the method described above for the
cylindrical coordinate case: set ∂/∂φ = 0 in the equations above, organize the resulting
six equations into TE and TM modes, and then discretize as we have done throughout
this chapter. In the interest of brevity, we will skip the 2D spherical coordinate case and
jump directly into the 3D discretization process.
The derivation of the FDTD algorithm in spherical coordinates was first presented
by R. Holland in 1983 [8]. Figure 4.18 shows the spherical coordinate system and the
unit cell for FDTD. Noting that in spherical coordinates, the (i, j, k) subscripts refer
to (r, θ, φ), respectively, the field vectors are all drawn pointing toward the increasing
ϕ ri+1,θj,ϕk
ri+1,θj,ϕk+1
ri+1
r
θj ri r θ
θ
ri,θj,ϕk+1
θj+1 ϕ
ϕk ϕ ri+1,θj+1,ϕ
ϕk+1 θ
r
θ r
ri,θj+1,ϕk
ri,θj+1,ϕk+1 ϕ
Figure 4.18 Spherical coordinate system used in FDTD. The left panel shows how the coordinate
system is drawn; the right panel shows an expanded view of a unit cell. Adapted from [8].
coordinate direction. Then, the locations of field components follow the same logic
from Table 4.1; i.e., the n components fall on the half-integer grid location for n
and the integer grid location for the other two coordinates, and vice versa for the
field components. For example, r is located at (i + 1/2, j, k) and r is located at
(i, j + 1/2, k + 1/2). Note, however, that our cell size is not strictly r × φ × θ ;
those last two do not even have spatial units. Rather, the grid cell has dimensions r ,
r θ , and r sin θ φ. Thus, if we use constant r , θ , and φ, the grid cells grow with
increasing r , and, to no surprise, we have another nonuniform grid.
As in the cylindrical coordinate case, we deal with the r and θ terms on the right
side of Equations (4.43) by simply using the values of r and θ at the location that the
equation is discretized. Similarly, the partial derivatives that include r or sin θ must be
treated in the same manner as Equation (4.31) earlier.
Applying the methods from this chapter to Equation (4.43a), and assuming σ = 0 for
simplicity, we find its update equation:
n+1 n

r = r
i+1/2, j,k i+1/2, j,k
⎡ n+1/2 n+1/2

φ sin θ j+1/2 − φ sin θ j−1/2
t ⎢ i+1/2, j+1/2,k i+1/2, j−1/2,k
+ ⎣
ri+1/2 sin θ j θ
n+1/2 n+1/2 ⎤

θ −θ
i+1/2, j,k+1/2 i+1/2, j,k−1/2 ⎥
− ⎦.
φ
(4.44)
Notice that this equation requires the value of r at a half-step, coincident in the i index
with the r field; and it also requires the values of θ at j and j ± 1/2. Of course, these
values are readily available and are simply defined by the geometry of the simulation
space. For a regular grid, ri = ir and θ j = jθ . The other update equations follow
from similar discretization of Equations (4.43b–f).
Now, recall that in polar coordinates, we needed to take care of the 2π discontinuity
in φ; the same issue exists in spherical coordinates, and must be accounted for in the
same manner. In addition, recall that in the 2D polar and 3D cylindrical cases, special
attention was given to field values at r = 0; this point is known as a singularity, since
some component of the update equations has a zero in the denominator. In spherical
coordinates, we have the same singularity at r = 0, which again requires special attention,
which we shall deal with in a moment; in addition, however, there are problems at θ = 0
and θ = π , i.e., along the z-axis. Observe that in Equation (4.44), the term sin θ appears
in the denominator; since sin(0) = sin(π ) = 0, we again have a singularity. As one might
anticipate, this singularity also exists in the discretized equivalents of Equations (4.43b,
d, and e), where sin θ appears in the denominator.
Holland’s [8] original formulation described a fix to the θ = 0 ( j = 1) and θ = π
( j = Nθ ) problem. Note that the only field component that presents a problem on this
axis is the r component. The solution is to apply the same Ampère’s law contour
integral around a small contour C enclosing the z-axis, either near the north pole or the
south pole of the sphere in Figure 4.18, and evaluating the integral in φ as a summation,4
we find the following corrected update equations at θ = 0 and θ = π :
n+1 n t sin θ1/2 (Nφ n+1/2

r = r + φ φ (4.45a)
i+1/2,1,k i+1/2,1,k 2π ri+1/2 (1 − cos θ1/2 ) p=1 i+1/2,1, p
n+1 n t sin θ Nθ −1/2 (Nφ n+1/2

r = r + φ φ
i+1/2,Nθ ,k i+1/2,Nθ ,k 2π ri+1/2 (1 + cos θ Nθ −1/2 ) p=1 i+1/2,Nθ , p
(4.45b)
where we note that θ1/2 = θ/2 and θ Nθ −1/2 = π − θ/2.
Now, dealing with the origin (r = 0 or i = 0) is somewhat more complicated. If one
imagines the unit cell in Figure 4.18 collapsed down to the origin, we end up with a
small trapezoidal grid cell with an apex at the origin. Depending on where we put the
origin, we may either have the triplet (r , θ , φ ) defined at the origin, or taking it a
half-step further out in r , we have the complementary triplet (r , θ , φ ).
4 Note in this case that the area of integration on the right-hand side of Ampère’s law is the small “cap” on
the top of a sphere, and not a simple circle. The area of a small cap on a sphere at radius ri and angle θ j is:
2π θj
ri2 sin θ dθ dφ = 2πri2 (1 − cos θ j ).
0 0
A similar integral for the south pole cap of the sphere yields an analogous expression, except with a
(1 + cos θ j ) term. However, if we observe that cos(π − θ j ) = −cos θ j , we find that the denominators in
Equations (4.45) are the same.
4.7 Summary 107
In the rare cases where the origin is not part of the source radiator, the methods
used to deal with this singularity involve solving a small region around the origin in
Cartesian coordinates, and then mapping those fields back to spherical coordinates on
the boundary of the small cube. The analysis of this method, however, is beyond the
scope of our current discussion. We refer interested readers to the papers by Holland [8]
and G. Liu et al. [9].
4.7 Summary
In this chapter, the interleaved leapfrog method from Chapter 3 was used to derive the
FDTD algorithm in 1D, 2D, and 3D. In 1D, the FDTD algorithm for the TM mode is
given by:
n+1/2 n−1/2 n
t n
y = y + z − z
i+1/2, j i+1/2, j μi+1/2, j x i+1, j i, j
n+1 n (1D TM mode)
t n+1/2 n+1/2
z = z + y − y ,
i, j i, j i, j x i+1/2, j i−1/2, j
with similar equations for the TE mode. In 2D, the FDTD algorithm for the TM mode
is given by:
n+1/2 n−1/2 n
t n
x = x − z − z (2D TM mode)
i, j+1/2 i, j+1/2 μi, j+1/2 y i, j+1 i, j
n+1/2 n−1/2 n
t n
y = y + z − z
i+1/2, j i+1/2, j μi+1/2, j x i+1, j i, j
⎡ n+1/2 n+1/2 n+1/2 n+1/2 ⎤

n+1 n y − y x − x
t ⎢ i+1/2, j i−1/2, j i, j+1/2 i, j−1/2 ⎥
z = z + ⎣ − ⎦,
i, j i, j i, j x y
with similar equations for the TE mode. In 3D, there are no TE versus TM modes; the
FDTD algorithm involves six update equations, given in Equations (4.18) and (4.19).
Update equations for lossy media were derived in Section 4.4; it was shown that the
update equations for lossy media can be transcribed from the lossless equations, by
simply replacing the coefficients C1 and C2 as given in Table 4.2.
Next, it was shown in Section 4.5 that the FDTD algorithm is divergence-free, i.e.,
it inherently satisfies Gauss’s law that ∇ · = 0. This is an important and appealing
feature of the FDTD algorithm; zero divergence is not a feature that is inherent to any
particular discretization scheme.
In Section 4.6 we introduced other orthogonal coordinate systems, namely cylindrical
and spherical coordinates, and the discretization process that leads to FDTD update
equations in these coordinate systems. We first introduced 2D polar and cylindrical
coordinates, which are both commonly used where there is longitudinal or azimuthal
symmetry in a problem. The discretization of Maxwell’s equations in these coordinate
systems again leads to TE and TM modes, but terms involving r or 1/r need to be
handled carefully.
Similarly, the 3D cylindrical and spherical coordinate systems require special handling
of r and sin θ terms which appear in the update equations. In all cases, these terms lead
to singularities at the origin or axis (where r = 0) and, in spherical coordinates, on the
axis where θ = 0 (i.e., the z-axis). These singularities are handled in the equations above
by directly applying Ampère’s law to a small contour around the z-axis.
4.8 Problems
4.1. A lossy high-speed chip-to-chip interconnect. Consider applying the inter-

leaved leapfrog algorithm to solve for the voltage variation on a high-speed lossy
interconnect, treating it as a lossy transmission line with G = 0, and R = 150 -
(cm)−1 . Determine and plot the voltage distribution on the lossy line at t = 0.5
ns. Select appropriate values of x and t. Assume other line parameters and
initial voltage distribution to be as in Problem 3.4 and use the same absorbing
boundary condition specified therein. Comment on your results. Do you think
this particular microstrip line is a good one to use for chip-to-chip transmission
of 0.1 ns pulses over a 5-cm long interconnect on a circuit board? Repeat for
R = 50 -(cm)−1 and R = 10 -(cm)−1 ; comment on the results.
Note: You are given 0 , but no initial condition for . You will need to use
another method to determine 1/2 , the initial condition for the current, over the
entire transmission line. This initial condition for can have a considerable effect
on the results.
4.2. High-speed GaAs coplanar strip interconnect. The per-unit-length parameters
of a high-speed GaAs coplanar strip interconnect operating at 10 GHz are mea-
sured to be R 25 -(cm)−1 , C 0.902 pF-(cm)−1 , L 9.55 nH-(cm)−1 , and
G = 0. We shall study the propagation and attenuation on this line of a Gaussian
pulse given by
= e−(n−3τ ) /τ 2
2
i=1
n
where τ t = 0.1 ns.

(a) Use an appropriate algorithm and write a program to propagate the pulse on
this line and measure the attenuation rate (in dB-m−1 ). For time-harmonic
signals, it is well known that the attenuation rate on a lossy transmission
√
line is given by α =
e{ (R + jωL)( jωC)} (in np-m−1 ), where ω is the
frequency of operation. Compare the attenuation rate you measured with the
time-harmonic attenuation rate (for the frequency range of your Gaussian
pulse) and comment on your results.
(b) Assuming the coplanar strip interconnect to be 10-cm long, terminate it
with a short circuit and run your simulation until after the Gaussian pulse
is reflected and just before its front arrives back at the source point. Plot
4.8 Problems 109
the variation with time of the propagating pulse as observed at fixed points
on your grid, i.e., at i = i max /4, i = i max /2, and i = 3i max /4, and compare
its shape (by plotting them on top of one another) with that of the original
Gaussian pulse. Note that the pulse will be observed twice at each of these
locations, once before reflection and then again after reflection. Is the shape
of the Gaussian pulse changing as it propagates? If so, why? Comment on
your results.
(c) Repeat parts (a) and (b), but excite the simulation with the derivative of the
Gaussian above. What effect does this have on results?
4.3. 1D wave incident on a dielectric slab. Consider the propagation of a one-
dimensional electromagnetic wave in a medium and its reflection from and trans-
mission through a dielectric slab. Construct a 1D FDTD space with a lossless
dielectric slab (r 2 = 2, μr 2 = 1) at the center of your grid, free space (r 1 = 1,
μr 1 = 1) to the left of the slab, and another dielectric (r 3 = 4, μr 3 = 1) to the
right of the slab. You should allow for the length of the slab to be adjustable.
Launch a Gaussian pulse from the left boundary of your FDTD space, which is
given by
n

z = e−(n−3τ ) /τ
2 2
i=1
√
where τ t = 0.8 ns. Choose x ≥ λ0 /[20 r ], where λ0 = 0.3 m, t =
x/v p , and i max x = 3 m. Once the tail end of the pulse leaves the source
(i.e., when n = 6τ ), replace the source end with a simple boundary condition:
n+1 n v t − x n+1 n
p
z = z + z −z .
1 2 v p t + x 2 1
Take the rightmost end of your FDTD space to be as far away as needed so that
there are no reflections from that boundary, or implement a similar radiation
boundary at i = i max .
√
(a) Taking a slab width of d = λ0 /[2 r 2 ], write a program to propagate the
wave on your FDTD space and observe the behavior of the pulse as it
interacts with the slab. Can you measure the “reflection coefficient”? What
is its value? What is the transmission coefficient? Plot the reflection and
transmission coefficients as a function of frequency and compare to the
analytical expression. Comment on the differences observed. Considering
the finite precision of floating-point numbers on a computer, what is the
effective “numerical bandwidth” of your pulse?
√
(b) Repeat (a) for a slab width of d = λ0 /[4 r 2 ]. Comment on your results.
4.4. 2D FDTD solution of Maxwell’s equations. Write a program implementing the
FDTD algorithm for the 2D TM mode. Assume square unit cells (i.e., x =√ y)
and free space everywhere in the grid, and use a time step of t = x/(v p 2).
Assume the outer boundary of your grid to be surrounded by perfect electrically
conducting sheets and terminate your mesh in electric field components set to
be zero at all times. Excite a radially outgoing wave in the grid by implementing
a hard source for a single electric field component at the center of the grid. For
both of the sources specified below, display the electric and magnetic fields of
the wave distribution across the grid at a number of time snapshots before and
after the wave reaches the outermost grid boundary.
(a) First consider a Gaussian source waveform given as
n

= e−(nt−1.8 ns) /(0.6 ns) .
2 2
z
imax/2,jmax/2
Use a 50 × 50 grid with x = y = 0.1λmin , where λmin is the minimum

wavelength represented by the Gaussian source. How would you determine
λmin ?
(b) Now consider a sinusoidal source waveform which is instantly turned on at
time t = 0, having a frequency of f 0 = 30 GHz. Let x = y = λ0 /10,
and once again use a 50 × 50 grid. Why is the amplitude of the outgoing
wave not constant?
(c) Repeat part (a) for the TE mode. To emulate PEC at the boundaries, note
that the tangential field must be zero, and that this will be either x or y
depending on the boundary in question.
Note that the two codes developed in this problem, for the TM and TE modes,
will be useful in later problems in this and later chapters. The student will find
it useful to write the code including parameters σ and σm , so that losses can be
easily included later.
4.5. 3D rectangular resonator. Write a 3D FDTD algorithm using Equations (4.18)
and (4.19) in a 100 × 50 × 20 m box with perfectly conducting walls. Use
x = y = z = 1 m.
(a) Excite this resonator with a sinusoidal z line source at the bottom center
of the space with a frequency of 1.5 MHz. Make this source 10 grid cells
long in z. Plot the magnitude of the field in the plane perpendicular to
the source excitation, and in a plane parallel to the excitation. What do you
notice about the field evolution with time?
(b) Repeat part (a), but at frequencies of 3 MHz and 10 MHz. Compare to
part (a).
(c) This time, excite the box with the derivative of a Gaussian source that has a
pulsewidth τ equal to five time steps. Add the source to the update equations
as a current source so that it does not interfere with reflected waves. Find
the resonance frequencies for this resonator and compare to the theoretical
resonances, which for the TEmnp and TMmnp modes are given by [3]:
&
1 m2π 2 n2π 2 p2 π 2
ωmnp = √ + +
μ a2 b2 d2
for a resonator with dimensions a, b, and d.
4.6. Antenna modeling. In this problem, we will use the 3D code from the previous
problem to model two types of antennas. For each of these examples, measure
4.8 Problems 111
the radiation pattern by monitoring points at a constant radius from the center
of the antenna. Make the simulation space 100 × 100 × 100 grid cells, and use
20 grid cells per wavelength; note that this implies the simulation boundaries are
only 2.5 wavelengths from the antenna, so the results will be near-field radiation
patterns.
(a) Model a half-wave dipole, similar to Figure 4.9, by making five grid cells
perfect conductors above and below the source point. Drive the antenna
with a sinusoidal source at the appropriate frequency.
(b) Model a monopole antenna above a ground plane by placing the source just
above the bottom of the simulation space, and five grid cells of PEC above
the source. Compare the radiation pattern to that found in (a).
In both cases above, watch out for reflections from the simulation boundaries:
you will need to compute the radiation pattern before reflections corrupt the
solution. Why does this problem not work with a small dipole, one grid cell in
length? We will revisit and solve this problem in Chapter 7.
4.7. 2D polar coordinates. Write a 2D TM mode simulation in polar coordinates that
emulates the simulation space in Problem 4.4 above. Excite the simulation space
with a sinusoidal z source at 30 GHz; compare the outgoing waves to those in
Problem 4.4.
4.8. 2D cylindrical coordinates. Create a 2D cylindrical simulation, choosing appro-
priate parameters r , z, and t to model the same wavelength as in Problem
4.6 above. Drive the simulation with the same two sources as in Problem 4.6,
and compare the radiation patterns. How do the results differ from the previous
problem? How does the simulation time compare?
4.9. 3D circular waveguide in cylindrical coordinates. Create a simulation space
for a circular waveguide in 3D cylindrical coordinates. Let the radius of the wave-
guide be 2.7 cm, and design it for a frequency of 10 GHz. Excite the waveguide at
one end with a hard source at a frequency of 10 GHz, where the source field is
defined along a line through the cross-section of the waveguide. Plot a snapshot
of the field patterns a few wavelengths down the waveguide; what modes exist in
this waveguide? Repeat for frequencies of 8 GHz and 12 GHz and compare the
results.
4.10. Optical fiber. Modify the code from the previous problem so that the wave-
guide consists of a dielectric with n = 2 with radius 2.0 cm, and a region of
free space outside the dielectric. Terminate the simulation space outside the
free space region with PEC, at a radius large enough so that evanescent fields
outside the dielectric are small. Excite this fiber at one end as in the problem
above, and plot snapshots of the field pattern a few wavelengths down the fiber.
4.11. Cylindrical resonator. Modify your 3D cylindrical coordinate code from the
previous problems to model a cylindrical waveguide with radius 2 cm and thick-
ness 2 cm (the z-dimension). Give the resonator a metallic outer shell, defined
by a lossy material with σ = 104 S/m a few grid cells thick. Excite the res-
onator at 10 GHz with an z source at the center of the resonator, spanning the
z-dimension, similar to Problem 4.5 above. The quality factor of a resonator is
defined by the ratio of the energy stored in the resonant cavity to the energy dis-
sipated (per cycle) in the lossy walls. Can you find a way to measure the quality
factor Q of this resonator?
References
[1] K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations
in isotropic media,” IEEE Trans. on Ant. and Prop., vol. 14, no. 3, pp. 302–207, 1966.
[2] A. Taflove and S. Hagness, Computational Electrodynamics: The Finite-Difference Time-
Domain Method, 3rd edn. Artech House, 2005.
[3] U. S. Inan and A. S. Inan, Electromagnetic Waves. Prentice-Hall, 2000.
[4] Y. Chen, R. Mittra, and P. Harms, “Finite-difference time-domain algorithm for solving
Maxwell’s equations in rotationally symmetric geometries,” IEEE Trans. Microwave Theory
and Tech., vol. 44, pp. 832–839, 1996.
[5] H. Dib, T. Weller, M. Scardelletti, and M. Imparato, “Analysis of cylindrical transmission lines
with the finite difference time domain method,” IEEE Trans. Microwave Theory and Tech.,
vol. 47, pp. 509–512, 1999.
[6] W. L. Stutzman and G. A. Thiele, Antenna Theory and Design. Wiley, 1998.
[7] H. Dib, T. Weller, and M. Scardelletti, “Analysis of 3-D cylindrical structures using the finite
difference time domain method,” IEEE MTT-S International, vol. 2, pp. 925–928, 1998.
[8] R. Holland, “THREDS: A finite-difference time-domain EMP code in 3D spherical coordi-
nates,” IEEE Trans. Nuc. Sci., vol. NS-30, pp. 4592–4595, 1983.
[9] G. Liu, C. A. Grimes, and K. G. Ong, “A method for FDTD computation of field values at
spherical coordinate singularity points applied to antennas,” IEEE Microwave and Opt. Tech.
Lett., vol. 20, pp. 367–369, 1999.

The FDTD Grid and The Yee Algorithm PDF

Uploaded by

Copyright:

Available Formats

The FDTD Grid and The Yee Algorithm PDF

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

The FDTD Grid and The Yee Algorithm PDF

Uploaded by

Copyright:

Available Formats

4 The FDTD grid and the Yee algorithm

Leapfrog in time n+1

Rearranging terms, we have the time-marching procedure:

4.1 Maxwell’s equations in one dimension

After rearranging, we find:

of the simulation is highly dependent on the first time step 1/2

Sinusoidal Source at x = 0, Snapshot at t = 50 ns

25.2 ns 10.8 ns t = 3.6 ns 18.0 ns t = 32.4 ns

4.1.1 Example 1D simulations

4.2 Maxwell’s equations in two dimensions

We now proceed with discretization of Maxwell’s equations based on the staggered

j+1/2 (i+1, j+1)

4.2.1 Transverse electric (TE) mode

n+1 n  n+1/2 n+1/2 

4.2.2 Transverse magnetic (TM) mode

4.2.3 Example 2D simulations

4 Gaussian Source Sinusoidal Source

4 Gaussian Source Sinusoidal Source

Snapshots at t = m × 20 Δ t Snapshot at t = 200 Δ t

0.2 r–1/2 0.8

Table 4.1 Field component locations in the Yee cell

Field x location y location z location

4.3 FDTD expressions in three dimensions

component Equations (4.4) with second-order centered differences:

4.3.1 Example 3D simulation

4.4 FDTD algorithm for lossy media

4.4.1 1D waves in lossy media: waves on lossy transmission lines

Discretized at (i + 1/2, n):

Discretized at (i, n + 1/2):

Interleaved leapfrog for lossy line

The corresponding FDTD algorithms for 1D TE and 1D TM electromagnetic waves

25.2 ns 10.8 ns t = 3.6 ns 18.0 ns t = 32.4 ns

Lossy, σ = 0.001 S/m

4.4.2 2D and 3D waves in lossy media

Table 4.2 Coefficients for Lossy Media

Update Equation Lossless C1 Lossless C2 Lossy C1 Lossy C2

n+1 = C1 n + C2 ∇FD n+1/2

4.5 Divergence-free nature of the FDTD algorithm

where, for example, Bz simply denotes the z component

4.6 The FDTD method in other coordinate systems

4.6.1 2D polar coordinates

(4.30c), must be treated as if the “field” r φ is to be discretized; hence, our second-order

Example polar simulation

4.6.2 2D cylindrical coordinates

(i, k) z (i+1, k) (i, k) r (i+1, k)

⎡ n+1/2 n+1/2 n+1/2 n+1/2 ⎤

and for the TE mode:

Noting that ri=1/2 = r/2, we find the update equation for z at r = 0:

On-axis cell Regular cell

Example cylindrical simulation

Log10 (z) Amplitude

Since we have chosen a dipole radiator at the center of what is effectively a 3D

4.6.3 3D cylindrical coordinates

rules as for the 2D cylindrical case, we find the update equation:

4.6.4 3D spherical coordinates

n+1 n t sin θ Nθ −1/2 (Nφ n+1/2

4.1. A lossy high-speed chip-to-chip interconnect. Consider applying the inter-

where τ t = 0.1 ns.

n+1 n n+1/2 n+1/2

n+1 = C1 n + C2 ∇FD n+1/2

⎡ n+1/2 n+1/2 n+1/2 n+1/2 ⎤

n+1 n t sin θ Nθ −1/2 (Nφ n+1/2