简介一套基于C语言实现的时域有限差分法FDTD计算器源码面向电磁仿真学习者、课程设计学生与数值计算开发者聚焦二维真空中电磁波传播采用正弦波激励并融入完美匹配层PML边界条件。压缩包共15个文件核心源码为fetd_32.cpp同时包含可直接运行的exe、调试期生成的obj/pdb/ilk、以及dsp/dsw/opt/ncb/plg等工程配置文件整体仅231KB轻量且便于研读。该资源已有103人学习下载。通过学习源码可掌握FDTD在Yee网格上的场量更新流程、PML吸收边界的实现思路配合程序输出的Ez.txt电场数据能直观观察电磁波在计算域中的传播与吸收效果还能学习项目在VC6.0下的组织方式与调试信息生成结构既适合初学者快速建立数值仿真认知也便于进阶者在此基础上扩展复杂介质、调整激励源或做并行优化。1. FDTD_32 是什么用 C 语言实现的电磁场计算器FDTD_32 是一个用 C 语言写的二维时域有限差分计算器。它不做四则运算而是把空间网格上的电场 Ez、磁场 Hx、Hy 反复更新模拟电磁波在介质中的传播、反射和耦合。名字里的 32 指默认的 32×32 网格尺寸也暗示这是一个学习型小工程适合 C 语言课、毕业设计或者射频岗位笔试前的速成练习。你只需要一个本地 gcc 或 VS Code 的 C 语言开发环境输入网格参数、波源位置和仿真步数运行后就能得到 CSV 格式的场快照。这篇文章会从麦克斯韦方程的差分形式一路写到 Makefile中间所有代码都可以照着敲。2. FDTD 内核怎么建Yee 网格、更新方程与 CFL 条件FDTD 的核心是把麦克斯韦 curl 方程在时间域做有限差分。空间被均匀切成长方形网格时间以 dt 步进电场和磁场交错更新。先更新磁场再更新电场边走边模拟。每一轮的计算量只取决于网格数量没有矩阵求逆所以特别适合 C 语言实现。2.1 Yee 网格把场量错开半个网格点Yee 网格是 FDTD 能长时间稳定计算的骨架。电场放在整数格点上磁场放在半格点上中心差分天然具有二阶精度。二维 TMz 模型只需要电场 Ez 和两个磁场 Hx、Hy变量少网格小适合入门。下面是我常用的存放约定。场量实际坐标更新它时要读的邻近场Ez(i, j)Hy(i-1, j)、Hy(i1, j)、Hx(i, j-1)、Hx(i, j1)Hx(i, j-1/2)Ez(i, j-1)、Ez(i, j)Hy(i-1/2, j)Ez(i-1, j)、Ez(i, j)这套偏移关系在代码里不需要真的使用半格下标。只需要记住 Hx 的数组下标比 Ez 低半格Hy 的下标也比 Ez 低半格更新时统一偏移一位即可。使用一维连续数组存储时换算公式固定为index j * nx i读写速度比二维指针数组更快。2.2 二维 TMz 波的电场磁场更新公式真空中无损耗介质的旋度方程展开后得到三个标量方程∂Hx/∂t -μ⁻¹ ∂Ez/∂y∂Hy/∂t μ⁻¹ ∂Ez/∂x∂Ez/∂t ε⁻¹ (∂Hy/∂x - ∂Hx/∂y)用中心差分替换空间导数后C 语言更新函数可以写成下面这样。/* 均匀网格dx dy 时系数可以合并 */ const double c1 dt / (mu * dy); /* Hx 更新系数 */ const double c2 dt / (mu * dx); /* Hy 更新系数 */ const double c3 dt / (eps * dx); /* Ez 中 Hy 差分系数 */ const double c4 dt / (eps * dy); /* Ez 中 Hx 差分系数 */ void update_h(FDTD2D *g) { for (int j 1; j g-ny; j) { for (int i 1; i g-nx; i) { g-hx[j * g-nx i] c1 * ( g-ez[j * g-nx i] - g-ez[(j - 1) * g-nx i]); g-hy[j * g-nx i] c2 * ( g-ez[j * g-nx i] - g-ez[j * g-nx (i - 1)]); } } } void update_e(FDTD2D *g) { for (int j 1; j g-ny - 1; j) { for (int i 1; i g-nx - 1; i) { double dhy g-hy[j * g-nx (i 1)] - g-hy[j * g-nx i]; double dhx g-hx[(j 1) * g-nx i] - g-hx[j * g-nx i]; g-ez[j * g-nx i] c3 * dhy - c4 * dhx; } } }磁场更新依赖前后两个 Ez 的差值电场更新依赖两个 Hy 和一个 Hx 的差值。为了让边界上的点不越界循环索引从 1 开始网格最外一圈不在更新范围内而由吸收边界单独处理。细心的读者会发现 Hx 的更新式里没有负号Hy 的更新式也是正号。这取决于坐标轴方向和 Yee 网格的偏移方向。如果你从另一本教材里抄来的公式带负号不要急着改程序先统一坐标轴定义。FDTD 本身没有方向偏好只要符号与坐标一致最后结果就是相同的。2.3 CFL 条件dt 不能随意设显式 FDTD 的数值波速必须大于等于物理光速否则误差波会追上真实波导致结果发散。二维均匀网格的稳定性条件为dt 1 / (c * sqrt(1 / dx² 1 / dy²))在dx dy时这个上限简化为dt dx / (c * sqrt(2))。实际工程里我习惯取上限的 0.9留出安全余量。后面如果加入介质层光速变小CFL 上限会略微放宽但调试阶段不要一上来就用接近极限的值。先跑通真空模型再慢慢加介质。3. 用 C 语言把 FDTD_32 组织成完整源码FDTD 计算器的代码结构比算法本身更容易劝退人。我习惯用一个结构体把网格参数和三个场数组全部包起来所有更新函数只接收一个FDTD2D *指针。这样后续加 PML、加金属边界、加介质层都不需要改动主循环的函数签名。3.1 用结构体承载计算器参数先看核心结构体。typedef struct { int nx, ny; /* 网格尺寸 */ double dx, dy, dt; /* 空间步长和时间步长 */ double t; /* 当前仿真时间 */ double *ez, *hx, *hy; /* 场数组一维连续内存 */ } FDTD2D;选择一维数组而不是二维数组是为了减少malloc的次数并且让同一行的场数据在缓存里连续。nx、ny后面会频繁参与下标换算如果写死 32以后就要全文件搜索修改。把它们放进结构体后用指针传递是 C 语言项目的标准做法也锻炼了c语言结构体和c语言指针的配合能力。3.2 malloc 与 free网格数组的内存管理创建和销毁函数如下。FDTD2D *grid_create(int nx, int ny, double dx, double dy) { FDTD2D *g (FDTD2D *)malloc(sizeof(FDTD2D)); if (g NULL) { perror(malloc FDTD2D); exit(EXIT_FAILURE); } g-nx nx; g-ny ny; g-dx dx; g-dy dy; g-t 0.0; g-ez (double *)calloc((size_t)nx * ny, sizeof(double)); g-hx (double *)calloc((size_t)nx * ny, sizeof(double)); g-hy (double *)calloc((size_t)nx * ny, sizeof(double)); return g; } void grid_free(FDTD2D *g) { free(g-ez); free(g-hx); free(g-hy); free(g); }calloc会把整块内存清零比malloc加memset更顺手。数组长度用了(size_t)nx * ny避免在 32 位系统上做int乘法时溢出。创建和释放写成独立函数比在main里手工维护安全得多。这也是网上那批c语言必背100代码里最常见的题型构造函数配析构函数不让内存泄漏。3.3 波源设置与简单吸收边界波源通常使用高斯脉冲脉冲宽度由时间常数tau控制。void add_gaussian_source(FDTD2D *g, int sx, int sy, double t0, double tau, double amp) { double pulse amp * exp(-pow((g-t - t0) / tau, 2)); g-ez[sy * g-nx sx] pulse; }波源加在电场数组上物理上是一个电流源。如果加在磁场数组上就变成磁流源两者辐射方向场不同。学习程序里不建议直接放正弦源因为正弦源会立刻把边界打爆初学者根本分不清是边界反射还是算法发散。高斯脉冲从零平滑上升到最高点再降下来边界吸收的压力小很多。吸收边界我用简化版的海绵层不实现完整的 PML。把电场在网格边缘按线性因子缩小让外行波被慢慢吃掉。void apply_absorbing_taper(FDTD2D *g, int thickness) { for (int j 0; j g-ny; j) { for (int i 0; i g-nx; i) { double factor 1.0; if (i thickness) factor (double)i / thickness; else if (i g-nx - thickness) factor (double)(g-nx - 1 - i) / thickness; if (j thickness) factor (double)j / thickness; else if (j g-ny - thickness) factor (double)(g-ny - 1 - j) / thickness; g-ez[j * g-nx i] * factor; } } }这个海绵层会吸收大部分向外传播的能量但也会造成轻微反射。对 32×32 的演示网格足够真正研究回波损耗还是要上 PML。thickness一般取 4 到 8太小没有衰减效果太大则压缩有效计算域。3.4 把场快照写入 CSV 文件输出用纯文本 CSV方便后续交给 gnuplot 或 Python 出图。void save_ez_csv(FDTD2D *g, const char *path) { FILE *fp fopen(path, w); if (fp NULL) { perror(path); return; } for (int j 0; j g-ny; j) { for (int i 0; i g-nx; i) { fprintf(fp, %g%c, g-ez[j * g-nx i], (i g-nx - 1) ? \n : ,); } } fclose(fp); }fprintf使用%g而不是%f避免输出一堆没有物理意义的零。每个 CSV 文件的行数等于ny列数等于nx。文件读写操作里一定要检查fopen的返回值否则磁盘满或目录不存在时会出现静默失败。这一步踩过的坑和c语言文件读写操作代码的热搜里描述的问题基本一致。3.5 用 Makefile 一键编译运行一段最小的 Makefile 足够在 Linux、macOS 和 Windows 的 MSYS2 里使用。CC gcc CFLAGS -stdc99 -O2 -Wall LDLIBS -lm main: main.c $(CC) $(CFLAGS) -o main main.c $(LDLIBS) clean: rm -f main在 VS Code 配置 C 语言环境时tasks.json里把编译命令指定成make main再配置一个 build task就能一键编译。优化级别-O2很关键FDTD 是循环密集型程序-O0下同一套参数可能要慢两三倍。接下来是主循环的大致样子。int main(void) { FDTD2D *g grid_create(32, 32, 1e-3, 1e-3); double dt_max 1.0 / (C0 * sqrt(1.0 / (g-dx * g-dx) 1.0 / (g-dy * g-dy))); g-dt dt_max * 0.9; for (int n 0; n 500; n) { update_h(g); update_e(g); add_gaussian_source(g, 16, 16, 5e-11, 1e-11, 1.0); apply_absorbing_taper(g, 4); g-t g-dt; if (n % 50 0) save_ez_csv(g, field.csv); } grid_free(g); return 0; }主循环里电磁场交替更新波源在电场更新之后叠加。如果你把源叠加放在电场更新之前等价于把电流密度积分顺序换了一下波形上会出现半个时间步的偏移。先按这个顺序跑通再考虑优化。4. 参数整定与排错覆盖 CFL、网格尺寸和内存陷阱上一章的源码能跑但直接运行 32×32 网格如果参数不对屏幕上很快就会刷出 NaN。参数整定是 FDTD 从玩具变成实用工具的关键。4.1 常用参数表与推荐起点下面是一组能直接跑的起点参数。参数推荐值说明网格尺寸 nx, ny32 × 32默认值先跑通再放大空间步长 dx, dy1e-3 m每格 1 mm适合分米波演示时间步长 dt2.1e-12 s由 CFL 上限乘 0.9高斯源 t05 * tau让脉冲从零平滑上升高斯源 tau10 * dt覆盖 10 个左右时间步源位置16, 16网格正中间吸收层厚度432 网格下取 4 较合适空间步长决定了空间解析度。经验法则是每个波长至少要有 10 到 20 个网格。32×32 不适合分析精细结构但足够清楚展示波前传播。这样也能和二维波的解析解对比检查数值色散是否在合理范围内。4.2 一次完整的 CFL 计算例子以下代码展示了如何在main里计算安全dt。#define C0 3e8 double dx 1e-3, dy 1e-3; double dt_max 1.0 / (C0 * sqrt(1.0 / (dx * dx) 1.0 / (dy * dy))); double dt dt_max * 0.9; printf(dt_max %g s, dt %g s\n, dt_max, dt);把dx1e-3, dy1e-3代入1/dx²等于 1e6两项加起来是 2e6开根后约 1414。dt_max约等于 2.36e-12 秒取 0.9 倍后是 2.1e-12 秒。这个量级在编码时经常因为单位写错而差出十个数量级所以我初始化时会先打印一次dt再进入主循环。4.3 计算器不对时的检查顺序如果场值发散按下面顺序排查。检查dt是否满足 CFL。把dt调成上限的 0.5如果立刻正常问题就是时间步长太大。检查边界。apply_absorbing_taper必须在每次电场更新后调用漏掉的话能量会在边界积累。检查源的位置。源必须落在 1 到nx-2、1 到ny-2的范围内否则源点本身处于不更新区域能量无法进入求解域。检查数组类型。场数组必须是double*不要用int*整型数组会把场值截断成 0 或 1画出来不是波形而是噪点。用c语言的valgrind跑小步数定位越界读写。FDTD 数组下标换算最容易在循环边界差一拍常见错误是把i nx写成i nx这会读到ez数组后面的邻居内存。这五步覆盖大部分初学问题。如果还不对就回到最小模型去掉源只保留初始场看总能量是否逐步下降。这样能区分算法问题和数据问题。5. 验证技巧给 FDTD_32 加能量探针并用 gnuplot 成图最后一个技巧是给计算器加一个总能量探针。每次电场更新完后扫描全网格把电场和磁场能量加起来。double total_energy(FDTD2D *g) { double sum 0.0; for (int j 0; j g-ny; j) { for (int i 0; i g-nx; i) { double e g-ez[j * g-nx i]; double hx g-hx[j * g-nx i]; double hy g-hy[j * g-nx i]; sum 0.5 * e * e 0.5 * (hx * hx hy * hy); } } return sum; }在main循环里每 10 步写一行time energy到energy.dat。如果海绵吸收边界有效能量曲线会单调下降但不会出现突兀跳变如果指数上涨大概率是 CFL 被破坏。接下来用 gnuplot 直接看热力图。gnuplot -e set terminal png; set output ez.png; \ set view map; splot field.csv matrix with imagematrix关键字告诉 gnuplotCSV 是按网格行组织的二维数组。跑完后如果能看到圆形波前从中心扩散并在边缘被吸收层压平说明从公式到源码的整条链路是通的。把能量探针再扩展一步还能在固定点记录 Ez 随时间变化用离散傅里叶变换扫出幅度谱进而找谐振频率。FDTD_32 作为计算器最实用的下一步就是这个频域响应分析。本文还有配套的精品资源点击获取 
长河君 · 2026/9/15 17:31:35 · 阅读 1.2k 用C语言实现二维FDTD电磁场模拟：从麦克斯韦方程到代码实践
简介一套基于C语言实现的时域有限差分法FDTD计算器源码面向电磁仿真学习者、课程设计学生与数值计算开发者聚焦二维真空中电磁波传播采用正弦波激励并融入完美匹配层PML边界条件。压缩包共15个文件核心源码为fetd_32.cpp同时包含可直接运行的exe、调试期生成的obj/pdb/ilk、以及dsp/dsw/opt/ncb/plg等工程配置文件整体仅231KB轻量且便于研读。该资源已有103人学习下载。通过学习源码可掌握FDTD在Yee网格上的场量更新流程、PML吸收边界的实现思路配合程序输出的Ez.txt电场数据能直观观察电磁波在计算域中的传播与吸收效果还能学习项目在VC6.0下的组织方式与调试信息生成结构既适合初学者快速建立数值仿真认知也便于进阶者在此基础上扩展复杂介质、调整激励源或做并行优化。1. FDTD_32 是什么用 C 语言实现的电磁场计算器FDTD_32 是一个用 C 语言写的二维时域有限差分计算器。它不做四则运算而是把空间网格上的电场 Ez、磁场 Hx、Hy 反复更新模拟电磁波在介质中的传播、反射和耦合。名字里的 32 指默认的 32×32 网格尺寸也暗示这是一个学习型小工程适合 C 语言课、毕业设计或者射频岗位笔试前的速成练习。你只需要一个本地 gcc 或 VS Code 的 C 语言开发环境输入网格参数、波源位置和仿真步数运行后就能得到 CSV 格式的场快照。这篇文章会从麦克斯韦方程的差分形式一路写到 Makefile中间所有代码都可以照着敲。2. FDTD 内核怎么建Yee 网格、更新方程与 CFL 条件FDTD 的核心是把麦克斯韦 curl 方程在时间域做有限差分。空间被均匀切成长方形网格时间以 dt 步进电场和磁场交错更新。先更新磁场再更新电场边走边模拟。每一轮的计算量只取决于网格数量没有矩阵求逆所以特别适合 C 语言实现。2.1 Yee 网格把场量错开半个网格点Yee 网格是 FDTD 能长时间稳定计算的骨架。电场放在整数格点上磁场放在半格点上中心差分天然具有二阶精度。二维 TMz 模型只需要电场 Ez 和两个磁场 Hx、Hy变量少网格小适合入门。下面是我常用的存放约定。场量实际坐标更新它时要读的邻近场Ez(i, j)Hy(i-1, j)、Hy(i1, j)、Hx(i, j-1)、Hx(i, j1)Hx(i, j-1/2)Ez(i, j-1)、Ez(i, j)Hy(i-1/2, j)Ez(i-1, j)、Ez(i, j)这套偏移关系在代码里不需要真的使用半格下标。只需要记住 Hx 的数组下标比 Ez 低半格Hy 的下标也比 Ez 低半格更新时统一偏移一位即可。使用一维连续数组存储时换算公式固定为index j * nx i读写速度比二维指针数组更快。2.2 二维 TMz 波的电场磁场更新公式真空中无损耗介质的旋度方程展开后得到三个标量方程∂Hx/∂t -μ⁻¹ ∂Ez/∂y∂Hy/∂t μ⁻¹ ∂Ez/∂x∂Ez/∂t ε⁻¹ (∂Hy/∂x - ∂Hx/∂y)用中心差分替换空间导数后C 语言更新函数可以写成下面这样。/* 均匀网格dx dy 时系数可以合并 */ const double c1 dt / (mu * dy); /* Hx 更新系数 */ const double c2 dt / (mu * dx); /* Hy 更新系数 */ const double c3 dt / (eps * dx); /* Ez 中 Hy 差分系数 */ const double c4 dt / (eps * dy); /* Ez 中 Hx 差分系数 */ void update_h(FDTD2D *g) { for (int j 1; j g-ny; j) { for (int i 1; i g-nx; i) { g-hx[j * g-nx i] c1 * ( g-ez[j * g-nx i] - g-ez[(j - 1) * g-nx i]); g-hy[j * g-nx i] c2 * ( g-ez[j * g-nx i] - g-ez[j * g-nx (i - 1)]); } } } void update_e(FDTD2D *g) { for (int j 1; j g-ny - 1; j) { for (int i 1; i g-nx - 1; i) { double dhy g-hy[j * g-nx (i 1)] - g-hy[j * g-nx i]; double dhx g-hx[(j 1) * g-nx i] - g-hx[j * g-nx i]; g-ez[j * g-nx i] c3 * dhy - c4 * dhx; } } }磁场更新依赖前后两个 Ez 的差值电场更新依赖两个 Hy 和一个 Hx 的差值。为了让边界上的点不越界循环索引从 1 开始网格最外一圈不在更新范围内而由吸收边界单独处理。细心的读者会发现 Hx 的更新式里没有负号Hy 的更新式也是正号。这取决于坐标轴方向和 Yee 网格的偏移方向。如果你从另一本教材里抄来的公式带负号不要急着改程序先统一坐标轴定义。FDTD 本身没有方向偏好只要符号与坐标一致最后结果就是相同的。2.3 CFL 条件dt 不能随意设显式 FDTD 的数值波速必须大于等于物理光速否则误差波会追上真实波导致结果发散。二维均匀网格的稳定性条件为dt 1 / (c * sqrt(1 / dx² 1 / dy²))在dx dy时这个上限简化为dt dx / (c * sqrt(2))。实际工程里我习惯取上限的 0.9留出安全余量。后面如果加入介质层光速变小CFL 上限会略微放宽但调试阶段不要一上来就用接近极限的值。先跑通真空模型再慢慢加介质。3. 用 C 语言把 FDTD_32 组织成完整源码FDTD 计算器的代码结构比算法本身更容易劝退人。我习惯用一个结构体把网格参数和三个场数组全部包起来所有更新函数只接收一个FDTD2D *指针。这样后续加 PML、加金属边界、加介质层都不需要改动主循环的函数签名。3.1 用结构体承载计算器参数先看核心结构体。typedef struct { int nx, ny; /* 网格尺寸 */ double dx, dy, dt; /* 空间步长和时间步长 */ double t; /* 当前仿真时间 */ double *ez, *hx, *hy; /* 场数组一维连续内存 */ } FDTD2D;选择一维数组而不是二维数组是为了减少malloc的次数并且让同一行的场数据在缓存里连续。nx、ny后面会频繁参与下标换算如果写死 32以后就要全文件搜索修改。把它们放进结构体后用指针传递是 C 语言项目的标准做法也锻炼了c语言结构体和c语言指针的配合能力。3.2 malloc 与 free网格数组的内存管理创建和销毁函数如下。FDTD2D *grid_create(int nx, int ny, double dx, double dy) { FDTD2D *g (FDTD2D *)malloc(sizeof(FDTD2D)); if (g NULL) { perror(malloc FDTD2D); exit(EXIT_FAILURE); } g-nx nx; g-ny ny; g-dx dx; g-dy dy; g-t 0.0; g-ez (double *)calloc((size_t)nx * ny, sizeof(double)); g-hx (double *)calloc((size_t)nx * ny, sizeof(double)); g-hy (double *)calloc((size_t)nx * ny, sizeof(double)); return g; } void grid_free(FDTD2D *g) { free(g-ez); free(g-hx); free(g-hy); free(g); }calloc会把整块内存清零比malloc加memset更顺手。数组长度用了(size_t)nx * ny避免在 32 位系统上做int乘法时溢出。创建和释放写成独立函数比在main里手工维护安全得多。这也是网上那批c语言必背100代码里最常见的题型构造函数配析构函数不让内存泄漏。3.3 波源设置与简单吸收边界波源通常使用高斯脉冲脉冲宽度由时间常数tau控制。void add_gaussian_source(FDTD2D *g, int sx, int sy, double t0, double tau, double amp) { double pulse amp * exp(-pow((g-t - t0) / tau, 2)); g-ez[sy * g-nx sx] pulse; }波源加在电场数组上物理上是一个电流源。如果加在磁场数组上就变成磁流源两者辐射方向场不同。学习程序里不建议直接放正弦源因为正弦源会立刻把边界打爆初学者根本分不清是边界反射还是算法发散。高斯脉冲从零平滑上升到最高点再降下来边界吸收的压力小很多。吸收边界我用简化版的海绵层不实现完整的 PML。把电场在网格边缘按线性因子缩小让外行波被慢慢吃掉。void apply_absorbing_taper(FDTD2D *g, int thickness) { for (int j 0; j g-ny; j) { for (int i 0; i g-nx; i) { double factor 1.0; if (i thickness) factor (double)i / thickness; else if (i g-nx - thickness) factor (double)(g-nx - 1 - i) / thickness; if (j thickness) factor (double)j / thickness; else if (j g-ny - thickness) factor (double)(g-ny - 1 - j) / thickness; g-ez[j * g-nx i] * factor; } } }这个海绵层会吸收大部分向外传播的能量但也会造成轻微反射。对 32×32 的演示网格足够真正研究回波损耗还是要上 PML。thickness一般取 4 到 8太小没有衰减效果太大则压缩有效计算域。3.4 把场快照写入 CSV 文件输出用纯文本 CSV方便后续交给 gnuplot 或 Python 出图。void save_ez_csv(FDTD2D *g, const char *path) { FILE *fp fopen(path, w); if (fp NULL) { perror(path); return; } for (int j 0; j g-ny; j) { for (int i 0; i g-nx; i) { fprintf(fp, %g%c, g-ez[j * g-nx i], (i g-nx - 1) ? \n : ,); } } fclose(fp); }fprintf使用%g而不是%f避免输出一堆没有物理意义的零。每个 CSV 文件的行数等于ny列数等于nx。文件读写操作里一定要检查fopen的返回值否则磁盘满或目录不存在时会出现静默失败。这一步踩过的坑和c语言文件读写操作代码的热搜里描述的问题基本一致。3.5 用 Makefile 一键编译运行一段最小的 Makefile 足够在 Linux、macOS 和 Windows 的 MSYS2 里使用。CC gcc CFLAGS -stdc99 -O2 -Wall LDLIBS -lm main: main.c $(CC) $(CFLAGS) -o main main.c $(LDLIBS) clean: rm -f main在 VS Code 配置 C 语言环境时tasks.json里把编译命令指定成make main再配置一个 build task就能一键编译。优化级别-O2很关键FDTD 是循环密集型程序-O0下同一套参数可能要慢两三倍。接下来是主循环的大致样子。int main(void) { FDTD2D *g grid_create(32, 32, 1e-3, 1e-3); double dt_max 1.0 / (C0 * sqrt(1.0 / (g-dx * g-dx) 1.0 / (g-dy * g-dy))); g-dt dt_max * 0.9; for (int n 0; n 500; n) { update_h(g); update_e(g); add_gaussian_source(g, 16, 16, 5e-11, 1e-11, 1.0); apply_absorbing_taper(g, 4); g-t g-dt; if (n % 50 0) save_ez_csv(g, field.csv); } grid_free(g); return 0; }主循环里电磁场交替更新波源在电场更新之后叠加。如果你把源叠加放在电场更新之前等价于把电流密度积分顺序换了一下波形上会出现半个时间步的偏移。先按这个顺序跑通再考虑优化。4. 参数整定与排错覆盖 CFL、网格尺寸和内存陷阱上一章的源码能跑但直接运行 32×32 网格如果参数不对屏幕上很快就会刷出 NaN。参数整定是 FDTD 从玩具变成实用工具的关键。4.1 常用参数表与推荐起点下面是一组能直接跑的起点参数。参数推荐值说明网格尺寸 nx, ny32 × 32默认值先跑通再放大空间步长 dx, dy1e-3 m每格 1 mm适合分米波演示时间步长 dt2.1e-12 s由 CFL 上限乘 0.9高斯源 t05 * tau让脉冲从零平滑上升高斯源 tau10 * dt覆盖 10 个左右时间步源位置16, 16网格正中间吸收层厚度432 网格下取 4 较合适空间步长决定了空间解析度。经验法则是每个波长至少要有 10 到 20 个网格。32×32 不适合分析精细结构但足够清楚展示波前传播。这样也能和二维波的解析解对比检查数值色散是否在合理范围内。4.2 一次完整的 CFL 计算例子以下代码展示了如何在main里计算安全dt。#define C0 3e8 double dx 1e-3, dy 1e-3; double dt_max 1.0 / (C0 * sqrt(1.0 / (dx * dx) 1.0 / (dy * dy))); double dt dt_max * 0.9; printf(dt_max %g s, dt %g s\n, dt_max, dt);把dx1e-3, dy1e-3代入1/dx²等于 1e6两项加起来是 2e6开根后约 1414。dt_max约等于 2.36e-12 秒取 0.9 倍后是 2.1e-12 秒。这个量级在编码时经常因为单位写错而差出十个数量级所以我初始化时会先打印一次dt再进入主循环。4.3 计算器不对时的检查顺序如果场值发散按下面顺序排查。检查dt是否满足 CFL。把dt调成上限的 0.5如果立刻正常问题就是时间步长太大。检查边界。apply_absorbing_taper必须在每次电场更新后调用漏掉的话能量会在边界积累。检查源的位置。源必须落在 1 到nx-2、1 到ny-2的范围内否则源点本身处于不更新区域能量无法进入求解域。检查数组类型。场数组必须是double*不要用int*整型数组会把场值截断成 0 或 1画出来不是波形而是噪点。用c语言的valgrind跑小步数定位越界读写。FDTD 数组下标换算最容易在循环边界差一拍常见错误是把i nx写成i nx这会读到ez数组后面的邻居内存。这五步覆盖大部分初学问题。如果还不对就回到最小模型去掉源只保留初始场看总能量是否逐步下降。这样能区分算法问题和数据问题。5. 验证技巧给 FDTD_32 加能量探针并用 gnuplot 成图最后一个技巧是给计算器加一个总能量探针。每次电场更新完后扫描全网格把电场和磁场能量加起来。double total_energy(FDTD2D *g) { double sum 0.0; for (int j 0; j g-ny; j) { for (int i 0; i g-nx; i) { double e g-ez[j * g-nx i]; double hx g-hx[j * g-nx i]; double hy g-hy[j * g-nx i]; sum 0.5 * e * e 0.5 * (hx * hx hy * hy); } } return sum; }在main循环里每 10 步写一行time energy到energy.dat。如果海绵吸收边界有效能量曲线会单调下降但不会出现突兀跳变如果指数上涨大概率是 CFL 被破坏。接下来用 gnuplot 直接看热力图。gnuplot -e set terminal png; set output ez.png; \ set view map; splot field.csv matrix with imagematrix关键字告诉 gnuplotCSV 是按网格行组织的二维数组。跑完后如果能看到圆形波前从中心扩散并在边缘被吸收层压平说明从公式到源码的整条链路是通的。把能量探针再扩展一步还能在固定点记录 Ez 随时间变化用离散傅里叶变换扫出幅度谱进而找谐振频率。FDTD_32 作为计算器最实用的下一步就是这个频域响应分析。本文还有配套的精品资源点击获取