简介面向数值分析与偏微分方程求解方向的学习者这份资源聚焦二维声波方程的间断有限元DG数值实现将空间离散与三阶龙格库塔时间推进相结合适合已具备Matlab基础、希望深入掌握高阶数值方法的中高级读者。压缩包共6个文件以5个m脚本和1个asv备份文件为主涵盖主程序、时间推进、基函数与初始条件等模块整体约4KB体量轻便便于逐行研读。已有1566人学习下载说明该实现路线在相关课程与科研练习中具有一定参考价值。读者可从中获得从网格划分、界面条件处理到时间步迭代的完整代码框架并借助可视化观察正弦初值下声波传播的演化过程为理解DG方法的稳定性与精度提供可运行的实践样本。1. 声波方程数值模拟为什么我最后选了间断有限元这套 MATLAB 实现第一次接触声波方程数值模拟很多人是从有限差分起步的——网格一画、差分格式一套跑起来确实快。但真到了复杂介质、强变速、甚至带裂缝的模型上有限差分那套规则网格就开始翻车界面处精度掉得厉害想加密网格又得全局加密算力直接爆炸。我后来转向间断有限元Discontinuous Galerkin简称 DG核心原因就一个——它把“局部高阶精度”和“复杂几何适应性”这两件本来矛盾的事捏到了一起。这套 MATLAB 实现把声波方程的一阶速度-应力方程组用 DG 空间离散、低存储显式时间推进网格单元之间靠数值通量耦合天然适合并行和自适应。它适合谁适合正在做地震波场模拟、声学成像、或者想拿一个干净代码骨架入门 DG 的研究生和工程开发——不是那种调个库就出图的工具而是你能真正拆开看每个矩阵怎么组装的教学级工程代码。2. 架构拆解从一阶方程组到时间推进的完整链路2.1 控制方程与变量布局声波方程在 DG 里通常不直接离散二阶形式而是拆成一阶速度-应力方程组这样每个单元内部只涉及一阶导数数值通量写起来干净。二维情况下% 一阶声波方程组二维 % dvx/dt (1/rho) * dp/dx % dvy/dt (1/rho) * dp/dy % dp/dt rho * c^2 * (dvx/dx dvy/dy) % 其中 vx, vy 为质点速度p 为声压rho 密度c 声速代码里把三个物理量都存成按单元组织的矩阵每个单元内是 N 阶多项式基函数的系数。常见做法是用 Legendre 多项式做基节点取 Gauss-Legendre 或 Gauss-Lobatto 点这样质量矩阵是对角的求逆就是逐点除法省掉大规模线性求解。变量布局上我一般把vx、vy、p各存一个[Np, Ne]的矩阵Np是单元内节点数Ne是单元总数这样后面通量计算可以整矩阵向量化避免 for 循环拖慢 MATLAB。2.2 单元刚度矩阵与数值通量DG 和连续有限元最大的区别就在通量项。单元内部用分部积分得到体积分项和边界通量项边界上的不连续量必须用数值通量粘起来。声波方程常用的是迎风型通量upwind flux按特征线方向取上游值% 迎风通量计算以 x 方向为例 % nx 为单元面外法向 x 分量Z rho*c 为阻抗 % 面上左右两侧状态分别为 (vxL,pL) 和 (vxR,pR) function [Fx_v, Fx_p] upwind_flux_x(vxL, pL, vxR, pR, nx, Z) % 特征变量分解 wL_p pL Z*vxL*nx; % 右行特征 wL_m pL - Z*vxL*nx; % 左行特征 wR_p pR Z*vxR*nx; wR_m pR - Z*vxR*nx; % 迎风取上游 wp 0.5*(wL_p wR_p) 0.5*(wL_p - wR_p); % 简化示意 wm 0.5*(wL_m wR_m) - 0.5*(wL_m - wR_m); p_flux 0.5*(wp wm); vx_flux 0.5*(wp - wm)/Z; Fx_v vx_flux * nx; Fx_p p_flux * nx; end这段是通量的骨架实际代码里会把左右状态从相邻单元索引里取出来按面循环或向量化拼装。参数说明nx是面法向分量Z是介质阻抗迎风通量的物理意义就是让信息只从上游传向下游保证格式稳定。这里最容易出错的是法向方向——相邻两个单元对同一个面的法向是相反的拼装全局通量时必须统一符号否则波场会出现莫名其妙的反射。2.3 时间推进与 CFL 条件空间离散完时间上用低存储 Runge-KuttaLSRK显式推进。DG 的稳定性受 CFL 条件约束时间步长和单元尺寸、多项式阶数、波速都相关% 低存储 RK 时间推进四阶 for stage 1:5 % 计算当前 stage 的右端项 rhs compute_rhs(vx, vy, p, mesh, mat); % LSRK 系数四阶五级 if stage 5 vx vx dt * a(stage) * rhs.vx; vy vy dt * a(stage) * rhs.vy; p p dt * a(stage) * rhs.p; end enddt的选取要满足dt CFL * h / (c_max * (2N1))其中h是单元特征尺寸N是多项式阶数。我一般先估一个保守的 CFL 系数0.1 到 0.2跑通后再往上调。这里有个血泪经验多项式阶数一高稳定域收得很快很多人按二阶格式的经验设 CFL结果跑几十步就发散还以为是通量写错了。2.4 目录结构与模块职责这套代码的目录组织一般长这样按职责分层方便二次开发时定位project/ ├── main.m % 主入口参数配置、初始化、时间循环 ├── mesh/ │ ├── build_mesh.m % 生成结构化/非结构化网格 │ └── connect.m % 建立单元邻接关系 ├── basis/ │ ├── legendre.m % Legendre 基函数与导数 │ └── nodes.m % Gauss 节点与权重 ├── solver/ │ ├── compute_rhs.m % 右端项体积分 通量 │ ├── upwind_flux.m % 数值通量 │ └── time_step.m % LSRK 推进 ├── medium/ │ └── model.m % 速度、密度模型赋值 └── post/ └── plot_wave.m % 波场快照与地震记录main.m只做编排不写物理solver/里是纯数值逻辑换方程只改这里medium/和post/是外围改模型和出图不动核心。这种分层是我拆过一堆科研代码后最认的做法——核心求解器和业务逻辑分开二次开发时不会牵一发动全身。3. 跑起来从零配置到第一张波场快照3.1 环境准备与参数配置MATLAB 版本建议 R2018b 以上用到的都是基础矩阵运算和sparse不需要额外工具箱。主脚本开头集中放参数我习惯把物理参数和数值参数分块% 物理参数 c0 3000; % 背景声速 (m/s) rho0 2500; % 背景密度 (kg/m^3) f0 15; % 震源主频 (Hz) % 数值参数 N 4; % 多项式阶数 CFL 0.15; % 稳定性系数 T 0.5; % 模拟总时长 (s) order 4; % 时间格式阶数N决定单元内节点数Np (N1)*(N2)/2二维三角形或(N1)^2四边形阶数越高单单元越贵但全局网格可以更粗。CFL先给保守值跑通再调。T和f0决定你能看到几个波长的传播别一上来就设很大先跑 0.2 秒看波前形状对不对。3.2 网格生成与邻接关系结构化四边形网格最容易上手build_mesh.m里按nx × ny划分记录每个单元的节点坐标和四个面的邻居索引% 生成结构化四边形网格 nx 80; ny 80; Lx 2000; Ly 2000; [xg, yg] meshgrid(linspace(0,Lx,nx1), linspace(0,Ly,ny1)); % 单元顶点坐标与邻居表 elem zeros(nx*ny, 4); neigh zeros(nx*ny, 4); % 四个面左、右、下、上 for j 1:ny for i 1:nx eid (j-1)*nx i; % 顶点索引、邻居索引填充边界置 0 表示无邻居 % ... end end邻接表是通量计算的关键neigh(eid, face)存对面单元编号边界单元存 0通量函数里遇到 0 就走边界条件自由表面或吸收边界。这里翻车的点在于索引顺序——MATLAB 列优先(j-1)*nxi和(i-1)*nyj混用会导致邻居错位波场会出现网格状伪影排查半天才发现是编号约定不统一。3.3 震源注入与边界处理震源一般用 Ricker 子波在某个空间点按时间加源项。边界处理有两种常见选择自由表面声压反射和吸收边界PML 或 sponge 层。教学代码里常用 sponge 层在边界几个单元内加阻尼% Ricker 子波 t (0:dt:T); src (1 - 2*pi^2*f0^2*(t-1/f0).^2) .* exp(-pi^2*f0^2*(t-1/f0).^2); % sponge 阻尼边界层内衰减 damp zeros(Ne,1); for e 1:Ne dist distance_to_boundary(e, mesh); if dist sponge_width damp(e) damp_max * (1 - dist/sponge_width)^2; end end % 时间推进中每步乘衰减因子 p p .* exp(-damp*dt);sponge_width一般取 5 到 10 个单元damp_max要试太小反射压不住太大边界附近波场失真。我一般先不加吸收边界跑一遍看自由表面反射是否符合预期确认内部格式没问题再加 sponge 调参数。3.4 运行与结果检查跑main.m输出波场快照和接收点地震记录。检查顺序建议先看单点源在均匀介质里的波前是不是标准圆形再看两层模型界面处是否有非物理反射最后看能量是否随时间单调衰减无源区域。如果波前出现菱形或网格状畸变八成是通量法向或邻接索引错了如果能量暴涨先查 CFL 和通量符号。这套检查流程是我每次改完求解器必走的比盯着代码看有效得多。4. 避坑与排查DG 声波模拟里最容易翻车的五件事4.1 波场出现网格状伪影现象波前不是光滑圆弧而是沿网格线出现锯齿或菱形纹路。原因通常是相邻单元法向符号不统一或者邻接表索引错位导致通量在面上不守恒。解决写一个单元测试构造两个相邻单元、给一个解析解检查通量项在面上左右残差是否抵消再打印几个面的neigh索引和法向人工核对方向。4.2 时间推进几十步后发散现象前几步正常之后数值指数增长。原因多为 CFL 超限尤其是高阶多项式时稳定域收缩。解决把CFL降到 0.05 重跑若稳定则逐步上调找到临界值同时确认dt计算里用了(2N1)因子漏掉它是最常见的低级错误。4.3 界面处出现非物理反射现象两层介质模型里除了真实界面反射还多出几道弱反射。原因可能是介质参数在单元内被平均了而 DG 要求参数在单元内常数或按基函数展开。解决把rho、c按单元赋值而不是按节点插值界面处让单元边界与介质界面重合别让一个单元跨两种介质。4.4 边界反射压不住现象加了 sponge 层边界仍有明显反射。原因阻尼层太薄或damp_max太小或者阻尼只加在声压上没加速度。解决加宽 sponge 到 10 个单元damp_max按c/h量级估三个变量同步加阻尼也可以换 PML但实现复杂度高教学阶段 sponge 够用。4.5 MATLAB 循环太慢现象网格一细跑一次要几十分钟。原因通量计算里用了三重 for 循环没向量化。解决把面状态用索引矩阵一次性取出通量整矩阵运算compute_rhs里避免动态扩容预分配rhs.vx等数组。向量化后通常能快一个数量级这是 MATLAB 写数值代码的基本功。5. 二次开发换方程、加并行、做验证的实操思路拿到这套代码最常见的改造方向有三个换物理方程、提升性能、做收敛性验证。换方程最直接——DG 框架是通用的把compute_rhs.m里的通量和体积分换成弹性波或电磁波的形式基函数和网格模块不用动。我一般先在均匀介质里跑解析解对比确认新方程的通量符号和特征速度写对再上复杂模型。性能上MATLAB 的瓶颈在通量循环除了向量化还可以把compute_rhs用mex编译或者把面循环改成arrayfun/pagefun但后者对老版本兼容性差我倾向先做纯矩阵化。收敛性验证是科研代码的后悔药固定 CFL把网格从粗到细跑几组算L2误差对h的斜率理论上是N1阶斜率对不上说明实现有隐藏 bug。下面这个验证脚本我每次改完核心都会跑% 收敛性验证均匀介质解析解对比 h_list [1/20, 1/40, 1/80, 1/160]; err zeros(size(h_list)); for k 1:length(h_list) mesh build_mesh(h_list(k)); [p_num, p_exact] run_case(mesh, N, CFL, T); err(k) sqrt(mean((p_num(:)-p_exact(:)).^2)); end % 拟合斜率 slope polyfit(log(h_list), log(err), 1); fprintf(收敛阶 %.2f (理论 %d)\n, -slope(1), N1);h_list至少四组run_case里保证除网格外所有参数一致p_exact用均匀介质的解析 Green 函数或已知行波解。斜率低于理论值一阶以上优先查通量实现和边界处理而不是怀疑基函数。从那以后我每次动求解器核心都强制先跑一遍这个收敛脚本再谈结果——它救过我至少三次有一次是迎风通量里一个0.5系数写成了1.0波场看着正常收敛阶直接掉到一阶肉眼根本发现不了。希望这套拆解能帮到你少走几段我踩过的弯路。本文还有配套的精品资源点击获取 
长河君 · 2026/10/29 6:02:47 · 阅读 1.2k 间断有限元DG声波方程数值模拟MATLAB实现与避坑指南
简介面向数值分析与偏微分方程求解方向的学习者这份资源聚焦二维声波方程的间断有限元DG数值实现将空间离散与三阶龙格库塔时间推进相结合适合已具备Matlab基础、希望深入掌握高阶数值方法的中高级读者。压缩包共6个文件以5个m脚本和1个asv备份文件为主涵盖主程序、时间推进、基函数与初始条件等模块整体约4KB体量轻便便于逐行研读。已有1566人学习下载说明该实现路线在相关课程与科研练习中具有一定参考价值。读者可从中获得从网格划分、界面条件处理到时间步迭代的完整代码框架并借助可视化观察正弦初值下声波传播的演化过程为理解DG方法的稳定性与精度提供可运行的实践样本。1. 声波方程数值模拟为什么我最后选了间断有限元这套 MATLAB 实现第一次接触声波方程数值模拟很多人是从有限差分起步的——网格一画、差分格式一套跑起来确实快。但真到了复杂介质、强变速、甚至带裂缝的模型上有限差分那套规则网格就开始翻车界面处精度掉得厉害想加密网格又得全局加密算力直接爆炸。我后来转向间断有限元Discontinuous Galerkin简称 DG核心原因就一个——它把“局部高阶精度”和“复杂几何适应性”这两件本来矛盾的事捏到了一起。这套 MATLAB 实现把声波方程的一阶速度-应力方程组用 DG 空间离散、低存储显式时间推进网格单元之间靠数值通量耦合天然适合并行和自适应。它适合谁适合正在做地震波场模拟、声学成像、或者想拿一个干净代码骨架入门 DG 的研究生和工程开发——不是那种调个库就出图的工具而是你能真正拆开看每个矩阵怎么组装的教学级工程代码。2. 架构拆解从一阶方程组到时间推进的完整链路2.1 控制方程与变量布局声波方程在 DG 里通常不直接离散二阶形式而是拆成一阶速度-应力方程组这样每个单元内部只涉及一阶导数数值通量写起来干净。二维情况下% 一阶声波方程组二维 % dvx/dt (1/rho) * dp/dx % dvy/dt (1/rho) * dp/dy % dp/dt rho * c^2 * (dvx/dx dvy/dy) % 其中 vx, vy 为质点速度p 为声压rho 密度c 声速代码里把三个物理量都存成按单元组织的矩阵每个单元内是 N 阶多项式基函数的系数。常见做法是用 Legendre 多项式做基节点取 Gauss-Legendre 或 Gauss-Lobatto 点这样质量矩阵是对角的求逆就是逐点除法省掉大规模线性求解。变量布局上我一般把vx、vy、p各存一个[Np, Ne]的矩阵Np是单元内节点数Ne是单元总数这样后面通量计算可以整矩阵向量化避免 for 循环拖慢 MATLAB。2.2 单元刚度矩阵与数值通量DG 和连续有限元最大的区别就在通量项。单元内部用分部积分得到体积分项和边界通量项边界上的不连续量必须用数值通量粘起来。声波方程常用的是迎风型通量upwind flux按特征线方向取上游值% 迎风通量计算以 x 方向为例 % nx 为单元面外法向 x 分量Z rho*c 为阻抗 % 面上左右两侧状态分别为 (vxL,pL) 和 (vxR,pR) function [Fx_v, Fx_p] upwind_flux_x(vxL, pL, vxR, pR, nx, Z) % 特征变量分解 wL_p pL Z*vxL*nx; % 右行特征 wL_m pL - Z*vxL*nx; % 左行特征 wR_p pR Z*vxR*nx; wR_m pR - Z*vxR*nx; % 迎风取上游 wp 0.5*(wL_p wR_p) 0.5*(wL_p - wR_p); % 简化示意 wm 0.5*(wL_m wR_m) - 0.5*(wL_m - wR_m); p_flux 0.5*(wp wm); vx_flux 0.5*(wp - wm)/Z; Fx_v vx_flux * nx; Fx_p p_flux * nx; end这段是通量的骨架实际代码里会把左右状态从相邻单元索引里取出来按面循环或向量化拼装。参数说明nx是面法向分量Z是介质阻抗迎风通量的物理意义就是让信息只从上游传向下游保证格式稳定。这里最容易出错的是法向方向——相邻两个单元对同一个面的法向是相反的拼装全局通量时必须统一符号否则波场会出现莫名其妙的反射。2.3 时间推进与 CFL 条件空间离散完时间上用低存储 Runge-KuttaLSRK显式推进。DG 的稳定性受 CFL 条件约束时间步长和单元尺寸、多项式阶数、波速都相关% 低存储 RK 时间推进四阶 for stage 1:5 % 计算当前 stage 的右端项 rhs compute_rhs(vx, vy, p, mesh, mat); % LSRK 系数四阶五级 if stage 5 vx vx dt * a(stage) * rhs.vx; vy vy dt * a(stage) * rhs.vy; p p dt * a(stage) * rhs.p; end enddt的选取要满足dt CFL * h / (c_max * (2N1))其中h是单元特征尺寸N是多项式阶数。我一般先估一个保守的 CFL 系数0.1 到 0.2跑通后再往上调。这里有个血泪经验多项式阶数一高稳定域收得很快很多人按二阶格式的经验设 CFL结果跑几十步就发散还以为是通量写错了。2.4 目录结构与模块职责这套代码的目录组织一般长这样按职责分层方便二次开发时定位project/ ├── main.m % 主入口参数配置、初始化、时间循环 ├── mesh/ │ ├── build_mesh.m % 生成结构化/非结构化网格 │ └── connect.m % 建立单元邻接关系 ├── basis/ │ ├── legendre.m % Legendre 基函数与导数 │ └── nodes.m % Gauss 节点与权重 ├── solver/ │ ├── compute_rhs.m % 右端项体积分 通量 │ ├── upwind_flux.m % 数值通量 │ └── time_step.m % LSRK 推进 ├── medium/ │ └── model.m % 速度、密度模型赋值 └── post/ └── plot_wave.m % 波场快照与地震记录main.m只做编排不写物理solver/里是纯数值逻辑换方程只改这里medium/和post/是外围改模型和出图不动核心。这种分层是我拆过一堆科研代码后最认的做法——核心求解器和业务逻辑分开二次开发时不会牵一发动全身。3. 跑起来从零配置到第一张波场快照3.1 环境准备与参数配置MATLAB 版本建议 R2018b 以上用到的都是基础矩阵运算和sparse不需要额外工具箱。主脚本开头集中放参数我习惯把物理参数和数值参数分块% 物理参数 c0 3000; % 背景声速 (m/s) rho0 2500; % 背景密度 (kg/m^3) f0 15; % 震源主频 (Hz) % 数值参数 N 4; % 多项式阶数 CFL 0.15; % 稳定性系数 T 0.5; % 模拟总时长 (s) order 4; % 时间格式阶数N决定单元内节点数Np (N1)*(N2)/2二维三角形或(N1)^2四边形阶数越高单单元越贵但全局网格可以更粗。CFL先给保守值跑通再调。T和f0决定你能看到几个波长的传播别一上来就设很大先跑 0.2 秒看波前形状对不对。3.2 网格生成与邻接关系结构化四边形网格最容易上手build_mesh.m里按nx × ny划分记录每个单元的节点坐标和四个面的邻居索引% 生成结构化四边形网格 nx 80; ny 80; Lx 2000; Ly 2000; [xg, yg] meshgrid(linspace(0,Lx,nx1), linspace(0,Ly,ny1)); % 单元顶点坐标与邻居表 elem zeros(nx*ny, 4); neigh zeros(nx*ny, 4); % 四个面左、右、下、上 for j 1:ny for i 1:nx eid (j-1)*nx i; % 顶点索引、邻居索引填充边界置 0 表示无邻居 % ... end end邻接表是通量计算的关键neigh(eid, face)存对面单元编号边界单元存 0通量函数里遇到 0 就走边界条件自由表面或吸收边界。这里翻车的点在于索引顺序——MATLAB 列优先(j-1)*nxi和(i-1)*nyj混用会导致邻居错位波场会出现网格状伪影排查半天才发现是编号约定不统一。3.3 震源注入与边界处理震源一般用 Ricker 子波在某个空间点按时间加源项。边界处理有两种常见选择自由表面声压反射和吸收边界PML 或 sponge 层。教学代码里常用 sponge 层在边界几个单元内加阻尼% Ricker 子波 t (0:dt:T); src (1 - 2*pi^2*f0^2*(t-1/f0).^2) .* exp(-pi^2*f0^2*(t-1/f0).^2); % sponge 阻尼边界层内衰减 damp zeros(Ne,1); for e 1:Ne dist distance_to_boundary(e, mesh); if dist sponge_width damp(e) damp_max * (1 - dist/sponge_width)^2; end end % 时间推进中每步乘衰减因子 p p .* exp(-damp*dt);sponge_width一般取 5 到 10 个单元damp_max要试太小反射压不住太大边界附近波场失真。我一般先不加吸收边界跑一遍看自由表面反射是否符合预期确认内部格式没问题再加 sponge 调参数。3.4 运行与结果检查跑main.m输出波场快照和接收点地震记录。检查顺序建议先看单点源在均匀介质里的波前是不是标准圆形再看两层模型界面处是否有非物理反射最后看能量是否随时间单调衰减无源区域。如果波前出现菱形或网格状畸变八成是通量法向或邻接索引错了如果能量暴涨先查 CFL 和通量符号。这套检查流程是我每次改完求解器必走的比盯着代码看有效得多。4. 避坑与排查DG 声波模拟里最容易翻车的五件事4.1 波场出现网格状伪影现象波前不是光滑圆弧而是沿网格线出现锯齿或菱形纹路。原因通常是相邻单元法向符号不统一或者邻接表索引错位导致通量在面上不守恒。解决写一个单元测试构造两个相邻单元、给一个解析解检查通量项在面上左右残差是否抵消再打印几个面的neigh索引和法向人工核对方向。4.2 时间推进几十步后发散现象前几步正常之后数值指数增长。原因多为 CFL 超限尤其是高阶多项式时稳定域收缩。解决把CFL降到 0.05 重跑若稳定则逐步上调找到临界值同时确认dt计算里用了(2N1)因子漏掉它是最常见的低级错误。4.3 界面处出现非物理反射现象两层介质模型里除了真实界面反射还多出几道弱反射。原因可能是介质参数在单元内被平均了而 DG 要求参数在单元内常数或按基函数展开。解决把rho、c按单元赋值而不是按节点插值界面处让单元边界与介质界面重合别让一个单元跨两种介质。4.4 边界反射压不住现象加了 sponge 层边界仍有明显反射。原因阻尼层太薄或damp_max太小或者阻尼只加在声压上没加速度。解决加宽 sponge 到 10 个单元damp_max按c/h量级估三个变量同步加阻尼也可以换 PML但实现复杂度高教学阶段 sponge 够用。4.5 MATLAB 循环太慢现象网格一细跑一次要几十分钟。原因通量计算里用了三重 for 循环没向量化。解决把面状态用索引矩阵一次性取出通量整矩阵运算compute_rhs里避免动态扩容预分配rhs.vx等数组。向量化后通常能快一个数量级这是 MATLAB 写数值代码的基本功。5. 二次开发换方程、加并行、做验证的实操思路拿到这套代码最常见的改造方向有三个换物理方程、提升性能、做收敛性验证。换方程最直接——DG 框架是通用的把compute_rhs.m里的通量和体积分换成弹性波或电磁波的形式基函数和网格模块不用动。我一般先在均匀介质里跑解析解对比确认新方程的通量符号和特征速度写对再上复杂模型。性能上MATLAB 的瓶颈在通量循环除了向量化还可以把compute_rhs用mex编译或者把面循环改成arrayfun/pagefun但后者对老版本兼容性差我倾向先做纯矩阵化。收敛性验证是科研代码的后悔药固定 CFL把网格从粗到细跑几组算L2误差对h的斜率理论上是N1阶斜率对不上说明实现有隐藏 bug。下面这个验证脚本我每次改完核心都会跑% 收敛性验证均匀介质解析解对比 h_list [1/20, 1/40, 1/80, 1/160]; err zeros(size(h_list)); for k 1:length(h_list) mesh build_mesh(h_list(k)); [p_num, p_exact] run_case(mesh, N, CFL, T); err(k) sqrt(mean((p_num(:)-p_exact(:)).^2)); end % 拟合斜率 slope polyfit(log(h_list), log(err), 1); fprintf(收敛阶 %.2f (理论 %d)\n, -slope(1), N1);h_list至少四组run_case里保证除网格外所有参数一致p_exact用均匀介质的解析 Green 函数或已知行波解。斜率低于理论值一阶以上优先查通量实现和边界处理而不是怀疑基函数。从那以后我每次动求解器核心都强制先跑一遍这个收敛脚本再谈结果——它救过我至少三次有一次是迎风通量里一个0.5系数写成了1.0波场看着正常收敛阶直接掉到一阶肉眼根本发现不了。希望这套拆解能帮到你少走几段我踩过的弯路。本文还有配套的精品资源点击获取