JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

1

A Dissipative Systems Theory for FDTD with
Application to Stability Analysis and Subgridding

arXiv:1606.08761v1 [] 28 Jun 2016

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

Abstract—This paper establishes a far-reaching connection
between the Finite-Difference Time-Domain method (FDTD) and
the theory of dissipative systems. The FDTD equations for a
rectangular region are written as a dynamical system having
the magnetic and electric fields on the boundary as inputs and
outputs. Suitable expressions for the energy stored in the region
and the energy absorbed from the boundaries are introduced,
and used to show that the FDTD system is dissipative under
a generalized Courant-Friedrichs-Lewy condition. Based on the
concept of dissipation, a powerful theoretical framework to
investigate the stability of FDTD methods is devised. The new
method makes FDTD stability proofs simpler, more intuitive, and
modular. Stability conditions can indeed be given on the individual components (e.g. boundary conditions, meshes, embedded
models) instead of the whole coupled setup. As an example of
application, we derive a new subgridding method with material
traverse, arbitrary grid refinement, and guaranteed stability. The
method is easy to implement and has a straightforward stability
proof. Numerical results confirm its stability, low reflections, and
ability to handle material traverse.
Index Terms—finite-difference time-domain, stability, energy,
dissipation, subgridding.

I. I NTRODUCTION
The Finite-Difference Time-Domain (FDTD) method is
widely used to solve Maxwell’s equations numerically in microwave and antenna engineering, photonics and physics [1],
[2]. FDTD is versatile, easy to implement, and has a low
computational cost per time step. Through explicit update
equations, FDTD recursively computes the electric and magnetic field in the region of interest without requiring the
solution of linear systems. This update process is stable if time
step ∆t satisfies the Courant-Friedrichs-Lewy (CFL) stability
limit [2]
1
∆t < q q
,
(1)
1
1
1
1
εµ
∆x2 + ∆y 2 + ∆z 2
where ∆x, ∆y, ∆z denote cell size, ε denotes permittivity
and µ denotes permeability. While many breakthroughs have
been achieved in FDTD since Yee’s original algorithm [3],
the efficiency of FDTD for multiscale problems remains an
open problem. The simultaneous presence of large and small
geometrical features can dramatically reduce FDTD efficiency
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).

because of two factors. First, the mesh has to be refined,
at least locally, to properly resolve small features, which
increases the number of unknowns and the cost per iteration.
Second, because of (1), a mesh refinement imposes a smaller
time step, which further increases computational cost. For
example, a 3X refinement of the entire FDTD mesh increases
computational cost by 81 times. These issues are unfortunate
since multiscale problems abound in practice.
Numerous solutions have been proposed to mitigate this
issue, including sophisticated boundary conditions to mimic
open spaces [2], local grid refinement [4], [5], [6] (commonly
known as subgridding), thin wire models [1], lumped elements [7], [8], and hybridizations with model order reduction [9], [10], [11], finite elements [12], integral equations [13],
ray tracing [14] and implicit schemes like ADI-FDTD [15],
[16]. From a system theory viewpoint, many of these methods
consist of the interconnection of subsystems, such as models,
algorithms, boundary conditions. For example, in [13], an
FDTD model, used to describe an inhomogeneous scatterer,
is coupled to the time-domain method of moments, used to
model a thin-wire antenna. In this way, one can leverage the
respective strengths of different algorithms. However, while
the stability properties of the individual algorithms may be
well understood, ensuring the stability of their combination
can be a formidable task. Ensuring stability is not trivial even
in the relatively simple case of FDTD subgridding, where
one just couples a coarse and a fine FDTD grid through
an interpolation rule to relate fields at the grid transition.
The lack of a systematic approach to ensure the stability of
advanced FDTD methods is a major issue, that limits the
development of new methods and their adoption by industry,
where guaranteed stability is mandatory. This issue is the
main motivation for this paper, that proposes a new theoretical
framework to investigate and ensure the stability of both
simple and advanced FDTD methods.
Several techniques are available to analyze the stability
of FDTD-like methods. Von Neumann analysis [1] is the
simplest, but is only applicable to uniform meshes and homogeneous materials. The iteration method [2] checks if the
eigenvalues of the matrix that relates the current and next
field solution are all below one in magnitude. In the energy
method [8] instead, one writes an expression for the total
energy stored in the simulation domain, and then checks if
the algorithm satisfies a discrete equivalent of the principle of
energy conservation. Since energy conservation prevents an
unphysical growth of the solution, this implies stability [17].
The iteration and energy methods are general, but they can lead
to long derivations, since they require the analysis of the whole

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

West
boundary

North boundary

(Nx + 1, Ny + 1)

∆x

Ey |1, 3
2

Hz |1, 3

2

(1, 1)

∆y

Hz | 3 , 3

East boundary

coupled scheme, which may consist of many subsystems such
as FDTD meshes, different boundary conditions, reduced order
models, and so on. For example, in a subgridding scenario,
one must derive the iteration matrix or energy function of the
whole scheme, taking simultaneously into account coarse and
fine meshes, interpolation rules and boundary conditions. This
issue makes stability analysis quite involved. Moreover, it does
not yield stability conditions on the individual subsystems.
Consequently, if a subsystem is changed (for example, a different boundary condition is introduced), the iteration matrix or
energy function of the whole problem must be derived again.
In this paper, we propose a new stability framework for
FDTD. The framework is based on the theory of dissipative
systems [18], and generalizes the energy method. First, starting
from FDTD update equations, we develop a self-contained
mathematical model for a region with arbitrary permittivity,
permeability and conductivity. The model is in the form of a
discrete time dynamical system. The magnetic field tangential
to the boundary is taken as input, while the electric field
tangential to the boundary is taken as output. Through this
form, we reveal that an FDTD model can be interpreted as
a dynamical system which is dissipative when ∆t satisfies a
generalized CFL condition. If this condition is not met, the
system can generate energy on its own, leading to an unstable
simulation. We thus establish a connection between FDTD and
the elegant theory of dissipative systems, which is a novel
result. We believe that this connection will greatly benefit the
FDTD community, since the theory of dissipative systems has
been extremely successful in control theory for ensuring the
stability of interconnected systems. This key result sets the
basis for a powerful FDTD stability theory, where each part
of a given FDTD setup (standard meshes, boundary conditions,
interpolation schemes, ...) is interpreted as a subsystem, and is
required to be dissipative. Since the connection of dissipative
systems is dissipative [18], [17], this will ensure the stability
of the FDTD algorithm resulting from the connection of the
subsystems. The proposed theory has numerous advantages. It
simplifies stability proofs, since conditions can be imposed
on each subsystem individually, rather than on the whole
coupled algorithm. Stability proofs are thus made modular
and “reusable”: once a given FDTD model (e.g., an advanced
boundary condition) has been deemed to be dissipative, it can
be combined to any other dissipative FDTD subsystem with
guaranteed stability. The proposed approach also naturally
provides the CFL stability limit of the resulting scheme,
which will be the most restrictive CFL limit of the individual
subsystems. Finally, the theory is intuitive, since it is based on
the concept of energy, familiar to most scientists. The proposed
theory is presented in 2D, for the sake of clarity. An extension
to 3D is feasible and is currently under development.
As an example of application, the proposed theory is
used to derive the subgridding algorithm which is stable
by construction, and has several desirable features. Material
traverse is supported for both dielectrics and highly conductive
materials, which is a limitation of other stable subgridding
methods [19]. The proposed method has low reflections, and
avoids non-rectangular cells [6], finite element concepts [20]
and Withney forms [21], [22]. Corners are natively supported

2

y
z

x

2 2

Hz | 3 ,1 Ex | 3 ,1
2

2

South boundary

Figure 1: Graphical representation of the 2D region considered
in Sec. II. The hanging variables introduced on the four
boundaries are denoted in green.
without any special treatment [5] nor L-shaped cells [6].
Ultimately, the proposed algorithm just consists of a compact
FDTD-like update equation for the edges between coarse and
fine mesh, and is thus easy to implement. This update equation
is provided explicitly for an arbitrary integer refinement ratio
r, while several previous works [5], [6] provide the update
weights only for specific refinements (typically r = 2 or 3),
leaving the derivation of other cases to the Reader.
The paper is organized as follows. In Sec. II, we cast the
FDTD update equations for a rectangular region into the form
of a dynamical system with suitable inputs and outputs. In
Sec. III, we show that FDTD equations can be interpreted as
a dissipative system, and propose the new stability theory. The
theory is applied to derive a stable subgridding algorithm in
Sec. IV, followed by numerical results in Sec. V.
II. D ISCRETE T IME DYNAMICAL M ODEL FOR A 2D FDTD
R EGION
The goal of this section is to cast the FDTD equations for a
2D region into the form of a discrete time dynamical model.
The model shall be self-contained, involving only field samples
from the nodes belonging to the region. This goal will be
achieved by introducing suitable magnetic field samples at the
boundaries.
We consider the 2D rectangular region shown in Fig. 1,
operating in a TE mode with components Ex , Ey , and Hz .
The region is discretized with a uniform rectangular grid with
Nx × Ny cells, of width ∆x and height ∆y. In addition to
the field samples used by standard FDTD, we also sample the
H field on the four boundaries of the region. These additional
samples will be referred to as hanging variables [22], and will
allow us to:
1) develop a self-contained model for the region, which
does not involve field samples beyond its boundaries;
2) derive an expression for the energy absorbed from each
boundary;
3) connect the FDTD grid to other subsystems while maintaining stability.
To keep the notation compact, we collect all Ex and Ey
samples into column vectors Enx and Eny , of size NEx =
Nx (Ny + 1) and NEy = (Nx + 1)Ny , respectively. The Hz

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

3

samples at the internal nodes are collected into column vector
n+ 1
Hz 2 of size NHz = Nx Ny . The hanging variables on the
South boundary of the region are collected into the Nx × 1
vector
h
iT
n+ 1
n+ 1
n+ 1
(2)
HS 2 = Hz |1+ 12,1 . . . Hz |N +2 1 ,1 .
x

2

A. State Equation for Each Node
The dynamical model for the region consists of an update
equation for each E and H sample, excluding the hanging
variables. Those variables will indeed be eliminated once the
region is connected to the surrounding subsystems or to some
boundary conditions.
1) Internal Hz Nodes: For these nodes, we use a standard
FDTD update [2]
µ
µ
n+ 1
n− 1
Hz |i+ 12,j+ 1 = ∆x∆y Hz |i+ 12,j+ 1
2
2
2
2
∆t
∆t
− ∆xEx |ni+ 1 ,j + ∆xEx |ni+ 1 ,j+1
2

2

2

(3)

2

for i = 1, . . . , Nx and j = 1, . . . , Ny . In (3), µ denotes
n+ 1
the average permittivity on the edge where Hz |i+ 12,j+ 1 is
2
2
sampled. Subscripts are omitted from µ in order to simplify
the notation. Equations (3) can be written in matrix form as
Dµ n+ 12
Dµ n− 12
Hz
= DA
Hz
+ Gy Dlx Enx − Gx Dly Eny ,
∆t
∆t
(4)
where DA is a diagonal matrix containing the area of the
primary cells, Dµ is a diagonal matrix containing the average
permittivity on each edge of the secondary grid. Diagonal
matrices
DA

Dly = ∆yINEy

(5)

contain the length of the x- and y-directed edges of the primary
grid, respectively. With Im , we denote the m × m identity
matrix. Matrix Gx is the discrete derivative operator along x,
which can be written as
Gx = INy ⊗ WNx ,

(6)

where ⊗ is the Kronecker’s product [23] and WNx is the
Nx × (Nx + 1) matrix [24]


−1 +1


−1 +1


Wn = 
(7)
.
..
..


.
.
−1 +1
Similarly,
Gy = WNy ⊗ INx .

2

2

2

for i = 1, . . . , Nx and j = 2, . . . , Ny . In (9), εx and σx denote,
respectively, the average permittivity and average conductivity
on the corresponding edge. While the common factor ∆x
could be eliminated from (9), it is kept as it will be useful
later.
For the Ex nodes on the South boundary, the standard
FDTD update equation involves a Hz sample outside the
region. In order to avoid this, we apply the finite difference
approximation to the half step of the
 secondary grid between
nodes i + 12 , 1 + 12 and i + 12 , 1 , making use of the hanging variables. In this way, we obtain an update equation that
involves only field samples from the considered region
σx 
∆y  εx
+
Ex |n+1
=
∆x
i+ 21 ,1
2 ∆t
2

σx 
∆y εx
−
Ex |ni+ 1 ,1
∆x
2
2 ∆t
2
n+ 12
n+ 1
+ ∆xHz |i+ 1 ,1+ 1 − ∆xHz |i+ 12,1 , (10)
2

+ ∆yEy |ni,j+ 1 − ∆yEy |ni+1,j+ 1

Dlx = ∆xINEx

2

2

Similarly, the hanging variables on the North, East and West
n+ 1
n+ 1
n+ 1
sides are cast into column vectors HN 2 , HE 2 and HW 2 ,
respectively.

∆x∆y

2) Ex Nodes: For the Ex nodes that fall strictly inside the
region, we use a standard FDTD update [2]
ε
ε
σx 
σx 
x
x
∆x∆y
Ex |n+1
=
∆x∆y
Ex |ni+ 1 ,j
+
−
1
i+ 2 ,j
2
∆t
2
∆t
2
n+ 21
n+ 21
+ ∆xHz |i+ 1 ,j+ 1 − ∆xHz |i+ 1 ,j− 1 (9)

(8)

2

2

for i = 1, . . . , Nx . In (10), εx and σx denote the material

properties of the half cell between nodes i + 12 , 1 + 12 and
i + 12 , 1 . The update equation for the Ex nodes on the North
boundary is obtained similarly, and reads
∆y  εx
σx 
∆x
=
+
Ex |n+1
i+ 21 ,Ny +1
2 ∆t
2


∆y εx
σx
∆x
−
Ex |ni+ 1 ,Ny +1
2
2 ∆t
2
n+ 21
n+ 1
+ ∆xHz |i+ 1 ,N +1 − ∆xHz |i+ 12,N + 1 , (11)
2

y

2

y

2

for i = 1, . . . , Nx . Relations (9), (10) and (11) can be
compactly written as


D εx
Dσx
Dlx Dly0
+
En+1
=
x
∆t
2


Dσx
Dεx
n+ 1
−
Enx − Dlx GTy Hz 2
Dlx Dly0
∆t
2
"
1#
2

 Hn+
S
+ Dlx BS Dlx BN
, (12)
n+ 1
HN 2
where
• Dl0 is an NEx × NEx diagonal matrix containing the
y
length of the y-directed edges of the secondary grid,
including the half-edges of length ∆y/2 that intersect
the North and South boundaries;
• Dεx and Dσx are diagonal matrices storing the permittivity and conductivity on each x-directed primary edge,
respectively;
• BS has all entries set to zero, except for a −1 at
the intersection of each row associated with a South
boundary edge and the column of the corresponding
n+ 1
hanging variable in HS 2 ;

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

4

BN is defined similarly to BS , but has a +1 on each
entry associated with an edge of the North boundary.
3) Ey Nodes: The update equations for the Ey nodes are
derived with a similar procedure. A standard FDTD update
equation is written for the Ey nodes that fall strictly inside the
region. For the Ey nodes on the West and East boundaries,
instead, finite differences are applied on a half grid step,
similarly to (10) and (11). This process leads to


Dσy
D εy
+
En+1
=
Dly Dlx0
y
∆t
2


Dσy
D εy
n+ 1
Dly Dlx0
−
Eny + Dly GTx Hz 2
∆t
2
#
"
n+ 21

 HW
, (13)
+ Dly BW Dly BE
n+ 1
HE 2
•

where
• Dl0 is an NEy × NEy diagonal matrix containing the
x
length of the x-directed edges of the secondary grid,
including the half-edges of length ∆x/2 that intersect
the East and West boundaries;
• Dεy and Dσy are diagonal permittivity and conductivity
matrices associated with the y-directed edges of the
primary grid, respectively;
• BW has all entries set to zero, except for a +1 at
the intersection of each row associated with a West
boundary edge and the column of the corresponding
n+ 1
hanging variable in HW 2 ;
• BE is analogous to BW , but has a −1 on each entry
associated with an East boundary edge.
B. Descriptor System Formulation
The update equations derived in the previous sections form
a complete dynamical model for the rectangular region
(R + F)xn+1
y

n

=

1

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

n

= L x .

(14a)

equation (14b) extracts these values from the state vector xn .
The coefficients matrices in (14a) and (14b) read


Dεx
1
T
Dlx Dly0 ∆t
0
2 Dlx Gy


Dε
R=
0
Dly Dlx0 ∆ty − 12 Dly GTx  , (17)
D
1
− 12 Gx Dly
DA ∆tµ
2 Gy Dlx


D
1
T
Dlx Dly0 2σx
0
D
G
l
y
x
2


Dσ
F=
0
Dly Dlx0 2 y − 21 Dly GTx  , (18)
1
− 12 Gy Dlx
0
2 Gx Dly


Dlx BS Dlx BN
0
0
0
Dly BW Dly BE  ,
B= 0
(19)
0
0
0
0


−BS BN
0
0
0
BW −BE  .
L= 0
(20)
0
0
0
0
Equations (14a)-(14b) define an FDTD-like model for a
2D lossy, source-free region, with possibly non-uniform permittivity and permeability. From a control perspective, this
model is a discrete time descriptor system [25], also known
as generalized state-space. From an electrical standpoint, the
model is an impedance-type description of the region, because
it gives the electric field tangential to the boundary in terms
of the tangential magnetic field. In [24], a special case of
the proposed representation was derived for a lossy region
enclosed by PEC boundaries. The proposed developments
generalize [24] in two directions:
1) the assumption of a PEC termination is removed. Suitable input-output variables are introduced in such a
way that model (14a)-(14b) can be connected to other
subsystems, such as different FDTD grids, boundary
conditions, or reduced models;
2) the proposed formulation will allow us to prove, in the
next section, that FDTD equations can be interpreted as
a dissipative dynamical system.

(14b)
III. FDTD AS A D ISSIPATIVE D ISCRETE T IME S YSTEM

The first equation (14a) is formed by update equations (4),
(12) and (13), and updates the state vector
 n 
Ex
n
xn =  Ey  ,
(15)
n− 12
Hz
which consists of all electric and magnetic field samples in the
1
region, excluding hanging variables. The input vector un+ 2
n
and output vector y are given by
 n+ 1 
 n
H 2
ES
 Sn+ 21 
n 



1
H
E
N
 yn =  nN  .
un+ 2 = 
(16)
1
 n+ 2 
EW 
HW 
n
EE
n+ 1
HE 2
The input vector contains all hanging variables, i.e. all magnetic field samples on the region boundaries. The output
vector is made by the E samples at the same nodes, which
are collected into vectors EnS , EnN , EnW and EnE . Output

In this section, we define suitable expressions for the energy
stored in the 2D region and the energy absorbed from its
four boundaries. These expressions are then used to derive
the conditions under which system (14a)-(14b) is dissipative.
A. Dissipation Inequality
A discrete time system like (14a)-(14b) is dissipative when
it satisfies the following condition [26].
Definition 1. Dynamic system (14a)-(14b) is said to be
1
dissipative with supply rate s(yn , un+ 2 ) if there exists a
n
nonnegative function E(x ) with E(0) = 0, called storage
function, such that
1

E(xn+1 ) − E(xn ) ≤ s(yn , un+ 2 )
n+ 12

for all u

(21)

and all n.

The storage function E(xn ) can be interpreted as the energy stored in the system at time n, while the supply rate

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

5

1

s(yn , un+ 2 ) is the energy absorbed by the system from its
boundaries between time n and n+1. Clearly, only dissipative
systems can satisfy the dissipation inequality (21). Indeed,
their stored energy can increase between time n and n + 1
by at most the energy absorbed from the outside world. If
this limit is violated, the system is considered active, since it
can generate energy on its own. For most practical systems,
inequality (21) is satisfied strictly because of the presence of
losses.
B. Storage Function and Supply Rate
For system (14a)-(14b), we choose as candidate storage
function
∆t n T
(x ) Rxn ,
(22)
E(xn ) =
2
and as supply rate
n+ 12

n

s(y , u

yn + yn+1
) = ∆t
2

T

1

LT Bun+ 2 .

Substituting (27) into (23), we obtain
(EnS + En+1
)T n+ 12
S
HS
2
n+1 T
n
n
T
(E + EN )
(E + En+1
n+ 1
n+ 1
W )
+∆t∆x N
HN 2 +∆t∆y W
HW 2
2
2
T
(EnE + En+1
n+ 1
E )
HE 2 , (28)
− ∆t∆y
2
and see that the supply rate is the sum of the energy absorbed
by the region from each boundary between time n and n + 1.
Signs in (28) are consistent with the direction of the Poynting
vector on each boundary.
1

s(yn , un+ 2 ) = −∆t∆x

C. Dissipativity Conditions
Using the proposed storage function (22) and supply
rate (23), we can derive simple dissipativity conditions on the
coefficients matrices R, F, B and L in (14a)-(14b).
Theorem 1. If

(23)

R = RT > 0 ,

(29a)

Before showing that these functions satisfy (21), we investigate
their physical meaning. Substituting (17) into (22), we get

F + FT ≥ 0 ,

(29b)

B = LLT B ,

(29c)

1 n T
(E ) Dlx Dly0 Dεx Enx
2 x
T
1
Eny Dly Dlx0 Dεy Eny
+
2


∆t
Dµ n− 21
n
n
Hz
+
Gy Dlx Ex − Gx Dly Ey + DA
. (24)
2
∆t

E(xn ) =

Since the term between square brackets is the right hand side
of (4), we can reduce (24) to
1
T
E(x ) = (Enx ) Dlx Dly0 Dεx Enx
2
T
1
Eny Dly Dlx0 Dεy Eny
+
2
1  n− 12 T
n+ 1
Hz
DA Dµ Hz 2 . (25)
+
2
n

Z 
A


1 2
1
1 2
2
εE (t) + εEy (t) + µHz (t) dA ,
2 x
2
2

(26)

the continuous expression for the electromagnetic energy per
unit height of a 2D TE mode in a region of area A. The first
term of (25) is the energy stored in the x component of the
electric field in the region. Indeed, the diagonal matrix Dlx Dly0
is the area of the half cells around the Ex edges. Similarly,
the second and third term in (25) represent the energy stored
in the Ey and Hz components, respectively.
We now discuss the physical meaning of supply rate (23).
Direct inspection reveals that


−∆xINx
0
0
0


0
+∆xINx
0
0
.
LT B = 


0
0
+∆yINy
0
0
0
0
−∆yINy
(27)

then system (14a)-(14b) is dissipative according to Definition 1, with (22) as storage function and (23) as supply rate.
Proof. Condition (29a) makes storage function (22) nonnegative for all xn , as required by Definition 1. Next, we show
that if (29b) and (29c) hold, the dissipation inequality (21) will
hold as well. Substituting (22) and (23) into (21), we obtain
T
T
1
T
xn+1 Rxn+1 −(xn ) Rxn − yn + yn+1 LT Bun+ 2 ≤ 0 .
This inequality can be rewritten as
T

T
xn+1 R xn+1 − xn + xn+1 − xn Rxn
T
1
− yn + yn+1 LT Bun+ 2 ≤ 0

(30)

From (14a), we have that


1
R xn+1 − xn = −F xn+1 + xn − Bun+ 2
(31)


T
Using (31), the terms R xn+1 − xn and xn+1 − xn R
in (30) can be expressed in terms of F
T

T
− xn+1 F xn+1 + xn − xn+1 + xn FT xn

T
T
1
1
+ xn+1 Bun+ 2 + un+ 2
BT xn
T
1
1
T
− (xn ) LLT Bun+ 2 − xn+1 LLT Bun+ 2 ≤ 0 (32)
Under (29c), we can finally rewrite (32) as
T


xn+1 + xn
F + FT xn+1 + xn ≥ 0

(33)

which is clearly satisfied if (29b) holds.
We now investigate the physical meaning of the three
dissipativity conditions, starting from (29c), which can be
directly verified by combining (27) and (20). This condition
holds because the region inputs and outputs in (16) are
sampled at the same nodes, and is reminiscent of a similar

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

6

relation which holds for linear circuits under the impedance
representation [27]. Condition (29b) is related to losses, and
reads


Dlx Dly0 Dσx
0
0
F + FT = 
0
Dly Dlx0 Dσy 0 ≥ 0 . (34)
0
0
0
Since Dlx , Dly , Dlx0 and Dly0 are diagonal with positive
elements only, (34) boils down to
D σx

≥ 0,

(35)

Dσy

≥ 0.

(36)

As one may expect, the average conductivity on each primary
edge must be non-negative in a dissipative system. Condition (29a) is the most interesting, and is a generalized CFL
criterion. From (17), we see that R = RT by construction.
With the Schur’s complement [28], we can transform (29a)
into
Dlx Dly0 Dεx > 0 ,

(37)

Dly Dlx0 Dεy > 0 ,
4
IN ,
SST <
∆t2 Hz

(38)
(39)

where
h
1
−1
−1
−1
−1
S = DA 2 Dµ 2 Gy Dl2x Dl0 2 Dεx2

i
1
−1
−1
−Gx Dl2y Dl0 2 Dεy2 .
y
x
(40)
Inequalities (37) and (38) are clearly satisfied by construction.
If we denote with sk the singular values of (40), inequality (39)
will hold if
2
∆t <
∀k .
(41)
sk
This condition sets an upper bound on the FDTD time step,
and is analogous to the generalized CFL constraint derived
in [24] for an FDTD grid terminated on PEC boundaries.
From (40), we indeed see that the time step upper limit
depends on grid size, permittivity and permeability. Using a
time step which violates (41) has two consequences. First, the
discretization of Maxwell’s equations leads to an active discrete time model, even if the real physical system has positive
losses everywhere, and should therefore be dissipative. The
numerical model is thus inconsistent with the actual physics.
Second, being able to generate energy on its own, model (14a)(14b) can lead to divergent transient simulations. Even when
such model is connected to other dissipative subsystems, its
ability to generate energy on its own can destabilize the
whole simulation. The proposed theory provides a new, deeper
understanding on the root causes of FDTD instability, and is
able to pinpoint which part of an FDTD setup is responsible
for it. It is hoped that this new explanation will facilitate
the development of powerful FDTD schemes with guaranteed
stability.
D. Relation to CFL Stability Limit
The relation between (41) and the CFL limit can be further
understood if we apply (41) to a single cell. The purpose
of (41) is to make storage function (22) non-negative. A

sufficient condition for this is to require the energy stored in
each primary cell to be positive. This condition can be derived
by applying (39) to a single primary cell. We consider the
primary cell extending from node (i, j) to node (i + 1, j + 1).
The coefficient matrices for this special case can be obtained
from the formulas in Sec. II by setting Nx = Ny = 1
Dlx = ∆xI2
∆x
I2
Dlx0 =
2

Gx = Gy = −1

1

Dly = ∆yI2
∆y
Dly0 =
I2
2
DA = ∆x∆y



(42)
(43)
(44)

The permeability matrix Dµ = µ|i+ 12 ,j+ 12 is the average
permeability on the secondary edge, while the permittivity
matrices




0
0
εx |i+ 12 ,j
εy |i,j+ 12
Dεx =
Dεy =
0
εx |i+ 21 ,j+1
0
εy |i+1,j+ 21
(45)
contain the average permittivity on the four primary edges.
Substituting (42)-(45) into (39), we obtain
"
!
1
1
1
∆t <
+
2∆x2 µ|i+ 12 ,j+ 12 εy |i,j+ 12
εy |i+1,j+ 12
! #− 21
1
1
1
, (46)
+
+
2∆y 2 µ|i+ 21 ,j+ 12 εx |i+ 12 ,j
εx |i+ 21 ,j+1
which is a generalized CFL condition. For a uniform medium
with permittivity ε and permeability µ, (46) reduces to
∆t < q

1
1
εµ

q

1
∆x2

+

1
∆y 2

,

(47)

the CFL limit of 2D FDTD. This derivation confirms that (29a)
is a generalized CFL condition, here reinterpreted in the
context of dissipation.
E. Application to Stability Analysis
The proposed dissipation theory can be effectively used
to investigate and enforce the stability of FDTD algorithms.
Most FDTD setups consist of an interconnection of FDTD
subsystems. In the simplest scenario, a uniform FDTD grid
is connected to some boundary conditions. In most advanced
scenarios, one may want to couple a main FDTD grid to
refined grids, reduced models, lumped elements, or models
from other numerical techniques, such as finite elements,
integral equations or ray tracing.
Ensuring stability of these hybrid schemes can be very
challenging, since stability is a property of the overall scheme,
rather than of its individual subsystems. By invoking the concept of dissipation, we can instead achieve stability in an easy
and modular way. Each part is seen as an FDTD subsystem
and required to satisfy dissipativity conditions (29a)-(29c).
Since the connection of dissipative systems is dissipative by
construction [17], the overall method will be guaranteed to be
stable. The proposed theoretical framework generalizes the socalled energy method [8], and has numerous advantages over
the state of the art:

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

1) non-uniform problems can be handled, unlike in the von
Neumann analysis [1];
2) stability conditions can be given on each subsystem
separately, unlike in the iteration method [2], which
requires the analysis of the iteration matrix of the whole
scheme. This makes stability analysis modular and thus
simpler;
3) once some given FDTD models have been proven dissipative, they can be arbitrarily interconnected without
having to carry out further stability proofs. With the
iteration method, when a single part of a coupled scheme
changes, the whole proof must be revised;
4) the CFL limit of the resulting scheme can be easily
determined by applying (37), (38) or (39) to each
subsystem and taking the most restrictive CFL limit;
5) the stability framework is intuitive, since it is based on
the fundamental physical concept of energy dissipation.

7

∆x
r

Ĥz |ı̂+ 12 ,̂+ 12

Ĥz |ı̂+ 32 ,̂+ 12

∆y
r

(ı̂, ̂) Ĥz |ı̂+ 21 ,̂ Êx |ı̂+ 12 ,̂ Ĥz |ı̂+ 32 ,̂ Êx |ı̂+ 32 ,̂
Hz |i+ 1 ,j

(i, j)

Ex |i+ 1 ,j

2

y

2

z

Hz |i+ 1 ,j− 1
2

x

2

Figure 2: Subgridding scenario considered in Sec. IV for the
case of r = 2. For clarity, a virtual gap has been inserted
between the two grids. This virtual gap is closed when the
two grids are connected.
Fine mesh
Ĥz |ı̂+ 1 ,̂
2

Êx |ı̂+ 3 ,̂ Ĥz |ı̂+ 3 ,̂
2

2

Êx |ı̂+ 1 ,̂
2

IV. A PPLICATION : S TABLE FDTD S UBGRIDDING
We demonstrate the proposed theory by deriving a subgridding algorithm which is stable by construction, easy to
implement and supports an arbitrary grid refinement ratio. The
goal is to derive stable update equations for a setup where one
or more fine grids are embedded in a main coarse grid. Without
loss of generality, we consider the case where a coarse grid
with cell size ∆x × ∆y hosts a single fine grid with cell size
∆x/r × ∆y/r, where r is an arbitrary integer. The algorithm
will ultimately consist of conventional FDTD equations to
update the fields that fall strictly inside the two grids, and
a special update equation to update the fields on the edges
at the grid transition. To derive the method, it is sufficient
to consider the interface between a single coarse cell and the
corresponding fine cells, as shown in Fig. 2. In the figure a
virtual gap has been opened between the two grids for clarity.
Without loss of generality, we consider a refinement in the
positive x direction. The other three cases can be derived in
the same way. The coarse cell under consideration is centered
at node (i+ 12 , j − 12 ) and the corresponding fine cells at nodes
(ı̂ + 12 , ̂ + 21 ), ... , (ı̂ + r − 12 , ̂ + 21 ), where coordinates (ı̂, ̂)
correspond to the same physical location as (i, j) in the coarse
grid. Superscript “ˆ” denotes variables related to the fine grid.
A subgridding algorithm can be interpreted as a the result
of the connection of the three subsystems, as shown in
Fig. 3. Two of those subsystems correspond to the coarse and
fine grids to be coupled. The third subsystem represents the
interpolation rule which is used to relate the fields on the
boundaries of the two grids, that are sampled with different
resolution.
A. State Equations for the Coarse and Fine Cells on the
Boundary
For the coarse and fine grids, we adopt the formulation of
Sec. II, introducing hanging variables on the two boundaries
to be connected. The purpose of the hanging variables is
to facilitate the coupling of the two meshes and the proof
of stability. These extra variables will be eliminated when
deriving the update equation for the fields at the interface.

Interpolation rule
Hz |i+ 1 ,j
2

Ex |i+ 1 ,j
2

Coarse mesh

Figure 3: Interpretation of the subgridding method as the
connection of three dynamical systems, representing the coarse
grid, the fine grid, and the interpolation rule.
On the North boundary of the coarse cell, we introduce the
hanging variable
n+ 12

HN

n+ 1

= Hz |i+ 12,j .

(48)

2

Similarly, on the South fine cell boundaries we introduce the
hanging variables
h
iT
n+ 1
n+ 1
n+ 12
ĤS 2 = Ĥz |ı̂+ 12,̂ . . . Ĥz |ı̂+r−
,
(49)
1
,̂
2

2

as shown in Fig. 2.
From (11), we obtain the following state equation for the
coarse E-field sample at the interface and the coarse hanging
variable
σx  n+1
∆y  εx
+
EN =
2 ∆t
2

∆y εx
σx  n
n+ 1
n+ 1
−
EN + HN 2 − Hj− 12 , (50)
2
2 ∆t
2
where
n+ 1

n+ 1

Hj− 12 = Hz |i+ 12,j− 1 ,
2

2

2

n
EN
= Ex |ni+ 1 ,j ,

(51)

2

and εx and σx are the permittivity and conductivity on the
primary edge of the coarse cell below the interface. Similarly,
the state equation for the r fine cells can be written as
!
∆y D̂εx
D̂σx
+
Ên+1
=
S
2r
∆t
2
!
∆y D̂εx
D̂σx
n+ 1
n+ 1
−
ÊnS + Ĥ̂+ 12 − ĤS 2 , (52)
2
2r
∆t
2

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

8

where
n+ 1
Ĥ̂+ 12
2

h
iT
n+ 1
n+ 21
= Ĥz |ı̂+ 12,̂+ 1 . . . Ĥz |ı̂+r−
,
1
1
2
2
2 ,̂+ 2
iT
h
ÊnS = Êx |nı̂+ 1 ,̂ . . . Êx |nı̂+r− 1 ,̂ ,
2

2

(53)
(54)

and D̂εx and D̂σx are r × r diagonal matrices containing the
values of permittivity and conductivity above the interface for
the South edges where the fine Êx fields are sampled.
B. Interpolation Rule
The interpolation rule [29] relates the field sampled on
the coarse and fine grids. From the boundary condition for
tangential electric fields, we have that
n
ÊnS = EN
T

∀n ,

(55)

where T is an r × 1 matrix of ones. Condition (55) forces
the coarsely- and finely-sampled E fields to be equal at all
times.
On the magnetic fields at the boundary, we impose a
constraint reciprocal to (55)
n+ 12

n+ 1
TT ĤS 2

∀n .
(56)
r
We will see in Sec. IV-D that the reciprocity between the E
and H interpolation rules is required to ensure stability.
HN

=

C. Explicit Update Equation for the Interface
Interpolation conditions (55) and (56) can now be used
to combine (50) and (52) in order to derive an explicit
update equation for the fields at the coarse-fine interface, and
eliminate hanging variables.
Substituting (55) into (52), and multiplying the obtained
equation by TT /r on the left yields
!
∆y TT D̂εx T TT D̂σx T
n+1
+
EN
=
2r
r∆t
2r
!
∆y TT D̂εx T TT D̂σx T
n
−
EN
2r
r∆t
2r
n+ 1

TT Ĥ̂+ 12

n+ 1

TT ĤS 2
+
−
. (57)
r
r
For the simplicity of notation we define symbols for the
average permittivity and conductivity of the r South boundary
fine cells
2

TT D̂σx T
TT D̂εx T
, σ̂x =
.
r
r
Equation (57) can now be added to (50), yielding
!
σx + σ̂rx
∆y εx + ε̂rx
n+1
+
EN
=
2
∆t
2
!
σx + σ̂rx
∆y εx + ε̂rx
n
−
EN
2
∆t
2
ε̂x =

n+ 12

+ HN

n+ 1

− Hj− 12 +
2

n+ 1
TT Ĥ̂+ 12
2

r

(58)

TT ĤS
r

Rearranging (60), we get the following explicit update equation for EN in terms of the neighboring magnetic fields
!−1
!
σ̂x
ε̂x
σ̂x
ε̂x
σ
+
ε
+
σ
+
ε
+
x
x
x
x
n+1
n
r
r
r
r
EN
+
−
EN
=
∆t
2
∆t
2

!−1  T n+ 21
T Ĥ̂+ 1
1
εx + ε̂rx
σx + σ̂rx
2
n+
2

+
+
− Hj− 12  .
2
∆y
∆t
2
r
(61)
The fine interface electric fields are then updated using (55).
It should be noted that, when r = 1, equation (61) reduces to
the standard FDTD update equation.
The overall subgridding algorithm can be summarized as
follows:
1) Calculate the magnetic field samples everywhere at time
n + 12 using conventional FDTD update equations.
2) Use standard FDTD update equations to compute the E
fields at time n + 1 on the edges that are strictly inside
the coarse and fine grids.
n+1
3) Compute EN
, the coarsely-sampled electric field at
the interface, using (61).
4) Update the finely-sampled Êx fields at the interface
using (55).
The computational overhead of this scheme is minimal, since
the coefficients in (61) can be pre-computed before the update
iterations.
The proposed method is thus simple to implement, since it
consists of conventional FDTD update equations inside the two
meshes and a modified update equation for the edges at the
interface. In comparison to previous subgridding methods, we
avoid non-rectangular cells [6], finite element concepts [20]
and Withney forms [21], [22]. The proposed update equation
can be also used at corners with no modifications, unlike in
previous works that require special treatment [5] or L-shaped
cells [6]. Finally, we remark that in Sec. II, the FDTD update
equations have been given in matrix form in order to reveal the
dissipative nature of FDTD systems. This form, however, does
not have to be used in the practical implementation, which can
use conventional for loops or, in languages like MATLAB,
vectorized operations.
D. Proof of Stability

n+ 21

−

With the interpolation rule (56), we cancel the hanging variables to obtain
!
σx + σ̂rx
∆y εx + ε̂rx
n+1
+
EN
=
2
∆t
2
!
n+ 12
T
T
Ĥ
1
σx + σ̂rx
∆y εx + ε̂rx
̂+ 12
n+
n
EN
− Hj− 12 +
−
.
2
2
∆t
2
r
(60)

. (59)

The proposed dissipation theory makes it straightforward to
prove that the subgridding algorithm is stable under the CFL
limit of the fine grid. For stability, all three subsystems in
Fig. 3 need to be dissipative, which requires one to use the

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

9

more restrictive fine grid time step. With the time step chosen
correctly, in order to guarantee the overall stability we need
to only ensure dissipativity of the interpolation rule.
Analogously to (23) the supply rate for the interpolation
subsystem is defined as
n+1
EN

+
n+ 1
) = −∆t∆x
HN 2
s(yn , u
2
)T n+ 21
∆x (ÊnS + Ên+1
S
+ ∆t
ĤS . (62)
r
2
Substituting (55) and (56) into (62), we have
n+ 12

20 mm

∆x = 1 mm
∆y = 2 mm

40 mm

n
EN

PEC
Probe
y

∆x̂ = ∆x/4
∆ŷ = ∆y/4

z

x

40 mm

Source

60 mm
Figure 4: Layout of the PEC cavity considered in Sec. V-A.

1

s(yn , un+ 2 ) =
n+ 12

TT ĤS
r

−

n+ 1
HN 2

!
= 0 . (63)

Therefore, the proposed interpolation rule is a lossless system
that does not dissipate nor absorb any energy. Physically, this
result makes sense, since the connection system corresponds to
an infinitely thin region where no energy dissipation can take
place. In conclusion, since the proposed subgridding method
can be seen as the connection of three dissipative systems, it
is overall dissipative, and thus stable.

0.5
Hz (A/m)

n+1 T
(E n + EN
)
∆t∆x N
2

0
−0.5
0

2

4

6
Timestep

8

10
5

x 10

Figure 5: Magnetic field at the probe for the the empty cavity
with subgridding of Sec. V-A, computed for 106 time steps.

V. N UMERICAL E XAMPLES
The following sections provide the results of FDTD simulations that were done to verify the proposed theory. The
subgridding algorithm was implemented in Matlab and tests
were performed in order to check its stability, ability to handle
material traverse, its accuracy and speedup capability.
A. Stability Verification
Stability was verified by simulating an empty cavity with
perfect electric conductor (PEC) walls with a centrally placed
subgridding region for 106 time steps. The layout of the
simulation is shown in Fig. 4. The cavity was excited using
a modulated Gaussian magnetic current source with central
frequency of 3.75 GHz and half-width at half-maximum of
0.74 GHz. Magnetic field was recorded at a probe placed
inside the cavity. The time step was set 1% below the CFL
limit of the fine grid.
The resulting waveform in Fig. 5 shows that no instability
occurred after 106 time steps. Stable behavior after such a
large number of time steps verifies the correctness of the
proposed stability enforcement technique, especially since no
lossy materials were present to dissipate any spurious energy
artificially created by the algorithm.
B. Material Traverse
The ability of the proposed method to produce meaningful
results when objects traverse the subgridding interface was
tested using the setup in Fig. 6. A 16 × 16 mm slab of
material was simulated for three different placements of the
subgridding region: enclosing, traversing and away from the
slab. The test was done for copper and a lossy dielectric
with conductivity of 5 S/m and relative permittivity of 2.

As a reference, uniformly discretized all-coarse and all-fine
simulations were performed at the fine time step, in addition
to the subgridding simulations. Coarse and fine meshes in
uniformly discretized and subgridding runs were chosen as
1 mm and 0.2 mm respectively. 15 mm-thick perfectly matched
layer (PML) terminated the simulation region. Modulated
Gaussian magnetic current excitation was used at 15.0 GHz
central frequency with 8.82 GHz half-width at half-maximum
bandwidth. The time step was chosen as 0.467 ps in all test
cases.
The magnetic field waveforms at the probe recorded in the
different subgridding scenarios are shown in Fig. 7, and are in
excellent agreement among each other. This result confirms
that the proposed subgridding method can properly handle
material traverse, for both very good conductors and for lossy
dielectrics.
C. Application to Simulating Scatterer Reflections
In order to investigate the accuracy of the proposed scheme,
we looked at waveguide reflections from the scatterer shown in
Fig. 8, which consisted of four copper rods with 1 mm radius.
We have also investigated the reflections from the subgridding
interface, in order to assess the quality of the subgridding
scheme. The coarse cell was chosen to be 1 mm, which was
1
10 of the minimum wavelength of interest that corresponded
to 30.0 GHz. The fields of the incident wave were computed
by running a reference simulation without the scatterer and
without subgridding. Reflected wave fields were found by
subtracting the incident wave fields from the total fields in
simulations with the scatterer. 15 mm-thick PML boundary
was chosen to terminate the two sides of the waveguide. The

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

10

PEC

40 mm

x

z
Source
subgrids

y

Jy current
Probes
8×8 mm

y

Probe

PML

PML

x

z
PML

20 mm

2 mm

∆x = 1 mm
∆y = 1 mm

Figure 6: Layout used for the material traverse test of Sec. V-B,
and the three different placements of the subgridding region.

2 mm
Scatterer

66 mm

Figure 8: Layout of the four-rod reflection simulation discussed in Sec. V-C. The dashed line shows the location of
the subgrid in the subgridding run.
Lossy dielectric

−5
With scatterer

−10
0

11

|S | (dB)

Hz (A/m)

0.02

−0.02

−15
−20
−25
−30
−35

Copper

Sub−grid interface only

−40
0

|S | (dB)

Hz (A/m)

0.02

11

−0.02
0.2

0.25
Time (ns)

0.3

0.35

Figure 7: Time-domain magnetic field at the probe recorded for
the three locations of the subgrid in Sec. V-B: enclosed (
),
outside (
) and traversing the object (
). Waveforms
from the uniformly discretized all-coarse (  ) and allfine (
) simulations are also shown.

time steps in the subgridding runs were set 1% below the CFL
limit of the refined region. Uniformly discretized simulations
were run 1% below the CFL limits of the respective grids.
The resulting reflections are shown in Fig. 9 for the subgridding case with different refinement ratios, as well as the
reference run with full refinement by a factor of 6. The
simulation times are shown in Table I. It can be seen that
the local refinement of the grid around the scatterer with
the proposed subgridding method can improve the accuracy
substantially compared to the coarse grid run. The larger
choice of refinement ratio makes the solution very close to
the reference all-fine solution. Moreover, very good speedup
- almost by a factor of 11 - is achieved even for the grid
refinement of 6 when the redundant high resolution of the grid
in air is eliminated. The reflections from subgridding interface
were significantly lower than those from the scatterer, further
demonstrating the accuracy of the proposed method.

−60
−80
−100
0

5

10
15
20
Frequency (GHz)

25

Figure 9: Reflected power with respect to the incident for the
example of Sec. V-C. Top panel: reflections from the four-rod
scatterer for different global discretization: all-coarse (  )
and all-fine with r = 6 (
); and for the subgridding runs:
r=2(
), r = 4 (
) and r = 6 (
). Bottom panel:
reflections from the subgrid interface only.

D. Application to Exposure Studies
We show the possibility of applying the method for multiscale human exposure simulations. The chosen setup is shown
in Fig. 10. A transverse cross-section of the head of the model
of Ella from IT’IS Virtual Population V1.x [30] was used
to assign the material properties [31] to the FDTD cells. A
900 MHz source was placed approximately 3 meters away
from the human head. The simulation region was terminated
with 20 cm PML. The reference (fine) resolution in the tissues
was set to 2 mm, based on the mesh size chosen by [32] in a
radiation exposure study at that frequency. Specific absorption
rate, or SAR, was evaluated according to the formula in [33],
which was used as follows for sinusoidal excitation
SAR =

σ(Ex2p + Ey2p )
2ρ

,

(64)

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

Table I: Simulation times for different meshes in four-rod
scatterer simulations in Sec. V-C. Simulation times were
recorded for the total field run only - not for the reference
run.
Method
All-fine (r = 6)
All-coarse
Subgridding (r = 2)
Subgridding (r = 4)
Subgridding (r = 6)

Simulation time (s)
427.4
3.2
10.7
23.3
39.5

3.05m

Head
slice

Method
All-fine
All-coarse
Subgridding

All-fine

z

Error in SAR integral
Not applicable
59.5%
-3.1%

Subgridding

Simulation time (s)
1596.4
6.9
34.5

All-coarse

x

PML
4.01m

Figure 10: Layout of the simulation in Sec. V-D with human
head cross-section placed approximately 3 meters away from
a point source.
where subscript “p” denotes the peak absolute value of a
field component and ρ corresponds to tissue density. SAR
was calculated for each of the primary cells. The values of
the electric field components at the primary cell centers were
found by averaging the nearest known samples at the cell
edges. The peak values of the fields were found for the time
interval from 25.6 ns to 26.8 ns, which gave the wave sufficient
time to reach the head and penetrate inside it.
In the subgridding run, the empty space was coarsened to
1
1 cm, which corresponded to 33.3
of the wavelength. The
reference all-fine run was performed at 2 mm resolution, along
with the all-coarse run where the entire simulation region was
discretized at 1 cm. The volumetric integral of SAR over the
tissues was used as an accuracy metric
XX
SAR|i+ 12 ,j+ 12 ∆x∆y ,
(65)
i

Table II: Error in the integral of SAR and simulation times in
the simulations discussed in Sec. V-D.

y

subgrid
Source

11

j

where ∆x and ∆y are FDTD cell dimensions in the tissues.
Time step of 4.67 ps was chosen for the all-fine simulation
and for the subgridding simulation. The coarse grid case was
run at 23.11 ps.
The resulting SAR maps are shown in Fig. 11. Table II
shows simulation times and percent error in the total SAR
with respect to the all-fine simulation. The results show that
with subgridding a speedup of 46X can be obtained with only
3.1% loss in accuracy. No noticeable difference can be seen
on the SAR maps in Fig. 11 between the subgridding and
all-fine simulations. When, instead, the coarse resolution was
chosen for the entire grid, the head tissues were not sufficiently
resolved and the integral SAR differed from the reference by
59.5%, showing the necessity of local grid refinement.
VI. C ONCLUSION
This paper proposed a dissipative systems theory for FDTD,
recognizing that FDTD equations can be seen as a dynamical

Figure 11: SAR maps obtained in Sec. V-D with all-fine
(2 mm), all-coarse (1 cm) discretization and with subgridding
(2 mm in the tissues and 1 cm in the air).

system which is dissipative under a generalized CourantFriedrichs-Lewy condition. The theory provides a new, powerful framework to make FDTD stability analysis simpler
and modular. Stability conditions can indeed be given on
the individual components (e.g. boundary conditions, meshes,
thin-wire models) rather than on the whole coupled FDTD
setup. The theory is intuitive since rooted on the familiar
concept of energy dissipation, and sheds new light on the
root mechanisms behind FDTD instability. As an example of
application, a simple yet effective subgridding algorithm is
derived, with straightforward stability proof. The proposed algorithm allows material traverse, is simpler to implement than
existing solutions, and supports an arbitrary grid refinement
ratio. Numerical results confirm its stability and accuracy.
Speedups of up to 46X were observed with only 3.1% error
with respect to standard FDTD.
R EFERENCES
[1] A. Taflove and S. C. Hagness, Computational electrodynamics. Artech
house, 2005.
[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] 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.
[4] 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.
[5] 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.
[6] 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.
[7] W. Sui, D. A. Christensen, and C. H. Durney, “Extending the twodimensional FDTD method to hybrid electromagnetic systems with
active and passive lumped elements,” IEEE Trans. Microw. Theory
Techn., vol. 40, no. 4, pp. 724–730, 1992.
[8] 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.

JOURNAL OF LATEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012

[9] 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.
[10] X. Li, C. D. Sarris, and P. Triverio, “Structure-Preserving Reduction of
Finite-Difference Time-Domain Equations with Controllable Stability
Beyond the CFL Limit,” IEEE Trans. Microw. Theory Techn., vol. 62,
no. 12, pp. 3228–3238, 2014.
[11] X. Li and P. Triverio, “Stable FDTD Simulations with Subgridding at the
Time Step of the Coarse Grid: a Model Order Reduction Approach,” in
IEEE MTT-S Int. Conf. on Numerical Electromagnetic and Multiphysics
Modeling and Optimization, Ottawa, Canada, August 11–14 2015.
[12] J.-F. Lee, R. Lee, and A. Cangellaris, “Time-domain finite-element
methods,” IEEE Trans. Antennas Propag., vol. 45, no. 3, pp. 430–442,
1997.
[13] 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.
[14] Y. Wang, S. Safavi-Naeini, and S. K. Chaudhuri, “A hybrid technique
based on combining ray tracing and FDTD methods for site-specific
modeling of indoor radio wave propagation,” IEEE Trans. Antennas
Propag., vol. 48, no. 5, pp. 743–754, 2000.
[15] T. Namiki, “A new FDTD algorithm based on alternating-direction
implicit method,” IEEE Microw. Wireless Compon. Lett., vol. 47, no. 10,
pp. 2003–2007, 1999.
[16] F. Zheng, Z. Chen, and J. Zhang, “A Finite-Difference time-domain
method without the Courant stability conditions,” IEEE Microwave
Guided Wave Lett., vol. 9, no. 11, pp. 441–443, Nov 1999.
[17] P. Triverio, S. Grivet-Talocia, M. S. Nakhla, F. Canavero, R. Achar,
“Stability, causality, and passivity in electrical interconnect models,”
IEEE Trans. Adv. Packag., vol. 30, no. 4, pp. 795–808, 2007.
[18] J. C. Willems, “Dissipative dynamical systems part i: General theory,”
Archive for rational mechanics and analysis, vol. 45, no. 5, pp. 321–351,
1972.
[19] Y. Wang, S. Langdon, and C. Penney, “Analysis of accuracy and stability
of FDTD subgridding schemes,” in 2010 European Microwave Conf.
IEEE, 2010, pp. 1297–1300.
[20] F. Collino, T. Fouquet, and P. Joly, “Conservative space-time mesh
refinement methods for the FDTD solution of Maxwell’s equations,”
J. Comput. Phys., vol. 211, no. 1, pp. 9–35, 2006.
[21] R. A. Chilton and R. Lee, “Conservative and provably stable FDTD
subgridding,” IEEE Trans. Antennas Propag., vol. 55, no. 9, pp. 2537–
2549, 2007.
[22] 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.
[23] J. Brewer, “Kronecker products and matrix calculus in system theory,”
IEEE Trans. Circuits Syst., vol. 25, no. 9, pp. 772–781, 1978.
[24] B. Denecker, F. Olyslager, L. Knockaert, and D. De Zutter, “A new statespace-based algorithm to assess the stability of the Finite-Difference
time-domain method for 3D finite inhomogeneous problems,” AEU-Int.
J. Electron. C, vol. 58, no. 5, pp. 339 – 348, 2004.
[25] L. Dai, “Singular control systems,” 1989.
[26] 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.
[27] S. Grivet-Talocia and B. Gustavsen, Passive Macromodeling: Theory
and Applications. Hoboken, NJ: Wiley, 2015.
[28] 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.
[29] 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.
[30] A. Christ et al., “The Virtual Family - Development of surface-based
anatomical models of two adults and two children for dosimetric
simulations,” Physics in Medicine and Biology, vol. 55, no. 2, pp. N23–
N38, 2010.
[31] IT’IS Foundation, “Overview - Database of Tissue Properties.”
[Online]. Available: http://www.itis.ethz.ch/virtual-population/tissueproperties/overview/ [Aug. 19, 2015].
[32] S. Kühn et al., “MMF/GSMA phase 2: scientific basis for base station
exposure compliance standards,” IT’IS Foundation, Tech. Rep. 21, June
2009.
[33] D. A. Sánchez-Hernández, High Frequency Electromagnetic Dosimetry.
Boston, MA: Artech House, 2009.

12

