要彻底弄懂FDTD,首先就要从基本的麦克斯韦方程组出发,一步步推导至我们编程需要的公式。
典型的时域麦克斯韦方程组如下:
{
∂
H
z
∂
y
−
∂
H
y
∂
z
=
ϵ
∂
E
x
∂
t
+
σ
E
x
∂
H
x
∂
z
−
∂
H
z
∂
x
=
ϵ
∂
E
y
∂
t
+
σ
E
y
∂
H
y
∂
x
−
∂
H
x
∂
y
=
ϵ
∂
E
z
∂
t
+
σ
E
z
\begin{cases} \frac {\partial H_z}{\partial y}-\frac {\partial H_y}{\partial z}=\epsilon\frac {\partial E_x}{\partial t}+\sigma E_x\\[2ex] \frac {\partial H_x}{\partial z}-\frac {\partial H_z}{\partial x}=\epsilon\frac {\partial E_y}{\partial t}+\sigma E_y\\[2ex] \frac {\partial H_y}{\partial x}-\frac {\partial H_x}{\partial y}=\epsilon\frac {\partial E_z}{\partial t}+\sigma E_z \end{cases}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧∂y∂Hz−∂z∂Hy=ϵ∂t∂Ex+σEx∂z∂Hx−∂x∂Hz=ϵ∂t∂Ey+σEy∂x∂Hy−∂y∂Hx=ϵ∂t∂Ez+σEz
以及
{
∂
E
z
∂
y
−
∂
E
y
∂
z
=
−
μ
∂
H
x
∂
t
−
σ
m
H
x
∂
E
x
∂
z
−
∂
E
z
∂
x
=
−
μ
∂
H
y
∂
t
−
σ
m
H
y
∂
E
y
∂
x
−
∂
E
x
∂
y
=
−
μ
∂
H
z
∂
t
−
σ
m
H
z
\begin{cases} \frac {\partial E_z}{\partial y}-\frac {\partial E_y}{\partial z}=-\mu\frac {\partial H_x}{\partial t}-\sigma_m H_x\\[2ex] \frac {\partial E_x}{\partial z}-\frac {\partial E_z}{\partial x}=-\mu\frac {\partial H_y}{\partial t}-\sigma_m H_y\\[2ex] \frac {\partial E_y}{\partial x}-\frac {\partial E_x}{\partial y}=-\mu\frac {\partial H_z}{\partial t}-\sigma_m H_z \end{cases}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧∂y∂Ez−∂z∂Ey=−μ∂t∂Hx−σmHx∂z∂Ex−∂x∂Ez=−μ∂t∂Hy−σmHy∂x∂Ey−∂y∂Ex=−μ∂t∂Hz−σmHz
为了后续推导的方便,我们首先定义
f
(
x
,
y
,
z
,
t
)
=
f
(
i
Δ
x
,
j
Δ
y
,
k
Δ
z
,
n
Δ
t
)
=
f
n
(
i
,
j
,
k
)
f(x, y, z, t) = f(i\Delta x, j\Delta y, k\Delta z, n\Delta t)=f^n(i, j, k)
f(x,y,z,t)=f(iΔx,jΔy,kΔz,nΔt)=fn(i,j,k)
然后对麦克斯韦方程组中的一节偏导数取中心差分近似,有
{
∂
f
(
x
,
y
,
z
,
t
)
∂
x
∣
x
=
i
Δ
x
≈
f
n
(
i
+
1
/
2
,
j
,
k
)
−
f
n
(
i
−
1
/
2
,
j
,
k
)
Δ
x
∂
f
(
x
,
y
,
z
,
t
)
∂
y
∣
x
=
j
Δ
y
≈
f
n
(
i
,
j
+
1
/
2
,
k
)
−
f
n
(
i
,
j
−
1
/
2
,
k
)
Δ
y
∂
f
(
x
,
y
,
z
,
t
)
∂
z
∣
x
=
k
Δ
z
≈
f
n
(
i
,
j
,
k
+
1
/
2
)
−
f
n
(
i
,
j
,
k
−
1
/
2
)
Δ
z
∂
f
(
x
,
y
,
z
,
t
)
∂
t
∣
x
=
n
Δ
t
≈
f
n
+
1
/
2
(
i
,
j
,
k
)
−
f
n
−
1
/
2
(
i
,
j
,
k
)
Δ
t
\begin{cases} \frac {\partial f(x, y, z, t)}{\partial x}|_{x=i\Delta x}\approx \frac { f^n(i+1/2, j, k) - f^n(i-1/2, j, k)}{\Delta x}\\[2ex] \frac {\partial f(x, y, z, t)}{\partial y}|_{x=j\Delta y}\approx \frac { f^n(i, j + 1/2, k) - f^n(i, j - 1/2, k)}{\Delta y}\\[2ex] \frac {\partial f(x, y, z, t)}{\partial z}|_{x=k\Delta z}\approx \frac { f^n(i, j, k + 1/2) - f^n(i, j, k - 1/2)}{\Delta z}\\[2ex] \frac {\partial f(x, y, z, t)}{\partial t}|_{x=n\Delta t}\approx \frac { f^{n + 1/2}(i, j, k) - f^{n - 1/2}(i, j, k)}{\Delta t} \end{cases}
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧∂x∂f(x,y,z,t)∣x=iΔx≈Δxfn(i+1/2,j,k)−fn(i−1/2,j,k)∂y∂f(x,y,z,t)∣x=jΔy≈Δyfn(i,j+1/2,k)−fn(i,j−1/2,k)∂z∂f(x,y,z,t)∣x=kΔz≈Δzfn(i,j,k+1/2)−fn(i,j,k−1/2)∂t∂f(x,y,z,t)∣x=nΔt≈Δtfn+1/2(i,j,k)−fn−1/2(i,j,k)
然后便可以得到著名的Yee元胞:
由上述推导以及Yee元胞中可以看到,在FDTD中电场与磁场的采样位置和采样时间皆相差半个步长。如此,当给定电磁问题的初始值以及边界条件后,便可以根据差分公式一步步地推导空间任意点处的电场与磁场。
下面便是根据差分公式推导出的麦克斯韦方程组的FDTD形式(无耗各向同性媒质中):
{
E
x
n
+
1
(
i
+
1
/
2
,
j
,
k
)
=
C
A
⋅
E
x
n
(
i
+
1
/
2
,
j
,
k
)
+
C
B
⋅
[
H
z
n
+
1
/
2
(
i
+
1
/
2
,
j
+
1
/
2
,
k
)
−
H
z
n
+
1
/
2
(
i
+
1
/
2
,
j
−
1
/
2
,
k
)
Δ
y
−
H
y
n
+
1
/
2
(
i
+
1
/
2
,
j
,
k
+
1
/
2
)
−
H
y
n
+
1
/
2
(
i
+
1
/
2
,
j
,
k
−
1
/
2
)
Δ
z
]
E
y
n
+
1
(
i
,
j
+
1
/
2
,
k
)
=
C
A
⋅
E
y
n
(
i
,
j
+
1
/
2
,
k
)
+
C
B
⋅
[
H
x
n
+
1
/
2
(
i
,
j
+
1
/
2
,
k
+
1
/
2
)
−
H
x
n
+
1
/
2
(
i
,
j
+
1
/
2
,
k
−
1
/
2
)
Δ
z
−
H
z
n
+
1
/
2
(
i
+
1
/
2
,
j
+
1
/
2
,
k
)
−
H
z
n
+
1
/
2
(
i
−
1
/
2
,
j
+
1
/
2
,
k
)
Δ
x
]
E
z
n
+
1
(
i
,
j
,
k
+
1
/
2
)
=
C
A
⋅
E
z
n
(
i
,
j
,
k
+
1
/
2
)
+
C
B
⋅
[
H
y
n
+
1
/
2
(
i
+
1
/
2
,
j
,
k
+
1
/
2
)
−
H
y
n
+
1
/
2
(
i
−
1
/
2
,
j
,
k
+
1
/
2
)
Δ
x
−
H
x
n
+
1
/
2
(
i
,
j
+
1
/
2
,
k
+
1
/
2
)
−
H
x
n
+
1
/
2
(
i
,
j
−
1
/
2
,
k
+
1
/
2
)
Δ
y
]
\begin{cases} {E^{n + 1}_{x}(i + 1/2, j, k)} = CA \cdot {E^{n}_{x}(i + 1/2, j, k)} + CB \cdot \left [ \frac { H^{n + 1/2}_z(i + 1/2, j + 1/2, k) - H^{n + 1/2}_z(i + 1/2, j - 1/2, k)}{\Delta y} - \frac { H^{n + 1/2}_y(i + 1/2, j, k + 1/2) - H^{n + 1/2}_y(i + 1/2, j, k - 1/2)}{\Delta z} \right ] \\[2ex] {E^{n + 1}_{y}(i, j + 1/2, k)} = CA \cdot {E^{n}_{y}(i, j + 1/2, k)} + CB \cdot \left [ \frac { H^{n + 1/2}_x(i, j + 1/2, k + 1/2) - H^{n + 1/2}_x(i, j + 1/2, k - 1/2)}{\Delta z} - \frac { H^{n + 1/2}_z(i + 1/2, j + 1/2, k) - H^{n + 1/2}_z(i - 1/2, j + 1/2, k)}{\Delta x} \right ] \\[2ex] {E^{n + 1}_{z}(i, j, k + 1/2)} = CA \cdot {E^{n}_{z}(i, j, k + 1/2)} + CB \cdot \left [ \frac { H^{n + 1/2}_y(i + 1/2, j, k + 1/2) - H^{n + 1/2}_y(i - 1/2, j, k + 1/2)}{\Delta x} - \frac { H^{n + 1/2}_x(i, j + 1/2, k + 1/2) - H^{n + 1/2}_x(i, j - 1/2, k + 1/2)}{\Delta y} \right ] \end{cases}
⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧Exn+1(i+1/2,j,k)=CA⋅Exn(i+1/2,j,k)+CB⋅[ΔyHzn+1/2(i+1/2,j+1/2,k)−Hzn+1/2(i+1/2,j−1/2,k)−ΔzHyn+1/2(i+1/2,j,k+1/2)−Hyn+1/2(i+1/2,j,k−1/2)]Eyn+1(i,j+1/2,k)=CA⋅Eyn(i,j+1/2,k)+CB⋅[ΔzHxn+1/2(i,j+1/2,k+1/2)−Hxn+1/2(i,j+1/2,k−1/2)−ΔxHzn+1/2(i+1/2,j+1/2,k)−Hzn+1/2(i−1/2,j+1/2,k)]Ezn+1(i,j,k+1/2)=CA⋅Ezn(i,j,k+1/2)+CB⋅[ΔxHyn+1/2(i+1/2,j,k+1/2)−Hyn+1/2(i−1/2,j,k+1/2)−ΔyHxn+1/2(i,j+1/2,k+1/2)−Hxn+1/2(i,j−1/2,k+1/2)]
以及
{
H
x
n
+
1
/
2
(
i
,
j
+
1
/
2
,
k
+
1
/
2
)
=
C
P
⋅
H
x
n
−
1
/
2
(
i
,
j
+
1
/
2
,
k
+
1
/
2
)
−
C
Q
⋅
[
E
z
n
(
i
,
j
+
1
,
k
+
1
/
2
)
−
E
z
n
(
i
,
j
,
k
+
1
/
2
)
Δ
y
−
E
y
n
(
i
,
j
+
1
/
2
,
k
+
1
)
−
E
y
n
(
i
,
j
+
1
/
2
,
k
)
Δ
z
]
H
y
n
+
1
/
2
(
i
+
1
/
2
,
j
,
k
+
1
/
2
)
=
C
P
⋅
H
y
n
−
1
/
2
(
i
+
1
/
2
,
j
,
k
+
1
/
2
)
−
C
Q
⋅
[
E
x
n
(
i
+
1
/
2
,
j
,
k
+
1
)
−
E
x
n
(
i
+
1
/
2
,
j
,
k
)
Δ
z
−
E
z
n
(
i
+
1
,
j
,
k
+
1
/
2
)
−
E
z
n
(
i
,
j
,
k
+
1
/
2
)
Δ
x
]
H
z
n
+
1
/
2
(
i
+
1
/
2
,
j
+
1
/
2
,
k
)
=
C
P
⋅
H
z
n
−
1
/
2
(
i
+
1
/
2
,
j
+
1
/
2
,
k
)
−
C
Q
⋅
[
E
y
n
(
i
+
1
,
j
+
1
/
2
,
k
)
−
E
y
n
(
i
,
j
+
1
/
2
,
k
)
Δ
x
−
E
x
n
(
i
+
1
/
2
,
j
+
1
,
k
)
−
E
x
n
(
i
+
1
/
2
,
j
,
k
)
Δ
y
]
\begin{cases} {H^{n + 1/2}_{x}(i, j + 1/2, k + 1/2)} = CP \cdot {H^{n - 1/2}_{x}(i, j + 1/2, k + 1/2)} - CQ \cdot \left [ \frac { E^{n}_z(i, j + 1, k + 1/2) - E^{n}_z(i, j, k + 1/2)}{\Delta y} - \frac { E^{n}_y(i, j + 1/2, k + 1) - E^{n}_y(i, j + 1/2, k)}{\Delta z} \right ] \\[2ex] {H^{n + 1/2}_{y}(i + 1/2, j, k + 1/2)} = CP \cdot {H^{n - 1/2}_{y}(i + 1/2, j, k + 1/2)} - CQ \cdot \left [ \frac { E^{n}_x(i + 1/2, j, k + 1) - E^{n}_x(i + 1/2, j, k)}{\Delta z} - \frac { E^{n}_z(i + 1, j, k + 1/2) - E^{n}_z(i, j, k + 1/2)}{\Delta x} \right ] \\[2ex] {H^{n + 1/2}_{z}(i + 1/2, j + 1/2, k)} = CP \cdot {H^{n - 1/2}_{z}(i + 1/2, j + 1/2, k)} - CQ \cdot \left [ \frac { E^{n}_y(i + 1, j + 1/2, k) - E^{n}_y(i, j + 1/2, k)}{\Delta x} - \frac { E^{n}_x(i + 1/2, j + 1, k) - E^{n}_x(i + 1/2, j, k)}{\Delta y} \right ] \end{cases}
⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧Hxn+1/2(i,j+1/2,k+1/2)=CP⋅Hxn−1/2(i,j+1/2,k+1/2)−CQ⋅[ΔyEzn(i,j+1,k+1/2)−Ezn(i,j,k+1/2)−ΔzEyn(i,j+1/2,k+1)−Eyn(i,j+1/2,k)]Hyn+1/2(i+1/2,j,k+1/2)=CP⋅Hyn−1/2(i+1/2,j,k+1/2)−CQ⋅[ΔzExn(i+1/2,j,k+1)−Exn(i+1/2,j,k)−ΔxEzn(i+1,j,k+1/2)−Ezn(i,j,k+1/2)]Hzn+1/2(i+1/2,j+1/2,k)=CP⋅Hzn−1/2(i+1/2,j+1/2,k)−CQ⋅[ΔxEyn(i+1,j+1/2,k)−Eyn(i,j+1/2,k)−ΔyExn(i+1/2,j+1,k)−Exn(i+1/2,j,k)]
上式中
C
A
=
ϵ
Δ
t
−
σ
2
ϵ
Δ
t
+
σ
2
=
1
−
σ
Δ
t
2
ϵ
1
+
σ
Δ
t
2
ϵ
CA = \frac {\frac {\epsilon}{\Delta t} - \frac {\sigma}{2}}{\frac {\epsilon}{\Delta t} + \frac {\sigma}{2}} = \frac {1 - \frac {\sigma \Delta t}{2 \epsilon}}{1 + \frac {\sigma \Delta t}{2 \epsilon}}
CA=Δtϵ+2σΔtϵ−2σ=1+2ϵσΔt1−2ϵσΔt
C
B
=
1
ϵ
Δ
t
+
σ
2
=
Δ
t
ϵ
1
+
σ
Δ
t
2
ϵ
CB = \frac {1}{\frac {\epsilon}{\Delta t} + \frac {\sigma}{2}} = \frac {\frac {\Delta t}{\epsilon}}{1 + \frac {\sigma \Delta t}{2 \epsilon}}
CB=Δtϵ+2σ1=1+2ϵσΔtϵΔt
C
P
=
μ
Δ
t
−
σ
m
2
μ
Δ
t
+
σ
m
2
=
1
−
σ
m
Δ
t
2
μ
1
+
σ
m
Δ
t
2
μ
CP = \frac {\frac {\mu}{\Delta t} - \frac {\sigma_m}{2}}{\frac {\mu}{\Delta t} + \frac {\sigma_m}{2}} = \frac {1 - \frac {\sigma_m \Delta t}{2 \mu}}{1 + \frac {\sigma_m \Delta t}{2 \mu}}
CP=Δtμ+2σmΔtμ−2σm=1+2μσmΔt1−2μσmΔt
C
Q
=
1
μ
Δ
t
+
σ
m
2
=
Δ
t
μ
1
+
σ
m
Δ
t
2
μ
CQ = \frac {1}{\frac {\mu}{\Delta t} + \frac {\sigma_m}{2}} = \frac {\frac {\Delta t}{\mu}}{1 + \frac {\sigma_m \Delta t}{2 \mu}}
CQ=Δtμ+2σm1=1+2μσmΔtμΔt
接下来讲解具体的编程思路,在上述FDTD公式中包含很多“半步长”的式子,而这些“半步长”的公式在C++中是比较难直接实现的。所以为了简化编程的复杂度,我们把这些“半步长”阉割掉,同时六个节点都用相同的整数步长。牢记,对任意给定的一组整数步长,电场节点向其所指的方向移动半个步长,磁场节点向其未指的方向移动半个步长,例如FDTD的电场公式:
{
E
x
n
+
1
(
i
+
1
/
2
,
j
,
k
)
=
C
A
⋅
E
x
n
(
i
+
1
/
2
,
j
,
k
)
+
C
B
⋅
[
H
z
n
+
1
/
2
(
i
+
1
/
2
,
j
+
1
/
2
,
k
)
−
H
z
n
+
1
/
2
(
i
+
1
/
2
,
j
−
1
/
2
,
k
)
Δ
y
−
H
y
n
+
1
/
2
(
i
+
1
/
2
,
j
,
k
+
1
/
2
)
−
H
y
n
+
1
/
2
(
i
+
1
/
2
,
j
,
k
−
1
/
2
)
Δ
z
]
E
y
n
+
1
(
i
,
j
+
1
/
2
,
k
)
=
C
A
⋅
E
y
n
(
i
,
j
+
1
/
2
,
k
)
+
C
B
⋅
[
H
x
n
+
1
/
2
(
i
,
j
+
1
/
2
,
k
+
1
/
2
)
−
H
x
n
+
1
/
2
(
i
,
j
+
1
/
2
,
k
−
1
/
2
)
Δ
z
−
H
z
n
+
1
/
2
(
i
+
1
/
2
,
j
+
1
/
2
,
k
)
−
H
z
n
+
1
/
2
(
i
−
1
/
2
,
j
+
1
/
2
,
k
)
Δ
x
]
E
z
n
+
1
(
i
,
j
,
k
+
1
/
2
)
=
C
A
⋅
E
z
n
(
i
,
j
,
k
+
1
/
2
)
+
C
B
⋅
[
H
y
n
+
1
/
2
(
i
+
1
/
2
,
j
,
k
+
1
/
2
)
−
H
y
n
+
1
/
2
(
i
−
1
/
2
,
j
,
k
+
1
/
2
)
Δ
x
−
H
x
n
+
1
/
2
(
i
,
j
+
1
/
2
,
k
+
1
/
2
)
−
H
x
n
+
1
/
2
(
i
,
j
−
1
/
2
,
k
+
1
/
2
)
Δ
y
]
\begin{cases} {E^{n + 1}_{x}(i + 1/2, j, k)} = CA \cdot {E^{n}_{x}(i + 1/2, j, k)} + CB \cdot \left [ \frac { H^{n + 1/2}_z(i + 1/2, j + 1/2, k) - H^{n + 1/2}_z(i + 1/2, j - 1/2, k)}{\Delta y} - \frac { H^{n + 1/2}_y(i + 1/2, j, k + 1/2) - H^{n + 1/2}_y(i + 1/2, j, k - 1/2)}{\Delta z} \right ] \\[2ex] {E^{n + 1}_{y}(i, j + 1/2, k)} = CA \cdot {E^{n}_{y}(i, j + 1/2, k)} + CB \cdot \left [ \frac { H^{n + 1/2}_x(i, j + 1/2, k + 1/2) - H^{n + 1/2}_x(i, j + 1/2, k - 1/2)}{\Delta z} - \frac { H^{n + 1/2}_z(i + 1/2, j + 1/2, k) - H^{n + 1/2}_z(i - 1/2, j + 1/2, k)}{\Delta x} \right ] \\[2ex] {E^{n + 1}_{z}(i, j, k + 1/2)} = CA \cdot {E^{n}_{z}(i, j, k + 1/2)} + CB \cdot \left [ \frac { H^{n + 1/2}_y(i + 1/2, j, k + 1/2) - H^{n + 1/2}_y(i - 1/2, j, k + 1/2)}{\Delta x} - \frac { H^{n + 1/2}_x(i, j + 1/2, k + 1/2) - H^{n + 1/2}_x(i, j - 1/2, k + 1/2)}{\Delta y} \right ] \end{cases}
⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧Exn+1(i+1/2,j,k)=CA⋅Exn(i+1/2,j,k)+CB⋅[ΔyHzn+1/2(i+1/2,j+1/2,k)−Hzn+1/2(i+1/2,j−1/2,k)−ΔzHyn+1/2(i+1/2,j,k+1/2)−Hyn+1/2(i+1/2,j,k−1/2)]Eyn+1(i,j+1/2,k)=CA⋅Eyn(i,j+1/2,k)+CB⋅[ΔzHxn+1/2(i,j+1/2,k+1/2)−Hxn+1/2(i,j+1/2,k−1/2)−ΔxHzn+1/2(i+1/2,j+1/2,k)−Hzn+1/2(i−1/2,j+1/2,k)]Ezn+1(i,j,k+1/2)=CA⋅Ezn(i,j,k+1/2)+CB⋅[ΔxHyn+1/2(i+1/2,j,k+1/2)−Hyn+1/2(i−1/2,j,k+1/2)−ΔyHxn+1/2(i,j+1/2,k+1/2)−Hxn+1/2(i,j−1/2,k+1/2)]
在实际编程时实现时,我们会把它简化成:
{
E
x
n
+
1
(
i
,
j
,
k
)
=
C
A
⋅
E
x
n
(
i
,
j
,
k
)
+
C
B
⋅
[
H
z
n
(
i
,
j
,
k
)
−
H
z
n
(
i
,
j
−
1
,
k
)
Δ
y
−
H
y
n
(
i
,
j
,
k
)
−
H
y
n
(
i
,
j
,
k
−
1
)
Δ
z
]
E
y
n
+
1
(
i
,
j
,
k
)
=
C
A
⋅
E
y
n
(
i
,
j
,
k
)
+
C
B
⋅
[
H
x
n
(
i
,
j
,
k
)
−
H
x
n
(
i
,
j
,
k
−
1
)
Δ
z
−
H
z
n
(
i
,
j
,
k
)
−
H
z
n
(
i
−
1
,
j
,
k
)
Δ
x
]
E
z
n
+
1
(
i
,
j
,
k
)
=
C
A
⋅
E
z
n
(
i
,
j
,
k
)
+
C
B
⋅
[
H
y
n
(
i
,
j
,
k
)
−
H
y
n
(
i
−
1
,
j
,
k
)
Δ
x
−
H
x
n
(
i
,
j
,
k
)
−
H
x
n
(
i
,
j
−
1
,
k
)
Δ
y
]
\begin{cases} {E^{n + 1}_{x}(i, j, k)} = CA \cdot {E^{n}_{x}(i, j, k)} + CB \cdot \left [ \frac { H^{n}_z(i, j, k) - H^{n}_z(i, j - 1, k)}{\Delta y} - \frac { H^{n}_y(i, j, k) - H^{n}_y(i, j, k - 1)}{\Delta z} \right ] \\[2ex] {E^{n + 1}_{y}(i, j, k)} = CA \cdot {E^{n}_{y}(i, j, k)} + CB \cdot \left [ \frac { H^{n}_x(i, j, k) - H^{n}_x(i, j, k - 1)}{\Delta z} - \frac { H^{n}_z(i, j, k) - H^{n}_z(i - 1, j, k)}{\Delta x} \right ] \\[2ex] {E^{n + 1}_{z}(i, j, k)} = CA \cdot {E^{n}_{z}(i, j, k)} + CB \cdot \left [ \frac { H^{n}_y(i, j, k) - H^{n}_y(i - 1, j, k)}{\Delta x} - \frac { H^{n}_x(i, j, k) - H^{n}_x(i, j - 1, k)}{\Delta y} \right ] \end{cases}
⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧Exn+1(i,j,k)=CA⋅Exn(i,j,k)+CB⋅[ΔyHzn(i,j,k)−Hzn(i,j−1,k)−ΔzHyn(i,j,k)−Hyn(i,j,k−1)]Eyn+1(i,j,k)=CA⋅Eyn(i,j,k)+CB⋅[ΔzHxn(i,j,k)−Hxn(i,j,k−1)−ΔxHzn(i,j,k)−Hzn(i−1,j,k)]Ezn+1(i,j,k)=CA⋅Ezn(i,j,k)+CB⋅[ΔxHyn(i,j,k)−Hyn(i−1,j,k)−ΔyHxn(i,j,k)−Hxn(i,j−1,k)]
磁场的FDTD则简化为:
{
H
x
n
+
1
(
i
,
j
,
k
)
=
C
P
⋅
H
x
n
(
i
,
j
,
k
)
−
C
Q
⋅
[
E
z
n
+
1
(
i
,
j
+
1
,
k
)
−
E
z
n
+
1
(
i
,
j
,
k
)
Δ
y
−
E
y
n
+
1
(
i
,
j
,
k
+
1
)
−
E
y
n
+
1
(
i
,
j
,
k
)
Δ
z
]
H
y
n
+
1
(
i
,
j
,
k
)
=
C
P
⋅
H
y
n
(
i
,
j
,
k
)
−
C
Q
⋅
[
E
x
n
+
1
(
i
,
j
,
k
+
1
)
−
E
x
n
+
1
(
i
,
j
,
k
)
Δ
z
−
E
z
n
+
1
(
i
+
1
,
j
,
k
)
−
E
z
n
+
1
(
i
,
j
,
k
)
Δ
x
]
H
z
n
+
1
(
i
,
j
,
k
)
=
C
P
⋅
H
z
n
(
i
,
j
,
k
)
−
C
Q
⋅
[
E
y
n
+
1
(
i
+
1
,
j
,
k
)
−
E
y
n
+
1
(
i
,
j
,
k
)
Δ
x
−
E
x
n
+
1
(
i
,
j
+
1
,
k
)
−
E
x
n
+
1
(
i
,
j
,
k
)
Δ
y
]
\begin{cases} {H^{n + 1}_{x}(i, j, k)} = CP \cdot {H^{n}_{x}(i, j, k)} - CQ \cdot \left [ \frac { E^{n + 1}_z(i, j + 1, k) - E^{n + 1}_z(i, j, k)}{\Delta y} - \frac { E^{n + 1}_y(i, j, k + 1) - E^{n + 1}_y(i, j, k)}{\Delta z} \right ] \\[2ex] {H^{n + 1}_{y}(i, j, k)} = CP \cdot {H^{n}_{y}(i, j, k)} - CQ \cdot \left [ \frac { E^{n + 1}_x(i, j, k + 1) - E^{n + 1}_x(i, j, k)}{\Delta z} - \frac { E^{n + 1}_z(i + 1, j, k) - E^{n + 1}_z(i, j, k)}{\Delta x} \right ] \\[2ex] {H^{n + 1}_{z}(i, j, k)} = CP \cdot {H^{n}_{z}(i, j, k)} - CQ \cdot \left [ \frac { E^{n + 1}_y(i + 1, j, k) - E^{n + 1}_y(i, j, k)}{\Delta x} - \frac { E^{n + 1}_x(i, j + 1, k) - E^{n + 1}_x(i, j, k)}{\Delta y} \right ] \end{cases}
⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧Hxn+1(i,j,k)=CP⋅Hxn(i,j,k)−CQ⋅[ΔyEzn+1(i,j+1,k)−Ezn+1(i,j,k)−ΔzEyn+1(i,j,k+1)−Eyn+1(i,j,k)]Hyn+1(i,j,k)=CP⋅Hyn(i,j,k)−CQ⋅[ΔzExn+1(i,j,k+1)−Exn+1(i,j,k)−ΔxEzn+1(i+1,j,k)−Ezn+1(i,j,k)]Hzn+1(i,j,k)=CP⋅Hzn(i,j,k)−CQ⋅[ΔxEyn+1(i+1,j,k)−Eyn+1(i,j,k)−ΔyExn+1(i,j+1,k)−Exn+1(i,j,k)]
注意,上式中的
H
n
H^n
Hn和
E
n
E^n
En是相差半个时间步长的,如下图:
这样,磁场和电场就可以随着时间交替求解。