追赶法(三对角矩阵算法)模拟器 返回
数值分析

追赶法(三对角矩阵算法)模拟器

体验以 O(n) 求解三对角线性系统的追赶法(TDMA)。它将一维泊松边值问题 −u''=f 离散化并求解,让你实时观察从前向消元到回代的过程,以及与解析解的误差。

参数设置
内部网格点数 n
三对角系统的未知数个数(分辨率)
左端边界值 u(0)
狄利克雷边界条件(x=0 处解的值)
右端边界值 u(1)
狄利克雷边界条件(x=1 处解的值)
源项大小 f
方程 −u''=f 右端的强度
源分布
右端 f(x) 的空间形状
计算结果
网格点数 n
运算次数(≈8n)
中央 u 值
最大 u 值
最大误差(对比解析解・均匀分布时)
复杂度等级
算法可视化 — 前向消元与回代

在绘制解 u(x) 曲线的同时,高亮标记从左→右(前向消元)、再从右→左(回代)扫过。右上角小矩阵仅点亮三对角的三条对角线。

数值解 vs 解析解
运算次数对比 — 追赶法 O(n) vs 高斯消元 O(n³)
理论与主要公式

$$a_i\,u_{i-1} + b_i\,u_i + c_i\,u_{i+1} = d_i$$

三对角系统的第 i 行。a:下对角,b:主对角,c:上对角,d:右端项。中心差分离散时 a=c=−1、b=2。

$$c'_i=\frac{c_i}{b_i-a_i c'_{i-1}},\qquad d'_i=\frac{d_i-a_i d'_{i-1}}{b_i-a_i c'_{i-1}}$$

前向消元(forward sweep)。从 i=1 的 c'_1=c_1/b_1、d'_1=d_1/b_1 开始,依次推进 i=2..n。

$$u_n=d'_n,\qquad u_i=d'_i-c'_i\,u_{i+1}$$

回代(back substitution),按逆序 i=n−1..1 求解。总运算量为 O(n),相比一般矩阵的高斯消元 O(n³) 快得多。

什么是追赶法

🙋
"追赶法"我还是第一次听说。我知道它是解联立方程的方法,但它和普通的高斯消元有什么不同呢?
🎓
简单说,它就是高斯消元的"特快版"。追赶法针对的是三对角矩阵——只在主对角线,以及紧挨它上方、下方这三条对角线上有数字,其余全是零的特殊矩阵。一般高斯消元什么矩阵都能解,但计算量按 n³ 爆炸增长。而对三对角矩阵,零的部分完全不用碰,所以只需 O(n),大约 8n 次运算就能解出来。
🙋
只有三条对角线……这么凑巧的矩阵,实际中真会出现吗?
🎓
出现得多到让你吃惊。这个工具求解的问题就是。把一维微分方程 −u''=f 用中心差分这种标准方法离散后,每个节点的方程里出现的未知数只有三个:"自己"、"左邻"和"右邻"。写成矩阵就自动变成三对角的。你在左边改变网格点数 n 时,工具每次都用追赶法求解那个 n×n 的三对角系统。
🙋
原来如此!那具体是怎么求解的呢?右上角的动画里有箭头从左往右动。
🎓
分两个阶段。首先是"前向消元"——从最左边的方程开始,逐个消去下对角项。每行一边计算 c'_i = c_i/(b_i − a_i·c'_{i-1}) 一边向右推进,这就是动画里左→右的扫描。到了最后一行,那一行只剩一个未知数,u_n 立即确定。接着是"回代",从右→左:用 u_i = d'_i − c'_i·u_{i+1},借助右侧刚求出的值,逐个向左求解。这就是动画里右→左的扫描。
🙋
解出来的答案有多准确呢?"最大误差"显示的是 1e-13 之类极小的数字。
🎓
好问题。对于默认的均匀源,这个问题有精确的解析解 u(x) = 50·x·(1−x)。追赶法不是近似,它是"精确地"求解那个离散系统,所以除去舍入误差,误差几乎为零——只有约 1e-13,即浮点数的极限。中央的 u 应当如理论所示正好是 12.5。但若把源切换为"中央集中"或"线性",解析解的公式就变了,误差栏便显示"—"。均匀分布这种情形,正是用来确认追赶法本身是精确的。
🙋
既然又快又准,那是不是全部都用追赶法就好了……我会这么想。
🎓
心情可以理解,但它只能在"矩阵是三对角时"才能用。而且追赶法不做选主元,所以除非矩阵对角占优(对角元素不小于上下对角之和),否则消元途中分母可能趋近于零,误差就会失控。幸好这个工具的矩阵对角是 2、上下对角是 −1,属于对角占优,可以放心用。实务中只要凑齐"三对角+对角占优",就毫不犹豫地用追赶法,这是常规做法。

常见问题

追赶法是专门求解三对角线性系统 A·u = d 的算法,三对角矩阵指仅在主对角线及其上下相邻对角线上有非零元素的矩阵。它是把一般高斯消元针对这种带状结构特化而来的,由前向消元和回代两个阶段组成。一般高斯消元约需 (2/3)n³ 次运算,而追赶法仅需约 8n 次,即 O(n)。三对角系统在工程中无处不在——一维边值问题、三次样条插值、偏微分方程的隐式时间推进——因此它是数值计算中极其重要的基础算法。
在一般高斯消元中,每一步消元都要更新其下方所有行和列,因此计算量按 n³ 增长。而三对角矩阵每行最多只有三个非零元素,消去某一行只影响相邻的一行。因此前向消元每行只需常数次运算,总量与 n 成正比。回代同样每行只需一次减法和一次乘法,所以也是 O(n)。合计约 8n 次浮点运算,内存也只需保存三条对角线,即 O(n)。
追赶法不做选主元,因此当矩阵是对角占优(每行对角元素的绝对值不小于上下对角之和)或对称正定时数值稳定。本工具求解的一维泊松问题矩阵对角元素为 2、上下对角为 −1,属于对角占优,可稳定求解。而对于非对角占优的一般三对角矩阵,前向消元过程中分母 m = b_i − a_i·c'_{i-1} 可能趋近于零,从而放大误差;这时需要使用带选主元的带状矩阵求解器。
最典型的是一维边值问题:像本工具求解的 −u''=f 那样,用中心差分离散二阶导数必然得到三对角矩阵。此外还有三次样条插值的系数确定、热传导或扩散方程隐式(如克兰克-尼科尔森法)时间积分的每一步、二维问题 ADI 法的分裂、三对角结构的结构分析等,应用非常广泛。一旦出现三对角系统,首先考虑能否用追赶法以 O(n) 求解,是常规做法。

实际应用

偏微分方程的隐式时间积分:在时间方向推进热传导方程或扩散方程时,为保证稳定性而采用隐式方法(如克兰克-尼科尔森法),就需要在每个时间步求解一次三对角系统。一万步的模拟意味着要用追赶法求解一万次。这里若改用 O(n³) 的高斯消元,计算时间会膨胀数千倍,因此追赶法是让隐式方法变得实用的关键。

三次样条插值:在 CAD、图形学与数据可视化中用于绘制光滑曲线的三次样条,需要求解一组方程来确定各分段的系数。由样条"二阶导数连续"这一条件产生的方程组,恰好是三对角的。即使是穿过数百个点的光滑曲线,用追赶法也能瞬间确定系数。

一维结构与传热分析:梁的挠度、肋片的温度分布、地下热传导等,一维边值问题正是追赶法的主场。本工具求解的 −u''=f 正是这一典型,相当于泊松方程、拉普拉斯方程的一维版本。在 CAE 教学与验证中,通常先用这个简单问题确认离散化与求解器的正确性,再推进到二维、三维。

作为大型稀疏求解器的部件:即便是二维、三维问题,像 ADI 法(交替方向隐式法)那样把问题逐方向分裂,内部也归结为求解大量三对角系统。此外,块三对角求解器以及预处理迭代法的预处理部分也用到了追赶法的思想。它不仅作为独立求解器使用,也是大型求解器的核心。

常见误解与注意事项

首先常见的是"追赶法能安全求解任何三对角矩阵"这种误解。追赶法完全不做选主元(不交换行),所以一旦前向消元的分母 m = b_i − a_i·c'_{i-1} 趋近于零,误差就会瞬间被放大。只有当矩阵对角占优或对称正定时,才能保证不发生这种情况。本工具的一维泊松矩阵(对角 2、上下对角 −1)属于对角占优,因而安全;但把一般三对角矩阵原样灌进去是有风险的。非对角占优时,请使用带选主元的带状矩阵求解器(如 LAPACK 的 dgtsv)。

其次,误以为"因为是 O(n),所以误差也为零"。追赶法(除去舍入误差)精确地求解离散后的线性系统,但这意味着"离散方程的答案"是精确的——这与"原微分方程的答案"并不相同。中心差分的离散有 O(h²) 的截断误差,网格点数 n 越大,越接近真解。本工具在均匀源情形下误差几乎为零,是因为对这个特定问题,离散解与解析解恰好一致的特殊情形。一般而言,必须区分"算法的舍入误差"与"离散化的截断误差"。

最后,"只要是三对角,追赶法就一定最快"并不尽然。作为直接法,追赶法确实是最优的,但若系数矩阵要在多个右端项中复用,只要先做一次 LU 分解(同样是 O(n)),就能反复执行前向、回代代入。此外,由周期边界条件在"角上"出现元素的循环三对角矩阵,纯粹的追赶法无法求解,需要 Sherman-Morrison 公式之类的修正。仔细观察问题的结构,在追赶法、块追赶法与循环版之间灵活选用,才是实务中正确的态度。

使用指南

  1. 设置网格点数n(推荐10~100),划分求解区间[0,1]为n-1个等距单元
  2. 输入边界条件:左端点ul与右端点ur的Dirichlet值
  3. 选择源项函数f(x)类型(如f=1、f=x(1-x)等),确定泊松方程-u''=f的右端项
  4. 点击求解,模拟器执行追赶法三对角分解:先向前消元生成L、U因子,再逐步回代求解
  5. 对比图表显示数值解曲线与解析解偏差,误差指标反映离散化精度

具体计算示例

设n=21,区间[0,1]网格步长h=0.05。边界条件u(0)=0、u(1)=0,源项f(x)=1。离散化后三对角方程组系数:主对角线a=400(=1/h²),上下对角线b=c=-200。追赶法需约8×21=168次乘除法运算。数值求解得中央点u(0.5)≈0.0625,解析解u(x)=0.5x(1-x)在x=0.5处值为0.125,最大误差约0.005(h²阶收敛)。

实务注意事项

  1. 网格密集化时步长h减半,三对角系数随h²变化,追赶法数值稳定性依赖条件数κ(A),对角占优矩阵(|a|>|b|+|c|)保证无需主元交换
  2. 热传导、结构梁弯曲等1D偏微分方程离散后均产生三对角系统,追赶法相比LU分解节省87.5%存储空间
  3. 源项f(x)非光滑(如分段常数)时,数值解收敛阶降至一阶,宜增加网格点数至50+补偿
  4. 边界条件从Dirichlet改为Neumann时,需修改最后一行方程,追赶法算法框架不变