3. 一维FDTD方法基础

时域有限差分（FDTD）方法的核心思想非常直接：把麦克斯韦旋度方程中的空间偏导和时间偏导都用中心差分近似，然后在离散网格上交替推进电场与磁场。一维情形是理解整套方法的最佳入口——它保留了Yee格式的全部关键特征（交错采样、蛙跳推进、稳定性约束），却把计算量压缩到可以手算验证的程度。本节以自由空间中沿 x 方向传播的均匀平面波为例，建立一维FDTD的完整框架，并给出可直接运行的程序结构与验证方法。

Yee网格：电场与磁场的空间交错

在一维问题中，假设场量只随 x 变化，且电场沿 y 方向、磁场沿 z 方向（TEM波），则旋度方程退化为两个标量方程：

∂E_y/∂t = -(1/ε₀) · ∂H_z/∂x
∂H_z/∂t = -(1/μ₀) · ∂E_y/∂x

Yee网格的关键在于：E_y 与 H_z 不在同一空间位置采样，也不在同一时刻采样。电场节点位于整数网格点 x = iΔx，磁场节点位于半整数点 x = (i+1/2)Δx；时间上电场取整数步 t = nΔt，磁场取半整数步 t = (n+1/2)Δt。这种交错使每个中心差分恰好跨越一个完整的网格间距，二阶精度得以保证，同时天然满足离散形式的法拉第定律与安培定律。

图中蓝色圆点为电场采样点，橙色方块为磁场采样点。空间上二者相差半个网格，时间上相差半个时间步，这正是“蛙跳”（leapfrog）格式的几何来源。

时间步进更新公式

对上述两个方程分别做中心差分，得到一维FDTD的标准更新式。为书写简洁，记 E_iⁿ 表示 E_y(iΔx, nΔt)，H_i+1/2^n+1/2 表示 H_z((i+1/2)Δx, (n+1/2)Δt)：

H[i+1/2]^(n+1/2) = H[i+1/2]^(n-1/2) - (Δt/(μ₀Δx)) · (E[i+1]^n - E[i]^n)

E[i]^(n+1) = E[i]^n - (Δt/(ε₀Δx)) · (H[i+1/2]^(n+1/2) - H[i-1/2]^(n+1/2))

注意两个系数的量纲：Δt/(μ₀Δx) 与 Δt/(ε₀Δx) 分别对应磁场和电场的更新增益。在自由空间中，若定义 Courant 数 S = c₀Δt/Δx，则两系数可统一写成 S/(Z₀Δx) 形式（Z₀ = √(μ₀/ε₀) ≈ 377 Ω），这也是程序中常用的归一化写法。更新顺序必须是：先用第 n 步的电场更新第 n+1/2 步的磁场，再用新磁场更新第 n+1 步的电场，如此循环。

Courant稳定性条件

显式差分格式并非无条件稳定。对一维Yee格式做von Neumann分析可得，数值解不发散的条件为：

S = c₀ · Δt / Δx ≤ 1

这就是一维Courant（CFL）条件。其物理含义是：在一个时间步内，信息传播的距离不能超过一个空间网格，否则差分模板“追不上”真实波前，数值解会指数增长。工程上通常取 S = 0.99 或 0.5：取接近1可减少数值色散，但余量小；取0.5则更稳健，代价是时间步数翻倍。三维情形对应 c₀Δt ≤ 1/√(1/Δx²+1/Δy²+1/Δz²)，一维是其最简特例。

数值色散： Courant 数如何影响波形

即使满足稳定性条件，离散网格上的波速也与真实光速略有差异，且随频率变化，这就是数值色散。对一维Yee格式可推导出数值色散关系：

sin(ωΔt/2) = S · sin(kΔx/2)

当 Δx → 0 且 Δt → 0 时该式退化为 ω = c₀k。有限网格下，高频分量（每波长采样点少）的相速偏差更大。实践中的经验判据是：每个波长至少取 10～20 个网格点，此时数值色散引起的波形展宽通常可控制在百分之几以内。当 S = 1 时，一维格式对沿轴向传播的波恰好无色散（这是1D的特殊性质），这也是许多一维验证程序偏爱取 S = 1 的原因；但实际多维问题中不能取到该极限。

程序实现：自由空间高斯脉冲传播

下面给出一个可直接运行的一维FDTD核心循环（Python/NumPy风格伪代码，参数为教学示例值）。设 Δx = 1 mm，c₀ ≈ 3×10⁸ m/s，取 S = 0.9，则 Δt = S·Δx/c₀ = 3×10⁻¹² s。计算域取 800 个网格，两端各留 100 格作为吸收区（简单一阶Mur吸收边界或PML均可，本例用Mur一阶）。

import numpy as np

c0 = 3e8; dx = 1e-3; S = 0.9
dt = S*dx/c0
eps0 = 8.854e-12; mu0 = 4e-7*np.pi
N = 800; T = 1200

E = np.zeros(N); H = np.zeros(N)   # H[i] 表示 i+1/2 处磁场
# 高斯脉冲源：中心频率对应约 20 格/波长
t0 = 30*dt; tau = 12*dt
for n in range(T):
    # 1) 更新磁场（内部点）
    H[:-1] -= (dt/(mu0*dx)) * (E[1:] - E[:-1])
    # 2) 更新电场（内部点）
    E[1:] -= (dt/(eps0*dx)) * (H[1:] - H[:-1])
    # 3) 软源注入（在 i=150 处加高斯脉冲）
    E[150] += np.exp(-((n*dt - t0)/tau)**2)
    # 4) 一阶Mur吸收边界（简化示意）
    E[0]  = E[1]  + (S-1)/(S+1)*(E[1]-E[0])
    E[-1] = E[-2] + (S-1)/(S+1)*(E[-2]-E[-1])

运行后应观察三个现象：第一，脉冲以接近 c₀ 的速度向右传播，可用两个时刻的峰值位置差除以时间差估算数值相速；第二，脉冲形状基本保持，若把 Δx 增大到每波长仅 5 格，会看到明显的拖尾与展宽，这就是数值色散；第三，把 S 改为 1.05 重跑，场值会迅速发散，验证CFL条件的必要性。

验证清单与判据

为了让“验证数值色散与稳定性”落到实处，建议按以下三项做定量检查，而不是仅凭肉眼观察波形：

• 稳定性验证：分别取 S = 0.5、0.9、1.0、1.05，记录 1200 步后场值最大绝对值。前三者应保持在 O(1) 量级，S = 1.05 应出现指数增长。
• 相速验证：追踪脉冲峰值到达两个固定探针（如 i=300 与 i=600）的时刻差 Δt_p，数值相速 v_p = 300Δx/Δt_p，与 c₀ 比较，相对误差通常应在 1% 以内（S 接近1时更小）。
• 色散验证：固定 S = 0.9，把每波长采样数从 20 降到 10、5，测量脉冲宽度（半高全宽）的相对展宽比例，应随采样数减少而单调增大。

总结本节要点：Yee网格用空间与时间的双重交错换取二阶精度和天然守恒性；蛙跳更新式决定了必须先场后场、交替推进；CFL条件 S ≤ 1 是显式格式的硬约束；数值色散则要求每个波长有足够采样。把这三件事在一维高斯脉冲程序中亲手复现一遍，比记住任何公式都更有助于后续多维仿真的正确实施。