高斯消去法的改進形式為Gauss-Jordan Elimination Method，要求每一行的主元素所在列元素全部消去為0，除了主元素本身，區別如下：
目錄：1 Gauss-Jordan法的程式、2 考慮病態方程的情況（某兩行系數的比例幾乎完全相同）
代碼實作如下：
# -*- coding: utf-8 -*- # @Author : ZhaoKe # @Time : 2022-09-05 23:34 from typing import List # input a augmented matrix, output its simpler form def GaussJordanMethod(matrix: List[List]) -> List[List]: PREC = 3 m, n = len(matrix), len(matrix[0]) # from left to right for j in range(n): # 主元素為0的話，要交換行令這一行主元素不為0 if j >= m: break if matrix[j][j] == 0: for k in range(j + 1, m): if matrix[k][j] == 0: continue else: matrix[j], matrix[k] = matrix[k], matrix[j] print(matrix) break # from above to bottom # 為了實作Gauss-Jordan，pivot位置的元素應該置為1 fac = matrix[j][j] for l in range(0, n): matrix[j][l] = matrix[j][l] / fac # 不再拘泥于對角線下方消除，整整一列都要消除 for i in range(m): # 主元素不可以消去，直接跳過該行 if i == j: continue # 當前行的該列元素為0的話，不必執行消去步驟，跳過即可 if matrix[i][j] == 0: continue # replace the jth equation by a combination of itself plus a multiple of the ith equation coef = matrix[i][j] for k in range(n): matrix[i][k] = round(matrix[i][k] - coef * matrix[j][k], PREC) print(matrix[i]) # elimination end print(matrix) # solution as follow for i in range(m - 1, -1, -1): for j in range(n - 1, -1, -1): if matrix[i][j] == 0: continue else: print("x" + str(i+1), "=", matrix[i][-1]) break return matrix if __name__ == '__main__': # input_m0 = [[2, 1, 1, 1], [6, 2, 1, -1], [-2, 2, 1, 7]] # 結果正確 # input_m1 = [[0, 1, -1, 3], [-2, 4, -1, 1], [-2, 5, -4, -2]] # 結果正確 input_m2 = [[2, 2, 6, 4], [2, 1, 7, 6], [-2, -6, -7, -1]] # 結果正確 GaussJordanMethod(input_m2) # 病態方程： # input_m3 = [[47, 28, 19], [89, 53, 36]] # 該程式暫不適合病態方程的數值計算，另開隨筆研究該問題 # GaussJordanMethod(input_m3)
經驗證該程式正確，但是對于病態方程可能會出現例外：
Traceback (most recent call last): matrix[j][l] = matrix[j][l] / fac ZeroDivisionError: float division by zero
修改如下：
# -*- coding: utf-8 -*- # @Author : ZhaoKe # @Time : 2022-09-05 23:34 from typing import List # input a augmented matrix, output its simpler form def GaussJordanMethod(matrix: List[List], prec=3) -> List[List]: m, n = len(matrix), len(matrix[0]) # from left to right for j in range(n): # 主元素為0的話，要交換行令這一行主元素不為0 if j >= m: break if matrix[j][j] == 0: for k in range(j + 1, m): if matrix[k][j] == 0: continue else: matrix[j], matrix[k] = matrix[k], matrix[j] print("exchange ", matrix) break # from above to bottom # 為了實作Gauss-Jordan，pivot位置的元素應該置為1 print("row", j, matrix[j]) # 如果是病態方程，會出現消去之后整行都為0的情況 # 判斷如果pivot element為0，那么直接結束程式 fac = matrix[j][j] if fac == 0: break for l in range(0, n): matrix[j][l] = matrix[j][l] / fac print("row", j, matrix[j]) # 不再拘泥于對角線下方消除，整整一列都要消除 for i in range(m): # 主元素不可以消去，直接跳過該行 if i == j: continue # 當前行的該列元素為0的話，不必執行消去步驟，跳過即可 if matrix[i][j] == 0: continue # replace the jth equation by a combination of itself plus a multiple of the ith equation coef = matrix[i][j] # / matrix[j][j] matrix[j][j]必為1（根據Gauss-Jordan法） # if coef == 0: # continue print("coef:", coef) for k in range(n): matrix[i][k] = round(matrix[i][k] - coef * matrix[j][k], prec) print("-> ", matrix[i]) # elimination end print(matrix) # solution as follow for i in range(m - 1, -1, -1): for j in range(n - 1, -1, -1): if matrix[i][j] == 0: continue else: print("x" + str(i+1), "=", matrix[i][-1]) break return matrix if __name__ == '__main__': # input_m0 = [[2, 1, 1, 1], [6, 2, 1, -1], [-2, 2, 1, 7]] # 結果正確 # input_m1 = [[0, 1, -1, 3], [-2, 4, -1, 1], [-2, 5, -4, -2]] # 結果正確 # input_m2 = [[2, 2, 6, 4], [2, 1, 7, 6], [-2, -6, -7, -1]] # 結果正確 # GaussJordanMethod(input_m1) # 病態方程： # input_m3 = [[47, 28, 19], [89, 53, 36]] input_m4 = [[0.835, 0.667, 0.168], [0.333, 0.266, 0.067]] GaussJordanMethod(input_m4, prec=5) print("=================") input_m4 = [[0.835, 0.667, 0.168], [0.333, 0.266, 0.067]] GaussJordanMethod(input_m4, prec=6)
有效數字為5的時候，經過高斯消去的矩陣為：
[[1.0, 0.798802395209581, 0.20119760479041918], [0.0, -0.0, 0.0]]
顯然無解
有效數字為6的時候，消去得：
[[1.0, 0.0, 1.0], [-0.0, 1.0, -1.0]]
x2 = -1.0
x1 = 1.0
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/504443.html
標籤：其他