一维到三维FDTD：MATLAB电磁仿真代码、避坑与加速实战

简介这是一套基于时域有限差分FDTD方法的电磁场数值模拟程序集使用矩阵实验室Matlab语言编写面向电磁场仿真入门者、光电通信及微波工程相关专业的学生与工程师能够满足一维、二维、三维电磁波传播问题的仿真需求。压缩包内共十份文件以脚本为主包含六个脚本程序分别实现一维、二维、三维的差分计算另有备份文件、图形界面与说明文档整体大小仅八十七千字节87KB轻量便携适合快速运行和二次开发。已有七百零三人下载学习热度较高。程序不仅包含核心计算逻辑还配有动态图形呈现运行后可直接观察电磁场在空间中的演化过程有助于直观理解吸收边界条件、数值色散等关键概念同时学习矩阵实验室的数值编程技巧对入门和进阶均有实用价值。1. FDTD从一维到三维这份MATLAB代码包到底在帮你解决什么问题电磁仿真做到一定阶段几乎绕不开 FDTD时域有限差分法。它的思路很朴素把空间切成 Yee 网格把时间切成步长每一步显式更新电场和磁场直接模拟电磁波在结构里的传播、反射、散射和透射。围绕这个标题的 MATLAB 代码包就是把一维、二维、三维 FDTD 的实现放在一起方便按需取用。对做光子晶体、天线辐射、电磁兼容或超表面的人来说这类代码最大的价值是让你在理论推导之后第一时间亲眼看到场分布和频谱结果不用盲目依赖商用软件的黑匣子。适合两类人一类是刚入门的研究生想从一维开始把 FDTD 的数值逻辑亲手跑通另一类是手上有矩形波导、柱状散射体这类规则结构需要快速算出反射率和透射率的工程师。下面我按一维、二维、三维的顺序把代码结构、参数选择和常见坑一次讲透。2. 一维FDTD从麦克斯韦旋度方程到第一段能跑的MATLAB代码2.1 一维FDTD为什么是入门最快的路径一维是理解 FDTD 的最小系统。假设电磁波沿 x 方向传播电场只有 Ez 分量磁场只有 Hy 分量麦克斯韦旋度方程在这个限制下只剩两个标量方程∂Ez/∂t (1/ε)(∂Hy/∂x)∂Hy/∂t (1/μ)(∂Ez/∂x)离散时让 Ez 和 Hy 在空间上错开半个网格、在时间上错开半个步长这就是 Yee 网格加蛙跳格式。这么做的直接好处是显式推进每一步先用上一时刻的 Ez 算新的 Hy再用新的 Hy 算新的 Ez全程不需要解线性方程组迭代几百步都不会累积矩阵误差。代价是时间步长必须满足 CFL 条件否则数值误差按指数增长很快爆掉。一维的 CFL 条件是 dt ≤ dx/c。我一般取 dt dx/(2c)哪怕后面改了网格参数这个裕量也能兜住大部分手误。从物理上看CFL 条件表达的是因果性每个时间步内信息最多只能跨过一个网格。离散格式如果不尊重这个约束等于是让波在一步之内跳过多个网格数值色散会放大到不可控。所以你会发现 FDTD 的调参顺序永远是先定 dx再定 dt顺序不能颠倒。dx 决定空间分辨率dt 决定时间步数两者一起决定总计算量。2.2 一维FDTD最小代码高斯脉冲在自由空间传播% 一维FDTD高斯脉冲在自由空间中传播 clear; clc; % 空间与网格参数 nx 600; % 空间网格数 dx 1e-3; % 空间步长dx 1 mm c 3e8; % 真空光速 eps0 8.854e-12; % 真空介电常数 mu0 pi*4e-7; % 真空磁导率 dt dx / (2*c); % 时间步长满足CFL条件并留裕量 nt 800; % 总时间步数 Ez zeros(1, nx); % 电场分量z方向 Hy zeros(1, nx); % 磁场分量y方向 % 高斯脉冲参数 t0 40; % 脉冲峰值对应时间步 spread 12; % 脉冲半宽度时间步 for n 1:nt % 先更新磁场Hy 由相邻 Ez 的差分决定 for i 1:nx-1 Hy(i) Hy(i) - (dt/(mu0*dx)) * (Ez(i1) - Ez(i)); end % 再更新电场Ez 由相邻 Hy 的差分决定 for i 2:nx Ez(i) Ez(i) - (dt/(eps0*dx)) * (Hy(i) - Hy(i-1)); end % 在网格第20个点注入高斯脉冲源 Ez(20) Ez(20) exp(-((n - t0)/spread)^2); end % 画最终波形 plot(1:nx, Ez, b-, LineWidth, 1.2); xlabel(网格坐标); ylabel(Ez (V/m));这段代码是硬源注入脉冲会同时向左右两个方向传播这是正常现象。更新顺序是先 Hy 后 Ez对应蛙跳的时间关联n 时刻的 Hy 与 n-1/2 时刻的 Ez 同层n 时刻的 Ez 与 n-1/2 时刻的 Hy 同层。索引上Hy(i) 看 Ez(i) 和 Ez(i1)Ez(i) 看 Hy(i) 和 Hy(i-1)正好是半个网格的错位照抄即可。两个更新公式里的系数 dt/(mu0dx) 和 dt/(eps0dx) 是缩放因子它们的准确性直接影响波速所以 eps0 和 mu0 必须用真实值不能随手填 1。参数上dx 决定最高可分辨频率。波动性问题里 dx 一般不超过最小波长的 1/20这里 dx1mm 对应最高有效频率约 150GHz演示足够。t0 和 spread 控制脉冲形状spread 越大频谱越窄波包在空间里拉得越长模拟周期结构时 spread 通常要设到目标波周期的若干倍。想看某个位置随时间变化的信号就在循环里加一个探针数组把 Ez(probeIdx) 存下来probeIdx 300; % 探针位置 probeEz zeros(1, nt); % 预分配探针记录 for n 1:nt % ... 前面两段更新代码 ... probeEz(n) Ez(probeIdx); end plot(1:nt, probeEz); % 画时域波形这里预分配 probeEz 是关键习惯MATLAB 里在循环内不断扩数组会让性能断崖式下降。探针法也是后面二维、三维仿真里做时域监视的标准手段一维阶段养成习惯会省很多事。2.3 一维边界为什么是第一个要补的洞上面代码跑完你会看到波形到达数组两端后反弹回来。这不是 bug是默认的 PEC 边界在起作用。做反射谱或透射谱时反弹波会污染结果所以边界处理是 FDTD 里绕不过去的一关。最简单的补法是加一阶 Mur 吸收边界在左右两个端点分别用前一时刻的场来近似外向波% 左边界i1的一阶Mur吸收 Ez(1) Ez(2) (c*dt - dx)/(c*dt dx) * (Ez_old(2) - Ez(1)); % 右边界inx的一阶Mur吸收 Ez(nx) Ez(nx-1) (c*dt - dx)/(c*dt dx) * (Ez_old(nx-1) - Ez(nx)); % 注意每次更新前先把 Ez_old Ez 存起来这段代码里 Ez_old 是进入主循环前保存的上一时刻电场Mur 边界本质上是用一维外向波方程做外推实现简单但吸收效果有限。更实用的做法是从一开始就上 CPML一维 CPML 只需要多维护两个辅助变量代码量只多十几行但能把边界反射压到 -40dB 以下。我自己的习惯是一维阶段就养成源、更新、边界三段式结构后面升到二维、三维只改更新段和边界段的维度源与输出的骨架不动这个习惯帮我省掉了大量排错时间。3. 二维FDTDTE/TM模怎么选点源辐射代码与网格参数搭配3.1 先定TE还是TM极化方向决定用哪组场分量二维 FDTD 第一个选择是 TE 还是 TM 模。TM 模的场分量为 Ez、Hx、Hy电场只有 z 分量适合模拟电场垂直入射平面的情况比如无限长介质柱的散射、光子晶体波导的横截面TE 模的分量为 Hz、Ex、Ey电场在平面内适合磁场垂直入射平面的场景。两组公式完全对称差别只是把电场更新换成磁场更新、把介电常数和磁导率对调。选择依据很简单看激励源极化方向是垂直还是平行于结构轴向。绝大多数光电器件仿真用 TM 模更多因为 z 方向的电场用一个标量 Ez 就能表示代码省一半。如果结构里有金属选择更要谨慎对 PEC 边界TE 模需要强制电场切向分量为零TM 模则要处理磁场法向分量换模时边界条件要跟着改。3.2 二维TM模点源辐射完整的可运行代码% 二维FDTDTM模Ez/Hx/Hy点源辐射演示 clear; clc; nx 201; ny 201; % 网格数量 dx 1e-3; dy 1e-3; % 空间步长 c 3e8; % 光速 eps0 8.854e-12; mu0 pi*4e-7; % 二维CFL条件取90%安全系数 dt 0.9 / (c * sqrt(1/dx^2 1/dy^2)); nt 300; Ez zeros(nx, ny); % 电场 z 分量 Hx zeros(nx, ny); % 磁场 x 分量 Hy zeros(nx, ny); % 磁场 y 分量 t0 30; spread 8; % 脉冲参数 for n 1:nt % Hx 更新依赖 y 方向的 Ez 差分 for j 1:ny-1 for i 1:nx Hx(i,j) Hx(i,j) - (dt/(mu0*dy)) * (Ez(i,j1) - Ez(i,j)); end end % Hy 更新依赖 x 方向的 Ez 差分 for j 1:ny for i 1:nx-1 Hy(i,j) Hy(i,j) (dt/(mu0*dx)) * (Ez(i1,j) - Ez(i,j)); end end % Ez 更新由 Hx 和 Hy 的差商共同决定 for j 2:ny-1 for i 2:nx-1 Ez(i,j) Ez(i,j) (dt/eps0) * ... ((Hy(i,j) - Hy(i-1,j))/dx - (Hx(i,j) - Hx(i,j-1))/dy); end end % 中心点源注入硬源 Ez(101,101) Ez(101,101) exp(-((n - t0)/spread)^2); % 每40步可视化一次 if mod(n, 40) 0 imagesc(Ez); axis equal; colorbar; drawnow; end end这段代码对应 TM 模核心是三个更新循环顺序不能乱。先算 Hx 和 Hy再算 Ez这保证时间上保持半格交错。二维 CFL 公式比一维紧dt 取了 90% 的安全系数原因是斜向传播时波的等效相速度更高CFL 约束更严格。imagesc(Ez) 里的转置容易被忽略imagesc 把每一列画成一行转置后 x 轴对应 i、y 轴对应 j坐标方向才和物理空间一致不转置会让你看到电场分布像被翻转了 90 度。参数上nx 和 ny 取奇数好用点源放在 (101,101) 正好在中心。dx 和 dy 可以不等但 CFL 必须同时用两个方向步长参与计算。spread8 的高斯脉冲在空间上大约覆盖 16 个网格如果这个值太小比如小于 3频谱太宽网格分辨率跟不上你会看到明显的各向异性伪影波形呈方形扩散而不是圆形。检查这种伪影最好的办法是看可视化结果各向同性介质中点源辐射应该是标准圆形波前一旦出现菱形或方形轮廓不是 dx 太大就是 CFL 系数给得太满。3.3 网格分辨率怎么取dx不是越小越好dx 的选择要同时权衡计算量、数值色散和结构精度。经验法则是 dx ≤ λ_min / 20λ_min 是模拟频段内的最短波长。可见光波段若用中心波长 550nmdx 取 25nm 以下微波 3GHzdx 取 5mm 以下即可。数值色散是更隐蔽的约束当 dx 相对波长不够小时不同传播方向的波会有不同的相速度导致波形在传播中畸变脉冲越走越宽甚至出现拖尾。仿真场景推荐dxdt参考c3e8典型网格量二维微波 3GHz5 mm约 8 ps很小几十万量级太赫兹 1THz15 μm约 25 fs中等百万量级可见光 550nm25 nm约 40 as大千万量级注意 dx 减半二维网格量和内存涨 4 倍计算时长也约涨 4 倍所以不要无脑加密网格。斜边或圆弧结构用矩形网格拟合时还有阶梯近似误差常见做法是对边界网格的介电常数做加权平均或者用共形 FDTD。至少心里要有数斜边界结构的仿真结果会随 dx 变化收敛速度比直边界慢得多。二维 FDTD 改成 TE 模也不复杂把 Ez/Hx/Hy 换成 Hz/Ex/Ey更新公式里 eps 和 mu 互换源类型从电偶极子换成磁偶极子代码改动量大概 20 行。4. 三维FDTD内存墙与计算量才是真正的门槛4.1 三维FDTD的六分量更新规则三维 FDTD 不再有极化方向的简化电场和磁场各有三个分量Ex、Ey、Ez、Hx、Hy、Hz总共 6 个三维数组。每个场分量由另外两个磁场或电场分量的空间差分驱动。比如 Ex 由 Hz 在 y 方向上的差分和 Hy 在 z 方向上的差分共同驱动写法是Ex 的更新量 (dt/ε) * ( ∂Hz/∂y - ∂Hy/∂z )Ey 由 Hx/z 和 Hz/x 驱动Ez 由 Hy/x 和 Hx/y 驱动磁场更新则是完全对称的对偶关系。在 MATLAB 里组织这个问题最常用的方式是维护 6 个三维数组主循环里写 6 段更新代码每段都是三重循环。二维里先磁场后电场的顺序在三维里变成先三个磁场分量再三个电场分量同样保持蛙跳格式的半时间步交错。三维的索引规则比二维更复杂我建议在代码旁注释每个偏导对应的索引差否则调试时非常容易看晕。4.2 跑一个200×200×200的网格要吃多少内存这是三维 FDTD 最现实的问题。6 个 double 型数组每个尺寸 Nx×Ny×Nz总内存就是 6 × Nx × Ny × Nz × 8 字节。我列出几个典型规模你感受一下增长速度网格规模6场分量总内存约实际可行性100×100×10048 MB轻松200×200×200384 MB普通台式机可以300×300×3001.29 GB需要留意建议单精度500×500×3003.6 GB纯MATLAB基本跑不动从这个表能看出三维比二维不是多一个维度而是多一到两个数量级的内存和计算量。很多人第一次在三维 FDTD 翻车不是公式写错而是内存直接耗尽或者每个时间步慢到无法忍受。常见做法是先把网格降到单精度能省一半内存或者干脆把研究问题压回二维只有当结构在 z 方向确实有变化时才上三维。若真要跑大网格建议把核心更新代码单独写成 MEX 或 C 函数MATLAB 只负责调度和可视化。标题里提到的 fdtd 方向通常指这种面向对象封装的 C FDTD 实现原理与 MATLAB 版完全一致只是把数组操作换成指针和高性能循环。4.3 从二维升维三维改造代码的三个关键步骤从二维 TM 升三维不是简单加一个 z 方向循环就完事。我一般按三步走。第一步把二维数组改成三维数组给所有场分量建立从 1 到 N 的完整索引zeros 预分配一个都不能少。第二步把更新公式里每个偏导数写成对应方向的差分比值这一步最容易错建议先在纸上把 Ex、Ey、Ez 对应的两个偏导项写全再对照 Hz、Hy 的实际索引导入代码。第三步处理 z 方向上的网格边界原来的 PEC 边界改到 z 方向的两个端面侧向四边保留吸收边界。% 三维FDTD电场更新示意只列Ex与Hz、Hy的耦合 % 数组尺寸Ex(nx,ny,nz)其余分量同 for k 2:nz-1 for j 2:ny-1 for i 1:nx-1 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 ); end end end % 注意Hz的y方向差分和Hy的z方向差分对应Ex的两个驱动项 % 磁场更新与之对称顺序上先算Hz、Hy再算Ex这个片段里最容易错的地方是索引Hz 的分母是 dy因为它是 y 方向差分Hy 的分母是 dz因为它是 z 方向差分。一旦搞反仿真结果会非常怪异。常见错误是把二维代码里的j2:ny-1直接抄成k2:nz-1却忘掉 y 方向的更新仍需要 j 循环还有的代码漏掉某个磁场分量的更新导致能量不守恒跑几百步后场分布明显畸变。验证三维代码是否写对的快速手段是平面波正入射一个已知厚度的介质平板与传输矩阵计算结果对比只要在 z 方向取一条线上的场做 FFT看透射系数的幅度和相位是否一致这个测试能过滤掉绝大部分索引错误。5. FDTD避坑指南数值爆炸、边界反射与结果对不上的排查记录5.1 数值爆炸电场值跑到1e50现象仿真跑到一两百步场值突然变成 1e50 甚至 NaN波形完全失控。原因最常见的是时间步长 dt 超过 CFL 上限。FDTD 显式格式的稳定性完全依赖 CFL 条件违反它误差按指数增长。另一个常见原因是介质参数有问题介电常数被设成了 0 或负值或者 PML 参数里电导率为负。第三个原因是源注入时把硬源赋值加在了已经被边界更新过的点上当场破坏物理场。解决先按 dt 0.5 / (c * sqrt(sum(1./dx.^2))) 重取时间步长把所有系数打到原来的一半再跑如果立刻稳定说明是 CFL 问题再打印 eps 和 mu 分布检查有没有非物理数值最后确认源的位置没有落在更新区域的边界网格上。这个排查顺序我用了很多次90% 的爆炸都是前两个原因。5.2 边界反射透射谱上出现周期波纹现象计算透射谱时频谱上叠加周期性波纹波纹间隔与仿真区域尺寸相关。用一维代码观察更直观高斯脉冲到达末端后反弹产生一串等间距的小波包。原因没加吸收边界端面默认 PEC电磁波全反射回来。PML 厚度不够比如小于 8 层或者电导率渐变过猛也会造成寄生反射。还有一种隐蔽情况是 PML 与结构贴得太近近场分量还没展开就撞上 PML反射量显著上升。解决至少用一阶 Mur 或 PML 把边界包住对光谱类仿真CPML 厚度取 1015 层电导率用多项式渐变并在结构边缘与 PML 之间留出 46 个网格的缓冲区。判断边界是否合格的方法是测反射系数在边界前放探针比较反射波峰值与入射波峰值的比值通常要求 -40dB 以下。5.3 透射率峰位对不上解析解或实验值现象计算的谐振波长偏了几百纳米或者透射率曲线整体左移右移怎么调参数都差一截。原因单位制混乱占大头。有人用纳米定义几何但 dx 按米写成了 1e-9结构尺寸却用微秒数直接乘 dx结果尺寸错了几百倍。另一个常见错误是把折射率当相对介电常数用把 1.5 直接代入 eps导致波速算错所有谐振峰集体偏移。解决在代码开头统一单位定义 scale 1e-6 表示微米dx 0.05*scale所有几何尺寸都用格点数 × dx表达。折射率 n1.5 对应 eps_r 2.25μ_r 通常取 1。跑完结构后先用解析平板透射率做一次校零确认单位换算无误再进入参数扫描。血泪经验单位问题排查起来最耗时间因为它不会报错只会在结果上悄悄偏一点点。5.4 仿真速度慢三维网格跑不动现象200×200×200 网格每个时间步循环要十几秒一次仿真跑到需要大半天的量级。原因MATLAB 的嵌套 for 循环解释执行开销巨大没有预分配数组会导致动态扩容每步都在重新分配内存两者叠加性能可以差几十倍。还有一个隐蔽原因是可视化代码加在了主循环里每步都执行 imagesc把 GPU 或 CPU 数据不断拷贝刷新严重拖慢进度。解决先给所有场数组做 zeros 预分配把主循环里的可视化改为每 40 步或每 100 步一次再把三重循环改成向量化的矩阵切片运算。改完还慢就把核心数组用 gpuArray 搬到显卡上配合 gather 取结果。向量化之前直接上 GPU 没有意义嵌套循环在 GPU 上依旧慢。5.5 能量缓慢增长格式稳定但物理不对现象监视全部网格的总能量即 sum(epsE.^2/2 muH.^2/2)发现能量随时间缓慢上升甚至越涨越快。原因最常见的是边界没处理好反弹波回到源区与入射波相干叠加形成驻波能量被持续泵进仿真区域。另一个常见原因是硬源本身在持续注入能量高斯脉冲的时间窗已经过了但源还在周期性赋值等价于往系统里不断加功率。解决把脉冲源设置成一次性用时间窗包住或者改用软源把源写成电流密度项加在更新公式里而不是直接覆盖场值。软源的形式是 Ez(i) Ez(i) exp(...)它通过方程注入能量比硬源物理得多。最后再检查边界 PML 是否真正吸收了尾部能量这个检查做一次就能发现前面所有问题。6. 进阶CPML吸收边界与GPU加速让FDTD从能跑到跑得快6.1 CPML边界参数与使用要点CPML卷积完全匹配层是目前最常用的吸波边界比手写 Mur 一阶吸收边界干净得多。参数经验值层数取 1015 层电导率从内到外按多项式渐变p 取 34最大电导率约取 sigma_max 0.8*(p1)/(150pidx)。CPML 的辅助变量要单独开数组数量与场分量相同内存占比不小用单精度存储即可。比较容易被忽略的是缓冲区结构边缘与 PML 之间留 46 个网格让波在到达 PML 前基本成为行波吸收效果会明显更好。6.2 GPU加速三行代码让二维代码提速数倍% 把核心场数组搬上GPU Ez gpuArray(Ez); Hx gpuArray(Hx); Hy gpuArray(Hy); % 主循环照常运行MATLAB自动在GPU上执行向量化操作 % 结果取回CPU用gather Ez_cpu gather(Ez);gpuArray 不是魔法嵌套 for 循环搬到 GPU 上依旧慢必须先做向量化或者至少把最内层循环改成 arrayfunGPU 才有明显收益。MATLAB 在 GPU 上执行时会自动管理显存但要注意显存容量200×200×200 的双精度 6 分量数组已经接近 6GB家用显卡很可能放不下。实践上我通常把二维和三维的数组都转成 single 再给 GPU速度更快内存压力也更小。向量化加单精度这一套组合下来二维点源辐射代码提速 510 倍是很正常的。6.3 一个值得养成的验证习惯每次写新结构前先跑一个已知解的最小模型。一维平板透射、二维圆柱散射、三维球体 Mie 散射都可以算好后与解析解对比误差在几个百分点内再继续下一步。这个习惯成本很低但能把代码 bug和物理参数错误清晰分开省下的调试时间非常可观。我现在跑新结构第一步永远是写一个五分钟能跑完的最小模型验证边界和源没有问题这个习惯救了我很多次。希望帮到你。本文还有配套的精品资源点击获取