93 lines
2.5 KiB
Python
93 lines
2.5 KiB
Python
|
||
# 列主元高斯消元法
|
||
def SovleRowMain(A,b):
|
||
ks = 0.00000001
|
||
n = len(A)
|
||
if len(A[0]) != n:
|
||
raise ValueError("A要为方阵")
|
||
if len(b) != n:
|
||
raise ValueError("b与A的行数不匹配")
|
||
p = list(range(n))
|
||
for i in range(n):
|
||
row_max = abs(A[i][i])
|
||
row_max_index = i
|
||
for j in range(i + 1, n):
|
||
if abs(A[j][i]) > row_max:
|
||
row_max = abs(A[j][i])
|
||
row_max_index = j
|
||
A[i], A[row_max_index] = A[row_max_index], A[i]
|
||
b[i], b[row_max_index] = b[row_max_index], b[i]
|
||
p[i], p[row_max_index] = p[row_max_index], p[i]
|
||
|
||
print(f"第{i+1}次交换，交换行{i+1}和行{row_max_index+1}，A矩阵为：")
|
||
for row in A:
|
||
print(row)
|
||
print(f"b向量为：{b}\n")
|
||
if abs(A[i][i]) < ks:
|
||
raise ValueError("A矩阵奇异，无法进行高斯消元")
|
||
for j in range(i + 1, n):
|
||
m = A[j][i] / A[i][i]
|
||
A[j][i] = m
|
||
for k in range(i + 1, n):
|
||
A[j][k] -= m * A[i][k]
|
||
b[j] -= m * b[i]
|
||
print(f"系数为{-1*m}用加号")
|
||
print("消元后的A矩阵：")
|
||
for row in A:
|
||
print(row)
|
||
print(f"消元后的b向量：{b}\n")
|
||
if abs(A[n - 1][n - 1]) < ks:
|
||
raise ValueError("A矩阵奇异，无法进行高斯消元")
|
||
|
||
# 回代求解
|
||
b[n - 1] /= A[n - 1][n - 1]
|
||
for i in range(n - 2, -1, -1):
|
||
for j in range(i + 1, n):
|
||
b[i] -= A[i][j] * b[j]
|
||
b[i] /= A[i][i]
|
||
|
||
# 得到L,U和P矩阵
|
||
L = [[0 for i in range(n)] for j in range(n)]
|
||
U = [[0 for i in range(n)] for j in range(n)]
|
||
P = [[0 for i in range(n)] for j in range(n)]
|
||
for i in range(n):
|
||
for j in range(n):
|
||
if i == j:
|
||
L[i][j] = 1
|
||
U[i][j] = A[i][j]
|
||
elif i < j:
|
||
U[i][j] = A[i][j]
|
||
else:
|
||
L[i][j] = A[i][j]
|
||
|
||
P[i][p[i]] = 1
|
||
return P,L,U,b
|
||
#打印矩阵
|
||
def prettyPrintMatrix(matrix):
|
||
for row in matrix:
|
||
print(row)
|
||
|
||
|
||
if __name__ == "__main__":
|
||
# A = np.array(A, dtype=float)
|
||
# b = np.array(b, dtype=float)
|
||
#把矩阵A和b改成题干要求的#####################################
|
||
A = [
|
||
[0, 3, 4],
|
||
[1, -1, 1],
|
||
[2, 1, 2]
|
||
]
|
||
|
||
b = [1, 2, 3]
|
||
|
||
P,L,U,x = SovleRowMain(A, b)
|
||
print("P:")
|
||
prettyPrintMatrix(P)
|
||
print("L:")
|
||
prettyPrintMatrix(L)
|
||
print("U:")
|
||
prettyPrintMatrix(U)
|
||
print("x:")
|
||
print(x)
|
||
|
||
|