三对角线性方程组(tridiagonal systems of equations)的求解

ゝ一纸荒年。 2022-05-25 05:58 206阅读 0赞

## 三对角线性方程组(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的读者，也不必担心，我们还会在我们的博客中介绍这类矩阵。现在，我们只要记住，该算法对于部分系数矩阵  A A 是可以求解的。

## 算法详述 ##

追赶法或者Thomas算法的具体步骤如下：

1.  创建新系数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 
2.  改写原先的方程组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 ∗ \] 
3.  计算解向量x  x ，如下：  
    
    
    xn=d∗n,xi=d∗i−c∗ixi\+1,i=n−1,n−2,...,2,1 x n=d n∗,x i=d i∗−c i∗xi \+ 1,i = n − 1 , n − 2 , . . . , 2 , 1

以上算法得到的解向量x  x 即为原方程Ax=d A x = d的解。  
  下面，我们来证明该算法的正确性，只需要证明该算法保持原方程组的形式不变。  
  首先，当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][https_github.com_percent4_Numeric_Analysis_blob_master_TDMA.py] .  
  最后再次声明，追赶法或者Thomas算法并不是对所有的三对角矩阵都是有效的，只是部分三对角矩阵可行。

## 参考文献 ##

1.  [https://www.quantstart.com/articles/Tridiagonal-Matrix-Solver-via-Thomas-Algorithm][https_www.quantstart.com_articles_Tridiagonal-Matrix-Solver-via-Thomas-Algorithm]
2.  [https://en.wikipedia.org/wiki/Tridiagonal\_matrix\_algorithm][https_en.wikipedia.org_wiki_Tridiagonal_matrix_algorithm]
3.  [https://wenku.baidu.com/view/336bafa3daef5ef7ba0d3ccc.html][https_wenku.baidu.com_view_336bafa3daef5ef7ba0d3ccc.html]

[https_github.com_percent4_Numeric_Analysis_blob_master_TDMA.py]: https://github.com/percent4/Numeric_Analysis/blob/master/TDMA.py
[https_www.quantstart.com_articles_Tridiagonal-Matrix-Solver-via-Thomas-Algorithm]: https://www.quantstart.com/articles/Tridiagonal-Matrix-Solver-via-Thomas-Algorithm
[https_en.wikipedia.org_wiki_Tridiagonal_matrix_algorithm]: https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
[https_wenku.baidu.com_view_336bafa3daef5ef7ba0d3ccc.html]: https://wenku.baidu.com/view/336bafa3daef5ef7ba0d3ccc.html