追赶法求解三对角方程组
2024-07-061. 来源和背景
对于一个（主）三对角方程组，我们常用“追赶法”来进行求解. 而三对角方程组常常出现于微分方程的数值求解，例如热传导方程的边值问题
{y′′(x)=f(x,y,y′), a≤x≤by(a)=η1, y(b)=η2
当
f(x,y,y′)
是一个线性函数时，对该边值问题的数值解转化为一个典型的三对角方程组求解.
“追赶法”目前比较可靠的来源是下面的文章： Thomas, L.H., Elliptic Problems in Linear Differential Equations over a Network. Watson Science Computer Laboratory Report, 1949. 其中的一个依据是，在国外的文章和教材中，“追赶法”被称为“Thomas算法”.
2. 追赶法的基本原理
追赶法的基本原理是矩阵的LU分解，即将矩阵
A
分解为
A=LU 其中，
L
为一个对角线上元素为1的下三角矩阵，
U
为一个上三角矩阵. 容易验证，一个三对角矩阵作LU分解以后，得到一个下二对角矩阵与一个上二对角矩阵的乘积，即
A=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢a11a21a12a22a32a23a33⋱a34⋱⋱an−1,n−2an−1,n−1an,n−1an−1,nan−1,n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
L=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢1ℓ211ℓ321⋱⋱ℓn−1,n−21ℓn,n−11⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
U=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢u11u12u22u23u33u34⋱⋱un−1,n−1un−1,nun−1,n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
三对角矩阵
A
的LU分解计算过程如下：
for i = 2 to n
A(i,i-1) = A(i,i-1)/A(i-1,i-1);
A(i,i) = A(i,i) - A(i-1,i) * A(i,i-1);
end
在计算过程中，将下三角矩阵L和上三角矩阵
U
的值保存在原矩阵A中. 计算结束以后，矩阵
A
中的元素为
⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢u11ℓ21u12u22ℓ32u23u33⋱u34⋱⋱ℓn−1,n−2un−1,n−1ℓn,n−1un−1,nun−1,n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
注： 三对角矩阵
A
做LU分解以后，严格上三角部分的元素没有发生变化，即上三角矩阵U中的元素
ui,i+1=ai,i+1, i=1,2,…,n−1
3. 追赶法求解三对角方程组
使用LU分解的求解线性方程组时，不需要存储下三角矩阵，而上三角矩阵将被用于回代求解.
3.1 “追”的过程：分解
对于
n
阶的三对角方程组
Ax=b 我们先用LU分解得到
Ux=y
注： 由
Ax=LUx=b
，得
Ux=L−1b
记
y=L−1b
，即得到方程组
Ux=y
.
计算过程如下：
for i = 2 to n
A(i,i-1) = A(i,i-1)/A(i-1,i-1);
A(i,i) = A(i,i) - A(i-1,i) * A(i,i-1);
b(i) = b(i) - b(i-1) * A(i,i-1);
end
循环里面的前两行与LU分解完全相同，第三行负责对常数项做相应的变换. 在计算过程中，上三角矩阵
U
的值保存在原矩阵A中，变换后的常数
y=L−1b
保存在
b
中.
3.2 “赶”的过程：回代
接着，我们用回代法求解上三角形方程组. 从三对角矩阵得到的上三角形方程组如下：
⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢u11u12u22u23u33u34⋱⋱un−1,n−1un−1,nun−1,n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎛⎝⎜⎜⎜⎜⎜⎜⎜x1x2⋮xn−1xn⎞⎠⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜y1y2⋮yn−1yn⎞⎠⎟⎟⎟⎟⎟⎟⎟ 注意在前面的计算过程中，我们将上三角矩阵
U
保存在A中，常数项
y
保存在b中. 因此，我们得到如下的回代过程：
x(n) = b(n) / A(i,i);
for i = n-1 to 1
x(i) = (b(i) - A(i,i+1) * x(i+1)) / A(i,i);
end
4. 实用的程序代码
在三对角矩阵中，三对角线以外的元素均为
0
，为了提高存储的效率，我们只需存储三对角线上的元素即可. 因此，对于前面的矩阵A，我们只存储三个向量：
d=[A(1,1),A(2,2),...,A(n,n)];
u=[A(1,2),A(2,3),...,A(n-1,n)];
l=[A(2,1),A(3,2),...,A(n,n-1)];
这三个向量分别为矩阵
A
三条对角线上的元素. 假定常数向量为
b=[b(1),b(2),...,b(n)];
则实用的追赶法（亦称为“Thomas算法”）求解三对角方程组的过程如下：
% 追
for i = 2 to n
l(i-1) = l(i-1)/d(i-1);
d(i,i) = d(i,i) - u(i-1) * l(i-1);
b(i) = b(i) - b(i-1) * l(i-1);
end
% 赶
x(n) = b(n) / d(i);
for i = n-1 to 1
x(i) = (b(i) - u(i) * x(i+1)) / d(i);
end