线性代数方程组的直接解法
调包
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))