简介一套基于Matlab的时域有限差分FDTD仿真源码包覆盖一维、二维和三维电磁场模拟场景并配有动态图形界面适合电磁场数值计算初学者、高校相关专业学生以及需要快速搭建FDTD模型进行科研验证的工程师。压缩包共10个文件以6个.m主程序为核心分别实现从一维到三维的差分迭代计算另含.fig图形用户界面、.asv自动备份及说明文档整体仅87KB轻量易用可直接运行。目前已有703人学习下载。通过运行这套代码读者可以逐步理解FDTD的网格划分、边界条件与时间步进逻辑动态观察电磁场在简单波导、平面结构及复杂三维模型中的传播与演化为后续研究光子学器件、天线设计或微波电路打下扎实基础。1. 从一维到三维的 FDTD 仿真这份 Matlab 代码包能帮你少走三个月弯路做电磁场仿真的人早晚都会撞上 FDTD时域有限差分这道坎。无论是天线设计、波导分析还是超材料结构仿真FDTD 都是绕不开的基础工具。但网上散落的 FDTD 代码大多只覆盖单一维度而且注释稀碎想从一维推到三维往往要自己啃好几个月的论文。这份 fdtd.rar 资源包的价值就在于它把一维、二维、三维的 FDTD 实现都放在了一起而且是用 Matlab 写的打开就能跑改改参数就能用。对于刚接触 FDTD 的学生或者需要快速验证电磁方案的工程师来说这份代码能直接帮你跳过从理论到代码实现之间那段最容易放弃的路。接下来我按维度拆开讲每一段都会给出能直接复现的代码和参数说明把我踩过的坑也一并交代清楚。2. 一维 FDTD从麦克斯韦方程到第一行可运行代码2.1 一维 FDTD 的核心迭代逻辑一维 FDTD 是所有维度里最接近教科书的一版它只考虑电磁波沿单一方向传播电场和磁场交替更新。核心思路是 Yee 网格电场分量在整数时间步上更新磁场分量在半时间步上更新空间上也错开半个网格。这种蛙跳格式保证了数值稳定性也让它实现起来特别直白。这里我用自由空间传播做一个最简示例代码里没有加入源项和边界条件让你看清楚迭代内核到底长什么样。% 1D FDTD 自由空间传播核心迭代 % 物理参数 dx 0.01; % 空间步长单位m dt dx / (2 * 3e8); % 时间步长满足CFL条件dt dx / c Nx 500; % 网格总数 Nz 600; % 时间步数 % 场初始化 Ez zeros(1, Nx); % 电场 z 分量 Hy zeros(1, Nx); % 磁场 y 分量 % 材料参数真空 eps0 8.854e-12; mu0 1.2566e-6; ca ones(1, Nx); % 电场更新系数自由空间为1 cb dt / (eps0 * dx) * ones(1, Nx); % 电场更新中磁场贡献的系数 da ones(1, Nx); % 磁场更新系数 db dt / (mu0 * dx) * ones(1, Nx); % 磁场更新中电场贡献的系数 % 主循环 for n 1:Nz % 更新磁场时间步 n0.5 for i 1:Nx-1 Hy(i) da(i) * Hy(i) db(i) * (Ez(i1) - Ez(i)); end % 更新电场时间步 n1 for i 2:Nx Ez(i) ca(i) * Ez(i) cb(i) * (Hy(i) - Hy(i-1)); end % 在网格中心加一个高斯脉冲源 Ez(250) Ez(250) exp(-0.5 * ((n - 60) / 12)^2); end代码逻辑不复杂但有两个地方要注意。磁场更新循环里索引从 1 到 Nx-1因为计算 Hy(i) 需要用到 Ez(i1)最右边的电场没有右邻居所以磁场在最右端不更新。电场更新循环里索引从 2 到 Nx因为 Ez(i) 需要用到 Hy(i-1)最左边的电场没有左邻居。这两个边界截断在后面会专门讨论现在先记住这个结构。参数选择上spatial step dx 决定空间分辨率一般取最小波长的 1/10 到 1/20。时间步 dt 必须满足 CFL 条件一维下就是 dt dx / c这里取了 0.5 倍的安全系数。系数 ca、cb、da、db 是从麦克斯韦旋度方程离散化推导出来的如果仿真区域里有介质ca 和 da 就不再是 1需要按相对介电常数和磁导率重新计算。2.2 激励源与吸收边界的配合一维 FDTD 跑起来之后你会发现波传到边界会反射回来这在大多数场景里都不是你想要的结果。解决办法是加吸收边界最简单的是一阶 Mur 边界。% 一阶 Mur 吸收边界实现接在上面的主循环内 % 在更新电场之前先保存左、右边界的前一步值 Ez_left_old Ez(2); Ez_right_old Ez(Nx-1); % ... 电场更新循环 ... % 应用 Mur 边界 Ez(1) Ez_left_old (dt * 3e8 - dx) / (dt * 3e8 dx) * (Ez(2) - Ez(1)); Ez(Nx) Ez_right_old (dt * 3e8 - dx) / (dt * 3e8 dx) * (Ez(Nx-1) - Ez(Nx));Mur 边界的本质是把外向行波的条件直接离散化让边界处的场值根据内侧邻居的当前值和上一步值推算出来。实现上要注意保存时序Ez_left_old必须是在电场更新之前保存的旧值如果你先更新了电场再存边界的吸收效果就废了。这个顺序问题我初学时踩过边界反射率直接飙到 5% 以上肉眼可见地离谱。激励源的位置和类型也影响仿真结果。高斯脉冲适合看宽频响应连续正弦波适合做窄带分析。上面代码里的源是硬源直接把场值赋进去它会同时向两个方向辐射。如果要模拟只往一个方向传播的波需要用总场散射场TF/SF方法一维实现就是在某个网格点左边和右边分别注入代码量会大一些但这份资源包里有现成实现直接用就行。3. 二维 FDTDTMz 模的完整实现与网格区域划分3.1 二维 Yee 网格与 TMz 模式的场分布二维 FDTD 比一维多了一个空间维度场分量变成三个Ez、Hx、Hy对应 TMz 偏振电场只有 z 方向分量磁场在 xy 平面内。这是二维 FDTD 最常处理的模式因为很多平面波导和光子晶体结构都是这个偏振。Yee 网格在二维的布局是Ez 在整数网格点上Hx 和 Hy 在半网格偏移的位置上具体拓扑关系决定了迭代公式里空间差分的方向。% 2D FDTD TMz 模式初始化 Nx 200; Ny 200; % 网格数 dx 0.005; dy 0.005; % 空间步长单位m dt dx / (2 * 3e8); % 2D CFLdt 1 / (c * sqrt(1/dx^2 1/dy^2)) steps 400; % 场分量 Ez zeros(Nx, Ny); Hx zeros(Nx, Ny); Hy zeros(Nx, Ny); % 材料图1为真空2为介质 material ones(Nx, Ny); material(80:120, 80:120) 2; % 在中心放一个正方形介质块 % 参数场 eps_r ones(Nx, Ny); eps_r(material 2) 4.0; % 介质块相对介电常数 mu_r ones(Nx, Ny); ca ones(Nx, Ny); cb dt / (eps0 * dx) ./ eps_r; da ones(Nx, Ny); db dt / (mu0 * dx) ./ mu_r;这段初始化里material矩阵用来标记不同区域后续可以根据区域切换材料参数。eps_r是相对介电常数分布介质块的折射率就是它的平方根。在实际仿真光通信器件时我会把折射率分布从外部文件读进来比如用一个二维数组存着器件截面这样改结构不用改代码。二维 CFL 条件和一维不同它要满足 dt 1 / (c * sqrt(1/dx^2 1/dy^2))如果 dx 和 dy 相等简化为 dt dx / (c * sqrt(2))。代码里用的 dt dx / (2c) 比上限更保守稳定性没问题但时间步数会变多。实际项目中我一般取 CFL 数的 0.9 倍兼顾精度和效率。3.2 主迭代循环磁场双分量更新% 2D TMz 主迭代循环 for n 1:steps % 更新 Hx依赖 Ez 在 y 方向的差分 for i 1:Nx for j 1:Ny-1 Hx(i, j) Hx(i, j) db(i, j) * (Ez(i, j) - Ez(i, j1)); end end % 更新 Hy依赖 Ez 在 x 方向的差分 for i 1:Nx-1 for j 1:Ny Hy(i, j) Hy(i, j) db(i, j) * (Ez(i1, j) - Ez(i, j)); end end % 更新 Ez依赖 Hx 和 Hy 的差分 for i 2:Nx for j 2:Ny Ez(i, j) ca(i, j) * Ez(i, j) cb(i, j) * (Hy(i, j) - Hy(i-1, j) Hx(i, j-1) - Hx(i, j)); end end % 添加高斯源 Ez(100, 100) Ez(100, 100) exp(-0.5 * ((n - 40) / 10)^2); % 可视化每10步显示一次 if mod(n, 10) 0 imagesc(abs(Ez)); axis equal; colorbar; title(sprintf(Step %d, n)); drawnow; end end注意看三个更新循环里索引范围各不相同这个对应的是 Yee 网格的半格偏移拓扑。Hx 更新中Ez(i,j) - Ez(i,j1)是负的 y 方向差分这是由旋度方程里电场和磁场的右手关系决定的。如果你把差分方向写反波就会往错误方向传播这在调试时会表现为场分布左右颠倒尤其诡异。可视化部分用imagesc(abs(Ez))转置是因为 Matlab 的矩阵是按列存储的转置后图像方向和空间坐标对得上。如果你不转置看到的仿真区域就是上下翻转的对结果判断会形成干扰。二维仿真的性能瓶颈在这个三重循环上NxNy200 时一个时间步要算 4 万次迭代400 步就是 1600 万次Matlab 跑起来大概需要几十秒到几分钟。这个资源包里用的是纯 Matlab 循环没有用 C MEX 加速好处是移植性极好任何机器上有 Matlab 就能跑。真要提速可以把磁场更新写成矩阵运算一次完成代码会短很多但物理过程的逐点推演逻辑就不那么直观了。4. 三维 FDTD六个场分量的更新矩阵与散射计算4.1 三维 Yee 网格的场分量排布三维 FDTD 是完整版的 FDTD六个电磁场分量Ex、Ey、Ez、Hx、Hy、Hz全部参与更新每一对电场和磁场在空间上相差半个网格。这意味着三维 FDTD 的内存消耗呈立方级增长200 的立方网格就要 800 万个场点每个点六个分量占内存差不多 400MB这还没算材料参数数组。所以三维仿真里网格大小和仿真区域的设计直接影响可行性。三维电场更新的核心公式是Ex(i,j,k) Ex(i,j,k) (dt/eps) * [(Hz(i,j,k) - Hz(i,j-1,k)) / dy - (Hy(i,j,k) - Hy(i,j,k-1)) / dz]这个公式看着长但结构很有规律每个电场分量由两个磁场分量的旋度差分驱动。Matlab 里实现时索引关系容易错我建议先手动画出 Yee 网格图标清每个场分量在哪条边和哪个面上。资源包里三段独立的更新循环我拆开写虽然代码行长但每段对应一个物理分量出了错也好定位。% 3D FDTD 主更新循环电场部分 % 假设已初始化 Ex, Ey, Ez, Hx, Hy, Hz 六分量三维数组尺寸 (Nx,Ny,Nz) % 材料参数 eps_r, mu_r 同尺寸 % 取局部变量加速访问 ce1 dt / eps0; % 通用系数 % 更新 Ex依赖 Hz 的 y 差分 和 Hy 的 z 差分 for k 2:Nz for j 2:Ny for i 1:Nx Ex(i,j,k) Ex(i,j,k) (ce1 / eps_r(i,j,k)) * ( ... (Hz(i,j,k) - Hz(i,j-1,k)) / dy - ... (Hy(i,j,k) - Hy(i,j,k-1)) / dz ); end end end % 更新 Ey依赖 Hx 的 z 差分 和 Hz 的 x 差分 for k 2:Nz for j 1:Ny for i 2:Nx Ey(i,j,k) Ey(i,j,k) (ce1 / eps_r(i,j,k)) * ( ... (Hx(i,j,k) - Hx(i,j,k-1)) / dz - ... (Hz(i,j,k) - Hz(i-1,j,k)) / dx ); end end end % 更新 Ez依赖 Hy 的 x 差分 和 Hx 的 y 差分 for k 1:Nz for j 2:Ny for i 2:Nx Ez(i,j,k) Ez(i,j,k) (ce1 / eps_r(i,j,k)) * ( ... (Hy(i,j,k) - Hy(i-1,j,k)) / dx - ... (Hx(i,j,k) - Hx(i,j-1,k)) / dy ); end end end这段代码的关键点是每个电场分量的索引范围不一样。Ex 在 k 和 j 方向从 2 开始i 方向从 1 开始这是因为 Ex 分量所在的 Yee 网格位置在 i 方向是整数点在 j 和 k 方向是半格点。如果所有循环都写成从 2 到 N计算结果在边界处就是错的而且这种错不会崩只是场值慢慢发散特别难排查。三维代码里eps_r出现在除法分母上如果某个网格点设成 0直接会出现 inf 或 NaN整个仿真就毁了。我遇到过一种诡异情况介质的几何结构用 CAD 模型导进来边界处产生了一层厚度为零的网格单元eps_r被设成 0然后仿真到几千步时全部场值变成 NaN。排查到最后是拿三维切片图逐层看才发现那一层网格点颜色完全不对。从那以后我每次导入材料分布都会做一次非零检查。4.2 三维散射计算近场到远场变换三维仿真的最终结果往往要算远场方向图比如天线辐射方向或者雷达散射截面RCS。FDTD 直接算出来的是近场分布需要用等效原理把近场变换到远场。这个变换在代码包里单独放了一个函数文件核心流程在包围目标的闭合面上记录电场和磁场的切向分量然后通过等效电流和等效磁流在频域里积分得到远场方向图。% 近场到远场变换示意代码 % 记录面上的等效电流 Js n × H等效磁流 Ms -n × E % surface 是包围目标的记录面n 是外法线方向 % Js_tang, Ms_tang 为记录到的时域切向分量 % 傅里叶变换到频域这里用 DFT 累积方式 freq_points 100; E_far zeros(freq_points, 181); % 每个频点、每个角度的远场 M size(Js_tang, 1); % 记录面网格数 for f_idx 1:freq_points f f_min (f_idx-1) * df; for m 1:M phase exp(1j * 2 * pi * f * (r(m) - r_hat * R / 3e8)); E_far(f_idx, :) E_far(f_idx, :) ... (Js_tang(m, f_idx) * phase) * cos_theta_factor; end end这段代码给了整体框架但实际实现里相位因子的计算很容易出偏差。r(m)是记录面上每个小面元到原点的距离向量r_hat是远场观察方向单位向量两者点乘后的结果要除以光速才能得到时间延迟。很多初学的人把单位弄错算出来的方向图相位是乱的幅度却看不出问题结果被人一问数据就露馅。我的经验是近场记录面至少离目标半个波长太近会有倏逝波污染记录数据。记录面要完全封闭目标不能有开口否则等效源不成立远场结果会缺一块。这两个约束条件在资源包的注释里有写明但很容易被忽略。5. 边界条件与稳定性避坑三种吸收边界的对比和五个血泪排查记录5.1 三种吸收边界对比Mur、PML、CPML边界条件是 FDTD 仿真的最大坑之一。边界设不好反射波会在仿真区域里来回弹污染整个场分布。三种常见方案我尽量把选型边界说清楚。Mur 边界实现最简单二阶 Mur 在一维二维都容易写但斜入射时吸收效果急剧下降适合波近似垂直入射的临时仿真。PML完美匹配层是经典方案引入损耗介质层让入射波在里面衰减掉。实现上需要额外更新 PML 区域内的辅助场变量代码量比 Mur 多几倍吸收效果在宽频和斜入射下都很稳定。CPML卷积完美匹配层是 PML 的改进版用递推卷积替代复杂的分场更新内存占用更低通用性更好是当前主流使用的方案。这份资源包里三种都有实现我个人的用法调试阶段先用 Mur 快速看波形确认源和结构没问题后再切换到 CPML 跑正式数据。Mur 是出问题概率最低的选择一旦 CPML 参数没配对反射曲线会出现奇怪的毛刺但你很难判断是结构本身的响应还是边界引起的伪影。5.2 五个具体踩坑记录坑一CFL 数设置过大导致场值爆炸。现象仿真跑到几百步场值突然变成 inf或者数量级跳到 1e100。原因时间步长超过稳定性极限误差指数增长。解决把 dt 减小到原来的 0.8 倍同时检查 dx 是否太小——有时候空间分辨率提得过高时间步必须跟着降两者是绑定的。坑二边界区域出现高反射场图里能看到明显的驻波条纹。现象波传播到边界后反弹在边界附近形成明暗交替的图案。原因PML 层的厚度不够或者 PML 的损耗系数分布设置得不合理。解决PML 层厚度至少 8 到 12 个网格损耗系数从内到外按多项式渐变通常是二次或三次多项式。我一般会把 PML 层的损耗系数画出来看一眼确保从 0 平滑增大到最大值中间没跳变。坑三介质边界处场值异常跳动。现象两个介质的界面处场值突然出现尖锐突起时间上还在不停震荡。原因材料参数赋值和场分量的 Yee 网格位置不对应——介电常数应该取在两个半格之间的平均值而你直接用了整数点上的值。解决在给 ca、cb 赋材料参数时先对齐 Yee 网格的采样位置。二维里 Ez 是在整数网格点上而 Hx 和 Hy 在半格点所以需要分别对 eps_r 做插值不能共用一套。坑四可视化图像随时间步数增加整体变亮或者场图整体漂移。现象动画播放时背景颜色越来越亮场分布形状没变但数值整体变大。原因直流分量没有消除。高斯脉冲里如果带有零频分量它不会被传播走而是在仿真区域里累积。解决把激励源改成调制高斯脉冲乘上一个正弦载波或者在可视化时显示场值的差分而不是绝对场值。宽带仿真时更推荐调制高斯脉冲中心频率设在实际关注的范围频带宽度按需求调。坑五三维仿真长时间运行后内存耗尽程序崩溃。现象Matlab 弹窗提示内存不足或者系统开始疯狂用虚拟内存导致卡死。原因场分量数组和材料参数数组全部放在内存里三维时尺寸太大。解决把材料参数从 double 转成 single可以省一半内存。再用单精度浮点跑场分量FDTD 精度要求没那么高单精度损失可以接受。实在不够就把仿真区域拆块一次只算一个子区域把边界场值存到硬盘上下次继续算但这属于高级操作了。6. 从代码包到自己的仿真工具CPML 加载和频域响应提取的进阶操作这份资源包真正体现价值的场景是你能把一维里验证过的物理过程快速推广到二维、三维结构里拿到可信的数据。以我过去在项目中用这套代码包的经验最后补一个高频操作把 CPML 正确加载进现有代码并在仿真结束后提取频域响应。CPML 的加载不能直接在主迭代外一两行搞完。它需要在初始化阶段额外分配 psi 数组辅助卷积变量数量和场分量一样多每个场分量的吸收边界更新公式都要在迭代循环里单独加一段。资源包里 CPML 代码的更新顺序是这样先用普通 FDTD 公式更新场值再用 CPML 项修正。修正项本质上是对场梯度做一次有损耗的递推滤波滤波系数由位置决定——接近外层边界时吸收越强。% CPML 电场修正段示意 % psi_Ez_x 和 psi_Ez_y 为 CPML 辅助变量 % kappa_x, alpha_x, a_x, b_x 为 CPML 位置相关参数 for i 2:Nx for j 2:Ny % 修正 Hy 差分项 psi_Ez_x(i,j) b_x(i) * psi_Ez_x(i,j) a_x(i) * (Hy(i,j) - Hy(i-1,j)); psi_Ez_y(i,j) b_y(j) * psi_Ez_y(i,j) a_y(j) * (Hx(i,j-1) - Hx(i,j)); Ez(i,j) Ez(i,j) cb(i,j) * (psi_Ez_x(i,j) / kappa_x(i) psi_Ez_y(i,j) / kappa_y(j)); end end这里 a_x 和 b_x 的推导依赖稳定参数 sigma、alpha、kappa代码注释里给了推荐初值但用的时候要检查一件事PML 区域和外层边界的间距必须大于等于 3 个网格否则边界上的离散误差会把 CPML 的吸收效果削弱。这是我在一次微带线仿真里验证出来的——把 PML 贴着外层边界设反射系数比正常情况下高了一个量级。仿真结束后提取频域响应通用做法是在时域记录某个探针点的场值序列然后做快速傅里叶变换得到频谱再和入射场的频谱相除得到传输系数。如果做的是宽频带仿真频域分辨率由时间步数和时间步长决定频域分辨率 1 / (Nt * dt)。要分辨 0.1 GHz 的精细结构就需要至少 10ns 的总仿真时长换算成步数就是几千到几万步不等。这个换算关系在项目里提前算好不然跑完才发现频率分辨率不够又要加大步数重跑一等就是几小时。另外提一句这份代码是我在 A 同学的项目文档里看到的当时他做的某种特定结构中 CPML 参数直接拿初值用就出了问题换了入射角度之后反射误差大得没法看。调整之后把 sigma 的厚度分布改成非对称的吸收效果才稳定下来。从那以后我每次搭建新几何模型都强制走一遍固定流程先用一维验证时间步长和激励源再在二维里检查 CPML 反射水平最后迭代三维模型。哪怕只是改了一个介质块的位置也把每一维度都过一遍。这个习惯看着费时间但至少能确保最后跑出来的远场数据拿得出手。希望帮到你动手跑一遍比看十遍文章都有用。本文还有配套的精品资源点击获取