三对角线性方程组(tridiagonal systems of equations)
三对角线性方程组,对于熟悉数值分析的同学来说，并不陌生，它经常出现在微分方程的数值求解和三次样条函数的插值问题中。三对角线性方程组可描述为以下方程组：
\]
其中\(1\leq i \leq n, a_{1}=0, c_{n}=0.\) 以上方程组写成矩阵形式为\(Ax=d\)，即：
{b_{1}}&{c_{1}}&{}&{}&{0}\\
{a_{2}}&{b_{2}}&{c_{2}}&{}&{}\\
{}&{a_{3}}&{b_{3}}&\ddots &{}\\
{}&{}&\ddots &\ddots &{c_{n-1}}\\
{0}&&&{a_{n}}&{b_{n}}\\
\end{bmatrix}}
{\begin{bmatrix}{x_{1}}\\{x_{2}}\\{x_{3}}\\\vdots \\{x_{n}}\\\end{bmatrix}}={\begin{bmatrix}{d_{1}}\\{d_{2}}\\{d_{3}}\\\vdots \\{d_{n}}\\\end{bmatrix}}
\]
三对角线性方程组的求解采用追赶法或者Thomas算法，它是Gauss消去法在三对角线性方程组这种特殊情形下的应用，因此，主要思想还是Gauss消去法，只是会更加简单些。我们将在下面的算法详述中给出该算法的具体求解过程。
当然，该算法并不总是稳定的，但当系数矩阵\(A\)为严格对角占优矩阵（Strictly D iagonally Dominant, SDD）或对称正定矩阵（Symmetric Positive Definite, SPD）时，该算法稳定。对于不熟悉SDD或者SPD的读者，也不必担心，我们还会在我们的博客中介绍这类矩阵。现在，我们只要记住，该算法对于部分系数矩阵\(A\)是可以求解的。
算法详述
追赶法或者Thomas算法的具体步骤如下：
- 创建新系数\(c_{i}^{*}\)和\(d_{i}^{*}\)来代替原先的\(a_{i},b_{i},c_{i}\),公式如下：
\begin{array}{lr}
\frac{c_1}{b_1} & ; i = 1\\
\frac{c_i}{b_i - a_i c^{*}_{i-1}} & ; i = 2,3,...,n-1
\end{array}
\right.\\
d^{*}_i = \left\{
\begin{array}{lr}
\frac{d_1}{b_1} & ; i = 1\\
\frac{d_i- a_i d^{*}_{i-1}}{b_i - a_i c^{*}_{i-1}} & ; i = 2,3,...,n-1
\end{array}
\right.
\]
- 改写原先的方程组\(Ax=d\)如下：
1 & c^{*}_1 & 0 & 0 & ... & 0 \\
0 & 1 & c^{*}_2 & 0 & ... & 0 \\
0 & 0 & 1 & c^{*}_3 & 0 & 0 \\
. & . & & & & . \\
. & . & & & & . \\
. & . & & & & c^{*}_{n-1} \\
0 & 0 & 0 & 0 & 0 & 1 \\
\end{bmatrix} \begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
.\\
.\\
.\\
x_k \\
\end{bmatrix} = \begin{bmatrix}
d^{*}_1 \\
d^{*}_2 \\
d^{*}_3 \\
.\\
.\\
.\\
d^{*}_n \\
\end{bmatrix}
\]
- 计算解向量\(x\)，如下：
\]
以上算法得到的解向量\(x\)即为原方程\(Ax=d\)的解。
下面，我们来证明该算法的正确性，只需要证明该算法保持原方程组的形式不变。
首先，当\(i=1\)时，
\]
当\(i>1\)时，
(b_{i} - a_{i} c^{*}_{i-1})x_{i}+c_{i}x_{i+1}=d_{i}- a_{i} d^{*}_{i-1}\]
结合\(a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i}\)，只需要证明\(x_{i-1}+c_{i-1}^{*}x_{i}=d_{i-1}^{*}\),而这已经在该算法的第（3）步的中的计算\(x_{i-1}\)中给出。证明完毕。
Python实现
我们将要求解的线性方程组如下：
4&1&{0}&{0}&{0}\\
{1}&{4}&{1}&{0}&{0}\\
{0}&{1}&{4}&{1}&{0}\\
{0}&{0}&{1}&{4}&{1}\\
{0}&{0}&{0}&{1}&{4}\\
\end{bmatrix}}
{\begin{bmatrix}{x_{1}}\\{x_{2}}\\{x_{3}}\\{x_{4}} \\{x_{5}}\\\end{bmatrix}}={\begin{bmatrix}{1\\0.5\\ -1\\3\\2}\\\end{bmatrix}}
\]
接下来，我们将用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算法并不是对所有的三对角矩阵都是有效的，只是部分三对角矩阵可行。