# Modified Splitting FDTD Methods for Two-Dimensional Maxwell’s Equations

Authors: Liping Gao; Shouhui Zhai
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

In this paper, we develop a new method to reduce the error in the splitting finite-difference method of Maxwell’s equations. By this method two modified splitting FDTD methods (MS-FDTDI, MS-FDTDII) for the two-dimensional Maxwell equations are proposed. It is shown that the two methods are second-order accurate in time and space and unconditionally stable by Fourier methods. By energy method, it is proved that MS-FDTDI is second-order convergent. By deriving the numerical dispersion (ND) relations, we prove rigorously that MS-FDTDI has less ND errors than the ADI-FDTD method and the ND errors of ADI-FDTD are less than those of MS-FDTDII. Numerical experiments for computing ND errors and simulating a wave guide problem and a scattering problem are carried out and the efficiency of the MS-FDTDI and MS-FDTDII methods is confirmed.


---

## Body

## 1. Introduction

The finite-difference time-domain (FDTD) method for Maxwell’s equations, which was first proposed by Yee (see [1], also called Yee’s scheme) in 1966, is a very efficient numerical algorithm in computational electromagnetism (see [2]) and has been applied in a broad range of practical problems by combining absorbing boundary conditions (see [3–7] and the references therein). It is well known from [8] that the Yee Scheme is stable when time and spatial step sizes (Δ
t, Δ
x, and Δ
y for 2D case) satisfy the Courant-Friedrichs-Lewy (CFL) condition c
Δ
t
≤
[
1
/
(
Δ
x
)
2
+
1
/
(
Δ
y
)
2
]
-
1
/
2, where c is the wave velocity. To overcome the restriction of the CFL condition there are many research works on this topic; for example, see [9–17] and the references therein. In [15], two unconditionally stable FDTD methods (named as S-FDTDI and S-FDTDII) were proposed by using splitting of the Maxwell equations and reducing of the perturbation error, where S-FDTDII, based on S-FDTDI (first-order accurate), is second-order accurate and has less numerical dispersion (ND) error than S-FDTDI. However, the second convergence of S-FDTDII was not proved by the energy method.In this letter, by introducing a new method to reduce the error caused by splitting of equations [15] (other methods of reducing perturbation error caused by splitting of differential equations can be seen in [18]), we propose two modified splitting FDTD methods (called MS-FDTDI and MS-FDTDII) for the 2D Maxwell equations. It is proved by the energy method that MS-FDTDI with the perfectly electric conducting boundary conditions is second-order convergent in both time and space. By Fourier method we derive the amplification factors and ND relations of MS-FDTDI and MS-FDTDII. Then, we prove that these two methods are unconditionally stable and that MS-FDTDI has less ND errors than S-FDTDII (or ADI-FDTD [10, 11]). Numerical experiments to compute numerical dispersion errors and convergence orders and to simulate a scattering problem are carried out. Computational results confirm the analysis of MS-FDTDI and MS-FDTDII.


## 2. Modified Splitting FDTD Method for the Maxwell Equations

### 2.1. Maxwell Equations

Consider the two-dimensional Maxwell equations in a lossless and homogeneous medium:(1)
∂
E
x
∂
t
=
1
ε
∂
H
z
∂
y
,
∂
E
y
∂
t
=
-
1
ε
∂
H
z
∂
x
,
∂
H
z
∂
t
=
1
μ
∂
E
x
∂
y
-
∂
E
y
∂
x
,where ε and μ are the electric permittivity and magnetic permeability of the medium and E
=
(
E
x
(
x
,
y
,
t
)
,
E
y
(
x
,
y
,
t
)
) and H
z
=
H
z
(
x
,
y
,
t
) for (
x
,
y
)
∈
Ω
=
[
0
,
a
]
×
[
0
,
b
] and t
∈
0
,
T denote the electric and magnetic fields, respectively. We assume that the spatial domain Ω is surrounded by perfectly electric conductor (PEC). Then the PEC boundary condition below is satisfied:(2)
E
,
0
×
n
→
,
0
=
0
on
0
,
T
×
∂
Ω
,where ∂
Ω denotes the boundary of Ω and n
→ is the outward normal vector on ∂
Ω. The initial conditions are assumed to be (3)
E
x
,
y
,
0
=
E
0
x
,
y
,
H
z
x
,
y
,
0
=
H
z
0
x
,
y
,where E
0
(
x
,
y
)
=
(
E
x
0
(
x
,
y
)
,
E
y
0
(
x
,
y
)
).


### 2.2. Partition of the Domains and Notations

LetΩ be partitioned as Yee’s staggered grids [1]: {
(
x
α
,
y
β
)
∣
α
=
i
,
i
+
1
/
2
,
β
=
j
,
j
+
1
/
2
}, and let [
0
,
T
] be divided into equidistant subintervals, [
t
n
,
t
n
+
1
], where(4)
x
i
=
i
Δ
x
,
x
i
+
1
/
2
=
x
i
+
1
2
Δ
x
,
x
I
=
I
Δ
x
=
a
,
y
j
=
j
Δ
y
,
y
j
+
1
/
2
=
j
+
1
2
Δ
y
,
y
J
=
J
Δ
y
=
b
,
t
n
=
n
Δ
t
,
t
n
+
1
/
2
=
t
n
+
1
2
Δ
t
,
t
N
=
N
Δ
t
=
T
,where Δ
x and Δ
y are the spatial step sizes, Δ
t is the time increment, and I, J, and N are positive integers. For a function F
(
x
,
y
,
t
) and u
,
v
=
x or y, we define(5)
F
α
,
β
m
=
F
α
Δ
x
,
β
Δ
y
,
m
Δ
t
,
δ
x
F
α
,
β
m
=
F
α
+
1
/
2
,
β
m
-
F
α
-
1
/
2
,
β
m
Δ
x
,
δ
y
F
α
,
β
m
=
F
α
,
β
+
1
/
2
m
-
F
α
,
β
-
1
/
2
m
Δ
y
,
δ
u
δ
v
F
α
,
β
m
=
δ
u
δ
v
F
α
,
β
m
.


### 2.3. Modified Splitting FDTD Methods

Denote by(
E
x
i
¯
,
j
m
,
E
y
i
,
j
¯
m
) and H
z
i
¯
,
j
¯
m the approximations to (
E
x
(
x
i
¯
,
y
j
,
t
m
), E
y
(
x
i
,
y
j
¯
,
t
m
)
) and H
z
(
x
i
¯
,
y
j
¯
,
t
m
), respectively, where, and in what follows, i
¯
=
i
+
1
/
2, j
¯
=
j
+
1
/
2. Based on the S-FDTDII scheme (see [15]) and the idea of reducing the splitting error, we propose a modified splitting FDTD method (called MS-FDTDI) for (1)–(3).Stage 1.
(6)
E
y
i
,
j
¯
n
+
1
-
E
y
i
,
j
¯
n
Δ
t
=
-
1
2
ε
δ
x
H
z
i
,
j
¯
∗
+
H
z
i
,
j
¯
n
-
Δ
t
2
μ
ε
δ
x
δ
y
E
x
i
,
j
¯
n
;
H
z
i
¯
,
j
¯
∗
-
H
z
i
¯
,
j
¯
n
Δ
t
=
-
1
2
μ
δ
x
E
y
i
¯
,
j
¯
n
+
1
+
E
y
i
¯
,
j
¯
n
.Stage 2.
(7)
E
x
i
¯
,
j
n
+
1
-
E
x
i
¯
,
j
n
Δ
t
=
1
2
ε
δ
y
H
z
i
¯
,
j
n
+
1
+
H
z
i
¯
,
j
n
+
Δ
t
4
μ
ε
δ
x
δ
y
E
y
i
¯
,
j
n
+
1
-
E
y
i
¯
,
j
n
;
(8)
H
z
i
¯
,
j
¯
n
+
1
-
H
z
i
¯
,
j
¯
∗
Δ
t
=
1
2
μ
δ
y
E
x
i
¯
,
j
¯
n
+
1
+
E
x
i
¯
,
j
¯
n
.The boundary conditions for (6)–(8) obtained from (2) are(9)
E
x
i
+
1
/
2,0
m
=
E
x
i
+
1
/
2
,
J
m
=
E
y
0
,
j
+
1
/
2
m
=
E
y
I
,
j
+
1
/
2
m
=
0
,where m
=
n or n
+
1, i
=
0,1
,
…
,
I
-
1, j
=
0
,
…
,
J
-
1.
The initial values for (6)–(8) are E
x
i
¯
,
j
0
=
E
x
0
(
x
i
¯
,
y
j
), and (10)
E
y
i
,
j
¯
0
=
E
y
0
x
i
,
y
j
¯
,
H
z
i
¯
,
j
¯
0
=
H
z
0
x
i
¯
,
y
j
¯
.
In the implementation of MS-FDTDI, Stage 1 (or Stage 2) can be reduced into a tridiagonal system of linear equations forE
y
i
,
j
+
1
/
2
n
+
1 with i
=
1
,
…
,
I
-
1 (or E
x
i
+
1
/
2
,
j
n
+
1 with j
=
1
,
…
,
J
-
1) and a formula for H
z
∗ (or H
z
n
+
1), which can be solved directly.Remark 1.
(1) In order to see the difference between MS-FDTDTI and the S-FDTDII method in [15], we give the equivalent forms of the two methods:(11)
E
x
i
¯
,
j
n
+
1
-
E
x
i
¯
,
j
n
Δ
t
=
1
2
ε
δ
y
H
z
i
¯
,
j
n
+
1
+
H
z
i
¯
,
j
n
+
P
,
(12)
E
y
i
,
j
¯
n
+
1
-
E
y
i
,
j
¯
n
Δ
t
=
-
1
2
ε
δ
x
H
z
i
,
j
¯
n
+
1
+
H
z
i
,
j
¯
n
+
Δ
t
4
μ
ε
δ
x
δ
y
E
x
i
,
j
¯
n
+
1
-
E
x
i
,
j
¯
n
,
(13)
H
z
i
¯
,
j
¯
n
+
1
-
H
z
i
¯
,
j
¯
n
Δ
t
=
1
2
μ
δ
y
E
x
i
¯
,
j
¯
n
+
1
+
E
x
i
¯
,
j
¯
n
-
δ
x
E
y
i
¯
,
j
¯
n
+
1
+
E
y
i
¯
,
j
¯
n
,where (11)–(13) with P
=
0 being the equivalent form of S-FDTDII (Stage 1 of S-FDTDII is the same as (6); Stage 2 of S-FDTDII is (7)-(8) with the last term on the right hand side of (7) removed); (11)–(13) with the case (14)
P
=
Δ
t
4
μ
ε
δ
x
δ
y
E
y
i
+
1
/
2
,
j
n
+
1
-
E
y
i
+
1
/
2
,
j
nis the equivalent form of MS-FDTDI.
By these forms we see that MS-FDTDI is different from the S-FDTDII and ADI-FDTD methods (see [10, 11], where splitting of the equations is not used; however, the equivalent form of S-FDTDII is the same as that of 2D ADI-FDTD).
(2) MS-FDTDI has similar perturbation term as the D’yakonov scheme (see [19]). The equivalent form of this scheme is (11)–(13) with (15)
P
=
Δ
t
4
μ
ε
δ
x
δ
y
E
y
i
+
1
/
2
,
j
n
+
1
-
E
y
i
+
1
/
2
,
j
nand the perturbation term on the right hand side of (12) being removed. In the comparison of these equivalent forms we see that the perturbation term and its location of MS-FDTDI are different from those of the D’yakonov’s scheme. This implies that they are different.Remark 2.
Based on S-FDTDII, we propose another modified splitting FDTD method (denoted by MS-FDTDII).
Stage 1 of MS-FDTDII.(16)
E
y
i
,
j
¯
n
+
1
-
E
y
i
,
j
¯
n
Δ
t
=
-
1
2
ε
δ
x
H
z
i
,
j
¯
∗
+
H
z
i
,
j
¯
n
-
Δ
t
2
μ
ε
δ
x
δ
y
E
x
i
,
j
¯
n
;
H
z
i
¯
,
j
¯
∗
-
H
z
i
¯
,
j
¯
n
Δ
t
=
-
1
2
μ
δ
x
E
y
i
¯
,
j
¯
n
+
1
+
E
y
i
¯
,
j
¯
n
.
Stage 2 of MS-FDTDII.(17)
E
x
i
¯
,
j
n
+
1
-
E
x
i
¯
,
j
n
Δ
t
=
1
2
ε
δ
y
H
z
i
¯
,
j
n
+
1
+
H
z
i
¯
,
j
∗
+
Δ
t
2
μ
ε
δ
x
δ
y
E
y
i
¯
,
j
n
,
H
z
i
¯
,
j
¯
n
+
1
-
H
z
i
¯
,
j
¯
∗
Δ
t
=
1
2
μ
δ
y
E
x
i
¯
,
j
¯
n
+
1
+
E
x
i
¯
,
j
¯
n
.The boundary and initial conditions of MS-FDTDII are the same as MS-FDTDI.
The equivalent form of MS-FDTDII is (11)–(13) with(18)
P
=
-
Δ
t
4
μ
ε
δ
x
δ
y
E
y
i
+
1
/
2
,
j
n
+
1
-
E
y
i
+
1
/
2
,
j
n
.By these equivalent forms we see that MS-FDTDI and MS-FDTDII are of second-order accuracy.


## 2.1. Maxwell Equations

Consider the two-dimensional Maxwell equations in a lossless and homogeneous medium:(1)
∂
E
x
∂
t
=
1
ε
∂
H
z
∂
y
,
∂
E
y
∂
t
=
-
1
ε
∂
H
z
∂
x
,
∂
H
z
∂
t
=
1
μ
∂
E
x
∂
y
-
∂
E
y
∂
x
,where ε and μ are the electric permittivity and magnetic permeability of the medium and E
=
(
E
x
(
x
,
y
,
t
)
,
E
y
(
x
,
y
,
t
)
) and H
z
=
H
z
(
x
,
y
,
t
) for (
x
,
y
)
∈
Ω
=
[
0
,
a
]
×
[
0
,
b
] and t
∈
0
,
T denote the electric and magnetic fields, respectively. We assume that the spatial domain Ω is surrounded by perfectly electric conductor (PEC). Then the PEC boundary condition below is satisfied:(2)
E
,
0
×
n
→
,
0
=
0
on
0
,
T
×
∂
Ω
,where ∂
Ω denotes the boundary of Ω and n
→ is the outward normal vector on ∂
Ω. The initial conditions are assumed to be (3)
E
x
,
y
,
0
=
E
0
x
,
y
,
H
z
x
,
y
,
0
=
H
z
0
x
,
y
,where E
0
(
x
,
y
)
=
(
E
x
0
(
x
,
y
)
,
E
y
0
(
x
,
y
)
).


## 2.2. Partition of the Domains and Notations

LetΩ be partitioned as Yee’s staggered grids [1]: {
(
x
α
,
y
β
)
∣
α
=
i
,
i
+
1
/
2
,
β
=
j
,
j
+
1
/
2
}, and let [
0
,
T
] be divided into equidistant subintervals, [
t
n
,
t
n
+
1
], where(4)
x
i
=
i
Δ
x
,
x
i
+
1
/
2
=
x
i
+
1
2
Δ
x
,
x
I
=
I
Δ
x
=
a
,
y
j
=
j
Δ
y
,
y
j
+
1
/
2
=
j
+
1
2
Δ
y
,
y
J
=
J
Δ
y
=
b
,
t
n
=
n
Δ
t
,
t
n
+
1
/
2
=
t
n
+
1
2
Δ
t
,
t
N
=
N
Δ
t
=
T
,where Δ
x and Δ
y are the spatial step sizes, Δ
t is the time increment, and I, J, and N are positive integers. For a function F
(
x
,
y
,
t
) and u
,
v
=
x or y, we define(5)
F
α
,
β
m
=
F
α
Δ
x
,
β
Δ
y
,
m
Δ
t
,
δ
x
F
α
,
β
m
=
F
α
+
1
/
2
,
β
m
-
F
α
-
1
/
2
,
β
m
Δ
x
,
δ
y
F
α
,
β
m
=
F
α
,
β
+
1
/
2
m
-
F
α
,
β
-
1
/
2
m
Δ
y
,
δ
u
δ
v
F
α
,
β
m
=
δ
u
δ
v
F
α
,
β
m
.


## 2.3. Modified Splitting FDTD Methods

Denote by(
E
x
i
¯
,
j
m
,
E
y
i
,
j
¯
m
) and H
z
i
¯
,
j
¯
m the approximations to (
E
x
(
x
i
¯
,
y
j
,
t
m
), E
y
(
x
i
,
y
j
¯
,
t
m
)
) and H
z
(
x
i
¯
,
y
j
¯
,
t
m
), respectively, where, and in what follows, i
¯
=
i
+
1
/
2, j
¯
=
j
+
1
/
2. Based on the S-FDTDII scheme (see [15]) and the idea of reducing the splitting error, we propose a modified splitting FDTD method (called MS-FDTDI) for (1)–(3).Stage 1.
(6)
E
y
i
,
j
¯
n
+
1
-
E
y
i
,
j
¯
n
Δ
t
=
-
1
2
ε
δ
x
H
z
i
,
j
¯
∗
+
H
z
i
,
j
¯
n
-
Δ
t
2
μ
ε
δ
x
δ
y
E
x
i
,
j
¯
n
;
H
z
i
¯
,
j
¯
∗
-
H
z
i
¯
,
j
¯
n
Δ
t
=
-
1
2
μ
δ
x
E
y
i
¯
,
j
¯
n
+
1
+
E
y
i
¯
,
j
¯
n
.Stage 2.
(7)
E
x
i
¯
,
j
n
+
1
-
E
x
i
¯
,
j
n
Δ
t
=
1
2
ε
δ
y
H
z
i
¯
,
j
n
+
1
+
H
z
i
¯
,
j
n
+
Δ
t
4
μ
ε
δ
x
δ
y
E
y
i
¯
,
j
n
+
1
-
E
y
i
¯
,
j
n
;
(8)
H
z
i
¯
,
j
¯
n
+
1
-
H
z
i
¯
,
j
¯
∗
Δ
t
=
1
2
μ
δ
y
E
x
i
¯
,
j
¯
n
+
1
+
E
x
i
¯
,
j
¯
n
.The boundary conditions for (6)–(8) obtained from (2) are(9)
E
x
i
+
1
/
2,0
m
=
E
x
i
+
1
/
2
,
J
m
=
E
y
0
,
j
+
1
/
2
m
=
E
y
I
,
j
+
1
/
2
m
=
0
,where m
=
n or n
+
1, i
=
0,1
,
…
,
I
-
1, j
=
0
,
…
,
J
-
1.
The initial values for (6)–(8) are E
x
i
¯
,
j
0
=
E
x
0
(
x
i
¯
,
y
j
), and (10)
E
y
i
,
j
¯
0
=
E
y
0
x
i
,
y
j
¯
,
H
z
i
¯
,
j
¯
0
=
H
z
0
x
i
¯
,
y
j
¯
.
In the implementation of MS-FDTDI, Stage 1 (or Stage 2) can be reduced into a tridiagonal system of linear equations forE
y
i
,
j
+
1
/
2
n
+
1 with i
=
1
,
…
,
I
-
1 (or E
x
i
+
1
/
2
,
j
n
+
1 with j
=
1
,
…
,
J
-
1) and a formula for H
z
∗ (or H
z
n
+
1), which can be solved directly.Remark 1.
(1) In order to see the difference between MS-FDTDTI and the S-FDTDII method in [15], we give the equivalent forms of the two methods:(11)
E
x
i
¯
,
j
n
+
1
-
E
x
i
¯
,
j
n
Δ
t
=
1
2
ε
δ
y
H
z
i
¯
,
j
n
+
1
+
H
z
i
¯
,
j
n
+
P
,
(12)
E
y
i
,
j
¯
n
+
1
-
E
y
i
,
j
¯
n
Δ
t
=
-
1
2
ε
δ
x
H
z
i
,
j
¯
n
+
1
+
H
z
i
,
j
¯
n
+
Δ
t
4
μ
ε
δ
x
δ
y
E
x
i
,
j
¯
n
+
1
-
E
x
i
,
j
¯
n
,
(13)
H
z
i
¯
,
j
¯
n
+
1
-
H
z
i
¯
,
j
¯
n
Δ
t
=
1
2
μ
δ
y
E
x
i
¯
,
j
¯
n
+
1
+
E
x
i
¯
,
j
¯
n
-
δ
x
E
y
i
¯
,
j
¯
n
+
1
+
E
y
i
¯
,
j
¯
n
,where (11)–(13) with P
=
0 being the equivalent form of S-FDTDII (Stage 1 of S-FDTDII is the same as (6); Stage 2 of S-FDTDII is (7)-(8) with the last term on the right hand side of (7) removed); (11)–(13) with the case (14)
P
=
Δ
t
4
μ
ε
δ
x
δ
y
E
y
i
+
1
/
2
,
j
n
+
1
-
E
y
i
+
1
/
2
,
j
nis the equivalent form of MS-FDTDI.
By these forms we see that MS-FDTDI is different from the S-FDTDII and ADI-FDTD methods (see [10, 11], where splitting of the equations is not used; however, the equivalent form of S-FDTDII is the same as that of 2D ADI-FDTD).
(2) MS-FDTDI has similar perturbation term as the D’yakonov scheme (see [19]). The equivalent form of this scheme is (11)–(13) with (15)
P
=
Δ
t
4
μ
ε
δ
x
δ
y
E
y
i
+
1
/
2
,
j
n
+
1
-
E
y
i
+
1
/
2
,
j
nand the perturbation term on the right hand side of (12) being removed. In the comparison of these equivalent forms we see that the perturbation term and its location of MS-FDTDI are different from those of the D’yakonov’s scheme. This implies that they are different.Remark 2.
Based on S-FDTDII, we propose another modified splitting FDTD method (denoted by MS-FDTDII).
Stage 1 of MS-FDTDII.(16)
E
y
i
,
j
¯
n
+
1
-
E
y
i
,
j
¯
n
Δ
t
=
-
1
2
ε
δ
x
H
z
i
,
j
¯
∗
+
H
z
i
,
j
¯
n
-
Δ
t
2
μ
ε
δ
x
δ
y
E
x
i
,
j
¯
n
;
H
z
i
¯
,
j
¯
∗
-
H
z
i
¯
,
j
¯
n
Δ
t
=
-
1
2
μ
δ
x
E
y
i
¯
,
j
¯
n
+
1
+
E
y
i
¯
,
j
¯
n
.
Stage 2 of MS-FDTDII.(17)
E
x
i
¯
,
j
n
+
1
-
E
x
i
¯
,
j
n
Δ
t
=
1
2
ε
δ
y
H
z
i
¯
,
j
n
+
1
+
H
z
i
¯
,
j
∗
+
Δ
t
2
μ
ε
δ
x
δ
y
E
y
i
¯
,
j
n
,
H
z
i
¯
,
j
¯
n
+
1
-
H
z
i
¯
,
j
¯
∗
Δ
t
=
1
2
μ
δ
y
E
x
i
¯
,
j
¯
n
+
1
+
E
x
i
¯
,
j
¯
n
.The boundary and initial conditions of MS-FDTDII are the same as MS-FDTDI.
The equivalent form of MS-FDTDII is (11)–(13) with(18)
P
=
-
Δ
t
4
μ
ε
δ
x
δ
y
E
y
i
+
1
/
2
,
j
n
+
1
-
E
y
i
+
1
/
2
,
j
n
.By these equivalent forms we see that MS-FDTDI and MS-FDTDII are of second-order accuracy.


## 3. Analysis of Stability and Numerical Dispersion Error

In this section we first derive the amplification factors and numerical dispersion (ND) relations of MS-FDTDI and MS-FDTDII and then we analyze the stability and ND error.


### 3.1. Stability Analysis

Let the trial time-harmonic solution of the Maxwell equations be(19)
E
x
i
+
1
/
2
,
j
n
=
E
1
ξ
n
e
-
i
m
k
x
i
+
1
/
2
Δ
x
+
k
y
j
Δ
y
,
E
y
i
,
j
+
1
/
2
n
=
E
2
ξ
n
e
-
i
m
k
x
i
Δ
x
+
k
y
j
+
1
/
2
Δ
y
,
H
z
i
+
1
/
2
,
j
+
1
/
2
n
=
H
0
ξ
n
e
-
i
m
k
x
i
+
1
/
2
Δ
x
+
k
y
j
+
1
/
2
Δ
y
,where i
m
=
-
1 is the unit of complex numbers, E
1, E
2, and H
0 are the amplitudes, k
x and k
y are the wave numbers along the x-axis and y-axis, and ξ is the amplification factor.Substituting the above expressions into the equivalent form of MS-FDTDI and evaluating the determinant of the coefficient matrix of the resulting system of equations forE
1, E
2, and H
0, we get a quadratic equation of ξ. Solving this equation yields the amplification factors for MS-FDTDI:(20)
ξ
1
=
-
d
1
+
i
m
d
0
2
-
d
1
2
d
0
,
ξ
2
=
-
d
1
-
i
m
d
0
2
-
d
1
2
d
0
,where the coefficients are(21)
d
0
=
1
+
Δ
t
2
μ
ε
a
x
2
+
b
y
2
+
Δ
t
4
μ
ε
2
a
x
b
y
2
,
d
1
=
-
1
+
Δ
t
2
μ
ε
a
x
2
+
b
y
2
+
3
Δ
t
4
μ
ε
2
a
x
b
y
2
,
a
x
=
sin
⁡
0.5
k
x
Δ
x
Δ
x
,
b
y
=
sin
⁡
0.5
k
y
Δ
y
Δ
y
.The modulus ofξ
1 or ξ
2 is equal to one, implying that MS-FDTDI is unconditionally stable and nondissipative.Similarly, we obtain the amplification factors of MS-FDTDII:(22)
ξ
~
1
=
-
d
~
1
+
i
m
d
0
2
-
d
~
1
2
d
0
,
ξ
~
2
=
-
d
~
1
-
i
m
d
0
2
-
d
~
1
2
d
0
,where d
0 is the same as that in (20), and d
~
1 is (23)
d
~
1
=
-
1
+
Δ
t
2
μ
ε
a
x
2
+
b
y
2
-
Δ
t
4
μ
ε
2
a
x
b
y
2
.Thatξ
~
1
=
ξ
~
2
=
1 implies that MS-FDTDII is also unconditionally stable and nondissipative.Remark 3.
The amplification factors of S-FDTDII, which are the same as those of ADI-FDTD (the derivation is seen in [15]), are(24)
ξ
¯
1
=
-
d
¯
1
+
i
m
d
0
2
-
d
¯
1
2
d
0
,
ξ
¯
2
=
-
d
¯
1
-
i
m
d
0
2
-
d
¯
1
2
d
0
,where d
0 is the same as that in (20), and d
¯
1 is (25)
d
¯
1
=
-
1
+
Δ
t
2
μ
ε
a
x
2
+
b
y
2
+
Δ
t
4
μ
ε
2
a
x
b
y
2
.


### 3.2. Numerical Dispersion Analysis

Letc
=
1
/
μ
ε be the wave speed. Substituting ξ
=
e
i
ω
Δ
t into (20), we obtain the ND relation of MS-FDTDI:(26)
1
-
c
Δ
t
4
a
x
b
y
2
sin
2
⁡
1
2
ω
Δ
t
=
cos
2
⁡
1
2
ω
Δ
t
·
c
Δ
t
2
a
x
2
+
b
y
2
+
2
c
Δ
t
2
a
x
b
y
2
,where a
x and b
y are defined under (20).Similarly, the ND relation of MS-FDTDII is(27)
1
+
c
Δ
t
4
a
x
b
y
2
sin
2
⁡
1
2
ω
Δ
t
=
c
Δ
t
2
a
x
2
+
b
y
2
cos
2
⁡
1
2
ω
Δ
t
.Remark 4.
The ND relation of S-FDTDII is the same as that of ADI-FDTD (see [15]), which is(28)
sin
2
⁡
1
2
ω
Δ
t
=
cos
2
⁡
1
2
ω
Δ
t
·
c
Δ
t
2
a
x
2
+
b
y
2
+
c
Δ
t
2
a
x
b
y
2
.
By using the Taylor expansions ofsin
⁡
(
x
) and cos
⁡
(
x
) and the continuous dispersion relation: ω
2
=
c
2
[
(
k
x
)
2
+
(
k
y
)
2
], we derive the main truncation errors of the ND relations of MS-FDTDI, MS-FDTDII, and S-FDTDII, denoted by T
m
i, T
m
i
i, and T
s
i
i, which are(29)
T
m
i
=
-
Δ
t
4
48
ω
4
-
6
c
4
k
x
k
y
2
+
c
2
Δ
t
2
48
k
x
4
+
k
y
4
h
2
;
T
m
i
i
=
ω
4
24
Δ
t
4
+
c
2
Δ
t
2
48
k
x
4
+
k
y
4
h
2
;
T
s
i
i
=
ω
4
24
-
c
4
16
k
x
k
y
2
Δ
t
4
+
c
2
Δ
t
2
48
k
x
4
+
k
y
4
h
2
.By the second and third terms of truncation errors we see that T
s
i
i
<
T
m
i
i, implying that the ND error of S-FDTDII or ADI-FDTD is less than that of MS-FDTDII. Noting that (30)
T
m
i
-
T
s
i
i
=
-
c
4
Δ
t
4
32
k
x
2
-
k
y
2
2
+
k
x
4
+
k
y
4
<
0
,we obtain that the ND error of MS-FDTDI is less than that of S-FDTDII (or ADI-FDTD).


## 3.1. Stability Analysis

Let the trial time-harmonic solution of the Maxwell equations be(19)
E
x
i
+
1
/
2
,
j
n
=
E
1
ξ
n
e
-
i
m
k
x
i
+
1
/
2
Δ
x
+
k
y
j
Δ
y
,
E
y
i
,
j
+
1
/
2
n
=
E
2
ξ
n
e
-
i
m
k
x
i
Δ
x
+
k
y
j
+
1
/
2
Δ
y
,
H
z
i
+
1
/
2
,
j
+
1
/
2
n
=
H
0
ξ
n
e
-
i
m
k
x
i
+
1
/
2
Δ
x
+
k
y
j
+
1
/
2
Δ
y
,where i
m
=
-
1 is the unit of complex numbers, E
1, E
2, and H
0 are the amplitudes, k
x and k
y are the wave numbers along the x-axis and y-axis, and ξ is the amplification factor.Substituting the above expressions into the equivalent form of MS-FDTDI and evaluating the determinant of the coefficient matrix of the resulting system of equations forE
1, E
2, and H
0, we get a quadratic equation of ξ. Solving this equation yields the amplification factors for MS-FDTDI:(20)
ξ
1
=
-
d
1
+
i
m
d
0
2
-
d
1
2
d
0
,
ξ
2
=
-
d
1
-
i
m
d
0
2
-
d
1
2
d
0
,where the coefficients are(21)
d
0
=
1
+
Δ
t
2
μ
ε
a
x
2
+
b
y
2
+
Δ
t
4
μ
ε
2
a
x
b
y
2
,
d
1
=
-
1
+
Δ
t
2
μ
ε
a
x
2
+
b
y
2
+
3
Δ
t
4
μ
ε
2
a
x
b
y
2
,
a
x
=
sin
⁡
0.5
k
x
Δ
x
Δ
x
,
b
y
=
sin
⁡
0.5
k
y
Δ
y
Δ
y
.The modulus ofξ
1 or ξ
2 is equal to one, implying that MS-FDTDI is unconditionally stable and nondissipative.Similarly, we obtain the amplification factors of MS-FDTDII:(22)
ξ
~
1
=
-
d
~
1
+
i
m
d
0
2
-
d
~
1
2
d
0
,
ξ
~
2
=
-
d
~
1
-
i
m
d
0
2
-
d
~
1
2
d
0
,where d
0 is the same as that in (20), and d
~
1 is (23)
d
~
1
=
-
1
+
Δ
t
2
μ
ε
a
x
2
+
b
y
2
-
Δ
t
4
μ
ε
2
a
x
b
y
2
.Thatξ
~
1
=
ξ
~
2
=
1 implies that MS-FDTDII is also unconditionally stable and nondissipative.Remark 3.
The amplification factors of S-FDTDII, which are the same as those of ADI-FDTD (the derivation is seen in [15]), are(24)
ξ
¯
1
=
-
d
¯
1
+
i
m
d
0
2
-
d
¯
1
2
d
0
,
ξ
¯
2
=
-
d
¯
1
-
i
m
d
0
2
-
d
¯
1
2
d
0
,where d
0 is the same as that in (20), and d
¯
1 is (25)
d
¯
1
=
-
1
+
Δ
t
2
μ
ε
a
x
2
+
b
y
2
+
Δ
t
4
μ
ε
2
a
x
b
y
2
.


## 3.2. Numerical Dispersion Analysis

Letc
=
1
/
μ
ε be the wave speed. Substituting ξ
=
e
i
ω
Δ
t into (20), we obtain the ND relation of MS-FDTDI:(26)
1
-
c
Δ
t
4
a
x
b
y
2
sin
2
⁡
1
2
ω
Δ
t
=
cos
2
⁡
1
2
ω
Δ
t
·
c
Δ
t
2
a
x
2
+
b
y
2
+
2
c
Δ
t
2
a
x
b
y
2
,where a
x and b
y are defined under (20).Similarly, the ND relation of MS-FDTDII is(27)
1
+
c
Δ
t
4
a
x
b
y
2
sin
2
⁡
1
2
ω
Δ
t
=
c
Δ
t
2
a
x
2
+
b
y
2
cos
2
⁡
1
2
ω
Δ
t
.Remark 4.
The ND relation of S-FDTDII is the same as that of ADI-FDTD (see [15]), which is(28)
sin
2
⁡
1
2
ω
Δ
t
=
cos
2
⁡
1
2
ω
Δ
t
·
c
Δ
t
2
a
x
2
+
b
y
2
+
c
Δ
t
2
a
x
b
y
2
.
By using the Taylor expansions ofsin
⁡
(
x
) and cos
⁡
(
x
) and the continuous dispersion relation: ω
2
=
c
2
[
(
k
x
)
2
+
(
k
y
)
2
], we derive the main truncation errors of the ND relations of MS-FDTDI, MS-FDTDII, and S-FDTDII, denoted by T
m
i, T
m
i
i, and T
s
i
i, which are(29)
T
m
i
=
-
Δ
t
4
48
ω
4
-
6
c
4
k
x
k
y
2
+
c
2
Δ
t
2
48
k
x
4
+
k
y
4
h
2
;
T
m
i
i
=
ω
4
24
Δ
t
4
+
c
2
Δ
t
2
48
k
x
4
+
k
y
4
h
2
;
T
s
i
i
=
ω
4
24
-
c
4
16
k
x
k
y
2
Δ
t
4
+
c
2
Δ
t
2
48
k
x
4
+
k
y
4
h
2
.By the second and third terms of truncation errors we see that T
s
i
i
<
T
m
i
i, implying that the ND error of S-FDTDII or ADI-FDTD is less than that of MS-FDTDII. Noting that (30)
T
m
i
-
T
s
i
i
=
-
c
4
Δ
t
4
32
k
x
2
-
k
y
2
2
+
k
x
4
+
k
y
4
<
0
,we obtain that the ND error of MS-FDTDI is less than that of S-FDTDII (or ADI-FDTD).


## 4. Error Estimates and Convergence of MS-FDTDI

LetE
w
α
,
β
n
=
e
w
(
t
n
,
x
α
,
y
β
)
-
E
w
α
,
β
n and H
z
α
,
β
n
=
h
z
(
t
n
,
x
α
,
y
β
)
-
H
z
α
,
β
n, where e
w
(
t
n
,
x
α
,
y
β
) with w
=
x
,
y and h
z
(
t
n
,
x
α
,
y
β
) denote the values of the exact solution of the Maxwell equations (1)–(3) and E
w
α
,
β
n with w
=
x
,
y and H
z
α
,
β
n denote the solution of the MS-FDTDI scheme (6)–(8).Subtracting the equivalent form of MS-FDTDI (11)–(13) from the discretized Maxwell equations (whose form is like (11)–(13) with extra truncation errors), we obtain the following error equations:(31)
E
x
i
¯
,
j
n
+
1
-
E
x
i
¯
,
j
n
Δ
t
=
1
2
ε
δ
y
H
z
i
¯
,
j
n
+
1
+
H
z
i
¯
,
j
n
+
Δ
t
4
μ
ε
δ
x
δ
y
E
y
i
¯
,
j
n
+
1
-
E
y
i
¯
,
j
n
+
ξ
x
i
¯
,
j
n
+
1
/
2
,
E
y
i
,
j
¯
n
+
1
-
E
y
i
,
j
¯
n
Δ
t
=
-
1
2
ε
δ
x
H
z
i
,
j
¯
n
+
1
+
H
z
i
,
j
¯
n
+
Δ
t
4
μ
ε
δ
x
δ
y
E
x
i
,
j
¯
n
+
1
-
E
x
i
,
j
¯
n
+
ξ
y
i
,
j
¯
n
+
1
/
2
,
H
z
i
¯
,
j
¯
n
+
1
-
H
z
i
¯
,
j
¯
n
Δ
t
=
-
1
2
μ
δ
y
E
x
i
¯
,
j
¯
n
+
1
+
E
x
i
¯
,
j
¯
n
-
δ
x
E
y
i
¯
,
j
¯
n
+
1
+
E
y
i
¯
,
j
¯
n
+
η
z
i
¯
,
j
¯
n
+
1
/
2
,where ξ
w
n
+
1
/
2 with w
=
x
,
y and η
z
n
+
1
/
2 are the truncation errors, which can be derived by using Taylor formula and discretizing of Maxwell equations. These local truncation error terms are bounded by (32)
ξ
x
i
¯
,
j
n
+
1
/
2
≤
C
μ
ε
M
Δ
t
2
+
Δ
y
2
,
ξ
y
i
,
j
¯
n
+
1
/
2
≤
C
μ
ε
M
Δ
t
2
+
Δ
x
2
,
η
z
i
¯
,
j
¯
n
+
1
/
2
≤
C
μ
ε
M
Δ
t
2
+
Δ
x
2
+
Δ
y
2
,where C
μ
ε
=
1
/
24
+
1
/
8
ε
+
1
/
8
μ
+
1
/
8
μ
ε and M is a constant dependent on L
2 norms of the derivatives of the solution of (1)–(3).Multiplying both sides of (31) by ε
(
E
x
i
¯
,
j
n
+
1
+
E
x
i
¯
,
j
n
)
Δ
v, ε
(
E
y
i
,
j
¯
n
+
1
+
E
y
i
,
j
¯
n
)
Δ
v, and μ
(
H
z
i
¯
,
j
¯
n
+
1
+
H
z
i
¯
,
j
¯
n
+
1
)
Δ
v
(
Δ
v
=
Δ
t
Δ
x
Δ
y
), respectively, and applying the summation by parts and the Schwarz inequality we have(33)
α
x
E
x
N
E
x
2
+
α
y
E
y
N
E
y
2
+
1
-
Δ
t
2
H
z
N
H
z
2
≤
C
E
x
0
E
x
2
+
E
y
0
E
y
2
+
H
z
0
H
z
2
+
Δ
t
4
+
Δ
x
4
+
Δ
y
4
+
Δ
t
∑
k
=
0
N
-
1
E
x
k
E
x
2
+
E
y
k
E
y
2
+
H
z
k
H
z
2
,where for u
=
x
,
y and m
=
0
,
N
(34)
α
u
=
1
-
Δ
t
2
-
Δ
t
2
μ
ε
Δ
u
2
,
E
x
m
E
x
2
=
∑
i
=
0
I
-
1
∑
j
=
1
J
-
1
ε
E
x
m
2
Δ
x
Δ
y
,
E
y
m
E
y
2
=
∑
i
=
1
I
-
1
∑
j
=
0
J
-
1
ε
E
y
m
2
Δ
x
Δ
y
,
H
z
m
H
z
2
=
∑
i
=
0
I
-
1
∑
j
=
0
J
-
1
μ
H
z
m
2
Δ
x
Δ
y
.Moreover, if the initial conditions and step sizes satisfy(35)
E
x
0
E
x
2
+
E
y
0
E
y
2
+
H
z
0
H
z
2
≤
C
Δ
t
4
+
Δ
x
4
+
Δ
y
4
,
c
Δ
t
Δ
x
<
1
-
Δ
t
2
,
c
Δ
t
Δ
y
<
1
-
Δ
t
2
,then, by the discrete Growall’s lemma, we have(36)
E
x
N
E
x
2
+
E
y
N
E
y
2
+
H
z
N
H
z
2
≤
C
Δ
t
4
+
Δ
x
4
+
Δ
y
4
.Remark 5.
(1) The convergence of MS-FDTDI requires thatc
Δ
t
/
Δ
u
≤
(
1
-
Δ
t
/
2
) with u
=
x
,
y, which is weaker than Courant stability condition: c
Δ
t
1
/
(
Δ
x
)
2
+
1
/
(
Δ
y
)
2
<
1.
(2) By the similar method to the above it can not be proved that MS-FDTDII is convergent since the perturbation terms in this scheme are not controlled.


## 5. Numerical Experiments

We do some experiments to compute the ND errors of MS-FDTDI and MS-FDTDII, to solve a wave guide problem, and to simulate a scattering problem by the two methods.


### 5.1. Computation of Numerical Dispersion Errors

Letλ be the wave length, Δ
x
=
Δ
y
=
h, and N
λ
=
λ
/
h be the number of points per wavelength, and S
=
(
c
Δ
t
)
/
h is a multiple of the CFL number (CFL number equals 2
c
Δ
t
/
h in this case); ϕ is the wave propagation angle. Then, by k
x
=
k
cos
⁡
(
ϕ
), k
y
=
k
sin
⁡
(
ϕ
), k
2
=
k
x
2
+
k
y
2, k
=
2
π
/
λ, and the expressions of a
x and b
y(defined in Section 3.1), we see that the amplification or stability factor ξ is a function of S, ϕ, and N
λ; that is, ξ
=
ξ
(
S
,
ϕ
,
N
λ
).The ND errors of MS-FDTDI, MS-FDTDII, and S-FTTDII are computed by the following formula (see [20]): (37)
v
p
c
=
1
c
k
Δ
t
tan
-
1
⁡
I
ξ
R
ξ
=
N
λ
2
π
S
tan
-
1
⁡
I
ξ
R
ξ
,where I
(
ξ
) and R
(
ξ
) denote the imaginary and real parts of the amplification factor ξ. We plot the normalized phase velocity v
p
/
c with respect to S
,
ϕ, and N
λ (see Figures 1-2).Figure 1
Normalized phase velocities of MS-FDTDI, MS-FDTDII, and S-FDTDII against wave propagation angleϕ with S
=
3.5 and N
λ
=
40.Figure 2
Normalized phase velocities of MS-FDTDI, MS-FDTDII, and S-FDTDII against numbers of points per wavelengthN
λ with ϕ
=
3
5
∘ and S
=
1.5
.Figure1 shows the variation of v
p
/
c against the wave propagation ϕ with N
λ
=
40 and S
=
3.5 for MS-FDTDI, MS-FDTDII, and S-FDTDII. From the curves we see that v
p
/
c for MS-FDTDI is more close to 1 than that of S-FDTDII (or ADI-FDTD; S-FDTDII is equivalent to ADI-FDTD), and the latter is more close to 1 than MS-FDTDII. This means that the ND error of MS-FDTDI is less than that of S-FDTDII (or ADI-FDTD) and that the ND error of S-FDTDII (or ADI-FDTD) is less than that of MS-FDTDII.Figures2 and 3 give the graphs of v
p
/
c against N
λ with ϕ
=
3
5
∘ and S
=
1.5 and against S with ϕ
=
3
5
∘ and N
λ
=
60, respectively. From these curves we see the same conclusion as that drawn from Figure 1.Figure 3
Normalized phase velocities of MS-FDTDI, MS-FDTDII, and S-FDTDII against multiple of CFL numberS with ϕ
=
3
5
∘ and N
λ
=
60.


### 5.2. Computation of a Wave Guide Problem

Consider the normalized model problem (1)–(3) with ε
=
μ
=
1, Ω
=
[
0,1
]
×
[
0,1
], and the initial conditions derived from the exact solution of the problem, (
e
x
,
e
y
)
,
h
z, where(38)
e
x
=
-
cos
⁡
2
π
t
cos
⁡
π
x
sin
⁡
π
y
,
e
y
=
cos
⁡
2
π
t
sin
⁡
π
x
cos
⁡
π
y
,
h
z
=
-
2
sin
⁡
2
π
t
cos
⁡
π
x
cos
⁡
π
y
.In computation of this problem, we take the step sizesΔ
x
=
Δ
y
=
Δ
t
=
h
=
0.01 and use MS-FDTDI, MS-FDTDI, and S-FDTDII to find the numerical solutions, errors, and convergence orders.LetE
x
i
¯
,
j
n
=
e
x
(
x
i
¯
,
y
j
,
t
n
)
-
E
x
i
¯
,
j
n, E
y
i
,
j
¯
n
=
e
y
(
x
i
,
y
j
¯
,
t
n
)
-
E
y
i
,
j
¯
n, and H
z
i
¯
,
j
¯
n
=
h
z
(
x
i
¯
,
y
j
¯
,
t
n
)
-
H
z
i
¯
,
j
¯
n be the absolute errors, where E
x
i
¯
,
j
n, E
y
i
,
j
¯
n, and H
z
i
¯
,
j
¯
n are the approximate solutions. Then, the relative errors in the discrete L
2 norm, denoted by Re-Err-E and Re-Err-H, are written as Re-Err-E
=
Err-E
/
e
n and Re-Err-H
=
Err-H
/
h
z
n, where e
n and h
z
n are the L
2 norms of the two functions, e
n
=
(
e
x
(
x
,
y
,
t
n
)
,
e
y
(
x
,
y
,
t
n
)
) and h
z
n
=
h
z
(
x
,
y
,
t
n
), and(39)
Err-E
=
E
x
n
E
x
2
+
E
y
n
E
y
2
1
/
2
,
Err-H
=
H
z
n
H
z
2
1
/
2
.The order of convergence is calculated by the formulaOrder
=
log
2
⁡
(
E
r
r
o
r
(
h
,
Δ
t
)
/
E
r
r
o
r
(
h
/
2
,
Δ
t
/
2
)
), where l
o
g is the logarithmic function.The experimental results are shown in Tables1-2, where the drive routines are written in Fortran, and the computation was run on a 2.53 GHz PC having 2.0 GB RAM and Windows 7 operating system.Table 1
Performance of MS-FDTDI, MS-FDTDII, and S-FDTDII with different sizes of stepsΔ
x
=
Δ
y
=
Δ
t
=
h.
Scheme
h
Re-Err-E
Order
Re-Err-H
Order
MS-FDTDI
0.02
0.5286
e
-
2
0.8977
e
-
3
0.01
0.1323
e
-
2
1.9997
0.2243
e
-
3
2.0002
0.005
0.3307
e
-
3
1.9999
0.5608
e
-
4
2.0000
MS-FDTDII
0.02
0.1402
e
-
1
0.1015
e
-
2
0.01
0.3508
e
-
2
1.9991
0.2527
e
-
3
2.0058
0.005
0.8773
e
-
3
1.9997
0.6311
e
-
4
2.0015
S-FDTDII
0.02
0.9543
e
-
2
0.9556
e
-
3
0.01
0.2388
e
-
2
1.9989
0.2385
e
-
3
2.0026
0.005
0.5970
e
-
3
1.9997
0.5959
e
-
4
2.0007Table 2
Performance of MS-FDTDI, MS-FDTDII, and S-FDTDII withh
=
0.01, Δ
t
=
h, and different time length T.
Scheme
T
Re-Err-E
Re-Err-H
CPU
MS-FDTDI
10
0.1743
e
-
2
0.7760
e
-
2
0.2044
e
+
1
20
0.9051
e
-
2
0.6031
e
-
2
0.3994
e
+
1
40
0.6692
e
-
1
0.2966
e
-
2
0.8128
e
+
1
MS-FDTDII
10
0.4331
e
-
2
0.1912
e
-
1
0.2246
e
+
1
20
0.2253
e
-
1
0.1487
e
-
1
0.3838
e
+
1
40
0.1677
e
+
0
0.7322
e
-
2
0.9095
e
+
1
S-FDTDII
10
0.3041
e
-
2
0.1343
e
-
1
0.1810
e
+
1
20
0.1580
e
-
1
0.1043
e
-
1
0.3588
e
+
1
40
0.1172
e
+
0
0.5204
e
-
2
0.7301
e
+
1Table1 gives the relative errors Re-Err-E and Re-Err-H and convergence orders of the approximate electric and magnetic fields computed by MS-FDTDI, MS-FDTDII, and S-FDTDII at time t
=
1 with Δ
x
=
Δ
y
=
Δ
t
=
h
=
0.02,0.01,0.005 (the results computed by ADI-FDTD are the same as S-FDTDII). By the relative errors we see that MS-FDTDI is more accurate than S-FDTDII (or ADI-FDTD) and that S-FDTDII (or ADI-FDTD) is more accurate than MS-FDTDII. The computed convergence orders of these three methods are approximately equal to 2, implying that MS-FDTDI and MS-FDTDII are second methods. This is consistent with the analysis in theory.To see the long time behavior, Table2 lists the relative errors and CPU time (1
/
c seconds) in the cases Δ
x
=
Δ
y
=
Δ
t
=
0.01 and different time lengths T
=
10,20,40. From this table we see that MS-FDTDI is better than S-FDTDI, that S-FDTDII is better than MS-FDTDII in a long time computation, and that the CPU time for the three methods is of a little difference.


### 5.3. Computation of a Scattering Problem

The scattering problem is produced by a source at the point(
0,0.5
) in the strip S
=
{
-
1
≤
x
≤
1
,
-
∞
≤
y
≤
∞
}, where in the upper part {
-
1
≤
x
≤
1
,
0
≤
y
≤
∞
}, μ
=
1, and ε
=
1 and in the lower part of S, μ
=
1
,
ε
=
4, and the periodic boundary conditions are assumed at both sides x
=
±
1.Letr be the distance from a point to (
0,0.5
). The point source used to be the initial fields is defined by(40)
E
x
0
=
-
y
-
0.5
H
z
0
,
E
y
0
=
x
H
z
0
,
H
z
0
=
f
r
,
if
r
≤
0.3
,
0
,
otherwise
,where f
(
r
)
=
1
/
3
+
5
/
12
c
o
s
(
10
π
r
/
3
)
+
1
/
6
c
o
s
(
20
π
r
/
3
)
+
1
/
12
c
o
s
(
10
π
r
).To compute this problem we use the perfectly matched sponge layers (see [7]) to be placed in the upper (
2
≤
y
≤
2.2
) and lower (-
0.2
≤
y
≤
0) parts of the open domain with electric and magnetic losses (see [20]): (41)
σ
=
0
,
0
≤
y
≤
2
,
σ
m
y
0.2
2
,
-
0.2
≤
y
≤
0
,
σ
m
y
-
2
0.2
2
,
2
<
y
≤
2.2
,and the Maxwell equations in the sponge layers [7] are(42)
∂
E
x
∂
t
+
σ
E
x
=
1
ε
∂
H
z
∂
y
,
ε
∂
E
y
∂
t
=
-
∂
H
z
∂
x
-
σ
∫
0
t
∂
H
z
∂
x
d
t
,
∂
H
z
∂
t
+
σ
H
z
=
1
μ
∂
E
x
∂
y
-
∂
E
y
∂
x
.The implementation of MS-FDTDI for (42) in the sponge layer is given in the following two stages.Stage 1.
(43)
E
y
n
+
1
-
E
y
n
Δ
t
=
-
1
2
ε
δ
x
H
z
∗
+
H
z
n
-
σ
j
+
1
/
2
ε
∑
k
=
0
n
-
1
δ
x
H
z
k
+
H
z
k
+
1
Δ
t
2
+
δ
x
H
z
n
+
H
z
∗
Δ
t
4
-
Δ
t
2
μ
ε
1
+
Δ
t
2
σ
j
+
1
/
2
δ
x
δ
y
E
x
n
i
,
j
+
1
/
2
;
H
z
∗
-
H
z
n
Δ
t
+
σ
j
+
1
/
2
H
z
∗
+
H
z
n
2
=
-
1
2
μ
δ
x
E
y
n
+
1
+
E
y
n
-
Δ
t
2
μ
σ
j
+
1
/
2
δ
y
E
x
n
i
+
1
/
2
,
j
+
1
/
2
.Stage 2.
(44)
E
x
n
+
1
-
E
x
n
Δ
t
+
σ
j
E
x
n
+
1
+
E
x
n
2
=
1
2
ε
δ
y
H
z
n
+
1
+
H
z
n
+
Δ
t
4
μ
ε
δ
x
δ
y
E
y
n
+
1
-
E
y
n
i
+
1
/
2
,
j
,
H
z
n
+
1
-
H
z
∗
Δ
t
=
1
2
μ
δ
y
E
x
n
+
1
+
E
x
n
i
+
1
/
2
,
j
+
1
/
2
,where the subscripts of the fields for the spatial indexes are omitted for the simplicity in notation.The MS-FDTDII scheme for the Maxwell equations (42) is given in the following two stages.Stage 1.
Stage  1 is the same as Stage 1 of MS-FDTDI.Stage 2.
(45)
E
x
n
+
1
-
E
x
n
Δ
t
+
σ
j
E
x
n
+
1
+
E
x
n
2
=
1
2
ε
δ
y
H
z
n
+
1
+
H
z
n
+
Δ
t
2
μ
ε
δ
x
δ
y
E
y
n
i
+
1
/
2
,
j
,
H
z
n
+
1
-
H
z
∗
Δ
t
=
1
2
μ
δ
y
E
x
n
+
1
+
E
x
n
i
+
1
/
2
,
j
+
1
/
2
.The implementation of S-FDTDII and Yee Scheme in the sponge layer is seen in [15]. We take the step sizes Δ
x
=
Δ
y
=
Δ
t
=
0.01 and do the computation by MS-FDTDI, MS-FDTDII, S-FDTDII, and Yee Scheme. The contours of the numerical magnetic fields H
z
n are plotted in Figures 4-5.Figure 4
Contours ofH
z
n with t
n
=
0.2,0.5,1 by MS-FDTDI (a) and MS-FDTDII (b) with Δ
x
=
Δ
y
=
Δ
t
=
0.01.
(a)
(b)Figure 5
Contour ofH
z
n with t
n
=
0.2,0.5,1, by S-FDTDII (a) with Δ
x
=
Δ
y
=
Δ
t
=
0.01, and by the Yee Scheme (b) with Δ
x
=
Δ
y
=
0.01 and Δ
t
=
0.002.
(a)
(b)Figures4 and 5 give the contours of H
z
n with t
n
=
0.2,0.5 and 1 obtained by MS-FDTDI, MS-FDTDII, S-FDTDII, and Yee’s scheme. Comparing these figures, we find that the numerical solutions computed by MS-FDTDI and MS-FDTDII are in agreement with those by S-FDTDII and Yee’s scheme.Also clearly shown in the contours atT
=
0.2 and 0.5 of Figures 4-5 are the reflected and transmitted wavefronts at the dielectric interface and the absorption in the PMLs layer in the bottom of the rectangle for the four methods. The symmetry of the curves on the left and right hand sides reflects the use of periodic boundary conditions.This confirms that MS-FDTDI and MS-FDTDII are effective in solving the scattering problem.


## 5.1. Computation of Numerical Dispersion Errors

Letλ be the wave length, Δ
x
=
Δ
y
=
h, and N
λ
=
λ
/
h be the number of points per wavelength, and S
=
(
c
Δ
t
)
/
h is a multiple of the CFL number (CFL number equals 2
c
Δ
t
/
h in this case); ϕ is the wave propagation angle. Then, by k
x
=
k
cos
⁡
(
ϕ
), k
y
=
k
sin
⁡
(
ϕ
), k
2
=
k
x
2
+
k
y
2, k
=
2
π
/
λ, and the expressions of a
x and b
y(defined in Section 3.1), we see that the amplification or stability factor ξ is a function of S, ϕ, and N
λ; that is, ξ
=
ξ
(
S
,
ϕ
,
N
λ
).The ND errors of MS-FDTDI, MS-FDTDII, and S-FTTDII are computed by the following formula (see [20]): (37)
v
p
c
=
1
c
k
Δ
t
tan
-
1
⁡
I
ξ
R
ξ
=
N
λ
2
π
S
tan
-
1
⁡
I
ξ
R
ξ
,where I
(
ξ
) and R
(
ξ
) denote the imaginary and real parts of the amplification factor ξ. We plot the normalized phase velocity v
p
/
c with respect to S
,
ϕ, and N
λ (see Figures 1-2).Figure 1
Normalized phase velocities of MS-FDTDI, MS-FDTDII, and S-FDTDII against wave propagation angleϕ with S
=
3.5 and N
λ
=
40.Figure 2
Normalized phase velocities of MS-FDTDI, MS-FDTDII, and S-FDTDII against numbers of points per wavelengthN
λ with ϕ
=
3
5
∘ and S
=
1.5
.Figure1 shows the variation of v
p
/
c against the wave propagation ϕ with N
λ
=
40 and S
=
3.5 for MS-FDTDI, MS-FDTDII, and S-FDTDII. From the curves we see that v
p
/
c for MS-FDTDI is more close to 1 than that of S-FDTDII (or ADI-FDTD; S-FDTDII is equivalent to ADI-FDTD), and the latter is more close to 1 than MS-FDTDII. This means that the ND error of MS-FDTDI is less than that of S-FDTDII (or ADI-FDTD) and that the ND error of S-FDTDII (or ADI-FDTD) is less than that of MS-FDTDII.Figures2 and 3 give the graphs of v
p
/
c against N
λ with ϕ
=
3
5
∘ and S
=
1.5 and against S with ϕ
=
3
5
∘ and N
λ
=
60, respectively. From these curves we see the same conclusion as that drawn from Figure 1.Figure 3
Normalized phase velocities of MS-FDTDI, MS-FDTDII, and S-FDTDII against multiple of CFL numberS with ϕ
=
3
5
∘ and N
λ
=
60.


## 5.2. Computation of a Wave Guide Problem

Consider the normalized model problem (1)–(3) with ε
=
μ
=
1, Ω
=
[
0,1
]
×
[
0,1
], and the initial conditions derived from the exact solution of the problem, (
e
x
,
e
y
)
,
h
z, where(38)
e
x
=
-
cos
⁡
2
π
t
cos
⁡
π
x
sin
⁡
π
y
,
e
y
=
cos
⁡
2
π
t
sin
⁡
π
x
cos
⁡
π
y
,
h
z
=
-
2
sin
⁡
2
π
t
cos
⁡
π
x
cos
⁡
π
y
.In computation of this problem, we take the step sizesΔ
x
=
Δ
y
=
Δ
t
=
h
=
0.01 and use MS-FDTDI, MS-FDTDI, and S-FDTDII to find the numerical solutions, errors, and convergence orders.LetE
x
i
¯
,
j
n
=
e
x
(
x
i
¯
,
y
j
,
t
n
)
-
E
x
i
¯
,
j
n, E
y
i
,
j
¯
n
=
e
y
(
x
i
,
y
j
¯
,
t
n
)
-
E
y
i
,
j
¯
n, and H
z
i
¯
,
j
¯
n
=
h
z
(
x
i
¯
,
y
j
¯
,
t
n
)
-
H
z
i
¯
,
j
¯
n be the absolute errors, where E
x
i
¯
,
j
n, E
y
i
,
j
¯
n, and H
z
i
¯
,
j
¯
n are the approximate solutions. Then, the relative errors in the discrete L
2 norm, denoted by Re-Err-E and Re-Err-H, are written as Re-Err-E
=
Err-E
/
e
n and Re-Err-H
=
Err-H
/
h
z
n, where e
n and h
z
n are the L
2 norms of the two functions, e
n
=
(
e
x
(
x
,
y
,
t
n
)
,
e
y
(
x
,
y
,
t
n
)
) and h
z
n
=
h
z
(
x
,
y
,
t
n
), and(39)
Err-E
=
E
x
n
E
x
2
+
E
y
n
E
y
2
1
/
2
,
Err-H
=
H
z
n
H
z
2
1
/
2
.The order of convergence is calculated by the formulaOrder
=
log
2
⁡
(
E
r
r
o
r
(
h
,
Δ
t
)
/
E
r
r
o
r
(
h
/
2
,
Δ
t
/
2
)
), where l
o
g is the logarithmic function.The experimental results are shown in Tables1-2, where the drive routines are written in Fortran, and the computation was run on a 2.53 GHz PC having 2.0 GB RAM and Windows 7 operating system.Table 1
Performance of MS-FDTDI, MS-FDTDII, and S-FDTDII with different sizes of stepsΔ
x
=
Δ
y
=
Δ
t
=
h.
Scheme
h
Re-Err-E
Order
Re-Err-H
Order
MS-FDTDI
0.02
0.5286
e
-
2
0.8977
e
-
3
0.01
0.1323
e
-
2
1.9997
0.2243
e
-
3
2.0002
0.005
0.3307
e
-
3
1.9999
0.5608
e
-
4
2.0000
MS-FDTDII
0.02
0.1402
e
-
1
0.1015
e
-
2
0.01
0.3508
e
-
2
1.9991
0.2527
e
-
3
2.0058
0.005
0.8773
e
-
3
1.9997
0.6311
e
-
4
2.0015
S-FDTDII
0.02
0.9543
e
-
2
0.9556
e
-
3
0.01
0.2388
e
-
2
1.9989
0.2385
e
-
3
2.0026
0.005
0.5970
e
-
3
1.9997
0.5959
e
-
4
2.0007Table 2
Performance of MS-FDTDI, MS-FDTDII, and S-FDTDII withh
=
0.01, Δ
t
=
h, and different time length T.
Scheme
T
Re-Err-E
Re-Err-H
CPU
MS-FDTDI
10
0.1743
e
-
2
0.7760
e
-
2
0.2044
e
+
1
20
0.9051
e
-
2
0.6031
e
-
2
0.3994
e
+
1
40
0.6692
e
-
1
0.2966
e
-
2
0.8128
e
+
1
MS-FDTDII
10
0.4331
e
-
2
0.1912
e
-
1
0.2246
e
+
1
20
0.2253
e
-
1
0.1487
e
-
1
0.3838
e
+
1
40
0.1677
e
+
0
0.7322
e
-
2
0.9095
e
+
1
S-FDTDII
10
0.3041
e
-
2
0.1343
e
-
1
0.1810
e
+
1
20
0.1580
e
-
1
0.1043
e
-
1
0.3588
e
+
1
40
0.1172
e
+
0
0.5204
e
-
2
0.7301
e
+
1Table1 gives the relative errors Re-Err-E and Re-Err-H and convergence orders of the approximate electric and magnetic fields computed by MS-FDTDI, MS-FDTDII, and S-FDTDII at time t
=
1 with Δ
x
=
Δ
y
=
Δ
t
=
h
=
0.02,0.01,0.005 (the results computed by ADI-FDTD are the same as S-FDTDII). By the relative errors we see that MS-FDTDI is more accurate than S-FDTDII (or ADI-FDTD) and that S-FDTDII (or ADI-FDTD) is more accurate than MS-FDTDII. The computed convergence orders of these three methods are approximately equal to 2, implying that MS-FDTDI and MS-FDTDII are second methods. This is consistent with the analysis in theory.To see the long time behavior, Table2 lists the relative errors and CPU time (1
/
c seconds) in the cases Δ
x
=
Δ
y
=
Δ
t
=
0.01 and different time lengths T
=
10,20,40. From this table we see that MS-FDTDI is better than S-FDTDI, that S-FDTDII is better than MS-FDTDII in a long time computation, and that the CPU time for the three methods is of a little difference.


## 5.3. Computation of a Scattering Problem

The scattering problem is produced by a source at the point(
0,0.5
) in the strip S
=
{
-
1
≤
x
≤
1
,
-
∞
≤
y
≤
∞
}, where in the upper part {
-
1
≤
x
≤
1
,
0
≤
y
≤
∞
}, μ
=
1, and ε
=
1 and in the lower part of S, μ
=
1
,
ε
=
4, and the periodic boundary conditions are assumed at both sides x
=
±
1.Letr be the distance from a point to (
0,0.5
). The point source used to be the initial fields is defined by(40)
E
x
0
=
-
y
-
0.5
H
z
0
,
E
y
0
=
x
H
z
0
,
H
z
0
=
f
r
,
if
r
≤
0.3
,
0
,
otherwise
,where f
(
r
)
=
1
/
3
+
5
/
12
c
o
s
(
10
π
r
/
3
)
+
1
/
6
c
o
s
(
20
π
r
/
3
)
+
1
/
12
c
o
s
(
10
π
r
).To compute this problem we use the perfectly matched sponge layers (see [7]) to be placed in the upper (
2
≤
y
≤
2.2
) and lower (-
0.2
≤
y
≤
0) parts of the open domain with electric and magnetic losses (see [20]): (41)
σ
=
0
,
0
≤
y
≤
2
,
σ
m
y
0.2
2
,
-
0.2
≤
y
≤
0
,
σ
m
y
-
2
0.2
2
,
2
<
y
≤
2.2
,and the Maxwell equations in the sponge layers [7] are(42)
∂
E
x
∂
t
+
σ
E
x
=
1
ε
∂
H
z
∂
y
,
ε
∂
E
y
∂
t
=
-
∂
H
z
∂
x
-
σ
∫
0
t
∂
H
z
∂
x
d
t
,
∂
H
z
∂
t
+
σ
H
z
=
1
μ
∂
E
x
∂
y
-
∂
E
y
∂
x
.The implementation of MS-FDTDI for (42) in the sponge layer is given in the following two stages.Stage 1.
(43)
E
y
n
+
1
-
E
y
n
Δ
t
=
-
1
2
ε
δ
x
H
z
∗
+
H
z
n
-
σ
j
+
1
/
2
ε
∑
k
=
0
n
-
1
δ
x
H
z
k
+
H
z
k
+
1
Δ
t
2
+
δ
x
H
z
n
+
H
z
∗
Δ
t
4
-
Δ
t
2
μ
ε
1
+
Δ
t
2
σ
j
+
1
/
2
δ
x
δ
y
E
x
n
i
,
j
+
1
/
2
;
H
z
∗
-
H
z
n
Δ
t
+
σ
j
+
1
/
2
H
z
∗
+
H
z
n
2
=
-
1
2
μ
δ
x
E
y
n
+
1
+
E
y
n
-
Δ
t
2
μ
σ
j
+
1
/
2
δ
y
E
x
n
i
+
1
/
2
,
j
+
1
/
2
.Stage 2.
(44)
E
x
n
+
1
-
E
x
n
Δ
t
+
σ
j
E
x
n
+
1
+
E
x
n
2
=
1
2
ε
δ
y
H
z
n
+
1
+
H
z
n
+
Δ
t
4
μ
ε
δ
x
δ
y
E
y
n
+
1
-
E
y
n
i
+
1
/
2
,
j
,
H
z
n
+
1
-
H
z
∗
Δ
t
=
1
2
μ
δ
y
E
x
n
+
1
+
E
x
n
i
+
1
/
2
,
j
+
1
/
2
,where the subscripts of the fields for the spatial indexes are omitted for the simplicity in notation.The MS-FDTDII scheme for the Maxwell equations (42) is given in the following two stages.Stage 1.
Stage  1 is the same as Stage 1 of MS-FDTDI.Stage 2.
(45)
E
x
n
+
1
-
E
x
n
Δ
t
+
σ
j
E
x
n
+
1
+
E
x
n
2
=
1
2
ε
δ
y
H
z
n
+
1
+
H
z
n
+
Δ
t
2
μ
ε
δ
x
δ
y
E
y
n
i
+
1
/
2
,
j
,
H
z
n
+
1
-
H
z
∗
Δ
t
=
1
2
μ
δ
y
E
x
n
+
1
+
E
x
n
i
+
1
/
2
,
j
+
1
/
2
.The implementation of S-FDTDII and Yee Scheme in the sponge layer is seen in [15]. We take the step sizes Δ
x
=
Δ
y
=
Δ
t
=
0.01 and do the computation by MS-FDTDI, MS-FDTDII, S-FDTDII, and Yee Scheme. The contours of the numerical magnetic fields H
z
n are plotted in Figures 4-5.Figure 4
Contours ofH
z
n with t
n
=
0.2,0.5,1 by MS-FDTDI (a) and MS-FDTDII (b) with Δ
x
=
Δ
y
=
Δ
t
=
0.01.
(a)
(b)Figure 5
Contour ofH
z
n with t
n
=
0.2,0.5,1, by S-FDTDII (a) with Δ
x
=
Δ
y
=
Δ
t
=
0.01, and by the Yee Scheme (b) with Δ
x
=
Δ
y
=
0.01 and Δ
t
=
0.002.
(a)
(b)Figures4 and 5 give the contours of H
z
n with t
n
=
0.2,0.5 and 1 obtained by MS-FDTDI, MS-FDTDII, S-FDTDII, and Yee’s scheme. Comparing these figures, we find that the numerical solutions computed by MS-FDTDI and MS-FDTDII are in agreement with those by S-FDTDII and Yee’s scheme.Also clearly shown in the contours atT
=
0.2 and 0.5 of Figures 4-5 are the reflected and transmitted wavefronts at the dielectric interface and the absorption in the PMLs layer in the bottom of the rectangle for the four methods. The symmetry of the curves on the left and right hand sides reflects the use of periodic boundary conditions.This confirms that MS-FDTDI and MS-FDTDII are effective in solving the scattering problem.


## 6. Conclusions and Remarks

In this letter we proposed two FDTD methods (MS-FDTDI and MS-FDTDII) for the 2D Maxwell’s equations by introducing two new methods to reduce the perturbation error caused by splitting of Maxwell equations. It was shown that the two methods are second-order accurate. By energy method MS-FDTDI was proved to be second-order convergent. By Fourier method the unconditional stability of the two methods was proved and the numerical dispersion (ND) relations were derived. By analyzing the truncation errors of the ND relations, we proved rigorously that MS-FDTDI has less ND error than S-FDTDII or ADI-FDTD and that the ND error of S-FDTDII or ADI-FDTD is less than that of MS-FDTDII. Numerical experiments were carried out, and the analysis on stability and ND error as well as the efficiency of MS-FDTDI and MS-FDTDII in computing a wave guide problem and a scattering problem was confirmed. The new methods for reducing splitting errors and analyzing ND errors could be used for construction of other methods and ND analysis.



---

*Source: 6063176-2017-02-22.xml*