克兰克—尼科尔森格式

克兰克—尼科尔森格式（Crank–Nicolson scheme）将空间离散算子在新旧两个时刻的作用取平均，是热方程等演化问题常用的二阶隐式时间格式。它对线性扩散问题具有无条件能量稳定性，但高频误差在大步长下可能换号振荡；无条件稳定并不等于无条件非负或任意步长都准确。

只有一个内部网格点的热传导

考虑 $0 < x < 1$ 上的热方程 $u_{t} = κ u_{x x}$ ，边界温度为零， $κ > 0$ 。令网格间距 $Δ x = 1 / 2$ ，只有中点一个未知温度 $U (t)$ 。二阶中心空间差分给出 $U^{'} = κ \frac{0 - 2 U + 0}{(1 / 2)^{2}} = - 8 κ U .$ 取 $κ = 1, U (0) = 1, Δ t = 0.1$ ，新旧时刻平均导数得到 $\frac{U^{n + 1} - U^{n}}{0.1} = - 4 (U^{n + 1} + U^{n}), U^{n + 1} = \frac{3}{7} U^{n} .$ 所以前两步为 $U^{1} = 3 / 7, U^{2} = 9 / 49$ 。用来比较时间误差的精确对象是空间离散系统的解 $e^{- 8 t}$ ；它仍与连续热方程有空间离散误差，不能混作一个误差来源。

一般离散与二阶精度

在均匀网格的内部节点定义 $(Δ_{h} U)_{i} = (U_{i - 1} - 2 U_{i} + U_{i + 1}) / Δ x^{2}$ 。格式为 $\frac{U^{n + 1} - U^{n}}{Δ t} = \frac{κ}{2} Δ_{h} (U^{n + 1} + U^{n}) .$ 令 $r = κ Δ t / (2 Δ x^{2})$ ，逐行写成 $\begin{matrix} - r U_{i - 1}^{n + 1} + (1 + 2 r) U_{i}^{n + 1} - r U_{i + 1}^{n + 1} & = r U_{i - 1}^{n} \\ + (1 - 2 r) U_{i}^{n} + r U_{i + 1}^{n} . \end{matrix}$ 零边界值直接省去；非零或随时间变化的边界还应加入两时刻的相应项。每一步需解一个三对角线性系统。

对光滑半离散真解 $Y^{'} = A Y$ ，区间积分使用梯形求积： $\begin{matrix} Y (t + Δ t) - Y (t) & = Δ t [A Y (t) + A Y (t + Δ t)] / 2 \\ + O (Δ t^{3}) \end{matrix}$ 。这说明单步状态缺陷为三阶，稳定后固定时间的全局时间误差为二阶；中心空间差分在四阶空间导数有界等条件下给出二阶空间误差。非光滑初值、初边值不相容或缺少所需正则性时，不能直接沿用同一误差常数和阶数。

该时间离散也是亚当斯—莫尔顿法的梯形公式应用到空间离散方程。数值 PDE 教材常以热方程比较它与显式、后向 Euler 方法。^[1]

无条件能量稳定性的证明

取离散内积 $⟨ U, V ⟩_{h} = Δ x \sum_{i} U_{i} V_{i}$ ，并将边界零值补齐。设中点状态 $W = (U^{n + 1} + U^{n}) / 2$ 。把格式与 $W$ 作内积，左边利用差平方得到 $\frac{‖ U^{n + 1} ‖_{h}^{2} - ‖ U^{n} ‖_{h}^{2}}{2 Δ t} = κ ⟨ Δ_{h} W, W ⟩_{h} .$ 将二阶差分展开，移位相邻求和并用端点为零，得离散分部求和恒等式 $⟨ Δ_{h} W, W ⟩_{h} = - Δ x \sum_{i = 0}^{M} {(\frac{W_{i + 1} - W_{i}}{Δ x})}^{2} .$ 右边非正，故任意正时间步长下 $‖ U^{n + 1} ‖_{h} \leq ‖ U^{n} ‖_{h}$ 。证明针对这个对称扩散算子和边界条件；非线性、反应项、非对称离散或近似线性求解器需要另作分析。

高频振荡与步长选择

若离散特征模式满足 $Y^{'} = - λ Y$ 、 $λ > 0$ ，其放大因子为 $R = \frac{1 - λ Δ t / 2}{1 + λ Δ t / 2} .$ 模始终小于一，但当 $λ Δ t > 2$ 时为负； $λ Δ t \to \infty$ 时趋于 $- 1$ ，很快的真实衰减模式在数值上只缓慢交替。

小步长梯形推进平滑衰减，大步长结果正负交替但绝对值下降 — 单内部节点的 λ=8。小步长 Δt=0.1 的因子为 3/7；大步长 Δt=1 的因子为 −3/5。能量下降仍允许温度出现负值。

上文模型若改用 $Δ t = 1$ ，首步就是 $- 3 / 5$ ，而真实解为正。这是完整的有解变式，也给出误用“稳定”的反例。对初始不连续或尖峰数据，可先用几个较小的后向 Euler 子步消减高频，再使用二阶梯形推进；具体起步与误差容限仍须核查。

三对角实现

下面对零 Dirichlet 边界求解一次时间步。输入内部节点值、扩散率、空间步长和时间步长，返回新列表，不修改输入。Thomas 消元在本问题的严格对角占优正定系统中没有零主元。每步时间与空间均为 $O (M)$ ；固定参数时可预先分解，避免重复消元。

Python 3

import math


def crank_nicolson_heat(interior, kappa, dx, dt):
    u = list(map(float, interior))
    kappa, dx, dt = map(float, (kappa, dx, dt))
    if not u or not all(math.isfinite(v) for v in u + [kappa, dx, dt]):
        raise ValueError("nonempty finite data required")
    if kappa <= 0 or dx <= 0 or dt <= 0:
        raise ValueError("positive parameters required")
    r = kappa * dt / (2 * dx * dx)
    if not math.isfinite(r) or not math.isfinite(1+2*r):
        raise OverflowError("coefficient outside floating-point range")
    n = len(u)
    old = [0.0] + u + [0.0]
    rhs = [r*old[i] + (1-2*r)*old[i+1] + r*old[i+2] for i in range(n)]
    diag = [1+2*r] * n
    for i in range(1, n):
        multiplier = -r / diag[i-1]
        diag[i] += multiplier*r
        rhs[i] -= multiplier*rhs[i-1]
    result = [0.0] * n
    result[-1] = rhs[-1] / diag[-1]
    for i in range(n-2, -1, -1):
        result[i] = (rhs[i] + r*result[i+1]) / diag[i]
    if not all(math.isfinite(v) for v in result):
        raise OverflowError("nonfinite linear solve")
    return result


u = crank_nicolson_heat([1], 1, 0.5, 0.1)
print(u)  # [0.4285714285714286]
print(crank_nicolson_heat(u, 1, 0.5, 0.1))  # [0.1836734693877551]
print(crank_nicolson_heat([1], 1, 0.5, 1))    # [-0.6]

参考资料

↑ Tobin A. Driscoll, Richard J. Braun. Fundamentals of Numerical Computation: Stiffness，扩散的空间离散与时间积分。

[fnc-1] Tobin A. Driscoll, Richard J. Braun. Fundamentals of Numerical Computation: Stiffness，扩散的空间离散与时间积分。

[1]