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

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

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

aixi1+bixi+cixi+1=diaixi−1+bixi+cixi+1=di

其中1in,a1=0,cn=0.1≤i≤n,a1=0,cn=0. 以上方程组写成矩阵形式为Ax=dAx=d,即:
b1a20c1b2a3c2b3an0cn1bnx1x2x3xn=d1d2d3dn[b1c10a2b2c2a3b3⋱⋱⋱cn−10anbn][x1x2x3⋮xn]=[d1d2d3⋮dn]

  三对角线性方程组的求解采用追赶法或者Thomas算法,它是Gauss消去法在三对角线性方程组这种特殊情形下的应用,因此,主要思想还是Gauss消去法,只是会更加简单些。我们将在下面的算法详述中给出该算法的具体求解过程。
  当然,该算法并不总是稳定的,但当系数矩阵AA为严格对角占优矩阵(Strictly D iagonally Dominant, SDD)或对称正定矩阵(Symmetric Positive Definite, SPD)时,该算法稳定。对于不熟悉SDD或者SPD的读者,也不必担心,我们还会在我们的博客中介绍这类矩阵。现在,我们只要记住,该算法对于部分系数矩阵AA是可以求解的。

算法详述

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

  1. 创建新系数cici∗didi∗来代替原先的ai,bi,ciai,bi,ci,公式如下:
    ci={c1b1cibiaici1;i=1;i=2,3,...,n1di=d1b1diaidi1biaici1;i=1;i=2,3,...,n1ci∗={c1b1;i=1cibi−aici−1∗;i=2,3,...,n−1di∗={d1b1;i=1di−aidi−1∗bi−aici−1∗;i=2,3,...,n−1
  2. 改写原先的方程组Ax=dAx=d如下:
    100...0c110...00c21000c30......00000..cn11x1x2x3...xk=d1d2d3...dn[1c1∗00...001c2∗0...0001c3∗00........cn−1∗000001][x1x2x3...xk]=[d1∗d2∗d3∗...dn∗]
  3. 计算解向量xx,如下:
    xn=dn,xi=dicixi+1,i=n1,n2,...,2,1xn=dn∗,xi=di∗−ci∗xi+1,i=n−1,n−2,...,2,1

  以上算法得到的解向量xx即为原方程Ax=dAx=d的解。
  下面,我们来证明该算法的正确性,只需要证明该算法保持原方程组的形式不变。
  首先,当i=1i=1时,

1x1+c1x2=d11x1+c1b1x2=d1b1b1x1+c1x2=d11∗x1+c1∗x2=d1∗⇔1∗x1+c1b1x2=d1b1⇔b1∗x1+c1x2=d1

  当i>1i>1时,
1xi+cixi+1=di1xi+cibiaici1xi+1=diaidi1biaici1(biaici1)xi+cixi+1=diaidi11∗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∗

结合aixi1+bixi+cixi+1=diaixi−1+bixi+cixi+1=di,只需要证明xi1+ci1xi=di1xi−1+ci−1∗xi=di−1∗,而这已经在该算法的第(3)步的中的计算xi1xi−1中给出。证明完毕。

Python实现

  我们将要求解的线性方程组如下:

4100014100014100014100014x1x2x3x4x5=10.5132[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算法并不是对所有的三对角矩阵都是有效的,只是部分三对角矩阵可行。

参考文献

  1. https://www.quantstart.com/articles/Tridiagonal-Matrix-Solver-via-Thomas-Algorithm
  2. https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
  3. https://wenku.baidu.com/view/336bafa3daef5ef7ba0d3ccc.html
上一篇:jQuery技术内幕:深入解析jQuery架构设计与实现原理. 2.8 静态属性和方法


下一篇:企业级业务架构设计方法与“中台”概念的比较