简介面向数值分析与偏微分方程求解方向的学习者这份资源聚焦二维声波方程的间断有限元DG数值实现将空间离散与三阶龙格库塔时间推进相结合适合已具备Matlab基础、希望深入理解间断伽辽金方法的中高级读者参考。压缩包共6个文件以5个m脚本为主辅以1个asv自动备份文件整体约4KB体量轻便便于直接阅读与二次修改。资源围绕矩形域上正弦初值条件下的声波传播问题涵盖网格划分、线性基函数选取、界面通量处理、边界条件设置与时间步迭代等核心环节代码结构对应空间离散与时间推进两条主线可帮助读者理清DG方法从公式到程序的落地路径。目前已有1567人学习下载适合作为课程设计、科研入门或算法复现的参考素材通过对照脚本理解数值稳定性与高阶精度的实现细节。1. 声波方程数值模拟为什么我最后选了间断有限元这套 MATLAB 实现第一次接触声波方程的数值求解是在一个地震波场正演的小项目里。当时用常规有限差分跑二维均匀介质没问题一换到复杂速度模型网格步长和稳定性条件就开始互相拉扯改一处崩一处。后来换成间断有限元Discontinuous Galerkin简称 DG思路才发现它在处理复杂介质、局部网格加密和非结构网格上确实顺手——单元之间只通过数值通量耦合每个单元内部可以独立选阶数这对声波这种双曲型方程特别友好。这份 MATLAB 实现就是围绕二维声波方程搭的一套间断有限元求解器从网格、基函数、数值通量到时间推进都有对应代码适合正在做波场模拟、想搞懂 DG 到底怎么落地的人也适合拿它当二次开发的起点。2. 项目结构与技术选型从目录看这套求解器怎么组织2.1 整体目录与模块划分拿到源码先别急着跑先把目录扫一遍。这套实现按数值方法的逻辑分层常见做法是拆成下面几块dg_acoustic/ ├── main.m % 入口脚本参数配置、调用求解、出图 ├── mesh/ │ ├── build_mesh.m % 生成结构化/非结构网格 │ └── mesh_info.m % 单元邻接、面信息、法向 ├── basis/ │ ├── legendre_basis.m % 一维 Legendre 基函数 │ ├── basis_2d.m % 二维张量积基函数 │ └── vandermonde.m % Vandermonde 矩阵与节点 ├── solver/ │ ├── assemble_rhs.m % 组装右端项体积分面积分 │ ├── numerical_flux.m % 数值通量迎风/Lax-Friedrichs │ └── time_step.m % 时间推进RK ├── postprocess/ │ └── plot_wavefield.m % 波场快照与动画 └── utils/ └── gauss_quad.m % 高斯积分节点与权重这个划分不是随便摆的它对应 DG 求解的四个核心动作建网格、选基函数、算通量、推时间。新手最容易犯的错是把所有逻辑塞进一个脚本改一个参数要翻半天这套结构至少让你知道改数值通量该去numerical_flux.m改时间格式该去time_step.m。2.2 技术选型为什么是 MATLAB Legendre 基 显式 RKMATLAB 做这套东西的优势在于矩阵运算和可视化原型阶段改公式、看波场都快。基函数选 Legendre 多项式是因为它在参考单元[-1,1]上正交质量矩阵接近对角求逆成本低二维用张量积构造basis_2d.m里就是两个一维基的外积。时间推进用显式 Runge-Kutta因为 DG 半离散后是一个常微分方程组dU/dt L(U)显式 RK 实现简单配合 CFL 条件控制步长即可。提示显式时间推进对 CFL 数敏感网格越小、阶数越高允许的时间步越短。跑之前先确认dt和h、p的关系否则波场会直接发散。参数配置一般集中在main.m顶部常见字段如下参数含义典型取值p多项式阶数14nx, ny网格数50200c声速15004000T总模拟时间0.52.0cflCFL 数0.10.3flux_type通量类型upwind/LF3. 核心模块拆解网格、基函数、通量与时间推进3.1 网格构建与单元信息DG 的网格不仅要给出节点坐标还要给出每个单元的面、法向和邻接关系。build_mesh.m负责生成规则矩形网格mesh_info.m把单元-面-邻接整理成结构体。核心是给每个面存左右单元编号和单位外法向后面算面积分时直接查表。function mesh build_mesh(nx, ny, Lx, Ly) % 生成 nx×ny 的结构化矩形网格 % 返回单元顶点坐标、单元中心、面信息 dx Lx / nx; dy Ly / ny; [X, Y] meshgrid(linspace(0, Lx, nx1), linspace(0, Ly, ny1)); mesh.dx dx; mesh.dy dy; mesh.nx nx; mesh.ny ny; % 单元中心坐标 mesh.xc X(1:end-1,1:end-1) dx/2; mesh.yc Y(1:end-1,1:end-1) dy/2; % 面法向左右面 (±1,0)上下面 (0,±1) mesh.normals [ -1 0; 1 0; 0 -1; 0 1 ]; end这段代码只搭了几何骨架真正让 DG 跑起来的是mesh_info.m里补的邻接表。每个内部面要记录左右单元索引边界面的外侧单元记为 0后面通量计算遇到 0 就走边界条件分支。参数nx, ny决定分辨率Lx, Ly是物理域尺寸改分辨率时记得同步检查dt因为dx变小会压缩稳定步长。3.2 Legendre 基函数与二维张量积基函数模块是整个求解器的地基。一维 Legendre 基在参考单元上定义legendre_basis.m返回基函数值和导数二维通过张量积组合basis_2d.m输出每个单元内所有基函数的取值矩阵。function [phi, dphi] legendre_basis(p, xi) % 计算 0~p 阶 Legendre 基函数及其导数 % xi 为参考单元 [-1,1] 上的积分点 phi zeros(length(xi), p1); dphi zeros(length(xi), p1); phi(:,1) 1; dphi(:,1) 0; if p 1 phi(:,2) xi; dphi(:,2) 1; end for k 2:p % 三项递推公式 phi(:,k1) ((2*k-1)*xi.*phi(:,k) - (k-1)*phi(:,k-1)) / k; dphi(:,k1) dphi(:,k-1) (2*k-1)*phi(:,k); end end递推公式是 Legendre 多项式的标准构造p控制阶数xi是高斯积分点。阶数越高单单元内能表达的波场细节越多但计算量和存储也上去。实际项目里p2或p3是常见折中再高就要考虑质量矩阵求逆的代价。二维版本在basis_2d.m里对两个方向做外积得到(p1)^2个基函数。3.3 数值通量与右端项组装DG 和连续有限元最大的区别就在面通量。numerical_flux.m实现迎风格式或 Lax-Friedrichs 通量assemble_rhs.m把体积分和面积分拼成右端项。function F numerical_flux(uL, uR, n, c, type) % 计算面上的数值通量 % uL/uR: 左右单元在面上的值n: 外法向c: 声速 if strcmp(type, upwind) % 迎风按特征方向取上游值 vn c * (n(:,1) n(:,2)); % 简化示意 F 0.5*(c*n).*(uL uR) 0.5*abs(vn).*(uL - uR); elseif strcmp(type, LF) % Lax-Friedrichs加数值耗散 alpha max(abs(c)); F 0.5*(c*n).*(uL uR) - 0.5*alpha*(uR - uL); end end迎风格式精度高但依赖特征方向Lax-Friedrichs 更稳但耗散大。assemble_rhs.m里对每个单元循环先算体积分∫ φ_i * f(u) dx再遍历四个面算面积分∫ φ_i * F * ds最后相加。这里最容易翻车的是法向方向搞反导致左右单元通量符号不一致波场会出现非物理振荡。3.4 时间推进与稳定性控制半离散后得到dU/dt L(U)time_step.m用显式 RK 推进。常见做法是低存储 RK 或经典四阶 RK前者省内存后者精度高。function U time_step(U, dt, rhs, order) % 显式 Runge-Kutta 时间推进 switch order case 1 % 前向 Euler U U dt * rhs(U); case 2 % 中点法 k1 rhs(U); k2 rhs(U 0.5*dt*k1); U U dt * k2; case 4 % 经典 RK4 k1 rhs(U); k2 rhs(U 0.5*dt*k1); k3 rhs(U 0.5*dt*k2); k4 rhs(U dt*k3); U U dt/6 * (k1 2*k2 2*k3 k4); end enddt由 CFL 条件给出常见估计是dt cfl * h / (c * (2p1))h取最小网格尺度。order选 2 或 4 看精度需求原型阶段 RK2 够用出图前再换 RK4。跑的时候如果波场突然爆掉先查dt是不是超了稳定上限再查通量法向。4. 跑起来从零到第一张波场快照4.1 环境准备与入口脚本MATLAB 版本建议 R2018b 以上用到的都是基础矩阵运算和meshgrid、imagesc这类函数不需要额外工具箱。把整个目录加到路径打开main.m先改顶部参数块% main.m 顶部参数 p 2; % 多项式阶数 nx 100; ny 100; c 2000; % 声速 m/s T 0.8; % 总时间 s cfl 0.2; flux LF; % 先用 Lax-Friedrichs 保稳 src [0.5, 0.5]; % 震源位置归一化坐标第一次跑建议就用这套保守参数低阶、中等网格、LF 通量。确认能出图后再逐步加阶数、换迎风、加密网格。4.2 执行流程与出图main.m的执行顺序是建网格 → 初始化波场和震源 → 循环时间步 → 定期存快照 → 出图。核心循环大致如下mesh build_mesh(nx, ny, 1, 1); mesh mesh_info(mesh); U zeros((p1)^2, nx*ny); % 每个单元 (p1)^2 个自由度 t 0; it 0; while t T dt cfl * min(mesh.dx, mesh.dy) / (c*(2*p1)); U time_step(U, dt, (u) assemble_rhs(u, mesh, p, c, flux), 2); t t dt; it it 1; if mod(it, 20) 0 plot_wavefield(U, mesh, p); drawnow; end endassemble_rhs用匿名函数包一层把网格、阶数、声速、通量类型都传进去这样time_step不用关心物理参数。mod(it,20)控制出图频率太频繁会拖慢运行。第一张图出来如果看到清晰的圆形波前说明网格、通量、时间步基本对上了。4.3 结果验证均匀介质对照解析解跑通不等于跑对。均匀介质下声波方程有解析解常见做法是在同一点接收把数值解和解析解叠在一张图上比波形。如果相位偏移明显多半是时间步或阶数不够如果幅值差一个量级查一下震源归一化和积分权重。这个对照步骤别省后面改复杂介质时它就是你的基准。5. 避坑与排查这套 DG 代码最容易翻车的几个地方5.1 波场几步就发散现象跑几十步后数值爆炸NaN或Inf出现。原因多半是dt超过显式格式稳定上限或者通量法向符号写反导致能量注入。解决先把cfl降到 0.1确认能稳跑再逐面打印法向和左右单元值检查内部面两侧通量是否满足守恒。5.2 波前出现锯齿或非物理振荡现象波前不光滑高频毛刺明显。原因通常是阶数低但网格粗或者迎风通量的特征方向判断有误。解决先换LF通量看是否消失若消失则是迎风实现问题若仍在提高p或加密nx, ny同时按新h重算dt。5.3 边界处波场反射异常现象边界附近出现不该有的强反射。原因边界面的外侧单元索引为 0通量分支没处理好等价于硬壁反射。解决在numerical_flux.m里对边界单独走吸收边界或特征边界条件别让 0 号单元参与正常通量计算。5.4 改阶数后结果对不上现象p从 2 改到 3波场形态明显变化。原因基函数阶数变了但积分点数量、质量矩阵、dt估计没同步更新。解决把p、积分点数、dt公式绑在一起改改完先用均匀介质对照解析解验证。5.5 运行慢到无法调试现象网格一加密单步时间暴涨。原因assemble_rhs里用了多层循环逐单元逐面计算MATLAB 循环开销大。解决把单元循环向量化面通量用索引矩阵批量算或者先用粗网格调逻辑最后再上细网格。6. 二次开发把这套求解器改成你自己的波场工具这套代码的价值不在跑出一个漂亮波前而在它把 DG 的每个环节都摊开给你看。二次开发我一般从三个方向下手。第一是换物理声波方程改成弹性波或声-弹耦合核心改动在通量的特征分解和变量个数numerical_flux.m里把标量通量换成矩阵形式assemble_rhs的自由度维度跟着变。第二是换网格把结构化矩形换成三角形非结构网格build_mesh和mesh_info要重写但基函数、通量、时间推进的框架不动这也是 DG 相对有限差分的好处。第三是加吸收边界在边界面上实现 PML 或特征边界改的是边界分支不影响内部单元。验证改造是否正确的习惯我一般固定三步先用均匀介质对照解析解看波形和相位再把网格加密一档看结果是否收敛最后换两种通量类型看解对通量选择的敏感度是否在合理范围。这三步走完基本能排除大部分实现错误。改造方向主要改动文件验证重点弹性波numerical_flux.m、assemble_rhs.m特征速度、变量维度非结构网格build_mesh.m、mesh_info.m面法向、邻接一致性吸收边界numerical_flux.m边界分支边界反射强度有个血泪教训早期我改完通量直接上复杂模型结果波场乱了半天找不到原因最后发现是边界分支里法向没归一化。从那以后我每次改完核心模块都强制先跑一遍均匀介质对照再动复杂模型。这套 MATLAB 实现适合当你的起点把框架吃透之后改起来比从零搭快得多。希望帮到你。本文还有配套的精品资源点击获取
