简介：面向数值分析与偏微分方程求解方向的学习者，这份资源聚焦二维声波方程的间断有限元（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-Kutta（LSRK）显式推进。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 * (2N+1))，其中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 = (N+1)*(N+2)/2（二维三角形）或(N+1)^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,nx+1), linspace(0,Ly,ny+1)); % 单元顶点坐标与邻居表 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)*nx+i和(i-1)*ny+j混用会导致邻居错位，波场会出现网格状伪影，排查半天才发现是编号约定不统一。
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计算里用了(2N+1)因子，漏掉它是最常见的低级错误。
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的斜率，理论上是N+1阶，斜率对不上说明实现有隐藏 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), N+1);h_list至少四组，run_case里保证除网格外所有参数一致，p_exact用均匀介质的解析 Green 函数或已知行波解。斜率低于理论值一阶以上，优先查通量实现和边界处理，而不是怀疑基函数。从那以后我每次动求解器核心，都强制先跑一遍这个收敛脚本再谈结果——它救过我至少三次，有一次是迎风通量里一个0.5系数写成了1.0，波场看着正常，收敛阶直接掉到一阶，肉眼根本发现不了。希望这套拆解能帮到你，少走几段我踩过的弯路。
本文还有配套的精品资源，点击获取
