【数值分析实验】线性代数方程组的直接解法：列主元高斯消去法、LU分解法、LU分解法求矩阵的逆（python）

调包

import numpy as np

消元函数和矩阵求解

#消元
def KillYuan(A,i):
    try:
        A[i] /= A[i][i]
    except:
        print("矩阵奇异")
    A[i] /= A[i][i]
    for j in range(i + 1,len(A)):
        A[j] -= A[j][i] * A[i]
    return A

#矩阵求解
def MatrixSolve(A):
    A = np.copy(A)
    b = A[:,-1]
    a = np.delete(A,-1,1)
    try:
        np.linalg.solve(a,b)
    except:
        print("矩阵奇异")
    return np.linalg.solve(a,b)

高斯消去法

#高斯消元法
def  GaussianElimination(A):
    A = np.copy(A)
    m = len(A)
    n = len(A[0])
    for i in range(m):
        A = KillYuan(A,i)
    return MatrixSolve(A),A

列主元高斯消去法

#列主元高斯消元法
def MainYuanGaussianElimination(A):
    A = np.copy(A)
    m = len(A)
    n = len(A[0])
    for i in range(m):
        maxYuan = 0
        for j in range(i,m):
            if abs(A[j][i]) > maxYuan:
                maxYuan = j
        A[[i,maxYuan],:] = A[[maxYuan,i],:]    
        A = KillYuan(A,i)
    return MatrixSolve(A),A

LU分解法

#LU分解法
def LU(A):
    n = len(A[0])
    L = np.zeros([n, n])
    U = np.zeros([n, n])
    for i in range(n):
        L[i][i] = 1
        if i == 0:
            for j in range(n):
                U[i][j] = A[i][j]
                L[j][i] = A[j][i] / U[i][i]
        else:
            for j in range(i, n):  
                temp = 0
                for k in range(i):
                    temp = temp + L[i][k] * U[k][j]
                U[i][j] = A[i][j] - temp
            for j in range(i + 1, n):  
                temp = 0
                for k in range(i):
                    temp = temp + L[j][k] * U[k][i]
                L[j][i] = (A[j][i] - temp) / U[i][i]
    return L,U

LU分解求矩阵的逆

下三角矩阵求逆

#下三角矩阵求逆
def DownTriangle_Inv(L):
    n = len(L)
    Li = np.zeros([n,n])
    for i in range(n):
        for j in range(n):
            if i < j:
                Li[i][j] = 0
            if i == j:
                Li[i][j] = 1 / L[i][j]
            if i > j:
                s = 0
                for k in range(j,i):
                    s += L[i][k] * Li[k][j]
                Li[i][j] = -1 / L[i][i] * s
    return Li

LU分解求逆

#LU分解法矩阵求逆
def A_Inv(A):
    L, U = LU(A)
    Li = DownTriangle_Inv(L)
    Ui = DownTriangle_Inv(U.T).T
    return np.dot(Ui,Li)

演示

Ae2 = np.array([[6, 2, -1, 2], 
                [5, 0, 4, -10], 
                [-9, -4, 2, 0], 
                [6, 8, 2, -10]])
print("LU分解法求矩阵的逆：")
print(np.around(A_Inv(Ae2), 2))
print("numpy验证：")
print(np.around(np.linalg.inv(Ae2), 2))