三对角线性方程组(tridiagonal systems of equations)的求解 |
您所在的位置：网站首页 › 解方程组的解法有 › 三对角线性方程组(tridiagonal systems of equations)的求解 |
|
三对角线性方程组(tridiagonal systems of equations)
三对角线性方程组,对于熟悉数值分析的同学来说，并不陌生，它经常出现在微分方程的数值求解和三次样条函数的插值问题中。三对角线性方程组可描述为以下方程组： aixi−1+bixi+cixi+1=di 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 , a 1 = 0 , c n = 0. 以上方程组写成矩阵形式为 Ax=d A x = d ，即： ⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢b1a20c1b2a3c2b3⋱⋱⋱an0cn−1bn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢x1x2x3⋮xn⎤⎦⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢d1d2d3⋮dn⎤⎦⎥⎥⎥⎥⎥⎥⎥ [ b 1 c 1 0 a 2 b 2 c 2 a 3 b 3 ⋱ ⋱ ⋱ c n − 1 0 a n b n ] [ x 1 x 2 x 3 ⋮ x n ] = [ d 1 d 2 d 3 ⋮ d n ] 三对角线性方程组的求解采用追赶法或者Thomas算法，它是Gauss消去法在三对角线性方程组这种特殊情形下的应用，因此，主要思想还是Gauss消去法，只是会更加简单些。我们将在下面的算法详述中给出该算法的具体求解过程。 当然，该算法并不总是稳定的，但当系数矩阵 A A 为严格对角占优矩阵（Strictly D iagonally Dominant, SDD）或对称正定矩阵（Symmetric Positive Definite, SPD）时，该算法稳定。对于不熟悉SDD或者SPD的读者，也不必担心，我们还会在我们的博客中介绍这类矩阵。现在，我们只要记住，该算法对于部分系数矩阵AA是可以求解的。 算法详述追赶法或者Thomas算法的具体步骤如下： 创建新系数 c∗i c i ∗ 和 d∗i d i ∗ 来代替原先的 ai,bi,ci a i , b i , c i ,公式如下： 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−1 c i ∗ = { c 1 b 1 ; i = 1 c i b i − a i c i − 1 ∗ ; i = 2 , 3 , . . . , n − 1 d i ∗ = { d 1 b 1 ; i = 1 d i − a i d i − 1 ∗ b i − a i c i − 1 ∗ ; i = 2 , 3 , . . . , n − 1 改写原先的方程组 Ax=d A x = d 如下： ⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢100...0c∗110...00c∗21000c∗30......00000..c∗n−11⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢x1x2x3...xk⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢d∗1d∗2d∗3...d∗n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥ [ 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 ] [ x 1 x 2 x 3 . . . x k ] = [ d 1 ∗ d 2 ∗ d 3 ∗ . . . d n ∗ ] 计算解向量 x x ，如下： 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以上算法得到的解向量 x x 即为原方程Ax=dAx=d的解。 下面，我们来证明该算法的正确性，只需要证明该算法保持原方程组的形式不变。 首先，当 i=1 i = 1 时， 1∗x1+c∗1x2=d∗1⇔1∗x1+c1b1x2=d1b1⇔b1∗x1+c1x2=d1 1 ∗ x 1 + c 1 ∗ x 2 = d 1 ∗ ⇔ 1 ∗ x 1 + c 1 b 1 x 2 = d 1 b 1 ⇔ b 1 ∗ x 1 + c 1 x 2 = d 1 当 i>1 i > 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−1 1 ∗ x i + c i ∗ x i + 1 = d i ∗ ⇔ 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 ∗ ⇔ ( 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=di a i x i − 1 + b i x i + c i x i + 1 = d i ，只需要证明 xi−1+c∗i−1xi=d∗i−1 x i − 1 + c i − 1 ∗ x i = d i − 1 ∗ ,而这已经在该算法的第（3）步的中的计算 xi−1 x i − 1 中给出。证明完毕。 Python实现我们将要求解的线性方程组如下： ⎡⎣⎢⎢⎢⎢⎢⎢4100014100014100014100014⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢x1x2x3x4x5⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢10.5−132⎤⎦⎥⎥⎥⎥ [ 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 ] [ x 1 x 2 x 3 x 4 x 5 ] = [ 1 0.5 − 1 3 2 ] 接下来，我们将用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算法并不是对所有的三对角矩阵都是有效的，只是部分三对角矩阵可行。 参考文献 https://www.quantstart.com/articles/Tridiagonal-Matrix-Solver-via-Thomas-Algorithmhttps://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithmhttps://wenku.baidu.com/view/336bafa3daef5ef7ba0d3ccc.html |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |