利用高斯消去法計算下面聯立方程式

例題6-1
利用高斯消去法計算下面聯立方程式


-3.0 x1  + 2.0  x2 -6.0 x3  = 6.0
 5.0 x1  + 7.0 x2  -5.0 x3  = 6.0
,1.0 x1  + 4.0 x2  -2.0 x3  = 8.0

-1.0 x1 + 1.0  x2  + 2.0  x3 = 2.0
 3.0 x1  - 1.0  x2  + 1.0  x3 = 6.0
-1.0 x1 +  3.0 x2  + 4.0  x3 = 4.0

 2.0 x1  - 1.5 x2  + 3.0  x3  = 1.0
-1.0 x1  + 0.0 x2 + 2.0  x3 = 3.0
 4.0 x1  - 4.5  x2 + 5.0  x3  = 1.0

 4.01  x1 + 1.23 x2   + 1.43 x3  -  0.73  x4 = 5.94
 1.23  x1 + 7.41 x2   + 2.41 x3  +  3.02 x4 = 14.07
 1.43  x1 + 2.41 x2   + 5.79  x3 -  1.11  x4 = 8.52
-0.73  x1 + 3.02 x2   -  1.11  x3 +  6.41 x4 = 7.59


程式
=============================================
def myGauss(m):
    #eliminate columns
    for col in range(len(m[0])):
        for row in range(col+1, len(m)):
            r = [(rowValue * (-(m[row][col] / m[col][col]))) for rowValue in m[col]]
            m[row] = [sum(pair) for pair in zip(m[row], r)]
    #now backsolve by substitution
    ans = []
    m.reverse() #makes it easier to backsolve
    for sol in range(len(m)):
            if sol == 0:
                ans.append(m[sol][-1] / m[sol][-2])
            else:
                inner = 0
                #substitute in all known coefficients
                for x in range(sol):
                    inner += (ans[x]*m[sol][-2-x])
                #the equation is now reduced to ax + b = c form
                #solve with (c - b) / a
                ans.append((m[sol][-1]-inner)/m[sol][-sol-2])
    ans.reverse()
    return ans

a= myGauss([[-3.0,2.0,-6.0,6.0],
                        [5.0,7.0,-5.0,6.0],
                        [1.0,4.0,-2.0,8.0]])

print("x1=",round(a[0],4), '\t',"x2=",round(a[1],4), '\t',"x3=",round(a[2],4) )


a= myGauss([[-1.0 , 1.0 , 2.0 , 2.0],
                        [3.0 , -1.0 , 1.0 , 6.0],
                        [-1.0 , 3.0 , 4.0 , 4.0]])
print("x1=",round(a[0],4), '\t',"x2=",round(a[1],4), '\t',"x3=",round(a[2],4) )
               
a= myGauss( [[2.0 , -1.5 , 3.0 , 1.0] ,
                         [-1.0 ,  0.0 , 2.0 , 3.0] ,
                         [ 4.0  ,-4.5 , 5.0 , 1.0] ])
print("x1=",round(a[0],4), '\t',"x2=",round(a[1],4), '\t',"x3=",round(a[2],4) )
           
a= myGauss( [[4.01 , 1.23 ,  1.43 , -0.73 ,  5.94],
                         [1.23 , 7.41  , 2.41 ,  3.02 , 14.07],
                         [1.43 , 2.41 ,  5.79 , -1.11 ,  8.52],
                         [-0.73 ,  3.02 , -1.11 ,   6.41 ,  7.59]])
print("x1=",round(a[0],4), '\t',"x2=",round(a[1],4), '\t',"x3=",round(a[2],4), '\t',"x3=",round(a[3],4) )

輸出結果
======= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX6-1-1.py ==========
x1= -2.0 x2= 3.0 x3= 1.0
x1= 1.0 x2= -1.0 x3= 2.0
x1= -1.0 x2= -0.0 x3= 1.0
x1= 1.0 x2= 1.0 x3= 1.0 x3= 1.0
>>>