1. 来源和背景

对于一个(主)三对角方程组,我们常用“追赶法”来进行求解. 而三对角方程组常常出现于微分方程的数值求解,例如热传导方程的边值问题 

 

是一个线性函数时,对该边值问题的数值解转化为一个典型的三对角方程组求解.

 

“追赶法”目前比较可靠的来源是下面的文章: 
Thomas, L.H., Elliptic Problems in Linear Differential Equations over a Network. Watson Science Computer Laboratory Report, 1949. 
其中的一个依据是,在国外的文章和教材中,“追赶法”被称为“Thomas算法”.


2. 追赶法的基本原理

追赶法的基本原理是矩阵的LU分解,即将矩阵分解为 

 

其中,为一个对角线上元素为的下三角矩阵,为一个上三角矩阵. 容易验证,一个三对角矩阵作LU分解以后,得到一个下二对角矩阵与一个上二对角矩阵的乘积,即 
 

 

 

 

三对角矩阵的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
  • 1
  • 2
  • 3
  • 4

在计算过程中,将下三角矩阵和上三角矩阵的值保存在原矩阵中. 计算结束以后,矩阵中的元素为 

 

 

注: 三对角矩阵做LU分解以后,严格上三角部分的元素没有发生变化,即上三角矩阵中的元素 

 

 


3. 追赶法求解三对角方程组

使用LU分解的求解线性方程组时,不需要存储下三角矩阵,而上三角矩阵将被用于回代求解.

3.1 “追”的过程:分解

对于阶的三对角方程组 

 

我们先用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);     b(i) = b(i) - b(i-1) * A(i,i-1); end
  • 1
  • 2
  • 3
  • 4
  • 5

循环里面的前两行与LU分解完全相同,第三行负责对常数项做相应的变换. 在计算过程中,上三角矩阵的值保存在原矩阵中,变换后的常数保存在中.

3.2 “赶”的过程:回代

接着,我们用回代法求解上三角形方程组. 从三对角矩阵得到的上三角形方程组如下: 

 

注意在前面的计算过程中,我们将上三角矩阵保存在中,常数项保存在中. 因此,我们得到如下的回代过程:

 

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
  • 1
  • 2
  • 3
  • 4

4. 实用的程序代码

在三对角矩阵中,三对角线以外的元素均为,为了提高存储的效率,我们只需存储三对角线上的元素即可. 因此,对于前面的矩阵,我们只存储三个向量:

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)];
  • 1
  • 2
  • 3

这三个向量分别为矩阵三条对角线上的元素. 假定常数向量为

b=[b(1),b(2),...,b(n)];
  • 1

则实用的追赶法(亦称为“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