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

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

三对角线性方程组,对于熟悉数值分析的同学来说，并不陌生，它经常出现在微分方程的数值求解和三次样条函数的插值问题中。三对角线性方程组可描述为以下方程组：

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,a1=0,cn=0. 以上方程组写成矩阵形式为

Ax=dAx=d，即：

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ b 1 a 2 0 c 1 b 2 a 3 c 2 b 3 ⋱ ⋱ ⋱ a n 0 c n - 1 b n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ x 1 x 2 x 3 ⋮ x n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ d 1 d 2 d 3 ⋮ d n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

三对角线性方程组的求解采用追赶法或者Thomas算法，它是Gauss消去法在三对角线性方程组这种特殊情形下的应用，因此，主要思想还是Gauss消去法，只是会更加简单些。我们将在下面的算法详述中给出该算法的具体求解过程。

当然，该算法并不总是稳定的，但当系数矩阵

AA为严格对角占优矩阵（Strictly D iagonally Dominant, SDD）或对称正定矩阵（Symmetric Positive Definite, SPD）时，该算法稳定。对于不熟悉SDD或者SPD的读者，也不必担心，我们还会在我们的博客中介绍这类矩阵。现在，我们只要记住，该算法对于部分系数矩阵

AA是可以求解的。

算法详述

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

创建新系数c∗i和d∗i来代替原先的ai,bi,ci,公式如下：

$c * i = {c 1 b 1 c i b i - a i c * i - 1; i = 1; i = 2, 3, . . ., n - 1 d * i = ⎧ ⎩ ⎨ ⎪ ⎪ d 1 b 1 d i - a i d * i - 1 b i - a i c * i - 1; i = 1; i = 2, 3, . . ., n - 1$

改写原先的方程组Ax=d如下：

$⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 100 . . . 0 c * 1 10 . . . 00 c * 2 1000 c * 3 0 . . . . . . 00000 . . c * n - 1 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ x 1 x 2 x 3 . . . x k ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ d * 1 d * 2 d * 3 . . . d * n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥$

计算解向量x，如下：

$x n = d * n, x i = d * i - c * i x i + 1, i = n - 1, n - 2, . . ., 2, 1$

以上算法得到的解向量x即为原方程Ax=d的解。

下面，我们来证明该算法的正确性，只需要证明该算法保持原方程组的形式不变。

首先，当

i=1i=1时，

1 * x 1 + c * 1 x 2 = d * 1 \Leftrightarrow 1 * x 1 + c 1 b 1 x 2 = d 1 b 1 \Leftrightarrow b 1 * x 1 + c 1 x 2 = d 1

当

i>1i>1时，

1 * x i + c * i x i + 1 = d * i \Leftrightarrow 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 \Leftrightarrow (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=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 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ x 1 x 2 x 3 x 4 x 5 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ 1 0.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_algorithmimport 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=ddef 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 edef 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地址为： .

最后再次声明，追赶法或者Thomas算法并不是对所有的三对角矩阵都是有效的，只是部分三对角矩阵可行。