高斯消元法——解有唯一解的方程组

import numpy as np
def gauss_elimination(a,b):
    rows , columns = np.shape(a)
    matrix = np.concatenate((a,b) ,axis=1)
    matrix = matrix.astype("float64")
    #部分选主元，从第一列到倒数第二列
    for p in range(rows):
        #找到列中绝对值最大的行
        max_row = p + np.argmax(np.abs(matrix[p:, p]))
        # 交换最大行和当前行
        matrix[[p, max_row]] = matrix[[max_row, p]]
        #消元部分
    for m in range(rows):
        temp =matrix[m,m]#找到了每一行的主对角线元素
        for n in range(m+1 , columns):#在第j中取出第i行的元素
            l = matrix[n,m]/temp
            #变化第n即(m+1)行元素，为将a(i+1)j元素变化为0
            matrix[n,:] = matrix[n,:] - l*matrix[m,:]
        print(matrix)
        print("\n")
    #回代部分
    #解的初始化，列表类型，方便追加各个分量
    result=[]
    #将增广矩阵的类型从numpy数组转化为list
    matrix=list(matrix)
    matrix.reverse()#对增广矩阵进行逆序索引
    for k in range(rows):#用k遍历增广矩阵的行
        if k==0:#由于矩阵翻转，代表为最后一行
            result.append(matrix[0][-1]/matrix[0][-2])#xn的值
        else:
            sum = 0
            for x in range(k):
            #计算第k个方程中与第k个变量无关的其他项
                sum = sum + result[x]*matrix[k][-x-2]
            result.append((matrix[k][-1]-sum)/matrix[k][-k-2])
    result.reverse()#将列表倒序表示
    print("\n")
    print(result)

a =[[1,1,0,-4],[-1,1,1,3],[1,3,5,-4],[0,1,2,-1]]
b = [[1],[-2],[-4],[-2]]
gauss_elimination(a,b)