CHAPTER 1

Introduction to FDTD

Computational electromagnetics (CEM) has evolved rapidly during the past decade to a
point where now extremely accurate predictions can be given for a variety of electro-
magnetic problems, including the scattering cross-section of radar targets and the precise
design of antennas and microwave devices. In general, commonly used CEM methods today
can be classified into two categories. The first is based on differential equation (DE) meth-
ods, whereas the second is based on integral equation (IE) methods. Both IE and DE solution
methods are based on the applications of Maxwell’s equations and the appropriate boundary
conditions associated with the problem to be solved. The IE methods in general provide
approximations for IEs in terms of finite sums, whereas the DE methods provide approx-
imations for DEs as finite differences.
In previous years, most numerical electromagnetic analysis has taken place in the
frequency domain where time-harmonic behavior is assumed. Frequency domain was
favored over time domain because a frequency-domain approach is more suitable for
obtaining analytical solutions for canonical problems, which are used to verify the numerical
results obtained as a first step before depending on a newly developed numerical method for
generating data for real-world applications. Furthermore, the experimental hardware avail-
able for making measurements in past years was largely confined to the frequency-domain
approach.
The recent development of faster and more powerful computational resources allowed for
more advanced time-domain CEM models. More focus is directed toward DE time-domain
approaches as they are easier to formulate and to adapt in computer simulation models
without complex mathematics. They also provide more physical insight to the characteristics
of the problems.
Therefore, an in-depth analysis and implementation of the commonly used time-domain
DE approach, namely, the finite-difference time-domain (FDTD) method for CEM applica-
tions, is covered in this book, along with applications related to antenna designs, microwave
filter designs, and radar cross-section analysis of three-dimensional targets.
The FDTD method has gained tremendous popularity in the past decade as a tool for
solving Maxwell’s equations. It is based on simple formulations that do not require complex
asymptotic or Green’s functions. Although it solves the problem in time, it can provide
frequency-domain responses over a wide band using the Fourier transform. It can easily
handle composite geometries consisting of different types of materials including dielectric,

1
2 CHAPTER 1 ● Introduction to FDTD

magnetic, frequency-dependent, nonlinear, and anisotropic materials. The FDTD technique

is easy to implement using parallel computation algorithms. These features of the FDTD
method have made it the most attractive technique of CEM for many microwave devices and
antenna applications.
FDTD has been used to solve numerous types of problems arising while studying many
applications, including the following:
● scattering, radar cross-section
● microwave circuits, waveguides, fiber optics
● antennas (radiation, impedance)
● propagation
● medical applications
● shielding, coupling, electromagnetic compatibility (EMC), electromagnetic pulse (EMP)
protection
● nonlinear and other special materials
● geological applications
● inverse scattering
● plasma

1.1 The finite-difference time-domain method

basic equations
The starting point for the construction of an FDTD algorithm is Maxwell’s time-domain
equations. The differential time-domain Maxwell’s equations needed to specify the field
behavior over time are
!
@D !
!
rH ¼ þ J; ð1:1aÞ
@t
!
! @B !
rE ¼  M; ð1:1bÞ
@t
!
r  D ¼ re ; ð1:1cÞ
!
r  B ¼ rm ; ð1:1dÞ
! !
where E is the electric field strength vector
!
in volts per meter, D is the electric displacement
vector in! coulombs per square meter, H is the magnetic field strength vector!in amperes per
meter, B is the magnetic flux density vector in webers !
per square meter, J is the electric
current density vector in amperes per square meter, M is the magnetic current density vector
in volts per square meter, re is the electric charge density in coulombs per cubic meter, and
rm is the magnetic charge density in webers per cubic meter.
Constitutive relations are necessary to supplement Maxwell’s equations and characterize
the material media. Constitutive relations for linear, isotropic, and nondispersive materials
can be written as
! !
D ¼ eE; ð1:2aÞ
! !
B ¼ mH ; ð1:2bÞ
 1.1 ● The finite-difference time-domain method basic equations 3

where e is the permittivity, and m is the permeability of the material. In free space

e ¼ e0  8:854  1012 farad=meter;

m ¼ m0 ¼ 4p  107 henry=meter:
We only need to consider curl equations (1.1a) and (1.1b) while deriving FDTD
equations because the divergence equations can be !
satisfied by the developed FDTD
updating!equations
!
[1]. The electric current density J is
!
the !sum !of the! conduction current
density J c ¼ se E and the !impressed
!
current
!
density
!
J i as J! ¼ J c þ J i : Similarly, for the
magnetic current density, M ¼ M c þ M i ; where M c ¼ sm H: Here se is the electric con-
ductivity in siemens per meter, and sm is the magnetic conductivity in ohms per meter. Upon
decomposing the current densities in (1.1) to conduction and impressed components and by
using the constitutive relations (1.2) we can rewrite Maxwell’s curl equations as
!
! @E ! !
rH ¼e þ se E þ J i ; ð1:3aÞ
@t
!
! @H ! !
r  E ¼ m  sm H  M i : ð1:3bÞ
@t
! ! ! !
This formulation treats only the electromagnetic fields E and H and not the fluxes D and B:
All four constitutive parameters e, m, se, and sm are present so that any linear isotropic material
can be specified. Treatment of electric and magnetic sources is included through the impressed
currents. Although only the curl equations are used and the divergence equations are not part of
the FDTD formalism, the divergence!
equations
! !
can! be used as a test on !the predicted
!
field
response, !so that!after forming D ¼ eE and B ¼ mH from the predicted E and H fields, the
resulting D and B must satisfy the divergence equations.
Equation (1.3) is composed of two vector equations, and each vector equation can be
decomposed into three scalar equations for three-dimensional space. Therefore, Maxwell’s
curl equations can be represented with the following six scalar equations in a Cartesian
coordinate system (x, y, z):
 
@Ex 1 @Hz @Hy
¼   sx Ex  J ix ;
e
ð1:4aÞ
@t ex @y @z
 
@Ey 1 @Hx @Hz
¼   sy Ey  J iy ;
e
ð1:4bÞ
@t ey @z @x
 
@Ez 1 @Hy @Hx
¼   sz Ez  J iz ;
e
ð1:4cÞ
@t ez @x @y
 
@Hx 1 @Ey @Ez
¼   sx Hx  Mix ;
m
ð1:4dÞ
@t mx @z @y
 
@Hy 1 @Ez @Ex
¼   sy Hy  Miy ;
m
ð1:4eÞ
@t my @x @z
 
@Hz 1 @Ex @Ey
¼   sm H
z z  M iz : ð1:4f Þ
@t mz @y @x
4 CHAPTER 1 ● Introduction to FDTD

The material parameters ex, ey, and ez are associated with electric field components Ex,
Ey, and Ez through constitutive relations Dx ¼ exEx, Dy ¼ eyEy, and Dz ¼ ezEz, respectively.
Similarly, the material parameters mx, my, and mz are associated with magnetic field compo-
nents Hx, Hy, and Hz through constitutive relations Bx ¼ mxHx, By ¼ myHy, and Bz ¼ mzHz,
respectively. Similar decompositions for other orthogonal coordinate systems are possible,
but they are less attractive from the applications point of view.
The FDTD algorithm divides the problem geometry into a spatial grid where electric and
magnetic field components are placed at certain discrete positions in space, and it solves
Maxwell’s equations in time at discrete time instances. This can be implemented by first
approximating the time and space derivatives appearing in Maxwell’s equations by finite
differences and next by constructing a set of equations that calculate the values of fields at a
future time instant from the values of fields at a past time instant, therefore constructing a
time marching algorithm that simulates the progression of the fields in time [2].

1.2 Approximation of derivatives by finite differences

An arbitrary continuous function can be sampled at discrete points, and the discrete function
becomes a good approximation of the continuous function if the sampling rate is sufficient
relative to the function’s variation. Sampling rate determines the accuracy of operations
performed on the discrete function that approximates the operations on the continuous
functions as well. However, another factor that determines the accuracy of an operation on a
discrete function is the choice of the discrete operator. Most of the time it is possible to use
more than one way of performing an operation on a discrete function. Here we will consider
the derivative operation.
Consider the continuous function given in Figure 1.1(a–c), sampled at discrete points.
The expression for the derivative of f ðxÞ at point x can be written as
f ðx þ DxÞ  f ðxÞ
f 0ðxÞ ¼ lim : ð1:5Þ
Dx!0 Dx
However, since Dx is a nonzero fixed number, the derivative of f ðxÞ can be approxi-
mately taken as
f ðx þ DxÞ  f ðxÞ
f 0ðxÞ  : ð1:6Þ
Dx
The derivative of f ðxÞ is the slope of the dashed line as illustrated in Figure 1.1(a).
Equation (1.6) is called the forward difference formula since one forward point f (x þ Dx) is
used to evaluate f 0ðxÞ together with f ðxÞ.
It is evident that another formula for an approximate f 0ðxÞ can be obtained by using a
backward point f ðx  DxÞ rather than the forward point f ðx þ DxÞ as illustrated in Figure 1.1(b),
which can be written as
f ðxÞ  f ðx  DxÞ
f 0ðxÞ  : ð1:7Þ
Dx
This equation is called the backward difference formula due to the use of the backward
point fðx  DxÞ.
 1.2 ● Approximation of derivatives by finite differences 5

f (x + Δ x)

f (x – Δ x)

f (x)

(a) x – Δx x x + Δx

f (x + Δ x)

f (x – Δ x)
f (x)

(b) x – Δx x x + Δx

f (x +Δ x)

f (x – Δ x)

f (x)

(c) x – Δx x x + Δx

Figure 1.1 (a) Approximation of the derivative of f ðxÞ at x by finite differences: forward
difference; (b) approximation of the derivative of f ðxÞ at x by finite differences:
backward difference; (c) approximation of the derivative of f ðxÞ at x by finite
differences: central difference.
6 CHAPTER 1 ● Introduction to FDTD

The third way of obtaining a formula for an approximate f 0 ðxÞ is by averaging the for-
ward difference and backward difference formulas, such that

f ðx þ DxÞ  f ðx  DxÞ
f 0 ðxÞ  : ð1:8Þ
Dx

Equation (1.8) is called the central difference formula since both the forward and back-
ward points around the center are used. The line representing the derivative of f ðxÞ calcu-
lated using the central difference formula is illustrated in Figure 1.1(c). It should be noted
that the value of the function f ðxÞ at x is not used in central difference formula.
Examination of Figure 1.1 immediately reveals that the three different schemes yield
different values for f 0 ðxÞ, with an associated amount of error. The amount of error introduced
by these difference formulas can be evaluated analytically by using the Taylor series
approach. For instance, the Taylor series expansion of f ðx þ DxÞ can be written as

ðDxÞ2 00 ðDxÞ3 000 ðDxÞ4 0000

f ðx þ DxÞ ¼ f ðxÞ þ Dxf 0 ðxÞ þ f ðxÞ þ f ðxÞ þ f ðxÞ þ ⋯: ð1:9Þ
2 6 24

This equation gives an exact expression for f ðx þ DxÞ as a series in terms of Dx and
derivatives of f ðxÞ, if f ðxÞ satisfies certain conditions and infinite number of terms, theore-
tically, are being used. Equation (1.9) can be rearranged to express f 0 ðxÞ as

f ðx þ DxÞ  f ðxÞ Dx 00 ðDxÞ2 000 ðDxÞ3 0000

f 0 ðxÞ ¼  f ðxÞ  f ðxÞ  f ðxÞ  ⋯: ð1:10Þ
Dx 2 6 24

Here it can be seen that the first term on the right-hand side of (1.10) is the same as the
forward difference formula given by (1.6). The sum of the rest of the terms is the difference
between the approximate derivative given by the forward difference formula and the exact
derivative f 0 ðxÞ, and hence is the amount of error introduced by the forward difference
formula. Equation (1.10) can be rewritten as

f ðx þ DxÞ  f ðxÞ
f 0 ðxÞ ¼ þ OðDxÞ; ð1:11Þ
Dx

where O(Dx) represents the error term. The most significant term in O(Dx) is Dx/2, and the
order of Dx in this most significant term is one. Therefore, the forward difference formula is
first-order accurate. The interpretation of first-order accuracy is that the most significant
term in the error introduced by a first-order accurate formula is proportional to the
sampling period. For instance, if the sampling period is decreased by half, the error reduces
by half.
A similar analysis can be performed for evaluation of the error of the backward formula
starting with the Taylor series expansion of f ðx  DxÞ:

ðDxÞ2 00 ðDxÞ3 000 ðDxÞ4 0000

f ðx  DxÞ ¼ f ðxÞ  Dxf 0 ðxÞ þ f ðxÞ  f ðxÞ þ f ðxÞ þ ⋯: ð1:12Þ
2 6 24
 1.2 ● Approximation of derivatives by finite differences 7

This equation can be rearranged to express f 0 ðxÞ as

f ðxÞ  f ðx  DxÞ Dx 00 ðDxÞ2 000 ðDxÞ3 0000

f 0 ðxÞ ¼ þ f ðxÞ  f ðxÞ þ f ðxÞ  ⋯: ð1:13Þ
Dx 2 6 24

The first term on the right-hand side of (1.13) is the same as the backward difference
formula and the sum of the rest of the terms represents the error introduced by the backward
difference formula to the exact derivative of f ðxÞ, such that

f ðxÞ  f ðx  DxÞ
f 0ðxÞ ¼ þ OðDxÞ: ð1:14Þ
Dx

The order of Dx in the most significant term of O(Dx) is one; hence the backward dif-
ference formula is first-order accurate.
The difference between the Taylor series expansions of f ðx þ DxÞ and f ðx  DxÞ can be
expressed using (1.9) and (1.12) as

2ðDxÞ3 000
f ðx þ DxÞ  f ðx  DxÞ ¼ 2Dxf 0 ðxÞ þ f ðxÞ þ ⋯: ð1:15Þ
6

This equation can be rearranged to express f 0 ðxÞ as

f ðx þ DxÞ  f ðx  DxÞ ðDxÞ2 000 f ðx þ DxÞ  f ðx  DxÞ  

f 0 ðxÞ ¼  f ðxÞ þ ⋯ ¼ þ O ðDxÞ2 ;
2Dx 6 2Dx
ð1:16Þ

where the first term on right-hand side is the same as the central difference formula given in
(1.8) and the order of Dx in the most significant term of the error O((Dx)2) is two; hence the
central difference formula is second-order accurate. The interpretation of second-order
accuracy is that the most significant term in the error introduced by a second-order accurate
formula is proportional to the square of sampling period. For instance, if the sampling period
is decreased by half, the error is reduced by a factor of four. Hence, a second-order accurate
formula such as the central difference formula is more accurate than a first-order accurate
formula.
For an example, consider a function f ðxÞ ¼ sin(x)e0.3x as displayed in Figure 1.2(a).
The exact first-order derivative of this function is

f 0ðxÞ ¼ cosðxÞe0:3x  0:3 sinðxÞe0:3x :

This function f ðxÞ is sampled with a sampling period Dx ¼ p/5, and approximate deri-
vatives are calculated for f ðxÞ using the forward difference, backward difference, and central
difference formulas. The derivative of the function f 0 ðxÞ and its finite difference approx-
imations are plotted in Figure 1.2(b). The errors introduced by the difference formulas,
which are the differences between f 0 ðxÞ and its finite difference approximations, are plotted
in Figure 1.2(c) for the sampling interval Dx ¼ p/5. It is evident that the error introduced by
the central difference formula is smaller than the errors introduced by the forward difference
8 CHAPTER 1 ● Introduction to FDTD

0.8

0.6

0.4

0.2
f (x)

–0.2

–0.4

–0.6

–0.8

–1
0 5 10 15
(a) x

1
exact
0.8 forward difference
backward difference
0.6 Δ x = p /5
central difference

0.4

0.2
f ′ (x)

–0.2

–0.4

–0.6

–0.8

–1
0 5 10 15
(b) x

Figure 1.2 (a) f ðxÞ, f 0 ðxÞ, and differences between f 0 ðxÞ and finite difference approximations of
f 0 ðxÞ for Dx ¼ p/5 and Dx ¼ p/10: f ðxÞ ¼ sin(x)e0.3x; (b) f ðxÞ, f 0 ðxÞ, and differences between
f 0 ðxÞ and finite difference approximations of f 0 ðxÞ for Dx ¼ p/5 and Dx ¼ p/10: f 0 ðxÞ ¼
cos(x)e0.3x – 0.3 sin(x)e0.3x and finite difference approximations of f 0 ðxÞ for Dx ¼ p/5;
(c) f ðxÞ, f 0 ðxÞ, and differences between f 0 ðxÞ and finite difference approximations of f 0 ðxÞ
for Dx ¼ p/5 and Dx ¼ p/10: error(f 0 ðxÞ) for Dx ¼ p/5; (d) f ðxÞ, f 0 ðxÞ, and differences between
f 0 ðxÞ and finite difference approximations of f 0 ðxÞ for Dx ¼ p/5 and Dx ¼ p/10: error(f 0 ðxÞ)
for Dx ¼ p/10.
 1.2 ● Approximation of derivatives by finite differences 9

0.2
forward difference
0.15 Δ x = p /5 backward difference
central difference
0.1

0.05
Error [ f ' (x)]

–0.05

–0.1

–0.15

–0.2
0 5 10 15
(c) x
0.2
forward difference
0.15 Δ x = p /10 backward difference
central difference
0.1

0.05
Error [ f ' (x)]

–0.05

–0.1

–0.15

–0.2
0 5 10 15
(d) x

Figure 1.2 (Continued )

and backward difference formulas. Furthermore, the errors introduced by the difference
formulas for the sampling period Dx ¼ p/10 are plotted in Figure 1.2(d). It can be realized
that as the sampling period is halved, the errors of the forward difference and backward
difference formulas are halved as well, and the error of the central difference formula is
reduced by a factor of four.
The MATLAB code calculating f ðxÞ and its finite difference derivatives, and generat-
ing the plots in Figure 1.2 is shown in Listing 1.1.
10 CHAPTER 1 ● Introduction to FDTD

Listing 1.1 MATLAB code generating Figure 1.2(a–d)

1.2 ● Approximation of derivatives by finite differences 11
12 CHAPTER 1 ● Introduction to FDTD

Values of the function f ðxÞ at two neighboring points around x have been used to obtain a
second-order accurate central difference formula to approximate f 0(x). It is possible to obtain
formulas with higher orders of accuracy by including a larger number of neighboring
points in the derivation of a formula for f 0 ðxÞ. However, although there are FDTD
formulations developed based on higher-order accurate formulas, the conventional FDTD is
based on the second-order accurate central difference formula, which is found to be suffi-
ciently accurate for most electromagnetics applications and simple in implementation and
understanding.
It is possible to obtain finite-difference formulas for approximating higher-order
derivatives as well. For instance, if we take the sum of the Taylor series expansions of
f (x þ Dx) and f (x – Dx) using (1.9) and (1.12), we obtain

ðDxÞ4 0000
f ðx þ DxÞ þ f ðx  DxÞ ¼ 2f ðxÞ þ ðDxÞ2 f 00 ðxÞ þ f ðxÞ þ ⋯: ð1:17Þ
12

After rearranging the equation to have f 00(x) on the left-hand side, we get

f ðx þ DxÞ  2f ðxÞ þ f ðx  DxÞ ðDxÞ2 0000

f 00 ðxÞ ¼ 2
 f ðxÞ þ ⋯
ðDxÞ 12

f ðx þ DxÞ  2f ðxÞ þ f ðx  DxÞ

¼ 2
þ O ðDxÞ2 : ð1:18Þ
ðDxÞ

Using (1.18) we can obtain a central difference formula for the second-order derivative
f 00 (x) as

f ðx þ DxÞ  2f ðxÞ þ f ðx  DxÞ

f 00 ðxÞ  ; ð1:19Þ
ðDxÞ2

which is second-order accurate due to O(Dx)2.

Similarly, some other finite difference formulas can be obtained for the first- and
second-order derivatives with different orders of accuracy based on different sampling
points. A list of finite difference formulas is given for the first- and second-order derivatives
in Table 1.1 as a reference (Figure 1.3).
 1.3 ● FDTD updating equations for three-dimensional problems 13

Table 1.1 Finite difference formulas for the first- and second-order derivatives where the function
f with the subscript i is an abbreviation of f ðxÞ. Similar notation can be implemented
for f (x þ Dx), f (x þ 2Dx), etc., as shown in Figure 1.3. FD, forward difference; BD,
backward difference; CD, central difference.
2
Derivative ∂f=∂x Derivative ∂ f=∂x2

Difference scheme Type error Difference scheme Type error

fiþ1  fi FD OðD xÞ fiþ2  2fiþ1 þ fi FD OðD xÞ
Dx ðDxÞ2
fi  fi1 BD OðD xÞ fiþ2  2fi1 þ fi2 BD OðD xÞ
Dx ðDxÞ2
fiþ1  fi1 CD OððD xÞ2 Þ fiþ1  2fi þ fi1 CD OððD xÞ2 Þ
2D x ðDxÞ2
fiþ2 þ 4fiþ1  3fi FD OððD xÞ2 Þ CD OððD xÞ4 Þ
fiþ2 þ 16fiþ1  30fi þ 16fi1  fi2
2D x
12ðDxÞ2
3fiþ1  4fi1 þ fi2 BD OððD xÞ2 Þ
2D x
CD OððD xÞ4 Þ
fiþ2 þ 8fiþ1  8fi1 þ fi2
12D x

fi−2 fi−1 fi fi+1 fi+2

x − 2Δx x − Δx x x + Δx x + 2Δx

Figure 1.3 Sample points of f ðxÞ.

1.3 FDTD updating equations for three-dimensional

problems
In 1966, Yee originated a set of finite-difference equations for the time-dependent
Maxwell’s curl equations system [2]. These equations can be represented in discrete form,
both in space and time, employing the second-order accurate central difference formula. As
mentioned before, the electric and magnetic field components are sampled at discrete posi-
tions both in time and space. The FDTD technique divides the three-dimensional problem
geometry into cells to form a grid. Figure 1.4 illustrates an FDTD grid composed of
(Nx  Ny  Nz) cells. A unit cell of this grid is called a Yee cell. Using rectangular Yee
cells, a stepped or ‘‘staircase’’ approximation of the surface and internal geometry of the
structure of interest is made with a space resolution set by the size of the unit cell.
The discrete spatial positions of the field components have a specific arrangement in the
Yee cell, as demonstrated in Figure 1.5. The electric field vector components are placed at
14 CHAPTER 1 ● Introduction to FDTD

(Nx, Ny, Nz)

(i, j, k)

y (1, 1, 1)
x

Figure 1.4 A three-dimensional FDTD computational space composed of (Nx  Ny  Nz)

Yee cells.

Δx
node(i + 1, j + 1, k + 1)

Δz

Ez(i, j, k) Hz(i, j, k) Hy(i, j, k)

Hz(i, j, k)
Ey(i, j, k)
node(i, j, k) Δy
z Ex(i, j, k)
y
x

Figure 1.5 Arrangement of field components on a Yee cell indexed as (i, j, k).

the centers of the edges of the Yee cells and oriented parallel to the respective edges, and the
magnetic field vector components are placed at the centers of the faces of the Yee cells and
are oriented normal to the respective faces. This provides a simple picture of three-
dimensional space being filled by an interlinked array of Faraday’s law and Ampere’s law
contours. It can be easily noticed in Figure 1.5 that each magnetic field vector is surrounded
 1.3 ● FDTD updating equations for three-dimensional problems 15

by four electric field vectors that are curling around the magnetic field vector, thus simu-
lating Faraday’s law. Similarly, if the neighboring cells are also added to the picture, it
would be apparent that each electric field vector is surrounded by four magnetic field vectors
that are curling around the electric field vector, thus simulating Ampere’s law.
Figure 1.5 shows the indices of the field components, which are indexed as (i, j, k),
associated with a cell indexed as (i, j, k). For a computational domain composed of uniform
Yee cells having dimension Dx in the x direction, Dy in the y direction, and Dz in the z
direction, the actual positions of the field components with respect to an origin coinciding
with the position of the node (1, 1, 1) can easily be calculated as

Ex ði; j; kÞ ) ðði  0:5ÞDx; ðj  1ÞDy; ðk  1ÞDzÞ;

Ey ði; j; kÞ ) ðði  1ÞDx; ðj  0:5ÞDy; ðk  1ÞDzÞ;
Ez ði; j; kÞ ) ðði  1ÞDx; ðj  1ÞDy; ðk  0:5ÞDzÞ;
Hx ði; j; kÞ ) ðði  1ÞDx; ðj  0:5ÞDy; ðk  0:5ÞDzÞ;
Hy ði; j; kÞ ) ðði  0:5ÞDx; ðj  1ÞDy; ðk  0:5ÞDzÞ;
Hz ði; j; kÞ ) ðði  0:5ÞDx; ðj  0:5ÞDy; ðk  1ÞDzÞ:

The FDTD algorithm samples and calculates the fields at discrete time instants; how-
ever, the electric and magnetic field components are not sampled at the same time instants.
For a time-sampling period Dt, the electric field components are sampled at time instants
0, Dt, 2Dt, . . . , nDt, . . . ; however, the magnetic field components are sampled at time
instants 12 Dt; ð1 þ 12ÞDt; . . . ; ðn þ 12ÞDt; . . . . Therefore, the electric field components are
calculated at integer time steps, and magnetic field components are calculated at half-integer
time steps, and they are offset from each other by Dt/2. The field components need to be
referred not only by their spatial indices which indicate their positions in space, but also by
their temporal indices, which indicate their time instants. Therefore, a superscript notation is
adopted to indicate the time instant. For instance, the z component of an electric field vector
positioned at ((i – 1)Dx, ( j – 1)Dy, (k – 0.5)Dz) and sampled at time instant nDt is referred to
as Ezn ði; j; kÞ: Similarly, the y component of a magnetic field vector positioned at ((i – 0.5)Dx,
( j – 1)Dy, (k – 0.5)Dz) and sampled at time instant ðn þ 12ÞDt is referred to as Hynþ1=2 ði; j; kÞ:
The material parameters (permittivity, permeability, electric, and magnetic con-
ductivities) are distributed over the FDTD grid and are associated with field components;
therefore, they are indexed the same as their respective field components. For instance,
Figure 1.6 illustrates the indices for the permittivity and permeability parameters.
The electric conductivity is distributed and indexed the same as the permittivity, and the
magnetic conductivity is distributed and indexed the same as the permeability.
Having adopted an indexing scheme for the discrete samples of field components in both
time and space, Maxwell’s curl equations (1.4) that are given in scalar form can be expressed
in terms of finite differences. For instance, consider again (1.4a):
 
@Ex 1 @Hz @Hy
¼   sex Ex  J ix :
@t ex @y @z
16 CHAPTER 1 ● Introduction to FDTD

ex(i, j + 1, k + 1)

ey(i, j, k + 1)
mz(i, j, k + 1)
ey(i + 1, j, k + 1)
ex(i, j, k + 1)
uy(i, j + 1, k)
ez(i, j + 1, k) e (i + 1, j + 1, k)
mx(i + 1, j, k) z
ez(i, j, k) mz(i, j, k)
my(i, j, k)
ez(i + 1, j, k)
ex(i, j + 1, k)

ey(i, j, k) mz(i, j, k)
ey(i + 1, j, k)
node(i, j, k)
z ex(i, j, k)

y
x

Figure 1.6 Material parameters indexed on a Yee cell.

The derivatives in this equation can be approximated by using the central difference
formula with the position of Ex(i, j, k) being the center point for the central difference
formula in space and time instant ðn þ 12ÞDt as being the center point in time. Considering the
field component positions given in Figure 1.7, we can write
nþ12 nþ12
Exnþ1 ði; j; kÞ  Exn ði; j; kÞ 1 Hz ði; j; kÞ  Hz ði; j  1; kÞ
¼
Dt ex ði; j; kÞ Dy
nþ1 nþ1
1 Hy 2 ði; j; kÞ  Hy 2 ði; j; k  1Þ

ex ði; j; kÞ Dz
sex ði; j; kÞ nþ12 1 nþ1
 Ex ði; j; kÞ  J ix 2 ði; j; kÞ: ð1:20Þ
ex ði; j; kÞ ex ði; j; kÞ

It has already been mentioned that the electric field components are defined at integer
time steps; however, the right-hand side of (1.20) includes an electric field term at time
instant ðn þ 12ÞDt, that is, Exnþ1=2 ði; j; kÞ. This term can be written as the average of the terms
at time instants (n þ 1)Dt and nDt, such that

nþ1 Exnþ1 ði; j; kÞ þ Exn ði; j; kÞ

Ex 2 ði; j; kÞ ¼ : ð1:21Þ
2

Using (1.21) in (1.20) and arranging the terms such that the future term Exnþ1 ði; j; kÞ is
kept on the left side of the equation and the rest of the terms are moved to the right-hand side
of the equation, we can write
 1.3 ● FDTD updating equations for three-dimensional problems 17

z
Δx
y
x node(i + 1, j + 1, k + 1)
Δy

Hy(i, j, k)
Δz

node(i, j, k) Hz(i, j, k)
node(i, j − 1, k) Hz(i, j − 1, k)
Ex(i, j, k)

Hy(i, j, k − 1)

node(i, j, k − 1)

Figure 1.7 Field components around Ex(i, j, k).

2ex ði; j; kÞ þ Dtsex ði; j; kÞ nþ1 2ex ði; j; kÞ  Dtsex ði; j; kÞ n

Ex ði; j; kÞ ¼ Ex ði; j; kÞ
2ex ði; j; kÞ 2ex ði; j; kÞ

Dt  
nþ1 nþ1
þ Hz 2 ði; j; kÞ  Hz 2 ði; j  1; kÞ
ex ði; j; kÞDy

Dt  
nþ1 nþ1
 Hy 2 ði; j; kÞ  Hy 2 ði; j; k  1Þ
ex ði; j; kÞDz

Dt nþ1
 J 2 ði; j; kÞ
ex ði; j; kÞ ix
ð1:22Þ

After some manipulations, we get

2ex ði; j; kÞ  Dtsex ði; j; kÞ n

Exnþ1 ði; j; kÞ ¼ E ði; j; kÞ
2ex ði; j; kÞ þ Dtsex ði; j; kÞ x

2Dt  
nþ12 n þ12
þ
Hz ði; j; kÞ  Hz ði; j  1; kÞ
2ex ði; j; kÞ þ Dtsex ði; j; kÞ Dy

2Dt  
nþ1 nþ1

Hy 2 ði; j; kÞ  Hy 2 ði; j; k  1Þ
2ex ði; j; kÞ þ Dtsx ði; j; kÞ Dz
e

2Dt nþ1 ð1:23Þ

 J ix 2 ði; j; kÞ:
2ex ði; j; kÞ þ Dtsx ði; j; kÞ
e
18 CHAPTER 1 ● Introduction to FDTD

y
x
Ey(i, j, k + 1)
Δx
node(i, j, k + 1)
Δy

Ez(i, j + 1, k)
Ez(i, j, k)
Hx(i, j, k)
Δz

node(i, j + 1, k)

node(i, j, k) Ey(i, j, k)

Figure 1.8 Field components around Hx(i, j, k).

The form of (1.23) demonstrates how the future value of an electric field component can
be calculated by using the past values of the electric field component, the magnetic field
components, and the source components. This form of an equation is called an FDTD
updating equation. Updating equations can easily be obtained for calculating Eynþ1 ði; j; kÞ
starting from (1.4b) and Eznþ1 ði; j; kÞ starting from (1.4c) following the same methodology
that has been used to obtain (1.23).
Similarly, updating equations can be obtained for magnetic field components following
the same methodology. However, while applying the central difference formula to the time
derivative of the magnetic field components, the central point in time shall be taken as nDt.
For instance, (1.4c), which is
 
@Hx 1 @Ey @Ez
¼   sx Hx  Mix ;
m
@t mx @z @y

can be approximated using finite differences based on the field positions (as shown in
Figure 1.8) as

nþ1 n1
Hx 2 ði; j; kÞ  Hx 2 ði; j; kÞ 1 Eyn ði; j; k þ 1Þ  Eyn ði; j; kÞ
¼
Dt mx ði; j; kÞ Dz

1 Ezn ði; j þ 1; kÞ  Ezn ði; j; kÞ


mx ði; j; kÞ Dy

sm
x ði; j; kÞ n 1
 H ði; j; kÞ  M n ði; j; kÞ: ð1:24Þ
mx ði; j; kÞ x mx ði; j; kÞ ix
 1.3 ● FDTD updating equations for three-dimensional problems 19

nþ1
After some manipulations, the future term Hx 2 ði; j; kÞ in (1.24) can be moved to the left-
hand side and the other terms can be moved to the right-hand side such that

nþ1 2mx ði; j; kÞ  Dtsm x ði; j; kÞ n2

1
Hx 2 ði; j; kÞ ¼ H x ði; j; kÞ
2mx ði; j; kÞ þ Dtsm x ði; j; kÞ
2Dt  
þ
E n
ði; j; k þ 1Þ  E n
ði; j; kÞ
2mx ði; j; kÞ þ Dtsm x ði; j; kÞ Dz
y y

2Dt


Ezn ði; j þ 1; kÞ  Ezn ði; j; kÞ
2mx ði; j; kÞ þ Dtsx ði; j; kÞ Dy
m

2Dt
 Mixn ði; j; kÞ: ð1:25Þ
2mx ði; j; kÞ þ Dtsm
x ði; j; kÞ
nþ1
This equation is the updating equation for Hx 2 ði; j; kÞ. Similarly, updating equations can
nþ1 nþ1
easily be obtained for Hy 2 ði; j; kÞ starting from (1.4e) and Hz 2 ði; j; kÞ starting from (1.4f)
following the same methodology used to obtain (1.25).
Finally, (1.4a)–(1.4f) can be expressed using finite differences and can be arranged to
construct the following six FDTD updating equations for the six components of electro-
magnetic fields by introduction of respective coefficient terms:

Exnþ1 ði; j; kÞ ¼ Cexe ði; j; kÞ  Exn ði; j; kÞ

 
nþ1 nþ1
þ Cexhz ði; j; kÞ  Hz 2 ði; j; kÞ  Hz 2 ði; j  1; kÞ
 
nþ1 nþ1
þ Cexhy ði; j; kÞ  Hy 2 ði; j; kÞ  Hy 2 ði; j; k  1Þ
nþ1
þ Cexj ði; j; kÞ  J ix 2 ði; j; kÞ; ð1:26Þ
where
2ex ði; j; kÞ  Dtsex ði; j; kÞ
Cexe ði; j; kÞ ¼ ;
2ex ði; j; kÞ þ Dtsex ði; j; kÞ
2Dt
Cexhz ði; j; kÞ ¼
;
2ex ði; j; kÞ þ Dtsex ði; j; kÞ Dy
ð1:27Þ
2Dt
Cexhy ði; j; kÞ ¼ 
;
2ex ði; j; kÞ þ Dtsex ði; j; kÞ Dz
2Dt
Cexj ði; j; kÞ ¼  :
2ex ði; j; kÞ þ Dtsex ði; j; kÞ

Eynþ1 ði; j; kÞ ¼ Ceye ði; j; kÞ  Eyn ði; j; kÞ

 
nþ1 nþ1
þ Ceyhx ði; j; kÞ  Hx 2 ði; j; kÞ  Hx 2 ði; j; k  1Þ
 
nþ1 nþ1
þ Ceyhz ði; j; kÞ  Hz 2 ði; j; kÞ  Hz 2 ði  1; j; kÞ
nþ1
þ Ceyj ði; j; kÞ  J iy 2 ði; j; kÞ;
20 CHAPTER 1 ● Introduction to FDTD

where

2ey ði; j; kÞ  Dtsey ði; j; kÞ

Ceye ði; j; kÞ ¼ ;
2ey ði; j; kÞ þ Dtsey ði; j; kÞ

2Dt
Ceyhx ði; j; kÞ ¼   ;
2ey ði; j; kÞ þ Dtsey ði; j; kÞ Dz

2Dt
Ceyhz ði; j; kÞ ¼    ;
2ey ði; j; kÞ þ Dtsey ði; j; kÞ Dx

2Dt ð1:28Þ
Ceyj ði; j; kÞ ¼  :
2ey ði; j; kÞ þ Dtsey ði; j; kÞ

Eznþ1 ði; j; kÞ ¼ Ceze ði; j; kÞ  Ezn ði; j; kÞ

 
nþ1 nþ1
þ Cezhy ði; j; kÞ  Hy 2 ði; j; kÞ  Hy 2 ði  1; j; kÞ
 
nþ1 nþ1
þ Cezhx ði; j; kÞ  Hx 2 ði; j; kÞ  Hx 2 ði; j  1; kÞ

nþ1
þ Cezj ði; j; kÞ  J iz 2 ði; j; kÞ;

where

2ez ði; j; kÞ  Dtsez ði; j; kÞ

Ceze ði; j; kÞ ¼ ;
2ez ði; j; kÞ þ Dtsez ði; j; kÞ

2Dt
Cezhy ði; j; kÞ ¼
;
2ez ði; j; kÞ þ Dtsez ði; j; kÞ Dx

2Dt
Cezhx ði; j; kÞ ¼ 
;
2ez ði; j; kÞ þ Dtsez ði; j; kÞ Dy

2Dt ð1:29Þ
Cezj ði; j; kÞ ¼  :
2ez ði; j; kÞ þ Dtsez ði; j; kÞ
n þ12 n1
Hx ði; j; kÞ ¼ Chxh ði; j; kÞ  Hx 2 ði; j; kÞ
 
þ Chxey ði; j; kÞ  Eyn ði; j; k þ 1Þ  Eyn ði; j; kÞ


þ Chxez ði; j; kÞ  Ezn ði; j þ 1; kÞ  Ezn ði; j; kÞ

þ Chxm ði; j; kÞ  Mixn ði; j; kÞ;

 1.3 ● FDTD updating equations for three-dimensional problems 21

where
2mx ði; j; kÞ  Dtsmx ði; j; kÞ
Chxh ði; j; kÞ ¼ ;
2mx ði; j; kÞ þ Dtsmx ði; j; kÞ
2Dt
Chxey ði; j; kÞ ¼
;
2mx ði; j; kÞ þ Dtsm x ði; j; kÞ Dz
2Dt
Chxez ði; j; kÞ ¼ 
;
2mx ði; j; kÞ þ Dtsm
x ði; j; kÞ Dy ð1:30Þ
2Dt
Chxm ði; j; kÞ ¼  :
2mx ði; j; kÞ þ Dtsm
x ði; j; kÞ
n þ12 n1
Hy ði; j; kÞ ¼ Chyh ði; j; kÞ  Hy 2 ði; j; kÞ þ Chyez ði; j; kÞ


 Ezn ði þ 1; j; kÞ  Ezn ði; j; kÞ þ Chyex ði; j; kÞ


 Exn ði; j; k þ 1Þ  Exn ði; j; kÞ þ Chym ði; j; kÞ  Miyn ði; j; kÞ
where
2my ði; j; kÞ  Dtsm
y ði; j; kÞ
Chyh ði; j; kÞ ¼ ;
2my ði; j; kÞ þ Dtsm
y ði; j; kÞ
2Dt
Chyez ði; j; kÞ ¼   ;
2my ði; j; kÞ þ Dtsm
y ði; j; kÞ Dx

2Dt
Cbyex ði; j; kÞ ¼    ;
2my ði; j; kÞ þ Dtsm
y ði; j; kÞ Dz
ð1:31Þ
2Dt
Chym ði; j; kÞ ¼  :
2my ði; j; kÞ þ Dtsm
y ði; j; kÞ
nþ12 n12
Hz ði; j; kÞ ¼ Chzh ði; j; kÞ  Hz ði; j; kÞ


þ Chzex ði; j; kÞ  Exn ði; j þ 1; kÞ  Exn ði; j; kÞ
 
þ Chzey ði; j; kÞ  Eyn ði þ 1; j; kÞ  Eyn ði; j; kÞ
þ Chzm ði; j; kÞ  Mizn ði; j; kÞ;
where
2my ði; j; kÞ  Dtsmz ði; j; kÞ
Chzh ði; j; kÞ ¼ ;
2mz ði; j; kÞ þ Dtsz ði; j; kÞ
m

2Dt
Chzex ði; j; kÞ ¼
;
2mz ði; j; kÞ þ Dtsmz ði; j; kÞ Dy
2Dt
Chzey ði; j; kÞ ¼ 
;
2mz ði; j; kÞ þ Dtsm
z ði; j; kÞ Dx
2Dt
Chzm ði; j; kÞ ¼  :
2mz ði; j; kÞ þ Dtsm
z ði; j; kÞ
22 CHAPTER 1 ● Introduction to FDTD

Start

Set problem space and define parameters

Compute field coefficients

Update magnetic field components

at time instant (n+0.5)Δ t

Update electric field components

Output data at time instant (n+1)Δ t

Apply boundary conditions

Increment time step, n n+1

Yes No
Stop Last iteration?

Figure 1.9 Explicit FDTD procedure.

It should be noted that the first two subscripts in each coefficient refer to the
corresponding field component being updated. For three subscripts coefficients, the third
subscript refers to the type of the field or source (electric or magnetic) that this coefficient is
multiplied by. For four subscripts coefficients, the third and fourth subscripts refer to the
type of the field that this coefficient is multiplied by.
Having derived the FDTD updating equations, a time-marching algorithm can be
constructed as illustrated in Figure 1.9. The first step in this algorithm is setting up the
problem space – including the objects, material types, and sources – and defining any other
parameters that will be used during the FDTD computation. Then the coefficient terms
appearing in (1.26)–(1.31) can be calculated and stored as arrays before the iteration is
started. The field components need to be defined as arrays as well and shall be initialized
with zeros since the initial values of the fields in the problem space in most cases are zeros,
and fields will be induced in the problem space due to sources as the iteration proceeds. At
every step of the time-marching iteration the magnetic field components are updated for time
instant (n þ 0.5)Dt using (1.29)–(1.31); then the electric field components are updated for
time instant (n þ 1)Dt using (1.26)–(1.28). The problem space has a finite size, and specific
boundary conditions can be enforced on the boundaries of the problem space. Therefore, the
field components on the boundaries of the problem space are treated according to the type of
boundary conditions during the iteration. The types of boundary conditions and the
 1.4 ● FDTD updating equations for two-dimensional problems 23

techniques used to integrate them into the FDTD algorithm are discussed in detail in
Chapters 7 and 8. After the fields are updated and boundary conditions are enforced, the
current values of any desired field components can be captured and stored as output data, and
these data can be used for real-time processing or postprocessing to calculate some other
desired parameters. The FDTD iterations can be continued until some stopping criteria are
achieved.

1.4 FDTD updating equations for two-dimensional problems

The FDTD updating equations given in (1.26)–(1.31) can be used to solve three-dimensional
problems. In the two-dimensional case where there is no variation in the problem geometry
and field distributions in one of the dimensions, a simplified set of updating equations can be
obtained starting from the Maxwell’s curl equations system (1.4). Since there is no variation
in one of the dimensions, the derivative terms with respect to that dimension vanish. For
instance, if the problem is z dimension independent, equations (1.4) reduce to
 
@Ex 1 @Hz
¼  sx Ex  J ix ;
e
ð1:32aÞ
@t ex @y
 
@Ey 1 @Hz
¼   sy Ey  J iy ;
e
ð1:32bÞ
@t ey @x
 
@Ez 1 @Hy @Hx
¼   sez Ez  J iz ; ð1:32cÞ
@t ez @x @y
 
@Hx 1 @Ez
¼   sx Hx  Mix ;
m
ð1:32dÞ
@t mx @y
 
@Hy 1 @Ez
¼  sy Hy  Miy ;
m
ð1:32eÞ
@t my @x
 
@Hz 1 @Ex @Ey
¼   sz Hz  Miz :
m
ð1:32f Þ
@t mz @y @x

One should notice that equations (1.32a), (1.32b), and (1.32f) are dependent only on the
terms Ex, Ey, and Hz, whereas equations (1.32c), (1.32d), and (1.32e) are dependent only on
the terms Ez, Hx, and Hy. Therefore, the six equations (1.32) can be treated as two separate
sets of equations. In the first set (1.32a), (1.32b), and (1.32f) – all the electric field compo-
nents are transverse to the reference dimension z; therefore, this set of equations constitutes
the transverse electric to z case – TEz. In the second set – (1.32c), (1.32d), and (1.32e) – all
the magnetic field components are transverse to the reference dimension z; therefore, this set
of equations constitutes the transverse magnetic to z case – TMz. Most two-dimensional
problems can be decomposed into two separate problems, each including separate field
components that are TEz and TMz for the case under consideration. These two problems can
be solved separately, and the solution for the main problem can be achieved as the sum of the
two solutions.
24 CHAPTER 1 ● Introduction to FDTD

Ex(i − 1, j + 2) Ex(i, j + 2) Ex(i + 1, j + 2)

Ey(i − 1, j + 1)

Ey(i + 1, j + 1)

Ey(i + 2, j + 1)
Ey(i, j + 1)
Hz(i, j + 1) Hz(i + 1, j + 1)
Hz(i − 1, j + 1)

Ex(i − 1, j + 1) Ex(i, j + 1) Ex(i + 1, j + 1)

Ey(i − 1, j)

Ey(i + 1, j)

Ey(i + 2, j)
Ey(i, j)
Hz(i − 1, j) Hz(i, j) Hz(i + 1, j)

Ey(i + 1, j − 1)
Ex(i − 1, j) Ex(i, j) Ex(i + 1, j)

Ey(i + 2, j − 1)
Ey(i − 1, j − 1)

Ey(i, j − 1)

Δy
Hz(i − 1, j − 1) Hz(i, j − 1) Hz(i + 1, j − 1)

Ex(i − 1, j − 1) Ex(i, j − 1) Ex(i + 1, j − 1)

Δx

Figure 1.10 Two-dimensional TEz FDTD field components.

The FDTD updating equations for the TEz case can be obtained by applying the central
difference formula to the equations constituting the TEz case based on the field positions
shown in Figure 1.10, which is obtained by projection of the Yee cells in Figure 1.5 on the xy
plane in the z direction. The FDTD updating equations for the TEz case are therefore
obtained as
 
nþ1 nþ1
Exnþ1 ði; jÞ ¼ Cexe ði; jÞ  Exn ði; jÞ þ Cexhz ði; jÞ  Hz 2 ði; jÞ  Hz 2 ði; j  1Þ
nþ1
þ Cexj ði; jÞ  J ix 2 ði; jÞ; ð1:33Þ

where
2ex ði; jÞ  Dtsex ði; jÞ
Cexe ði; jÞ ¼ ;
2ex ði; jÞ þ Dtsex ði; jÞ
2Dt
Cexhz ði; jÞ ¼
;
2ex ði; jÞ þ Dtsex ði; jÞ Dy
2Dt
Cexj ði; jÞ ¼  :
2ex ði; jÞ þ Dtsex ði; jÞ
 
nþ1 nþ1
Eynþ1 ði; jÞ ¼ Ceye ði; jÞ  Eyn ði; jÞ þ Ceyhz ði; jÞ  Hz 2 ði; jÞ  Hz 2 ði  1; jÞ
nþ1
þ Ceyj ði; jÞ  J iy 2 ði; jÞ; ð1:34Þ
 1.4 ● FDTD updating equations for two-dimensional problems 25

where
2ey ði; jÞ  Dtsey ði; jÞ
Ceye ði; jÞ ¼ ;
2ey ði; jÞ þ Dtsey ði; jÞ
2Dt
Ceyhz ði; jÞ ¼    ;
2ey ði; jÞ þ Dtsey ði; jÞ Dx
2Dt ð1:35Þ
Ceyj ði; jÞ ¼  :
2ey ði; jÞ þ Dtsey ði; jÞ
nþ12 n 12

Hz ði; jÞ ¼ Chzh ði; jÞ  Hz ði; jÞ þ Cbzex ði; jÞ  Exn ði; j þ 1Þ  Exn ði; jÞ
 
þ Chzey ði; jÞ  Eyn ði þ 1; jÞ  Eyn ði; jÞ þ Chzm ði; jÞ  Mizm ði; jÞ;

where
2mz ði; jÞ  Dtsm z ði; jÞ
Chzh ði; jÞ ¼ ;
2mz ði; jÞ þ Dtsm z ði; jÞ
2Dt
Chzex ði; jÞ ¼
;
2mz ði; jÞ þ Dtsm z ði; jÞ Dy
2Dt
Chzey ði; jÞ ¼ 
;
2mz ði; jÞ þ Dtsm z ði; jÞ Dx
2Dt
Chzm ði; jÞ ¼  :
2mz ði; jÞ þ Dtsm z ði; jÞ

Since the FDTD updating equations for the three-dimensional case are readily available,
they can be used to derive (1.33), (1.34), and (1.35) by simply setting the coefficients
including 1/Dz to zero. Hence, the FDTD updating equations for the TMz case can be
obtained by eliminating the term Chxey(i, j, k) in (1.29) and eliminating the term Chyex(i, j, k)
in (1.30) based on the field positions shown in Figure 1.11, such that
 
nþ1 nþ1
Eznþ1 ði; jÞ ¼ Ceze ði; jÞ  Ezn ði; jÞ þ Cezhy ði; jÞ  Hy 2 ði; jÞ  Hy 2 ði  1; jÞ
 
nþ1 nþ1 nþ1
þ Cezhx ði; jÞ  Hx 2 ði; jÞ  Hx 2 ði; j  1Þ þ Cezj ði; jÞ  J iz 2 ði; jÞ; ð1:36Þ

where
2ez ði; jÞ  Dtsez ði; jÞ
Ceze ði; jÞ ¼ ;
2ez ði; jÞ þ Dtsez ði; jÞ
2Dt
Cezhy ði; jÞ ¼
;
2ez ði; jÞ þ Dtsez ði; jÞ Dx
2Dt
Cezhx ði; jÞ ¼ 
; ð1:37Þ
2ez ði; jÞ þ Dtsez ði; jÞ Dy
2Dt
Cezj ði; jÞ ¼  :
2ez ði; jÞ þ Dtsez ði; jÞ
nþ1 n1

Hx 2 ði; jÞ ¼ Chxh ði; jÞ  Hx 2 ði; jÞ þ Chxez ði; jÞ  Ezn ði; j þ 1Þ  Ezn ði; jÞ
þ Chxm ði; jÞ  Mixn ði; jÞ;
26 CHAPTER 1 ● Introduction to FDTD

2) 2) 2)
, j+ 2 )
, j+ , j+
−1 j+ +1 +2
i i, i i
E z( E z( E z( E z(

Hx(i − 1, j + 1) Hy(i − 1, j + 2) Hy(i, j + 2) Hy(i + 1, j + 2)

Hx(i + 1, j + 1)

Hx(i + 2, j + 1)
Hx(i, j + 1)
1) 1) 1)
, j+ 1 )
, j+ , j+
−1 j+ +1 +2
i i, i i
E z( E z( E z( E z(
Hy(i − 1, j + 1) Hy(i, j + 2) Hy(i + 1, j + 1)
Hx(i − 1, j)

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

j) j) j)
1, j) 1, 2,
i− i, i+ i+
E z( E z( E z( E z(
Hx(i − 1, j − 1)

Hy(i − 1, j) Hy(i, j) Hx(i + 1, j − 1) Hy(i + 1, j)

Hx(i + 2, j − 1)
Hx(i, j − 1)

Δy
1) 1) 1)
, j− 1) , j− j−
− +1 ,
i −1 i, j i i +2
E z( E z( E z( E z(

Hy(i − 1, j − 1) Hy(i, j − 1) Hy(i + 1, j − 1)

Δx

Figure 1.11 Two-dimensional TMz FDTD field components.

where
2mx ði; jÞ  Dtsmx ði; jÞ
Chxh ði; jÞ ¼ ;
2mx ði; jÞ þ Dtsmx ði; jÞ
2Dt
Chxez ði; jÞ ¼ 
;
2mx ði; jÞ þ Dtsm x ði; jÞ Dy
2Dt ð1:38Þ
Chxm ði; jÞ ¼  :
2mx ði; jÞ þ Dtsm
x ði; jÞ
nþ1 n1

Hy 2 ði; jÞ ¼ Chyh ði; jÞ  Hy 2 ði; jÞ þ Chyez ði; jÞ  Ezn ði þ 1; jÞ  Ezn ði; jÞ
þ Chym ði; jÞ  Miym ði; jÞ;

where
2my ði; jÞ  Dtsm
y ði; jÞ
Chyh ði; jÞ ¼ ;
2my ði; jÞ þ Dtsm
y ði; jÞ
2Dt
Chyez ði; jÞ ¼   ;
2my ði; jÞ þ Dtsm
y ði; jÞ Dx

2Dt
Chym ði; jÞ ¼  :
2my ði; jÞ þ Dtsm
y ði; jÞ
 1.5 ● FDTD updating equations for one-dimensional problems 27

1.5 FDTD updating equations for one-dimensional problems

In the one-dimensional case, there is no variation in the problem geometry and field dis-
tributions in two of the dimensions. For instance, if the y and z dimensions have no variation,
the derivative with respect to the y and z dimensions vanish in Maxwell’s curl equations.
Therefore, the two-dimensional curl equations (1.32a)–(1.32f) reduce to

@Ex 1

¼ sex Ex  J ix ; ð1:39aÞ
@t ex
 
@Ey 1 @Hz
¼   sy Ey  J iy ;
e
ð1:39bÞ
@t ey @x
 
@Ez 1 @Hy
¼  sz Ez  J iz ;
e
ð1:39cÞ
@t ez @x

@Hx 1

¼ smx Hx  Mix ; ð1:39dÞ
@t mx
 
@Hy 1 @Ez
¼  sy Hy  Miy ;
m
ð1:39eÞ
@t my @x
 
@Hz 1 @Ey
¼   sm H
z z  M iz : ð1:39f Þ
@t mz @x

It should be noted that (1.39a) and (1.39d) include time derivatives but not space deri-
vatives. Therefore, these equations do not represent propagating fields, hence the field
components Ex and Hx, which only exist in these two equations, do not propagate. The other
four equations represent propagating fields, and both the electric and magnetic field com-
ponents existing in these equations are transverse to the x dimension. Therefore, transverse
electric and magnetic to x (TEMx) fields exist and propagate as plane waves in the one-
dimensional case under consideration.
Similar to the two-dimensional case, the one-dimensional case as well can be decomposed
into two separate cases, since (1.39b) and (1.39f), which include only the terms Ey and Hz, are
decoupled from (1.39c) and (1.39e), which include only the terms Ez and Hy. FDTD updating
equations can be obtained for (1.39b) and (1.39f) using the central difference formula based on
the field positioning in one-dimensional space as illustrated in Figure 1.12, such that
 
nþ1 nþ1 nþ1
Eynþ1 ðiÞ ¼ Ceye ðiÞ  Eyn ðiÞ þ Ceyhz ðiÞ  Hz 2 ðiÞ  Hz 2 ði  1Þ þ Ceyj ðiÞ  J iy 2 ðiÞ; ð1:40Þ

Ey(i − 1) Hz(i − 1) Ey(i) Hz(i) Ey(i + 1) Hz(i + 1) Ey(i + 2)

Δx

Figure 1.12 One-dimensional FDTD – positions of field components Ey and Hz.

28 CHAPTER 1 ● Introduction to FDTD

where
2ey ðiÞ  Dtsey ðiÞ
Ceye ðiÞ ¼ ;
2ey ðiÞ þ Dtsey ðiÞ
2Dt
Ceyhz ðiÞ ¼    ;
2ey ðiÞ þ Dtsey ðiÞ Dx
2Dt
Ceyj ðiÞ ¼  ;
2ey ðiÞ þ Dtsey ðiÞ

and
 
nþ12 n  12
Hz ðiÞ ¼ Chzh ðiÞ  Hz ðiÞ þ Chzey ðiÞ  Eyn ði þ 1Þ  Eyn ðiÞ þ Chzm ðiÞ  Mizn ðiÞ; ð1:41Þ

where
2mz ðiÞ  Dtsmz ðiÞ
Chzh ðiÞ ¼ ;
2mz ðiÞ þ Dtsmz ðiÞ
2Dt
Chzey ðiÞ ¼ 
;
2mz ðiÞ þ Dtsmz ðiÞ Dx
2Dt
Chzm ðiÞ ¼  :
2mz ðiÞ þ Dtsm
z ðiÞ

Similarly, FDTD updating equations can be obtained for (1.39c) and (1.39e) using the
central difference formula based on the field positioning in one-dimensional space as illu-
strated in Figure 1.13, such that
 
nþ1 nþ1 nþ1
Eznþ1 ðiÞ ¼ Ceze ðiÞ  Ezn ðiÞ þ Cezhy ðiÞ  Hy 2 ðiÞ  Hy 2 ði  1Þ þ Cezj ðiÞ  J iz 2 ðiÞ; ð1:42Þ

where
2ez ðiÞ  Dtsez ðiÞ
Ceze ðiÞ ¼ ;
2ez ðiÞ þ Dtsez ðiÞ
2Dt
Cezhy ðiÞ ¼
;
2ez ðiÞ þ Dtsez ðiÞ Dx
2Dt
Cezj ðiÞ ¼  ;
2ez ðiÞ þ Dtsez ðiÞ

Ez(i − 1) Hy(i − 1) Ez(i) Hy(i) Ez(i + 1) Hy(i + 1) Ez(i + 2)

Δx

Figure 1.13 One-dimensional FDTD – positions of field components Ez and Hy.

 1.5 ● FDTD updating equations for one-dimensional problems 29

and
nþ1 n1

Hy 2 ðiÞ ¼ Chyh ðiÞ  Hy 2 ðiÞ þ Chyez ðiÞ  Ezn ði þ 1Þ  Ezn ðiÞ þ Chym ðiÞ  Miyn ðiÞ; ð1:43Þ

where

2my ðiÞ  Dtsm

y ðiÞ
Chyh ðiÞ ¼ ;
2my ðiÞ þ Dtsm
y ðiÞ
2Dt
Chyez ðiÞ ¼   ;
2my ðiÞ þ Dtsm
y ðiÞ Dx

2Dt
Chym ðiÞ ¼  :
2my ðiÞ þ Dtsm
y ðiÞ

A MATLAB code is given in Appendix A for a one-dimensional FDTD implementation

based on the updating equations (1.42) and (1.43). The code calculates electric and magnetic
field components generated by a z-directed current sheet, Jz, placed at the center of a pro-
blem space filled with air between two parallel, perfect electric conductor (PEC) plates
extending to infinity in the y and z dimensions. Figure 1.14 shows snapshots of Ez and Hy
within the FDTD computational domain demonstrating the propagation of the fields and
their reflection from the PEC plates at the left and right boundaries.
In the demonstrated problem, the parallel plates are 1 m apart. The one-dimensional
problem space is represented by cells having width Dx ¼ 1 mm. The current sheet at the
center excites a current with 1 A/m2 density and Gaussian waveform

2
t  21010
J z ðtÞ ¼ e 51011 :

The current sheet is excited over a one-cell cross-section (i.e., 1 mm wide), which
translates Jz into a surface current density Kz of magnitude 1  103 A/m. The surface
current density causes!
a discontinuity
!
in magnetic field and generates Hy, which satisfies the
boundary condition K ¼ n^  H : Therefore, two waves of Hy are generated on both sides of
the current sheet, each of which has magnitude 5  104 A/m. Since the fields are propa-
gating in free space, the magnitude of the generated electric field is h0  5  104  0.1885
V/m, where h0 is the intrinsic impedance of free space (h0  377).
Examining the listing of the one-dimensional FDTD code in Appendix A, one immedi-
ately notices that the sizes of the arrays associated with Ez and Hy are not the same. The one-
dimensional problem space is divided into nx intervals; therefore, there are nx þ 1 nodes in
the problem space including the left and right boundary nodes. This code is based on the
field positioning scheme given in Figure 1.13, where the Ez components are defined at node
positions and the Hy components are defined at the center positions of the intervals. There-
fore, arrays associated with Ez have the size nx þ 1, whereas those associated with Hy have
the size nx.
Another point that deserves consideration is the application of the boundary conditions.
In this problem the boundaries are PEC; thus, the tangential electric field component (Ez in
this case) vanishes on the PEC surface. Therefore, this condition can be enforced on the
electric field component on the boundaries, Ez(1) and Ez(nx þ 1), as can be seen in the
30 CHAPTER 1 ● Introduction to FDTD

Time step = 100, Time = 0.3 ns

Ez
Hy × h0

0.2

0.1
[V/m]

–0.1

–0.2
0.2
1
0
0.5

–0.2 0
(a) [A/m] x [m]

Time step = 300, Time = 0.9 ns

Ez
Hy × h0

0.2

0.1
[V/m]

–0.1

–0.2
0.2
1
0
0.5

–0.2 0
(b) [A/m] x [m]

Figure 1.14 Snapshots of a one-dimensional FDTD simulation: (a) fields observed after 100
time steps; (b) fields observed after 300 time steps; (c) fields observed after 615
time steps; and (d) fields observed after 700 time steps.
 1.5 ● FDTD updating equations for one-dimensional problems 31

Time step = 615, Time = 1.845 ns

Ez
Hy × h0

0.2

0.1
[V/m]

–0.1

–0.2
0.2
1
0
0.5

–0.2 0
(c) [A/m] x [m]

Time step = 700, Time = 2.1 ns

Ez
Hy × h0

0.2

0.1
[V/m]

–0.1

–0.2
0.2
1
0
0.5

–0.2 0
(d) [A/m] x [m]

Figure 1.14 (Continued )

32 CHAPTER 1 ● Introduction to FDTD

code listing. Observing the changes in the fields from time step 650 to 700 in Figure 1.14
reveals the correct behavior of the incident and reflected waves for a PEC plate. This is
evident by the reversal of the direction of the peak value of Ez after reflection from the PEC
boundaries, which represents a reflection coefficient equal to 1 while the behavior of the
reflected Hy component represents the reflection coefficient that is equal to þ1 as expected.

1.6 Exercises
1.1 Follow the steps presented in Section 1.2 to develop the appropriate approximations for
the first derivative of a function with second-order accuracy using the forward and
backward difference procedure and fourth-order accuracy using the central difference
procedure. Compare your derived expressions with those listed in Table 1.1.
1.2 Update the MATLAB program in Listing 1.1 to regenerate Figure 1.2 using the second-
order accurate forward and backward differences and the fourth-order accurate central
difference approximations for the same function.
1.3 Update the MATLAB program in Listing 1.1 to generate the corresponding figures to
Figure 1.2, but for the second derivatives of the function. Use the approximate
expressions in the right column of Table 1.1 while updating the program.
1.4 The steps of the development of the updating equation for the x component of the
electric field are demonstrated in Section 1.3. Follow the same procedure to show the
steps required to develop the updating equations for the y and z components for the
electric field. Pay sufficient attention to the selection of the indices of each individual
field component in your expressions. The Yee cell presented in Figure 1.5 should help
you in understanding the proper use of these indices.
1.5 Repeat Exercise 1.4, but this time for developing the updating equations for the mag-
netic field components.