PYthon 高斯列主元消去法求增广矩阵/方程组的解 Numpy模块
“对 A * X = B” 矩阵的阶数不限
可通过修改 a 的值来改变A
a = np.array([[2, 1, 1], [3, 1, 2], [1, 2, 2]],dtype=float)
可通过修改 b 的值来修改B
b = np.array([[4],[6],[5]],dtype=float)
模块导入
import numpy as np
直接上代码
# 导入 numpy 模块
import numpy as np
# 行交换
def swap_row(matrix, i, j):
m, n = matrix.shape
if i >= m or j >= m:
print('错误! : 行交换超出范围 ...')
else:
matrix[i],matrix[j] = matrix[j].copy(),matrix[i].copy()
return matrix
# 变成阶梯矩阵
def matrix_change(matrix):
m, n = matrix.shape
main_factor = []
main_col = main_row = 0
while main_row < m and main_col < n:
# 选择进行下一次主元查找的列
main_row = len(main_factor)
# 寻找列中非零的元素
not_zeros = np.where(abs(matrix[main_row:,main_col]) > 0)[0]
# 如果该列向下全部数据为零,则直接跳过列
if len(not_zeros) == 0:
main_col += 1
continue
else:
# 将主元列号保存在列表中
main_factor.append(main_col)
# 将第一个非零行交换至最前
if not_zeros[0] != [0]:
matrix = swap_row(matrix,main_row,main_row+not_zeros[0])
# 将该列主元下方所有元素变为零
if main_row < m-1:
for k in range(main_row+1,m):
a = float(matrix[k, main_col] / matrix[main_row, main_col])
matrix[k] = matrix[k] - matrix[main_row] * matrix[k, main_col] / matrix[main_row, main_col]
main_col += 1
return matrix,main_factor
# 回代求解
def back_solve(matrix, main_factor):
# 判断是否有解
if len(main_factor) == 0:
print('主元错误,无主元! ...')
return None
m, n = matrix.shape
if main_factor[-1] == n - 1:
print('无解! ...')
return None
# 把所有的主元元素上方的元素变成0
for i in range(len(main_factor) - 1, -1, -1):
factor = matrix[i, main_factor[i]]
matrix[i] = matrix[i] / float(factor)
for j in range(i):
times = matrix[j, main_factor[i]]
matrix[j] = matrix[j] - float(times) * matrix[i]
# 先看看结果对不对
return matrix
# 结果打印
def print_result(matrix, main_factor):
if matrix is None:
print('阶梯矩阵为空! ...')
return
m, n = matrix.shape
result = [''] * (n - 1)
main_factor = list(main_factor)
for i in range(n - 1):
# 如果不是主元列,则为自由变量
if i not in main_factor:
result[i] = 'x_' + str(i + 1) + '(free var)'
# 否则是主元变量,从对应的行,将主元变量表示成非主元变量的线性组合
else:
# row_of_main表示该主元所在的行
row_of_main = main_factor.index(i)
result[i] = str(matrix[row_of_main, -1])
for j in range(i + 1, n - 1):
ratio = matrix[row_of_main, j]
if ratio > 0:
result[i] = result[i] + '-' + str(ratio) + 'x_' + str(j + 1)
if ratio < 0:
result[i] = result[i] + '+' + str(-ratio) + 'x_' + str(j + 1)
print('方程的通解是:', )
for i in range(n - 1):
print('x_' + str(i + 1), '=', result[i])
# 得到简化的阶梯矩阵和主元列
def Handle(matrix_a, matrix_b):
# 拼接成增广矩阵
matrix_01 = np.hstack([matrix_a, matrix_b])
print('增广矩阵为:')
print(matrix_01)
matrix_01, main_factor = matrix_change(matrix_01)
print('阶梯矩阵为:')
print(matrix_01)
matrix_01 = back_solve(matrix_01, main_factor)
print('方程的简化阶梯矩阵:')
print(matrix_01)
print('方程的主元列为:')
print(main_factor)
print_result(matrix_01, main_factor)
return matrix_01, main_factor
if __name__ == '__main__':
#a = np.array([[0, 1, 1], [0, 1, 0], [1, 0, 0]])
a = np.array([[2, 1, 1], [3, 1, 2], [1, 2, 2]],dtype=float)
b = np.array([[4],[6],[5]],dtype=float)
Handle(a, b)
print('*' * 20)