列主元高斯的python實現

相對於順序高斯只是每次循環的時候增加了一個選擇列主元的過程。

選擇列主元也就是找到余下的列中最大的一行，並以此行為主元

代碼如下：

# coding: utf8
import numpy as np


def getInput1():
    matrix_a = np.mat([[0.726, 0.81, 0.9],
                       [1, 1, 1],
                       [1.331, 1.21, 1.1]
                       ], dtype=float)
    matrix_b = np.mat([0.6867, 0.8338, 1])
    # 答案：-2 0 1 1
    return matrix_a, matrix_b


# 設置矩陣
def getInput():
    matrix_a = np.mat([[2, 3, 11, 5],
                     [1, 1, 5, 2],
                     [2, 1, 3, 2],
                     [1, 1, 3, 4]], dtype=float)
    matrix_b = np.mat([2, 1, -3, -3])
    # 答案：-2 0 1 1
    return matrix_a, matrix_b


# 交換
def swap(mat, num):
    print(num)
    print("調換前")
    print(mat)
    maxid = num
    for j in range(num, mat.shape[0]):
        if mat[j, num] > mat[num, num]:
            maxid = j
    if maxid is not num:
        mat[[maxid,num],:] = mat[[num,maxid],:]
    else:pass
    print("調換后")
    print(maxid)
    print(mat)
    return mat


def SequentialGauss(mat_a):
    for i in range(0, (mat_a.shape[0])-1):
        swap(mat_a, i)
        if mat_a[i, i] == 0:
            print("終斷運算：")
            print(mat_a)
            break
        else:
            for j in range(i+1, mat_a.shape[0]):
                mat_a[j:j+1 , :] = mat_a[j:j+1,:] - \
                                                    (mat_a[j,i]/mat_a[i,i])*mat_a[i, :]
    return mat_a


# 回帶過程
def revert(new_mat):
    # 創建矩陣存放答案 初始化為0
    x = np.mat(np.zeros(new_mat.shape[0], dtype=float))
    number = x.shape[1]-1
    # print(number)
    b = number+1
    x[0,number] = new_mat[number,b]/new_mat[number, number]
    for i in range(number-1,-1,-1):
        try:
            x[0, i] = (new_mat[i,b]-np.sum(np.multiply(new_mat[i,i+1:b],x[0,i+1:b])))/(new_mat[i,i])
        except:
            print("錯誤")
    print(x)


if __name__ == "__main__":
    mat_a, mat_b = getInput()
    # 合並兩個矩陣
    print("原矩陣")
    print(np.hstack((mat_a, mat_b.T)))
    new_mat = SequentialGauss(np.hstack((mat_a, mat_b.T)))
    print("三角矩陣")
    print(new_mat)
    print("方程的解")
    revert(new_mat)