间断有限元求解声波方程：MATLAB实现与稳定性避坑指南

简介面向数值分析与偏微分方程求解方向的学习者这份资源聚焦二维声波方程的间断有限元DG数值实现将空间离散与三阶龙格库塔时间推进相结合适合已具备Matlab基础、希望深入掌握间断伽辽金方法的中高级读者。压缩包共6个文件以5个m脚本为主辅以1个asv备份文件整体约4KB体量轻便便于直接阅读与二次修改。内容围绕矩形域上给定正弦初值的声波传播问题涵盖网格划分、线性基函数选取、界面条件处理、边界条件设置与时间步迭代等关键环节读者可据此理解DG方法在复杂几何与高阶精度上的优势并借助Matlab绘图观察解的演化过程。目前已有1567人学习下载可作为数值方法课程设计或科研入门的参考脚本。1. 声波方程为什么值得用间断有限元重写一遍做波动方程数值模拟的人迟早会撞上一个尴尬局面连续有限元在网格加密时收敛漂亮可一旦介质速度反差大、边界形状刁钻或者要算长时间波场演化全局刚度矩阵一拼、一解内存和时间就开始不讲道理。间断有限元Discontinuous GalerkinDG恰好是冲着这个痛点来的——单元之间不强制连续靠数值通量在界面传递信息天然适合并行、适合变阶次、适合处理复杂边界。这份 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/(2N1)$ 量级缩小。这个权衡在调参时比什么都直接。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*N1这个因子来自 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*N1因子重算dt或者换稳定域更大的 RK 格式。还有一种情况是参考单元上的积分精度不够高阶基函数需要更多积分点否则质量矩阵算不准表现为低频漂移。5.4 现象介质速度反差大时界面处出现非物理振荡原因是界面波速取了算术平均而非最大值迎风通量的耗散不足以压制跳跃。把界面波速改成两侧最大值或者对速度反差极大的算例适当增加界面耗散。这是 DG 处理强非均匀介质的经典难点不是调参能完全消除的必要时得换更鲁棒的通量。5.5 现象并行或向量化后结果与串行不一致原因是面循环里对共享边的写入存在顺序依赖或者用了全局变量在单元间传递迹。DG 的面贡献必须按“先算通量、再统一累加”的两步走不能边算边写。向量化时尤其要注意每个面的左右单元索引是否被正确广播。6. 收敛性验证与阶次调参用制造解把代码钉死代码能跑不等于算对。DG 的收敛阶理论上应该是N1L2 范数下拿制造解验证是最省事的办法人为构造一个解析解代进方程算出源项再把源项加进右端项跑不同网格看误差随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_order1);src_term是把解析解代进方程后剩下的非齐次项compute_l2_error用质量矩阵加权算 L2 误差。实测阶次接近N1说明空间离散没问题明显偏低通常是积分精度不够或通量实现有误偏高则要怀疑误差还没进入渐近区网格不够密。我一般会跑三到四级网格看斜率是否稳定只看单级网格的误差没有意义。调参上阶次N和网格尺度h要配合着看目标精度固定时高阶配粗网格往往比低阶配细网格更省因为 DG 的误差常数随阶次下降很快。但高阶对网格质量和时间步更敏感介质复杂时未必划算。我的习惯是先用N2或N3把物理流程跑通确认波形合理再往上调阶次做精度对比而不是一上来就冲高阶。从那以后我每次改完通量或边界代码都强制先跑一遍制造解收敛表斜率不对就不往下走——这个习惯帮我拦下过好几次“看起来能跑、其实算错”的改动。希望帮到你。本文还有配套的精品资源点击获取