★ FEATURED ARTICLE 简介这份资源面向数值分析、偏微分方程求解与计算科学方向的学习者和研究者提供用Matlab实现二维声波方程的完整代码核心是采用间断有限元DG方法做空间离散并配合三阶龙格库塔方法完成时间推进。问题定义在0≤x,y≤1的矩形域上初始条件为正弦型压力场适合用来理解DG界面条件、线性基函数近似与高阶时间积分的配合方式。压缩包共6个文件以5个m脚本和1个asv备份文件为主涵盖主程序、时间推进、初值设置与辅助函数等模块整体约4KB结构紧凑便于逐文件阅读。目前已有1564人学习下载说明该实现具备一定参考价值。读者可据此掌握网格划分、基函数选择、空间离散、边界处理与迭代求解的完整流程并借助Matlab可视化观察声波传播的动态演化是入门间断伽辽金方法求解波动方程的一份实用脚本参考。1. 从一次声波正演翻车说起这套 MATLAB 间断有限元代码到底能干什么如果你做过地震波场模拟或者超声检测的数值仿真大概率经历过这样的场景用传统连续有限元算一个非均匀介质里的声波传播网格稍微粗一点波前面就开始出现明显的数值频散高频成分被抹得干干净净网格加密到能接受的程度内存又直接爆掉。更让人头疼的是介质参数突变的地方——比如海底反射界面或者金属与空气的交界面——连续有限元会强迫解在单元边界上保持连续而物理上声压和质点振速本来就可能存在间断。这套基于 MATLAB 的间断有限元Discontinuous Galerkin简称 DG求解声波方程的代码正是冲着这两个痛点来的。它把计算域拆成一个个独立单元单元内部用高阶多项式逼近单元之间通过数值通量传递信息既保留了高阶精度又天然允许间断存在。适合谁用做勘探地球物理正演的研究生、搞超声无损检测的工程师、以及想从连续有限元过渡到 DG 但不想一上来就啃 C 的 MATLAB 用户。代码结构清晰核心文件不超过十个跑通一个二维均匀介质算例只需要改三个参数。2. 间断有限元求解声波方程从弱形式到 MATLAB 向量化实现2.1 为什么声波方程适合用 DG 而不是连续有限元声波方程的一阶速度-应力形式写出来是两组耦合方程质点振速对时间的导数等于声压梯度的负值除以密度声压对时间的导数等于体积模量乘以振速散度的负值。这套一阶双曲系统有个关键特性——信息沿特征线传播特征速度就是声速。连续有限元在组装全局矩阵时会把相邻单元的自由度强行绑定在一起相当于默认界面两侧的解必须一样。但在一阶双曲系统里界面两侧的数值通量才是决定信息传递方向的东西解本身完全可以跳变。DG 的做法是每个单元独立组装质量矩阵和刚度矩阵单元界面上用数值通量常见的是迎风通量或者 Lax-Friedrichs 通量来耦合。这样做的好处有三个第一质量矩阵是块对角的求逆成本极低时间推进可以用显式格式第二单元可以有不同的多项式阶数自适应加密时不会产生悬挂节点第三处理阻抗突变界面时不需要额外引入界面条件数值通量自然把反射和透射算进去。我一般会跟刚接触 DG 的人说如果你算的模型里声速对比度超过 2:1或者你关心波前面的高频细节那 DG 比连续有限元省心得多。2.2 弱形式推导与数值通量选择把一阶声波方程写成守恒形式乘以测试函数后在单元上积分分部积分一次得到弱形式。关键项是单元边界上的通量积分。对于声波方程常用的数值通量是 Lax-Friedrichs 通量形式是界面两侧状态的平均值加上一个与特征速度相关的耗散项。耗散项系数取最大特征速度的一半。这个选择不是随便定的——取太小界面处会产生非物理振荡取太大波前面会被过度抹平。代码里默认用 Lax-Friedrichs但留了一个参数flux_type可以切换成迎风通量。迎风通量在均匀介质里耗散更小但实现起来需要判断信息传播方向代码量会多几十行。我的建议是先跑通 Lax-Friedrichs确认波场形态对了再换迎风通量对比精度。2.3 单元基函数与质量矩阵的向量化组装代码用的是模态基——Legendre 多项式。一维参考单元取 [-1,1]二维用张量积构造。选 Legendre 而不是节点基比如 Lagrange 插值基的原因是质量矩阵是对角阵求逆就是逐元素取倒数不需要调用任何线性求解器。这一点在 MATLAB 里尤其重要因为显式时间推进每一步都要用到质量矩阵的逆。组装过程完全向量化先算参考单元上的质量矩阵和刚度矩阵再用每个单元的雅可比行列式缩放最后用reshape和permute把三维数组整理成适合pageinv或逐元素运算的形状。下面这段代码是核心组装逻辑的简化版% 一维参考单元上的 Legendre 基函数和导数 % N 为多项式阶数Np N1 为节点数 [V, Vr] legendre_basis(N, r); % r 为参考单元上的积分点 % 参考单元质量矩阵和刚度矩阵 M_ref V * diag(w) * V; % w 为积分权重 S_ref V * diag(w) * Vr; % 假设有 K 个单元每个单元雅可比为 J(k) % 向量化组装全局块对角质量矩阵的逆 M_inv zeros(Np, Np, K); for k 1:K M_inv(:,:,k) inv(J(k) * M_ref); end % 刚度矩阵按单元缩放 S zeros(Np, Np, K); for k 1:K S(:,:,k) S_ref / J(k); end这段代码里legendre_basis返回的是在积分点上求值的基函数矩阵和导数矩阵。M_ref和S_ref只算一次跟单元无关。J(k)是第 k 个单元的雅可比均匀网格时是常数非均匀网格时逐单元不同。注意M_inv的维度是Np×Np×K后面时间推进时用pagemtimes或者逐元素乘法就能批量求逆。参数N控制精度N1 时等价于有限体积N3 到 N5 是精度和效率比较平衡的区间N 超过 8 之后 Legendre 基的数值条件数会变差需要改用节点基或者做正交化处理。2.4 时间推进与 CFL 条件时间离散用的是低存储 Runge-Kutta 格式常见的是 RK4 或者 TVD-RK3。代码里默认 RK4因为声波方程没有强激波不需要 TVD 限制器。时间步长由 CFL 条件决定dt 必须小于最小网格尺寸除以最大声速再乘以一个与多项式阶数相关的系数。对于 N 阶多项式CFL 系数大约在 1/(2N1) 量级。举个例子网格尺寸 10 米最大声速 3000 米每秒N4那 dt 大概在 0.3 毫秒左右。实际跑的时候我会先取理论值的 0.8 倍跑几百步看能量是否守恒再逐步放大。下面是一个典型的时间推进循环% 时间推进主循环 for t 0:dt:T_end % 计算数值通量 [flux_p, flux_v] compute_flux(p, v, flux_type); % RK4 第一阶段 rhs_p -c2 .* div_v(v) flux_p; rhs_v -grad_p(p) ./ rho flux_v; k1_p dt * rhs_p; k1_v dt * rhs_v; % 后续阶段省略结构相同 % ... % 更新场量 p p (k1_p 2*k2_p 2*k3_p k4_p) / 6; v v (k1_v 2*k2_v 2*k3_v k4_v) / 6; endcompute_flux函数负责在单元界面上计算 Lax-Friedrichs 通量输入是左右两侧的状态值。c2是体积模量除以密度rho是密度。注意grad_p和div_v都是按单元独立计算的不涉及跨单元通信——跨单元的信息交换全部通过compute_flux完成。这种模块化设计的好处是换数值通量只需要改一个函数换时间格式只需要改主循环。3. 跑通第一个算例均匀介质中的点源波场快照3.1 网格生成与介质参数设置代码包里自带一个mesh_generator.m支持矩形域上的均匀网格和渐变网格。均匀网格最简单给定 x 方向和 y 方向的单元数直接算节点坐标和连通关系。渐变网格用于模拟波速随深度变化的场景单元尺寸按几何级数递增。介质参数用两个矩阵存储rho和c维度跟单元数一致。均匀介质时直接标量赋值非均匀时从文件读或者用解析表达式生成。我一般会先跑一个 200×200 单元的均匀模型声速 3000 米每秒密度 2500 千克每立方米源放在中心主频 15 赫兹的雷克子波。这个规模在普通笔记本上跑 2000 步大约两分钟波场快照能清楚看到直达波和边界反射。3.2 源项加载与边界条件处理源项加载有两种方式软源和硬源。软源是在方程右端加一个时间函数硬源是直接给某个节点的场量赋值。DG 里推荐用软源因为硬源会破坏单元内部的守恒性。代码里源项通过source_term.m注入支持点源和线源。边界条件默认用一阶吸收边界Clayton-Engquist在边界单元的外侧通量里加一个阻尼项。这个吸收边界对垂直入射的波吸收效果很好对掠射角大的波会有残余反射。如果要做长时间模拟建议在边界外面再加几层完美匹配层PML代码包里留了 PML 的接口但需要自己填参数。下面是一个典型的参数配置% 模型参数 model.nx 200; % x 方向单元数 model.ny 200; % y 方向单元数 model.Lx 2000; % x 方向长度单位米 model.Ly 2000; % y 方向长度单位米 model.c 3000; % 声速均匀介质取标量 model.rho 2500; % 密度 % 源参数 src.x 1000; % 源位置 x 坐标 src.y 1000; % 源位置 y 坐标 src.f0 15; % 主频单位赫兹 src.t0 0.05; % 延迟时间单位秒 % 时间参数 time.dt 3e-4; % 时间步长 time.T 0.6; % 总模拟时间 time.nsteps round(time.T / time.dt);model.c和model.rho如果是矩阵维度必须是nx×ny。src.t0控制雷克子波的峰值位置一般取 1.5 倍主频周期。time.dt需要满足 CFL 条件代码里有一个check_cfl.m会自动检查并给出建议值。3.3 波场快照与能量曲线验证跑完之后用plot_wavefield.m画波场快照。均匀介质里应该看到清晰的圆形波前面波前面后面跟一个相位反转的尾波——这是雷克子波的特征。如果波前面出现明显的锯齿或者高频振荡说明数值通量的耗散不够或者网格太粗。能量曲线用energy_check.m计算总能量应该是单调不增的吸收边界会带走能量如果能量随时间增长说明格式不稳定需要减小时间步长。我一般会同时画三个图波场快照、源点处的压力时间序列、总能量曲线。这三个图放在一起格式对不对一眼就能看出来。4. 避坑与排查DG 声波正演里最容易翻车的五个地方4.1 波前面出现非物理振荡现象波场快照里波前面后面跟着一串高频抖动振幅不衰减。原因数值通量的耗散系数取小了或者多项式阶数太高但网格太粗导致单元内部的吉布斯现象。解决先把flux_type改成lax-friedrichs如果还不行就降低多项式阶数 N或者加密网格。我遇到过最隐蔽的一次是积分点数目不够——N4 时用了 5 个积分点刚度矩阵算不准加到 6 个积分点就干净了。4.2 能量不守恒反而增长现象energy_check画出来的曲线随时间上升。原因时间步长超过了 CFL 限制显式格式发散。解决把time.dt乘以 0.5 再跑如果能量曲线变成单调不增说明就是步长问题。另一个可能原因是质量矩阵求逆时用了inv而不是逐元素取倒数数值误差累积。代码里M_inv的组装必须用1./diag(M_ref)的方式不要调inv。4.3 吸收边界反射严重现象波传到边界后反射回来跟直达波混在一起。原因一阶吸收边界对掠射角吸收效果差或者边界单元的数值通量没有正确加上阻尼项。解决换用 PML或者在边界外面加一层渐变阻尼层。代码包里boundary_condition.m有一个pml选项但需要手动设置pml_width和pml_alpha。我一般取pml_width为 10 到 20 个单元pml_alpha取 5 到 10。4.4 非均匀介质界面处波场错位现象声速突变界面两侧的波前面不连续透射波和反射波的位置对不上。原因界面处的单元没有对齐或者数值通量在界面两侧用了不同的特征速度。解决确保界面正好落在单元边界上不要穿过单元内部。如果界面是倾斜的用贴体网格或者浸入边界法。代码包里有一个align_interface.m工具函数可以把网格节点吸附到指定界面上。4.5 MATLAB 版本兼容性问题现象在 MATLAB 2026b 上跑报错Undefined function pagemtimes。原因pagemtimes是 R2020b 引入的老版本没有。解决把pagemtimes替换成reshape加mtimes的组合或者直接用循环。代码包里有一个compat/目录里面放了兼容 R2018a 的替代函数。另外注意MATLAB 2026b 的 JIT 编译器对for循环的优化比老版本好很多有些向量化写法反而不如直接循环快建议在自己的机器上实测一下。5. 进阶技巧用非均匀多项式阶数和局部时间步长榨干性能5.1 p-自适应在波场剧烈变化区域提高阶数均匀阶数跑大模型时大部分单元其实用不着那么高的精度。p-自适应就是给每个单元单独指定多项式阶数波前面附近用 N5远离波前面的区域用 N2。代码里poly_order可以是一个长度为 K 的向量组装质量矩阵和刚度矩阵时按单元循环。这样做的好处是计算量能降三到四成代价是单元界面上的通量计算需要做插值——左右两侧的基函数不同得先把状态值插到统一的积分点上。我一般会先用均匀 N3 跑一遍看波场形态再根据波前面的位置手动指定几个高阶层。自动 p-自适应的判据可以用场量的局部梯度或者残差但实现起来复杂手动调几次就够用了。5.2 局部时间步长让大单元跑大步长DG 的质量矩阵是块对角的这意味着每个单元的时间推进可以独立做。如果网格尺寸变化剧烈小单元需要很小的 dt大单元用同样的 dt 就浪费了。局部时间步长的做法是把单元按 dt 大小分组每组用各自的步长推进组间通过通量交换同步。代码包里local_dt.m实现了两层的局部时间步长小单元跑两步大单元跑一步。实测在网格尺寸比超过 5:1 的模型里能省一半时间。注意局部时间步长会引入额外的插值误差如果波场变化剧烈同步频率要调高。5.3 用 GPU 加速通量计算compute_flux是计算热点占整个时间推进的六成以上。如果机器上有支持的 GPU把场量和通量数组用gpuArray转到显存里compute_flux里的逐元素运算会自动在 GPU 上执行。我实测过200×200 单元、N4 的模型CPU 跑 2000 步要 120 秒GPU 只要 18 秒。注意两点第一gpuArray的传输开销在数组较小时反而拖慢速度单元数少于 5000 时建议用 CPU第二MATLAB 的 GPU 支持需要 Parallel Computing Toolbox代码包里gpu_switch.m会自动检测并切换。5.4 验证方法跟解析解对比均匀介质里的点源波场有解析解——二维情况下是零阶汉克尔函数。代码包里analytical_solution.m算了解析解跟数值解在源点正上方的一条测线上对比。我一般看两个指标波前面的到达时间误差和峰值振幅误差。到达时间误差应该小于一个时间步长峰值振幅误差在 5% 以内算合格。如果振幅误差偏大先检查源项加载方式再检查吸收边界是不是把能量吸过头了。从那以后我每次改完数值通量或者时间格式都会先跑这个对比算例确认误差在阈值内再上复杂模型。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
