简介面向电磁仿真学习者与科研人员的时域有限差分FDTDMATLAB实现资源包覆盖天线设计、微波工程、光子学等典型场景可帮助读者从Yee网格离散、PML边界条件到麦克斯韦方程迭代更新快速上手。压缩包共426个文件以393个m脚本为主体另有br备份、txt说明、asv存档以及少量dll、exe、fig配套文件整体大小37.23MB目录分层清晰便于对照学习与二次开发。已有521人学习下载。资源包内不仅给出完整二维FDTD仿真流程还包含源项设置、电磁场采样、输出参数定义等模块可直接运行观察场分布也可根据实际需求调整网格尺寸、时间步长和材料参数扩展至不同电磁仿真项目。对于希望系统理解FDTD算法并用MATLAB完成仿真验证的初学者和工程师是一份兼具教学与实用价值的参考资料。1. 一份时域有限差分MATLAB代码包跑通和跑对是两回事A同学从某网盘拿了一份「时域有限差分的MATLAB代码.zip」解压后主程序只有一百来行双击运行波形确实出来了可拿去向文献结果一对谐振频率偏了百分之八边界上还有一股反射波来回弹。这个现象在计算电磁学里太典型了时域有限差分FDTD把麦克斯韦方程组直接在时域离散迭代算法骨架简单核心循环写出来不到二十行但真正决定结果对错的是Yee网格怎么排、稳定条件取多少、吸收边界和激励源参数怎么配合。这篇文章就按一份能用的MATLAB实现往下拆——先立住最小可运行的代码再把参数和坑逐个讲透让新手拿这份代码能跑出可信结果让熟手能快速定位自己那份代码为什么翻车。2. FDTD的核心离散逻辑与MATLAB的适用边界先搞清楚解法再谈代码2.1 Yee网格和蛙跳迭代到底在算什么FDTD的出发点是把时间和空间都切成离散网格然后用中心差分去近似麦克斯韦方程里的旋度运算。国外在1966年提出这个思路时最关键的一步是把电场和磁场在空间上交错半个网格步长时间上也相差半步。这样做的直接好处是电场分量和磁场分量的更新各自只需要周围一圈邻居的值不需要解联立方程组每一时间步都是显式推进。以二维TM波为例场分量只有Ez、Hx、Hy三个。计算区域被划分成矩形网格Ez定义在整数格点(i,j)Hx定义在(i,j1/2)y方向错半格Hy定义在(i1/2,j)x方向错半格。时间上电场在第nΔt时刻取值磁场在第(n1/2)Δt时刻取值两者交替更新。一个完整时间步分成两半先用当前电场值算出磁场的变化量再用更新后的磁场算出电场的变化量如此循环这就是“蛙跳”迭代的由来。这种交错排布不是巧合它让空间差分达到二阶精度。如果用同样的网格去写电场和磁场中心差分就会退化成向前或向后差分数值色散会大得多结果错得离谱。看代码时如果发现Hx、Hy数组尺寸比Ez小一圈或者切片偏移一位不用惊讶那正是在体现这个半网格偏移。2.2 教学验证和器件级预研场景下MATLAB比C和Python好在哪我见过不少人在网上争论FDTD到底该用什么语言写结论往往吵成一团。实际上这个问题的答案高度依赖场景。一份教学用途的FDTD代码核心价值不是跑得快而是让读代码的人能在二十分钟内理解每一步在做什么。MATLAB在这点上优势非常明显数组切片直接对应场更新公式不写循环也能完成整个区域的更新向量化运算让二维区域的单步更新几乎就是几行矩阵运算内置的绘图和FFT让时域波形、场分布图、频谱结果随手就能画出来。C语言在FDTD里的优势是性能但代价是内存管理和指针操作占掉大量精力。Python科学栈接近MATLAB的开发效率但在多层循环场景下如果不做矢量化和Numba优化性能还不如MATLAB的向量化。表格对比更直观维度MATLABC/CPython(纯NumPy)代码可读性高接近公式低需管理索引中需熟练NumPy切片二维场更新性能高向量化高中视写法而定三维大算例慢内存压力大快可控慢调试体验好变量随时可视化一般需工具链好典型用途教学、验证、预研生产级仿真教学、数据处理判断准则很简单如果网格规模在百万量级以下、要跑的是二维问题或小尺寸三维问题、重点是验证算法和理解物理过程那MATLAB是性价比最高的选择。如果目标是一套能支撑工程设计迭代的仿真平台网格动辄几千万甚至上亿那无论如何都要转向编译型语言。这份代码包属于前者所以下文所有讨论都建立在MATLAB向量化写法上。2.3 稳定性和数值色散CFL条件取多少网格多密才够FDTD的显式迭代有一个硬性稳定条件叫CFL条件。它说的是时间步长不能超过电磁波在一个空间步长内传播所需的时间否则数值解会指数增长直至发散。二维均匀网格下稳定上限是dt_max 1 / (c * sqrt(1/dx² 1/dy²))当dx dy时dt_max dx / (c * sqrt(2))。实际取值时我习惯乘一个安全系数Courant数取0.9左右也就是dt 0.9 * dt_max。这个系数不是越大越好它只保证稳定不保证精度也不是越小越好因为时间步太小意味着达到相同物理时间需要更多迭代步纯粹浪费时间。另一个容易被忽略的量是每波长网格数英文缩写PPWPoints Per Wavelength。假设波在介质中的波长为λ网格尺寸dx应满足dx ≤ λ/10经验上取λ/12到λ/20能得到比较可靠的结果。如果网格太粗数值色散会让波以错误的相速度传播表现就是谐振频率偏移、波形拖尾、脉冲展宽。这个误差是FDTD方法固有的加密网格可以减小但无法完全消除验证代码正确性时网格收敛测试是必须做的一步后面第6章会专门讲。3. 跑通2D TM波最小实现网格参数、主循环与软源注入3.1 一份典型FDTD代码包解压后通常能看到哪些文件拿到这种zip解压后第一件事不是找主脚本而是先看目录结构。一份组织良好的FDTD教学代码一般包含四类内容一个主脚本负责参数设定和调用流程若干个函数文件比如场更新函数、边界条件函数、激励源函数一个或几个后处理脚本负责画场图、画探针波形、做FFT可能还有一个README或注释文档说明运行方式。如果解压出来的东西只有一个孤零零的m文件全部逻辑堆在一起那要格外留意。不是说不可以而是说明作者没有把参数、核心循环、后处理分开这类代码改起来牵一发动全身。我拿到代码后通常先搜三个关键变量的赋值dx/dy、dt或courant、介质数组eps_r或epsilon。这三个值决定了网格分辨率和材料分布是后续所有结果的基础先把它们看懂再去看主循环。3.2 参数初始化网格尺寸、时间步长、材料分布怎么定下面这份是最小可运行的二维TM波主程序没有用PML第4章会补也没有介质损耗目的是先验证迭代算法本身。把代码完整复制进MATLAB编辑器直接运行就能看到探针处的波形。% fdtd2d_tm_demo.m % 二维 TM 波 FDTD 最小实现Ez、Hx、Hy 三分量 % 使用软源激励观测点记录时域波形 c0 3e8; mu0 4*pi*1e-7; eps0 1 / (c0^2 * mu0); % ---- 网格参数 ---- nx 120; % x 方向网格数 ny 100; % y 方向网格数 dx 1e-3; % 空间步长单位米 dy 1e-3; courant 0.9; % CFL 安全系数 % 二维均匀网格的 CFL 上限乘系数得到实际时间步 dt courant / (c0 * sqrt(1/dx^2 1/dy^2)); nt 800; % 迭代步数 % ---- 材料分布 ---- eps_r ones(nx, ny); % 相对介电常数默认空气 [xg, yg] meshgrid((1:nx)*dx, (1:ny)*dy); idx (xg - 70*dx).^2 (yg - 50*dy).^2 (20*dx)^2; % 半径20格的介质圆 eps_r(idx) 4; % 介质圆内介电常数设为4 % ---- 场数组Yee半网格偏移体现在数组尺寸上 ---- Ez zeros(nx, ny); Hx zeros(nx, ny-1); % Hx 定义在 (i, j1/2) Hy zeros(nx-1, ny); % Hy 定义在 (i1/2, j) % ---- 激励源软源放在 (20, 50)高斯脉冲 ---- sx 20; sy 50; tau 100 * dt; % 脉冲时间常数约0.21ns t0 4 * tau; % 峰值时刻 Probe zeros(1, nt); % 观测点 (100, 50) 的电场记录 for n 1:nt % 1) 更新 HxEz 对 y 的差分注意负号 Hx(:, :) Hx(:, :) - (dt/(mu0*dy)) * (Ez(:, 2:ny) - Ez(:, 1:ny-1)); % 2) 更新 HyEz 对 x 的差分 Hy(:, :) Hy(:, :) (dt/(mu0*dx)) * (Ez(2:nx, :) - Ez(1:nx-1, :)); % 3) 更新 Ez由 Hy、Hx 的散度项介质区除以 eps_r Ez(2:nx-1, 2:ny-1) Ez(2:nx-1, 2:ny-1) (dt/eps0) * ... ( (Hy(2:nx-1, 2:ny-1) - Hy(1:nx-2, 2:ny-1))/dx ... -(Hx(2:nx-1, 2:ny-1) - Hx(2:nx-1, 1:ny-2))/dy ) ... ./ eps_r(2:nx-1, 2:ny-1); % 4) 软源注入累加到源点场值上不是覆盖 t n * dt; Ez(sx, sy) Ez(sx, sy) exp(-((t - t0)/tau)^2); % 5) 记录观测点 (100, 50) 的电场 Probe(n) Ez(100, 50); end % 画时域波形 plot((1:nt)*dt*1e9, Probe, LineWidth, 1.2); xlabel(时间 (ns)); ylabel(Ez (V/m)); title(观测点电场时域波形未加吸收边界);3.3 主循环里的每一步分别对应什么物理量代码以时间步为单位循环每个循环内做五件事。第一步更新Hx用的是Ez在y方向的差分乘上一个系数dt/(mu0*dy)。这里的负号来自旋度方程顺序不能反反了等效于时间倒流。第二步更新Hy正号是从方程里自然推出来的。数组切片(Ez(:,2:ny)-Ez(:,1:ny-1))沿y方向做前向差分结果正好落在Hx所在的错半格位置上这是Yee网格设计的核心。第三步更新Ez这是整个循环里最容易被写错的一行。它把Hy在x方向的差分减去Hx在y方向的差分再除以介电常数。注意./ eps_r这里用的是逐元素除法因为介质圆内外的介电常数不一样真空区和介质区必须分别处理。计算区域最外一圈的Ez没更新是因为边界上没有足够的邻居值可用这在没有吸收边界的版本里会自然形成一个反射面暂时不管它。第四步是软源注入。软源和硬源的区别在于它是把脉冲值加到源点原有场值上而不是直接把源点场值覆盖成脉冲值。硬源会阻挡入射场穿过源点等效于在源点放了一个金属边界会产生没来由的反向散射。很多新手跑出奇怪波形第一反应怀疑边界条件其实问题出在硬源上。换成软源后源点既发射波形又允许从别处返回的场自由穿过物理上更合理。3.4 这段代码的预期行为和改进空间源点在(20,50)观测点在(100,50)间距80个网格。光在真空走完80格大约0.27ns折合约130个时间步脉冲峰值在t0400步附近发射所以观测点在530步左右能收到第一个波峰这段代码跑完800步正好能看到一个完整脉冲经过。因为还没加吸收边界波形尾部会跟着一些边界反射回来的杂波幅度不大但确实存在。想验证这段代码是否正确可以把eps_r全设为1再跑一遍此时观测点波形应该和源点波形形状一致只有时延和幅度衰减。如果波形出现畸变或数值振荡说明索引切片写错了用这个小测试能很快定位。跑通这段之后下一步就是把它升级成带PML吸收边界的版本这是第4章的活。4. PML吸收边界与激励源参数结果可信度的两个开关4.1 PML层数和电导率剖面怎么配反射率才够低第3章的代码用零边界等效波碰到计算区域边缘会原路反射回来这会污染观测点的时域波形。解决方法是PML完全匹配层——在计算区域外圈加一层有耗介质让入射波在层内指数衰减到最外层时幅度已经小到可以忽略。理想PML理论上对所有入射角度和频率都无反射但离散实现里只能做到接近。PML的参数有三个直接决定吸收效果层数NPML、电导率最大值σ_max、电导率剖面的多项式阶数m。σ_max太小波衰减不够太大介质参数突变本身会产生反射。工程上常用多项式渐变剖面σ从内边界向外边界按(m1)次幂增长。σ_max的经验公式为σ_max -(m1) * ln(R0) / (2 * η0 * NPML * dx)其中η0是自由空间波阻抗约377欧姆R0是理论反射系数通常取1e-4到1e-2。实际调试时先固定R0再看NPML取多少合适。我的常用配置是NPML10m3R01e-4。层数再多效果提升有限但计算区域会明显变大不值得。参数推荐范围说明NPML812层少于8层低频吸收不足多于12层收益递减m243最常见m越大内部渐变越平缓R01e-41e-2越小σ_max越大但过大会引起界面反射σ_max按公式算不要拍脑袋按上面公式计算后微调MATLAB里生成σ剖面很简单。沿x方向的一维剖面写成sigma_xy方向同理在更新方程里把哈达玛乘法逐元素乘积加进去。核心思路是在PML区域内原场更新结果每步额外乘一个衰减因子因子大小正比于该位置的σ值。调试时检查PML是否生效的办法是跑一个一维测试把观测点放在PML内边界前面看反射波是否比零边界版本低两个数量级以上。4.2 高斯脉冲宽度与频谱覆盖想要多宽频带脉宽就有多窄FDTD的优势是一次时域计算得到宽频带响应但宽频带不是白来的它取决于激励源本身的频谱范围。最常用的高斯脉冲形式为E(t) exp(-((t - t0) / τ)²)它的频谱近似高斯分布有效频带上限约f_max ≈ 0.5 / τ-3dB带宽约0.3/τ。如果要做0到5GHz的响应τ取100ps比较合适如果只关心2.4GHz附近的窄带响应τ取200ps就够了。τ越小脉冲越陡频谱越宽对网格分辨率的要求也越高。有一个容易被忽略的坑普通高斯脉冲带有直流分量频域结果在0Hz附近会有一个很大的峰干扰低频段响应。要避免这个问题改用调制高斯脉冲E(t) sin(2π * fc * (t - t0)) * exp(-((t - t0) / τ)²)其中fc是载波中心频率。这样频谱以fc为中心向两侧扩展去掉了直流分量。后面第5章讲频域毛刺排查时会再提到它。4.3 探针、时窗和收敛判断怎么看记录到的波形是否可信探针位置的选择直接影响结果可信度。我的习惯是至少放两个探针一个靠近源用来验证源注入正常一个放在关心区域比如介质圆后方用来提取透射波形。两个探针的波形在时延上要符合预期如果第二个探针比第一个还早收到信号明显是索引写错了。时窗长度要覆盖到源发射的脉冲彻底离开观测区域。判定标准很简单波形尾部衰减到峰值的1%以下。如果波形尾部还拖着尾巴就被截断后面做FFT时会因为吉布斯效应产生频域振荡这就是第5章要讲的毛刺问题。通常的做法是先跑一遍看波形趋势再把nt加大到波形衰减为止而不是一开始就猜迭代步数。另外一个值得养成的习惯是每改一次参数就记录下波形峰值和到达时间。这两个量在纯空气环境里可以用解析解估算如果对不上说明代码写错了某处继续往下跑只会浪费时间。5. FDTD常见翻车现象与排查路径发散、反射、毛刺和慢5.1 跑到一半数值变成NaN或inf先查CFL再查硬源现象代码前几百步正常某个时间步开始出现NaN刷屏停不下来。如果是inf则波形幅值直接冲到1e300量级。原因两个最典型的来源。第一是CFL条件不满足时间步长超过稳定上限误差指数放大第二是用了硬源源点被赋成脉冲值后阻挡了场传播在源点附近形成数值奇点。解决先把courant系数从0.9降到0.5跑一遍。如果不再发散说明是稳定性问题如果还是发散检查源是不是写成Ez(sx, sy) exp(...)这种覆盖式赋值改成Ez(sx, sy) Ez(sx, sy) exp(...)累加式。还有一个小概率原因是有耗介质的电导率项被显式处理更新方程里出现了除以σ的运算检查一下分母有没有可能为零。5.2 波形在边界处回弹观测点出现意料之外的反射峰现象零边界版本里观测波形在正常脉冲过后又冒出第二串幅度不小的振荡或PML版本里反射波衰减不明显波形尾部迟迟不归零。原因零边界版本出现反射是必然的不是bug是边界条件的属性。PML版本里反射明显通常是NPML层数太少或者σ_max取值偏小。另一种情况是PML更新方程里的分裂场项psi变量没有正确实现衰减只做了方向性处理部分角度的入射波会漏进吸收层。解决先做一维PML验证。把二维区域砍成一行网格跑一个平面波入射看反射系数是否达到-60dB以下。一维都过不了二维更不可能过。二维下如果仍然反射再加两层PML并把σ_max按第4章的公式从R01e-4重新算一遍。这个调试流程每次都管用比我见过的一上来就加厚PML层数要高效得多。5.3 FFT后频域曲线全是毛刺时窗截断和激励源带宽不匹配现象探针波形做FFT得到的频谱不是平滑曲线而是密密麻麻的锯齿振荡主峰旁边还有旁瓣。原因本质是时域信号被矩形窗截断波形尾部没有衰减到零就硬切断了。另一种情况是高斯脉冲τ选得太大频谱没有覆盖到关心的最高频率高频段信噪比极差。解决加大nt让波形自然衰减到接近零再停如果nt已经很大但波形仍在震荡说明边界反射在持续污染信号先解决PML问题。频域毛刺的另一个退路是对时域波形加窗我一般用汉宁窗主瓣展宽一点但旁瓣压得很低。窗口参数的配置逻辑是先保证物理上波形完整再用窗函数做最后一层保险顺序不能反。5.4 MATLAB跑得太慢一个算例等半小时循环和向量化的取舍现象算法正确但仿真时间长得无法接受网格稍微加密就没有耐心跑完。原因最常见的是把Yee网格更新写成了三层for循环嵌套。MATLAB处理循环的效率远低于对数组整体运算三层循环对二维网格就是nx*ny次内层操作数量级上就输给了切片运算。解决把循环改造成数组切片。第3章主循环里的写法Ez的更新就三行向量化表达式实际执行时MATLAB调用底层BLAS库速度比循环快一到两个量级。网格更新之外另一个性能陷阱是动态增长数组——每次循环往Probe数组后面追加元素MATLAB会反复分配内存。先Probe zeros(1, nt)预分配性能立刻提升。内存允许的情况下一次分配好所有场数组再迭代是最基础但收益最大的优化。5.5 谐振频率整体偏移百分之几怎么调都不对Yee网格的材料位置出错现象介质谐振结构的仿真结果比理论值偏低或偏高幅度也不对检查了CFL和PML都没问题。原因FDTD里eps_r和sigma不是定义在任意位置而是定义在Ez所在的整数格点上。Hx和Hy的更新需要的是磁导率mu虽然多数情况设为真空值但介质材料的相对介电常数必须放在Ez数组的相同索引位置。新手常犯的错误是把材料分布赋值在整个网格区域然后去更新所有场分量导致介质的电边界位置和理论位置差了半个网格。解决画一张Yee单元的示意图把Ez、Hx、Hy三个分量的位置标出来再对照自己的数组切片确认eps_r用的是Ez的下标。验证方法是把介质边界放到网格线的整数位置而不是半网格位置对比两次仿真的频移。如果频移量和半个网格对应的时间延迟相符说明就是这个问题。6. 进阶用一次仿真提取宽频带S参数并做网格收敛验证跑通时域波形只是第一步多数场景要的是频域响应。在FDTD里计算S参数不需要做多次频点扫描而是在时域仿真的同时对观测点信号做离散傅里叶变换一次激励获得整个频带的结果这就是DFT监视点的思路。代码上其实很轻量% 在时间循环外预先定义频率轴和DFT累加器 f linspace(fmin, fmax, 201); FDFT zeros(size(f)); for n 1:nt % ...主循环里的场更新... % 对观测点处的场做DFT累加 FDFT FDFT Ez(px, py) * exp(-1i * 2 * pi * f(:) * n * dt); end % 用入射波参考信号归一化得到传输系数 S21 abs(FDFT ./ FDFT_ref);这行的关键是exp(-1i * 2 * pi * f(:) * n * dt)做了向量化乘法频率点越多优势越明显。如果f有几百个点整个循环跑完也就百万级运算对MATLAB毫无压力。FDFT_ref的取值我习惯在去掉介质结构后单独跑一遍在相同位置记录参考场再拿有结构时的结果除以它这样能直接消除源频谱形状的影响。网格收敛验证是检验仿真可信度的硬指标。把dx减半、dy也减半时间步按CFL重算重新跑一遍对比两次S21曲线。如果两条曲线基本重合说明空间离散误差已经小于你的容忍度如果还有明显偏差继续加密直到收敛。这个验证方法不花太多额外时间但它能回答一个问题仿真结果到底能不能代表真实物理。如果只知道跑代码而不做收敛测试任何结果都只是代码的输出不是计算电磁学的结论。这些年我的习惯是拿到任何一份“时域有限差分的MATLAB代码.zip”第一件事不是去跑大算例而是先用一维或二维的解析解验证代码的正确性。固定CFL0.9、PPW不低于12是保底操作PML参数按公式算而不是靠猜最后一定要做一次网格收敛测试再下结论。这套习惯救过我很多次希望你也能用上。希望帮到你。本文还有配套的精品资源点击获取
