简介以浅水波传播、波浪模拟为典型场景一套面向海洋/水利工程领域学生与科研人员的ADI浅水波方程数值求解MATLAB源码围绕二维浅水波模型的交替方向隐式求解展开。代码清晰划分物理域参数、初始条件、时间步长与网格尺寸设定、ADI迭代以及波浪结果可视化等模块主程序qianshuibo.m可直接运行zhuigan.m则以追赶法辅助完成隐式格式下三对角线性方程组的求解两者配合便于快速验证算法流程。压缩包共2个m文件大小仅2KB整体结构精简适合入门阅读与二次开发。目前已有162人浏览学习可作为偏微分方程数值解、计算流体力学等课程的实践参考。通过逐段调试读者能掌握ADI交替迭代与追赶法的实现细节并迁移至其他抛物型方程或更复杂的浅水流动问题中。1. 二维浅水波方程里ADI 隐式求解到底赢在哪前阵子处理近岸波浪传播的数值实验老板给的算例从一维推到二维显式格式在网格加密后追不动。翻到 code_qianshuibo_ADI_浅水波方程_源码 这个 MATLAB 包里面只有 qianshuibo.m 和 zhuigan.m 两个脚本前者搭主流程后者实现追赶法却把二维浅水波方程的隐式求解跑得比想象中稳。ADI 的核心是交替方向隐式把大的二维问题拆成沿 x、y 方向的两步三对角求解避免直接处理五对角带状矩阵特别适合海洋工程、水利工程场景里需要长时间推进的波浪模拟。这篇顺着源码把离散格式、追赶法、边界条件和验证流程讲透新手能照着跑熟手可以直接改出自己要的网格策略。2. ADI 交替方向隐式与浅水波方程的二维离散2.1 从浅水波方程组到线性波动方程不考虑底摩擦、科氏力和对流非线性项时二维浅水波方程的连续方程与动量方程可以写成∂η/∂t H(∂u/∂x ∂v/∂y) 0∂u/∂t g ∂η/∂x 0∂v/∂t g ∂η/∂y 0其中 η 是自由面相对静水面的扰动H 是静水水深u、v 是垂向平均流速在 x、y 方向的分量g 是重力加速度。对第一个方程再取时间偏导并把后两个方程的 ∂u/∂t、∂v/∂t 代进去就能消去流速项得到一个只含 η 的二维标量波动方程∂²η/∂t² gH (∂²η/∂x² ∂²η/∂y²)这个消元步骤在很多教材里是一笔带过但实际写代码前最好自己推一遍。因为后面 ADI 的每一步隐式求解本质上都是在交替处理这个方程右端的 x 方向与 y 方向二阶导数项。如果直接从原始方程组做耦合离散每个未知点会同时牵连水位和两个流速分量未知量翻三倍ADI 的优势就体现不出来了。2.2 Peaceman-Rachford 分裂与两步半步格式对上面的波动方程如果直接做中心差分离散空间上会形成一个五对角结构虽然矩阵稀疏但工程上要自己维护行索引写起来麻烦。ADI 的做法是把一个完整时间步 [t^n, t^{n1}] 分成两个半步中间引入一个中间层 η^{n1/2}。第一步沿 x 方向隐式y 方向用已知的 n 层值显式处理(1 − r_x/2 · δx²) η^{n1/2} (1 r_y/2 · δy²) η^n第二步反过来y 方向隐式x 方向用刚算出的 n1/2 层值(1 − r_y/2 · δy²) η^{n1} (1 r_x/2 · δx²) η^{n1/2}这里的 δx²、δy² 是二阶中心差分算子r_x gH·dt²/dx²r_y gH·dt²/dy²。展开第一步后对内部网格点 (i, j) 具体写成−r_x/2 · η_{i-1,j}^{n1/2} (1 r_x) · η_{i,j}^{n1/2} − r_x/2 · η_{i1,j}^{n1/2} r_y/2 · η_{i,j-1}^{n} (1 − r_y) · η_{i,j}^{n} r_y/2 · η_{i,j1}^{n}可见每一行只涉及三个点形成三对角系统。第二步同理。这种分裂格式的精度是二阶的两个半步合并后与直接 Crank-Nicolson 格式的效果相当但每次只需解三对角矩阵计算成本和存储开销都小一个量级。2.3 三对角系统与无条件稳定性的工程含义对于线性化的浅水波方程上面的 PR 分裂格式是无条件稳定的也就是说时间步长不再受显式 CFL 条件的硬性约束。这一点在实际算例中非常有价值万岁·说人话就是同样的水槽算例显式格式在 dx 从 1 米加密到 0.25 米时dt 要跟着缩小四倍才能保证数值稳定而 ADI 可以把 dt 维持在相对较大的水平因为隐式格式内在耗散会吸收高频分量。但这不意味着可以无脑放大 dt。表格里给出的参数就说明了这个问题符号表达式工程含义r_xgH·dt²/dx²x 方向网格比控制离散格式的相位误差r_ygH·dt²/dy²y 方向网格比过大会导致波形失真dx, dy空间步长决定可解析的最小波长一般至少 8 个网格覆盖一个波长dt时间步长决定时间分辨率影响数值色散r_x 和 r_y 偏大时虽然计算不崩但波速会偏离物理值波峰位置提前或滞后波形出现非物理振荡。所以实际项目中我一般把 r_x、r_y 控制在 0.1 到 1.0 之间这是精度和计算量之间比较稳妥的折中区间。理解了这些系数的来龙去脉再看源码里 dt 的赋值逻辑就不觉得突兀了。3. qianshuibo.m 与 zhuigan.m 的源码级拆解3.1 主程序 qianshuibo.m 的物理域与初始条件这个包里 qianshuibo.m 是主程序代码开头通常是物理域定义和初始水位扰动设置。以典型的水槽算例为模板核心逻辑可以还原成下面这样% qianshuibo.m 主程序片段 Lx 100; Ly 50; nx 201; ny 101; dx Lx / (nx - 1); dy Ly / (ny - 1); H 5.0; % 静水水深 g 9.81; x (0:nx-1) * dx; y (0:ny-1) * dy; [X, Y] meshgrid(x, y); A 0.05; % 初始波幅 sigma 3.0; % 高斯波包宽度 eta A * exp(-((X-40).^2 (Y-25).^2) / (2*sigma^2));这里水深 H 和波幅 A 的量纲要保持一致A 取 0.05 米时对 5 米水深是小扰动符合第 2 章线性化的适用前提。如果 A 与 H 的比值超过 0.1线性化误差就会明显增大此时解出来的波形只能看趋势不能拿去做定量分析。网格数 201×101 对应 dx、dy 都是 0.5 米对 100 米×50 米的计算域来说网格分辨率适中。3.2 追赶法 zhuigan.m 的实现与调用约定三对角系统用追赶法求解是最自然的选择zhuigan.m 就是这个包里承担该任务的函数。追赶法本质是高斯消元法在三对角矩阵上的特化版本时间复杂度 O(N)不需要显式构造完整矩阵。一个标准实现如下function x zhuigan(a, b, c, d) % 追赶法求解三对角系统 % a: 下对角线, b: 主对角线, c: 上对角线, d: 右端项 n length(d); x zeros(n, 1); % 消元过程 for i 2:n w a(i) / b(i-1); b(i) b(i) - w * c(i-1); d(i) d(i) - w * d(i-1); end % 回代过程 x(n) d(n) / b(n); for i n-1:-1:1 x(i) (d(i) - c(i) * x(i1)) / b(i); end end注意 a、b、c、d 四个数组的长度与主对角线一致a(1) 在循环中不会被使用传入时置 0 即可。主程序调用时要按 ADI 第一步、第二步分别组装系数数组。第一步 x 方向的系数是 a −r_x/2、b 1 r_x、c −r_x/2右端项来自 n 时刻 y 方向的显式差分离散第二步 y 方向则把 r_x 换成 r_y右端项来自 n1/2 时刻 x 方向的差分结果。3.3 ADI 迭代循环与可视化输出有了追赶法函数主循环就可以按方向交替推进。为了减少内存开销我习惯预分配 η_half 和 η_new 两个中间变量避免循环内动态扩容% 时间步长 C 0.6; dt C * dx / sqrt(g * H); nt round(60 / dt); eta_old eta; eta_half eta; eta_new eta; % 构建 x 方向三对角系数 ax -rx/2 * ones(nx-2, 1); bx (1 rx) * ones(nx-2, 1); cx -rx/2 * ones(nx-2, 1); % 第一步x 方向隐式y 方向显式 for j 2:ny-1 rhs zeros(nx, 1); for i 2:nx-1 rhs(i) ry/2 * eta_old(j, i-1) (1 - ry) * eta_old(j, i) ry/2 * eta_old(j, i1); end eta_half(j, 2:nx-1) zhuigan(ax, bx, cx, rhs(2:nx-1)); end eta_half(:, 1) eta_half(:, 2); eta_half(:, nx) eta_half(:, nx-1);第二步的代码与第一步对称只是换成了沿 i 循环、组装 y 方向的系数数组求解后得到 η_new再更新 eta_old eta_new。这里最容易被忽视的是数组维度的对应关系如果初始赋值用的是 meshgrid(x, y)则 eta 的第一维对应 y第二维对应 x赋值时务必保持下标一致。可视化用 surf 或 pcolor 直接绘制即可pcolor 前记得转置figure; pcolor(x, y, eta_new); shading interp; axis equal; colorbar; caxis([-A A]); title(sprintf(t %.2f s, n * dt)); drawnow;pcolor 默认按矩阵行列为坐标不转置会看到波面沿对角线偏转 90 度。这个坑在二维波模拟里非常常见检查的时候先看一眼 x、y 尺寸和矩阵维度是否匹配别急着怀疑算法。4. 网格、时间步长与边界条件的实战调参4.1 稳定性约束与时间步长选取ADI 虽然线性无条件稳定但时间步长仍受精度约束。下面这组对照来自同一个高斯波包算例在相同的初始条件和物理域下改变 dx 与 dt 会得到差别明显的波面算例nx×nydx [m]dy [m]dt [s]r_xr_y波形特征粗网格201×1010.500.500.050.490.49波形平滑相位略提前加密网格401×2010.250.250.020.310.31与解析解吻合较好大步长201×1010.500.500.207.857.85仍稳定但振荡失真明显r_x 从 0.49 涨到 7.85 时计算仍能推进但在波峰位置会出现明显的过冲波前出现细微的锯齿状。这说明无条件稳定指的是误差不会指数爆炸不代表波形保真。如果算例本身有非线性项或者变水深Δt 还得进一步收紧一般取显式 CFL 限制的 3 到 5 倍以内比较安全。4.2 固壁反射边界与海绵吸收层浅水波模拟最常见的两种边界是固壁和开边界。固壁边界要求法向速度为零反映到水位场上就是 Neumann 零梯度条件即边界内外的 η 相等。在 ADI 每半步完成之后直接赋边值即可% 固壁反射边界 eta_half(:, 1) eta_half(:, 2); eta_half(:, nx) eta_half(:, nx-1); eta_half(1, :) eta_half(2, :); eta_half(ny, :) eta_half(ny-1, :);开边界如果要吸收往外传播的波常见做法是在边界附近加一个海绵层Sponge Layer让波能在几行网格内逐渐衰减掉。实现时给靠近边界的网格乘上一个衰减系数数值上等价于增加瑞利阻尼项% 海绵层边界外 20 个网格内逐步衰减 nsponge 20; beta 0.15; sponge exp(-beta * (nsponge-1:-1:0)); eta_half(1:nsponge, :) eta_half(1:nsponge, :) .* sponge; eta_half(end-nsponge1:end, :) eta_half(end-nsponge1:end, :) .* sponge(end:-1:1);这里的 sponge 数组沿 y 方向作用x 方向边界同理。衰减系数 beta 取 0.1 到 0.3 比较合适太大会引入新的反射太小则吸收不干净。调试时可以放一个初始波包观察波碰到边界后的反射振幅反射波小于入射波振幅的 1% 就可以认为海绵层参数合格了。4.3 性能瓶颈定位与向量化改进MATLAB 里跑 ADI 最怕的就是用纯 for 循环嵌套双层网格尤其第一步内层循环逐一访问 eta_old(j, i) 时元素是按列访问的内存跳步严重204 万的网格点规模下效率明显下降。常见优化是把内层循环改成向量切片操作% 向量化实现第一步的右端项 rhs(2:nx-1) ry/2 * eta_old(j, 1:nx-2) (1-ry) * eta_old(j, 2:nx-1) ry/2 * eta_old(j, 3:nx);这样每个时间步少两层循环计算量不变但常数时间大幅下降。再进一步可以把 x 方向所有 j 行的三对角系统合并成一次矩阵运算但代码可读性会明显下降我一般只在确认算法正确后做这一步。调性能时先用 profile 定位最耗时的行再决定是否向量化不要一开始就写满矩阵版。5. 收敛性检验与限制器改造让 ADI 结果更可信线性浅水波方程存在解析解这让收敛性验证变得非常直接。取二维谐波解 η A·cos(k_x·x)·cos(k_y·y)·cos(ω·t)其中 ω² gH(k_x² k_y²)。在同样的网格上用 ADI 推进若干个周期每个时间步结束后计算一次 L2 相对误差kx 2*pi/40; ky 2*pi/60; omega sqrt(g*H*(kx^2 ky^2)); eta_exact A * cos(kx*x) .* cos(ky*y) .* cos(omega*t_now); errL2 sqrt(sum(sum((eta_new - eta_exact).^2)) / sum(sum(eta_exact.^2)));把 dx、dy、dt 同时减半重新跑一遍L2 误差应大约缩小到原来的四分之一即二阶收敛。如果斜率掉到 1.x通常不是 ADI 的问题而是边界处理不准确或初始条件本身不光滑。检查固壁边界的赋值顺序确保每半步结束后边界都被刷新别只在主循环末尾更新一次。非线性修正方面当波幅偏大或地形有突变时η 可能出现过冲。简单有效的做法是给每步更新后的 η 做一次斜率限制max_slope 0.1 * dx; eta_new(2:end-1, :) ... min(max(eta_new(2:end-1, :), eta_new(1:end-2, :) - max_slope), ... eta_new(1:end-2, :) max_slope);这相当于一个粗糙的 TVD 限制器能抑制局部振荡但也会稍微抹平波峰。要保留更多细节可以使用 minmod 限制器比较相邻网格的三个斜率取绝对值最小者参与差分更新。最后把 dt 缩小到理论极限的 60% 再跑一遍对比两次结果的波峰位置偏差小于半个网格宽度时这个结果就可以放心交给下游分析了。本文还有配套的精品资源点击获取