2. 三对角方程组追赶法与通用消元法性能对比

问题设定与两种算法的实现路径

本节针对同一批三对角方程组 Ax = d，其中 A 仅在主对角线、上次对角线、下次对角线非零，分别用追赶法（Thomas 算法）与通用高斯消元（带部分主元）求解，并在相同规模 n 下比较时间复杂度、内存占用与数值误差。三对角矩阵的典型形式为：

| b1  c1   0   ... |   | x1 |   | d1 |
| a2  b2  c2   ... |   | x2 |   | d2 |
|  0  a3  b3  c3 ...| * | x3 | = | d3 |
| ...               |   | .. |   | .. |

追赶法利用三对角结构，把消元与回代压缩为两次线性扫描：前向消元只更新 b、d，回代只依赖相邻分量。通用高斯消元则把矩阵视为稠密 n×n 数组，逐列选主元、消去下方元素，再回代。两者在数学上等价，但计算路径与存储需求差异显著。

可运行的对比实现

下面给出 Python 实现，两种算法共用同一组三对角数据，并统一用残差无穷范数 ‖Ax − d‖∞ 与相对误差衡量数值精度。为便于观察规模效应，用 n = 100、500、1000 三档测试。

import numpy as np
import time

def thomas(a, b, c, d):
    """a: 下次对角(长度n, a[0]不用), b: 主对角, c: 上次对角(长度n, c[-1]不用), d: 右端"""
    n = len(b)
    cp = np.zeros(n); dp = np.zeros(n)
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]
    for i in range(1, n):
        m = b[i] - a[i] * cp[i-1]
        cp[i] = c[i] / m if i < n-1 else 0.0
        dp[i] = (d[i] - a[i] * dp[i-1]) / m
    x = np.zeros(n)
    x[-1] = dp[-1]
    for i in range(n-2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i+1]
    return x

def gauss_tridiag(a, b, c, d):
    """把三对角矩阵展开为稠密矩阵后做带部分主元的高斯消元"""
    n = len(b)
    A = np.zeros((n, n))
    for i in range(n):
        A[i, i] = b[i]
        if i > 0: A[i, i-1] = a[i]
        if i < n-1: A[i, i+1] = c[i]
    rhs = d.copy()
    for k in range(n-1):
        p = k + np.argmax(np.abs(A[k:, k]))
        if p != k:
            A[[k, p]] = A[[p, k]]
            rhs[[k, p]] = rhs[[p, k]]
        for i in range(k+1, n):
            f = A[i, k] / A[k, k]
            A[i, k:] -= f * A[k, k:]
            rhs[i] -= f * rhs[k]
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (rhs[i] - A[i, i+1:] @ x[i+1:]) / A[i, i]
    return x

def build_tridiag(n):
    a = np.full(n, -1.0); b = np.full(n, 4.0); c = np.full(n, -1.0)
    a[0] = 0.0; c[-1] = 0.0
    d = np.ones(n)
    return a, b, c, d

for n in (100, 500, 1000):
    a, b, c, d = build_tridiag(n)
    t0 = time.perf_counter(); x1 = thomas(a, b, c, d); t1 = time.perf_counter()
    x2 = gauss_tridiag(a, b, c, d); t2 = time.perf_counter()
    A = np.diag(b) + np.diag(a[1:], -1) + np.diag(c[:-1], 1)
    r1 = np.linalg.norm(A @ x1 - d, np.inf)
    r2 = np.linalg.norm(A @ x2 - d, np.inf)
    print(f"n={n:4d}  Thomas={t1-t0:.4f}s  Gauss={t2-t1:.4f}s  "
          f"残差 Thomas={r1:.2e}  Gauss={r2:.2e}")

运行后应观察到：追赶法耗时随 n 近似线性增长，高斯消元耗时近似按 n³ 增长；两者残差都在机器精度量级（约 1e-12 以下），说明在本例这类对角占优的三对角矩阵上，追赶法无需选主元即可保持稳定。

复杂度与内存的定量对照

以 n = 1000 为例，追赶法约 8000 次浮点运算，高斯消元约 6.7×10⁸ 次，相差约 8 万倍；存储上追赶法约 4000 个双精度数（约 32 KB），高斯消元约 10⁶ 个（约 8 MB）。这正是结构化算法相对通用算法的核心收益。

数值误差与稳定性判读

追赶法的稳定性依赖系数矩阵的对角占优性。若 |b[i]| ≥ |a[i]| + |c[i]| 且至少一处严格大于，则前向消元中分母不会趋近于零，误差不会被放大。反之，若某步 b[i] − a[i]·cp[i-1] 接近零，追赶法会出现严重误差甚至除零。此时应改用带主元的高斯消元，或对三对角矩阵做局部主元交换（部分选主元的追赶法）。