IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

1

A Dissipation Theory for Three-Dimensional
FDTD with Application to Stability
Analysis and Subgridding

arXiv:1705.02274v2 [] 22 Jul 2017

Fadime Bekmambetova, Student Member, IEEE, Xinyue Zhang Student Member, IEEE, and Piero Triverio, Senior
Member, IEEE

Abstract—The finite-difference time-domain (FDTD) algorithm
is a popular numerical method for solving electromagnetic
problems. FDTD simulations can suffer from instability due
to the explicit nature of the method. Stability enforcement
can be particularly challenging in scenarios where a setup is
composed of multiple components, such as grids of different
resolution, advanced boundary conditions, reduced-order models,
and lumped elements. We propose a dissipation theory for 3-D
FDTD inspired by the principle of energy conservation. We view
the FDTD update equations for a 3-D region as a dynamical
system, and show under which conditions the system is dissipative. By requiring each component of an FDTD-like scheme to
be dissipative, the stability of the overall coupled scheme follows
by construction. The proposed framework enables the creation
of provably stable schemes in an easy and modular fashion, since
conditions are imposed on the individual components, rather than
on the overall coupled scheme as in existing approaches. With
the proposed framework, we derive a new subgridding scheme
with guaranteed stability, low reflections, support for material
traverse and arbitrary (odd) grid refinement ratio.
Index Terms—dissipation, finite-difference time-domain, stability, subgridding

I. I NTRODUCTION
Stability is a major challenge in explicit numerical
schemes for solving differential equations. The conventional
finite-difference time-domain (FDTD) algorithm for solving
Maxwell’s equations [1] is stable when the iteration time
step is below the Courant-Friedrichs-Lewy (CFL) limit [2].
However, stability enforcement becomes much more difficult
in the case of more advanced FDTD schemes such as, for
example, those including locally-refined grids [3–5], reducedorder models [6], [7], and hybridizations of FDTD with
integral equation methods [8] or circuit simulators [9]. Many
of these advanced schemes can be seen as the coupling of
multiple blocks, such as FDTD grids of different resolution,
reduced-order models, lumped components, boundary conditions, and so on. While the stability of each individual block
may be well understood, analyzing and enforcing the stability
of the overall coupled scheme can be a daunting task.
Manuscript received ...; revised ...
This work was supported in part by the Natural Sciences and Engineering
Research Council of Canada (Discovery grant program) and in part by the
Canada Research Chairs program.
F. Bekmambetova, X. Zhang and P. Triverio are with
the Edward S. Rogers Sr. Department of Electrical and
Computer
Engineering,
University
of
Toronto,
Toronto,
M5S 3G4 Canada (email: fadime.bekmambetova@mail.utoronto.ca,
xinyuezhang.zhang@mail.utoronto.ca, piero.triverio@utoronto.ca).

Several techniques have been proposed for FDTD stability
analysis, such as von Neumann analysis [10], the iteration
method [2] and the energy method [11]. Unfortunately, von
Neumann analysis cannot be applied in presence of inhomogeneous materials. The iteration method [2], [12] analyzes
the eigenvalues of a matrix describing FDTD iterations. Even
though an FDTD setup can often be divided into several
subsystems, such as grids of different resolution, reduced-order
models, interpolation conditions, and other blocks, the iteration
matrix needs to be derived for the entire scheme, which can
lead to lengthy derivations and make stability analysis very
tedious. A similar issue arises in the case of the energy
method [11], where one must write an expression for the stored
energy in the entire domain, and show that the coupled scheme
satisfies the principle of energy conservation.
In this paper, we present a dissipation theory for modular
stability enforcement of complex 3-D FDTD systems, generalizing previous work in two dimensions [13]. The method
is based on the theory of dissipative dynamical systems [14].
First, we write the FDTD update equations for an arbitrary
inhomogeneous region in the form of a discrete-time dynamical system. The tangential magnetic and electric field are seen
as input and output, respectively. Suitable discrete expressions
for the energy stored inside the region, and for the energy
absorbed through the boundaries are proposed and used to
investigate under which conditions the system is dissipative.
This result becomes the basis for a powerful framework to
create new FDTD algorithms with guaranteed stability. By
imposing each subsystem to to be dissipative, the overall
coupled scheme formed by these subsystems is dissipative by
construction, and thus stable. The main virtue of this approach
is that the stability of the overall scheme follows automatically
from conditions imposed on each subsystem individually. This
makes stability analysis simpler and more modular, which is
a remarkable advantage over existing approaches that require
the analysis of the overall coupled scheme. This modularity
also facilitates the generation of new schemes, since proving
a given subsystem (e.g. a sophisticated boundary condition)
to be dissipative allows coupling it to any other dissipative
block, with no need for additional analysis. Finally, we use the
theory to derive a stable subgridding scheme for 3-D FDTD
that works for any odd refinement ratio and naturally supports
traversal of material boundaries by the subgridding interface.
This paper is structured as follows. In Sec. II, we show
how the update equations for an FDTD region can be cast

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

z

1

2

in Fig. 2a, is described using a conventional FDTD update
equation [2]

(1, Ny +1, Nz +1)

H n− 2

∆

x

Top

∆z

West
South

En
1

U n+ 2
x

∆y

y

North

σx 
n+1
Ex i+ 1 ,j,k
∆t 2
ε 2 σ 
n
x
x
= ∆x∆y∆z
Ex i+ 1 ,j,k
−
2
∆t 2
n+ 1
n+ 1
+ ∆x∆zHz i+ 12,j+ 1 ,k − ∆x∆yHy i+ 12,j,k+ 1

∆x∆y∆z

East
Bottom

(Nx +1, Ny +1, 1)

Fig. 1. Illustration of the 3-D FDTD region with the hanging variables
shown in green and the regular electric and magnetic samples in blue and
red, respectively. The secondary grid is shown with the dashed line.

into the form of a dynamical system with suitable inputs and
outputs. In Sec. III, we provide dissipativity conditions for the
system and show their relation to the CFL limit. In Sec. IV,
we describe the modular method for enforcing FDTD stability
using the proposed theory, and in Sec. V we apply this method
to derive a stable subgridding scheme. Numerical examples are
given in Sec. VI.

II. DYNAMICAL S YSTEM F ORMULATION OF A 3-D FDTD
R EGION
Consider the FDTD region shown in Fig. 1, which contains
Nx , Ny , and Nz primary cells in the x, y, and z directions,
respectively. Electric field is sampled on the primary grid
edges at the integer time points, n, and magnetic field is
sampled on the primary grid faces at the time points shifted
by half of a time step, n−(1/2) [2]. For simplicity, we assume
that the region is discretized uniformly with primary cells of
dimensions ∆x×∆y×∆z. As shown in Fig. 1, we use cardinal
directions in the xy-plane, with +y being the North. “Top” and
“Bottom” denote the +z and −z directions, respectively.
In addition to the regular electric and magnetic field samples, denoted by E n and H n−(1/2) , respectively, we define
the so-called hanging variables [15], U n+(1/2) , which are
magnetic field samples collocated with the electric fields on
the boundary of the region. The hanging variables, which are
tangential to the boundary and perpendicular to the corresponding electric field samples, are used to define the power
supplied to the region and to facilitate the derivation of a selfcontained dynamical model for the region based on FDTD
update equations.
In this section, we take the FDTD equations for the regular
field samples, E n and H n−(1/2) , and cast them in the form
of a dynamical system with the hanging variables, U n+(1/2) ,
as its inputs. The output consists of electric field samples on
the boundary of the region.

A. Equations at Each Node
The temporal evolution of any electric field sample strictly
inside the region, such as the Ex |ni+(1/2),j,k sample shown

ε

x

+

2

−

2

n+ 1
∆x∆zHz i+ 12,j− 1 ,k
2
2

2

+

2

n+ 1
∆x∆yHy i+ 12,j,k− 1
2
2

, (1)

where σx and εx are electrical conductivity and permittivity
values, respectively, on the primary edge where the electric
field is sampled. For clarity, we do not show the dependence
of material properties on the location, although the proposed
developments are valid in the most general case where all
material properties are inhomogeneous. The iteration time step
is denoted by ∆t. Although ∆x, the length of the edge where
the sample is located, can be canceled in (1), we keep it
for later derivations. The product ∆x∆y∆z is the volume
associated with Ex |ni+(1/2),j,k , as illustrated in Fig. 3. We refer
to the internal electric fields such as Ex |ni+(1/2),j,k as electric
fields of Type 1.
A conventional FDTD equation for the electric fields on
the region’s faces, such as the South Ex |ni+(1/2),1,k sample in
Fig. 2b, would involve some magnetic field samples that are
beyond the considered region. The use of the magnetic fields
outside the region would require assumptions on the nature of
the subsystems connected to the FDTD grid. Instead, we write
a modified equation for the South electric field samples using
the hanging variables on the South boundary
∆x

ε
σx 
∆y
n+1
x
∆z
+
Ex i+ 1 ,1,k
2
2
∆t 2
ε
∆y
σx 
n
x
= ∆x
∆z
−
Ex i+ 1 ,1,k
2
2
∆t 2
∆y
n+ 12
n+ 12
Hy i+ 1 ,1,k+ 1
+ ∆x∆zHz i+ 1 , 3 ,k − ∆x
2 2
2
2
2
1
∆y
n+ 2
n+ 21
− ∆x∆zUz i+ 1 ,1,k + ∆x
Hy i+ 1 ,1,k− 1 . (2)
2
2
2
2

Equation (2) is obtained from an FDTD-like approximation of
Maxwell-Ampère law on a half-cell containing Ex |ni+(1/2),1,k .
Using (2) instead of a conventional FDTD equation, we obtain
a self-contained FDTD-like model for the fields in the region
that does not involve field samples beyond its boundaries. This
feature is crucial for investigating under which conditions the
region is dissipative. It is also needed to obtain an FDTD
model that can be connected to other subsystems, such as a
grid of different resolution or a reduced model, where some
magnetic field samples may not be available for a conventional
FDTD equation. Samples like Ex |ni+(1/2),1,k are classified as
Type 2, which are the electric field samples on the faces of
the boundary.
In a similar way, we write a modified FDTD equation for the
samples located at the edges where two faces of the boundary
come together. Those equations involve two hanging variables

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

3

(a)

(b)

(c)

z
Hy

Hz

n+ 1
2

n+ 1
2

n
Ex j,k

j− 1
,k
2

Hy

Hy

x

1
j,k+ 2

Hz

z

z
y

n+ 1
2

n
Ex 1,k

j+ 1
,k
2

n+ 1
2

Uz

Hy

j,k− 1
2

y

n+ 1
2

1,k+ 1
2

x

n+ 1
2

n+ 1
2

1,k

Hz

y
x
Hy

3 ,k
2

n
Ex 1,1

n+ 1
2

1,k− 1
2

Uz

1
n+ 2

1, 3
2

n+ 1
2
1,1

n+ 1
Uy 1,1 2

Hz

n+ 1
2
3 ,1
2

Fig. 2. Variables involved in the recurrence relations for (a) internal Exn samples (Type 1), (b) Exn samples of Type 2 on the South face, and (c) the Exn
samples of Type 3 on the Bottom-South edge of the boundary. The common subscript “i+(1/2)” is omitted for clarity.

z

z

where µx is the magnetic permeability at the sampling location. Equation for the boundary magnetic field sample
n−(1/2)
Hx |1,j+(1/2),k+(1/2) , which is located at a secondary edge of
length ∆x/2 normal to the West boundary, involves a factor
of ∆x/2, as opposed to ∆x

Exn
1
n− 2

∆z/2

Hx

∆z
∆x

y

y

∆x/2
∆x
∆x/2

∆z/2

x

x

Fig. 3. Illustration of the volumes (shown in dashed lines), associated with
each electric (left panel) and magnetic (right panel) field sample in (1)–(5).

for each sample. For instance, in the case of the Bottom-South
sample shown in Fig. 2c, the recurrence relation reads
∆x

∆y ∆z  εx σx 
n+1
+
Ex i+ 1 ,1,1
2
2 2 ∆t 2
∆y ∆z  εx σx 
n
= ∆x
−
Ex i+ 1 ,1,1
2
2 2 ∆t 2
∆z
∆y
n+ 21
n+ 21
+ ∆x
Hz i+ 1 , 3 ,1 − ∆x
Hy i+ 1 ,1, 3
2 2
2
2
2
2
∆y
∆z
n+ 12
n+ 12
− ∆x
Uz i+ 1 ,1,1 + ∆x
Uy i+ 1 ,1,1 . (3)
2
2
2
2

∆x
∆x
µx
µx
n+ 1
n− 1
∆y∆z Hx 1,j+2 1 ,k+ 1 =
∆y∆z Hx 1,j+2 1 ,k+ 1
2
2
2
2
2
∆t
2
∆t
∆x
∆x
n
n
−
∆zEz 1,j+1,k+ 1 +
∆yEy 1,j+ 1 ,k+1
2
2
2
2
∆x
∆x
n
n
+
∆zEz 1,j,k+ 1 −
∆yEy 1,j+ 1 ,k , (5)
2
2
2
2
and similarly for the boundary faces on the East side. Recurn−(1/2)
n−(1/2)
are obtained in a
and Hz
rence relations for Hy
similar fashion.
The coefficients ∆x∆y∆z in (4) and (∆x/2)∆y∆z in (5)
correspond to the dimensions of the volume associated with
n−(1/2)
sample, as shown in Fig. 3.
each Hx
B. Compact Matrix Form
In order to shorten the notation, we write the equations (1)–
(5) for all field samples in a matrix form. Collecting all
n−(1/2)
recurrence relations for Hx
samples (4) and (5) into
matrix equations, we obtain

n+(1/2)

In (3), the y- and z-directed hanging variables, Uy |i+(1/2),1,1
n+(1/2)
Uz |i+(1/2),1,1 ,

and
respectively, are used to complete the line
integral of magnetic field around Ex |ni+(1/2),1,1 . The sample
Ex |ni+(1/2),1,1 is an example of a Type 3 electric field, which is
a field shared between two faces of the region’s boundary. The
three types of Eyn and Ezn fields and their recurrence relations
are defined analogously to Exn and to (1)–(3).
n−(1/2)
The recurrence relation for the magnetic Hx
samples
that are strictly inside the region is the conventional update
FDTD equation, which reads
µx
µx
n+ 1
n− 1
Hx i,j+2 1 ,k+ 1 = ∆x∆y∆z Hx i,j+2 1 ,k+ 1
2
2
2
2
∆t
∆t
n
n
− ∆x∆zEz i,j+1,k+ 1 + ∆x∆yEy i,j+ 1 ,k+1

∆x∆y∆z

2

+ ∆x∆zEz

n
i,j,k+ 21

2

− ∆x∆yEy

n
i,j+ 12 ,k

, (4)

Dlx0 DAx

Dµx n+ 12
Dµx n− 12
Hx
= Dlx0 DAx
Hx
∆t
∆t
n
+ Dlx0 Gzy Dly Ey − Dlx0 Gyz Dlz Enz , (6)
n−(1/2)

n−(1/2)

where vectors Hx
and Eny contain all Hx
and
n
Ey samples in the region, respectively. The coefficient matrix
Dµx is a diagonal matrix containing the magnetic permeability
values on the edges where the corresponding elements of
n−(1/2)
Hx
are sampled. Matrix Gzy , which consists of zeros,
+1’s and −1’s, is a discrete z-derivative operator for Eny .
Similarly, Gyz is a discrete y-derivative operator for Enz .
Matrix Dlx0 is a diagonal matrix containing the length of
n−(1/2)
each secondary edge associated with samples in Hx
,
namely ∆x for internal samples and ∆x/2 for the samples
on the East and West boundaries. Diagonal matrix DAx
n−(1/2)
contains the area of primary faces where Hx
is sampled,

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

4

which reduces to ∆y∆zINHx for a uniform discretization,
where Im denotes an m × m identity matrix and NHx is
n−(1/2)
the number of Hx
samples in the region. As a result,
the product Dlx0 DAx contains the volumes shown in the right
n−(1/2)
panel of Fig. 3, associated with each sample in Hx
.
Matrices Dly and Dlz contain the length of the primary edges
where the elements of Eny and Enz are sampled, respectively.
These matrices reduce to ∆yINEy and ∆zINEz for uniform
discretization.
In a similar fashion, we write in matrix form the recurrence
n−(1/2)
n−(1/2)
relations for Hy
and Hz
as
Dly0 DAy

Dlz0 DAz

Dµy n− 12
Dµy n+ 21
= Dly0 DAy
Hy
Hy
∆t
∆t
− Dly0 Gzx Dlx Enx + Dly0 Gxz Dlz Enz , (7)
Dµz n− 12
Dµz n+ 21
= Dlz0 DAz
Hz
Hz
∆t
∆t
+ Dlz0 Gyx Dlx Enx − Dlz0 Gxy Dly Eny , (8)

with the coefficients defined similarly to those in (6).
From (1), (2), and (3), the recurrence relation for the Exn
samples can be written in matrix form as




Dεx Dσx
Dεx Dσx
0
Dlx DA0x
+
En+1
=
D
D
−
Enx
l
A
x
x
x
∆t
2
∆t
2
n+ 12

+ Dlx GTzx Dly0 Hy

n+ 1

n+ 12

− Dlx GTyx Dlz0 Hz

n+ 1

+ Dlx Bx,B Dly0 ,B Uy,B2 + Dlx Bx,N Dlz0 ,N Uz,N2
n+ 1

n+ 1

− Dlx Bx,T Dly0 ,T Uy,T 2 − Dlx Bx,S Dlz0 ,S Uz,S 2 , (9)
where Dεx and Dσx are diagonal matrices with the values
of permittivity and electric conductivity, respectively, on the
x-directed primary edges. The diagonal matrix DA0x contains
the area of the secondary faces where Exn is sampled, which is
equal to ∆y∆z for Type 1 nodes, ∆y∆z/2 for Type 2 nodes,
and ∆y∆z/4 for Type 3 nodes. Matrix Bx,B consists of 0’s
n+(1/2)
and 1’s that extract from the vector Uy,B
the hanging
variables corresponding to a Type 2 or Type 3 sample in Enx .
Matrices Bx,T , Bx,S , and Bx,N for the Top, South, and North
faces serve an analogous purpose. In a similar manner, we
write the matrix equations for Eny and Enz .
Equations (6)–(8) describing the recurrence relation for all
regular magnetic field variables can be written as
Dl0 DA

Dµ n+ 1
Dµ n− 1
H 2 = Dl0 DA
H 2 − Dl0 CDl En , (10)
∆t
∆t

where vectors En and Hn−(1/2) collect all electric and regular
magnetic field variables, respectively
 n− 1 
 n
Hx 2
Ex
1
1

n
n−
n
2
E = Ey  , H 2 = Hn−
(11)
.
y
n− 12
Enz
H
z

The diagonal matrices Dl0 and DA contain the secondary
edge lengths and primary face areas associated with samples
in Hn−(1/2) . Matrix Dµ contains permeability on the edges

where samples in Hn−(1/2) are located. Coefficient matrix C
is a discrete curl operator


0
−Gzy
Gyz
0
−Gxz 
C =  Gzx
(12)
−Gyx
Gxy
0
and Dl is a diagonal matrix containing the length of the
primary edges corresponding to the samples in En .
Similarly, we obtain the following recurrence relation for
the electric samples in the region by collecting (9) and similar
equations for Eny and Enz

Dl DA0

Dε Dσ
+
∆t
2



En+1 = Dl DA0



Dε Dσ
−
∆t
2

1



En

1

+ Dl CT Dl0 Hn+ 2 + Dl QPDlU0 Un+ 2 , (13)
where DA0 contains the secondary face areas associated with
the samples in En . Matrices Dε , and Dσ contain permittivity
and electrical conductivity on the edges that correspond to the
samples in En . Matrix Q contains 1’s and 0’s that extract the
hanging variables from vector Un+(1/2) . For Type 1 samples
in En , the corresponding rows in Q contain only zeros. For a
Type 2 sample in En , the row of Q contains a single 1 in the
column corresponding to the hanging variable that is involved
in an equation such as (2) for that electric field sample. For a
Type 3 sample, the row of Q contains two 1’s in the columns
corresponding to the two hanging variables in the equation
like (3). Matrix P is a diagonal matrix that selects the signs
for the hanging variables according to (9). The product QP
has the following structure


Bx 0
0
(14)
QP =  0 By 0  ,
0
0 Bz
where


Bx = Bx,B Bx,N −Bx,T −Bx,S ,


By = By,W By,T −By,E −By,B ,


Bz = Bz,S Bz,E −Bz,N −Bz,W .

(15a)
(15b)
(15c)

Matrix DlU0 in (13) is a diagonal matrix containing the length
of the edges where the hanging variables are sampled, namely
∆x, ∆y, or ∆z for the hanging variables collocated with
Type 2 electric fields and ∆x/2, ∆y/2, or ∆z/2 for the
hanging variables associated with the electric fields of Type 3.
C. Dynamical System Formulation
Matrix equations (10) and (13) can be written in the form
of a dynamical system as follows
(R + F)xn+1

=

yn

=

1

(R − F)xn + Bun+ 2 ,
LT xn ,

(16a)
(16b)

where the state vector xn contains the regular FDTD samples
in the region
 n 
E
n
x =
.
(17)
1
Hn− 2

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

5

The hanging variables U n+(1/2) are not included in the state
vector, and instead are regarded as the input of the FDTD
region
 n+ 1 
UEx 2


1
n+ 12 
un+ 2 = 
(18)
UEy  ,
n+ 12
UEz
n+(1/2)
UEx

where vector
collects the variables associated with
n+(1/2)
the Exn samples, namely the Uy
variables tangential to
n+(1/2)
the Top and Bottom faces and Uz
variables for the North
n+(1/2)
n+(1/2)
and South faces. Vectors UEy
and UEz
are defined
similarly. The output vector yn is defined as
 n 
Y Ex
n 
y n =  YE
,
(19)
y
n
YE
z
n
YE
x

Exn

where
contains the
samples of Types 2 and 3
n+(1/2)
associated with the hanging variables in UEx
. Type 2
n
samples are included in YEx once, while each Type 3 sample
n
twice, at the positions of the two coris included in YE
x
n+(1/2)
n
responding hanging variables in UEx
. Vectors YE
and
y
n
n
YEz are defined similarly to YEx .
Matrices R, F, B, and LT in (16a)–(16b) are given by


ε
Dl DA0 D
− 12 Dl CT Dl0
∆t
R=
,
(20a)
D
− 1 Dl0 CDl
Dl0 DA ∆tµ
 2

D D 0 Dσ − 12 Dl CT Dl0
F= 1l A 2
,
(20b)
0
0
2 Dl CDl


Dl QPDlU0
B=
,
(20c)
0


LT = QT 0 .
(20d)

energy absorbed by the region through its boundaries between
time points n and n+1.
We define the storage function and the supply rate as
∆t n T
(x ) Rxn ,
2
T


yn + yn+1
1
n
n+ 21
= ∆t
Sun+ 2 .
s y ,u
2
The diagonal matrix S in (23) is given by
E(xn ) =

S = DlY DlU0 P ,

(22)
(23)

(24)

where DlY contains the lengths of the primary edges on which
the electric samples in yn are taken. Expressions (22) and (23)
generalize those from the 2-D case [13].
A. Physical Meaning of the Supply Rate and Storage Function
By substituting (17) and (20a) into (22), we can write the
storage function as
1
E(xn ) = (En )T Dl DA0 Dε En
2


∆t  n− 1 T
Dµ n− 1
+
H 2
Dl0 DA
H 2 − Dl0 CDl En , (25)
2
∆t
which, using (10) to simplify the term inside the square
brackets, reduces to

E(0) = 0 ,
(21b)



1
1
E xn+1 − E(xn ) ≤ s yn , un+ 2 , ∀un+ 2 ∀n . (21c)

1 n T
(E ) Dl DA0 Dε En
2
1
1  n− 1 T
H 2
Dl0 DA Dµ Hn+ 2 . (26)
+
2
This expression for the stored energy has been proposed
in [11] for FDTD/FIT. Equation (26) reveals the physical
meaning of the storage function (22) as a discrete analogy
of
˚

1
εx Ex2 + εy Ey2 + εz Ez2 + µx Hx2 + µy Hy2 + µz Hz2 dV ,
2
V
(27)
which is the energy stored in electric and magnetic fields inside
a given volume.
For the supply rate, we can substitute (24) into (23),
obtaining
T
yn + yn+1
1
n
n+ 21
s(y , u
) = ∆t
DlY DlU0 Pun+ 2 . (28)
2
The product DlY DlU0 is a diagonal matrix containing the area
of boundary faces associated with each hanging variable and
the corresponding electric field sample. Diagonal matrix P has
a +1 on the faces where the Poynting vector points into the
region and a −1 otherwise. As a result, (28) is a discrete
version of the Poynting integral over the region’s boundary
from time sample n to n+1
ˆ (n+1)∆t ‹
~ ×H
~ · n̂ dA dt ,
E
(29)

The storage function E(xn ) quantifies the energy stored
in

the FDTD region and the supply rate s yn , un+(1/2) is the

where n̂ is the inward unit normal vector. Therefore, supply
rate (23) is the total energy that enters the region from its
boundaries between time n and time n+1.

The definitions given in this section allow us to write the
FDTD equations for a 3-D region in the compact matrix
form (16a)–(16b). Remarkably, equations (16a)–(16b) have the
same structure as in the 2-D case [13]. This fact allows us to
generalize previous results valid in two dimensions to the most
general 3-D case.
III. D ISSIPATIVITY OF A 3-D FDTD R EGION
In this section, we derive dissipativity conditions for the 3D FDTD system (16a)–(16b) defined in Sec. II and explain
their physical significance.
Definition 1. According to the theory of dissipative systems [16], the discrete-time dynamical system (16a)–(16b) is
dissipative with the supply rate s yn , un+(1/2) if we can find
a function E that satisfies
E(xn ) ≥ 0 ,

∀xn ,

(21a)

E(xn ) =

n∆t

S

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

6

B. Dissipativity Conditions for a 3-D FDTD Region
By substituting (22) and (23) into (21c), we arrive at the
following theorem, which gives some simple conditions for
the dissipativity of (16a)–(16b).
Theorem 1. A system in the form (16a)–(16b) is dissipative
with E(xn ) in (22) as storage function and s yn , un+(1/2)
in (23) as supply rate if
R = RT > 0 ,

(30a)

T

F + F ≥ 0,

(30b)

LS = B .

(30c)

Proof. See [13].
The physical meaning of the first condition (30a) can be
better understood by applying the Schur complement [17] to
transform (30a) into

Dl DA0 Dε > 0
∆t
. (31)
D 0 D Dµ − ∆t D 0 CD D−1 D−1 CT D 0 > 0
0
l
l
A ∆t
l
l
ε
A
4
The first matrix inequality in (31) is true because lengths,
areas, and permittivities in the region are positive. Applying
a congruence with matrix 2(∆tDl0 DA Dµ )−(1/2) to the left
hand side of the second inequality in (31), we obtain an
equivalent condition to (31):
4
I N − ΣT Σ > 0 ,
∆t2 H
where

1

−1

−1

1

−1

(32)
−1

Σ = Dl2 Dε 2 DA02 CT Dl20 Dµ 2 DA 2 .

(33)

Condition (32) holds when

∆t < min
k

2
sk


,

(34)

where sk are the non-zero singular values of Σ. Thus, condition (30a) can be seen as a generalized CFL limit, since
it is applicable to the most general case of a lossy region
with nonuniform material properties. Moreover, with a strategy
similar to [11] we can show that a sufficient condition for (30a)
to hold is that the classical FDTD CFL limit [2] is met
√
µε
q
,
(35)
∆t <
1
1
1
∆x2 + ∆y 2 + ∆z 2
where ε is the smallest primary edge permittivity in the region
and µ is the smallest secondary edge permeability. When
the CFL limit is violated, the energy stored in some cells
is no longer bounded below by zero, which can make them
capable of supplying unlimited energy to the rest of the system,
potentially causing instability.
Condition (30b) expands into


Dl DA0 Dσ 0
F + FT =
≥ 0,
(36)
0
0
which holds when the conductivity on each edge is nonnegative.

The third condition (30c) is always true. In order to see this,
we expand LS as


QDlY PDlU0
LS =
.
(37)
0
Right-multiplication by DlY has the same effect on Q as leftmultiplication by Dl . Hence, QDlY = Dl Q, and as a result
LS = B.
In summary, Theorem (1) shows that the FDTD system (16a)–(16b) associated to an arbitrary region is dissipative
if
1) all conductivities are non-negative, as one may obviously
expect;
2) the CFL limit (30a) is respected.
If these two conditions are satisfied, the FDTD region can be
arbitrarily interconnected with any other dissipative subsystem
without violating stability.
IV. S YSTEMATIC M ETHOD FOR S TABILITY E NFORCEMENT
The developments in the previous sections provide a powerful approach to construct both simple and advanced FDTD
schemes with guaranteed stability. The method is general
since it is applicable to complicated setups involving multiple
subsystems, such as FDTD grids with different resolution,
multiple boundary conditions, circuit models, or reduced-order
models. The approach involves the following steps:
1) Hanging variables, which are seen as inputs, are defined
at the boundary of each block. The corresponding electric field samples are regarded as outputs, as done in
Sec. II for a 3D-FDTD region.
2) If subsystems cannot be connected directly, for example
because of different grid resolution, an interpolation rule
is created, and viewed as a separate subsystem.
3) Stability is enforced by ensuring that all subsystems are
dissipative. For dynamical systems of the form (16a)–
(16b) this can be done using (30a)–(30c). Since the
connection of dissipative systems is also dissipative [18],
the resulting scheme is stable by construction.
4) The most restrictive time step is taken to ensure that all
subsystems are dissipative.
The proposed approach significantly simplifies the creation
of new FDTD schemes with guaranteed stability, since dissipativity conditions can be imposed on individual subsystems. In
contrast, existing techniques to analyze and enforce stability,
such as the energy [11] and iteration [2], [12] methods, require
the analysis of the entire scheme, which can be a formidable
task. A remarkable feature of the proposed approach is its
modularity. Given a set of subsystems that have been proven
to be dissipative, any of their combination is guaranteed to be
stable. This feature is hoped to accelerate scientific research in
the FDTD area, since it allows researchers to create new stable
schemes without having to redo stability analysis for every
change. The proposed modular framework is relevant also for
commercial FDTD solvers, where users want to be able to
combine available FDTD subsystems in the way most suitable
to model their problem, while having an absolute guarantee of
stability.

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

y

7

z

PEC

z

Ezn

∆z

fine grid

∆ẑ

PEC

7

8

9

n+ 1
Ux 2

interpolation rule

n
Êz6

4

5

6

1

2

3

n+ 1
2

Ûx6

coarse grid
x

z

PEC

x
(i, j, k)

Fig. 4. Illustration of the subgridding scenario. Left: xy cross section of the
primary grid. Right: subsystem interpretation.

It should be noted that, although we use a matrix formulation to investigate the dissipativity of FDTD equations, the
final scheme can be implemented with scalar FDTD equations
for optimal efficiency.

Fig. 5. Coarse and fine sides of the subgridding interface in Sec. V.

n+(1/2)

and the fine hanging variables Ûx
x direction
n+ 12

Ûx1

n+ 1
Ûx4 2
n+ 1
Ûx7 2

V. A PPLICATION TO S TABLE 3-D FDTD S UBGRIDDING
As a demonstration of the stability framework in Sec. IV,
we derive a stable subgridding algorithm for 3-D FDTD. In
the proposed algorithm, a coarse grid is refined in selected
regions to better resolve fine geometrical details. We denote
the refinement factors as rx , ry , and rz in the x, y, and z
directions, respectively. The refinement factors are odd integers
greater than one.
For simplicity, we present the proposed method by considering the scenario in Fig. 4, where the interface between the
coarse and fine grids is planar. In the Appendix, we discuss
how corners can be similarly treated.
In order to develop a stable subgridding scheme with the
method in Sec. IV, we view a subgridding algorithm as
a connection of four subsystems: boundary conditions, the
coarse grid, the fine grid, and the interpolation rule that relates
fields at the interface of the two grids, as depicted in Fig. 4.
A. Interpolation Conditions
We first describe the interpolation rule established between
the coarse and fine grid fields, which is the core of every
subgridding algorithm. In the scenario of Fig. 4, the interpolation rule has to properly relate the fields tangential to the
North boundary of the coarse grid to the fields tangential to the
South boundary of the fine grid. We discuss the interpolation
n+(1/2)
rule in detail for the Ezn -Ux
pairs. The derivations for
n+(1/2)
Exn -Uz
pairs can be done analogously.
As shown in the left panel of Fig. 5, we consider the portion
of the interface surrounding one coarse electric field sample
(denoted as Ezn ) and the corresponding hanging variable
n+(1/2)
(denoted as Ux
). The fine samples that fall in the same
region in the case of rx = rz = 3 are shown in the right panel
of Fig. 5.
The fine samples are numbered form 1 to 9, as the nodes
where they are sampled. The fine electric fields Êzn are set
equal in the z direction
n
n
n
Êz1
= Êz4
= Êz7
,
n
Êz2
n
Êz3

=
=

n
Êz5
n
Êz6

=
=

n
Êz8
n
Êz9

,
,

∀n ,

∀n ,

∀n ,

(38a)
(38b)
(38c)

x
∆x̂

∆x/2

n+ 12

= Ûx2
=
=

n+ 1
Ûx5 2
n+ 1
Ûx8 2

n+ 12

= Ûx3
=
=

n+ 1
Ûx6 2
n+ 1
Ûx9 2

are set equal in the

,

∀n ,

(39a)

,

∀n ,

(39b)

,

∀n ,

(39c)

where the “ˆ” notation refers to the variables of the refined
region.
n+(1/2)
are forced to be
The coarse samples Ezn and Ux
equal to the average of the fine samples through the following
interpolation rules
Ezn =
n+ 12

Ux

=

1
T
n
rx Trx Êz ,
n+ 21
1
T
rz Trz Ûx

where
 n
Êz1
Ê n 
n
Êz =  z2  ,
n
Êz3

n+ 21

Ûx

,

∀n ,

(40)

∀n ,

(41)

 n+ 1 
Ûx1 2
 n+ 1 
2
=
Ûx4 
n+ 1
Ûx7 2

(42)

and Tm is an m × 1 matrix of ones.
Interpolation rules (38), (39), (40), and (41) draw inspiration
from reciprocity theory for stable subgridding [19]. In the z
direction, the pair of rules (38) and (41) is similar to [19],
except the magnetic fields are related directly at the interface,
as opposed to a plane that is offset from the interface as in [19].
In the x direction, the pair (39) and (40) is analogous, except
the magnetic and electric interpolation rules are switched.
B. Stability Proof
The stability of the proposed subgridding algorithm can be
proved with the dissipation theory in Sec. IV. As shown in
Fig. 4, the proposed scheme can be seen as the connection of
various subsystems, which must all be dissipative in order to
ensure stability. As proven in Sec. III, the coarse and fine grids
are dissipative under their respective CFL limits. The perfect
electric conductor (PEC) boundary condition is lossless, since
imposing a null electric field means that no energy will be
exchanged at any time.
What is left to investigate is whether the interpolation rule
is dissipative or not. We show that the interpolation rule in
Sec. V-A is actually lossless for any time step, by proving
that the supply rate (28) of the interpolation rule is zero. The

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

supply rate of the interpolation rule is the summation of terms,
each one associated to one hanging variable on its boundaries.
The terms associated with the x-directed hanging variables
n+(1/2)
n+(1/2)
Ux
and Ûx
read
E n + Ezn+1 n+ 21
n+ 1
Ux
(43)
sUx 2 = +∆t∆x∆z z
2
T

n
n+1 

∆x ∆z Trz⊗ Êz + Trz⊗ Êz
n+ 1
− ∆t
Ûx 2 ⊗ Trx .
rx rz
2
Symbol “⊗” in (43) denotes the Kronecker product [20]. The
supply rate term (43) accounts for the contribution of the
n+(1/2)
n+(1/2)
Ux
and Ûx
samples in Fig. 5 and the discussion
for z-directed hanging variables is analogous. Using the properties of the Kronecker product [20], we can rewrite (43) as
n+ 1

Ezn + Ezn+1 n+ 21
Ux
2

T
n
n+1
Ê
+
Ê
1
z
z
∆x ∆z T n+ 2
− ∆t
Trx . (44)
T Ûx
rx rz rz
2

sUx 2 = +∆t∆x∆z

Substituting (40) and (41) into (44), we see that the terms
corresponding to the coarse and fine grids cancel each other
n+ 1

Ezn + Ezn+1 n+ 21
Ux
2

n+1
n
∆x ∆z
n+ 1 rx Ez + rx Ez
= 0 . (45)
− ∆t
rz Ux 2
rx rz
2

sUx 2 = +∆t∆x∆z

Iterating the argument for all other coarse hanging variables
and the corresponding groups of rx rz hanging variables on
the fine side, we conclude that the total supply rate to the
interpolation rule is zero. Hence, the energy leaving the
fine grid during a time step is equal to the energy that the
interpolation rule supplies to the coarse grid. This makes the
interpolation rule a lossless, and hence dissipative, subsystem.
Since the coarse and fine grids are dissipative under their own
CFL limits, the entire scheme is stable under the CFL limit
of the fine grid, which is the most restrictive CFL limit.
C. Practical Implementation
In this section we discuss how the proposed subgridding
algorithm can be implemented. For updating all field samples
strictly inside the two grids, one can use conventional FDTD
equations. The update equation for the field samples on the
refinement interface is instead derived from the interpolation
rules (38), (39), (40), and (41).
As discussed in Sec. II-A, each electric field sample Ezn on
the North boundary of the coarse grid satisfies
∆y  εz σz  n+1
+
Ez
2 ∆t 2

∆y εz σz  n
n+ 1
= ∆x
−
Ez − ∆xUx 2
2 ∆t 2
∆y
∆y
n+ 1
n+ 1
n+ 1
−
Hy i− 12 + ∆xHx j− 12 +
Hy i+ 12 , (46)
2
2
2
2
2
where only those subscripts that are different from i, j, and
k +(1/2) are explicitly shown for clarity. Each electric field
∆x

8

sample at a fine node m̂ ∈ {1, 2, . . . , 9} with coordinates
(ı̂, ̂ = 1, k̂+(1/2)) satisfies

ε̂zm̂ σ̂zm̂
+
Êzn+1
m̂
∆t
2


∆x n+ 12
∆x ∆y ε̂zm̂ σ̂zm̂
−
Êznm̂ +
Û
=
rx 2ry ∆t
2
rx xm̂
∆x
∆y
∆y
n+ 1
n+ 1
n+ 1
+
Ĥym̂ ı̂+ 12 −
Ĥxm̂ ̂+ 12 −
Ĥym̂ ı̂− 12 , (47)
2
2
2
2ry
rx
2ry

∆x ∆y
rx 2ry



n+(1/2)

where a sample Ĥym̂ ı̂+(1/2) has the same coordinates as the
node m̂, except for the x coordinate, which is ı̂+(1/2). Similar
notation is used for the other samples in (47).
Using (38a) and averaging (47) over the k̂ indexes in the
shaded area of Fig. 5, we obtain the following equation for
n
the fine sample Êz1
X  ε̂zm̂ σ̂zm̂ 
n+1
Êz1
+
∆t
2
m̂=1,4,7


∆x ∆y 1 X
ε̂zm̂ σ̂zm̂
n
=
Êz1
−
rx 2ry rz
∆t
2
m̂=1,4,7
∆x 1 X n+ 12
∆y 1 X
n+ 1
+
Ûxm̂ +
Ĥym̂ ı̂+ 12
2
rx rz
2ry rz
m̂=1,4,7
m̂=1,4,7
∆x 1 X
∆y 1 X
n+ 1
n+ 1
−
Ĥxm̂ ̂+ 12 −
Ĥym̂ ı̂− 12 ,
2
2
rx rz
2ry rz

∆x ∆y 1
rx 2ry rz

m̂=1,4,7

m̂=1,4,7

(48)
n
and similarly for the other samples in Ênz , namely for Êz2
n
and Êz3 . Next, we use (41) to replace the average of the fine
n+(1/2)
hanging variables in (48) with Ux

X  ε̂zm̂ σ̂zm̂ 
n+1
+
Êz1
∆t
2
m̂=1,4,7


ε̂zm̂ σ̂zm̂
∆x n+ 12
∆x ∆y 1 X
n
Êz1
+
=
−
Ux
rx 2ry rz
∆t
2
rx
m̂=1,4,7
∆y 1 X
∆x 1 X
n+ 1
n+ 1
+
Ĥym̂ ı̂+ 12 −
Ĥxm̂ ̂+ 12
2
2
2ry rz
rx rz
m̂=1,4,7
m̂=1,4,7
X
1
∆y 1
n+
−
Ĥym̂ ı̂− 12 . (49)
2
2ry rz

∆x ∆y 1
rx 2ry rz

m̂=1,4,7

Because of (39), the same substitution can be done in the
n
n
and Êz3
.
equations for Êz2
n
n
in matrix
Writing (49) and the equations for Êz2
and Êz3
n
form, we obtain the following relation for the Ê vector
∆x ∆y
rx 2ry

D̂εz D̂σz
+
∆t
2

∆x ∆y
=
rx 2ry
+

!
Ên+1
z

D̂εz D̂σz
−
∆t
2

!
Ênz +

∆x
n+ 1
Trx Ux 2
rx

∆y n+ 12
∆x n+ 12
∆y n+ 21
Ĥy,ı̂+ 1 −
Ĥx,̂+ 1 −
Ĥ
(50)
1 ,
2
2
2ry
rx
2ry y,ı̂− 2

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

9

where
 P

n+ 1
Ĥx,̂+2 1
2

=

Ĥxm̂
m̂=1,4,7
1 
 P Ĥxm̂

rz m̂=2,5,8
 P
Ĥxm̂
m̂=3,6,9

 P
n+ 1
Ĥy,ı̂−2 1
2

=

1
rz

Ĥym̂

m̂=1,4,7
 P

Ĥym̂

m̂=2,5,8
 P
Ĥym̂
m̂=3,6,9

n+ 21
̂+ 12




n+ 12 

̂+ 12 

,

(51a)

,

(51b)



n+ 12 
̂+ 12
n+ 21
ı̂− 21




n+ 12 

ı̂− 21 
n+ 12
ı̂− 21




n+(1/2)

and similarly for Ĥy,ı̂+(1/2) . Diagonal matrix D̂εz in (50)
contains the permittivity of the half-cells adjacent to the fine
grid’s South boundary, averaged in a way similar to the
magnetic fields in (51a)–(51b). The conductivity matrix D̂σz
is defined analogously.
Substituting (40) into (46) and multiplying the result on the
left by Trx /rx , we obtain
∆x ∆y  εz
σz  1rx n+1
+
Ê
rx 2 ∆t
2
rx z

∆x ∆y εz
σz  1rx n
Tr n+ 1
=
−
Êz − ∆x x Ux 2
rx 2 ∆t
2
rx
rx
Trx
∆y Trx
∆y Trx
n+ 12
n+ 12
n+ 1
−
Hy i− 1 + ∆x
Hx j− 1 +
Hy i+ 12 ,
2
2
2
2 rx
rx
2 rx
(52)
where 1rx denotes a square matrix of ones of size rx × rx .
n+(1/2)
, we obtain an
Adding (50) to (52) to eliminate Ux
rx × rx matrix equation



∆x ∆y
∆y
Mεz
Mσz
Ên+1
+
+
z
rx
2
2ry
∆t
2



∆x ∆y
∆y
Mεz
Mσz
=
Ênz
+
−
rx
2
2ry
∆t
2
∆y Trx
Tr
∆y Trx
n+ 1
n+ 1
n+ 1
−
Hy i− 12 + ∆x x Hx j− 12 +
Hy i+ 12
2
2
2
2 rx
rx
2 rx
∆y n+ 12
∆x n+ 21
∆y n+ 21
+
Ĥ
Ĥ
Ĥ
(53)
1 −
1 −
1 ,
2ry y,ı̂+ 2
rx x,̂+ 2
2ry y,ı̂− 2
with the vector of unknowns Ên+1
, where matrices Mεz and
z
Mσz are given by

−1 

∆y
∆y 1rx
∆y
∆y
+
εz +
D̂εz , (54)
Mεz =
2
2ry
2 rx
2ry
−1 


∆y
∆y 1rx
∆y
∆y
Mσz =
+
σz +
D̂σz . (55)
2
2ry
2 rx
2ry
Equation (53) relates the electric field vector Êz to the
surrounding magnetic fields in both the coarse and the fine
grids. We can see that this last step eliminates the hanging
variables on the boundaries of the two grids, which serve as a
temporary means to connect the two grids. Equation (53) is a
recursive relation for the electric field samples at the interface,
that can be used to compute Ên+1
knowing the electric and
z
magnetic samples at the previous times. Although (53) is not

fully explicit in terms of Ên+1
, the matrix in front of this
z
unknown is very small, having size rx × rx . Moreover, it is a
constant. Therefore, this matrix has to be inverted only once,
leading to an explicit update equation for Ên+1
, which will
z
be used for all time steps.
Overall, the proposed subgridding scheme consists of the
following steps.
1) Starting from E n and H n+(1/2) , we update all electric
fields strictly inside fine and coarse regions using standard FDTD equations such as (1).
2) We update the electric fields on the interface between
the two grids using (53) and a similar equation for Enx
samples. Once (53) is solved, the fine electric samples
n+1
n+1
Êz4
, . . . , Êz9
can be found using the equalities (38a)–(38c) and the coarse electric field sample can
be updated using (40).
3) All regular magnetic field samples are updated using
standard FDTD equations, such as (4) and (5), to obtain
the values at n+(3/2).
The inclusion of material properties throughout the derivations rigorously guarantees stability in the case when arbitrary
permittivities and conductivities are assigned to cells on each
side of the subgridding interface. In contrast, the stability
of many existing schemes with material traverse can only
be verified numerically [21], [22] and some schemes exhibit
inaccuracy when objects intersect the interface [23].
VI. N UMERICAL E XAMPLES
The subgridding algorithm from Sec. V was implemented in
MATLAB in order to test the validity of the proposed theory
for FDTD stability and to assess the accuracy of the proposed
subgridding method. Time-consuming portions of the code,
such as the update equations for the fields inside each grid,
were written using vectorized operations.
A. Stability Verification
We verify the stability of the proposed algorithm by simulating the subgridding scenario shown in Fig. 6 for a million
time steps. The setup consists of a 12 cm × 12 cm × 12 cm
cavity with PEC walls. The cavity is discretized with a
uniform coarse mesh with ∆x = ∆y = ∆z = 1 cm. A
4 cm × 4 cm × 4 cm subgridding region with refinement
ratio r = 5 in all directions is placed at the center of the
cavity. In order to test the correctness of stability analysis in
the case when material properties vary, we randomly assign
permittivity values to cells in both grids between the free space
value ε0 and 3ε0 using the rand() function in MATLAB. In
another simulation, we test if losses are properly handled by
introducing, in addition to the varying permittivities, random
conductivities between zero and 5×10−5 S/m. We excite
the cavity using a Gaussian pulse with half-width at halfmaximum bandwidth of 3.53 GHz. Time step is set to 99% of
the fine grid’s CFL limit.
From the results, presented in Fig. 7, it can be seen that
no signs of instability occur during the million time steps
of the simulation, even if nonuniform materials traverse the
coarse-fine interface. This result validates the correctness of
the stability analysis framework presented in this paper.

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

z

10

Case 2

PEC

Case 3

Case 1

z
y

probe

Jz

x

r=5
probe

y

ε0 ∼ 3ε0

x

Jz (0 ∼ 50 µS/m)
Fig. 8. Setup of the material traverse test of Sec. VI-B.

Ez (V/m)

Fig. 6. Setup of the stability verification test in Sec. VI-A.

Ez (V/m)

0.1 Lossless

0
−20

0
Ez (V/m)

−0.1
0.1 Lossy
Ez (V/m)

Lossy dielectric

20

0

Copper

20
0
−20
0

1

2

3

4

5

6

7

Time (ns)

−0.1
0

2

4
6
Timestep

8

10
x 10

5

Fig. 7. The z component of electric field at the probe in the stability
verification test of Sec. VI-A.

B. Material Traverse
In this simulation, we test whether the proposed subgridding
method produces consistent results when the subgridding
interface traverses material boundaries, which is known to
be a challenging scenario for subgridding schemes [23]. We
use the setup shown in Fig. 8, where the subgridding region
is placed in three different ways: in Case 1 and Case 3, a
block of material traverses different sides of the interface.
In the reference Case 2, the block is fully enclosed by the
subgridding boundary. We perform the test for a copper block
and for a lossy dielectric block (εr = 3, σ = 0.05 S/m). A 15cm perfectly matched layer (PML) terminates the domain. The
coarse grid is set to ∆x = ∆y = ∆z =1 cm and the subgridding
region is refined by a factor of r = 3 in each direction. A
Gaussian source with half-width at half-maximum bandwidth
of 1.02 GHz is used. The iteration time step is 99% of the
CFL limit of the fine grid.
Fig. 9 shows the electric field recorded at the probe in
the three cases. The two cases with material traverse are in
very good agreement with the reference simulation, showing
that the proposed algorithm does not lose accuracy when
different materials traverse the subgridding boundary. The test
is successful for both lossy dielectric and for copper.

Fig. 9. The z component of electric field at the probe in material traverse
test of Sec. VI-B for Case 1 ( ∗ ), Case 2 (  ), and Case 3 ( • ).

C. Meta-Screen
In this example, we use subgridding to study the three-slot
meta-screen shown in Fig. 10. The setup for this simulation
is similar to [24], except for differences in discretization and
screen material. In particular, we set the plate conductivity to
1.3×106 S/m and thickness of the plate to 6 mils, based on
the value in [25]. The slot in the center is 13.2 mm by 1.2 mm
and the two satellite slots on the sides are 17 mm by 0.6 mm,
placed at a center-to-center distance of 3 mm away from the
central slot. These dimensions were designed to achieve subwavelength focusing at 10 GHz [25]. We place a line of probes
4.572 mm away from the meta-screen, which corresponds to
approximately 15% of the wavelength at 10 GHz – distance at
which measurements were performed in [25]. The simulation
region was 25.2 mm×19.2024 mm×43.4 mm including PMLs.
As in [24], we terminate the region before the meta-screen
with PEC walls on the East and the West sides, and with
perfect magnetic conductor walls on the Top and the Bottom
sides. A five-cell PML layer is added on all sides in order to
mimic an unbounded domain.
The satellite slots have widths of only 2% of the wavelength, which makes the problem intrinsically multiscale.
Moreover, [24] reports resonance effects that make it difficult
to resolve the structure near the design frequency. We choose
this example to investigate the performance of the proposed
subgridding scheme.
We set the coarse grid to ∆x = ∆z = 0.7 mm, and ∆y =

11

|Ex (10 GHz)| (nV/m)

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

8
6

subgridding
region

4

k̂
~ inc
E

2

~ inc
H

z

probes

y

0
−8

x

−6

−4

−2

0

2

4

6

8

Probe x coordinate (mm)

Fig. 10. Setup of the meta-screen example of Sec. VI-C.

0.1524 mm. In order to properly resolve the meta-screen structure and the fields, we refine a 9.8 mm×2.5908 mm×21 mm
area around the slots with rx = rz = 7 and ry = 3, using the
proposed subgridding method. As a reference, we perform a
conventional FDTD simulation with a nonuniform grid, where
a region of dimensions 9.8 mm×1.0668 mm×21 mm is refined
by rx = rz = 7 and ry = 3 by means of varying ∆x, ∆y, and
∆z. We also simulate the case where the entire structure is
discretized with the coarse mesh.
A modulated Gaussian waveform with 10 GHz central frequency and 8.24 GHz half-width at half-maximum bandwidth
is used to excite the structure with a sheet of x-directed
uniform current density on the South side. The proposed
method and FDTD with nonuniform discretization are run at
∆t = 0.1362 ps, whereas the coarse grid simulation has the
iteration time step of 0.4810 ps. These time steps correspond to
99% of the CFL limit of each mesh. The waveform is recorded
for 5.5 ns, as in [24].
Fig. 11 shows the amplitude of the x component of electric
field at 10 GHz, which was obtained with a Fourier transform
of the time-domain waveforms. The results from subgridding
are in very good agreement with the reference simulation. In
contrast, the coarse resolution was insufficient to perform an
accurate simulation of the meta-screen. This is not surprising,
since any mesh with ∆x larger than 0.3 mm or ∆z larger than
0.1 mm cannot correctly resolve the meta-screen structure, let
alone the fields around the slots.
Simulation times, as well as the number of primary FDTD
cells are shown in Table I. Subgridding almost doubles the
speed-up compared to the nonuniform grid refinement, which
is already much faster than the simulation with uniform grid.
This is a result of the lower number of unknowns associated
with subgridding. The guarantee of stability of the proposed
method allows for the reduction of simulation time without
sacrificing reliability. Moreover, the ability of the method to
handle material traverse allows placing the subgridding region
across the meta-screen.
VII. C ONCLUSION
We proposed a dissipation theory for 3-D FDTD inspired by
the theory of dissipative dynamical systems [14]. The theory
shows that the FDTD equations for the fields in a 3-D region
can be interpreted as a dynamical system that is dissipative
when time step satisfies the CFL limit.

Fig. 11. Magnitude of x-directed electric field along the line of probes in the
),
focusing meta-screen example of Sec. VI-C for uniform coarse grid (
nonuniform grid ( ), and the proposed method ( ).
Table I
C OMPUTATIONAL C OST FOR D IFFERENT D ISCRETIZATION IN THE
M ETA -S CREEN E XAMPLE OF S EC . VI-C.
Method
Coarse
Fine
Nonuniform
Proposed
a

Number of primary cells
281,232
41,341,104
4,065,600
1,323,672

Simulation time
0.18 hours
95.81 hoursa
11.46 hours
5.80 hours

Speed-up
524.8
reference
8.4
16.5

Estimated runtime from 30 time steps.

The proposed dissipation theory provides a powerful and
systematic method to create new FDTD schemes with guaranteed stability. By showing that each component of the
scheme is dissipative (such as grids, boundary conditions, or
embedded lumped components), the overall coupled scheme
is guaranteed to be dissipative, and thus stable. This approach
makes stability analysis simpler, more intuitive, and modular,
since dissipation conditions are imposed on each component
independently and are based on the familiar concept of energy.
With the new framework, we develop a stable 3-D subgridding scheme supporting material traverse. Numerical examples
show good performance of the scheme when accelerating
multiscale problems, as well as its ability to produce accurate and stable simulations when various materials traverse
the subgridding interface. The dissipation framework is not
restricted to subgridding algorithms, and could be applied to
developing more advanced schemes in the future. As a proof
of concept, we have applied the dissipation framework in 2-D
to guarantee the stability of FDTD coupled to reduced-order
models [18].
A PPENDIX
In this section, we handle the case when the fine grid
is surrounded by the coarse grid on all six sides. In this
scenario, the locations where two faces of the boundary come
together require a special treatment. Fig. 12 shows the portion
of the fine grid’s boundary on its South-West edge, which is
discussed in detail in this section.
A. Interpolation Conditions
Unlike in Sec. V, the edge case involves two sets of hanging
variables: those in the y direction on the West side and those in
the x direction on the South side. As in Sec. V, we interpolate

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

12

the fields on the shaded portion of the interface so that the zdirected electric field samples are equal in the z direction
n
n
n
Êz1
= Êz3
= Êz5
,
n
Êz2
n
Êz2
0

=
=

n
Êz4
n
Êz4
0

=
=

n
Êz6
n
Êz6
0

∀n ,

(56a)

∀n ,

,

(56c)

The fine y-directed hanging variables are set equal in the
y direction and the fine x-directed hanging variables are set
equal in the x direction
n+ 12

Ûy1

n+ 21

Ûy3

n+ 21

Ûy5

n+ 12

= Ûy2

n+ 12

= Ûy4

n+ 12

= Ûy6

n+ 21

,

Ûx1

,

Ûx3

n+ 21
n+ 21

,

Ûx5

n+ 1
n+ 1

∀n , (57b)

= Ûx40 2 ,
= Û

n
Ênz(W ) = Êz2
,

n
Ênz(S) = Êz2
0 ,

(58)

where vector Ênz(W ) contains the distinct samples on the West
half of the shaded region in Fig. 12, except for the SouthWest sample at (1, 1, k̂+(1/2)) and similarly for Ênz(S) on
n
the South side. Variable Êz(SW
) is always the single sample
at (1, 1, k̂+(1/2)), which is excluded from Ênz(W ) and Ênz(S) .
In general, Ênz(W ) and Ênz(S) are vectors of size (ry −1)/2 and
(rx −1)/2, respectively, although in the rx = ry = 3 example
they reduce to scalar quantities.
Similarly, the distinct hanging variables in (57a)–(57c) are
n+(1/2)
n+(1/2)
and Ûx
collected into vectors Ûy
h
iT
n+ 1
n+ 21
n+ 1
n+ 1
Ûy 2 = Ûy1
,
(59)
Ûy3 2 Ûy5 2
n+ 12

Ûx

h
n+ 21
= Ûx1

n+ 21

Ûx3

n+ 21

Ûx5

iT

.

(60)

Analogously to (40), we force the coarse electric field
sample Ezn to be equal to the average of Êzn over the West
portion of the shaded area in Fig. 12
∆y n
∆y n
∆y T n
+
,
E =
Ê
T Ê
2 z
2ry z(SW )
ry W z(W )

∀n ,

(61)

with the weighting coefficients chosen to correctly account for
the area of boundary faces associated with each Êzn sample.
The column vector of ones TW has size (ry − 1)/2. A similar
condition is applied on the South portion of the shaded area
in Fig. 12
∆x n
∆x T n
∆x n
E =
Ê
+
T Ê
,
2 z
2rx z(SW )
rx S z(S)

∀n .

(62)

As in the planar case, the coarse hanging variables are
forced to be equal to the average of the corresponding fine
hanging variables
n+ 1
Uy 2
n+ 12

Ux

=
=

1 T n+ 21
T Ûy ,
rz rz
1 T n+ 21
T Ûx ,
rz rz

y

2

60

5

40

3
n
Êz1

1

n+1
West
Û 2
interface y1

20

x

South
n+1
Ûx1 2 interface

Fig. 12. Fine side of the subgridding interface in the Appendix.

∀n . (57c)

,

We collect the distinct electric field samples into vectors as
follows
n
n
Êz(SW
) = Êz1 ,

4

∀n , (57a)

= Ûx20 2 ,
n+ 21
x60

6

(56b)

∀n .

,

z

∀n ,

(63)

∀n .

(64)

B. Stability Proof
In this section, we show that with the interpolation rule in
Sec. A for the edges of the boundary, the interpolation rule remains a lossless subsystem. Consider the portion of the coarsefine interface consisting of the South shaded area in Fig. 12
on the fine side and of the adjacent ∆x/2 × ∆z boundary face
on the coarse side. For this portion of the interface, the net
contribution of the x-directed hanging variables to the supply
rate from the fine and coarse grids to the interpolation rule
boils down to
n+ 1

E n + Ezn+1 n+ 12
∆x
∆z z
Ux
2
2
T

n+1
n
Ê
+
Ê
1
z(S)
z(S)
∆x ∆z T n+ 2
− ∆t
T Ûx
TS
rx rz rz
2
n+1
n
∆x ∆z T n+ 21 Êz(SW ) + Êz(SW )
− ∆t
, (65)
Trz Ûx
2rx rz
2

sUx 2 = +∆t

analogously to (44). Substituting (64), (65) becomes
∆x
E n + Ezn+1 n+ 12
∆z
n+ 1
n+ 1
∆z z
Ux
− ∆t
sUx 2 = +∆t
rz Ux 2
2
2
rz

T
!
n+1
n
n+1
n
∆x Êz(S) + Êz(S)
∆x Êz(SW ) + Êz(SW )
×
TS +
,
rx
2
2rx
2
(66)
which is zero because of (62). The proof for other hanging
variables near the edges of the boundary can be done with a
similar argument. This shows that even in presence of corners,
the interpolation rule is a lossless system, since its supply rate
is zero.

C. Practical Implementation
The update equations for the electric fields in the shaded
area in Fig. 12 are derived using the modified FDTD equations
on each side of the interface and the interpolation conditions (56a)–(56c), (57a)–(57c), (61), (62), (63), and (64). The

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

13

coarse Ezn field at (i, j, k+(1/2)) needs to satisfy the modified
FDTD equation
3

∆x ∆y  εz σz  n+1
Ez
+
2 2 ∆t 2
∆x ∆y  εz σz  n ∆y n+ 21
=3
Ez +
−
Uy
2 2 ∆t 2
2
∆x
n+ 1
n+ 1
n+ 1
−
Hx j+ 12 − ∆yHy i− 12 + ∆xHx j− 12
2
2
2
2
∆y
∆x n+ 12
n+ 12
+
Hy i+ 1 −
Ux
. (67)
2
2
2

Applying a similar averaging procedure as in Sec. V-C, we
obtain the following relations for the fine electric field vectors
defined in (58)


∆x ∆y ε̂z(SW ) σ̂z(SW )
n+1
+
Êz(SW
)
2rx 2ry
∆t
2


∆x ∆y ε̂z(SW ) σ̂z(SW )
∆x n+ 12
n
=
−
Êz(SW
Ux
)+
2rx 2ry
∆t
2
2rx
∆y n+ 12
∆x n+ 12
∆y n+ 21
, (68a)
+
Ĥy(SW ),ı̂+ 1 −
Ĥx(SW ),̂+ 1 −
Uy
2
2
2ry
2rx
2ry
D̂εz(W ) D̂σz(W )
+
∆t
2

∆x ∆y
2rx ry

!

∆x n+ 12
Ĥ
1
2rx x(W ),̂− 2

Ênz(W ) +

∆y n+ 21
∆x n+ 21
∆y
n+ 1
Ĥy(W ),ı̂+ 1 −
Ĥx(W ),̂+ 1 −
TW Uy 2 ,
2
2
ry
2rx
ry
(68b)

∆x ∆y
rx 2ry

D̂εz(S) D̂σz(S)
+
∆t
2

!
Ên+1
z(S)

D̂εz(S) D̂σz(S)
−
∆t
2

∆x ∆y
=
rx 2ry
−

Ên+1
z(W )

D̂εz(W ) D̂σz(W )
−
∆t
2

∆x ∆y
=
2rx ry
+

!

!
Ênz(S) +

∆y n+ 12
Ĥ
1
2ry y(S),ı̂+ 2

∆x n+ 12
∆y n+ 21
∆x
n+ 1
Ĥx(S),̂+ 1 −
Ĥy(S),ı̂− 1 +
TS Ux 2 , (68c)
2
2
rx
2ry
rx

where
n+ 1

2
Ĥx(SW
),̂+ 1

=

2

n+ 1

Ĥx(W2 ),̂− 1

=

2

n+ 1

Ĥx(W2 ),̂+ 1

=

2

n+ 1

2
Ĥx(S),̂+
1

2

=

1
rz
1
rz
1
rz
1
rz

X

Ĥxm̂

n+ 12
̂+ 12

,

(69a)

n+ 12
̂− 12

,

(69b)

n+ 12
̂+ 12

,

(69c)

n+ 12
̂+ 12

,

(69d)

m̂=1,3,5

X

Ĥxm̂

m̂=2,4,6

X

Ĥxm̂

m̂=2,4,6

X

Ĥxm̂

m̂=20 ,40 ,60

where (ı̂ = 1, ̂ = 1, k̂ + (1/2)) is the coordinate of node m̂.
n+ 1
n+(1/2)
In general, vectors Ĥx(W2 ),̂±(1/2) and Ĥy(W ),ı̂+(1/2) have
size (ry −1)/2, with each element corresponding to a samn+(1/2)
n+(1/2)
ple in Êz(W ) . Quantities Ĥy(SW ),ı̂+(1/2) , Ĥy(S),ı̂±(1/2) , and

n+(1/2)

Ĥy(W ),ı̂+(1/2) are defined similarly to (69a)–(69d). Permittivity and conductivity coefficients ε̂z(SW ) , σ̂z(SW ) , D̂εz(W ) ,
D̂σz(W ) , D̂εz(S) , and D̂σz(S) are obtained from averaging
the values in the half-cells adjacent to the boundary over k̂
indexes. Equations (68a)–(68c) serve a purpose similar to (50).
We now combine (67) for the z-directed electric field on the
coarse side, the recurrence relations (68a)–(68c) for the fields
on the fine side, and the interpolation rules (61) and (62), into
a linear system of the form
Azn+1 = bn ,
with an (rx /2 + ry /2 + 3) × 1 vector of unknowns
 n+1 
Ez
 n+1 
Êz(SW ) 


 Ên+1 

z(W ) 

zn+1 = 
 Ên+1  .
 z(S) 


 n+ 12 

 Ux
n+ 12
Uy

(70)

(71)

n+1
n+1
Solving (70) for zn+1 gives values of Ezn+1 , Êz(SW
) , Êz(W ) ,
n+1
n+1
and Ên+1
and Êz5
are
z(S) . With (56a)–(56c), samples Êz3
n+1
n+1
n+1
n+1
obtained from Êz1 , samples Êz4 and Êz6 from Êz2
,
n+1
n+1
n+1
and Êz4
and Êz6
from Êz2
0
0
0 . The size of system (70) is
only (rx /2 + ry /2 + 3) × (rx /2 + ry /2 + 3) and the resulting
overhead is typically small.

R EFERENCES
[1] K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas
Propag., vol. 14, no. 3, pp. 302–307, 1966.
[2] S. D. Gedney, Introduction to the Finite-Difference Time-Domain
(FDTD) Method for Electromagnetics, 1st ed. San Rafael, CA: Morgan
& Claypool Publishers, 2011.
[3] M. Okoniewski, E. Okoniewska, and M. Stuchly, “Three-dimensional
subgridding algorithm for FDTD,” IEEE Trans. Antennas Propag.,
vol. 45, no. 3, pp. 422–429, 1997.
[4] P. Thoma and T. Weiland, “A consistent subgridding scheme for the
finite difference time domain method,” Int. J. Numer. Model El., vol. 9,
no. 5, pp. 359–374, 1996.
[5] K. Xiao, D. J. Pommerenke, and J. L. Drewniak, “A three-dimensional
FDTD subgridding algorithm with separated temporal and spatial interfaces and related stability analysis,” IEEE Trans. Antennas Propag.,
vol. 55, no. 7, pp. 1981–1990, 2007.
[6] Ł. Kulas and M. Mrozowski, “Reduced-order models in FDTD,” IEEE
Microw. Compon. Lett., vol. 11, no. 10, pp. 422–424, 2001.
[7] B. Denecker, F. Olyslager, L. Knockaert, and D. De Zutter, “Generation
of FDTD subcell equations by means of reduced order modeling,” IEEE
Trans. Antennas Propag., vol. 51, no. 8, pp. 1806–1817, 2003.
[8] A. R. Bretones, R. Mittra, and R. G. Martı́n, “A hybrid technique
combining the method of moments in the time domain and FDTD,”
IEEE Microw. Guided Wave Lett., vol. 8, no. 8, pp. 281–283, 1998.
[9] I. Erdin, M. S. Nakhla, and R. Achar, “Circuit analysis of electromagnetic radiation and field coupling effects for networks with embedded
full-wave modules,” IEEE Trans. Electromagn. Compat., vol. 42, no. 4,
2000.
[10] A. Taflove and S. C. Hagness, Computational electrodynamics. Artech
house, 2005.
[11] F. Edelvik, R. Schuhmann, and T. Weiland, “A general stability analysis
of FIT/FDTD applied to lossy dielectrics and lumped elements,” International Journal of Numerical Modelling: Electronic Networks, Devices
and Fields, vol. 17, no. 4, pp. 407–419, 2004.
[12] R. F. Remis, “On the stability of the finite-difference time-domain
method,” J. Comput. Phys., vol. 163, pp. 249–261, 2000.

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

[13] F. Bekmambetova, X. Zhang, and P. Triverio, “A dissipative systems
theory for FDTD with application to stability analysis and subgridding,”
IEEE Trans. Antennas Propag., vol. 65, no. 2, pp. 751–762, 2017.
[14] J. C. Willems, “Dissipative dynamical systems part i: General theory,”
Archive for rational mechanics and analysis, vol. 45, no. 5, pp. 321–351,
1972.
[15] N. V. Venkatarayalu, R. Lee, Y.-B. Gan, and L.-W. Li, “A stable FDTD
subgridding method based on finite element formulation with hanging
variables,” IEEE Trans. Antennas Propag., vol. 55, no. 3, pp. 907–915,
2007.
[16] C. Byrnes and W. Lin, “Losslessness, feedback equivalence, and the
global stabilization of discrete-time nonlinear systems,” IEEE Trans.
Autom. Control, vol. 39, no. 1, pp. 83–98, 1994.
[17] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, Linear Matrix
Inequalities in System and Control Theory, ser. Studies in Applied
Mathematics. SIAM, 1994, vol. 15.
[18] X. Zhang, F. Bekmambetova, and P. Triverio, “A Stable FDTD
Method with Embedded Reduced-Order Models,” IEEE Trans. Antennas
Propag., 2016, (submitted).
[19] L. Kulas and M. Mrozowski, “Reciprocity principle for stable subgridding in the finite difference time domain method,” in EUROCON, 2007.
The International Conference on ”Computer as a Tool”, Sept 2007, pp.
106–111.
[20] W.-H. Steeb, Matrix Calculus and Kronecker Product with Applications
and C++ Programs. Singapore: World Scientific, 1997.
[21] Z. Ye, C. Liao, X. Xiong, and M. Zhang, “A novel FDTD subgridding
method with improved separated temporal and spatial subgridding
interfaces,” IEEE Antennas Wireless Propag. Lett., vol. 16, pp. 1011–
1015, 2017.
[22] G. Kim, E. Arvas, V. Demir, and A. Z. Elsherbeni, “A novel nonuniform
subgridding scheme for FDTD using an optimal interpolation technique,”
Prog. Electromagn. Res. B, vol. 44, pp. 137–161, 2012.
[23] Y. Wang, S. Langdon, and C. Penney, “Analysis of accuracy and stability
of FDTD subgridding schemes,” in 2010 European Microwave Conf.,
2010, pp. 1297–1300.
[24] A. Ludwig, G. V. Eleftheriades, and C. D. Sarris, “FDTD analysis of
meta-screens for sub-wavelength focusing,” in 2011 IEEE International
Symposium on Antennas and Propagation (APS/URSI), July 2011, pp.
673–676.
[25] L. Markley, A. M. H. Wong, Y. Wang, and G. V. Eleftheriades, “Spatially
shifted beam approach to subwavelength focusing,” Phys. Rev. Lett., vol.
101, p. 113901, 2008.

14

