交替方向隐式法

交替方向隐式法（alternating direction implicit method，ADI）将多维隐式问题拆成若干沿单一坐标方向的线性求解。在矩形张量网格上，一个二维扩散时间步可以分成逐行和逐列的三对角求解，减少直接处理整个二维矩阵的成本。拆分的次序、系数和边界处理共同决定精度与稳定性。

一个二维扩散点的两个子步

考虑单位正方形内的 $u_{t} = κ (u_{x x} + u_{y y})$ ，边界恒为零。设 $κ = 1$ ，两方向网格间距均为 $1 / 2$ ，只有中央一个内部点 $U$ 。两个方向的离散算子分别为 $A_{x} U = - 8 U, A_{y} U = - 8 U$ 。

采用 Peaceman–Rachford 次序，步长 $τ = 0.05$ ，先求 $(I - τ A_{x} / 2) U^{*} = (I + τ A_{y} / 2) U^{n},$ 再求 $(I - τ A_{y} / 2) U^{n + 1} = (I + τ A_{x} / 2) U^{*} .$ 初值为一时，首式为 $1.2 U^{*} = 0.8$ ，所以 $U^{*} = 2 / 3$ ；次式给 $U^{n + 1} = 4 / 9$ 。半离散精确值为 $e^{- 16 τ} = e^{- 0.8} \approx 0.449328964$ 。中间变量是拆分算法的辅助状态，不必把它解释成真实半时刻的高精度解。

行列分解的代数结构

一般矩形网格中， $A_{x} = κ D_{x x}$ 只连接同一行的邻点， $A_{y} = κ D_{y y}$ 只连接同一列的邻点。因此第一子步中，每一行单独解三对角系统；第二子步中，每一列单独解三对角系统。若两个方向各有 $M_{x}, M_{y}$ 个未知点，每步总工作量为 $O (M_{x} M_{y})$ ，存储也为 $O (M_{x} M_{y})$ ，而不需要构造稠密的二维逆矩阵。

在常系数、张量积网格、相容的齐次边界条件下， $A_{x} A_{y} = A_{y} A_{x}$ 。消去中间变量得到 $(I - τ A_{x} / 2) (I - τ A_{y} / 2) U^{n + 1} = (I + τ A_{x} / 2) (I + τ A_{y} / 2) U^{n} .$ 展开两边，除了克兰克—尼科尔森格式的项，还出现 $τ^{2} A_{x} A_{y} (U^{n + 1} - U^{n}) / 4$ 。对光滑解，状态增量为 $O (τ)$ ，所以它贡献单步 $O (τ^{3})$ 的拆分差异。这解释了二阶时间精度，也说明两个格式并非同一个代数公式。Netlib 的矩阵迭代资料介绍了 ADI 的方向分解和因子化思想。^[1]

两个方向共同衰减的模式中，ADI 与精确指数和完整梯形推进比较 — 单内部点模型的两方向特征值均为 −8。横轴为时间步长，纵轴为一步放大因子；ADI 是两个方向梯形因子的乘积。

一个模式的稳定性证明

对同时对角化的两个对称负定离散算子，取公共特征向量，令 $A_{x} v = - λ_{x} v, A_{y} v = - λ_{y} v$ ， $λ_{x}, λ_{y} > 0$ 。该模式的放大因子为 $R = \frac{1 - τ λ_{x} / 2}{1 + τ λ_{x} / 2} \frac{1 - τ λ_{y} / 2}{1 + τ λ_{y} / 2} .$ 对任意 $b > 0$ ， $| 1 - b | < 1 + b$ ，所以每个因子的模小于一，乘积的模也小于一。公共特征向量构成正交基，模式能量相加便得到离散二范数不增。这个证明把无条件稳定性限定在上述可分离对称扩散问题，不能直接推广到任意混合导数、变系数或不相容的边界离散。

大步长仍会影响准确性与耗散。若两方向因子均接近负一，乘积反而接近正一，真实应迅速消失的高频模式可能长时间残留。ADI 的方向分解节省线性求解成本，并不取消时间分辨率的要求。

有解变式仍取 $λ_{x} = λ_{y} = 8$ ，把步长改为 $0.1$ ，两个因子均为 $3 / 7$ ，总因子 $9 / 49 \approx 0.183673469$ 。整体梯形法的因子却为 $(1 - 0.8) / (1 + 0.8) = 1 / 9$ 。两者都稳定，但数值答案不同，不能把一个方法的误差表套给另一个。

零边界热方程的行列实现

下面用二维列表按“行对应 y、列对应 x”保存内部值。函数返回新矩阵，不修改输入；所有边界均为零。内部三对角解算器只用于本问题的正定系数。代码展示一次方向分裂，参数固定时还可以缓存两个方向的分解。

Python 3

import math


def adi_heat_step(values, kappa, dx, dy, dt):
    U = [list(map(float, row)) for row in values]
    if not U or not U[0] or any(len(row) != len(U[0]) for row in U):
        raise ValueError("nonempty rectangular array required")
    kappa, dx, dy, dt = map(float, (kappa, dx, dy, dt))
    if not all(math.isfinite(v) for row in U for v in row):
        raise ValueError("finite states required")
    if not all(math.isfinite(v) and v > 0 for v in (kappa, dx, dy, dt)):
        raise ValueError("positive finite parameters required")
    rx, ry = kappa*dt/(2*dx*dx), kappa*dt/(2*dy*dy)
    if not all(math.isfinite(1+2*r) for r in (rx, ry)):
        raise OverflowError("nonfinite coefficients")

    def solve(rhs, r):
        b, d = [1+2*r]*len(rhs), rhs[:]
        for i in range(1, len(d)):
            factor = -r/b[i-1]
            b[i] += factor*r
            d[i] -= factor*d[i-1]
        x = [0.0]*len(d)
        x[-1] = d[-1]/b[-1]
        for i in range(len(d)-2, -1, -1):
            x[i] = (d[i]+r*x[i+1])/b[i]
        return x

    ny, nx = len(U), len(U[0])
    half = []
    for j in range(ny):
        rhs = [(1-2*ry)*U[j][i] + ry*(U[j-1][i] if j else 0)
               + ry*(U[j+1][i] if j+1 < ny else 0) for i in range(nx)]
        half.append(solve(rhs, rx))
    answer = [[0.0]*nx for _ in range(ny)]
    for i in range(nx):
        rhs = [(1-2*rx)*half[j][i] + rx*(half[j][i-1] if i else 0)
               + rx*(half[j][i+1] if i+1 < nx else 0) for j in range(ny)]
        column = solve(rhs, ry)
        for j in range(ny):
            answer[j][i] = column[j]
    if not all(math.isfinite(v) for row in answer for v in row):
        raise OverflowError("nonfinite split solution")
    return answer


print(adi_heat_step([[1]], 1, 0.5, 0.5, 0.05))  # [[0.4444444444444445]]

非零时间依赖边界需要与拆分次序相容的边界修正。非矩形区域、强各向异性、交叉导数以及非线性扩散也可能破坏这里的简单行列结构，应根据相应算子重新设计和验证分裂。

参考资料

↑ Richard Barrett et al. Templates for the Solution of Linear Systems，方向分裂、预条件和迭代求解。

[netlib-1] Richard Barrett et al. Templates for the Solution of Linear Systems，方向分裂、预条件和迭代求解。

[1]