FDTD（时域有限差分）这名字听起来挺唬人，但把它拆开看：有限差分，就是把麦克斯韦方程组里那些偏微分换成加减乘除；时域，就是一格一格往前推进时间。用 MATLAB 写一个基于 FDTD 的二维电磁波传播模拟，本质上就是搭一个数值“试验台”，观察波在自由空间中怎么冒出来、往外跑、最后被边界“吞”掉。我这次做的这个小项目，正是把整个过程用最直接的方式跑通了：TMz 模式、二维自由空间、点源激励、正弦脉冲，边界用 Mur 一阶 ABC 吸收，可视化用 MATLAB 的 pcolor 逐帧刷新。
项目本身不复杂，但把“原理—离散—边界—源—可视化”这条链路完整走一遍之后，你对 FDTD 的理解会和只看公式完全不一样。想做电磁仿真入门、需要演示电磁波传播场景，或者想在课程设计里加一个有物理内涵的 MATLAB 程序的同学，这篇文章应该能帮你省不少时间。不需要什么高端硬件，一台普通笔记本，装好 MATLAB，按下面的代码跑一遍就能看到波从中心扩散、接触边界、被“吸收”的全过程。
1. 项目概述——FDTD 自由空间模拟到底在做什么
1.1 技术选型背后的逻辑
先看模式选择。二维问题里最简单也最常用的是 TMz 模式：电场只有一个垂直分量 (E_z)，磁场有两个水平分量 (H_x)、(H_y)。好处是场方程只有三个标量方程，不会引入矢量场的复杂耦合，出图也直观。三维 FDTD 的实现差异主要是分量多了，核心思想完全一致。所以用 TMz 入门，性价比非常高。
再看边界条件。标题里写的 ABC 是吸收边界条件，我选用的是 Mur 一阶公式。为什么不用理想导体边界？因为它会把出射波反射回来，模拟出来的就不是自由空间而是谐振腔。为什么不用更先进的 PML？因为 Mur 一阶只有一种经典的单向波公式，代码量小、物理概念清楚，作为学习 FDTD 的边界处理是非常合适的起点。我没让问题一开始就变得复杂，这也是很多初学者最容易犯的毛病。
1.2 正弦脉冲激励源的特点
源选正弦脉冲，不是说随便挑一个正则的 sin 函数。它的本质是一个有限持续时间的窄带信号：载波频率定为 1 GHz，持续 5 个周期，5 ns 后源自动关闭。相比高斯脉冲的宽频特性，正弦脉冲的主要成分集中在单一频率附近，非常方便观察单频波在自由空间的传播和相速度表现。做演示时，你能清清楚楚看到一圈一圈的等相位面往外扩散，这种视觉反馈对理解 FDTD 很有帮助。
这套组合方案的关键优势是“每一块都能独立验证”：笛卡尔网格下的场更新是标准教科书内容，Mur 边界是标准教科书内容，正弦脉冲点源也是标准教科书内容。三个标准内容拼在一起，就能组装出一个完整的开放边界电磁波传播演示系统。后续想加介质、加散射体、加损耗，只需要在局部格点上改参数，不用重构主循环。
2. FDTD 原理：麦克斯韦方程是怎么变成递推公式的
2.1 从旋度方程到三项更新式
自由空间中，麦克斯韦方程的时域形式可以写作：
[ \nabla \times E = -\mu_0 \frac{\partial H}{\partial t} ]
[ \nabla \times H = \varepsilon_0 \frac{\partial E}{\partial t} ]
在二维 TMz 模式下，电场只保留 (E_z)，磁场只保留 (H_x)、(H_y)，展开后得到三个标量方程：
[ \frac{\partial H_x}{\partial t} = -\frac{1}{\mu_0} \cdot \frac{\partial E_z}{\partial y} ]
[ \frac{\partial H_y}{\partial t} = \frac{1}{\mu_0} \cdot \frac{\partial E_z}{\partial x} ]
[ \frac{\partial E_z}{\partial t} = \frac{1}{\varepsilon_0} \cdot \left( \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} \right) ]
剩下的事情就是把偏导数用差分替代。这里需要引入 Yee 网格：电场 (E_z) 放在整数网格点上，磁场 (H_x) 放在 y 方向半格处，磁场 (H_y) 放在 x 方向半格处，时间上也错开半步。这样做的好处是，空间差分立刻变成相邻两点的差，时间递推则是标准的蛙跳格式。
用索引写出来，每一步就非常机械了。假设 (E_z(i,j)) 表示整数点，(H_x(i,j)) 表示 ((i, j+1/2)) 的位置，(H_y(i,j)) 表示 ((i+1/2, j)) 的位置，那么在一个时间步 (dt) 里，三个更新式分别是：
[ H_x^{n+1/2}(i,j) = H_x^{n-1/2}(i,j) - \frac{\Delta t}{\mu_0 \Delta y} \left[ E_z^n(i,j+1) - E_z^n(i,j) \right] ]
[ H_y^{n+1/2}(i,j) = H_y^{n-1/2}(i,j) + \frac{\Delta t}{\mu_0 \Delta x} \left[ E_z^n(i+1,j) - E_z^n(i,j) \right] ]
[ E_z^{n+1}(i,j) = E_z^n(i,j) + \frac{\Delta t}{\varepsilon_0} \left[ \frac{H_y^{n+1/2}(i,j) - H_y^{n+1/2}(i-1,j)}{\Delta x} - \frac{H_x^{n+1/2}(i,j) - H_x^{n+1/2}(i,j-1)}{\Delta y} \right] ]
注意公式里的正负号，这是最容易出错的地方。(H_x) 的更新带负号，(H_y) 带正号，(E_z) 内部的旋度项是 (H_y) 对 x 的偏导减去 (H_x) 对 y 的偏导。很多初学者把符号搞反，结果波不但不往外传播，反而原地发散，看上去就像系统爆炸。
2.2 稳定性条件：为什么时间步和空间步必须配好
FDTD 是显式格式，存在数值稳定性上限。二维均匀网格的 Courant 条件给的是：
[ \Delta t \le \frac{1}{c_0 \sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}}} ]
当 (\Delta x = \Delta y = d) 时，简化成：
[ \Delta t \le \frac{d}{c_0 \sqrt{2}} ]
我的代码里直接取等号，即dt = dx / (c0 * sqrt(2))。这是很多教材推荐的做法，等于完全利用稳定边界；实际工程里有人会留 10% 余量，取 0.9 倍，避免数值色散积累。
除稳定性之外，空间分辨率也必须考虑。一个波长至少要有 10 到 20 个网格点，否则波的形状会严重失真。我这里 (f_0 = 1\text{ GHz})，波长 (0.3\text{ m})，而 (dx = 0.01\text{ m})，每个波长有 30 个网格点，绰绰有余。如果频率提高十倍，(dx) 就要缩小十倍，时间步也跟着变小，计算量会指数增加。这就是高频模拟为什么贵的根本原因。
2.3 一阶 Mur 边界条件的物理意义
再回来看 ABC。所谓吸收边界，本质是在截断边界上构造一个“只让波出去、不让波回来”的数值条件。Mur 一阶假设边界上的外向波满足单向波动方程，比如左边界 (x=0) 处有：
[ \frac{\partial u}{\partial t} - c_0 \frac{\partial u}{\partial x} = 0 ]
离散化后，可以解出边界上的新值。设 (s = c_0 \cdot \Delta t / \Delta x)，令 (r = (s-1)/(s+1))，则左边界更新为：
[ E_z^{n+1}(1,j) = E_z^{n}(2,j) + r \cdot \left[ E_z^{n+1}(2,j) - E_z^{n}(1,j) \right] ]
右边界、下边界、上边界形式类似。系数 (r) 在这里是负的，因为 (s = 0.707 < 1)，所以 (r \approx -0.1716)。边界更新本质上不是简单的“把边界清零”，而是根据内部的新值和旧值外推，让波近似“走出去”。
必须承认，一阶 Mur 对正入射波吸收效果最好，斜入射时会有一定反射。后文我会专门讲怎么识别和评估这种反射。它毕竟不是 PML，不能指望在所有入射角下都完美无缺。
3. MATLAB 代码实现与关键细节
3.1 参数设置：每一格该怎么选
写代码的第一步是定义物理参数和网格。代码里我给出了一组可以直接跑的默认参数，但你要清楚每个数的含义，改的时候才不会乱。
%===================================================== % FDTD 2D TMz 模拟：自由空间电磁波传播 % 吸收边界：Mur 一阶 ABC % 激励源：正弦脉冲（点源） %===================================================== clear; clc; close all; %% 物理常数与网格参数 c0 = 3e8; % 光速 m/s eps0 = 8.854e-12; % 真空介电常数 F/m mu0 = pi*4e-7; % 真空磁导率 H/m dx = 0.01; % 空间步长 x (m) dy = dx; % 空间步长 y (m) dt = dx / (c0 * sqrt(2)); % 时间步长 (Courant 极限) Nx = 300; % x 方向网格数量 Ny = 300; % y 方向网格数量 Nt = 900; % 总时间步数 % 物理坐标矢量 x_vec = (0:Nx-1) * dx; y_vec = (0:Ny-1) * dy; %% 场分量初始化 Ez = zeros(Nx, Ny); % 电场 z 分量 Hx = zeros(Nx, Ny); % 磁场 x 分量 Hy = zeros(Nx, Ny); % 磁场 y 分量 %% 激励源参数 f0 = 1e9; % 载波频率 1 GHz T0 = 1 / f0; % 周期 ncycle = 5; % 脉冲持续 5 个周期 t_source = ncycle * T0; % 源持续时间 % 源位置：取网格几何中心 src_x = round(Nx / 2); src_y = round(Ny / 2); %% Mur ABC 系数 s = c0 * dt / dx; % 当前 Courant 数 = 1/sqrt(2) r = (s - 1) / (s + 1); % Mur 一阶边界系数c0、eps0、mu0是真空常数；dx决定空间分辨率；dt由 Courant 条件推出；Nx、Ny决定计算区域大小；Nt是总时间步数。源参数f0、ncycle、t_source定义了激励信号。src_x、src_y是源的位置，我放在网格中心，这样图看起来最对称。
注意一个映射细节：数组索引和物理坐标的对应关系是 (x = (i-1)\Delta x)，(y = (j-1)\Delta y)。如果你把索引直接当坐标，等画图时标注的横轴就会差一个 (dx)，虽然整体分布不变，但坐标刻度就是错的。这一点在分析结果时很关键。
3.2 场更新：主循环里每一步到底发生了什么
先看主循环里的 Hx、Hy、Ez 三个更新式，以及源的注入位置。
%% 主时间步循环 figure('Color', 'w', 'Position', [80 80 640 540]); for n = 1:Nt t_advance = (n - 0.5) * dt; % 半时间步时刻 % --- 更新 Hx：位置 (i, j+1/2) --- Hx(:, 1:end-1) = Hx(:, 1:end-1) - (dt / (mu0 * dy)) * ... (Ez(:, 2:end) - Ez(:, 1:end-1)); % --- 更新 Hy：位置 (i+1/2, j) --- Hy(1:end-1, :) = Hy(1:end-1, :) + (dt / (mu0 * dx)) * ... (Ez(2:end, :) - Ez(1:end-1, :)); % --- 保存更新前的 Ez，Mur 边界需要上一时刻边界值 --- Ez_old = Ez; % --- 更新 Ez 内部区域 --- Ez(2:end-1, 2:end-1) = Ez(2:end-1, 2:end-1) + (dt / eps0) * ... ( (Hy(2:end-1, 2:end-1) - Hy(1:end-2, 2:end-1)) / dx - ... (Hx(2:end-1, 2:end-1) - Hx(2:end-1, 1:end-2)) / dy ); % --- 注入正弦脉冲点源（加性源） --- if t_advance <= t_source Ez(src_x, src_y) = Ez(src_x, src_y) + sin(2 * pi * f0 * t_advance); end % --- Mur 一阶 ABC 四边更新 --- Ez(1, :) = Ez_old(2, :) + r * (Ez(2, :) - Ez_old(1, :)); Ez(end, :) = Ez_old(end-1, :) + r * (Ez(end-1, :) - Ez_old(end, :)); Ez(:, 1) = Ez_old(:, 2) + r * (Ez(:, 2) - Ez_old(:, 1)); Ez(:, end) = Ez_old(:, end-1) + r * (Ez(:, end-1) - Ez_old(:, end)); % --- 每隔 20 步绘制快照 --- if mod(n, 20) == 0 pcolor(x_vec, y_vec, Ez'); shading interp; caxis([-0.3 0.3]); colormap(jet); colorbar; xlabel('x / m'); ylabel('y / m'); title(sprintf('Ez 分布 t = %.2f ns', t_advance * 1e9)); drawnow; end endHx、Hy、Ez三个更新式要按顺序执行：先更新磁场，再更新电场，因为蛙跳格式中，电场需要用刚刚更新过的磁场。这是我这里的具体顺序，如果你以后用其他库或者换语言，也要保持这种时间交错的思路。
更新完内部 Ez 后，立刻注入激励源。注意源是在整个场更新之后才加入的，这是显式 FDTD 添加点源的常见做法之一。实际物理上源也可以视为某个格点上的场值被外部驱动，放在更新后不会破坏差分格式的时间一致性。这里有个细节：我没有用“硬源”，而是直接把正弦函数加到Ez(src_x, src_y)上，这叫加性源。硬源会直接覆盖场值，入射波和散射波在源位置无法叠加；加性源则允许场自由通过源点，更接近物理真实。
提示：场更新顺序一定要是先磁场后电场。如果你先更新 Ez 再更新 Hx、Hy，整个蛙跳时间结构就塌了，模拟结果会偏得离谱。
3.3 Mur ABC：四边边界更新为什么要在电场更新之后做
需要特别强调的是，Mur 边界的实现顺序不能乱。我用Ez_old保存了电场更新前的值；然后更新内部 Ez；紧接着，在边界上用Ez_old和更新后的内部值算出边界新值。为什么必须这样？因为 Mur 公式同时涉及上一时刻的边界值和当前时刻的内部值。如果你先更新边界再更新内部，边界上新值会被内部更新覆盖；如果你不保存旧值，直接拿更新后的内部值当“旧值”用，公式就错了。
左边界和右边界关注 x 方向，下边界和上边界关注 y 方向。四个角落的计算可能被两条边重复更新，我的做法是让后执行的边界决定角点值。严格来说角点应单独用一维外向波方程处理，但一阶方案下误差其实很小，多数演示场景可以忽略。
边界更新的四条语句里，索引变化很有规律：
- 左边界：取出原数组第二行作为旧内部值，用更新后的第二行做外推，边界行是第一行；
- 右边界：用倒数第二行外推边界行，边界行是最后一行；
- 下边界：用第二列外推边界列，边界列是第一列；
- 上边界：用倒数第二列外推边界列，边界列是最后一列。
写错索引最常见的后果是：边界上出现镜像条纹，或者波一碰到边界就产生剧烈振荡。如果看到这种症状，先把这四条语句的索引逐行打印出来，对照物理坐标图看一遍。
3.4 动态可视化需要注意什么
% 循环结束后的最终静态图 figure; subplot(1,2,1); pcolor(x_vec, y_vec, Ez'); shading interp; colorbar; axis equal tight; title('Ez 空间分布'); xlabel('x / m'); ylabel('y / m'); subplot(1,2,2); plot(x_vec, Ez(:, src_y), 'b-', 'LineWidth', 1.2); grid on; xlabel('x / m'); ylabel('Ez'); title('沿 y 中心线的切片');pcolor画的是伪色彩图，注意转置：因为Ez是按Nx * Ny存储的，x 对应第一维，y 对应第二维，pcolor的输入是 x、y、Z，其中 Z 是Ny * Nx，所以要写Ez'。我每 20 步刷新一次，否则每步都drawnow会让整个循环变得很慢。颜色范围固定下来，不要让它自动缩放，不然你会看到每一步的色标都在跳，根本没法观察波前的衰减。
如果你觉得逐帧动画还不过瘾，可以在循环内把Ez存成三维数组，也可以在drawnow同时用getframe录制成视频。不过这会让内存占用骤增，我的建议是先用快照方式观察，确认现象正常后再考虑录动画。
4. 模拟结果与现象解读
4.1 波前扩散的完整过程
如果把前 200 个时间步的数据保存下来回放，看到的应该是：源点首先出现一个逐渐增大的正负交替环形波前，波峰和波谷像同心圆一样向外扩散。因为是单频正弦脉冲，波前是一圈一圈的；而在脉冲结束后，源位置不再有新信号，外围波前继续向前传播。
计算一下时间尺度：(dt \approx 23.6\text{ ps})，200 步约 (4.7\text{ ns})。此时 1 GHz 波完成了大约 4.7 个周期，加上前导波和边界之间还有距离，波没有碰到边界。做到大约 400～500 步，波前开始接触计算区域边缘，这时最能检验 ABC 的效果。
需要提醒的是，正弦脉冲虽然是窄带信号，但它的瞬时开通仍然会带来一个“前导波头”。这个波头从零突然跳到正弦的某个相位，空间上表现为一个较陡的过渡带，它比纯正弦波更容易产生数值色散。如果你把dx调大，会明显看到波头变得圆钝、拖尾变长，这就是数值色散的直接证据。
4.2 边界反射怎么判断
判断 ABC 好不好，有一个很直接的观察法：波前碰到边界后，边界内侧不应该出现明显的“二次波纹”再往回流。如果看到一圈清晰的反射波从边界折回中心，那说明边界条件没有正常工作，或者实现有误。
我测试时用的是 300×300 的网格，900 个时间步足够让波到达边界并大部分吸收。尾部还有一点残余反射是正常现象——一阶 Mur 对斜入射本来就有局限，圆形波前到达矩形边界时，除正对源点的区域外，大部分都是斜入射。后文我会专门讲如何进一步削弱反射。
另一个容易误判的点：波的柱面扩散本身就会让振幅随距离衰减。二维自由空间中，理想柱面波的电场幅度按 (1/\sqrt{R}) 衰减。所以看到靠近边界的波峰比中心位置低，并不代表能量被吸收了，而是几何衰减的必然结果。很多人第一次看到这个现象会误以为边界有问题，先记下这个物理规律，能省去不少排查时间。
4.3 观察相速度和频率的关系
正弦单频波在均匀网格里传播时，数值相速度和真实相速度会有细微差别，这叫数值色散。二维 FDTD 的色散关系比较长，这里我不展开公式，只讲观察方法：你可以在源点附近某条固定直径方向放几个探针点，记录 (E_z) 随时间变化，然后求出相邻探针点之间的相位差，反推数值相速度。通常情况下，每波长 30 格时，数值相速度和 (c_0) 的偏差在 1% 以内，肉眼看图基本分辨不出来。每波长低于 10 格时，偏差会明显到让波前形状都开始扭曲。
5. 参数敏感性与常见问题排查实录
5.1 空间步长对精度的影响
如果dx等于波长的 1/5，数值色散会极其明显，波的相速度不再等于 (c_0)，而且不同方向传播的相位误差不一样，原本的圆波前会变成多边形。有人把这种现象叫“各向异性色散”。如果你看到波前变形，先查空间分辨率：每波长至少保证 10～20 个网格点。
这里有个经验值：想看清干涉结构或驻波特征时，每波长 20 格以上比较稳妥；如果只是演示波传播趋势，10 格也能看。但千万别认为 10 格就“够用”了，因为 FDTD 的误差是累积的，跑得越久，色散导致的相位误差越明显。
5.2 Courant 条件不满足的典型症状
把dt取成dx / c0（不带根号二）就会违反二维 Courant 条件。短时间看不出问题，几十步之后场值开始指数膨胀，图像出现高频闪烁，最后直接变成 NaN。这不是边界或者源的问题，而是时间步太大。遇到这种情况，立刻回头检查dt的计算公式。
在排查时可以把dt再缩小 10%，观察是否稳定。如果缩小后数值稳定、波形正常，说明就是 Courant 条件踩线了。显式格式的一个特点是“要么稳定，要么爆掉”，不会给你中间过渡态。
5.3 Mur 边界总是有一点反射怎么办
首先确认公式和更新顺序没有错。如果确认无误，还想进一步减少反射，可以升级到 Mur 二阶，或者换上 PML。PML 是匹配层思想，通过在边界区引入人工各向异性介质吸收波，效果比 Mur 好很多，但实现代码要多好几行。我建议先把 Mur 跑通，再去看 PML，这样你能真正理解为什么要做匹配层。
Mur 二阶比一阶多了一项关于切向二阶导数的修正，边界附近需要额外存储一行一列的数据来构造二阶差分。这个升级不复杂，但代码美观度会下降。如果你只是做自由空间演示，一阶 Mur 已经足够；如果要研究散射体，那确实应该考虑 PML。
5.4 常见问题速查表
| 现象 | 可能原因 | 排查思路 |
|---|---|---|
| 场值几个时间步就发散 | dt 超过 Courant 条件 | 检查 dt 是否满足二维不等式 |
| 波前明显变形不是圆 | dx 太大，数值色散严重 | 加密网格，保证每波长 ≥ 20 格 |
| 边界出现强反射回流 | Mur 实现顺序或符号错 | 重新核对四条边的索引和 Ez_old 使用位置 |
| 源点不出现波 | 源位置数组下标越界 | 确认 src_x、src_y 在 2 到 Nx-1 内 |
| 动画太卡 | 每步都 drawnow | 每隔 20 或 50 步刷新一次 |
| 图像上下翻转 | 忘记转置 Ez | pcolor 里用 Ez'，因为存储顺序是 Nx×Ny |
| 波碰到边界后明显炸开 | 边界更新被内部更新覆盖 | 确认 Mur 更新在 Ez 内部更新之后执行 |
5.5 调试用的三个小技巧
把Nx、Ny先设成 50×50 的小网格，跑通再放大。加一个探针点记录Ez随时间变化，直接观察波形。把边界换成 PEC 跑一遍对照组，再和 ABC 结果对比，立刻明白吸收边界的价值。这是我实际调试时最常用的三个手段。
具体来说，PEC 对照组就是把 Mur 边界那四条更新语句注释掉，让边界值保持为 0。跑同样的源和时间步，你会看到波碰到边界后产生非常干净的镜面反射，反射波和原波相互叠加，最后形成复杂的驻波图案。再看 ABC 版本几乎滚不出来反射波，这么一对比，边界条件的贡献就一目了然了。
还有个小技巧：检查边界吸收率时，不一定要用完整的大网格。可以做一个窄长条网格，比如 60×200，把源放在一端，观察波沿线传播到另一端后的剩余幅度。这样能清楚地看到一次边界反射的幅度，比看二维平面图更直观。
6. 进一步扩展方向
这个项目做完之后，往几个方向扩展都很顺手：把源从单点改成两个同相或反相点源，就能看到干涉条纹；把局部介电常数改成一组数值，就成了介质柱散射模拟；把二维升级成三维，核心代码结构基本不变，只是多出几个场分量。
我个人觉得最划算的扩展是加一个简单的介质圆盘：在网格区域内把某个圆内的介电常数改成 2.25（水或者玻璃的典型值），然后重新跑同样源。你会看到波在圆盘边界发生部分反射和透射，透射波在圆盘内部波长被压短，因为介质中波长 (\lambda = \lambda_0 / \sqrt{\varepsilon_r})。这个现象直接用肉眼就能验证，对理解介质中的电磁波非常有帮助。
我实际跑下来最大的体会是：FDTD 的参数选择和实现顺序同等重要。写代码之前先把时间步、空间步、波长关系算清楚，比写完再调参要高效得多。你如果想自己动手试，建议第一步就把Nx、Ny降到 100 以下，把绘图间隔调大，专注观察波前和边界反射，然后再慢慢加大规模。跑通一次之后，你会发现电磁场模拟并不是什么神秘的东西，无非就是把麦克斯韦方程“翻译”成循环和数组操作，剩下的事就是让计算机按部就班地算下去。
