三對角線性方程組(tridiagonal systems of equations)的求解

三對角線性方程組(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算法并不是對所有的三對角矩陣都是有效的，隻是部分三對角矩陣可行。

三對角線性方程組(tridiagonal systems of equations)的求解