简介：面向数值分析与偏微分方程求解方向的学习者，这份资源聚焦二维声波方程的数值模拟，采用间断有限元（间断伽辽金）方法进行空间离散，并配合三阶龙格库塔方法完成时间推进，适合已具备Matlab基础、希望深入理解DG方法实现细节的中高级读者。压缩包共6个文件，以5个m脚本为主，另含1个asv备份文件，整体约4KB，脚本分别承担主程序调用、时间推进、基函数与初始条件设定等职责，结构紧凑便于逐文件研读。目前已有1566人学习下载。资源完整呈现了从网格划分、基函数选取、界面条件处理到时间步迭代的求解链路，读者可据此复现矩形区域内正弦初值下的声波传播过程，并结合可视化观察解的演化，是掌握高阶数值方法工程落地的实用参考。
1. 间断有限元：为什么声波方程偏偏要“断开”来算
连续有限元在求解声波方程时，最让人头疼的不是精度不够，而是它天生“怕间断”。当介质速度突变、边界反射剧烈、或者波前出现陡峭梯度时，连续基函数会在界面附近产生全局性的振荡污染，一个单元的误差能顺着网格爬到整个求解域。间断有限元（Discontinuous Galerkin，DG）反其道而行：它允许单元之间的解完全不连续，靠数值通量在界面上“谈判”交换信息。这个看似违背直觉的设计，恰恰是它能处理复杂声学场景的根本原因。间断伽辽金方法把“连续性”从强约束降级为弱耦合，换来的是局部紧致性、任意高阶精度和天然并行性。搞懂它为什么能 work，比会调一个求解器参数重要得多——因为声波方程里那些让连续方法翻车的坑，DG 从架构层面就绕过去了。
2. 弱形式与数值通量：DG 到底在“断”什么
2.1 从强形式到单元局部弱形式
声波方程的一阶速度-压力形式是理解 DG 的起点。把二阶波动方程拆成一阶系统，速度向量和压力标量分别满足：
∂p/∂t + ρc² ∇·v = 0 ∂v/∂t + (1/ρ) ∇p = 0连续有限元直接对这个系统做全局弱形式，要求试探函数在整域上属于同一个连续空间。DG 的做法是把求解域剖分成单元 K，在每个单元内部单独做弱形式，试探函数和检验函数只在单元内定义，跨单元没有任何连续性要求。对单元 K 乘检验函数 φ 并分部积分，得到：
∫_K (∂p/∂t) φ dV = ρc² ∫_K v·∇φ dV - ρc² ∫_∂K (v·n) φ dS注意最后那个边界积分项——它正是单元之间唯一的耦合通道。连续有限元里这一项会被相邻单元“抵消”掉，因为法向速度在界面上连续；DG 里不抵消，而是用一个数值通量替代真实的界面通量。这个替换是 DG 全部魔力的来源，也是它和连续方法最本质的分水岭。
常见的实现思路是采用迎风型通量（upwind flux）。对声波系统，界面上的数值通量取两侧状态的加权组合，权重由特征线方向决定。物理上可以理解为：波在界面上传播时，信息只从上游来，下游单元不该“提前知道”上游的状态。迎风通量把这种因果性编码进了离散格式，所以 DG 在时间推进时天然稳定，不需要额外加人工耗散。
2.2 数值通量为什么能保证稳定与守恒
数值通量的设计有两个硬指标：守恒性和稳定性。守恒性要求通量满足“反对称”条件——从单元 A 流向 B 的通量和从 B 流向 A 的通量大小相等方向相反，这样全局积分时内部界面上的贡献自动抵消，只有物理边界留下净通量。稳定性则要求通量是单调的，或者说满足所谓的“单侧 Lipschitz 条件”，保证离散能量不增长。
以压力通量为例，迎风格式写成：
p* = (p⁻ + p⁺)/2 + (Z/2)(v⁻·n - v⁺·n)其中 Z = ρc 是声阻抗，上标负号表示单元内侧、正号表示相邻单元侧。第一项是中心平均，第二项是阻抗加权的速度跳跃惩罚。这个结构不是拍脑袋凑的：把通量代入离散能量方程，跳跃项贡献一个负定的耗散项，耗散量正比于界面两侧解的跳跃平方。换句话说，DG 的稳定性不是靠调参数调出来的，而是数值通量的代数结构直接保证的。这也是为什么 DG 在长时间声波传播模拟中不容易出现能量漂移——耗散是格式内建的，不是事后补的。
提示：数值通量的选择直接决定格式的耗散特性。迎风通量耗散适中、精度高，适合大多数声学问题；Lax-Friedrichs 通量更耗散但更鲁棒，适合强非线性或激波类场景。声波方程是线性的，迎风通量通常是首选。
2.3 基函数与质量矩阵：为什么 DG 能做成显式
DG 的另一个关键设计是基函数的选择。在每个单元内用正交多项式基（比如 Legendre 多项式），质量矩阵就变成对角或块对角。这意味着时间推进时不需要求解全局线性系统——每个单元的质量矩阵求逆是局部的、小规模的，可以直接显式更新。
对比连续有限元：全局质量矩阵是稀疏但非对角的，显式时间推进需要解一个大规模线性系统，或者做质量集中（lumping）近似，后者会损失精度。DG 用正交基把质量矩阵对角化，代价是引入了界面通量计算，但换来了真正的显式时间步进。这个取舍在声波模拟里非常划算，因为声波问题的时间尺度通常要求小步长，显式格式的并行效率远高于每步解线性系统。
常见的实现思路是在参考单元上用 Gauss-Legendre 节点做数值积分，配合节点型基函数，这样质量矩阵严格对角，通量计算也只需要节点上的值，不需要做额外的投影。这套组合在 MATLAB 里实现起来尤其干净——矩阵运算天然适合向量化，单元循环可以压成一次批量操作。
3. 声波方程的 DG 离散：从半离散到时间推进
3.1 半离散系统的矩阵组装
把单元局部弱形式里的各项用基函数展开，得到每个单元的半离散方程：
M_K du_K/dt = S_K u_K - F_K(u_K, u_neighbor)M_K 是单元质量矩阵，S_K 是刚度矩阵（来自体积分项），F_K 是界面通量项。由于基函数只在单元内有定义，M_K 和 S_K 都是小尺寸稠密矩阵，尺寸等于单元内基函数个数。全局组装时不需要把 M_K 拼成一个大矩阵——每个单元独立存储、独立求逆，这正是 DG 局部性的体现。
在 MATLAB 里，一个自然的实现是把所有单元的 M_K 堆成一个三维数组，用 pageinv 或者循环做批量求逆。刚度矩阵和通量算子的组装则依赖单元几何映射：从参考单元到物理单元的 Jacobian 矩阵参与积分变换，声速和密度在单元内可以是常数，也可以按节点插值。如果介质非均匀，通量计算时界面两侧的阻抗要分别取各自单元的值，不能简单平均——这是非均匀介质 DG 实现里最容易翻车的地方之一。
3.2 时间积分：显式 RK 与 CFL 约束
半离散系统是常微分方程组，时间推进通常用显式 Runge-Kutta 格式。经典的四阶 RK（RK4）在精度和稳定性之间平衡得不错，但 DG 的稳定域对 CFL 数有严格要求。CFL 条件大致是：
Δt ≤ C * h / ( (2N+1) * c_max )其中 h 是单元特征尺寸，N 是多项式阶数，c_max 是最大声速。注意分母里的 (2N+1)——阶数越高，允许的时间步长越小。这是 DG 的一个固有代价：高阶精度换来了更紧的稳定性约束。实际实现时，CFL 常数 C 取决于通量类型和 RK 格式，迎风通量配 RK4 通常取 0.1 到 0.3 之间比较安全。
一个血泪经验是：不要用理论 CFL 上限直接跑长时模拟。数值实验里，理论边界附近往往已经出现缓慢的能量增长，短时间看不出来，跑几千步之后波场就糊了。稳妥的做法是先做收敛性测试，把 Δt 减半看解是否收敛，再留 20% 到 30% 的余量。
3.3 边界条件的弱施加
声波问题的边界条件（硬壁、吸收、自由表面）在 DG 里不是强加在节点上的，而是通过边界通量弱施加。以硬壁边界为例，法向速度为零，实现时把边界外侧的虚拟状态设为内侧状态的镜像：压力相同、法向速度反号。代入迎风通量，速度跳跃项自动为零，压力通量退化为内侧压力，物理上正好对应全反射。
吸收边界（比如一阶 Engquist-Majda 条件）则需要在边界通量里引入阻抗匹配项，让出射波不被反射。这里有个容易忽略的细节：吸收条件的系数依赖入射角，垂直入射时完美吸收，斜入射时仍有残余反射。如果模拟区域靠近边界有强源，残余反射会污染近场解。常见的补救是加一层海绵层（sponge layer），在边界附近逐渐增大阻尼，把残余反射压下去。海绵层的厚度和阻尼剖面需要调，太薄效果不够，太厚又增加计算量。
4. 避坑与排查：DG 声波模拟里那些让人后悔药的细节
4.1 通量符号写反导致能量爆炸
现象：跑几十步之后波场幅值指数增长，CFL 明明取得很小。
原因：迎风通量里速度跳跃项的符号写反了，耗散项变成负定，格式从稳定变成反稳定。DG 的稳定性对通量符号极其敏感，一个正负号错了，整个格式的性质就翻转。
解决：先用均匀介质、单频源做短时测试，监控离散能量随时间的变化。正常情况能量应该缓慢衰减或守恒，如果看到增长，逐项检查通量表达式。建议把通量计算单独写成一个函数，输入两侧状态和法向，输出通量值，方便单元测试。
4.2 非均匀介质界面阻抗取错
现象：两种介质交界面上出现非物理的强反射，反射系数和理论值对不上。
原因：界面通量计算时，两侧单元的声阻抗没有分别取各自的值，而是用了全局平均或者单侧值。声波在阻抗突变界面上的反射系数是 (Z2-Z1)/(Z2+Z1)，通量里如果阻抗取错，等效于人为改变了界面匹配条件。
解决：通量函数必须同时接收左右单元的 ρ 和 c，分别算 Z_left 和 Z_right。界面法向统一指向一侧，另一侧取反。测试时用两层介质、平面波垂直入射，对比解析反射系数，误差应该在离散精度量级。
4.3 高阶基函数下积分精度不足
现象：提高多项式阶数后精度不升反降，或者出现高频振荡。
原因：体积分和界面通量的数值积分用了低阶求积规则，高阶基函数的非线性项（比如变系数项）积分欠采样，产生混叠误差。DG 的精度依赖积分精确度，阶数越高对求积规则要求越高。
解决：体积分用至少 (N+1) 点 Gauss 求积，界面通量用 (N+1) 点或更高。如果介质参数在单元内变化剧烈，考虑过积分（over-integration）或者把参数也做多项式展开。MATLAB 里可以用 gaussquad 类函数生成节点和权重，别手写低阶梯形规则凑合。
4.4 时间步长只看 CFL 忽略源项刚度
现象：加了点源或者阻尼项之后，原本稳定的时间步长突然发散。
原因：CFL 条件只约束了波传播的稳定性，源项或阻尼项引入了额外的刚度。如果源的时间变化率很高，或者阻尼系数很大，显式 RK 的稳定域会被源项撑破。
解决：把源项和阻尼项纳入稳定性分析，估算它们的特征时间尺度，取更严格的步长。实在不行就把源项做隐式处理（IMEX 格式），波传播部分显式、源项部分隐式，兼顾效率和稳定。
4.5 并行组装时的索引错位
现象：单核跑正常，多核或者向量化之后结果对不上，误差集中在单元边界。
原因：DG 的界面通量需要访问相邻单元的数据，如果单元编号和邻居索引的映射在并行或向量化时错位，通量就用错了邻居状态。MATLAB 的向量化写法里，这种错位往往表现为隐式的维度广播错误，不报错但结果错。
解决：组装前显式建立邻居表，每个界面存左右单元编号和局部面编号。向量化时用这个表做索引，不要依赖单元数组的隐式顺序。测试时用两个单元的最小算例，打印每个界面的通量输入输出，逐一手工核对。
5. 设计取舍：DG 的代价与它换来的东西
DG 不是没有代价的。相比连续有限元，它的未知量数量翻倍（每个界面两侧都有独立自由度），通量计算增加了额外开销，高阶情况下时间步长更紧。这些代价在声波模拟里是否值得，取决于问题的性质。如果求解域规则、介质均匀、边界简单，连续有限元配合质量集中可能更经济；但只要出现介质突变、复杂几何、或者需要局部自适应，DG 的局部紧致性就变成压倒性优势——加密某个单元不会牵动全局自由度，误差不会跨界面传播，并行时单元之间只需要交换界面通量，通信量可控。
我一般会这么理解 DG 的设计哲学：它把“连续性”从一个全局硬约束拆解成一组局部弱约束，用数值通量做仲裁。这个拆解让格式获得了模块化和局部性，代价是引入了通量设计的自由度——通量选得好，稳定性和精度都有保证；选得不好，格式可能比连续方法还脆弱。所以 DG 的实现质量，很大程度上取决于对通量机制的理解深度，而不是代码写得多花哨。
验证一个 DG 声波求解器是否靠谱，我习惯走三步：先用均匀介质平面波做收敛性测试，看误差随网格加密的下降阶数是否匹配多项式阶数；再用两层介质界面测反射系数，对比解析解；最后跑一个长时能量监控，确认没有缓慢漂移。这三步过了，才敢往复杂场景上放。从那以后我每次改通量或者换基函数，都强制走一遍这套验证流程，省下来的调试时间远比写测试花的多。希望帮到你。
本文还有配套的精品资源，点击获取
