间断有限元求解声波方程：从源码解析到二次开发实战

简介：面向数值分析与偏微分方程求解方向的学习者，这份资源聚焦二维声波方程的间断有限元（DG）数值实现，将空间离散与三阶龙格库塔时间推进相结合，适合已具备Matlab基础、希望深入理解间断伽辽金方法的中高级读者。压缩包共6个文件，以m脚本为主，辅以一个asv自动备份文件，整体约4KB，体量轻巧，便于直接阅读与二次修改。资源围绕矩形域上给定正弦初值的声波传播问题，涵盖网格划分、线性基函数选取、界面条件处理、边界条件设置与时间步迭代等核心环节，并可通过Matlab绘图观察解的演化过程。已有1567人学习下载，说明该实现思路在同类数值方法练习中具有一定参考价值。读者可从中获得一套可运行的DG求解框架，理解空间离散与时间积分如何协同工作，并据此扩展到更复杂的几何与边界情形。

1. 从一份“跑不起来”的源码说起：间断有限元求解声波方程到底能干什么

很多人第一次拿到“间断有限元求解声波方程”这类源码时，打开文件夹看到一堆.m文件，心里第一反应是：这玩意儿到底从哪跑？我当初也是这么翻车的——目录里没有main.m，函数名和文件名对不上，mesh变量在三个文件里各定义了一遍。后来才明白，这类项目不是“点一下运行”的 Demo，而是一套把间断伽辽金（DG）方法落到声波方程上的完整数值实验框架。它解决的核心问题是：在一维/二维介质里，用高阶多项式分片逼近声压场，让波在单元界面处通过数值通量耦合，从而模拟声波传播、反射与透射。适合谁？适合正在做波动方程数值模拟、想搞懂 DG 弱形式怎么写成代码、或者要拿一套能改边界条件和介质参数的骨架做二次开发的人。间断有限元的好处是单元局部性好、并行友好、能处理复杂介质，代价是实现细节多，稍不留神就发散或者精度掉阶。

2. 项目结构与技术选型：为什么是 MATLAB 而不是 C++

2.1 目录骨架与文件职责

这类项目通常没有花哨的工程化包装，目录结构朴素但职责清晰。我拆过的版本大致长这样：

soundwave_dg/ ├── main.m % 主驱动：参数、网格、时间推进 ├── mesh/ │ ├── gen_mesh_1d.m % 一维均匀/非均匀网格生成 │ └── gen_mesh_2d.m % 二维结构化四边形网格 ├── basis/ │ ├── legendre_basis.m % Legendre 多项式基与导数 │ └── vandermonde.m % 节点到模态的变换矩阵 ├── flux/ │ ├── lax_friedrichs.m % 局部Lax-Friedrichs数值通量 │ └── upwind.m % 迎风通量（可选） ├── solver/ │ ├── assemble_rhs.m % 组装右端项（体积分+面积分） │ └── time_advance.m % RK4 / 低存储RK 时间推进 └── utils/ ├── plot_field.m % 声压场可视化 └── compute_error.m % 与解析解对比的L2误差

这个骨架的关键在于：basis负责参考单元上的多项式运算，flux负责界面耦合，solver把两者拼成右端项，main只做调度。我第一次改的时候把通量函数直接写进assemble_rhs.m，结果换通量格式要动主循环，血泪经验就是——通量一定要独立成函数，后面调参数、对比格式全靠它。

2.2 为什么选 MATLAB 做 DG 原型

选 MATLAB 不是因为它快，而是因为它把“矩阵拼装”和“可视化”这两件事压到了最低成本。DG 的核心是每个单元上的质量矩阵、刚度矩阵和界面通量项，这些在 MATLAB 里就是几行M = V' * diag(w) * V的事，换成 C++ 要写一堆循环和内存管理。对于声波方程这种需要反复调网格、调多项式阶数、调时间步长的场景，MATLAB 的交互式调试优势明显。代价是：当单元数上万、阶数到 6 以上时，纯 MATLAB 循环会明显拖慢，常见做法是把assemble_rhs里的单元循环向量化，或者对热点用mex编译。我一般会先在 MATLAB 里把格式调对，再考虑性能，别一上来就优化。

2.3 声波方程的 DG 弱形式落点

声波方程写成一阶系统更利于 DG：

∂p/∂t + c² ∇·u = 0 ∂u/∂t + ∇p = 0

其中p是声压，u是质点速度，c是介质声速。在单元K上乘测试函数φ，分部积分后得到体积分项和边界通量项。代码里对应的是assemble_rhs.m中两块：单元内用legendre_basis算∇φ的积分，界面处用lax_friedrichs算数值通量。这里最容易翻车的是通量的符号和法向，一维里左右单元法向相反，二维里每个面的法向要按单元外法向取，写反了波就会“倒着走”。

3. 核心模块拆解：从基函数到时间推进的完整链路

3.1 Legendre 基与 Vandermonde 矩阵

DG 在参考单元[-1,1]上用正交多项式做基。项目里legendre_basis.m一般返回基函数值和导数：

function [phi, dphi] = legendre_basis(x, N) % x: 参考单元上的节点坐标（列向量） % N: 多项式阶数 % phi: 基函数值矩阵，size = [length(x), N+1] % dphi: 基函数导数矩阵 phi = zeros(length(x), N+1); dphi = zeros(length(x), N+1); phi(:,1) = 1; % P0 dphi(:,1) = 0; if N >= 1 phi(:,2) = x; % P1 dphi(:,2) = 1; end for k = 2:N % 三项递推：(k)P_k = (2k-1)xP_{k-1} - (k-1)P_{k-2} phi(:,k+1) = ((2*k-1)*x.*phi(:,k) - (k-1)*phi(:,k-1)) / k; dphi(:,k+1) = ((2*k-1)*(phi(:,k) + x.*dphi(:,k)) - (k-1)*dphi(:,k-1)) / k; end end

这段的关键是递推关系，k从 2 开始，phi(:,k+1)对应P_k。参数N就是常说的“几阶 DG”，一阶就是分片线性，三阶以上能明显看到高精度优势，但时间步长要相应缩小（CFL 条件更严）。vandermonde.m则把节点值转成模态系数：

function V = vandermonde(x, N) % 构造节点到模态的变换矩阵 [phi, ~] = legendre_basis(x, N); V = phi; % 若节点为Gauss-Lobatto点，V可逆 end

实际项目里V用来做u_nodes = V * u_modes的转换，质量矩阵M = inv(V)' * diag(w) * inv(V)。我一般会预计算inv(V)存起来，别在时间循环里反复求逆，那是纯浪费。

3.2 数值通量：Lax-Friedrichs 的写法与参数

界面通量是 DG 的灵魂。lax_friedrichs.m典型实现：

function F = lax_friedrichs(pL, pR, uL, uR, c, nx) % pL,pR: 左右单元界面处的声压 % uL,uR: 左右单元界面处的速度（法向分量） % c: 局部声速 % nx: 外法向（一维取 +1/-1） alpha = max(c); % 最大特征速度 p_avg = 0.5*(pL + pR); u_avg = 0.5*(uL + uR); % 声波系统通量：F_p = c^2 * u, F_u = p F_p = c^2 * u_avg + 0.5*alpha*(pL - pR)*nx; F_u = p_avg + 0.5*alpha*(uL - uR)*nx; F = [F_p; F_u]; end

alpha取局部最大声速，保证耗散足够压制界面振荡；nx是法向，一维里左单元取+1、右单元取-1，二维里按面法向。这里有个坑：alpha如果取全局最大c，在非均匀介质里会过度耗散，高阶精度会掉；常见做法是取界面两侧c的较大值。另外nx的符号必须和pL/pR的左右定义一致，否则通量方向反了，波会“撞墙反弹”而不是透射。

3.3 右端项组装：体积分与面积分怎么拼

assemble_rhs.m是计算量最大的地方，结构一般是“先算单元内，再算界面”：

function rhs = assemble_rhs(p, u, mesh, basis, flux_func) % p,u: 模态系数向量（按单元分块） % mesh: 含节点坐标、单元连接、面信息 % basis: 预计算的基函数、质量矩阵逆 rhs_p = zeros(size(p)); rhs_u = zeros(size(u)); for e = 1:mesh.nelem % 体积分：c^2 * ∇·u 和 ∇p idx = mesh.elem_dof(e,:); rhs_p(idx) = -mesh.c(e)^2 * (basis.D * u(idx)); rhs_u(idx) = -(basis.D * p(idx)); end % 面积分：遍历界面，加通量修正 for f = 1:mesh.nface eL = mesh.face_elem(f,1); eR = mesh.face_elem(f,2); % 取界面两侧的迹值，调用通量函数 [pL,pR,uL,uR] = trace_values(p, u, eL, eR, mesh, basis); F = flux_func(pL, pR, uL, uR, mesh.c([eL eR]), mesh.face_normal(f)); % 把通量修正加回左右单元的右端项 rhs_p(mesh.elem_dof(eL,:)) += basis.face_mass * F(1); rhs_u(mesh.elem_dof(eL,:)) += basis.face_mass * F(2); % 右单元符号相反 rhs_p(mesh.elem_dof(eR,:)) -= basis.face_mass * F(1); rhs_u(mesh.elem_dof(eR,:)) -= basis.face_mass * F(2); end end

basis.D是参考单元上的微分矩阵，basis.face_mass是界面上的质量项。注意右单元那两行的减号，这是法向相反的直接体现，漏了就会能量不守恒。trace_values负责从模态系数取界面节点值，一维里就是取两端点，二维里按面取边上的点。这个循环在单元多的时候是瓶颈，我一般会把elem_dof预存成稀疏索引，避免每次mesh结构体访问。

3.4 时间推进：RK4 与 CFL 约束

time_advance.m用经典四阶 Runge-Kutta：

function [p, u, t] = time_advance(p, u, dt, T, mesh, basis, flux_func) % dt: 时间步长，需满足 CFL % T: 总模拟时间 t = 0; while t < T [k1p,k1u] = assemble_rhs(p, u, mesh, basis, flux_func); [k2p,k2u] = assemble_rhs(p+0.5*dt*k1p, u+0.5*dt*k1u, mesh, basis, flux_func); [k3p,k3u] = assemble_rhs(p+0.5*dt*k2p, u+0.5*dt*k2u, mesh, basis, flux_func); [k4p,k4u] = assemble_rhs(p+dt*k3p, u+dt*k3u, mesh, basis, flux_func); p = p + dt/6*(k1p + 2*k2p + 2*k3p + k4p); u = u + dt/6*(k1u + 2*k2u + 2*k3u + k4u); t = t + dt; end end

dt不能随便给，DG 的 CFL 条件大致是dt <= C * h / (c * (2N+1))，h是单元尺度，N是阶数。阶数越高、网格越细，dt越小。我见过有人把dt设成0.1跑三阶 DG，结果第一步就NaN，还以为是通量写错了，其实是 CFL 炸了。稳妥做法是先算h和c，取C=0.1试跑，看能量是否稳定再放大。

4. 跑起来：参数配置、边界条件与结果验证

4.1 主驱动脚本的参数怎么改

main.m一般集中了所有可调参数，典型结构：

% 物理参数 c0 = 1.0; % 背景声速 T = 2.0; % 总时间 % 网格参数 nelem = 100; % 单元数 N = 3; % 多项式阶数 % 时间步长（按CFL估计） h = 2/nelem; dt = 0.1 * h / (c0*(2*N+1)); % 初始条件：高斯脉冲 p0 = @(x) exp(-100*(x-0.5).^2);

改nelem看收敛阶，改N看高阶精度，改c0看介质变化。p0是初始声压，常见做法是高斯脉冲或者正弦波包。注意dt那行是估计值，实际跑之前建议先跑 10 步看max(abs(p))有没有爆掉。

4.2 边界条件的处理方式

声波方程常见边界有硬壁（反射）、吸收（无反射）、周期。项目里一般用“虚拟单元”或者“边界通量”实现：

% 硬壁边界：法向速度置零，声压对称 if is_boundary(e) u_bnd = 0; p_bnd = p(e); % 对称延拓 end % 吸收边界：用一阶Engquist-Majda条件 p_bnd = p(e) - c * u(e) * dt / h;

硬壁最简单，反射完全；吸收边界要小心，一阶条件在大角度入射时仍有反射，常见做法是加一层 PML（完美匹配层），但那要改网格和方程，属于二次开发范畴。我一般先用硬壁验证格式，再加吸收。

4.3 用解析解验证收敛阶

验证是这类项目最容易被跳过的一步。做法是：给一个已知解析解（比如一维行波p = sin(k(x-ct))），跑不同nelem，算 L2 误差：

err = zeros(1, length(nelem_list)); for i = 1:length(nelem_list) [p, ~] = run_solver(nelem_list(i), N, T); err(i) = compute_error(p, exact_solution, mesh); end % 拟合收敛阶 order = polyfit(log(nelem_list), log(err), 1);

N阶 DG 理论上 L2 收敛阶是N+1，如果拟合出来只有 1 阶，八成是通量或者边界处理掉了精度。这个对比表我一般会存下来，改代码后跑一遍，防止“改一处崩一处”。

（表中数值为示意，实际以自己跑出来的为准。）看N=3那列，单元数翻倍误差降约 16 倍，对应 4 阶，符合预期。

5. 避坑与排查：那些让波“消失”或“爆炸”的细节

5.1 现象：跑几步就 NaN，波场全白

原因：CFL 超限，dt太大。DG 的稳定域比有限差分窄，高阶尤其明显。解决：按dt = 0.1*h/(c*(2N+1))重算，先跑 10 步看max(abs(p))是否单调增长。如果还炸，检查通量里的alpha是不是取了负值或者零。

5.2 现象：波传播方向反了，脉冲往左跑

原因：界面法向nx符号写反，或者右单元通量修正的加减号搞错。解决：一维里固定“左单元法向 +1、右单元 -1”，在assemble_rhs里右单元那两行必须是减号。二维里打印mesh.face_normal抽查几个面，确认指向单元外部。

5.3 现象：收敛阶只有 1 阶，高阶基没起作用

原因：vandermonde矩阵和节点不匹配，或者质量矩阵逆算错。常见是节点用了等距点但基是 Legendre，导致V病态。解决：改用 Gauss-Lobatto 节点，或者直接用模态形式避免求逆。检查basis.D对x^N的导数是否精确。

5.4 现象：边界处出现非物理振荡

原因：吸收边界条件太粗糙，或者硬壁边界用了单侧通量但没做对称延拓。解决：先用硬壁验证内部格式，再加吸收；吸收边界附近加密网格，或者上 PML。别在边界没调好的时候怀疑内部通量。

5.5 现象：改介质声速后能量不守恒

原因：c在界面处跳变，但通量里alpha取了全局最大值，过度耗散。解决：alpha取界面两侧max(cL,cR)，并在assemble_rhs里按面取局部c，别用单元平均。

6. 二次开发：换通量格式、加 PML 与向量化提速

拿到这套骨架后，最值当的改造方向有三个。第一是换数值通量：把lax_friedrichs换成upwind或者HLL，对比耗散和精度。改的时候只动flux/下的函数，assemble_rhs的接口不变，这是当初把通量独立出来的回报。第二是加 PML 吸收层：在计算域外围加一层衰减介质，把∂p/∂t改成∂p/∂t + σ(x)p，σ在层内从 0 平滑增长到最大值。实现上是在assemble_rhs的体积分里加一个-sigma.*p项，网格生成时标记 PML 单元。第三是向量化提速：把assemble_rhs里的单元循环改成按elem_dof的稀疏矩阵乘法，界面循环用accumarray聚合。我实测过，nelem=2000、N=4时向量化能快 3 到 5 倍，代价是内存涨一些。

验证改造是否成功，我习惯用“三件套”：解析解收敛阶、能量随时间的变化曲线、以及和旧版本的波场快照对比。收敛阶不掉、能量不增、快照对得上，才算改对了。从那以后我每次动通量或者边界，都强制先跑一遍收敛阶测试再往下写，不然调到最后根本分不清是格式错了还是参数错了。希望这套拆解能帮到你，把这份源码真正跑起来、改下去。

本文还有配套的精品资源，点击获取