三對角線性方程組,對於熟悉數值分析的同窗來講,並不陌生,它常常出如今微分方程的數值求解和三次樣條函數的插值問題中。三對角線性方程組可描述爲如下方程組:
$$a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i}$$
其中$1\leq i \leq n, a_{1}=0, c_{n}=0.$ 以上方程組寫成矩陣形式爲$Ax=d$,即:html
$$ {\begin{bmatrix} {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}} $$git
三對角線性方程組的求解採用追趕法或者Thomas算法,它是Gauss消去法在三對角線性方程組這種特殊情形下的應用,所以,主要思想仍是Gauss消去法,只是會更加簡單些。咱們將在下面的算法詳述中給出該算法的具體求解過程。
固然,該算法並不老是穩定的,但當係數矩陣$A$爲嚴格對角佔優矩陣(Strictly D iagonally Dominant, SDD)或對稱正定矩陣(Symmetric Positive Definite, SPD)時,該算法穩定。對於不熟悉SDD或者SPD的讀者,也沒必要擔憂,咱們還會在咱們的博客中介紹這類矩陣。如今,咱們只要記住,該算法對於部分系數矩陣$A$是能夠求解的。github
追趕法或者Thomas算法的具體步驟以下:算法
1.建立新系數$c_{i}^{*}$和$d_{i}^{*}$來代替原先的$a_{i},b_{i},c_{i}$,公式以下:函數
$$ c^{*}_i = \left\{ \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. $$3d
2.改寫原先的方程組$Ax=d$以下:code
$$ \begin{bmatrix} 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} $$htm
3.計算解向量$x$,以下:
$$ x_n = d^{*}_n, \qquad x_i = d^{*}_i - c^{*}_i x_{i+1}, \qquad i = n-1, n-2, ... ,2,1$$ip
以上算法獲得的解向量$x$即爲原方程$Ax=d$的解。
下面,咱們來證實該算法的正確性,只須要證實該算法保持原方程組的形式不變。
首先,當$i=1$時,
$$1*x_{1}+c_{1}^{*}x_{2}=d_{1}^{*} \Leftrightarrow 1*x_{1}+\frac{c_{1}}{b_{1}}x_{2}=\frac{d_{1}}{b_{1}}\Leftrightarrow b_{1}*x_{1}+c_{1}x_{2}=d_{1}$$
當$i>1$時,ci
$$ 1*x_{i}+c_{i}^{*}x_{i+1}=d_{i}^{*} \Leftrightarrow 1*x_{i}+\frac{c_{i}}{b_{i} - a_{i} c^{*}_{i-1}}x_{i+1}=\frac{d_{i}- a_{i} d^{*}_{i-1}}{b_{i} - a_{i} c^{*}_{i-1}} \Leftrightarrow (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}$中給出。證實完畢。
咱們將要求解的線性方程組以下:
$$ {\begin{bmatrix} 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/N... .
最後再次聲明,追趕法或者Thomas算法並非對全部的三對角矩陣都是有效的,只是部分三對角矩陣可行。