简介压缩包内为Matlab实现间断有限元DG求解二维声波方程的完整源码适合具备偏微分方程基础、正在学习高阶数值方法的研究生或仿真工程师。代码将空间离散采用线性基函数的DG格式时间推进使用三阶龙格库塔RK3并配备初始条件设置、通量计算、网格初始化等模块可直接运行观察正弦初始波场在方形区域内的传播演化。包内共6个文件包括5个m脚本与1个asv自动备份文件整体仅4KB文件精炼、结构紧凑便于逐行研读。已有1564人学习下载。对于想快速上手DG方法、理解界面数值通量构造与时间积分配合的读者这份代码提供了清晰的参考骨架可在此基础上进一步扩展至变系数声速、不规则网格或高阶基函数也可与其他有限元实现对比验证。1. 间断有限元求解声波方程一张网格上的波前捕捉实验用 MATLAB 做间断有限元求解声波方程听起来有点“自找麻烦”——毕竟 MATLAB 的循环性能常年被吐槽而 DG 方法又是出了名的自由度多、通量计算繁琐。但当你试过用连续有限元算一个含强间断介质的声波传播看到波前被数值振荡糊成一团时就会理解为什么间断有限元DG值得在 MATLAB 里折腾一遍它允许每个单元内的解不连续用数值通量去连接单元间的信息天然适合处理突变界面和强对流项在粗网格上就能给出锐利的波前。这篇笔记会带你从弱形式推导到一维完整代码再讲清楚边界条件、CFL 限制和常见翻车场景适合想用 MATLAB 快速验证 DG 算法、做声波/弹性波教学实验或小规模正演的读者。2. 间断有限元的数学基础本地重构与数值通量2.1 声波方程的守恒形式与弱形式一维声波方程可以写成一组一阶双曲守恒律∂p/∂t ρc² ∂u/∂x 0 ∂u/∂t (1/ρ) ∂p/∂x 0这里 p 是压力扰动u 是质点速度ρ 是密度c 是声速。写成向量形式就是∂U/∂t ∂F(U)/∂x 0其中 U [p; u]通量 F [ρc²u; p/ρ]。DG 方法的第一步是把计算域剖分为若干单元然后在每个单元上乘以一个测试函数 v做分部积分得到弱形式∫_K (∂U/∂t) v dx - ∫_K F(U) ∂v/∂x dx [F_hat * v]_{边界} 0关键在于最右侧的边界项周期内部单元交界面处的 F_hat 是“数值通量”它不完全等于 p 或 u 在边界上的物理值而是由左右两个单元的解共同决定。这给了我们设计空间——不同的数值通量会导致不同的数值色散和耗散特性。对于声波方程这种线性双曲系统常见的数值通量有中心通量和迎风通量。中心通量就是取左右平均值(F_hat (F(U_L) F(U_R)) / 2)。迎风通量则需要做特征分解根据波传播方向决定取左侧还是右侧的值。在声波方程里特征值为 ±c所以迎风通量的具体形式可以写成矩阵形式对每个特征方向分别选边。2.2 数值通量中心通量与迎风通量如何选我的经验是如果你只想快速跑通 DG用中心通量因为矩阵形式简单不会因为特征分解引入额外计算。但它的问题在于耗散偏小高频震荡在粗网格上更明显。迎风通量虽然多几个矩阵乘法但对声波这种双曲问题能提供更好的稳定性尤其在介质突变处比如声速从 1500 m/s 跳到 5000 m/s不会出现明显的伪振荡。具体选型步骤先确认你要模拟的是均匀介质还是分层介质。均匀介质下两者差别极小。若介质有强不连续优先迎风通量。把通量矩阵写成F_hat (F(U_L)F(U_R))/2 - (1/2) * |A| * (U_R - U_L)其中 A 是通量雅可比矩阵|A| 是对角化后的绝对值矩阵。调试初期用中心通量因为一旦结果发散问题更容易定位在时间积分而非通量上。2.3 基函数与局部质量矩阵间断带来的好处DG 最吸引人的特性之一是每个单元内的质量矩阵只与本单元的基函数相关不需要像连续有限元那样组装全局矩阵。如果你在单元内采用 N 阶多项式基函数那么每个单元的质量矩阵是 (N1)×(N1) 的小矩阵甚至可以预计算逆矩阵。这意味着你可以在不同单元用不同阶次的基函数p 自适应或者用不同大小的单元h 自适应而无需处理跨单元的全局耦合。常见的基函数是勒让德多项式或幂函数。我推荐幂函数基底1, x, x², ...因为实现最简单但要注意当基函数阶数较高时幂函数基底会导致质量矩阵条件数变差。用勒让德多项式基底则更稳但需要额外计算标准正交基。在 MATLAB 里用polyfit或legendreP能快速生成但为了效率我一般直接硬编码低阶基函数。3. 一维声波方程的MATLAB实现完整可跑通的最小代码3.1 网格与自由度一段代码生成均匀网格先定义一个简单的一维问题计算域 [0, 1]声速 c1密度 ρ1初始条件是一个高斯脉冲。我们把它剖成 100 个单元每个单元内用 2 阶多项式P2那么每个单元有 3 个自由度总自由度 300。生成网格的代码如下% 生成均匀网格 - 1D声波方程DG求解L1, N100单元P2阶基函数 L 1; numel 100; P 2; % 多项式阶数 nnode numel 1; edge linspace(0, L, nnode); % 单元边界 nodes zeros(numel, P1); % 存储每个单元内的物理坐标点 h L / numel; for e 1:numel x0 edge(e); x1 edge(e1); legx -1 2 * (0:P) / P; % 局部点 nodes(e,:) x0 (legx 1) * h / 2; % 映射到物理坐标 end % 初始高斯脉冲 xc 0.5; sigma 0.05; u0 zeros(numel, P1, 2); % 两个分量第一页压力第二页速度 u0(:,:,1) exp(-((nodes - xc).^2) / (2*sigma^2)); u0(:,:,2) 0;这段代码的关键在于每个单元内的nodes矩阵存储的是高斯积分点或插值点而基函数是定义在标准区间 [-1,1] 上的局部函数。用legx -1 2*(0:P)/P生成等距插值点虽然这不是高斯积分点但对线性双曲问题足够稳定。如果你想更精确可以用gaussquad函数替换但这里先跑通为主。参数说明numel控制空间分辨率一般声波问题至少需要 20 个单元覆盖一个波长。P是多项式阶数P1 是线性近似P2 是二次近似。P 越大单阶精度越高但自由度也更多时间步长也要相应缩小。3.2 组装刚度矩阵与质量阵接下来计算单元内的局部矩阵。标准 DG 矩阵需要质量矩阵M测试函数与基函数的内积。刚度矩阵K测试函数与基函数导数的内积。边界通量矩阵F_hat左右单元边界处的数值通量相关。在幂函数基底里矩阵元素有解析表达式直接用int或数值积分。这里用高斯积分法做数值积分避免推导公式出错% 计算单元内局部矩阵基于局部坐标[-1,1] % 高斯积分点与权重2点精度足够因为被积函数是多项式 gx [-0.5773502692, 0.5773502692]; gw [1, 1]; Mlocal zeros(P1, P1); % 质量矩阵 Klocal zeros(P1, P1); % 刚度矩阵 % 标准区间 [-1,1] 内基函数 phi_i(xi) xi^(i-1) for i 1:P1 for j 1:P1 for q 1:2 xq gx(q); wq gw(q); Mlocal(i,j) Mlocal(i,j) wq * xq^(i-1) * xq^(j-1); % 导数: d/dxi (xi^(i-1)) * d/dxi (xi^(j-1)) * dx/dxi 的逆 % 对于一维映射dx/dxi h/2刚度项需要乘上 (2/h)^2 吗见下 Klocal(i,j) Klocal(i,j) wq * (i-1)*xq^(i-2) * (j-1)*xq^(j-2); end end end % 注意物理坐标下的刚度矩阵需要乘以 (2/h)^2 不这里Klocal是标准区间导数内积。 % 实际贡献(1/h) * Klocal因为 d/dx (2/h) d/dxi Kphysical Klocal * (2/h); % 此处仅示意具体组合在通量中逻辑说明在标准区间 [-1,1] 上物理坐标 x 与局部坐标 ξ 的关系是 x x0 (ξ1)h/2所以 dx/dξ h/2。弱形式中测试函数导数的贡献需要乘以 (dξ/dx)也就是 2/h。上面的Kphysical实际上应该是(2/h) * Klocal但Klocal已经对 ξ 积分所以物理刚度矩阵是Klocal * (2/h)。质量矩阵与映射无关就是Mlocal * h/2。实际组装时我建议把所有单元矩阵预计算因为均匀网格下每个单元完全相同。如果后续做自适应每个单元的映射不同才需要在循环内重新计算。3.3 半离散格式与时间推进主循环DG 半离散格式可以写为M * dU/dt RHS(U)其中 RHS 包含体积分与边界通量。我们把 RHS 计算封装在一个函数里然后使用 RK2二阶龙格库塔推进时间。下面是核心循环% 时间推进参数 tend 2.0; CFL 0.3; % 针对P2的经验值 dt CFL * h / (c * (2*P1)); % 按DG的特征速度估计 nsteps ceil(tend / dt); dt tend / nsteps; % 预分配中间变量 u u0; % RK2积分器 for n 1:nsteps u1 rk2_stage1(u, nodes, h, c, rho); u2 rk2_stage2(u, u1, nodes, h, c, rho); u u2; end % 定义RK2 stage函数 function rhs compute_rhs(u, nodes, h, c, rho) numel size(u,1); P size(u,2)-1; rhs zeros(numel, P1, 2); % 计算每个单元内体积分在局部坐标上 for e 1:numel % 体积项 -F(U) * dphi/dx 在单元内积分 % 这里用高斯积分遍历省略具体累加 % 边界通量项在左右端点加数值通量 % 左侧边界使用外伸边界条件右侧同理 % 先不考虑边界假设周期边界 % ... end end上面的compute_rhs只是一个骨架完整实现需要填充体积项和边界通量项。因为代码较长建议你把单元操作向量化。在 MATLAB 中用 for 循环遍历单元在小可题几百个单元下可以接受但单元数超过 5000 时就该考虑向量化。我在实际项目中会保留循环但确保内层使用.*和预分配的数组而不是动态扩容。3.4 后处理与误差计算怎么确认算对了算完后需要在单元内重构连续解并画图。一种做法是直接在节点上插值但 DG 解的跳跃使得直接画会显得有“锯齿”。我更推荐在每个单元内单独画曲线% 绘制不同单元内的压力分布 figure; hold on; for e 1:numel x_plot linspace(nodes(e,1), nodes(e,end), 50); p_val polyval(u(e,:,1), linspace(-1,1,50)); % 注意需将x映射到局部坐标 plot(x_plot, p_val, b-); end如果与解析解对比高斯脉冲在周期性边界下会在 2 秒后回到初始位置。计算 L2 误差时需要把 DG 解在每个单元内做积分。我一般先用polyval加密采样再用trapz近似积分。这里有个常见错误不能只比较单元边界的值因为 DG 解在边界不连续必须做单元内积分。4. 边界条件与稳定性控制从CFL数到吸收边界4.1 反射边界与自由边界在DG里的写法声波方程摄测中我们需要区分硬边界法向速度为零压力自由反射和软边界压力为零速度自由。在 DG 方法里边界条件是通过修改边界处的数值通量实现的。以左边界 x0 为例假设硬边界则边界处法向速度为零那么通量中的速度分量为零压力连续性保持。具体实现引入一个“外部单元”的虚值。硬边界条件下外部单元的压力p_ext p_in但速度u_ext -u_in。这样代入中心通量速度项平均后为零自然得到反射。软边界则取p_ext -p_in速度u_ext u_in。而吸收边界无反射边界在声波方程里最常见。一维情况下可以简化为一阶吸收条件在边界上通量只取从边界向外传播的特征分量。代码上就是令边界处的数值通量只由内部解与向外传播的特征组合构成。这相当于特征迎风通量的边界版实现并不复杂是你检查波前消失反射的关键。4.2 RK时间积分为什么二阶RK就够用DG 的空间精度通过提高多项式阶数来提升时间积分则需要匹配相应的阶数。最常见的是 RKDG 方法其中 Runge-Kutta 阶数取空间阶数附近。对于 P1线性基二阶 RK 就足够对于 P2二次基严格来说希望用三阶 RKSSP RK3但实践中我用二阶 RK 配合更小的 CFL 也能稳定只是时间精确度略降。关键原因是声波方程的时间尺度不长而且我们往往更关注波前位置而非时间相位精度二阶 RK 在 CFL 取 0.3 时等效误差已经足够小。所以我的选型法则是P≤2 用 RK2P3 用 RK3。MATLAB 内置的ode45虽然方便但它是变步长的对双曲问题会产生过多的小步长浪费计算。我建议自己写固定步长的 RK 循环代码如下function [U] rk2_step(u, dt, rhs_func, ...) k1 rhs_func(u,...); u_temp u dt * k1; k2 rhs_func(u_temp,...); U u 0.5*dt*(k1 k2); end参数说明rhs_func是计算空间离散 RHS 的函数k1是初始时刻斜率k2是中间时刻斜率。固定步长配合 CFL 条件能确保稳定性。4.3 CFL条件与时间步长选取经验值DG 方法的稳定时间步长不仅依赖于单元尺度和平度还取决于多项式阶数。一个常用的保守估计是dt CFL * h / (c * (2P1))这里 P 是多项式阶数系数 2P1 来自迎风 DG 的谱半径估计。对于中心通量我一般取 CFL ≤ 0.5对于迎风通量可以放宽到 0.8但为了保险还是 0.4 起步。有经验的人可能直接用CFL0.25因为初型阶梯性问题的误差即使略大也不容易爆。实际调试时如果你发现时间推进后数值出现指数发散先检查 CFL 是否大于 1再检查通量矩阵的线性稳定性。不改通量、死减dt是个好办法但若减小 10 倍后依然发散那么问题在空间离散而非时间步长。5. 避坑指南DG求解声波方程最常见的5个翻车点5.1 现象高频震荡越算越乱你在均匀介质中算一个平滑高斯脉冲几个时间步后压力场出现空间高频振荡之后蔓延到全域。原因数值通量耗散不够或者多项式阶数 P 相对于网格分辨率过高。声波方程没有物理耗散中心通量的不连续捕获能力弱容易被高频模激起。解决改用迎风通量或者给边界通量添加一个小的人工粘性项例如在 RHS 中加(1/Re)*∂²p/∂x²但更是修改通量更优雅。我的习惯是先检查初始条件是否在单元边界处存在跳跃。高斯脉冲应跨越至少 3 个单元否则初始投影本身就会引入高频伪波。5.2 现象边界处出现虚假反射用周期边界或反射边界波传到边界后产生附加的数值反射波形畸变。原因边界通量没有按特征方向处理。我用周期边界做测试时最容易忘了把左右边界的通量“配对”导致通量不连续。若用吸收边界但只对压力做了无反射处理忽略速度分量也会反射。解决仔细写出边界通量的矩阵形式。对于左边界需要先计算边界上的特征变量然后明确哪些特征是入境的、哪些是出境的。吸收边界就是让入境特征为零出境特征保留。一段常用代码如下% 左边界吸收边界 c2 c*c; rho2 rho*rho; wl 0.5*(u_left u_ext) - 0.5*rho*c*(u_ext - u_left); % 左行波 % 在此处让左行波出境而不引入右行波 p_flux 0.5*wl; u_flux 0.5*wl/(rho*c);5.3 现象质量矩阵奇异组装局部质量矩阵后求逆失败或线性系统无解。原因使用幂函数基底时如果单元内选择的插值点与基函数阶数不匹配导致矩阵欠秩。比如 P2 时你却放 3 个等距点可能没问题但 P3 时等距点在某些边界条件下会导致矩阵奇异。解决改用勒让德基底或使用高斯-洛巴托积分点。在 MATLAB 中可以用legpoly函数生成正交多项式质量矩阵会变成对角阵彻底避开奇异问题。另一种保险做法是在组装前计算cond(M)若条件数大于 1e12 就换基底。5.4 现象高阶基函数下程序变慢P2 还能勉强跑P3 时 MATLAB 循环耗时暴涨几十倍。原因没有向量化且单元循环内反复调用polyval或polydiff这些函数内部有动态分配。简化数学后RHS 计算中每个自由度需要至少 2P1 次乘加运算P 高了复杂度自然上去但 MATLAB 的循环开销才是主因。解决将单元循环重写为矩阵运算。把所有单元的基函数值预先存储在一个 (P1)×(P1) 的矩阵中然后一次性计算所有单元的通量。对 10000 个单元向量化版本比循环快 30 倍。另一个技巧是用parfor替代for但要注意每次迭代中不能有随机写入冲突所以需要先做稀疏组装。5.5 现象MATLAB运行速度无法忍受当你扩展二维或增加网格密度到 1e6 自由度时等待时间以小时计。原因MATLAB 解释性语言的天生缺陷。DG 方法的组装和 RHS 计算如果全程脚本化会很难受。解决把最热的内循环写成 C MEX 文件或者直接在 MATLAB 中调用 GPU 阵列。对于一维/二维问题我最快的方案是用arrayfun加gpuArray在 4090 上能提速 20 倍。如果不想碰 GPU就尽量把所有数组变成single精度内存减半且计算吞吐量改善。另外把整个时间积分封装进一个函数用 MATLAB 的 JIT 编译老版本 MATLAB 需要开启feature(jit ,1)新版本默认开了。6. 进阶从一维到二维从均匀网格到自适应6.1 二维扩展张量积基函数与网格剖分二维声波方程在 Cartesian 网格上做 DG最自然的做法是采用张量积基函数在 x 方向用 P 阶基y 方向也用 P 阶基每个单元内的自由度就是 (P1)×(P1)。通量计算需要分别处理 x 方向和 y 方向的边界意味着每个单元有 4 个面每个面需要做 2 维积分。MATLAB 实现时我建议把数组组织成(numel, (P1)^2, 2)的形式然后利用reshape实现二维数值通量。这种方法的代价是存储量增加但传写速度反而比非结构网格快因为避免了稀疏矩阵的重排序。如果你需要处理复杂几何建议用四边形网格加等参映射每个单元可以变形成不规则形状但这样基函数在物理坐标下不再是简单的多项式需要额外计算雅可比行列式。6.2 自适应加密在局部突变处自动加密DG 的自适应加密有两种h 加密细分单元和 p 加密提升基函数阶数。对于声波问题我建议先做逐单元残差指示器计算每个单元内解的二阶导数范数或者用相邻单元界的跳跃值。当某个单元的压力梯度超过全局平均值的数倍时把该单元细分为两个子单元。MATLAB 里的实现思路是维护一个单元列表每个单元记录其边界坐标和基函数阶数。加密后需要重新组装邻接关系数值通量只涉及同一物理位置的左右单元所以重头并不大。需要注意加密后时间步长需要按最小单元重新计算否则会局部不稳定。6.3 验证方法用解析解和能量守恒说话要确认代码正确不能只看图形长得像。我的验证流程是先设计一个具有解析解的声波传播问题例如平面波或高斯脉冲在周期边界下回到初始位置然后计算 L2 误差随网格加密的收敛阶。理论上P 阶基函数的 DG 解应该以 P1 阶收敛。如果你测得的收敛阶明显低于理论值说明通量或时间积分有问题。另一个更严格的方法是检查能量守恒。声波方程的无阻尼解具有总能量守恒性质DG 离散后能量耗散取决于通量格式。中心通量几乎不耗散迎风通量有小耗散。在足够细的网格下能量误差应随时间保持稳定不会突然上涨。我用一个简单的energy sum(sum(0.5*(p.^2 c^2*u.^2)))来监测如果能量曲线在某一时刻出现尖刺基本可以定位到边界处理或通量实现有误。最后我想说一句我在 MATLAB 里调通 DG 声波求解器之后最大的教训是“验证比求解更重要”。每一次改动哪怕是换个边界条件我都会先跑一遍初始高斯脉冲的回到原点测试确认波形和幅值不畸变。这个习惯帮我避开了至少三次因通量符号写反而导致的返工。希望这些经验和踩坑记录能帮你在自己的 MATLAB 代码里少走弯路。本文还有配套的精品资源点击获取