企业官网建设流程全解析

简介：这是一套基于时域有限差分（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(i+1) - 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(i+1)，Ez(i) 看 Hy(i) 和 Hy(i-1)，正好是半个网格的错位，照抄即可。两个更新公式里的系数 dt/(mu0dx) 和 dt/(eps0dx) 是缩放因子，它们的准确性直接影响波速，所以 eps0 和 mu0 必须用真实值，不能随手填 1。

参数上，dx 决定最高可分辨频率。波动性问题里 dx 一般不超过最小波长的 1/20，这里 dx=1mm 对应最高有效频率约 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 吸收边界，在左右两个端点分别用前一时刻的场来近似外向波：

% 左边界（i=1）的一阶Mur吸收 Ez(1) = Ez(2) + (c*dt - dx)/(c*dt + dx) * (Ez_old(2) - Ez(1)); % 右边界（i=nx）的一阶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. 二维FDTD：TE/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模点源辐射：完整的可运行代码

% 二维FDTD：TM模（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,j+1) - 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(i+1,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 必须同时用两个方向步长参与计算。spread=8 的高斯脉冲在空间上大约覆盖 16 个网格，如果这个值太小，比如小于 3，频谱太宽，网格分辨率跟不上，你会看到明显的各向异性伪影，波形呈方形扩散而不是圆形。检查这种伪影最好的办法是看可视化结果：各向同性介质中点源辐射应该是标准圆形波前，一旦出现菱形或方形轮廓，不是 dx 太大就是 CFL 系数给得太满。

3.3 网格分辨率怎么取：dx不是越小越好

dx 的选择要同时权衡计算量、数值色散和结构精度。经验法则是 dx ≤ λ_min / 20，λ_min 是模拟频段内的最短波长。可见光波段若用中心波长 550nm，dx 取 25nm 以下；微波 3GHz，dx 取 5mm 以下即可。数值色散是更隐蔽的约束：当 dx 相对波长不够小时，不同传播方向的波会有不同的相速度，导致波形在传播中畸变，脉冲越走越宽甚至出现拖尾。

注意 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 字节。我列出几个典型规模，你感受一下增长速度：

从这个表能看出，三维比二维不是多一个维度，而是多一到两个数量级的内存和计算量。很多人第一次在三维 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 方向差分。一旦搞反，仿真结果会非常怪异。常见错误是把二维代码里的"j=2:ny-1"直接抄成"k=2: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 厚度取 10～15 层，电导率用多项式渐变，并在结构边缘与 PML 之间留出 4～6 个网格的缓冲区。判断边界是否合格的方法是测反射系数：在边界前放探针，比较反射波峰值与入射波峰值的比值，通常要求 -40dB 以下。

5.3 透射率峰位对不上解析解或实验值

现象：计算的谐振波长偏了几百纳米，或者透射率曲线整体左移右移，怎么调参数都差一截。

原因：单位制混乱占大头。有人用纳米定义几何，但 dx 按米写成了 1e-9，结构尺寸却用微秒数直接乘 dx，结果尺寸错了几百倍。另一个常见错误是把折射率当相对介电常数用，把 1.5 直接代入 eps，导致波速算错，所有谐振峰集体偏移。

解决：在代码开头统一单位，定义 scale = 1e-6 表示微米，dx = 0.05*scale，所有几何尺寸都用"格点数 × dx"表达。折射率 n=1.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 一阶吸收边界干净得多。参数经验值：层数取 10～15 层，电导率从内到外按多项式渐变，p 取 3～4，最大电导率约取 sigma_max = 0.8*(p+1)/(150pidx)。CPML 的辅助变量要单独开数组，数量与场分量相同，内存占比不小，用单精度存储即可。比较容易被忽略的是缓冲区：结构边缘与 PML 之间留 4～6 个网格，让波在到达 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 上依旧慢，必须先做向量化，或者至少把最内层循环改成 arrayfun，GPU 才有明显收益。MATLAB 在 GPU 上执行时会自动管理显存，但要注意显存容量：200×200×200 的双精度 6 分量数组已经接近 6GB，家用显卡很可能放不下。实践上我通常把二维和三维的数组都转成 single 再给 GPU，速度更快，内存压力也更小。向量化加单精度这一套组合下来，二维点源辐射代码提速 5～10 倍是很正常的。

6.3 一个值得养成的验证习惯

每次写新结构前，先跑一个已知解的最小模型。一维平板透射、二维圆柱散射、三维球体 Mie 散射都可以，算好后与解析解对比，误差在几个百分点内再继续下一步。这个习惯成本很低，但能把"代码 bug"和"物理参数错误"清晰分开，省下的调试时间非常可观。我现在跑新结构，第一步永远是写一个五分钟能跑完的最小模型验证边界和源没有问题，这个习惯救了我很多次。希望帮到你。

本文还有配套的精品资源，点击获取