ADI隐式交替法（P-R格式）求解二维热传导方程及MATLAB实现

简介ADI隐式交替法及其P-R差分格式的MATLAB实现面向数值计算、偏微分方程数值解领域的学习者与研究人员适合具备偏微分方程和MATLAB基础的读者用于课程设计、毕业设计或科研入门。该方法将二维抛物型方程的隐式求解拆分为两个方向的一维子问题结合罚函数与松弛迭代处理边界条件和非线性因素支持较大时间步长并保持数值稳定性。资源包内含1个m文件压缩包整体仅760B代码精炼便于逐行分析。已有494人学习/下载。通过运行该程序可直观理解交替方向法的矩阵构造、迭代流程及离散化过程并可作为模板改写到热传导、流体扩散等实际场景中。资源包以单个脚本呈现便于直接打开对照分析快速复现典型算例。对正在学习隐式差分格式、希望掌握ADI算法实现细节的读者是一份紧凑而实用的参考。1. ADI隐式交替法解决的是哪一类问题二维隐式差分在网格加密后很难跑动原因是系数矩阵带宽随网格数一起涨ADI隐式交替法Alternating Direction Implicit换了个思路把二维问题沿 x、y 两个方向拆成两个半步每个半步只解一个三对角系统稳定性接近全隐式成本却从大型稀疏求解降到线性量级。P-RPeaceman-Rachford就是这种交替方向法里最常用的一族时间推进格式常和“交替隐式”“隐式差分”一起出现在老代码里。如果你从某个 ADI.rar 里解压出来的 MATLAB 脚本没有注释这篇可以把格式重推一遍让你知道 rx、ry、边界条件填在哪结果拿什么来验。2. P-R交替方向法的数学推导与离散矩阵2.1 二维隐式差分为什么卡在矩阵带宽上一维热传导方程 ∂u/∂t α·∂²u/∂x² 用隐式差分比如 Crank-Nicolson时未知量是一根线上的点系数矩阵是三对角Thomas 算法每个时间步成本为 O(N)。换成二维问题∂u/∂t α(∂²u/∂x² ∂²u/∂y²)如果直接对 ∂²u/∂y² 也做隐式离散每个时间步会得到五对角系数矩阵带宽约等于 min(nx, ny)。虽然五对角带内的直接解法仍是多项式复杂度但稀疏 LU 分解的填充元会在带宽外继续扩展内存和时间都随网格变方明显恶化。最常见的处理是走算子分裂把二维算子按坐标方向拆开每个方向各隐式一遍这就是 ADI 隐式交替法要干的事。三种常见格式的取舍可以放在一起比较格式每个时间步成本稳定性限制备注显式 FTCSO(N)N 为网格点数Δt ≤ h²/(4α)实现最简单时间步被死死按住全隐式/CN 直接解稀疏 LU带宽随 1/h 增长无条件稳定填充元让成本无法随网格线性增长ADI-P-R每步两次三对角求解 O(N)无条件稳定常数系数热传导存储、速度、稳定性兼顾标准的 P-R 实现不会去构造整个二维系数矩阵而是循环解若干次一维三对角系统这一点从代码结构上很容易认出来。2.2 Peaceman-Rachford 格式半步拆分P-R 格式把每个时间步拆成两步。记 μ αΔt/2第一步只对 x 方向隐式y 方向沿用旧层数据(1 − μ·δx²) u^{n1/2} (1 μ·δy²) u^n第二步反过来对 y 方向隐式(1 − μ·δy²) u^{n1} (1 μ·δx²) u^{n1/2}把两个半步展开对内部点 (i,j) 得到两个三对角系统第一步x方向隐式 -rx·u(i-1,j)^{n1/2} (12rx)·u(i,j)^{n1/2} - rx·u(i1,j)^{n1/2} u(i,j)^n ry·( u(i,j-1)^n - 2·u(i,j)^n u(i,j1)^n ) 第二步y方向隐式 -ry·u(i,j-1)^{n1} (12ry)·u(i,j)^{n1} - ry·u(i,j1)^{n1} u(i,j)^{n1/2} rx·( u(i-1,j)^{n1/2} - 2·u(i,j)^{n1/2} u(i1,j)^{n1/2} )其中rx α·Δt / (2·Δx²) ry α·Δt / (2·Δy²)rx 和 ry 决定了整个编程模型第一步固定 j即固定一行沿 i 方向解三对角第二步固定 i固定一列沿 j 方向解三对角。注意 u^{n1/2} 只是中间层不是真实时刻的值没有物理意义老代码里一般叫ustar或ut。2.3 稳定性和精度的边界对常数系数的抛物线型热传导方程P-R 格式的放大因子为G [ (1 - 4·rx·sin²(ξ/2)) · (1 - 4·ry·sin²(η/2)) ] / [ (1 4·rx·sin²(ξ/2)) · (1 4·ry·sin²(η/2)) ]每个因子都是 (1−a)/(1a) 形式|G| ≤ 1 恒成立所以 P-R 格式无条件稳定。但这不表示 Δt 可以无限放大时间方向是二阶精度Δt 太大时截断误差会盖过空间精度算出来虽然“稳定”数值已经偏离真实解。一般先把 rx、ry 控制在 0.22 之间再按 rx 不变的原则缩放 Δt是比较稳妥的起点。对流项或非线性项出现时无条件稳定结论不再成立。对流占优问题需要配合迎风离散且算子分裂还会引入额外的分裂误差必须单独做数值实验验证这是手写 ADI 代码最常见的问题来源。2.4 每个半步解的是什么样的系统以第一步为例固定第 j 行未知数是该行 x 方向上的全部网格点内部点 i2…nx-1 的方程组成一个标准三对角系统左右边界值贡献到第一个和最后一个方程。所以代码里能不能直接用 Thomas 算法取决于提取出的子矩阵是不是标准三对角。如果实现者把整行连同边界点塞进一个矩阵再求逆说明走了弯路那样的代码既慢又难维护。3. MATLAB实现Thomas算法与ADI主循环3.1 先写好一个能复用的Thomas求解器Thomas 算法是三对角系统的 O(N) 直接法最稳妥的写法是“前代求临时数组再回代”。常数系数场景下可以先把分母数组算好循环里每次都省一次除法。function x thomas_const(a, b, c, d) % thomas_const: 常系数三对角系统求解器 % 输入: % a - 下次对角线常数 % b - 主对角线常数 % c - 超对角线常数 % d - 右端向量, 长度 n % 输出: % x - 解向量 (列向量) d d(:); % 强制转成列向量 n numel(d); den zeros(n, 1); den(1) b; for i 2:n den(i) b - a * c / den(i-1); end y zeros(n, 1); y(1) d(1) / den(1); for i 2:n y(i) (d(i) - a * y(i-1)) / den(i); end x zeros(n, 1); x(n) y(n); for i n-1:-1:1 x(i) y(i) - c * x(i1) / den(i); end end这个版本把den放在每次调用里重新计算小网格下没什么问题。第 5 章会给预分解优化方案。调用时a、b、c是标量d是长度为 n 的列向量函数内部用d d(:)强制列向量避免行向量传入导致后续赋值形状错乱。3.2 ADI主循环先扫x方向再扫y方向下面是一个可以直接运行的完整函数算单位正方形上的热扩散四边 Dirichlet 零边界初始条件 sin(πx)·sin(πy)。function u adi_pr_heat(nx, ny, alpha, dt, nt) % ADI-PR: 交替方向隐式求解二维热传导方程 % u adi_pr_heat(nx, ny, alpha, dt, nt) % 输入: % nx, ny - x/y方向网格数(含边界点) % alpha - 热扩散系数 % dt - 时间步长 % nt - 时间步数 % 输出: % u - 最终温度场, 尺寸 ny×nx Lx 1; Ly 1; dx Lx / (nx - 1); dy Ly / (ny - 1); x linspace(0, Lx, nx); y linspace(0, Ly, ny); rx alpha * dt / (2 * dx^2); ry alpha * dt / (2 * dy^2); [X, Y] meshgrid(x, y); u sin(pi * X) .* sin(pi * Y); % 初始条件, 边界自动为零 for n 1:nt uold u; ustar uold; % u^{n1/2}, 边界保持旧值 % 第一步: x方向隐式, 固定每一行 j, 沿列方向解三对角 for j 2:ny-1 rhs uold(j, 2:nx-1) ... ry * (uold(j-1, 2:nx-1) ... - 2 * uold(j, 2:nx-1) ... uold(j1, 2:nx-1)); ustar(j, 2:nx-1) thomas_const(-rx, 12*rx, -rx, rhs); end % 第二步: y方向隐式, 固定每一列 i, 沿行方向解三对角 for i 2:nx-1 rhs ustar(2:ny-1, i) ... rx * (ustar(2:ny-1, i-1) ... - 2 * ustar(2:ny-1, i) ... ustar(2:ny-1, i1)); u(2:ny-1, i) thomas_const(-ry, 12*ry, -ry, rhs); end end end两个半步的右端分别来自另一个方向的显式贡献第一步把 y 方向差分用旧层 u^n 计算第二步把 x 方向差分用中间层 u^{n1/2} 计算。注意ustar必须先完整复制uold再更新内部点否则第二步读取ustar(2:ny-1, i-1)时可能已经覆盖本时间步的数据典型的“原地更新顺序错误”结果会先从四角偏掉。运行示例u adi_pr_heat(41, 41, 1.0, 0.0025, 20); surf(linspace(0, 1, 41), linspace(0, 1, 41), u);3.3 参数在实际代码里怎么填rx、ry 不是直接输入的物理量而是由网格和时间步算出来的。常见做法是先定空间分辨率再反推时间步只想快速看趋势取 rx ry 2即 dt 4·dx²/α需要论文级精度固定 rx ry 1保证时间误差不压过空间误差边界条件一变三对角系统首末行就需要补边界修正项见 5.1。网格均匀时 rx ry两个方向的 Thomas 系数相同网格非均匀时必须分别构造 x、y 两个方向的三对角系数。建议先把均匀网格调通再推广到非均匀网格。4. 算例验证解析解对比与rx/ry参数选择4.1 用解析解直接做对比单位正方形上α 为常数时存在精确解u(x,y,t) exp(-2π²·α·t) · sin(πx) · sin(πy)验证脚本如下把数值解和解析解做 L∞ 误差对比x linspace(0, 1, 41); y x; [X, Y] meshgrid(x, y); u adi_pr_heat(41, 41, 1.0, 0.0025, 20); t 20 * 0.0025; u_exact exp(-2*pi^2*t) * sin(pi*X) .* sin(pi*Y); err max(abs(u(:) - u_exact(:))); fprintf(t %.4f, L_inf误差 %.3e\n, t, err);t 0.05 时衰减因子 exp(−2π²·0.05) ≈ 0.373中心峰值还有约 0.37数值结果应该在 1e-3 量级。中间层 u^{n1/2} 在这一时刻的取值不代表任何物理场不用拿它和解析解比。4.2 网格加密时的时间步缩放做收敛阶测试最常见的失误是固定 dt 只加密空间网格。P-R 是时间二阶格式dt 不动时时间误差会成为底噪收敛曲线会先出现平台。固定 rx 1 时时间步随 dx² 同步缩小推荐参数如下dxdtnt跑到 t0.05用途0.05005.0e-310粗网格快速冒烟0.02501.25e-340正常误差对比0.01253.12e-4160观察二阶收敛h 减半后L∞ 误差大约缩到原来的 1/4说明空间二阶、时间二阶同时成立。如果你跑出的误差只缩小不到 2 倍先查边界处理再查中间层是否被原地覆盖。4.3 显式步长上限与常见报错显式 FTCS 在同一套网格下的稳定上限约 Δt ≤ h²/(4α)。以 dx 0.025、α 1 为例显式上限约 1.56e-4而上面 ADI 的 dt 1.25e-3 已经放大了 8 倍仍然稳定。实际工程里放大 8100 倍都很正常前提是时间截断误差在可接受范围。报错时先看这些信号Index exceeds number of array dimensions八成是第一步循环里行、列索引写反。记住 u 的第一维是 y不是 x报行号错误比如 Line 9先检查调用thomas_const时 rhs 的长度和系数矩阵维度是否一致结果出现 NaN 且从四角蔓延初始化或边界上某一点没有被赋值NaN 会在连续时间步里扩散。4.4 时间方向要单独验证只缩 h 不缩 dt 的测试常被拿来演示“隐式方法无限制”但它只能证明稳定性证明不了精度。建议至少做两组对比一组 h、dt 按比例缩放观察总二阶一组固定 h、只将 dt 减半观察 L∞ 误差是否按约 4 倍下降。后者能快速暴露出时间离散上的错误比如漏掉了一个半步。5. 进阶Neumann边界、稀疏矩阵装配与批量扫描5.1 Neumann边界用虚拟点处理左边界 x 0 处给 ∂u/∂x g 时用虚拟点公式 u(0,j) u(2,j) − 2·dx·g 消去边界外点。代入第一步的第一个方程-rx·u(1,j) (12rx)·u(2,j) − rx·u(3,j) rhs(1)变成(1rx)·u(2,j) − rx·u(3,j) rhs(1) rx·dx·g也就是说主对角线第一个元素从 12rx 改成 1rx右端加上 rx·dx·g。改代码时这一步比想象中容易漏只改矩阵不动右端误差会在边界附近先出现。5.2 用spdiags预装配矩阵加速逐行调用 Thomas 直观但 MATLAB 里批量解三对角系统更省事。注意内存布局第一步是沿行方向列方向变量求解需要把右端矩阵转置后再用左除。Nx nx - 2; Ny ny - 2; A_x spdiags([-rx*ones(Nx,1), (12*rx)*ones(Nx,1), -rx*ones(Nx,1)], -1:1, Nx, Nx); A_y spdiags([-ry*ones(Ny,1), (12*ry)*ones(Ny,1), -ry*ones(Ny,1)], -1:1, Ny, Ny); % 第一步: 对每一行解列方向三对角, 需要转置 rhs_x uold(2:ny-1, 2:nx-1) ... ry*(uold(1:ny-2, 2:nx-1) - 2*uold(2:ny-1, 2:nx-1) uold(3:ny, 2:nx-1)); ustar(2:ny-1, 2:nx-1) (A_x \ rhs_x.).; % 第二步: 对每一列解行方向三对角, A_y 直接作用在列向量上 rhs_y ustar(2:ny-1, 2:nx-1) ... rx*(ustar(2:ny-1, 1:nx-2) - 2*ustar(2:ny-1, 2:nx-1) ustar(2:ny-1, 3:nx)); u(2:ny-1, 2:nx-1) A_y \ rhs_y;[L,U]lu(A)返回的 L 已经吸收了行交换后续可以直接用U\(L\b)批量回代。网格超过 200×200 后这个版本对比逐行 for 循环的优势会很明显。5.3 批量扫描rx并观察误差曲线扫描时间步对误差的影响时注意 dt 和 nt 的对应关系dts [5e-4, 1e-3, 2.5e-3, 5e-3]; errs zeros(size(dts)); for k 1:numel(dts) dtk dts(k); u adi_pr_heat(41, 41, 1.0, dtk, round(0.05/dtk)); u_exact exp(-2*pi^2*0.05) * sin(pi*X) .* sin(pi*Y); errs(k) max(abs(u(:) - u_exact(:))); end semilogy(dts, errs, -o); xlabel(\Delta t); ylabel(L_\infty error);四种步长下误差会呈现近似二阶斜率的下降rx 超过 10 后曲线向上翘说明时间截断误差已经压过空间误差。如果让 Codex 这类工具直接帮你生成 ADI 代码最容易翻车的三个点是循环索引方向写反、未复制旧场导致原地更新、rx 和 ry 抄混写进同一方向。即使 rx 大到解仍然“稳定”地算完也只能说明误差没有增长不代表误差不大。本文还有配套的精品资源点击获取