# Frequency-Dependent FDTD Algorithm Using Newmark’s Method

Authors: Bing Wei; Le Cao; Fei Wang; Qian Yang
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

According to the characteristics of the polarizability in frequency domain of three common models of dispersive media, the relation between the polarization vector and electric field intensity is converted into a time domain differential equation of second order with the polarization vector by using the conversion from frequency to time domain. Newmarkβ
γ difference method is employed to solve this equation. The electric field intensity to polarizability recursion is derived, and the electric flux to electric field intensity recursion is obtained by constitutive relation. Then FDTD iterative computation in time domain of electric and magnetic field components in dispersive medium is completed. By analyzing the  solution stability of the above differential equation using central difference method, it is proved that this method has more advantages in the selection of time step. Theoretical analyses and numerical results demonstrate that this method is a general algorithm and it has advantages of higher accuracy and stability over the algorithms based on central difference method.


---

## Body

## 1. Introduction

In 1966, the finite difference time domain (FDTD) method [1] is put forward by Yee. With the further development, the FDTD method probably becomes the most popular numerical technique for solving complex electromagnetic problems [1–8]. FDTD is a time domain algorithm but the constitutive relation of dispersive medium is given by a frequency domain expression. Hence, it is difficult to use FDTD method directly in analyzing electromagnetic characteristics of dispersive medium. According to the different treatments of constitutive relation, the FDTD methods used in dispersive medium can be classified as RC method [5, 6], JEC method [7], PLJERC method [8], ADE method [9], Z transform method [10–12], SO method [13–15], SARC method [16], and so on.In the case of nondispersive medium, the spatial discretization of Maxwell curl equation can be determined by the requirement from numerical dispersion to stability, and the time discrete can be determined by Courant stability condition. So the calculation involves just the first order partial differential equation. The time domain expression of constitutive relation often can be written as a second order partial differential equation in the case of dispersive medium. To obtain an accurate solution of the partial differential equation, time step must be very small. In this paper, the general differential equation of dispersive medium is given and the Newmark method [17–19] is introduced to the treatment of the constitutive relation in dispersive medium. This new method has a higher accuracy and stability than the central difference method. The detailed process is as follows.According to the polarizability characteristics in frequency domain of three common models (Debye, Drude, and Lorentz models) of dispersive medium, the relation in frequency domain betweenP and E is converted into a second order differential equation (has a general form to the three common models) by using the conversion j
ω
→
∂
/
∂
t from frequency to time domain. Then the Newmark β
γ method which is used in FETD is employed to solve the equation. Hence the E
→
P recursion is derived, and the D
→
E recursion is obtained by constitutive relation.


## 2. The General Newmark-FDTD Method

### 2.1. Maxwell Equation in Dispersive Medium and Constitutive Relation

Maxwell curl equation in linear isotropic dispersive medium is(1)
∂
D
∂
t
=
∇
×
H
,
∂
H
∂
t
=
-
1
μ
0
∇
×
E
.
The electric constitutive relation is
(2)
D
(
ω
)
=
ε
(
ω
)
E
(
ω
)
=
ε
0
[
ε
∞
+
χ
(
ω
)
]
E
(
ω
)
,
where ε
(
ω
) and χ
(
ω
) are permittivity coefficient and polarizability in frequency domain, respectively. ε
∞ is the relative permittivity when frequency is infinite.The FDTD iterative computation ofE
→
H
→
D can be accomplished by the discretization of (1) according to the standing Yee cell. And it is needed to deal with constitutive relation in the iterative computation of D
→
E.


### 2.2. The United Expression of Polarizability Function of Common Dispersive Medium Models

Three common dispersive medium models are Debye model, Drude model, and Lorentz model [2]. Their polarizabilities can be written as [20]
(3)
χ
Debye
(
ω
)
=
∑
l
=
1
L
χ
l
Debye
(
ω
)
=
∑
l
=
1
L
Δ
ε
l
1
+
j
ω
τ
0
,
l
,
χ
Drude
(
ω
)
=
∑
l
=
1
L
χ
l
Drude
(
ω
)
=
∑
l
=
1
L
-
ω
l
2
ω
2
-
j
ω
ν
c
,
l
,
χ
Lorentz
(
ω
)
=
∑
l
=
1
L
χ
l
Lorentz
(
ω
)
=
∑
l
=
1
L
Δ
ε
l
ω
0
,
l
2
ω
0
,
l
2
+
2
j
ω
ν
c
,
l
-
ω
2
,
where L is numbers of poles, l represents the lth pole, Δ
ε
l
=
ε
s
,
l
-
ε
∞, ε
s
,
l is the relative permittivity at zero frequency, τ
0
,
l is the relaxation time of pole, ν
c
,
l is collision frequency, ω
l is Drude frequency, and ω
0
,
l is the natural frequency of oscillator.Let(4)
ξ
l
(
ω
)
=
1
χ
l
(
ω
)
.
ξ
l
(
ω
) of the three dispersive medium models in (3) can be written as a united form:
(5)
ξ
l
(
ω
)
=
M
l
·
(
j
ω
)
2
+
C
l
·
(
j
ω
)
+
K
l
,
where
(6)
M
l
=
0
,
C
l
=
τ
0
,
l
Δ
ε
l
,
K
l
=
1
Δ
ε
l
Debye
M
l
=
1
ω
l
2
,
C
l
=
ν
c
,
l
ω
l
2
,
K
l
=
0
Drude
M
l
=
1
Δ
ε
l
ω
0
,
l
2
,
C
l
=
2
ν
c
,
l
Δ
ε
l
ω
0
,
l
2
,
K
l
=
1
Δ
ε
l
Lorentz
.


### 2.3. Constitutive Relations in Time Domain and Its Newmark’s Solution

Polarization vectorP is introduced:
(7)
P
(
ω
)
=
χ
(
ω
)
E
(
ω
)
.
The constitutive relation equation (2) can be rewritten as
(8)
D
(
ω
)
=
ε
0
ε
∞
E
(
ω
)
+
ε
0
P
(
ω
)
.
Let
(9)
P
l
(
ω
)
=
χ
l
(
ω
)
E
(
ω
)
.
Then from (3) and (7) we know that
(10)
P
(
ω
)
=
∑
l
=
1
L
P
l
(
ω
)
.
From (4) and (9), we have
(11)
E
(
ω
)
=
ξ
l
(
ω
)
P
l
(
ω
)
.
After substituting (5) into (11) and using conversion j
ω
→
∂
/
∂
t from frequency to time domain, we have
(12)
E
(
t
)
=
M
l
d
2
d
t
2
P
l
(
t
)
+
C
l
d
d
t
P
l
(
t
)
+
K
l
P
l
(
t
)
.
E
(
t
) in the above differential equation with second order time derivative is equivalent to a excitation source. In the year of 1959, a time-stepping algorithm for solving the differential equations like (12) is given by Newmark according to the time derivative difference approximate solution of Taylor series and the mean value theorem. This algorithm provides higher accuracy and stability [17]. By the deriving method of Zienkiewicz, solution of (12) is [18]
(13)
P
l
n
+
1
=
w
1
,
l
P
l
n
+
w
2
,
l
P
l
n
-
1
+
u
0
,
l
E
n
+
1
+
u
1
,
l
E
n
+
u
2
,
l
E
n
-
1
,
where
(14)
w
1
,
l
=
2
M
l
-
(
1
-
2
γ
)
Δ
t
C
l
-
(
(
1
/
2
)
+
γ
-
2
β
)
(
Δ
t
)
2
K
l
S
l
,
w
2
,
l
=
-
M
l
-
(
γ
-
1
)
Δ
t
C
l
-
(
(
1
/
2
)
-
γ
+
β
)
(
Δ
t
)
2
K
l
S
l
,
u
0
,
l
=
(
Δ
t
)
2
β
S
l
,
u
1
,
l
=
(
Δ
t
)
2
(
(
1
/
2
)
+
γ
-
2
β
)
S
l
,
u
2
,
l
=
(
Δ
t
)
2
(
(
1
/
2
)
-
γ
+
β
)
S
l
,
S
l
=
M
l
+
γ
Δ
t
C
l
+
β
(
Δ
t
)
2
K
l
.
Equation (13) is called two-step or Newmark β
γ method [19]. The calculation of P
l
n
+
1 involves the values of P
l
n and P
l
n
-
1. The coefficients in (13) include two variables β and γ. The stability and accuracy of (12) when selecting different β and γ are also discussed in [19] in detail. Usually, 0
≤
γ
≤
1 and 0
≤
β
≤
1
/
2. To assure the accuracy, it is needed to choose an appropriate γ and β
≥
0.25
(
0.5
+
γ
)
2. According to [19], unconditional convergence can be obtained by taking β
=
0.25, γ
=
0.5. The time-stepping calculation from E
n
+
1 to P
l
n
+
1 can be finished by (13).After converting the constitutive relation equation (8) into discrete time domain, we have
(15)
D
n
+
1
(
t
)
=
ε
0
(
ε
∞
E
n
+
1
(
t
)
+
∑
l
=
1
L
P
l
n
+
1
(
t
)
)
.
Substituting (13) into (15), we know that
(16)
(
ε
∞
+
∑
l
=
1
L
u
0
,
l
)
E
n
+
1
=
D
n
+
1
ε
0
-
∑
l
=
1
L
(
w
1
,
l
P
l
n
)
-
∑
l
=
1
L
(
w
2
,
l
P
l
n
-
1
)
-
∑
l
=
1
L
u
1
,
l
·
E
n
-
∑
l
=
1
L
u
2
,
l
·
E
n
-
1
.
Then the time-stepping calculation from D
n
+
1, P
l
n to E
n
+
1 can be accomplished.In conclusion, the general Newmark-FDTD method in dispersive medium can be reduced as follows.(1)
Use the difference discrete expressions in curl (1) to get E
→
H
→
D
.
(2)
Use (16) to get D
,
P
→
E
.
(3)
Use (13) to get E
→
P
.
(4)
Back to step 1.


## 2.1. Maxwell Equation in Dispersive Medium and Constitutive Relation

Maxwell curl equation in linear isotropic dispersive medium is(1)
∂
D
∂
t
=
∇
×
H
,
∂
H
∂
t
=
-
1
μ
0
∇
×
E
.
The electric constitutive relation is
(2)
D
(
ω
)
=
ε
(
ω
)
E
(
ω
)
=
ε
0
[
ε
∞
+
χ
(
ω
)
]
E
(
ω
)
,
where ε
(
ω
) and χ
(
ω
) are permittivity coefficient and polarizability in frequency domain, respectively. ε
∞ is the relative permittivity when frequency is infinite.The FDTD iterative computation ofE
→
H
→
D can be accomplished by the discretization of (1) according to the standing Yee cell. And it is needed to deal with constitutive relation in the iterative computation of D
→
E.


## 2.2. The United Expression of Polarizability Function of Common Dispersive Medium Models

Three common dispersive medium models are Debye model, Drude model, and Lorentz model [2]. Their polarizabilities can be written as [20]
(3)
χ
Debye
(
ω
)
=
∑
l
=
1
L
χ
l
Debye
(
ω
)
=
∑
l
=
1
L
Δ
ε
l
1
+
j
ω
τ
0
,
l
,
χ
Drude
(
ω
)
=
∑
l
=
1
L
χ
l
Drude
(
ω
)
=
∑
l
=
1
L
-
ω
l
2
ω
2
-
j
ω
ν
c
,
l
,
χ
Lorentz
(
ω
)
=
∑
l
=
1
L
χ
l
Lorentz
(
ω
)
=
∑
l
=
1
L
Δ
ε
l
ω
0
,
l
2
ω
0
,
l
2
+
2
j
ω
ν
c
,
l
-
ω
2
,
where L is numbers of poles, l represents the lth pole, Δ
ε
l
=
ε
s
,
l
-
ε
∞, ε
s
,
l is the relative permittivity at zero frequency, τ
0
,
l is the relaxation time of pole, ν
c
,
l is collision frequency, ω
l is Drude frequency, and ω
0
,
l is the natural frequency of oscillator.Let(4)
ξ
l
(
ω
)
=
1
χ
l
(
ω
)
.
ξ
l
(
ω
) of the three dispersive medium models in (3) can be written as a united form:
(5)
ξ
l
(
ω
)
=
M
l
·
(
j
ω
)
2
+
C
l
·
(
j
ω
)
+
K
l
,
where
(6)
M
l
=
0
,
C
l
=
τ
0
,
l
Δ
ε
l
,
K
l
=
1
Δ
ε
l
Debye
M
l
=
1
ω
l
2
,
C
l
=
ν
c
,
l
ω
l
2
,
K
l
=
0
Drude
M
l
=
1
Δ
ε
l
ω
0
,
l
2
,
C
l
=
2
ν
c
,
l
Δ
ε
l
ω
0
,
l
2
,
K
l
=
1
Δ
ε
l
Lorentz
.


## 2.3. Constitutive Relations in Time Domain and Its Newmark’s Solution

Polarization vectorP is introduced:
(7)
P
(
ω
)
=
χ
(
ω
)
E
(
ω
)
.
The constitutive relation equation (2) can be rewritten as
(8)
D
(
ω
)
=
ε
0
ε
∞
E
(
ω
)
+
ε
0
P
(
ω
)
.
Let
(9)
P
l
(
ω
)
=
χ
l
(
ω
)
E
(
ω
)
.
Then from (3) and (7) we know that
(10)
P
(
ω
)
=
∑
l
=
1
L
P
l
(
ω
)
.
From (4) and (9), we have
(11)
E
(
ω
)
=
ξ
l
(
ω
)
P
l
(
ω
)
.
After substituting (5) into (11) and using conversion j
ω
→
∂
/
∂
t from frequency to time domain, we have
(12)
E
(
t
)
=
M
l
d
2
d
t
2
P
l
(
t
)
+
C
l
d
d
t
P
l
(
t
)
+
K
l
P
l
(
t
)
.
E
(
t
) in the above differential equation with second order time derivative is equivalent to a excitation source. In the year of 1959, a time-stepping algorithm for solving the differential equations like (12) is given by Newmark according to the time derivative difference approximate solution of Taylor series and the mean value theorem. This algorithm provides higher accuracy and stability [17]. By the deriving method of Zienkiewicz, solution of (12) is [18]
(13)
P
l
n
+
1
=
w
1
,
l
P
l
n
+
w
2
,
l
P
l
n
-
1
+
u
0
,
l
E
n
+
1
+
u
1
,
l
E
n
+
u
2
,
l
E
n
-
1
,
where
(14)
w
1
,
l
=
2
M
l
-
(
1
-
2
γ
)
Δ
t
C
l
-
(
(
1
/
2
)
+
γ
-
2
β
)
(
Δ
t
)
2
K
l
S
l
,
w
2
,
l
=
-
M
l
-
(
γ
-
1
)
Δ
t
C
l
-
(
(
1
/
2
)
-
γ
+
β
)
(
Δ
t
)
2
K
l
S
l
,
u
0
,
l
=
(
Δ
t
)
2
β
S
l
,
u
1
,
l
=
(
Δ
t
)
2
(
(
1
/
2
)
+
γ
-
2
β
)
S
l
,
u
2
,
l
=
(
Δ
t
)
2
(
(
1
/
2
)
-
γ
+
β
)
S
l
,
S
l
=
M
l
+
γ
Δ
t
C
l
+
β
(
Δ
t
)
2
K
l
.
Equation (13) is called two-step or Newmark β
γ method [19]. The calculation of P
l
n
+
1 involves the values of P
l
n and P
l
n
-
1. The coefficients in (13) include two variables β and γ. The stability and accuracy of (12) when selecting different β and γ are also discussed in [19] in detail. Usually, 0
≤
γ
≤
1 and 0
≤
β
≤
1
/
2. To assure the accuracy, it is needed to choose an appropriate γ and β
≥
0.25
(
0.5
+
γ
)
2. According to [19], unconditional convergence can be obtained by taking β
=
0.25, γ
=
0.5. The time-stepping calculation from E
n
+
1 to P
l
n
+
1 can be finished by (13).After converting the constitutive relation equation (8) into discrete time domain, we have
(15)
D
n
+
1
(
t
)
=
ε
0
(
ε
∞
E
n
+
1
(
t
)
+
∑
l
=
1
L
P
l
n
+
1
(
t
)
)
.
Substituting (13) into (15), we know that
(16)
(
ε
∞
+
∑
l
=
1
L
u
0
,
l
)
E
n
+
1
=
D
n
+
1
ε
0
-
∑
l
=
1
L
(
w
1
,
l
P
l
n
)
-
∑
l
=
1
L
(
w
2
,
l
P
l
n
-
1
)
-
∑
l
=
1
L
u
1
,
l
·
E
n
-
∑
l
=
1
L
u
2
,
l
·
E
n
-
1
.
Then the time-stepping calculation from D
n
+
1, P
l
n to E
n
+
1 can be accomplished.In conclusion, the general Newmark-FDTD method in dispersive medium can be reduced as follows.(1)
Use the difference discrete expressions in curl (1) to get E
→
H
→
D
.
(2)
Use (16) to get D
,
P
→
E
.
(3)
Use (13) to get E
→
P
.
(4)
Back to step 1.


## 3. The Advantages of Algorithm in This Paper

For the second order differential equations like(17)
M
d
2
d
t
2
ϕ
+
C
d
d
t
ϕ
+
K
ϕ
=
0
.
After using central difference method, we have
(18)
(
M
+
C
2
Δ
t
)
ϕ
n
+
1
-
[
2
M
-
K
(
Δ
t
)
2
]
ϕ
n
+
[
M
-
C
2
Δ
t
]
ϕ
n
-
1
=
0
.The time domain solution stability of the above differential equation can be analyzed by using theZ transform. To a time series of x
(
n
T
)
(
n
=
0,1
,
2
.
.
.
), its Z transform is defined as
(19)
X
(
z
)
=
∑
n
=
0
∞
x
(
n
T
)
z
-
n
.
X
(
z
)  is the Z transform of the sampling sequence x
(
n
T
).Z transform is proposed to (18); we know that
(20)
(
M
+
C
2
Δ
t
)
Z
2
-
[
2
M
-
K
(
Δ
t
)
2
]
Z
+
[
M
-
C
2
Δ
t
]
=
0
.According to the final value theorem, time domain functionx
(
n
T
)
(
n
=
0,1
,
2
.
.
.
) has time stable solutions when the poles of X
(
z
)  located in a unit circle on complex plane in Z domain. An analysis of the solution stability of second order homogeneous differential equation without first order derivative is made in [21]. In this paper, the case of including first order derivative is analyzed. Equation (20) is considered as a quadratic equation with one unknown using Z as dependent variable. The equation can be solved and let -
1
≤
Z
≤
1; we know that
(21)
K
C
(
Δ
t
)
3
+
2
M
K
(
Δ
t
)
2
-
4
M
C
Δ
t
-
8
M
2
≤
0
.The corresponding cubic equation with one unknown to the above inequality has three solutions:(22)
Δ
t
1
=
-
2
M
C
,
Δ
t
2
=
-
2
M
K
,
Δ
t
3
=
2
M
K
.LetM
,
C
,
K be all positive quantities. The derivative of the cubic equation with one unknown is 8
M
C
+
8
M
M
K at Δ
t
3, which is an increasing function satisfying the requirements of the inequality. Hence, when central difference method is used to solve (17), the time-stepping is limited to
(23)
Δ
t
≤
2
M
K
.
According to (23), the requirements of time-stepping when central difference method is used to solve the second order time domain differential equations corresponding to the three models of dispersive medium in (6) are
(24)
Δ
t
Debye
≤
0
,
Δ
t
Drude
≤
∞
,
Δ
t
Lorentz
≤
2
ω
0
,
l
.
That is to say, for Debye medium, corrected results cannot be obtained whatever the value of the time-stepping takes when central difference method is used to solve (17). For Drude medium, corrected results can be obtained whatever the value of the time-stepping takes. But for Lorentz medium, the solution stability is relevant to its resonance frequency at poles.However, the solution of (12) is unconditionally stable when using Newmark difference algorithm. Hence, the selection of time-stepping is more flexible. Meanwhile, a higher accuracy is obtained because the Newmark β
γ method has a combination of Taylor expansions and means value theorem is employed.To analyze the difference between using Newmarkβ
γ method and central difference method, the results obtained from these two methods and Runge-Kutta solutions [20] are given. Figure 1 is the absolute error comparing the classical Runge-Kutta solutions of using these two methods, in which the solid line and the circle are the result from Newmark β
γ and central difference method, respectively. As it is shown in the figure, the time it takes to get steady state for Newmark method is shorter than central difference method. Let
(25)
K
=
7.96
×
1
0
7
,
M
=
2.95
×
1
0
-
14
,
C
=
2.65
×
1
0
-
3
,
f
=
sin
(
4
×
1
0
9
t
)
in (17).Figure 1
Comparison of stability and calculation error between central difference and Newmark method.The theoretical analysis and numerical results show that the algorithm presented in this paper has higher steady state and calculation accuracy than central difference method.


## 4. Numerical Results

To demonstrate the correctness of our algorithm, the backscatter of three different dispersive media sphere is given and compared with Mie series solution. Letβ
=
0.25, γ
=
0.5, Δ
t
=
δ
/
2
c, where δ is spatial discrete grid-scale of Yee cell. Consider
(26)
E
i
(
t
)
=
exp
⁡
[
-
4
π
(
t
-
t
0
)
2
τ
2
]
,
where τ
=
60
Δ
t and t
0 is the pulse peak time.Example 1 (the backscatter of Debye sphere).
Let the radius of the sphere be 0.25 m and filled with single pole Debye medium. The complex relative permittivity is(27)
ε
r
(
ω
)
=
ε
∞
+
χ
(
ω
)
=
ε
∞
+
ε
s
,
1
-
ε
∞
1
+
j
ω
τ
0,1
,
where ε
s
,
1
=
1.16, ε
∞
=
1.01, τ
0,1
=
6.497
×
1
0
-
10 s. The backscatter RCS is shown in Figure 2. The solid line and circle are the results from this method and Mie series solution, respectively. It is shown that these two results are in good agreement with each other. Consider δ
=
3.3
×
1
0
-
3 min FDTD calculation.Figure 2
The backscatter RCS of Debye sphere.Example 2 (the backscatter of plasma sphere).
The radius of the sphere is 3.75 mm and filled with single pole Drude medium. The complex relative permittivity is(28)
ε
r
(
ω
)
=
ε
∞
+
χ
(
ω
)
=
1
+
-
ω
1
2
ω
2
-
j
ω
ν
c
,
1
,
where ω
1
=
1.8
×
1
0
11 rad/s, ν
c
,
1
=
2.0
×
1
0
10 Hz. The backscatter RCS is shown in Figure 3. The solid line and circle are the result from this method and Mie series solution, respectively. It is shown that these two results are in good agreement with each other. Consider δ
=
5.0
×
1
0
-
2 mm in FDTD calculation.Figure 3
The backscatter RCS of plasma sphere.Example 3 (the backscatter of Lorentz sphere).
The radius of the sphere is15.0
×
1
0
-
9 m and filled with single pole Lorentz medium. The complex relative permittivity is
(29)
ε
r
(
ω
)
=
ε
∞
+
χ
(
ω
)
=
ε
∞
+
(
ε
s
,
1
-
ε
∞
)
ω
0,1
2
ω
0,1
2
+
2
j
ω
ν
c
,
1
-
ω
2
,
where ε
s
,
1
=
2.25, ε
∞
=
1.0, ν
c
,
1
=
0.28
×
1
0
16 Hz, ω
0,1
=
4.0
×
1
0
16 Hz. The backscatter is shown in Figure 4. The solid line and circle are the result from this method and Mie series solution, respectively. It is shown that these two results are in good agreement with each other. Consider δ
=
3.0
×
1
0
-
10m in FDTD calculation.Figure 4
The backscatter RCS of Lorentz sphere.


## 5. Conclusions

In this paper, a general Newmark-FDTD algorithm is given to deal with the electromagnetic problems in dispersive medium. This new method combines the Newmarkβ
γ difference method which is widely used in FETD calculation in dispersive medium. Theoretical analyses and numerical results demonstrate that this algorithm has advantages of higher accuracy and stability over the algorithms based on central difference method.



---

*Source: 216763-2014-10-13.xml*