FDTD_1

要彻底弄懂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​​−σm​Hx​∂z∂Ex​​−∂x∂Ez​​=−μ∂t∂Hy​​−σm​Hy​∂x∂Ey​​−∂y∂Ex​​=−μ∂t∂Hz​​−σm​Hz​​
为了后续推导的方便,我们首先定义
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ϵσΔt​1−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​Δt​1−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σm​​1​=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是相差半个时间步长的,如下图:
在这里插入图片描述

这样,磁场和电场就可以随着时间交替求解。