三对角线性方程组(tridiagonal systems of equations)
三对角线性方程组,对于熟悉数值分析的同学来说，并不陌生，它经常出现在微分方程的数值求解和三次样条函数的插值问题中。三对角线性方程组可描述为以下方程组：
aixi−1+bixi+cixi+1=diaixi−1+bixi+cixi+1=di
其中1≤i≤n,a1=0,cn=0.1≤i≤n,a1=0,cn=0. 以上方程组写成矩阵形式为Ax=dAx=d，即：
⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢b1a20c1b2a3c2b3⋱⋱⋱an0cn−1bn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢x1x2x3⋮xn⎤⎦⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢d1d2d3⋮dn⎤⎦⎥⎥⎥⎥⎥⎥⎥[b1c10a2b2c2a3b3⋱⋱⋱cn−10anbn][x1x2x3⋮xn]=[d1d2d3⋮dn]
三对角线性方程组的求解采用追赶法或者Thomas算法，它是Gauss消去法在三对角线性方程组这种特殊情形下的应用，因此，主要思想还是Gauss消去法，只是会更加简单些。我们将在下面的算法详述中给出该算法的具体求解过程。
当然，该算法并不总是稳定的，但当系数矩阵AA为严格对角占优矩阵（Strictly D iagonally Dominant, SDD）或对称正定矩阵（Symmetric Positive Definite, SPD）时，该算法稳定。对于不熟悉SDD或者SPD的读者，也不必担心，我们还会在我们的博客中介绍这类矩阵。现在，我们只要记住，该算法对于部分系数矩阵AA是可以求解的。
算法详述
追赶法或者Thomas算法的具体步骤如下：
- 创建新系数c∗ici∗和d∗idi∗来代替原先的ai,bi,ciai,bi,ci,公式如下： c∗i={c1b1cibi−aic∗i−1;i=1;i=2,3,...,n−1d∗i=⎧⎩⎨⎪⎪d1b1di−aid∗i−1bi−aic∗i−1;i=1;i=2,3,...,n−1ci∗={c1b1;i=1cibi−aici−1∗;i=2,3,...,n−1di∗={d1b1;i=1di−aidi−1∗bi−aici−1∗;i=2,3,...,n−1
- 改写原先的方程组Ax=dAx=d如下： ⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢100...0c∗110...00c∗21000c∗30......00000..c∗n−11⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢x1x2x3...xk⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢d∗1d∗2d∗3...d∗n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥[1c1∗00...001c2∗0...0001c3∗00........cn−1∗000001][x1x2x3...xk]=[d1∗d2∗d3∗...dn∗]
- 计算解向量xx，如下： xn=d∗n,xi=d∗i−c∗ixi+1,i=n−1,n−2,...,2,1xn=dn∗,xi=di∗−ci∗xi+1,i=n−1,n−2,...,2,1
以上算法得到的解向量xx即为原方程Ax=dAx=d的解。
下面，我们来证明该算法的正确性，只需要证明该算法保持原方程组的形式不变。
首先，当i=1i=1时，
1∗x1+c∗1x2=d∗1⇔1∗x1+c1b1x2=d1b1⇔b1∗x1+c1x2=d11∗x1+c1∗x2=d1∗⇔1∗x1+c1b1x2=d1b1⇔b1∗x1+c1x2=d1
当i>1i>1时，
1∗xi+c∗ixi+1=d∗i⇔1∗xi+cibi−aic∗i−1xi+1=di−aid∗i−1bi−aic∗i−1⇔(bi−aic∗i−1)xi+cixi+1=di−aid∗i−11∗xi+ci∗xi+1=di∗⇔1∗xi+cibi−aici−1∗xi+1=di−aidi−1∗bi−aici−1∗⇔(bi−aici−1∗)xi+cixi+1=di−aidi−1∗
结合aixi−1+bixi+cixi+1=diaixi−1+bixi+cixi+1=di，只需要证明xi−1+c∗i−1xi=d∗i−1xi−1+ci−1∗xi=di−1∗,而这已经在该算法的第（3）步的中的计算xi−1xi−1中给出。证明完毕。
Python实现
我们将要求解的线性方程组如下：
⎡⎣⎢⎢⎢⎢⎢⎢4100014100014100014100014⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢x1x2x3x4x5⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢10.5−132⎤⎦⎥⎥⎥⎥[4100014100014100014100014][x1x2x3x4x5]=[10.5−132]
接下来，我们将用Python来实现该算法，函数为TDMA，输入参数为列表a,b,c,d, 输出为解向量x，代码如下：
# use Thomas Method to solve tridiagonal linear equation
# algorithm reference: https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
import numpy as np
# parameter: a,b,c,d are list-like of same length
# tridiagonal linear equation: Ax=d
# b: main diagonal of matrix A
# a: main diagonal below of matrix A
# c: main diagonal upper of matrix A
# d: Ax=d
# return: x(type=list), the solution of Ax=d
def TDMA(a,b,c,d):
try:
n = len(d) # order of tridiagonal square matrix
# use a,b,c to create matrix A, which is not necessary in the algorithm
A = np.array([[0]*n]*n, dtype='float64')
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]
# new list of modified coefficients
c_1 = [0]*n
d_1 = [0]*n
for i in range(n):
if not i:
c_1[i] = c[i]/b[i]
d_1[i] = d[i] / b[i]
else:
c_1[i] = c[i]/(b[i]-c_1[i-1]*a[i])
d_1[i] = (d[i]-d_1[i-1]*a[i])/(b[i]-c_1[i-1] * a[i])
# x: solution of Ax=d
x = [0]*n
for i in range(n-1, -1, -1):
if i == n-1:
x[i] = d_1[i]
else:
x[i] = d_1[i]-c_1[i]*x[i+1]
x = [round(_, 4) for _ in x]
return x
except Exception as e:
return e
def main():
a = [0, 1, 1, 1, 1]
b = [4, 4, 4, 4, 4]
c = [1, 1, 1, 1, 0]
d = [1, 0.5, -1, 3, 2]
'''
a = [0, 2, 1, 3]
b = [1, 1, 2, 1]
c = [2, 3, 0.5, 0]
d = [2, -1, 1, 3]
'''
x = TDMA(a, b, c, d)
print('The solution is %s'%x)
main()运行该程序，输出结果为：
The solution is [0.2, 0.2, -0.5, 0.8, 0.3]
本算法的Github地址为： https://github.com/percent4/Numeric_Analysis/blob/master/TDMA.py .
最后再次声明，追赶法或者Thomas算法并不是对所有的三对角矩阵都是有效的，只是部分三对角矩阵可行。