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

简介：ADI隐式交替法及其P-R差分格式的MATLAB实现，面向数值计算、偏微分方程数值解领域的学习者与研究人员，适合具备偏微分方程和MATLAB基础的读者用于课程设计、毕业设计或科研入门。该方法将二维抛物型方程的隐式求解拆分为两个方向的一维子问题，结合罚函数与松弛迭代处理边界条件和非线性因素，支持较大时间步长并保持数值稳定性。资源包内含1个m文件，压缩包整体仅760B，代码精炼，便于逐行分析。已有494人学习/下载。通过运行该程序，可直观理解交替方向法的矩阵构造、迭代流程及离散化过程，并可作为模板改写到热传导、流体扩散等实际场景中。资源包以单个脚本呈现，便于直接打开对照分析，快速复现典型算例。对正在学习隐式差分格式、希望掌握ADI算法实现细节的读者，是一份紧凑而实用的参考。

1. ADI隐式交替法解决的是哪一类问题

二维隐式差分在网格加密后很难跑动，原因是系数矩阵带宽随网格数一起涨；ADI隐式交替法（Alternating Direction Implicit）换了个思路，把二维问题沿 x、y 两个方向拆成两个半步，每个半步只解一个三对角系统，稳定性接近全隐式，成本却从大型稀疏求解降到线性量级。P-R（Peaceman-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 隐式交替法要干的事。

三种常见格式的取舍可以放在一起比较：

标准的 P-R 实现不会去构造整个二维系数矩阵，而是循环解若干次一维三对角系统，这一点从代码结构上很容易认出来。

2.2 Peaceman-Rachford 格式：半步拆分

P-R 格式把每个时间步拆成两步。记 μ = αΔt/2，第一步只对 x 方向隐式，y 方向沿用旧层数据：

(1 − μ·δx²) u^{n+1/2} = (1 + μ·δy²) u^n

第二步反过来，对 y 方向隐式：

(1 − μ·δy²) u^{n+1} = (1 + μ·δx²) u^{n+1/2}

把两个半步展开，对内部点 (i,j) 得到两个三对角系统：

第一步（x方向隐式）： -rx·u(i-1,j)^{n+1/2} + (1+2rx)·u(i,j)^{n+1/2} - rx·u(i+1,j)^{n+1/2} = u(i,j)^n + ry·( u(i,j-1)^n - 2·u(i,j)^n + u(i,j+1)^n ) 第二步（y方向隐式）： -ry·u(i,j-1)^{n+1} + (1+2ry)·u(i,j)^{n+1} - ry·u(i,j+1)^{n+1} = u(i,j)^{n+1/2} + rx·( u(i-1,j)^{n+1/2} - 2·u(i,j)^{n+1/2} + u(i+1,j)^{n+1/2} )

其中：

rx = α·Δt / (2·Δx²) ry = α·Δt / (2·Δy²)

rx 和 ry 决定了整个编程模型：第一步固定 j（即固定一行），沿 i 方向解三对角；第二步固定 i（固定一列），沿 j 方向解三对角。注意 u^{n+1/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)/(1+a) 形式，|G| ≤ 1 恒成立，所以 P-R 格式无条件稳定。但这不表示 Δt 可以无限放大：时间方向是二阶精度，Δt 太大时截断误差会盖过空间精度，算出来虽然“稳定”，数值已经偏离真实解。一般先把 rx、ry 控制在 0.2～2 之间，再按 rx 不变的原则缩放 Δt，是比较稳妥的起点。

对流项或非线性项出现时，无条件稳定结论不再成立。对流占优问题需要配合迎风离散，且算子分裂还会引入额外的分裂误差，必须单独做数值实验验证，这是手写 ADI 代码最常见的问题来源。

2.4 每个半步解的是什么样的系统

以第一步为例，固定第 j 行，未知数是该行 x 方向上的全部网格点，内部点 i=2…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(i+1) / 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^{n+1/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(j+1, 2:nx-1)); ustar(j, 2:nx-1) = thomas_const(-rx, 1+2*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, i+1)); u(2:ny-1, i) = thomas_const(-ry, 1+2*ry, -ry, rhs); end end end

两个半步的右端分别来自另一个方向的显式贡献：第一步把 y 方向差分用旧层 u^n 计算，第二步把 x 方向差分用中间层 u^{n+1/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^{n+1/2} 在这一时刻的取值不代表任何物理场，不用拿它和解析解比。

4.2 网格加密时的时间步缩放

做收敛阶测试最常见的失误是固定 dt 只加密空间网格。P-R 是时间二阶格式，dt 不动时时间误差会成为底噪，收敛曲线会先出现平台。固定 rx = 1 时，时间步随 dx² 同步缩小，推荐参数如下：

h 减半后，L∞ 误差大约缩到原来的 1/4，说明空间二阶、时间二阶同时成立。如果你跑出的误差只缩小不到 2 倍，先查边界处理，再查中间层是否被原地覆盖。

4.3 显式步长上限与常见报错

显式 FTCS 在同一套网格下的稳定上限约 Δt ≤ h²/(4α)。以 dx = 0.025、α = 1 为例，显式上限约 1.56e-4，而上面 ADI 的 dt = 1.25e-3 已经放大了 8 倍，仍然稳定。实际工程里放大 8～100 倍都很正常，前提是时间截断误差在可接受范围。

报错时先看这些信号：

• 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) + (1+2rx)·u(2,j) − rx·u(3,j) = rhs(1)

变成：

(1+rx)·u(2,j) − rx·u(3,j) = rhs(1) + rx·dx·g

也就是说主对角线第一个元素从 1+2rx 改成 1+rx，右端加上 rx·dx·g。改代码时这一步比想象中容易漏：只改矩阵不动右端，误差会在边界附近先出现。

5.2 用spdiags预装配矩阵加速

逐行调用 Thomas 直观，但 MATLAB 里批量解三对角系统更省事。注意内存布局：第一步是沿行方向（列方向变量）求解，需要把右端矩阵转置后再用左除。

Nx = nx - 2; Ny = ny - 2; A_x = spdiags([-rx*ones(Nx,1), (1+2*rx)*ones(Nx,1), -rx*ones(Nx,1)], -1:1, Nx, Nx); A_y = spdiags([-ry*ones(Ny,1), (1+2*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 大到解仍然“稳定”地算完，也只能说明误差没有增长，不代表误差不大。

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