例題7-1 利用Thomas 理則 求解對角線聯立方程式

例題7-1 利用Thomas 理則 求解

2x1  -   x2                       =1
-x1   +2x2    -x3             = 0
         -  x2    +2x3  - x4  = 0
                     -x3    +2x4 =1


左上又下 斜線方向 對角線聯立方程式

a= [ 0 ,  0 , -1 , -1 , -1 ]  從a[2]開始
b= [ 0 ,  2 , 2  , 2  , 2   ]   從b[1]開始
c= [ 0,  -1 , -1  , -1 ]  從c[1]開始
d= [ 0 ,  1 , 0 , 0 1 ]   從d[1]開始

K is a tridiagonal matrix:

    K =   [   \   \   \       ]
             |    \   \   \      ]
             |     a   b   c   ]
             |      \   \   \    ]
             |       \   \   \   ]

    Where,
            a: lower diagonal
            b: main diagonal
            c: upper diagonal
            d: RHS or the input vector
            x: solution vector



'''
 Tri Diagonal Matrix Algorithm(a.k.a Thomas algorithm) solver
 Uses Thomas algorithm for solving a tridiagonal matrix for n unknowns.
      a, b, and  are a list of the matrix entries
      Matrix form of:
      [b1 c1        ] [x1]   [d1]
      |a2 b2 c2   | |x2|   |d2|
      | a3 b3 c3  | |x3|   |d3|
      |                  | |'  | = | '|
      |                  | |'  |     | '|
      [ an bn cn] |xn] [dn]

    # Finding the size of the matrix and determine n
    #print n # Used for debugging
    # Test the size of a and c

This code solves the system Ku = f by gaussian elimination.
This particular algorithm is often referred to as the 'Thomas Algorithm'.

K is a tridiagonal matrix:

    K =   [   \   \   \       ]
             |    \   \   \      ]
             |     a   b   c   ]
             |      \   \   \    ]
             |       \   \   \   ]

    Where,
            a: lower diagonal
            b: main diagonal
            c: upper diagonal
            d: RHS or the input vector
            x: solution vector
'''

def tridiagonal_solver( n , a, b,  c,  d , x):
    for i in range (2 , n+1):
        r=a[i]/b[i-1]
        b[i]=b[i]-r*c[i-1]
        d[i]=d[i]-r*d[i-1]

    # /* The answers are stored in x[i] */
    x[n]=d[n]/b[n];

    for i in range (n-1 , 0 , -1):
        x[i]=(d[i]-c[i]*x[i+1]) / b[i]

    return x;


 #++++++++++++++++++++++++++++++++++++++++++

n=4
a = [0.0 , 0.0 , -1.0 , -1.0 , -1.0]
b = [0.0 ,2.0 , 2.0 , 2.0 , 2.0 ]
c = [0.0 , -1.0 , -1.0 , -1.0 , 0.0 ]
d = [0.0 , 1.0 ,  0.0 , 0.0  , 1.0 ]

x= [0.0 for i in range(n+1)]     #a [n][n]//係數矩陣
s= [0.0 for i in range(n+1)]     #a [n][n]//係數矩陣

s=(tridiagonal_solver(n, a,b,c,d,x))
for i in range(1 ,n+1):
    print( "x[",i,"]=", round(s[i],4),'\t',end='')

========= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX7-1.py ============
x[ 1 ]= 1.0 x[ 2 ]= 1.0 x[ 3 ]= 1.0 x[ 4 ]= 1.0
>>>