简介：面向数值分析与偏微分方程求解方向的学习者，这份资源聚焦二维声波方程的间断有限元（DG）数值实现，将空间离散与三阶龙格库塔时间推进相结合，适合已具备Matlab基础、希望深入掌握间断伽辽金方法的中高级读者。压缩包共6个文件，以5个m脚本为主，辅以1个asv备份文件，整体约4KB，体量轻便，便于直接阅读与二次修改。内容围绕矩形域上给定正弦初值的声波传播问题，涵盖网格划分、线性基函数选取、界面条件处理、边界条件设置与时间步迭代等关键环节，读者可据此理解DG方法在复杂几何与高阶精度上的优势，并借助Matlab绘图观察解的演化过程。目前已有1567人学习下载，可作为数值方法课程设计或科研入门的参考脚本。
1. 声波方程为什么值得用间断有限元重写一遍
做波动方程数值模拟的人，迟早会撞上一个尴尬局面：连续有限元在网格加密时收敛漂亮，可一旦介质速度反差大、边界形状刁钻，或者要算长时间波场演化，全局刚度矩阵一拼、一解，内存和时间就开始不讲道理。间断有限元（Discontinuous Galerkin，DG）恰好是冲着这个痛点来的——单元之间不强制连续，靠数值通量在界面传递信息，天然适合并行、适合变阶次、适合处理复杂边界。这份 MATLAB 资源就是围绕声波方程搭的一套 DG 求解流程，把弱形式、数值通量、时间推进和边界处理串成可跑的脚本。它适合已经懂一点有限元、想亲手把 DG 从公式落到代码上的从业者，也适合拿它当教学 Demo 反复改参数的人。下面按“先立住原理、再动手复现、最后盯住坑”的顺序拆。
2. 声波方程的 DG 弱形式：从标量波动方程到半离散系统
2.1 一阶化与通量选择：为什么先降阶再离散
标量声波方程写成二阶形式是 $p_{tt} = c^2 \nabla^2 p$，直接对二阶时间项做 DG 离散，界面上的通量定义会变得别扭。常见做法是先引入速度变量 $\mathbf{v}$，把方程降成一阶双曲系统：
$$ \frac{\partial p}{\partial t} = -c^2 \nabla \cdot \mathbf{v}, \qquad \frac{\partial \mathbf{v}}{\partial t} = -\nabla p $$
这样每个单元内部都是标准的一阶守恒律形式，界面通量可以套用数值通量框架。选型上，声波这类线性双曲系统最常用的是迎风型通量（upwind flux），它根据特征线方向决定信息从哪一侧传入，稳定性有理论保证。资源里如果用的是中心通量，数值耗散会偏小，长时间算容易出现高频振荡，这一点在改代码前要先确认。
弱形式推导的关键一步是对通量项做分部积分，把散度从试探函数身上挪到测试函数身上，代价是每个单元边界多出一个面积分项。这个面积分项正是 DG 与连续有限元的本质区别：连续有限元在拼装时把界面贡献抵消掉，DG 则显式保留，靠数值通量把左右单元耦合起来。理解这一点，后面看代码里的face循环就不会觉得是多余开销。
2.2 基函数与质量矩阵：为什么 DG 的质量矩阵可以求逆
DG 的另一个好处是基函数只在单元内部有定义，跨单元为零。于是全局质量矩阵是块对角的，每个块的规模等于单元内基函数个数乘以场变量个数。对一阶系统来说，这个块通常很小，直接求逆完全可行，时间推进可以用显式格式而不必解大型线性系统。
资源里一般会用模态基（比如 Legendre 多项式）或节点基。模态基的好处是质量矩阵在参考单元上接近对角，求逆几乎零成本；节点基的好处是物理量直接落在节点上，后处理方便。我一般会先看基函数阶次N是怎么设的：N越大单单元精度越高，但显式时间步受 CFL 条件限制会收紧，步长大致按 $1/(2N+1)$ 量级缩小。这个权衡在调参时比什么都直接。
2.3 半离散系统与时间推进：显式 RK 的步长约束
空间离散完得到的是 $\frac{d\mathbf{U}}{dt} = \mathbf{L}(\mathbf{U})$，其中 $\mathbf{L}$ 包含体项积分和界面通量。时间推进常用低存储显式 Runge-Kutta，三阶或四阶都行。步长由 CFL 条件控制：
% 估算显式时间步：c 为最大波速，h 为最小单元尺度，N 为基函数阶次 c = max(cellfun(@(x) max(x(:)), c_elem)); % 逐单元取最大速度 h = min(h_elem); % 全局最小单元尺寸 N = poly_order; % 多项式阶次 CFL = 0.2; % 保守取值，先跑稳再放大 dt = CFL * h / (c * (2*N + 1));这段估算里CFL取 0.2 是保守起点，实际能放到多少取决于通量类型和网格质量，迎风通量通常比中心通量允许更大步长。2*N+1这个因子来自 DG 谱半径的经典估计，阶次越高步长越紧，这是显式 DG 的固有代价。跑之前先用小网格试几步，看能量是否单调不增，再放大 CFL，比一上来就调大参数稳妥得多。
3. 在 MATLAB 里把 DG 循环搭起来：网格、面循环与通量装配
3.1 网格与单元数据结构：先定好连接关系
MATLAB 里做 DG，网格结构决定了后面所有循环怎么写。一维情形最简单，单元就是区间，左右邻居靠索引直接找；二维三角网格则需要显式的边表和邻接表。资源里通常会给一个mesh结构体，至少包含节点坐标、单元连接、边连接和边界标记。我建议先把这几样打印出来核对一遍：
% 检查网格结构完整性：节点数、单元数、边数、边界边数 fprintf('nodes=%d, elems=%d, edges=%d, bnd_edges=%d\n', ... size(mesh.nodes,1), size(mesh.elems,1), ... size(mesh.edges,1), sum(mesh.edge_tag ~= 0)); % 检查每个单元的面积/长度是否为正，负值说明节点顺序反了 for e = 1:size(mesh.elems,1) idx = mesh.elems(e,:); area = polyarea(mesh.nodes(idx,1), mesh.nodes(idx,2)); if area <= 0 warning('单元 %d 面积为非正：%g', e, area); end end节点顺序反了是二维 DG 最常见的翻车点之一，面积算出来是负的，后面法向量方向全错，通量符号跟着错，波场会莫名其妙发散。先跑这段检查，比事后对着波形猜原因省时间。边表里edge_tag非零表示边界边，边界条件就挂在这些边上处理。
3.2 面循环与数值通量：迎风通量的实现细节
面循环是 DG 的核心。对每条内部边，要取左右两侧单元的解，按特征速度判断信息流向，再组装通量。一维声波系统的迎风通量可以写成矩阵形式，二维则按法向投影。下面是一维标量情形下的通量骨架，二维把标量换成向量、乘上法向分量即可：
% 一维迎风通量：pL/vL 为左侧迹，pR/vR 为右侧迹，c 为局部波速 function [Fp, Fv] = upwind_flux_1d(pL, vL, pR, vR, c) % 特征变量分解：右行波与左行波 w1L = pL - c*vL; w1R = pR - c*vR; % 右行特征 w2L = pL + c*vL; w2R = pR + c*vR; % 左行特征 % 迎风选取：右行取左侧，左行取右侧 w1 = w1L; w2 = w2R; % 还原通量 Fp = 0.5*(w1 + w2); Fv = 0.5*(w2 - w1)/c; end这里w1、w2是特征变量，迎风的意思就是“信息从哪边来就用哪边的值”。Fp、Fv分别是压力和速度方程在界面上的数值通量，回填到弱形式的面积分项里。注意c要取界面两侧的某种平均或最大值，介质不连续时取最大波速偏保守但稳定。二维情形下把v换成v·n，法向量n的方向必须与左右单元约定一致，否则通量会反号，这是第二个高频翻车点。
3.3 体项积分与矩阵预组装：把能提前算的都算掉
体项积分在每个时间步里都要用，但参考单元上的质量矩阵、刚度矩阵与时间无关，完全可以预组装。资源里如果每个 RK 子步都重新算一遍基函数积分，速度会慢一个量级。正确做法是把参考单元上的算子算一次存起来，时间循环里只做“单元局部解乘以预存矩阵再映射回全局”的乘法：
% 预组装：参考单元上的质量矩阵 M_ref 与刚度矩阵 K_ref [M_ref, K_ref] = assemble_ref_operators(poly_order); % 时间循环内只做局部乘法，不再重复积分 for step = 1:n_steps for e = 1:n_elems ue = U(e,:); % 当前单元局部自由度 rhs(e,:) = -c_elem(e)^2 * (K_ref*ue) / M_ref; % 体项贡献 end % 面循环叠加通量贡献（略） U = rk_step(U, rhs, dt); endM_ref和K_ref只依赖阶次，与网格无关，这是 DG 相对连续有限元在实现上的一个便利。rhs先装体项，面循环再往里加通量项，顺序不能反，否则边界贡献会被覆盖。rk_step用低存储 RK 时注意每个子步都要重新算通量，因为界面迹随子步变化。
4. 边界条件与稳定性：吸收边界、CFL 与能量监控
4.1 吸收边界怎么加：一阶 Engquist-Majda 的落地
声波模拟最怕边界反射把内部波场污染。资源里一般会给一阶吸收边界条件，本质是在边界上只允许出射特征波，入射特征置零。一维情形下就是在边界面上令入射特征变量为零，二维则按法向分解。实现上它和内部通量共用同一套特征分解代码，只是边界侧的特征值取零：
% 左边界吸收：只保留出射特征，入射特征置零 function [Fp, Fv] = abc_left(p_in, v_in, c) w_out = p_in - c*v_in; % 出射（右行）特征 w_in = 0; % 入射特征强制为零 Fp = 0.5*(w_out + w_in); Fv = 0.5*(w_in - w_out)/c; end一阶吸收边界在正入射时效果不错，斜入射会有残余反射，这是方法本身的局限，不是代码 bug。想要更好效果得上 PML，但 PML 在 DG 里要额外引入辅助变量，复杂度上一个台阶。先用一阶 ABC 把流程跑通，再决定要不要升级，是我一贯的顺序。
4.2 能量监控：判断发散的第一手证据
显式 DG 跑飞之前往往有征兆，最直接的就是离散能量。每个时间步算一下 $\sum_e \int (p^2/c^2 + |\mathbf{v}|^2)$，正常情况下应该单调不增（有耗散通量时）或近似守恒（中心通量时）。一旦看到能量指数上升，基本可以断定是 CFL 超限、通量符号错或边界处理漏了项：
% 每若干步记录一次离散能量，用于判断稳定性 energy = zeros(n_steps,1); for step = 1:n_steps e_now = 0; for e = 1:n_elems ue = U(e,:); e_now = e_now + ue * (M_ref * ue'); % 局部能量 end energy(step) = e_now; if mod(step, 50) == 0 && e_now > energy(max(1,step-50)) * 1.05 warning('第 %d 步能量上升超过 5%%，检查 CFL 与通量符号', step); end end这段监控不占多少计算量，却能省下大量事后排查时间。M_ref在这里当权重用，保证能量是离散内积意义下的量。阈值 5% 是经验值，网格粗、阶次低时可以放宽，但持续上升就一定要停下来查。
5. 避坑与排查：五个让 DG 跑飞的常见原因
5.1 现象：波场几步内指数发散，能量监控直线上升
原因多半是 CFL 超限或通量符号反了。先按 2.3 的公式把dt砍半再跑，如果还发散，检查面循环里左右单元的迹是否取反、法向量方向是否与单元编号约定一致。迎风通量对符号极其敏感，一个负号就能让整个系统变成“反耗散”。
5.2 现象：边界附近出现强反射，内部波场被污染
原因是一阶吸收边界的入射特征没有真正置零，或者边界边的edge_tag没被识别到，边界被当成内部面处理。检查边表里边界标记是否非零、边界分支是否真的进了abc函数。另一个可能是边界法向朝内，出射特征判断反了。
5.3 现象：提高多项式阶次后结果反而变差
原因是高阶基函数下时间步没跟着收紧，显式 RK 的稳定域被突破。按2*N+1因子重算dt，或者换稳定域更大的 RK 格式。还有一种情况是参考单元上的积分精度不够，高阶基函数需要更多积分点，否则质量矩阵算不准，表现为低频漂移。
5.4 现象：介质速度反差大时界面处出现非物理振荡
原因是界面波速取了算术平均而非最大值，迎风通量的耗散不足以压制跳跃。把界面波速改成两侧最大值，或者对速度反差极大的算例适当增加界面耗散。这是 DG 处理强非均匀介质的经典难点，不是调参能完全消除的，必要时得换更鲁棒的通量。
5.5 现象：并行或向量化后结果与串行不一致
原因是面循环里对共享边的写入存在顺序依赖，或者用了全局变量在单元间传递迹。DG 的面贡献必须按“先算通量、再统一累加”的两步走，不能边算边写。向量化时尤其要注意每个面的左右单元索引是否被正确广播。
6. 收敛性验证与阶次调参：用制造解把代码钉死
代码能跑不等于算对。DG 的收敛阶理论上应该是N+1（L2 范数下），拿制造解验证是最省事的办法：人为构造一个解析解，代进方程算出源项，再把源项加进右端项，跑不同网格看误差随h的下降斜率。下面是一维情形的验证骨架：
% 制造解验证：p_exact = sin(pi*x - c*pi*t)，对应源项解析可得 err = zeros(n_levels,1); hs = zeros(n_levels,1); for lv = 1:n_levels mesh = build_uniform_mesh(2^lv); % 逐级加密 U = run_dg_with_source(mesh, @src_term, T_end); err(lv) = compute_l2_error(U, mesh, T_end, @p_exact); hs(lv) = 1/2^lv; end % 拟合收敛阶：log(err) 对 log(h) 的斜率 p_fit = polyfit(log(hs), log(err), 1); fprintf('实测收敛阶约 %.2f（理论 %d）\n', -p_fit(1), poly_order+1);src_term是把解析解代进方程后剩下的非齐次项，compute_l2_error用质量矩阵加权算 L2 误差。实测阶次接近N+1说明空间离散没问题；明显偏低通常是积分精度不够或通量实现有误；偏高则要怀疑误差还没进入渐近区，网格不够密。我一般会跑三到四级网格，看斜率是否稳定，只看单级网格的误差没有意义。
调参上，阶次N和网格尺度h要配合着看：目标精度固定时，高阶配粗网格往往比低阶配细网格更省，因为 DG 的误差常数随阶次下降很快。但高阶对网格质量和时间步更敏感，介质复杂时未必划算。我的习惯是先用N=2或N=3把物理流程跑通，确认波形合理，再往上调阶次做精度对比，而不是一上来就冲高阶。从那以后我每次改完通量或边界代码，都强制先跑一遍制造解收敛表，斜率不对就不往下走——这个习惯帮我拦下过好几次“看起来能跑、其实算错”的改动。希望帮到你。
本文还有配套的精品资源，点击获取
