简介：一套基于Matlab的二维声波方程间断有限元（DG）求解实现，面向数值计算、偏微分方程数值解方向的研究生与工程师。资源采用DG方法进行空间离散，并以三阶龙格库塔格式推进时间积分，覆盖网格划分、线性基函数构造、界面通量处理、边界条件设置与迭代求解等完整流程。包体共6个文件，包含5个.m脚本（主程序、初始条件、通量计算、数值通量函数等）和1个.asv备份文件，压缩包仅4KB，代码量精简，便于逐行理解与二次开发。已有1564人学习下载，适合希望掌握DG方法基础实现、或需要参考Matlab编程框架来开展波动方程模拟的读者。通过阅读源码可直观对照二维声波方程从初始条件到数值解的完整过程，并能基于现有框架修改初边值与基函数，扩展至更复杂的计算场景。
1. 间断有限元求解声波方程：先想清楚这资源到底替你解决了什么
手里同时拿到一套间断有限元求解声波方程的Matlab程序和一段教学视频时，大多数人的第一反应是先把波场跑起来看看。我的建议正相反：先在纸上把间断伽辽金（DG）方法里几个核心量写明白——多项式阶数 N、单元数 K、时间步长 dt、数值通量到底是什么——再动手指。不然程序爆了，你连该往哪个函数里查都不知道。
这套资源解决的是声波方程正演模拟里的核心问题：用间断有限元把一阶双曲型声波方程离散，每个单元内部用多项式近似，单元之间允许解不连续，靠数值通量交换信息，时间上显式推进，最终输出压力场和速度场的动态波场。它适合两类人：一类是刚把连续有限元编程跑通、想往 DG 跳的研究生；另一类是做波场模拟但被网格连续性和边界反射折磨的从业者。下载之后拿到的不是一张说明书，而是一条能把每个矩阵、每个通量函数都拆开看的主线。这篇笔记就按这条主线拆一遍，顺便把参数怎么设、坑在哪讲透。
2. 间断伽辽金方法：先弄懂DG在算什么，声波方程才不白解
在点开 demo 之前，先花二十分钟把 DG 的离散思路过一遍。这一步省不掉，因为你后面看到的每个函数、每个矩阵，全是从这个思路里长出来的。
2.1 声波方程的两种写法：DG为什么偏好一阶双曲系统
声波方程最常见的写法是二阶标量方程：
[ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} ]
这个形式有限元也能做，但要求相邻单元之间解至少 C0 连续，否则强形式里对空间两次求导根本没有定义。这正是连续伽辽金方法的边界所在：它逼着你在单元交界面上做连续性约束。
DG 的做法完全不同：单元内部随便用多项式，单元之间允许跳变，连续性根本不强加。这么一来就不能再用二阶方程了，不然裂缝处的通量没法定义。标准做法是把声波方程降成一阶双曲系统：
[ \frac{\partial p}{\partial t} + \rho c^2 \frac{\partial v}{\partial x} = 0 ]
[ \frac{\partial v}{\partial t} + \frac{1}{\rho}\frac{\partial p}{\partial x} = 0 ]
其中 p 是压力扰动，v 是质点速度，ρ 是密度，c 是声速。写成向量形式就是：
[ \frac{\partial Q}{\partial t} + A \frac{\partial Q}{\partial x} = 0, \quad Q = \begin{bmatrix}p\v\end{bmatrix}, \quad A = \begin{bmatrix}0 & \rho c^2\1/\rho & 0\end{bmatrix} ]
这个 A 叫通量雅可比矩阵，它的特征值是 ±c，对应左行波和右行波。这就是一维声波的 Riemann 问题的全部基础，也是数值通量设计的出发点。为什么必须写成这个形式？因为只有一阶双曲系统里，单元交界面上才能定义数值通量：左侧单元给一个状态，右侧单元给一个状态，中间用什么规则合成一个通量，这就是 DG 的命门。二阶方程写法的命门是连续性条件，一阶系统的命门是通量算法，两条技术路线分道扬镳的起点就在这里。
2.2 单元自由度与基函数：N、K、网格三个数先设对
DG 网格和连续有限元网格最大的区别是自由度编号方式。连续有限元里相邻单元共享节点，全局自由度编号需要查邻居表；DG 里每个单元独立拥有自己的一组节点，自由度数等于 K×(N+1)，单元之间不共享任何自由度。这段代码把这个结构写清楚：
% 单元数与多项式阶数 K = 100; % 单元个数 N = 3; % 单元内多项式次数 % 物理网格端点 xL = 0; xR = 1; xe = linspace(xL, xR, K+1); % 参考单元节点，范围 [-1,1]，共 N+1 个点 r = lobatto_points(N); % 常见做法：取 Legendre-Gauss-Lobatto 的根 % 预分配网格结构体 mesh = struct('J', zeros(K,1), 'x', zeros(K,N+1)); for e = 1:K mesh.J(e) = (xe(e+1) - xe(e)) / 2; % 雅可比因子 mesh.x(e,:) = (xe(e) + xe(e+1)) / 2 + mesh.J(e) * r; end这里 lobatto_points 是自己写的一个小函数，作用是给出 [-1,1] 上的 N+1 个插值节点。如果你用的是一套现成资源包，它内部一般已经封装了这个函数，你只需确认它返回的是 Gauss-Lobatto 点还是 Legendre 点。两者差别主要在质量矩阵是否严格对角，下文第 4 章会展开。现在你只需要记住三个参数：K 决定空间分辨率的下限，N 决定每个单元内部的逼近能力，J 是仿射变换的雅可比因子，它把物理单元映射到参考单元，所有积分都在参考单元上做。
第 2 章里这些代码先不参与求解，只是把数据结构铺好。真正进入迭代循环之前，你还需要理解数值通量怎么算、时间步长怎么限制，这是下一章的事。
3. 把声波方程写成DG格式：数值通量与主循环的参数经验
整个 DG 程序的核心不在地图网格生成里，而在两个函数：半离散的右端项 dg_rhs 和时间推进循环。前者负责把波场在每个单元内部推进半步，后者负责用合适的步长把时间往前走。这两个函数配对正确，波场动画基本就出来了。
3.1 半离散DG格式：单元内部弱形式，边界数值通量
对一维声波方程做 DG 离散，每个单元上的半离散格式长这样：
[ M_e \frac{d Q_e}{dt} = S_e (A Q_e) - F^* \phi \bigg|{x_R} + F^* \phi \bigg|{x_L} ]
逐项拆开看：M_e 是单元质量矩阵，S_e 是单元刚度矩阵，A 就是 2.1 节里那个通量雅可比矩阵，F^* 是数值通量，φ 是边界处的试探函数值。右边第一项是单元内部的通量贡献，后两项是左右两个交界面上数值通量的贡献，符号差异来自外法向方向：左边界外法向沿 −x，右边界外法向沿 +x。
关键点在 F^* 怎么取。如果直接把左右单元的解代入物理通量 F(Q) = A Q，那叫中心通量。中心通量形式上是能量守恒的，但实际算起来在粗糙网格上极容易产生震荡，波峰后面跟一串毛刺，像锯齿一样。双曲问题里最省事、最稳妥的选择是 Lax-Friedrichs（LF）通量：
[ F^* = \frac{F(Q_L) + F(Q_R)}{2} - \frac{c}{2}(Q_R - Q_L) ]
左边那一项是中心通量，右边减掉的那一项是耗散项，c 是通量雅可比矩阵谱半径的最大特征值。对声波方程这个值是声速 c。LF 通量的核心思想很好理解：在两个状态中间取平均，同时加一点数值黏性把高频振荡压下去。代价是波场会有一点数值耗散，但声波模拟对振幅衰减没那么敏感，这点耗散换来的健壮性非常划算。常见的另一种选择是 Roe 通量和迎风通量，数学上更精致，也更贴近物理特征方向，但代码复杂度明显上去：要按传播方向做特征分解。教学和资源包里最常用的是 LF，先把整体框架跑通，再换迎风通量做精度对比，是更稳妥的路线。
3.2 时间推进与CFL：显式RK3的参数经验
空间离散完了，时间方向用显式 Runge-Kutta 推进。声波方程是标准的双曲型问题，用显式方法就必须受 CFL 条件约束。经验公式是：
% dt 由最小尺寸、声速、多项式阶数共同决定 dx_min = min(mesh.J); % 最小半单元宽度 dt = CFL * dx_min / (c * (N + 1)); % 一维DG的CFL经验公式这里的 CFL 是无量纲库朗数，对三阶 SSP-RK3 方法，经验取值在 0.18 ~ 0.3 之间。初跑时设 0.2 左右最保险，跑通了再往大调。N+1 这个因子经常被忽略：单元内多项式的最高次数越高，对时间步长的限制越严，因为它决定了单元内最短的有效空间尺度。主循环用广播式的 SSP-RK3，写法有讲究：
% SSP-RK3 时间推进主循环 for n = 1:Nt % 第一级 u1 = u0 + dt * dg_rhs(u0); % 第二级 u2 = 0.75 * u0 + 0.25 * (u1 + dt * dg_rhs(u1)); % 第三级 u0 = (1/3) * u0 + (2/3) * (u2 + dt * dg_rhs(u2)); % 周期性输出波场 if mod(n, 20) == 0 plot_result(u0, xe, t(n)); end end这句代码里每级都要调一次 dg_rhs，所以一个完整时间步等于三次右端项求值。SSP-RK3 的好处是它保证强稳定性质，配合 LF 通量这种自带耗散的空间离散，整个格式的稳健性非常高。常见翻车写法是把 RK3 的系数写错，比如把第二级系数写成 0.5 和 0.5，或者把第三级的 1/3、2/3 记反，算出来的波场振幅会在时间方向明显涨落。检查方法很简单：对一个已知解析解的行波跑几个时间步，每隔十步输出一次波峰峰值，正常应该在数值耗散作用下缓慢下降，而不是上下抖动。
3.3 数值通量的Matlab实现：LF函数和中心通量的差别
把 LF 通量写成函数只有短短几行，但方向性是最容易错的地方：
function Fs = lf_flux(Q_L, Q_R, c, rho) % Q_L: 交界面上左侧单元状态 [p; v] % Q_R: 交界面上右侧单元状态 [p; v] % Fs : 输出数值通量 [rho*c^2*v; p/rho] 在交界面的估计值 F_L = [rho * c^2 * Q_L(2); Q_L(1) / rho]; F_R = [rho * c^2 * Q_R(2); Q_R(1) / rho]; Fs = 0.5 * (F_L + F_R) - 0.5 * c * (Q_R - Q_L); end这个函数返回的是交界面上的数值通量，调用它的时候要注意：对每个单元的左边界，F_L 是本单元的值，F_R 是左邻居的值；对单元的右边界则反过来。LF 通量里最后那个耗散项里，Q_R 和 Q_L 的次序不能颠倒，反了以后耗散项就变成增幅项，波场会在交界面处指数增长，程序两三步就爆。这个细节在 5.5 节里还会重点说。中心通量实现时只需把最后那项删掉，但删掉以后波场会在间断处振铃，新手很容易把原因归到 CFL 上，其实源头是通量缺了耗散。
4. 从参考单元到全局组装：质量矩阵、刚度矩阵与可逆性
很多从连续有限元转过来的人，第一个不适应的点就是 DG 的矩阵结构。连续有限元的质量矩阵是全局的一大块带状稀疏矩阵，解一次线性方程组跑一步；DG 因为单元间不共享自由度，质量矩阵天然是块对角的，这一步快得让人惊喜。但快的前提是你要理解它为什么快，以及在哪里可以更快。
4.1 Legendre正交基下的质量矩阵：对角阵白捡
如果参考单元上的基函数选 Legendre 多项式 L_i(r)，利用正交性质：
[ \int_{-1}^{1} L_i(r) L_j(r) dr = \frac{2}{2i+1} \delta_{ij} ]
质量矩阵直接就是对角阵。这个福利是 DG 从理论上就带来的：你们用自由排列的多项式做基，我只需要在参考单元上预先算好一个对角阵，每个单元复制一份就行。
% 参考单元上的质量矩阵，Legendre 正交基时严格对角 M_ref = diag(2 ./ (2*(0:N) + 1)); % 它的逆也是对角阵，预计算一次，循环里直接用 M_inv_ref = diag((2*(0:N) + 1) / 2);注意这里的 N 必须与网格生成时的 N 保持一致。如果你换成节点基（Lagrange 插值基），质量矩阵就不再对角了，但它依然是一个块对角稀疏矩阵，每个块对应一个单元。区别不在于能不能求逆，而在于对角阵的求逆是 O(N) 的，密块求逆是 O(N³) 的。对 N=4 以内的问题两者差别不大，但 N 一高，Legendre 基的优势立刻显现。资源包如果用的是节点基，别慌，它不是错的，只是没吃到全部红利；如果用的是 Legendre 基，那质量矩阵的逆就是上面那行代码。
4.2 刚度矩阵组装：数值积分与参考单元
刚度矩阵在参考单元上的定义为：
[ S_{ij} = \int_{-1}^{1} L_i'(r) L_j(r) dr ]
这个积分没有简单的对角结构，但可以在参考单元上用 Gauss-Legendre 数值积分一次性算好，因为所有单元共享同一套参考单元，算一次就够：
% 参考单元刚度矩阵组装（Gauss-Legendre 数值积分） [rq, wq] = gauss_legendre(2*N+1); % 2N+1 个积分点足够精确 Lq = legendre_eval(N, rq); % 基函数在积分点处的值 dLq = legendre_deriv(N, rq); % 基函数导数在积分点处的值 S_ref = dLq' * diag(wq) * Lq; % 矩阵化写法，避免双重循环代码里的矩阵化写法本质上是把积分定义换成了向量内积：每一列代表一个积分点，计算 S_ij = Σ_k w_k L_i'(r_k) L_j(r_k)。gauss_legendre 函数需要给 2N+1 个积分点，因为被积函数是 N-1 次多项式乘以 N 次多项式，总次数 2N-1，2N+1 个点的 Gauss 积分完全精确。这里有个性能细节：积分点数量宁多勿少，取 2N+3 也不会浪费多少时间，但取少了会引入不可见的离散误差，收敛阶测试时会现原形。
4.3 全局组装还是逐单元推进：按K和N选架构
矩阵预计算完之后，剩下的问题是怎么组织循环。两个做法各有边界，按问题规模选。
% 方案A：全局稀疏矩阵组装（K 在一两万以内适用） M_global = kron(speye(K), M_ref); S_global = kron(speye(K), S_ref); % 方案B：逐单元推进（K 超过五万时首选） for e = 1:K dQe = M_inv_ref * (S_ref * (A * Qe) - ... ); end方案 A 的优点是代码短、容易调试，但你不要在循环里写M_global \ rhs，那等于每一时间步都做一次稀疏分解，白白浪费 DG 块对角的优势。正确做法是一开始就算好M_inv_global = M_global \ speye(size(M_global))，循环里只做矩阵乘法。方案 B 是内存友好型，核心是提前把 M_inv_ref、S_ref 这两个参考矩阵准备好，单元循环里直接乘。哪个更快取决于你机器内存带宽和稀疏矩阵运算效率，但实际体验是 K 超过五万时方案 B 明显更省内存，调试时也更容易定位哪一段逻辑出错。初学阶段用方案 A 跑小算例，生产阶段切方案 B，这是常年跑 DG 的人普遍的做法。
5. 避坑指南：间断有限元声波模拟的五个常见翻车现场
波场模拟类代码的调试完全不同于普通数值计算：中间过程全是动画和数据文件，NaN 爆出来时往往已经被恶性传播了一百个时间步。以下五条是 DG 求解声波方程最常见的翻车记录，每一条都是真实跑过的血泪经验。
5.1 一加网格就 NaN：CFL 没跟着 K 走
现象：程序在 K=100 时跑得好好的，把 K 改成 300 后，跑到几百个时间步突然整场爆成 NaN。原因：网格加密后最小单元尺寸 dx_min 变小，但 dt 还是照着 K=100 时算的，CFL 条件被破坏。显式格式的时间步长极限是网格尺寸的函数，网格加密、时间步长必须同步缩小。解决：强制每次运行都从网格重算 dt。我自己的习惯是在主脚本里用dt = CFL * min(mesh.J) / (c * (N+1))这句，禁止把 dt 写成固定数字存进参数文件。调参时先动 CFL，再动 K，顺序不能反。
5.2 波到了边界像撞了墙：边界条件没做对
现象：高斯脉冲在自由传播中，波前到达边界后反弹回来，产生一列假波，分析人员误以为介质内部存在反射目标。原因：边界默认是反射边界。DG 的边界处理比有限元更隐蔽，因为它不会报错，只是安安静静地反射。实际物理计算里你要的是无反射出射边界。解决：最简单有效的做法是特征边界条件：左边界令 p − Zv = 0，右边界令 p + Zv = 0，其中 Z = ρc 是声阻抗，等价于消除对外特征值方向的反向波。实现时只需在 dg_rhs 里对边界单元调用一次 with 边界状态的通量计算。注意特征边界只能把外向波放出去，对接近垂直入射的波效果好，斜入射波会被部分反射，这时得换 PML 吸收层，代价是代码复杂度上一个台阶。
5.3 质量矩阵求逆把程序拖死了
现象：小规模算例几秒钟跑完，K 加到十万后每跑一步都像卡住，内存也一路飙升。原因：在时间循环内部调用了M_global \ rhs或者inv(M_global)。即使 M_global 是块对角稀疏矩阵，每次反斜杠运算符都会触发一次分解计算，而inv(M_global)更是直接生成满秩稠密矩阵，内存直接爆掉。解决：所有与质量矩阵相关的求逆计算全部移到循环外，预计算M_inv_global = M_global \ speye(size(M_global))，循环里只乘这个预计算矩阵；或者干脆用逐单元方案 B，每个单元手动乘 M_inv_ref。
5.4 多项式阶 N 提高但收敛速率不动：时间误差拖了后腿
现象：N 从 1 加到 4，网格加密一倍，误差几乎不变，甚至上升。原因：空间误差和总误差是两回事。空间收敛阶由 N+1 决定，但总误差还包含时间误差。当你只加 N 不加时间分辨率时，时间误差成了瓶颈，整体收敛曲线走平。另一个常见原因是网格单元数太少，波长根本分辨不了。解决：做收敛阶测试时同时缩小 dt 和网格，观察误差变化。如果误差随 dt 线性下降，那说明时间方向是主导误差源。先把时间步长减半，看空间收敛阶是否回归 N+1；如果回归了，说明程序逻辑本身没问题，是测试参数没配对。
5.5 通量方向性错误：波场左右不对称
现象：一个对称的高斯初始脉冲，跑几百步后左右不对称，或者波峰峰值逐渐大于初始值，完全违背能量守恒。原因：LF 通量中耗散项的符号或左右状态次序写反。耗散项应该是 Q_R − Q_L，一旦写成 Q_L − Q_R，数值上就变成增幅项；边界通量调用时左右邻居取反也会造成方向性错误。解决：用一个标准 Riemann 问题做单元测试。最理想的做法是用初始条件为不连续的阶跃波做试算，然后与解析解对比；或者先只跑几个时间步，检查波场中通量是不是从正方向传播到负方向，排除方向错误后再跑长周期。我见过太多人把这个错误归咎于 CFL 或边界，其实就是一个负号的问题。
6. 验证与快速开始：高斯脉冲、收敛阶测试与资源里优先看哪些文件
代码从下载到出图，中间隔着一条河，河里全是“看起来能跑”和“真的对”之间的差距。我的经验是拿到资源包后先做三件事，按顺序来，每一步都花不了十分钟，但能替你省下后面几十次的返工。
6.1 先跑一个高斯脉冲，确认波场形态
第一个验证不算物理算例，只是形态检查。初始条件取压力场高斯脉冲，速度场置零：
% 初始条件：高斯压力脉冲 x0 = 0.5; sigma = 0.02; p0 = exp(-((mesh.x - x0) / sigma).^2); % 压力场初值 v0 = zeros(size(mesh.x)); % 速度场零初值 u0 = [p0(:); v0(:)]; % 全局解向量跑几百步，先看波场是否分裂成左右两个脉冲，再盯着波峰峰值。正常情况下峰值会略降，这是 LF 通量的数值耗散，属于正常现象；如果峰值升高，立刻去检查通量耗散项符号。对于声波方程，高斯脉冲足够简单，任何一个有经验的数值模拟者都能凭肉眼从波场动画里判断程序对不对。这一步不过关，后面一切算例都不要开始。
6.2 收敛阶测试：用数据说话
形态检查过了，还不能证明程序是高精度的。下一步是收敛阶测试：对同一个光滑行波解，把网格逐步加密，观察误差变化速率。使用精确行波解 g(x − ct)，在波碰到边界之前取误差。
% 收敛阶测试：空间加密，验证格式的真实收敛速率 Ks = [50 100 200 400]; N = 2; % 固定多项式阶 CFL = 0.2; % 固定库朗数 err = zeros(size(Ks)); for j = 1:numel(Ks) [u, x, t] = run_dg(Ks(j), N, CFL); err(j) = max(abs(u(1:2:end-1) - g(x - c*t))); % 压力场误差 end % 相邻两组成倍网格加密时，误差比值的 log2 就是收敛阶 rate = log(err(1:end-1) ./ err(2:end)) / log(2); disp(rate);对 N 次多项式，理论空间收敛阶是 N+1，所以 N=2 时速率应接近 3。如果你的测试结果是 1 点几或 2 点几，先查边界条件是否干净：边界上的反射最容易破坏收敛阶。收敛阶测试是判断 DG 程序是否正确的金标准，它同时检验了网格映射、数值通量、边界条件、时间推进四者的协同。严格说，正式论文里的所有 DG 算例都该附这样的收敛阶表，做不出收敛阶的程序基本都有隐藏 bug。
6.3 资源包里优先看哪几个文件
下载解压后不要急着跑 main，先打开这几个文件过一遍。入口基准脚本管参数设置和网格初始化，重点看 N 和 CFL 是否作为变量传入，如果写在函数体里就只能手动改了；右端项函数是所有算例的核心，看它是否区分内部单元和边界单元；数值通量函数看耗散项符号和左右次序；后处理代码看输出的是物理量还是保守变量，这决定读数据时要不要单位换算。
资源包里如果净是些看不出结构的文件，优先挑文件名含 main、dg、rhs、flux、mesh 的看。要是没有这几个名字，就只能用 6.1 节的高斯脉冲做黑匣子测试，先判断整套代码是否可信。
从那以后，我每次拿到一套新的 DG 代码，都强制自己先走一遍高斯脉冲和收敛阶测试，再去看动画。这套流程救过我太多次，希望帮到你。
本文还有配套的精品资源，点击获取
