三对角线性方程组(tridiagonal systems of equations)
三对角线性方程组,对于熟悉数值分析的同学来说,并不陌生,它经常出现在微分方程的数值求解和三次样条函数的插值问题中。三对角线性方程组可描述为以下方程组:
aixi−1+bixi+cixi+1=di
a
i
x
i
−
1
+
b
i
x
i
+
c
i
x
i
+
1
=
d
i
其中 1≤i≤n,a1=0,cn=0. 1 ≤ i ≤ n , a 1 = 0 , c n = 0. 以上方程组写成矩阵形式为 Ax=d A x = d ,即:
⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢b1a20c1b2a3c2b3⋱⋱⋱an0cn−1bn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢x1x2x3⋮xn⎤⎦⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢d1d2d3⋮dn⎤⎦⎥⎥⎥⎥⎥⎥⎥
[
b
1
c
1
0
a
2
b
2
c
2
a
3
b
3
⋱
⋱
⋱
c
n
−
1
0
a
n
b
n
]
[
x
1
x
2
x
3
⋮
x
n
]
=
[
d
1
d
2
d
3
⋮
d
n
]
三对角线性方程组的求解采用追赶法或者Thomas算法,它是Gauss消去法在三对角线性方程组这种特殊情形下的应用,因此,主要思想还是Gauss消去法,只是会更加简单些。我们将在下面的算法详述中给出该算法的具体求解过程。
当然,该算法并不总是稳定的,但当系数矩阵 A A 为严格对角占优矩阵(Strictly D iagonally Dominant, SDD)或对称正定矩阵(Symmetric Positive Definite, SPD)时,该算法稳定。对于不熟悉SDD或者SPD的读者,也不必担心,我们还会在我们的博客中介绍这类矩阵。现在,我们只要记住,该算法对于部分系数矩阵是可以求解的。
算法详述
追赶法或者Thomas算法的具体步骤如下:
- 创建新系数
c∗i
c
i
∗
和
d∗i
d
i
∗
来代替原先的
ai,bi,ci
a
i
,
b
i
,
c
i
,公式如下:
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−1 c i ∗ = { c 1 b 1 ; i = 1 c i b i − a i c i − 1 ∗ ; i = 2 , 3 , . . . , n − 1 d i ∗ = { d 1 b 1 ; i = 1 d i − a i d i − 1 ∗ b i − a i c i − 1 ∗ ; i = 2 , 3 , . . . , n − 1 - 改写原先的方程组
Ax=d
A
x
=
d
如下:
⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢100...0c∗110...00c∗21000c∗30......00000..c∗n−11⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢x1x2x3...xk⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢d∗1d∗2d∗3...d∗n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥ [ 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 ] [ x 1 x 2 x 3 . . . x k ] = [ d 1 ∗ d 2 ∗ d 3 ∗ . . . d n ∗ ] - 计算解向量
x
x
,如下:
以上算法得到的解向量
x
x
即为原方程的解。
下面,我们来证明该算法的正确性,只需要证明该算法保持原方程组的形式不变。
首先,当
i=1
i
=
1
时,
1∗x1+c∗1x2=d∗1⇔1∗x1+c1b1x2=d1b1⇔b1∗x1+c1x2=d1
1
∗
x
1
+
c
1
∗
x
2
=
d
1
∗
⇔
1
∗
x
1
+
c
1
b
1
x
2
=
d
1
b
1
⇔
b
1
∗
x
1
+
c
1
x
2
=
d
1
当 i>1 i > 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−1
1
∗
x
i
+
c
i
∗
x
i
+
1
=
d
i
∗
⇔
1
∗
x
i
+
c
i
b
i
−
a
i
c
i
−
1
∗
x
i
+
1
=
d
i
−
a
i
d
i
−
1
∗
b
i
−
a
i
c
i
−
1
∗
⇔
(
b
i
−
a
i
c
i
−
1
∗
)
x
i
+
c
i
x
i
+
1
=
d
i
−
a
i
d
i
−
1
∗
结合 aixi−1+bixi+cixi+1=di a i x i − 1 + b i x i + c i x i + 1 = d i ,只需要证明 xi−1+c∗i−1xi=d∗i−1 x i − 1 + c i − 1 ∗ x i = d i − 1 ∗ ,而这已经在该算法的第(3)步的中的计算 xi−1 x i − 1 中给出。证明完毕。
Python实现
我们将要求解的线性方程组如下:
⎡⎣⎢⎢⎢⎢⎢⎢4100014100014100014100014⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢x1x2x3x4x5⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢10.5−132⎤⎦⎥⎥⎥⎥
[
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
]
[
x
1
x
2
x
3
x
4
x
5
]
=
[
1
0.5
−
1
3
2
]
接下来,我们将用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算法并不是对所有的三对角矩阵都是有效的,只是部分三对角矩阵可行。