简介：这份资源面向数值分析、偏微分方程求解与计算科学方向的学习者和研究者，提供用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 = N+1 为节点数 [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控制精度：N=1 时等价于有限体积，N=3 到 N=5 是精度和效率比较平衡的区间，N 超过 8 之后 Legendre 基的数值条件数会变差，需要改用节点基或者做正交化处理。

2.4 时间推进与 CFL 条件

时间离散用的是低存储 Runge-Kutta 格式，常见的是 RK4 或者 TVD-RK3。代码里默认 RK4，因为声波方程没有强激波，不需要 TVD 限制器。时间步长由 CFL 条件决定：dt 必须小于最小网格尺寸除以最大声速再乘以一个与多项式阶数相关的系数。对于 N 阶多项式，CFL 系数大约在 1/(2N+1) 量级。举个例子：网格尺寸 10 米，最大声速 3000 米每秒，N=4，那 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; end

compute_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，或者加密网格。我遇到过最隐蔽的一次是积分点数目不够——N=4 时用了 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-自适应就是给每个单元单独指定多项式阶数：波前面附近用 N=5，远离波前面的区域用 N=2。代码里poly_order可以是一个长度为 K 的向量，组装质量矩阵和刚度矩阵时按单元循环。这样做的好处是计算量能降三到四成，代价是单元界面上的通量计算需要做插值——左右两侧的基函数不同，得先把状态值插到统一的积分点上。我一般会先用均匀 N=3 跑一遍，看波场形态，再根据波前面的位置手动指定几个高阶层。自动 p-自适应的判据可以用场量的局部梯度或者残差，但实现起来复杂，手动调几次就够用了。

5.2 局部时间步长：让大单元跑大步长

DG 的质量矩阵是块对角的，这意味着每个单元的时间推进可以独立做。如果网格尺寸变化剧烈，小单元需要很小的 dt，大单元用同样的 dt 就浪费了。局部时间步长的做法是：把单元按 dt 大小分组，每组用各自的步长推进，组间通过通量交换同步。代码包里local_dt.m实现了两层的局部时间步长：小单元跑两步，大单元跑一步。实测在网格尺寸比超过 5:1 的模型里，能省一半时间。注意局部时间步长会引入额外的插值误差，如果波场变化剧烈，同步频率要调高。

5.3 用 GPU 加速通量计算

compute_flux是计算热点，占整个时间推进的六成以上。如果机器上有支持的 GPU，把场量和通量数组用gpuArray转到显存里，compute_flux里的逐元素运算会自动在 GPU 上执行。我实测过：200×200 单元、N=4 的模型，CPU 跑 2000 步要 120 秒，GPU 只要 18 秒。注意两点：第一，gpuArray的传输开销在数组较小时反而拖慢速度，单元数少于 5000 时建议用 CPU；第二，MATLAB 的 GPU 支持需要 Parallel Computing Toolbox，代码包里gpu_switch.m会自动检测并切换。

5.4 验证方法：跟解析解对比

均匀介质里的点源波场有解析解——二维情况下是零阶汉克尔函数。代码包里analytical_solution.m算了解析解，跟数值解在源点正上方的一条测线上对比。我一般看两个指标：波前面的到达时间误差和峰值振幅误差。到达时间误差应该小于一个时间步长，峰值振幅误差在 5% 以内算合格。如果振幅误差偏大，先检查源项加载方式，再检查吸收边界是不是把能量吸过头了。从那以后我每次改完数值通量或者时间格式，都会先跑这个对比算例，确认误差在阈值内再上复杂模型。希望帮到你。

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