例題 7-2 THOMES 理則解 三個對角線的聯立方程式

例題 7-2 THOMES 理則解 三個對角線的聯立方程式

tridiagonal matrix

'''
 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;


 #++++++++++++++++++++++++++++++++++++++++++
'''
# a 要從 a[2] 開始 其他的 b[1] c[1] d[1]

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


4
0.35  0.25  1
0.5   0.8   1    -2
0.25  0.4   0.5
0.35  0.77 -0.5  -2.25
'''
n=4
a = [0.0 , 0.0  , 0.35 , 0.25 , 1.0 ]
b = [0.0 , 0.5 , 0.8  ,1.0  , -2.0 ]
c =  [0.0 , 0.25 , 0.4 , 0.5  ]
d = [0.0 , 0.35 , 0.77, -0.5 , -2.25 ]

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-0.py ==========
x[ 1 ]= 1.0 x[ 2 ]= -1.0 x[ 3 ]= 0.0
>>> 
========= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX7-1-0.py =========
x[ 1 ]= -0.0936 x[ 2 ]= 1.5872 x[ 3 ]= -1.1674 x[ 4 ]= 0.5413
>>>