Relaxation of the Courant Condition in the Explicit Finite-Difference Time-Domain Method With Higher-Degree Differential Terms

概要

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

1

Relaxation of the Courant Condition in the Explicit
Finite-Difference Time-Domain Method with
Higher-Degree Differential Terms
Harune Sekido and Takayuki Umeda, Member, IEEE

Abstract—A new explicit and non-dissipative FDTD method
in two and three dimensions is proposed for relaxation of the
Courant condition. The third-degree spatial difference terms with
second- and fourth-order accuracy are added with coefficients to
the time-development equations of FDTD(2,4). Optimal coefficients are obtained by a brute-force search of the dispersion
relations, which reduces phase velocity errors but satisfies the
numerical stability as well. The new method is stable with large
Courant numbers where the conventional FDTD methods are
unstable. The new method also has smaller numerical errors in
the phase velocity than conventional FDTD methods with small
Courant numbers.
Index Terms—FDTD, Courant condition, dispersion relation,
phase velocity.

I. I NTRODUCTION
HE Finite-Difference Time-Domain (FDTD) method is a
numerical method for solving the time development of
electromagnetic fields [1], [2]. The time-development equations used in the FDTD method is obtained by approximating
Maxwell’s equations with the finite difference of second-order
accuracy in both time and space. A staggered grid (Yee grid)
system is used in the spatial difference so that Gauss’s law for
both electric and magnetic fields is always satisfied. Owing to
this advantage and the simpleness of the numerical algorithm,
the FDTD method has been used as the standard numerical
method for electromagnetic fields for more than a half century.
Maxwell’s equations are written as follows:

T

∂B
+∇×E =0
∂t
∂E
1
ϵ
− ∇ × B = −J
∂t
µ
ρ
∇·E =
ϵ
∇·B =0

(1a)
(1b)
(1c)
(1d)

where E is electric field, B is magnetic field, J is current
density, ϵ is permeability, µ is permittivity, and ρ is charge
density.
Assuming a case in vacuum (J = 0, ρ = 0), (1) are written
in the three vector components in the rectangular coordinate
system as follows:
−

∂Ez
∂Ey
∂Bx
=
−
∂t
∂y
∂z

(2a)

H. Sekido and T. Umeda are with Institute for Space-Earth Environmental
Research, Nagoya University (e-mail: sekido@isee.nagoya-u.ac.jp).

∂By
∂Ex
∂Ez
=
−
(2b)
∂t
∂z
∂x
∂Bz
∂Ey
∂Ex
−
=
−
(2c)
∂t
∂x
∂y


∂Ex
1 ∂Bz
∂By
ϵ
=
−
(2d)
∂t
µ ∂y
∂z


∂Ey
1 ∂Bx
∂Bz
ϵ
=
−
(2e)
∂t
µ
∂z
∂x


1 ∂By
∂Bx
∂Ez
=
−
.
(2f)
ϵ
∂t
µ ∂x
∂y
The time-development equations are derived by approximating temporal and spatial differential terms in (2) with the finite
difference of second-order accuracy:


∆z
∆y
t+ ∆t
Bx 2 x, y +
,z +
2
2


∆t
∆y
∆z
t− 2
x, y +
= Bx
,z +
2
2


(3a)
∆y
∆z
1 t
− Dy Ez x, y +
,z +
2
2


∆z
∆y
,z +
+ Dz1 Eyt x, y +
2
2




∆x
∆x
t+∆t
t
Ex
x+
, y, z = Ex x +
, y, z
2
2


∆x
t− ∆t
+c2 Dy1 Bz 2 x +
, y, z
(3b)
2


∆x
t− ∆t
−c2 Dz1 By 2 x +
, y, z
2


∆t
∆z
∆x
t+
By 2 x +
, y, z +
2
2


∆t
∆x
∆z
t− 2
= By
x+
, y, z +
2
2


(3c)
∆x
∆z
1 t
− D z Ex x +
, y, z +
2
2


∆z
∆x
, y, z +
+ Dx1 Ezt x +
2
2




∆y
∆y
t+∆t
t
Ey
x, y +
, z = Ey x, y +
,z
2
2


∆y
t− ∆t
+c2 Dz1 Bx 2 x, y +
,z
(3d)
2


∆y
t− ∆t
−c2 Dx1 Bz 2 x, y +
,z
2
−

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION


∆x
∆y
x+
,y +
,z
2
2


∆x
∆y
t− ∆t
2
= Bz
x+
,y +
,z
2
2


∆x
∆y
1 t
− D x Ey x +
,y +
,z
2
2


∆x
∆y
+ Dy1 Ext x +
,y +
,z
2
2



∆z
= Ezt x, y, z +
Ezt+∆t x, y, z +
2

∆t
2 1 t− 2
+c Dx By
x, y, z +

t− ∆t
−c2 Dy1 Bx 2 x, y, z +
t+ ∆t
2

2



Bz

(3e)



∆z
2

∆z
2

∆z
2

(3f)

where Dxn is a nth-degree spatial difference operator. Note
that Dzn = 0 in p
two dimensions. Here, the speed of light c is
defined as c = 1/ϵµ.
The first-degree spatial difference operator in the x direction
with the second-order accuracy Dx1 is defined as follows:


∆x
∆y
Dx1 Eyt x +
,y +
,z
2
2
 

∆t
∆y
(4)
=
Eyt x + ∆x, y +
,z
∆x
2


∆y
,z
.
−Eyt x, y +
2
The FDTD method has a disadvantage that numerical
oscillations occur in discontinuous waveforms and even in
continuous waveforms with a large slope. The numerical
oscillations are caused by an error in the phase velocity, which
depends on frequency or wavenumber. To reduce the numerical
error in the phase velocity, higher-order finite differences are
used. The FDTD(2,4) method uses the fourth-order spatial
difference [3], [4].
The first-degree spatial difference operator in the x direction
with the fourth-order accuracy Dx1 is defined as follows:


∆x
∆y
1 t
D x Ey x +
,y +
,z
2
2



1 ∆t
∆y
t
=
−Ey x + 2∆x, y +
,z
24 ∆x
2




∆y
∆y
+ 27Eyt x + ∆x, y +
, z − 27Eyt x, y +
,z
2
2


∆y
+Eyt x − ∆x, y +
,z
.
2
(5)
The FDTD method using tth-order accuracy in time and
xth-order accuracy in space is generally referred to as
FDTD(t,x).
The numerical error in the phase velocity of FDTD(2,4)
is smaller than that of FDTD(2,2). However, the Courant
condition of FDTD(2,4) method is more restricted than that
of FDTD(2,2). In general, higher-order finite differences in
space with the second-order finite difference in time make the
Courant condition more restrictive, which requires smaller ∆t
and larger number of time steps.

The Courant condition is derived from dispersion relation.
Considering only x direction in one dimension, the dispersion
relation of FDTD(2,2) is derived by Fourier transform of (3)
and (4) as follows:
 
2 

2
ω∆t
k∆x
sin
= C sin
(6)
2
2
where C = c∆t/∆x is the Courant number. In the same way,
the dispersion relation for FDTD(2,4) is derived from (3) with
(5) as follows:
2
 
ω∆t
sin
2
(7)
  


2
k∆x
1
k∆x
= C sin
+ sin3
.
2
6
2
The right-hand sides of (6) and (7) are both maximized at
k∆x = π:


ω∆t
2
sin
= C2
2

 
2
ω∆t
7
2
sin
=
C .
2
6
Therefore, the right-hand side of (7) is more than 1 for C >
6/7 ∼ 0.857, which causes a numerical instability.
Zhou et al. [5] used different operators to update electric
and magnetic field. The Zhou schemes have almost the same
numerical dispersion as the standard FDTD method but reduce
the computational costs. However, the Courant conditions of
the Zhou schemes are as restricted as the standard FDTD
method.
The FDTD modified(2,4) (M(2,4)) method [6] corrects the
phase velocity of FDTD(2,4) by adding a difference term
which uses diagonal grid points. The second-order nonstandard (NS)-FDTD method [7], [8] reduces anisotropic errors
in the phase velocity by considering diagonal differences.
The numerical phase velocity is matched with the physical
speed of light at a specific frequency only in this method.
Furthermore, the Wideband (W)NS-FDTD method aiming for
analysis in a wide frequency band corrects the phase velocity
error of NS-FDTD by a post process [9]. In these methods, the
phase velocity error is reduced by adding various difference
terms with coefficients which corrects numerical dispersion.
However, it is not easy to determine appropriate coefficients
for the correction terms.
The Crank-Nicolson (CN)-FDTD method [10]–[12] and the
Alternating Direction Implicit (ADI)-FDTD method [13]–[16]
use implicit time-development equations which relaxes the
Courant condition. However, the implicit equations have larger
computational costs which are solved with the matrix inversion
or iterative convergence.
This paper aims to relax the Courant condition of the
FDTD(2,4) method in two and three dimensions. A new
explicit and non-dissipative method is developed by adding
third-degree spatial difference terms with coefficients to the
time-development equations of FDTD(2,4). The coefficients
are determined so that the phase velocity error is minimized
but numerical instabilities are suppressed as well.

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION

3

This paper is organized as follows. In section 2, the timedevelopment equations and the dispersion relation of the new
methods are shown. In section 3 and 4, the optimal coefficients, the phase velocity errors, and the results of numerical
tests using the new methods in two and three dimensions are
shown, respectively. In section 5, the conclusion is given.

II. F ORMULATION AND N UMERICAL D ISPERSION
R ELATION





∆y
∆y
Eyt+∆t x, y +
, z = Eyt x, y +
,z
2
2


∆t
∆y
t−
+ c2 Dz1 Bx 2 x, y +
,z
2


∆t
∆y
2 3 t− 2
+ αc Dz Bx
x, y +
,z
2


∆y
t− ∆t
− c2 Dx1 Bz 2 x, y +
,z
2


∆t
∆y
t−
,z
− αc2 Dx3 Bz 2 x, y +
2

(8d)

A. General Form
The time-development equations with third-degree spatial
difference terms are written as follows:


∆y
∆z
,z +
x, y +
2
2


∆y
∆z
t− ∆t
2
= Bx
x, y +
,z +
2
2


∆z
∆y
,z +
− Dy1 Ezt x, y +
2
2


∆z
∆y
− αDy3 Ezt x, y +
,z +
2
2


∆y
∆z
1 t
+ Dz Ey x, y +
,z +
2
2


∆z
∆y
,z +
+ αDz3 Eyt x, y +
2
2
t+ ∆t
2

Bx




∆x
∆x
t
, y, z = Ex x +
, y, z
x+
2
2


∆x
t− ∆t
+ c2 Dy1 Bz 2 x +
, y, z
2


∆t
∆x
2
3 t− 2
, y, z
+ c αDy Bz
x+
2


∆t
∆x
2 1 t− 2
− c Dz B y
x+
, y, z
2


∆x
t− ∆t
− c2 αDz3 By 2 x +
, y, z
2


∆y
∆x
,y +
,z
x+
2
2


∆x
∆y
t− ∆t
2
= Bz
x+
,y +
,z
2
2


∆x
∆y
− Dx1 Eyt x +
,y +
,z
2
2


∆y
∆x
,y +
,z
− αDx3 Eyt x +
2
2


∆x
∆y
1 t
+ D y Ex x +
,y +
,z
2
2


∆y
∆x
3 t
+ αDy Ex x +
,y +
,z
2
2

t+ ∆t
Bz 2



(8a)



(8e)



Ext+∆t




∆x
∆z
x+
, y, z +
2
2


∆t
∆x
∆z
t−
, y, z +
= By 2 x +
2
2


∆x
∆z
1 t
− D z Ex x +
, y, z +
2
2


∆x
∆z
3 t
− αDz Ex x +
, y, z +
2
2


∆x
∆z
+ Dx1 Ezt x +
, y, z +
2
2


∆x
∆z
+ αDx3 Ezt x +
, y, z +
2
2
t+ ∆t
2

(8b)





∆z
∆z
Ezt+∆t x, y, z +
= Ezt x, y, z +
2
2


∆t
∆z
t−
+ c2 Dx1 By 2 x, y, z +
2


∆t
∆z
t−
2 3
2
+ αc Dx By
x, y, z +
2


∆t
∆z
t−
2 1
2
− c Dy B x
x, y, z +
2


∆t
∆z
t−
− αc2 Dy3 Bx 2 x, y, z +
2

(8f)

By

(8c)

where a coefficient α is added to the third-degree difference
terms. The optimal coefficient α depends on the Courant
number C. The time-development (8) is based on the Taylor
series of the central difference equation in time. In the present
study, odd-degree difference terms are added only, because
even-degree difference terms lead to a numerical dissipation. ...