要彻底弄懂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是相差半个时间步长的，如下图：
这样，磁场和电场就可以随着时间交替求解。 我用C++实现了上述过程，点击此处查看。
