简介一套基于MATLAB的FDTD有限差分时域电磁场计算入门工具面向电磁仿真初学者、天线与微波工程方向的学生和工程师解决不知如何在MATLAB中落地FDTD算法的问题。压缩包内仅有1个MATLAB脚本文件.m体积约1KB代码精简聚焦Yee网格空间离散、电场磁场时间步交替更新、边界条件处理等FDTD核心环节尤其适合先读懂主干再逐步扩展的学习路径。目前已有314人学习下载。借助该脚本用户可快速搭建二维/三维简单FDTD仿真环境直观观察电磁波在网格中的传播与散射过程通过修改源项波形、材料参数或加入吸收边界还能延伸验证PML边界、色散媒质、谐振结构等进阶内容。对于想掌握如何用MATLAB实现FDTD的入门者这份代码比纯理论推导更具体可作为第一份可运行、可调试、可改写的参考实现为后续复杂电磁场仿真打下基础。1. 一个 FDTD.m 背后为什么 MATLAB 适合做时域电磁场计算很多人在搜“fdtd怎么matlab”时并不是想先啃完一整本电磁场数值方法而是手里有一个从 FDTD.rar 解压出来的 FDTD.m想知道这个脚本到底在算什么、怎么跑、改成什么参数能让结果看起来合理。这个脚本用有限差分时域FDTD方法求解麦克斯韦方程组核心逻辑其实只有十几行电场和磁场交替更新再加一个激励源和一组边界条件。它适合两类人一类是刚接触电磁仿真的学生想通过改网格步长和介电常数理解波的传播与反射另一类是用过 Lumerical FDTD 或 HFSS 的工程师需要不依赖 License、能在 MATLAB 里快速验证物理想法的参考实现。MATLAB 的优势在于数组运算和可视化让场更新更像是矩阵差分而不是写三重 for 循环。2. 从麦克斯韦到 Yee 网格读懂 FDTD.m 的核心更新公式2.1 交错网格与半步更新FDTD 的底层逻辑FDTD 的基本操作是把空间切成一个个网格然后在格点上近似旋度。Yee 网格的核心是让电场和磁场在空间上错开半个网格时间上也错开半个步长。这样一来麦克斯韦方程组中的时间偏导可以用中心差分写成显式更新不需要求解矩阵方程。以二维 TMz 波为例电场只有 Ez 分量磁场有 Hx 和 Hy 分量。电场节点在整数格点上磁场节点落在半网格位置。MATLAB 里没有“半个索引”但可以用数组移位来模拟Hx(i,j) 实际表示 (i, j0.5) 处的磁场Hy(i,j) 表示 (i0.5, j) 处的磁场。这种交错结构让每次更新都只涉及相邻节点天然适合向量化也让新手更容易把物理公式映射到矩阵索引上。先看一维 FDTD 的更新公式理解它之后再扩展到二维会顺很多% 一维FDTD电场ez沿x方向传播 for n 1:Nt % 更新磁场相邻电场差产生磁场 for i 1:nx-1 hy(i) hy(i) (dt/(mu0*dx)) * (ez(i) - ez(i1)); end % 更新电场相邻磁场差产生电场 for i 2:nx-1 ez(i) ez(i) (dt/(eps0*dx)) * (hy(i) - hy(i-1)); end % 源点激励源位置和波形可自行修改 ez(source_pos) ez(source_pos) amp * sin(2*pi*freq*n*dt); end这段代码里hy(i) 和 hy(i-1) 在空间上正好相差一个网格中心差分的系数是 dt/(mu0*dx)。注意电场更新时索引从 2 开始到 nx-1 结束这是为了不越界两个端点的电场值需要由边界条件单独处理而不是直接做差分。如果你拿到的 FDTD.m 里也有类似的双循环说明作者保留了最容易读懂的写法。2.2 把更新方程写成 MATLAB 数组运算上面的双层 for 循环在网格数超过 200 后速度会明显下降所以更常见的 FDTD.m 实现会改成向量化写法。对于一维问题可以用 diff 函数一次算出全部相邻差% 电场和磁场数组长度不同diff自动处理相邻差 hy hy - (dt/(mu0*dx)) * diff(ez); ez(2:end-1) ez(2:end-1) - (dt/(eps0*dx)) * diff(hy);diff(ez) 返回 ez(2)-ez(1)、ez(3)-ez(2)……长度比 ez 少 1恰好和 hy 的个数匹配。第二行的 diff(hy) 长度也比 hy 少 1所以电场只更新中间点。这个思路扩展到二维后就是把 diff 放在不同维度上Ez(2:end-1,2:end-1) Ez(2:end-1,2:end-1) ... (dt/(eps0*dx)) * (diff(Hy,1,1) - diff(Hx,1,2));diff(Hy,1,1) 表示在行方向做一阶差分diff(Hx,1,2) 表示在列方向做差分。两个结果形状不同但组合起来正好对齐电场矩阵的中间区域。如果你在 FDTD.m 里看到类似代码说明作者已经做过向量化优化如果你是第一次写优先保证逻辑正确再考虑性能。2.3 源、CFL 条件和材料参数改哪些会改变波形运行 FDTD.m 时最常改的三个地方是源、网格步长和材料参数。下表给出了常见参数的参考范围和它们对结果的影响参数常见值含义与影响dx0.001~0.01 m空间网格步长决定分辨率越小越精确但耗时指数增长dtdx/(2c)时间步长必须满足 CFL 条件才能稳定Nx、Ny100~500网格数量决定计算区域大小Nt500~5000总迭代步数决定模拟的物理时间长度epsr1~12相对介电常数决定波速和界面反射source_freq1e9~10e9激励频率中心频率要低于网格能够分辨的最高频率源的波形也很关键。高斯脉冲适合宽带计算因为只需跑一次就能获取多个频点的响应正弦波适合观察特定频率下的驻波和反射。在 MATLAB 中高斯脉冲可以这样加入ez(source_pos) ez(source_pos) exp(-((n-n0)/tau)^2);n0 是脉冲中心时间步tau 是脉冲宽度。tau 越小脉冲包含的高频成分越多空间步长 dx 就要越小否则会出现明显的数值色散。这个点经常被忽略很多人只调小 tau发现波形乱了其实是 dx 没有跟着缩小。3. 把 FDTD.m 跑起来参数设置、可视化与排查3.1 解压 FDTD.rar 后先做三件事解压得到 FDTD.m 后不要急着直接点运行。先打开文件确认三件事开头有没有清空变量的语句主循环外有没有定义 eps0、mu0、dx、dt边界条件写在哪里。很多参考脚本会在开头写 clear all; close all;这是为了避免上一次运行残留的变量影响当前结果。如果文件里没有可以在运行前手动加一行。确认参数后在 MATLAB 命令窗口进入解压目录并运行脚本cd(你的解压路径); run(FDTD.m);如果脚本没有报错工作区里会出现 Ez、Hx、Hy、dx、dt 等变量。用 whos 查看它们的尺寸whos Ez Hx Hy正常的 Ez 维度应该和网格数 Nx*Ny 对应。如果你看到 Ez 是 200x1 而不是二维矩阵说明这是一维模型坐标索引方式和二维不同。先看清数组维度再决定从哪里改参数。3.2 运行脚本并实时看场imagesc 与 getframe很多 FDTD.m 会只在循环结束后画一张最终场图但调试阶段更需要看波的传播过程。常见做法是把绘图放进时间循环每几十步刷新一次figure(1); for n 1:Nt % FDTD 更新代码 if mod(n, 50) 0 imagesc(squeeze(Ez)); colormap(jet); colorbar; caxis([-0.1 0.1]); title(sprintf(Time step %d, n)); drawnow; end endimagesc 把 Ez 矩阵显示为二维色图squeeze 用来去掉多余的单一维度。caxis 把色标范围固定住否则每帧都会自动缩放入射波和反射波混在一起时很难分辨幅度差异。drawnow 放在 if 语句里可以避免每步都刷新从而减少绘图对计算速度的影响。如果想输出动画可以把 drawnow 换成 getframe 并把帧保存到数组。注意 getframe 本身会触发图窗刷新所以不要再重复调用 drawnow否则画面会闪烁。3.3 波形发散、反射大、结果不对怎么办运行后最常见的三种异常是发散、边界反射和波形畸变。下表列出了排查方向现象可能原因检查方法场值变成 NaN 或 Infdt 超过 CFL 条件检查 dt 是否小于 dx/(sqrt(2)*c0)波形碰到边界后反弹缺少吸收边界加大计算区域或补吸收层脉冲变宽或拖尾网格步长太大每个波长至少 10 个网格点数值发散的判断很简单直接在循环外校验c0 1/sqrt(eps0*mu0); dt_max dx / (sqrt(2) * c0); if dt dt_max error(dt 太大请改用 dt dx/(3*c0) 再试); end边界反射大的表现是波形传到边界后反弹叠加到真实信号上。最快补救办法是扩大网格区域让反射波在观察时间内不回到源区。更完整的做法是用吸收边界或 PML。波形畸变的原因通常是源附近网格太粗。经验法则是每个最小波长内至少 10 个网格点。用lambda_min c0 / (freq_max)估算最小波长然后要求dx lambda_min / 10。4. 从一维到二维FDTD.m 如何扩展成面阵天线模型4.1 二维网格的搭建和索引对应一维 FDTD 只用一行数组二维 TMz 模型则要把电场定义成矩阵。常见做法是把 Ez 设成 Nx*Ny 的矩阵Hx 和 Hy 在不同方向错开半个网格Nx 200; Ny 200; Ez zeros(Nx, Ny); Hx zeros(Nx, Ny-1); % Hx位于(i,j0.5) Hy zeros(Nx-1, Ny); % Hy位于(i0.5,j) eps eps0 * ones(Nx, Ny); mu mu0 * ones(Nx, Ny);这里的 Hx 和 Hy 尺寸各少一列或一行正是为了匹配差分后的结果。更新 Hy 时用 diff 得到的数组自然对应前 Nx-1 行Hy(1:end-1,:) Hy(1:end-1,:) - (dt./mu(1:end-1,:)) .* diff(Ez,1,1);这里必须使用./而不是/因为 mu 是矩阵。如果你把介质从均匀真空改成多层介质这个点最容易出错。4.2 在同一套循环里支持真空和介质区域要在二维模型里放置介质块最直接的做法是维护一个介电常数矩阵 epsr把更新系数提前算好放在循环外C_E dt ./ (eps0 * epsr); % 电场更新系数逐点不同 C_H dt / (mu0 * dx); % 磁场更新系数真空区域相同然后在时间循环里直接乘以这些系数。这样做的好处是介质块边界不会因为源的位置而改变更新逻辑清晰。缺点是在介电常数变化剧烈的边界上电场更新应使用两侧介电常数的平均值简单方案会产生阶梯误差。对于初学先用逐点 epsr 也能得到合理波形之后再考虑边界平均。4.3 Mur一阶吸收边界二维扩展的第一步二维模型必须处理边界反射最简单的方案是 Mur 一阶吸收边界。以左边界为例电场更新可以写成Ez(1,2:end-1) Ez(2,2:end-1) ... (c0*dt - dx)/(c0*dt dx) * (Ez(2,2:end-1) - Ez(1,2:end-1));这个公式用相邻点的当前值推算边界外侧的场本质上是在假设波垂直边界向外传播。其余三条边界需要分别写公式结构相同只是方向不同。Mur 边界在斜入射时效果一般下表给出了几种边界处理方式的对比边界类型垂直入射反射斜入射效果实现复杂度固定边界100%100%最低Mur一阶5%~10%较差低渐变衰减层小于5%一般中标准 PML小于1%好高对于只观察短时段内的前几次反射Mur 边界已经够用如果要做驻波比或谐振分析还是得用 PML。4.4 观察近场记录探针时间波形二维场图适合看空间分布定量分析则需要放探针。探针点在每个时间步记录 Ez 值循环结束后画成波形probe_x 100; probe_y 120; probe_Ez zeros(1, Nt); for n 1:Nt % FDTD 更新代码 probe_Ez(n) Ez(probe_x, probe_y); end plot((1:Nt)*dt*1e9, probe_Ez); xlabel(时间 (ns)); ylabel(Ez (V/m));探针波形能直接读出入射波和反射波到达的时间差。如果知道介质厚度就能算出波速再反推介电常数设置是否正确。这和 Lumerical FDTD 里的 time monitor 是同一个思路只是实现简单很多。5. 从验证到进阶让 FDTD.m 拥有 PML 和远场外推二维 FDTD 能跑通后下一步是把边界反射进一步压下去。Mur 一阶边界只能应付接近垂直入射的波真正工程上常用的是 PML。PML 的核心思想是在边界外增加一层有损介质让进入的波指数衰减同时保证介质阻抗匹配不产生额外反射。MATLAB 里最容易验证衰减机制的做法是在边界区域设置渐变电导率。为了让吸收层和真空阻抗匹配还需要引入等效磁导率损耗项。电导率从内到外逐渐增大避免突变反射sigma_e 0.02 * ones(Nx, Ny); sigma_e(:, 1:10) 0.02 * ((10:-1:1)/10).^2; % 渐变导电率层 Ez Ez .* exp(-sigma_e * dt / eps0);这行代码并不是严格 PML因为真正的 PML 需要分裂场或坐标伸缩坐标代码量会翻倍。但指数衰减形式能帮你理解“损耗层为什么能压反射”。验证时放一个高斯脉冲记录边界处探针的返回峰值。不加吸收层时峰值可能高达入射波的 50% 以上加了渐变衰减层后应降到 5% 以下。远场外推则需要额外做一步在源周围取一个闭合矩形框记录框上的切向电场和磁场。循环结束后用 fft 提取目标频点的复振幅再按等效原理叠加得到远场方向图。MATLAB 的 fft 可以直接处理时间序列不需要自己推导格林函数。如果你身边有 Lumentum 或 Lumerical FDTD 的测试环境可以把同一个平板模型搬进去比较单点探针波形。两者波的传播趋势应该一致差别主要在数值色散和边界处理上。这个对照能帮你确认 MATLAB 脚本的介质参数和边界条件没有写错也方便后续把 FDTD.m 改成自己的天线或滤波器模型。本文还有配套的精品资源点击获取