步驟1 檢查A，b是否行數相同

步驟2 構造增廣矩陣Ab

步驟3 逐列轉換Ab為化簡行階梯形矩陣 中文維基連結

對於Ab的每一列（最後一列除外）
    當前列為列c
    尋找列c中 對角線以及對角線以下所有元素（行 c~N）的絕對值的最大值
    如果絕對值最大值為0
        那麼A為奇異矩陣，返回None (你可以在選做問題2.4中證明為什麼這裡A一定是奇異矩陣)
    否則
        使用第一個行變換，將絕對值最大值所在行交換到對角線元素所在行（行c） 
        使用第二個行變換，將列c的對角線元素縮放為1
        多次使用第三個行變換，將列c的其他元素消為0
 

步驟4 返回Ab的最後一列

注： 我們並沒有按照常規方法先把矩陣轉化為行階梯形矩陣，再轉換為化簡行階梯形矩陣，而是一步到位。如果你熟悉常規方法的話，可以思考一下兩者的等價性。

from decimal import Decimal
from fractions import Fraction
# TODO 實現 Gaussain Jordan 方法求解 Ax = b

""" Gaussian Jordan 方法求解 Ax = b.
    引數
        A: 方陣 
        b: 列向量
        decPts: 四捨五入位數，預設為4
        epsilon: 判讀是否為0的閾值，預設 1.0e-16
        
    返回列向量 x 使得 Ax = b 
    返回None，如果 A，b 高度不同
    返回None，如果 A 為奇異矩陣
"""

def gj_Solve(A, b, decPts=4, epsilon = 1.0e-16):
    if len(A)!= len(b) :
        return None
#     構造擴增矩陣
    result = augmentMatrix(A, b)
#     轉換為階梯型矩陣
    for j in range(len(result[0])-1):
        row,maxNum = 0,0  
        for i in range(len(result)):
            if i >= j:
                if abs(result[i][j]) > maxNum:
                    maxNum = abs(result[i][j])
                    row = i
#         返回奇異矩陣
        if(abs(maxNum) < epsilon):
            return None
        else:
#         找出最大放到最上面
            if(row != j):
               swapRows(result,j,row)
               row,maxNum = 0,0 
        scaleRow(result,j,Fraction(1,result[j][j]))
#         消除下面  
        for dis in range(len(result)):              
            if dis != j:
                if abs(-result[dis][j]) > epsilon:
                    addScaledRow(result,dis,j,-result[dis][j])
        newResult =[]
        for i in range(len(result)):
            row =[]
            row.append(result[i][len(result)])
            newResult.append(row)       
    return newResult
 

測試

    def test_gj_Solve(self):

        for _ in range(9999):
            r = np.random.randint(low=3,high=4)
            A = np.random.randint(low=-10,high=10,size=(r,r))
            b = np.arange(r).reshape((r,1))
            x = gj_Solve(A.tolist(),b.tolist(),epsilon=1.0e-8)
            
         
            if np.linalg.matrix_rank(A) < r:
                self.assertEqual(x,None,"Matrix A is singular")
            else:
                self.assertNotEqual(x,None,"Matrix A is not singular")
                self.assertEqual(np.array(x).shape,(r,1),"Expected shape({},1), but got shape{}".format(r,np.array(x).shape))
                Ax = np.dot(A,np.array(x))
                loss = np.mean((Ax - b)**2)
                self.assertTrue(loss<0.1,"Bad result.")