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

方程组(tridiagonal systems of equations)的求解"/>

三对角线性方程组(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算法的具体步骤如下：

1. 创建新系数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
2. 改写原先的方程组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.......−1∗000001][x1x2x3...xk]=[d1∗d2∗d3∗...dn∗]
3. 计算解向量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:  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 algorithmA = 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 coefficientsc_1 = [0]*nd_1 = [0]*nfor 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=dx = [0]*nfor 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 xexcept 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地址为： .py .
  最后再次声明，追赶法或者Thomas算法并不是对所有的三对角矩阵都是有效的，只是部分三对角矩阵可行。

参考文献

3. .html