企业官网建设流程全解析

简介：这是一份基于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) = Σ_{j=0}^{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, deg+1); dphi_dxi = zeros(num, deg+1); % 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(:,k+1) = ((2*k-1) * xi .* phi(:,k) ... - (k-1) * phi(:,k-1)) / k; dphi_dxi(:,k+1) = phi(:,k-1) ... + xi .* dphi_dxi(:,k) ... - (k-1)/k * dphi_dxi(:,k-1); end end

参数说明：deg=0 是分片常数格式，太耗散，一般不用；deg=1 是线性DG；deg=3 以上对光滑声波场能明显压制频散，但CFL步长会按 (2deg+1) 的比例收紧。注意递推公式里 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, N+1); % N+1个网格顶点 % 高斯积分点与权重，nq = deg+2 保证 2*deg+1 阶精确 nq = deg + 2; [xiq, wq] = gauss_legendre(nq); % 每个单元的物理坐标映射 xq = zeros(nq, N); for e = 1 : N xL_e = xNodes(e); xR_e = xNodes(e+1); xq(:,e) = xL_e + 0.5 * (xiq + 1) * h; end % 体积转化系数 detJ = h / 2;

参数说明：高斯积分节点阶数至少要比被积多项式阶数高，这里取 deg+2 点；点数过少会导致内部积分不准，表现是长时间推进后能量缓慢增长，而不是瞬间发散。gauss_legendre 函数可以直接用常用数值教材里的递推版本，精度到 1e-15 就够用。

3.2 状态评估与半离散右端项

在DG代码里，把模态系数转换为边界值和积分点值是最高频的操作。单独封装一个函数可以避免到处索引出错：

function val = eval_modal(ue, phi_at_xi) % ue: 长度为 deg+1 的模态系数 % phi_at_xi: 该参考点上的基函数值，行向量 val = ue(:)' * phi_at_xi(:); end

注意这里 ue 统一用列向量存储，phi 用行向量，这样ue' * phi就是标量。如果行向量和列向量反了，MATLAB会隐式扩展，直接得到一个矩阵，导致维度爆炸。这个问题在第5章还会遇到。

体积分在半离散右端项里是单元内通量和基函数导数的卷积。对线性DG，两个高斯点就能精确积分，但为了程序通用，我保留 deg+2 点。核心片段如下：

% 体积分：sum_j w_j * f(q(xi_j)) * dphi_i/dx dqe_vol = zeros(deg+1, 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，雅可比 detJ=h/2，所以dphi_dx = dphi_dxi / detJ。如果你忘了除以 detJ，短单元会表现出过大的刚度和过小的有效波速，数值结果会在几个周期后明显“发胖”。

面积分部分会在循环外层先计算所有界面通量，再分别加到左右单元上。LF通量函数已经在第2章定义，这里不需要改。

3.3 RK3时间推进与CFL

空间离散完成后，得到一组ODE：dq/dt = L(q)。显式时间推进用三阶强稳定RK3，虽然比RK4多一次求值，但对DG来说强稳定保号的特性更重要。标准步骤如下：

% q 是 [N, deg+1, 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*deg+1))，CFL取0.4左右比较稳。比如 N=100、deg=1、c=1500，h=0.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没问题再升阶。内部积分必须用高斯求积，积分点数不小于 deg+2。如果必须用高阶，把时间推进换成低耗散的RK4，并把CFL收紧到0.2左右。

5.2 时间步长不满足CFL导致NaN

现象：代码运行几十步后压力场出现成片数值震荡，随后能量变成NaN。

原因：显式时间推进稳定性条件不满足。一维下 dt 超过 h/(c*(2deg+1))；二维还要再乘方向数相关系数，通常再除以2。

解决：统一用dt = 0.4 * min(dx,dy) / (cmax * (2*deg+1))。注意多介质里 cmax 要取全局最大声速，而不是平均声速。我后来把CFL计算写成一个独立函数，每次改网格都自动重算，不再手动填。

5.3 边界反射明显：把通量边界当成零压力边界

现象：波传到计算域边缘后被完全反射，反射波能量和入射波同级，很难分辨目标区域的多次波。

原因：边界上通量没做处理，用 p=0 的Dirichlet边界，在DG里等于强制边界外侧压力为零，构成刚性反射面。

解决：最简单的方法是周期边界或使用大物理域，让反射波在观测时间内到不了区域中心。若要做吸收边界，可以走完美匹配层路线。本资源里一维验证用周期边界，二维算例把边界扩展了足够的距离，没有把零压力当吸收边界用。

5.4 多介质界面反射系数偏大：用单元平均参数算通量

现象：一个平面波垂直入射到两种介质界面，反射波幅度比解析值高，且入射波一侧出现拖尾振荡。

原因：通量函数里用了单元平均密度或平均声速，而不是界面两侧真正介质参数；界面处的面力矩不被满足。

解决：对每个界面，先根据左右单元编号取出实际 ρ 和 c，再将状态插值到界面同一点。垂直入射的解析反射系数是 (Z_right-Z_left)/(Z_right+Z_left)，Z=ρc，可以用这个公式快速检查通量实现是否正确。

5.5 维度不匹配：MATLAB隐式扩展把行向量变成矩阵

现象：eval_modal 函数明明应该返回标量，却报错，或者计算结果出一个 N×N 矩阵。

原因：模态系数用了行向量，基函数也用了行向量，ue(:)' * phi没问题；但如果 ue 是 n×1 数组，phi 是 1×m 数组，MATLAB会隐式扩展成矩阵，后续全部炸掉。

解决：在函数入口统一加ue = ue(:); phi = phi(:).';，保证内部一定是抖量。另外所有状态数组固定为[N, deg+1, 2]，不要混用cell或不同方向的行列排列。

6. 让代码跑得再快一点：向量化、并行与收敛阶验证

资源里的一维版本为了可读性用了多重循环，实际做二维模拟时建议逐步向量化。体积分循环其实是对所有单元同时做同一个线性变换，MATLAB里可以直接用矩阵乘法替代。比如把每个单元的高斯点状态排成子块，一次性算完所有单元的 dqe_vol，速度能提升一个数量级。更省事的是用parfor替代时间层外的单元循环，但要注意rhs函数里不能有全局变量依赖。

验证代码是否写对，比加速更重要。我每做完一个算例都会跑一段收敛阶测试：对固定模拟时间，比较不同网格步长 h 下的误差，看误差随 h 的下降斜率是否接近 deg+1 阶。下面是一维情况的检验思路：

% 使用解析解作为初值，在 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 接近 deg+1，说明空间离散正确；如果接近 deg，先查面积分符号；如果乱七八糟，直接查雅可比和通量方向。这套验证流程对我排查代码问题比任何调试器都管用。从那以后，我每次拿到别人的DG代码，第一件事不是看主循环，而是先跑一个已知解析解的收敛阶测试。这个习惯救了我很多次，也希望帮到你。

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