三对角方程组与追赶法：沿着相邻关系求解

把一根均匀细杆的左端维持在 ，右端维持在 。杆的侧面绝热，内部没有热源。等温度不再随时间变化时，在杆内等距选 4 个点，它们的温度分别是多少？

在这个理想模型里，每个内部点的温度等于左右邻点温度的平均。于是，虽然有 4 个未知温度，每条方程只牵涉相邻的两个或三个未知数。把点增加到一万个，这种局部关系也不改变。

先试着给这 4 个点分配温度，再逐点检查是否符合平均关系。可以按位置估计，也可以从一组相同的温度开始修改。需要找到的是让所有点同时平衡的分配；接下来再把这个任务整理成方程，并寻找不必反复试填的计算办法。

从相邻温度看三对角结构

把内部温度依次记为 ，下标表示从左到右的位置。两个端点温度已知，记为 、，它们不计入未知数。

“等于两侧平均”可以写成

这里假设材料和截面积均匀、节点等距，且没有侧向散热或内部热源。平衡时，相邻等长小段传递的热量相等，温度差也就相等，得到上面的关系。对于更一般的导热问题，系数和右端会改变，不能直接套用“取平均”。

在下面填一组温度，点击“检查当前试算”。读数是各点温度偏离左右邻点平均值的程度，单位仍为 ℃。可以先保留一次结果，再修改偏差较大的位置，比较整组温度是否更接近平衡；修改一个点也会影响邻点，因此只修好一处还不够。

如果暂时没有修改思路，可以逐级展开提示，或直接载入方程组的解；后面的正文也会完整算出结果，不需要先完成试算才能继续阅读。改变右端温度后，可以检验自己发现的分配方式还能否使用。

检查得到的是这组数是否满足局部规则的证据，并不能凭几次改进就保证某种试填策略总能成功。若观察到相邻温差相等，可以将这个猜想代回平均关系核对： 正好等价于上面的方程。对这根无热源的均匀细杆，4 个内部点把总温差分成 5 段，每段 20 ℃，于是得到 ℃。这给出一个可核对的解；当系数或热源改变时，还需要更通用的办法。

将已知边界移到等号右边，4 条方程为

第一行右端的 0 来自左端温度，最后一行右端的 100 来自右端温度。把系数排列成表，就是系数矩阵；将未知数和右端各排成一列，可以简写成 ：

从左上到右下的 构成主对角线，紧挨着它的上下两条斜线是副对角线。只有这三条对角线上允许出现非零元素，这样的矩阵叫三对角矩阵，对应的线性方程组叫三对角方程组。三条线上也可以有零，不要求每个位置都非零。

一般的输入长什么样

有 个未知数时，采用同一套记号：

是主对角系数， 是下对角系数， 是上对角系数， 是已知右端。约定 、，首尾方程省略不存在的未知数项。输入是这四列数，输出是满足所有方程的 ；求唯一解还需要矩阵可逆，并且本文不换行的算法需要后面列出的主元条件。

三对角矩阵不必对称， 和 可以不同。温度例子只是其中一种特别容易观察的输入。

先向前消掉，再向后算回来

先只看前两条方程：

不必同时猜 4 个温度。先尝试一个更局部的任务：把第一条方程乘以一个数再加到第二条上，能否让第二条里的 消失？可以先算出所需的倍数，再核对剩下的系数。这里只改变方程的表达方式，不要求改变未知温度。

第一行的 系数是 2，第二行是 −1，所以乘以 再相加， 恰好抵消。第二条变为

这叫消元：用一条已知方程消去另一条中的一个未知数。之所以不改变解，是因为新第二行等于“旧第二行加半倍第一行”，反过来再减去半倍第一行就能恢复。

接下来，用新第二行消掉第三行里的 ，再用新第三行消掉第四行里的 ：

处理到哪一行	用什么操作消去前一个未知数	处理后的该行
第 1 行	保留原式
第 2 行	加上第 1 行的 倍
第 3 行	加上新第 2 行的 倍
第 4 行	加上新第 3 行的 倍

到这里，每条方程中只剩当前未知数和它右边的一个未知数，最后一行更是只剩 。这就是“追”的方向。

下面是这套计算规则的分步演示。点击“下一步”，观察当前行左侧的系数变成 0，主对角和右端怎样更新。完成消元后继续单步，查看“赶”回来的回代过程；也可以切换到非对称例子，观察右端不为零时的更新。它用于核对每一步操作，与开头自行构造温度的任务相互补充。

从最后一行开始，计算顺序为

这叫回代：先算出末端未知数，再将它代入上一行。把 放回原方程，例如第三行是 ，满足要求；两端加上去，温度恰好每段升高 。

这种针对三对角结构的高斯消去法，称为追赶法，也叫 Thomas 算法。它是两趟有限步骤的直接解法，不是反复猜测直到收敛的迭代法。NPTEL 课程中的 Thomas 算法采用的是先消元、再回代的同一思路。

把手算变成递推公式

用 表示消元后第 行的主对角系数，用 表示它的新右端。原来的 与更新后的 分开记，避免把更新前后的数混在一起。

第一行不需要消元，因此

当第 行已经变成

第 行仍然含有 。只要 ，令

再做“第 行减去 倍的第 行”，便能消掉 。由此得到向前消元公式：

是为了抵消左邻系数所用的倍数。温度例子里 ，所以 ；“减去负的半倍”正好就是手算时“加上半倍”。上对角系数 不用更新。

消元结束后，第 行是 。因此向后回代公式为

有些教材会继续把每行除以 ，令 、，写成 。这与上面的公式等价，只是保存的中间量不同。

为什么消元不会越算越宽

一般矩阵消元时，原本为零的地方可能出现非零数，称为“填充”。这里，在不交换行的前提下，第 行消元后最多只含 ，拿它更新第 行，只会改变 和 的系数，影响不到更远的未知数。

每轮操作保持两件事：已处理的行都只含当前及右邻未知数；新方程组与原方程组有相同的解。向前循环处理完第 行后得到上双对角方程组。若所有 都非零，末行唯一确定 ，随后每一行依次唯一确定 。这给出了算法的正确性理由。

追赶法与LU分解的对应关系

把所有消元倍数 收集起来，可以将原矩阵写成 ：

的主对角线都是 1，只有紧邻的下方可能非零； 只有主对角线及紧邻上方可能非零。相乘后，下对角元素是 ，主对角元素是 ，上对角元素仍是 ，因此确实得到 。

求解 等价于先求 ，再求 。前一个方程给出 ，后一个就是回代。

若同一个 要配上很多组不同右端，可以保存 ，只做一次分解，再分别求解。每个新右端仍需线性工作量。

为什么计算量是线性的

令 为未知数个数，按固定精度计算，把一次加、减、乘、除各算作一次算术操作。

向前每处理一行，需要一次除法得到 ，两次乘法和两次减法更新 ，共 5 次；这样的行有 个。回代末行需要 1 次除法，其余每行各需要一次乘法、一次减法、一次除法，共 次。合起来是

这是本文这组公式的算术操作计数，未计入读写数组、检查输入等开销，也不等于运行毫秒数。它说明求一个右端的总时间为 ：规模翻倍，主要计算量约翻倍。若用普通稠密高斯消元而不利用零元素结构，经典算法的计算量为 ；差异来自是否利用稀疏结构，而非“追赶法与消元法”互不相关。

矩阵也无需保存成 的整张表。只保存三条对角线，去掉首尾占位的零，共需 个数；完整矩阵需要 个数。

拖动规模，比较保存同一个矩阵需要多少个数，以及本文公式的算术次数。这里不运行大矩阵计算，读数直接由存储长度和计数公式得出。

例如 ，三条对角线共 个数，稠密矩阵有 个数。仅按每数 8 字节估算，分别约为 MB 与 MB，均采用十进制 MB，且不含右端、输出和容器开销。

使用 等长度为 的数组，总空间仍为 。覆盖输入可以减少常数，但若要核对原方程的残差，就应保留或重建原系数。这里的线性空间是包含输出的量级，并不是“只需要常数空间”。

有唯一解，追赶法也可能走不下去

消元公式反复除以 ，回代还要除以 。这些作为除数的数叫主元。本文追赶法不交换方程的顺序，要求所有主元非零；矩阵可逆并不自动满足这个要求。

一个可手算的反例是

第一行直接给出 ，第二行给出 ，唯一解非常清楚。但追赶法的第一步主元是 ，接下来要算 ，无法执行。将两行交换后，第一主元变为 1，就能继续。

非零却很小的主元

把左上角的 0 改成一个小数 ，保持右端 不变：

当 时，精确解为

零主元的问题很明确；把它改成非零数以后，计算是否就足够可靠？先看默认的 ，再自行选择一个更小的非零数，或回到零主元对照。实验会保留上一次输入与误差，便于检查“能继续算”与“算得准确”是否总是一回事。每次比较时，两种方法使用相同的方程组。

面板中的解误差是计算结果与上面解析参考值的最大绝对差；残差是将计算结果代回原方程后，左右两侧差的最大绝对值。先观察两项读数，再比较两边的消元倍数；需要时可以展开提示。

若选到 或 ，本例不换行计算会出现明显的解误差。追溯中间量，不换行时 ，交换两行后则是 。前一种写法产生很大的中间量，随后相减可能丢失有效数字，再在回代中放大影响。这修正了“只要主元非零就足够可靠”的猜想；并不表示所有非零小主元都会产生同样大小的误差。

实验使用 JavaScript Number 的双精度浮点运算，显示的小数位数不参与算法。以解析式在双精度下计算的值作为参考；这里 ，解析式不需要除以很小的数，足以分辨示例中的大误差。接近双精度分辨率的差异，只能解释为浮点量级的差异。

这里的解误差是 ，其中 是计算值， 是解析参考值；残差是 ，衡量结果满足原方程的程度。两者不是同一个量。小残差一般不能单独保证小解误差，还要看问题对扰动的敏感程度。

本例在 时，矩阵并不接近奇异矩阵：其逆为

各元素仍有界。因此，这个例子主要揭示不合适的主元造成的算法问题。实际软件会根据矩阵性质选算法；例如 LAPACK 的 DGTSV为一般三对角系统使用带部分选主元的高斯消元。它仍利用带状结构，但换行可能让 多出一条上对角线，不能直接照搬本文不换行的存储和递推。

哪些条件可以保证主元非零

两个常用的充分条件是：

• 按行严格对角占优：每行都有 ，缺失的首尾系数按 0 计。主对角绝对值严格大于同一行其他系数绝对值之和。
• 实对称正定：矩阵对称，且对任意非零实向量 ，都有 。这个二次式可以理解为由矩阵定义的一种严格为正的“能量”。

在精确算术中，这两类矩阵都能保证本文消元所需的主元非零。它们是充分条件，不是必要条件；不满足时要继续检查，不能直接断言无解或算法必定失败。仅有非严格的“”也不够，例如 就是奇异矩阵。

两个充分条件的证明，以及温度矩阵为什么适用

严格对角占优。 首行给出 。若已有 ，则

因此每个 都非零，包括 时的末行。这是可逐行检查的归纳证明。

实对称正定。 每个左上角 子矩阵也正定，其行列式 为正。令 ，按三对角行列式展开，

将 代入消元递推，得到 ，所以没有零主元。这里使用了正定矩阵的标准性质：所有左上角主子式为正。

温度例子的内部行满足 ，并非每行都严格占优。但令 ，可以直接展开验证

只要内部向量不全为零，就不可能每个相邻差都为零，所以矩阵正定。还可以直接代入递推得到 ，温度例子确实能一路消元到底。

非零主元的数学保证，与浮点误差最终有多大，是不同层次的问题。整体尺度、舍入和问题条件性仍会影响计算；将任意小主元硬改成一个固定小数，会改变原方程，不能作为通用修复方式。

一个可以直接运行的实现

下面用 Python 实现正文的不换行版本。数学下标从 1 开始，代码下标从 0 开始；四个输入数组都长为 ，首尾占位为 a[0]=0、c[-1]=0。函数复制输入，保留原系数便于验算。

import math


def thomas(a, b, c, d):
    """求解三对角系统；要求全部消元主元非零，不执行选主元。"""
    a, b, c, d = [list(map(float, v)) for v in (a, b, c, d)]
    n = len(b)
    if n < 1 or any(len(v) != n for v in (a, c, d)):
        raise ValueError("四个数组必须等长且非空")
    if not all(math.isfinite(t) for v in (a, b, c, d) for t in v):
        raise ValueError("输入必须是有限数")
    if a[0] != 0 or c[-1] != 0:
        raise ValueError("要求 a[0]=0、c[-1]=0")

    u, r = b.copy(), d.copy()

    def check_pivot(i):
        if u[i] == 0 or not math.isfinite(u[i]):
            raise ArithmeticError(f"第 {i + 1} 个主元为零或非有限数")

    check_pivot(0)
    for i in range(1, n):
        m = a[i] / u[i - 1]
        u[i] -= m * c[i - 1]
        r[i] -= m * r[i - 1]
        check_pivot(i)
        if not math.isfinite(r[i]):
            raise ArithmeticError("右端计算上溢")

    x = [0.0] * n
    x[-1] = r[-1] / u[-1]
    for i in range(n - 2, -1, -1):
        x[i] = (r[i] - c[i] * x[i + 1]) / u[i]
    if not all(math.isfinite(t) for t in x):
        raise ArithmeticError("回代出现非有限数")
    return x


a = [0, -1, -1, -1]
b = [2, 2, 2, 2]
c = [-1, -1, -1, 0]
d = [0, 0, 0, 100]
x = thomas(a, b, c, d)
print(x)  # 约为 [20.0, 40.0, 60.0, 80.0]

# 用原始方程核对最大绝对残差。
residual = max(
    abs(b[i] * x[i]
        + (a[i] * x[i - 1] if i > 0 else 0)
        + (c[i] * x[i + 1] if i + 1 < len(x) else 0)
        - d[i])
    for i in range(len(x))
)
print("最大绝对残差:", residual)

这里检测零主元和非有限数，是为了明确报告无法继续的算术步骤；它没有诊断所有小主元风险，也不提供精度保证。对一般输入，应使用带选主元和误差诊断的成熟求解器，而不是仅凭“程序返回了一个数组”判断可靠。

从一维相邻关系到其他问题

回到细杆，改变两个条件，先决定原来的分配规律与求解方法各需怎样调整。判断时可以列出受影响的那一行，不必先算出整个新解。

加入局部热源。 两端仍为 0 ℃、100 ℃，仍有 4 个内部点，但某个已知热源使第二条离散方程改为 ，其余三条不变。这里的 10 是将热源和导热参数折算后的方程右端量，单位为 ℃。原来的等温差分配还满足方程吗？改变的是矩阵、右端，还是两者？需要重新选择求解算法吗？

核对：局部热源改变了什么

将原解代入第二行，左边仍为 0，已经不等于新的右端 10，因此等温差分配失效。但系数矩阵完全没变，仍是已经证明正定的三对角矩阵，可以继续使用追赶法。右端改为 ，所得解为 。逐行代入可得 ，可以独立核对。如果已经保存 LU 分解，只需为新右端重新前代、回代。

改成周期连接。 再考虑由至少 3 个节点组成的闭环：第一个节点同时与第二个和最后一个相邻。这是新的周期模型，不再保留细杆两端固定温度的条件。尝试写出第一行和最后一行，检查原来存三条对角线的数组是否还能容纳全部关系。

核对：周期连接为何需要另外处理

若仍使用相邻平均关系，第一行含 ，最后一行含 。系数矩阵右上角、左下角多出非零元素，形成循环三对角结构，不能直接交给本文的普通三对角实现。并且，没有固定温度或其他约束时，所有节点取任意同一常数都满足平均关系，解不唯一。实际求解还需要检查约束与可解性，不能只把两个角上的系数漏掉。

一维差分方程、某些三次样条的系数系统，也可能产生三对角方程组。判断能否采用本文的方法，要同时检查相邻关系形成的结构，以及主元条件。二维网格通常会产生更宽的带状或块结构，也不再是每行只有这三个标量系数。

对于本文的细杆模型，解析温度本来就是一条直线，等距节点上的平均关系恰好精确成立；这不意味着所有连续问题离散后都没有误差。追赶法负责解给定的离散方程组，离散方程是否准确描述物理问题，需要另外判断。

当相邻关系确实形成普通三对角结构，且主元条件合适时，每一行只向下一行传递少量信息，两趟扫描就能求解。这也是从矩阵的零元素分布中节省计算和存储的一个具体例子。

参考资料

1. NPTEL：Thomas Algorithm for Tridiagonal and Block Tridiagonal Matrices。用于核对算法机制；课程使用归一化形式，本文使用未归一化的 。
2. LAPACK：DGTSV 官方文档与源代码。用于核对一般三对角系统的部分选主元策略及第二条上对角线存储。
3. ScaLAPACK Users’ Guide：Linear Equations。对照一般带状、可不选主元和正定三对角系统的不同分解；其中并行实现还包含额外置换，不等同于本文的串行两趟扫描。