简介这是一份基于Matlab的间断有限元DG求解二维声波方程的资源包面向需要数值求解偏微分方程、关注高阶精度方法的研究生和工程计算人员。代码采用线性基函数进行空间离散并用三阶龙格库塔法推进时间完整实现了从网格划分、基函数定义、矩阵组装、界面通量处理到时间迭代的DG求解流程。压缩包共6个文件包括5个.m脚本和1个.asv备份文件分别承担初值设置、通量计算、主循环推进等分工代码体量精简便于对照论文逐步理解DG方法的核心步骤。目前已有1564人学习下载适合作为入门间断伽辽金方法的参考实现。运行后可直观看到二维声波传播的数值模拟效果尤其能帮助体会DG方法在界面不连续处理与高精度时间积分方面的特点。1. 间断有限元求解声波方程为什么用MATLAB从一维开始而不直接上二维间断有限元间断伽辽金求解声波方程第一反应是“这不就是普通有限元吗”真跑起来才发现单元之间自由度完全不连续每个单元独立做多项式展开跨界面全靠数值通量连接。我最早用有限差分算声波两种介质界面处总出现人为反射换成间断有限元后界面通量反而干净。这份MATLAB资源从一维P1DG起步一直到二维多介质算例覆盖网格生成、LF通量、RK3时间推进和可视化适合正在做声波模拟的研究生和想试新格式的工程师。不需要把DG理论背熟拿着代码改参数就能跑出波场快照。2. 数学模型与DG离散给声波方程写出一套分片多项式的半离散格式2.1 一阶双曲形式的声波方程声波方程在压力形式里是二阶偏微分方程但DG框架里更常用一阶双曲守恒律。把压力 p 和质点速度 v 作为两个未知量写成一阶系统∂p/∂t ρc² ∂v/∂x 0∂v/∂t (1/ρ) ∂p/∂x 0其中 ρ 是介质密度c 是声速ρc² 正好是体积模量。这样写的好处是守恒律理论可以直接用线性双曲性质清晰特征波速 ±c 也容易求。后续的LF数值通量、CFL条件全部建立在这一阶形式上。我一般先把介质参数独立放成一个结构体这样换介质时不需要改主程序% 一维模型参数rho 密度c 声速K 为体积模量 rho 1.0; % 单位 kg/m^3 c 1500.0; % 单位 m/s K rho * c^2; % 单位 Pa % 高斯脉冲中心位置与半宽 x0 0.3; sigma 0.02; % 初始压力场初始速度为0 p_init (x) exp(-((x - x0)/sigma).^2); % 状态变量统一用 [p; v] 表示 q_init (x) [p_init(x); 0.0];参数说明密度和声速必须严格为正K 在厘米克秒制下会很大导致通量量级差异明显。如果模拟海水介质ρ 取 1000、c 取 1500模拟气体则要换小几个数量级否则LF通量中的 C 会偏大时间步长被压得很小。这个初始条件只是测试用的高斯脉冲后面做二维时建议改用Ricker子波有更好的频谱行为。2.2 间断伽辽金弱形式与LF数值通量在每个单元上把 p 和 v 用分片多项式展开比如单元 i 上的值写成q_h(x,t) Σ_{j0}^{deg} q_{i,j}(t) φ_j(x)φ_j 是参考单元 [-1,1] 上的基函数。用同一个基函数做检验函数在单元内积分分部积分后会出现单元边界上的通量项。由于相邻单元的 q_h 在界面处不一样这个通量不能两边同时取必须用一个数值通量代替。这就是“间断”二字的关键。最简单的稳定选择是Lax-Friedrichs通量也叫Rusanov通量。对一维线性声波方程通量函数 f(q)[ρc² v; p/ρ]两侧值分别代入再减去界面最大特征速度乘以状态跳跃的一半function Fstar LF_flux(qL, qR, matL, matR) % 输入 qL,qR: 左右单元在界面处的状态 [p; v] % matL,matR: 左右介质参数 [rho; c] rhoL matL(1); cL matL(2); rhoR matR(1); cR matR(2); FL [rhoL * cL^2 * qL(2); qL(1) / rhoL]; FR [rhoR * cR^2 * qR(2); qR(1) / rhoR]; % 界面最大波速估计 C max(cL, cR); Fstar 0.5 * (FL FR) - 0.5 * C * (qR - qL); end参数说明这里 C 取 max(cL,cR) 而不是平均是为了加密熵稳定。如果你图省事改用中心通量即去掉后半部分非均匀网格或强弱间断处容易产生振荡。我最初就这么干过结果在界面处出现明显的锯齿波。LF通量虽然耗散大一些但胜在稳新手阶段先用它不会错。2.3 正交基函数与对角质量矩阵一维上最省事的基函数是勒让德多项式。它在 [-1,1] 上正交质量矩阵自动退化为对角阵。这意味着每个单元内的自由度更新互相独立不需要组装全局矩阵、不需要求逆。对MATLAB来说数组维度就是 单元数×(阶数1)×2内存也直观。我习惯把基函数求值和对应的高斯积分封装在一起。下面这段生成参考点上的基函数值和导数function [phi, dphi_dxi] legendre_basis(deg, xi) % 计算勒让德基函数在参考点 xi 处的函数值和导数值 % xi 可以是列向量 num length(xi); phi zeros(num, deg1); dphi_dxi zeros(num, deg1); % 0 阶 phi(:,1) 1.0; dphi_dxi(:,1) 0.0; if deg 0 return; end % 1 阶 phi(:,2) xi; dphi_dxi(:,2) 1.0; % 递推关系 for k 2 : deg phi(:,k1) ((2*k-1) * xi .* phi(:,k) ... - (k-1) * phi(:,k-1)) / k; dphi_dxi(:,k1) phi(:,k-1) ... xi .* dphi_dxi(:,k) ... - (k-1)/k * dphi_dxi(:,k-1); end end参数说明deg0 是分片常数格式太耗散一般不用deg1 是线性DGdeg3 以上对光滑声波场能明显压制频散但CFL步长会按 (2deg1) 的比例收紧。注意递推公式里 dphi_dxi 的表达式容易写错我就在这里踩过坑如果直接对递推公式求导少算一项导数高频分量会莫名放大。有了对角质量矩阵半离散方程就是每个自由度的常微分方程。右端项包含两部分单元内部基函数与通量梯度的积分以及界面数值通量的修正。具体代码放到第3章因为那里要配合完整的数据结构单独讲右端项很难说清楚。还需要记住一点单元映射到参考区间后物理坐标与参考坐标的雅可比是 h/2这个系数会出现在所有体积分和面积分里。3. MATLAB实现一维DG求解的完整流程3.1 网格生成与高斯积分节点一维DG最忌讳用等距节点做高次多项式积分一定要用高斯求积。每个单元的物理坐标 x 映射到参考坐标 ξx xL (ξ1)/2 * h。雅可比是 h/2。网格生成时只需要单元左右端点剩下交给基函数。% 求解域 [0,1]100个单元线性DG N 100; deg 1; xL 0; xR 1; h (xR - xL) / N; xNodes linspace(xL, xR, N1); % N1个网格顶点 % 高斯积分点与权重nq deg2 保证 2*deg1 阶精确 nq deg 2; [xiq, wq] gauss_legendre(nq); % 每个单元的物理坐标映射 xq zeros(nq, N); for e 1 : N xL_e xNodes(e); xR_e xNodes(e1); xq(:,e) xL_e 0.5 * (xiq 1) * h; end % 体积转化系数 detJ h / 2;参数说明高斯积分节点阶数至少要比被积多项式阶数高这里取 deg2 点点数过少会导致内部积分不准表现是长时间推进后能量缓慢增长而不是瞬间发散。gauss_legendre 函数可以直接用常用数值教材里的递推版本精度到 1e-15 就够用。3.2 状态评估与半离散右端项在DG代码里把模态系数转换为边界值和积分点值是最高频的操作。单独封装一个函数可以避免到处索引出错function val eval_modal(ue, phi_at_xi) % ue: 长度为 deg1 的模态系数 % phi_at_xi: 该参考点上的基函数值行向量 val ue(:) * phi_at_xi(:); end注意这里 ue 统一用列向量存储phi 用行向量这样ue * phi就是标量。如果行向量和列向量反了MATLAB会隐式扩展直接得到一个矩阵导致维度爆炸。这个问题在第5章还会遇到。体积分在半离散右端项里是单元内通量和基函数导数的卷积。对线性DG两个高斯点就能精确积分但为了程序通用我保留 deg2 点。核心片段如下% 体积分sum_j w_j * f(q(xi_j)) * dphi_i/dx dqe_vol zeros(deg1, 2); for j 1 : nq p_j eval_modal(qe(:,1), phi(j,:)); v_j eval_modal(qe(:,2), phi(j,:)); f1 mat(1) * mat(2)^2 * v_j; f2 p_j / mat(1); dphi_dx dphi_dxi(j,:) / detJ; % 物理微分 dqe_vol(:,1) dqe_vol(:,1) wq(j) * f1 * dphi_dx; dqe_vol(:,2) dqe_vol(:,2) wq(j) * f2 * dphi_dx; end这里 dqe_vol 是本单元模态的时导。因为在参考坐标系里φ_i 对 x 的导数等于对 ξ 的导数除以雅可比。物理单元长度 h雅可比 detJh/2所以dphi_dx dphi_dxi / detJ。如果你忘了除以 detJ短单元会表现出过大的刚度和过小的有效波速数值结果会在几个周期后明显“发胖”。面积分部分会在循环外层先计算所有界面通量再分别加到左右单元上。LF通量函数已经在第2章定义这里不需要改。3.3 RK3时间推进与CFL空间离散完成后得到一组ODEdq/dt L(q)。显式时间推进用三阶强稳定RK3虽然比RK4多一次求值但对DG来说强稳定保号的特性更重要。标准步骤如下% q 是 [N, deg1, 2] 状态数组 % rhs 是半离散右端项函数句柄 q1 q dt * rhs(q); q2 0.75 * q 0.25 * q1 0.25 * dt * rhs(q1); q3 1/3 * q 2/3 * q2 2/3 * dt * rhs(q2); q q3;参数说明RK3的系数是强稳定型的适合双曲问题。dt 必须满足CFL条件一维线性声波的经验上限是dt CFL * h / (c * (2*deg1))CFL取0.4左右比较稳。比如 N100、deg1、c1500h0.01那么 dt 大约是 0.40.01/(15003)8.9e-7 秒。如果 dt 大于1e-6大概率跑几步就NaN。我通常先按解析解估算 dt再在时间循环里打印能量来确认。为了快速验证可以每10步记录一次全局能量sum(p.^2 / (rho*c^2) rho*v.^2)。无源情况下这个能量应逐步下降因为LF通量有耗散但不会跳跃。如果能量突然暴涨第一步查CFL第二步查边界通量第三步查雅可比。4. 二维声波模拟多介质界面与波场快照4.1 二维矩形网格上的张量积基直接从一维跳到二维最关键的区别是每个单元多了 y 方向自由度。如果使用矩形网格基函数可以取一维勒让德的张量积 φ_i(x) φ_j(y)积分节点也是高斯节点的乘积。单元上的质量矩阵仍然是块对角每个阶数组合一个对角块。我先把网格参数存成结构体% 二维矩形网格参数 Nx 100; Ny 100; deg 1; xmin 0; xmax 2; ymin 0; ymax 2; hx (xmax-xmin)/Nx; hy (ymax-ymin)/Ny; % 高斯积分节点一维 nq deg 2; [xiq, wq] gauss_legendre(nq); % 张量积积分权重 w2d wq(:) * wq(:); xi2d [kron(ones(nq,1), xiq), kron(xiq, ones(nq,1))];参数说明张量积网格的好处是基函数简单界面通量计算可以拆成 x 方向和 y 方向分别处理。代价是矩形网格在弯曲边界上需要阶梯近似对于光滑曲面的散射问题误差较大。如果研究对象是任意形状异常体建议换成三角形网格但三角形网格上正交基构造要麻烦一些至少要处理非对角的单元质量矩阵。4.2 多介质参数与界面通量处理二维声波方程的状态变量变成 [p; vx; vy]通量有两个方向F_x [ρc² vx; p/ρ; 0]F_y [ρc² vy; 0; p/ρ]界面处仍然用方向分解后的LF通量。不同介质在单元交界处可能参数不同通量计算必须用左右两侧的 ρ、c。我会预先把每个单元的介质编号存成二维矩阵界面通量循环时查这个矩阵% mediumIdx: [Nx, Ny] 整型矩阵1 或 2 mat1 [1.0; 1500.0]; % 水介质 rho,c mat2 [2.0; 2000.0]; % 沉积层 % 取界面两侧介质参数 function mat get_mat(idx) if idx 1 mat mat1; else mat mat2; end end注意多介质界面最忌讳只在单元中心取介质参数。声波方程的通量是边界问题界面上的物理通量应该由两侧状态和两侧介质参数共同决定。如果简单取平均密度或平均声速界面反射系数就完全错了。二维情况下物理通量函数可以写成function F physical_flux(q, rho, c, dir) % q [p; vx; vy] % dir 1 表示 x 方向2 表示 y 方向 switch dir case 1 F [rho*c^2*q(2); q(1)/rho; 0]; case 2 F [rho*c^2*q(3); 0; q(1)/rho]; end end界面法向通量同样用两侧 F 的平均减波速乘跳跃。先单独把每个方向算清楚再合到同一个更新循环里比写一个通用的向量通量函数更好调试。4.3 震源加载与波场可视化二维中常用Ricker子波作为震源。可以用一个高斯空间分布和Ricker时间函数相乘加到某个单元的模态系数上。如果是压力源就只在 p 的自由度上累加% Ricker 子波时间函数 function src ricker(t, f0) % f0 是主频 tau pi * f0 * (t - 1.0/f0); src (1 - 2*tau.^2) .* exp(-tau.^2); end % 在时间循环中往某个压力自由度加源 q(source_cell, 1, 1) q(source_cell, 1, 1) src(t) * dt;参数说明f0 要避开网格的最高解析频率。一般要求每个波长至少覆盖5个单元如果 f0 太高源本身就有接近网格色散关系的能量会产生尾波。加源时如果忘记乘 dt相当于提高子波幅度看不出格式错误但能量曲线会突然变得不合理。可视化用pcolor或surf显示压力场。注意 pcolor 输入格点值需要把模态系数映射到网格顶点直接用单元中心平均也可以% 取每个单元压力平均作为当前波场快照 Pmean zeros(Nx, Ny); for e 1 : Nx*Ny % 从结构体或数组中取该单元模态系数算零阶模态 Pmean(e) q_modal(e, 1, 1); end pcolor(xc2d, yc2d, Pmean); shading interp; axis equal; colorbar;如果你希望看到高次DG的细节用单元中心平均会抹平高频分量建议把每个单元内部的高斯点值也画出来用scatter而不是pcolor。我踩过这个坑一幅看起来很光滑的云图实际上已经丢掉了P1DG的高频信息。5. 避坑与常见问题频散、CFL与边界反射的五个排查点DG方法理论漂亮但代码里翻车的地方往往不在理论而在实现细节。下面五条是我在模拟项目X里真实踩过的坑每条都按现象、原因、解决给出照着排查能省下半天时间。5.1 高频振荡多项式阶数越高反而越毛糙现象使用P2或P3格式时光滑初始波场跑到几十步后出现明显的高频“毛刺”而且加密网格后毛刺依然存在只是位置变化。原因DG的色散关系在不同阶数下行为不同高阶格式在中高频段可能出现与物理不符的群速度另一个常见原因是内部积分用了等距点而不是高斯点导致伪高阶。解决先退回P1验证模型逻辑确认P1没问题再升阶。内部积分必须用高斯求积积分点数不小于 deg2。如果必须用高阶把时间推进换成低耗散的RK4并把CFL收紧到0.2左右。5.2 时间步长不满足CFL导致NaN现象代码运行几十步后压力场出现成片数值震荡随后能量变成NaN。原因显式时间推进稳定性条件不满足。一维下 dt 超过 h/(c*(2deg1))二维还要再乘方向数相关系数通常再除以2。解决统一用dt 0.4 * min(dx,dy) / (cmax * (2*deg1))。注意多介质里 cmax 要取全局最大声速而不是平均声速。我后来把CFL计算写成一个独立函数每次改网格都自动重算不再手动填。5.3 边界反射明显把通量边界当成零压力边界现象波传到计算域边缘后被完全反射反射波能量和入射波同级很难分辨目标区域的多次波。原因边界上通量没做处理用 p0 的Dirichlet边界在DG里等于强制边界外侧压力为零构成刚性反射面。解决最简单的方法是周期边界或使用大物理域让反射波在观测时间内到不了区域中心。若要做吸收边界可以走完美匹配层路线。本资源里一维验证用周期边界二维算例把边界扩展了足够的距离没有把零压力当吸收边界用。5.4 多介质界面反射系数偏大用单元平均参数算通量现象一个平面波垂直入射到两种介质界面反射波幅度比解析值高且入射波一侧出现拖尾振荡。原因通量函数里用了单元平均密度或平均声速而不是界面两侧真正介质参数界面处的面力矩不被满足。解决对每个界面先根据左右单元编号取出实际 ρ 和 c再将状态插值到界面同一点。垂直入射的解析反射系数是 (Z_right-Z_left)/(Z_rightZ_left)Zρc可以用这个公式快速检查通量实现是否正确。5.5 维度不匹配MATLAB隐式扩展把行向量变成矩阵现象eval_modal 函数明明应该返回标量却报错或者计算结果出一个 N×N 矩阵。原因模态系数用了行向量基函数也用了行向量ue(:) * phi没问题但如果 ue 是 n×1 数组phi 是 1×m 数组MATLAB会隐式扩展成矩阵后续全部炸掉。解决在函数入口统一加ue ue(:); phi phi(:).;保证内部一定是抖量。另外所有状态数组固定为[N, deg1, 2]不要混用cell或不同方向的行列排列。6. 让代码跑得再快一点向量化、并行与收敛阶验证资源里的一维版本为了可读性用了多重循环实际做二维模拟时建议逐步向量化。体积分循环其实是对所有单元同时做同一个线性变换MATLAB里可以直接用矩阵乘法替代。比如把每个单元的高斯点状态排成子块一次性算完所有单元的 dqe_vol速度能提升一个数量级。更省事的是用parfor替代时间层外的单元循环但要注意rhs函数里不能有全局变量依赖。验证代码是否写对比加速更重要。我每做完一个算例都会跑一段收敛阶测试对固定模拟时间比较不同网格步长 h 下的误差看误差随 h 的下降斜率是否接近 deg1 阶。下面是一维情况的检验思路% 使用解析解作为初值在 tEnd 与解析解比较 L2 误差 hList [0.1, 0.05, 0.025, 0.0125]; errList zeros(1, length(hList)); for k 1 : length(hList) % 重新生成网格与初值 [errList(k), ~] run_dg_1d(hList(k), deg, tEnd); end % 误差斜率 order diff(log10(errList)) / diff(log10(hList));如果 order 接近 deg1说明空间离散正确如果接近 deg先查面积分符号如果乱七八糟直接查雅可比和通量方向。这套验证流程对我排查代码问题比任何调试器都管用。从那以后我每次拿到别人的DG代码第一件事不是看主循环而是先跑一个已知解析解的收敛阶测试。这个习惯救了我很多次也希望帮到你。本文还有配套的精品资源点击获取