简介：面向数值分析与偏微分方程求解方向的学习者，这份资源聚焦二维声波方程的数值模拟，采用间断有限元（间断伽辽金）方法完成空间离散，并配合三阶龙格库塔格式推进时间步，适合已具备Matlab基础、希望深入理解DG方法实现细节的高年级本科生与研究生参考。压缩包共6个文件，以5个m脚本和1个asv备份文件为主，涵盖主程序、时间推进、初值设定与辅助函数等模块，整体约4KB，体量轻便、便于逐行阅读与二次修改。已有1567人学习下载，说明该实现思路在同类问题中具有一定参考价值。读者可从中获得从网格划分、基函数选取、界面条件处理到边界条件设置的完整代码框架，并借助可视化观察正弦初值下声波在矩形域内的传播演化，为掌握高阶数值方法提供可运行的实践样本。
1. 间断有限元：为什么声波方程偏偏要“断开”来算
连续有限元给人的直觉是网格越连续越精确，可一旦拿它去解声波方程，工程师很快会撞上一个尴尬现象：网格稍微加密，波形前沿就开始出现非物理振荡，长时间推进后能量莫名其妙地漂移。间断有限元（Discontinuous Galerkin，简称 DG）反其道而行，允许相邻单元的函数值在界面上“断开”，靠数值通量把两边缝起来。这个看似退步的选择，恰恰是它能同时拿到高阶精度、局部守恒和并行友好的根源。搞懂 DG 为什么能 work，比会调一个求解器重要得多——它解释了高阶波动模拟里那些玄学般的稳定性问题到底从哪来。这份围绕 MATLAB 实现声波方程间断有限元的资源，正好把这条链路从弱形式一路铺到时间推进，值得逐层拆开看。
2. 弱形式与数值通量：界面“断开”之后靠什么缝合
2.1 从强形式到单元局部弱形式
声波方程的一阶形式通常写成一阶双曲系统：压力对时间的导数等于速度散度乘以体积模量，速度对时间的导数等于压力梯度除以密度。DG 不去全局组装一个大矩阵，而是在每个单元内部独立做分部积分，把二阶空间导数降成一阶，代价是边界上多出一个通量项。这个通量项就是整个方法的命门：单元内部完全解耦，单元之间只通过界面通量交换信息。
常见的实现思路是先把方程乘上测试函数，在单元上积分，再对散度项做分部积分，得到含边界积分的弱形式。由于测试函数在单元外为零，边界积分只在单元界面上出现，而界面两侧的函数值一般不相等，于是必须人为定义一个“唯一”的界面通量来替代物理通量。这一步决定了格式的守恒性和稳定性，也是 DG 与连续有限元最本质的分野。
2.2 数值通量为什么必须“迎风”
界面两侧各有一套解，物理上只有一个真实通量，数值上必须选一个规则把两者合成一个。中心通量看似对称优雅，但它对波动问题不引入任何耗散，高频误差会不断累积，长时间推进必然翻车。工程上普遍采用迎风型通量（如 Lax-Friedrichs 或 Roe 型），按特征速度方向给两侧加权，本质是往格式里注入恰好够用的数值耗散，压住非物理振荡。
% 一维声波系统界面通量（Lax-Friedrichs 型） % qL, qR 为界面左右两侧的守恒量 [p; u] % A 为通量雅可比矩阵，alpha 为最大特征速度上界 function F = numerical_flux(qL, qR, A, alpha) FL = A * qL; % 左侧物理通量 FR = A * qR; % 右侧物理通量 F = 0.5*(FL + FR) ... % 中心部分 - 0.5*alpha*(qR - qL); % 迎风耗散项，alpha 需覆盖谱半径 end这段代码里alpha的取值是关键参数：它必须不小于局部雅可比矩阵谱半径的上界，取小了耗散不足、振荡压不住，取大了耗散过量、波形被抹平。实际实现中常按单元局部声速动态估计，而不是全局取一个常数，否则在非均匀介质里要么局部失稳要么过度耗散。通量函数返回后，会被乘上界面长度累加回相邻两个单元的右端项，这一步保证了“断开”的单元之间仍然严格守恒——一侧流出多少，另一侧就流入多少。
2.3 质量矩阵的块对角结构带来的红利
连续有限元组装出的是全局稀疏但耦合的刚度矩阵，而 DG 的单元质量矩阵是块对角的：每个单元一块小矩阵，块与块之间没有耦合。这意味着求时间导数时可以逐单元做小规模求逆，甚至预先算出逆矩阵反复使用。对显式时间推进来说，这是巨大的性能优势，也是 DG 天然适合并行和 GPU 加速的结构性原因。代价是自由度总数比连续有限元多（界面上的值重复存储），但换来的局部性和高阶精度通常更划算。
3. 基函数与高阶精度：单元内部怎么逼近才不失真
3.1 节点基与模态基的取舍
单元内部的多项式逼近有两种主流选择：节点型基函数（如 Lagrange 插值多项式，系数就是节点上的物理值）和模态型基函数（如正交多项式，系数是模态分量）。节点基直观、边界值直接可取，做通量计算方便；模态基在质量矩阵上更接近对角，正交性好，做投影和滤波更自然。声波问题里界面通量频繁用到边界值，节点基往往更顺手，但高阶时节点分布必须选得好。
等距节点在高阶下会触发 Runge 现象，插值边缘剧烈振荡，直接毁掉精度。血泪经验是：阶数超过三四阶就改用 Gauss-Legendre 或 Chebyshev-Lobatto 节点，把节点向单元两端聚集，压低 Lebesgue 常数。这一步不做，后面再怎么调通量都救不回来。
3.2 数值积分精度与欠积分隐患
弱形式里的积分一般用 Gauss 求积近似。多项式阶数为 p 时，质量矩阵被积函数是 2p 次，刚度项含通量后次数更高，求积点数不足会引入“欠积分”。欠积分不总是坏事——某些格式故意欠积分来引入耗散稳定高频——但对声波这类要求低耗散的问题，欠积分会让高阶模态悄悄失稳。
% 单元内 Gauss-Legendre 求积配置（阶数 p，点数 p+1 起步） p = 4; [xi, w] = gauss_legendre(p+1); % xi 为参考单元节点，w 为权重 % 质量矩阵 M_ij = sum_k w_k * phi_i(xi_k) * phi_j(xi_k) * J % J 为单元雅可比（物理长度/参考长度）参数说明：p+1个求积点对 2p 次被积函数是精确的，但含通量的项次数更高，实践中常取p+2甚至p+3个点留余量。J是几何映射的雅可比，非均匀网格下逐单元不同，必须在组装时带上，否则守恒性会被几何误差破坏。我一般会这么理解：求积点数是精度和稳定之间的一个旋钮，宁可多一两个点，也别在长时程模拟里赌欠积分不出事。
3.3 高阶带来的时间步限制
DG 的显式时间步受 CFL 条件约束，而 CFL 数随多项式阶数升高而收紧，大致按 p 的平方量级恶化。这意味着追求高阶精度时，时间步会变得很小，计算量陡增。常见的应对是配对高阶时间积分（如低存储 Runge-Kutta）或局部时间步进，让不同尺寸的单元各走各的步长。这个取舍在后面的稳定性章节还会再碰到。
4. 时间推进与稳定性：显式格式为什么总在长时程翻车
4.1 显式 Runge-Kutta 与 CFL 的硬约束
空间离散完成后，半离散系统是一组常微分方程，最直接的做法是显式 Runge-Kutta 推进。显式格式每步只做矩阵向量乘，无需解线性系统，配合 DG 的块对角质量矩阵极其高效。但显式格式的条件稳定性是硬约束：时间步必须小于某个由最大特征速度和最小单元尺寸决定的上界，否则误差指数放大。
% 低存储 RK 推进一阶系统 dq/dt = L(q) dt = cfl * hmin / (alpha_max * (2*p+1)); % 经验 CFL 缩放 q = q0; for n = 1:nt k1 = L(q); k2 = L(q + dt*k1); q = q + 0.5*dt*(k1 + k2); % 二阶 RK，示意 endcfl一般取 0.1 到 0.3 之间，2*p+1是对高阶格式的保守缩放因子，hmin是最小单元特征尺寸。这个公式不是精确理论界，而是工程上反复验证过的经验缩放，实际还要按具体通量和基函数微调。翻车往往就出在这里：网格里混进一个畸变的小单元，hmin骤降，全局时间步被拖垮，或者反过来 CFL 取太大，短波模态先失稳。
4.2 色散与耗散误差的物理表现
数值格式对不同的波数响应不同，短波（高波数）成分的相速度会偏离真实值，这就是数值色散；振幅被压低则是数值耗散。声波模拟里色散表现为波形后拖的伪振荡，耗散表现为脉冲被抹宽。DG 的高阶特性让它在解析良好的波数范围内色散极小，但接近网格分辨率极限的短波仍然会失真。理解这一点，才能判断结果是物理还是数值伪影。
4.3 边界条件怎么加才不反射
外边界处理是另一个易翻车点。简单地把边界值置零会引入强反射，污染整个计算域。常用做法是特征边界条件：按出入射特征分解，只对入射特征赋物理值，出射特征让它自由流出。实现上就是在边界界面上套用与内部一致的迎风通量逻辑，把外部状态设为满足无反射条件的状态。这一步做对了，开放边界的伪反射能压到很低。
5. 避坑与排查：DG 声波模拟最常见的几类翻车
5.1 现象：波形前沿出现高频毛刺，越跑越密
原因多半是通量耗散不足或 CFL 偏大，短波模态没被压住。先检查alpha是否覆盖了局部最大特征速度，再确认 CFL 是否在经验范围内。解决方式是适度增大迎风耗散或减小时间步，同时排查网格是否有畸变单元拉低了hmin。
5.2 现象：能量随时间单调增长，最终溢出
典型是欠积分或边界条件实现错误导致能量注入。检查求积点数是否足够、边界通量是否把出射特征也强行赋了值。用能量范数逐时刻打印是最快的定位手段，能量该守恒（无源无耗）时却漂移，问题一定在通量或边界。
5.3 现象：高阶精度没体现，误差随阶数不降反升
多半是节点分布等距引发 Runge 现象，或几何雅可比算错。换成 Gauss-Legendre 节点，逐单元核对J的符号和量级，非均匀网格下尤其容易把雅可比方向搞反。
5.4 现象：并行后结果和串行对不上
DG 单元解耦，并行边界上的通量需要两侧单元各算一次并保证一致。若两侧用了不同的alpha估计或不同的状态插值，通量就不守恒。统一通量计算函数、在界面两侧用同一份数据，是避免这类不一致的基本纪律。
6. 设计取舍：局部性与精度的平衡点在哪
DG 的全部设计都围绕一个核心张力：单元越独立，并行和守恒越好，但界面通量和自由度冗余的代价越高。选高阶还是多单元、选显式还是隐式、选节点基还是模态基，本质上都是在“局部性红利”和“精度成本”之间找平衡。高阶少单元能压低色散，但时间步被 CFL 拖累；低阶多单元时间步宽松，但界面通量计算量上升、色散变差。没有普适最优，只有针对具体波数范围和计算规模的最优。
验证一套 DG 声波实现是否可信，我习惯做三件事：一是收敛性测试，固定网格加密阶数、固定阶数加密网格，看误差是否按理论阶下降；二是能量守恒测试，无源无耗场景下能量漂移应随分辨率提高而收敛；三是单脉冲传播测试，看波形是否只有可控的色散而无异常增长。这三关过了，才敢往复杂介质和长时程上推。
从那以后我每次搭 DG 求解器，都强制先跑一遍单脉冲加能量监控，再谈精度调优——跳过这步的代价，往往是几天后才发现结果早就不可信。希望帮到你。
本文还有配套的精品资源，点击获取
