问题描述
为求解一个线性方程组，首先构造增广矩阵[A|B]，采用偏序选主元策略的高斯消去法变换成上三角矩阵，再执行回代过程得到解。
输入形式
在屏幕上依次输入方阵阶数n，系数矩阵A和常数矩阵B。
输出形式
首先输出上三角矩阵（变换后的增广矩阵），然后每一行输出一个根
样例1输入
4
1 2 1 4
2 0 4 3
4 2 2 1
-3 1 3 2
13
28
20
6
样例1输出
[[ 4. 2. 2. 1. 20. ]
[ 0. 2.5 4.5 2.75 21. ]
[ 0. 0. 4.8 3.6 26.4 ]
[ 0. 0. 0. 3.75 7.5 ]]
[[ 3.]
[-1.]
[ 4.]
[ 2.]]
样例1说明
输入：第1行为方阵阶数4，第2行至5行为系数矩阵A，第6行至9行为常数矩阵B。
输出：首先输出上三角矩阵（变换后的增广矩阵），然后每行依次输出方程解：x1, x2, x3, x4。
代码
# 高斯消元法
import numpy as np
def Input():
n = int(input())
A = np.zeros([n, n], dtype=np.double())
for r in range(n):
A[r:] = np.array(input().split(), dtype=np.double)
B = np.zeros([n, 1], dtype=np.double)
for r in range(n):
B[r:] = np.array(input(), dtype=np.double)
return A, B, n
def Gaussian(Aug):
n = Aug.shape[0]
epslion = np.finfo(np.float32).eps
for p in range(0, n):
max_index = np.argmax(np.fabs(Aug[p:n, p]))
# 交换p, max_row行
max_row = max_index + p
Aug[[p, max_row], :] = Aug[[max_row, p], :]
if abs(Aug[p, p]) < epslion :
print('Array a was singular. No unique solution')
break
for k in range(p + 1, n):
m = Aug[k, p] / Aug[p, p]
Aug[k, p:n + 1] = Aug[k, p:n + 1] - m * Aug[p, p:n + 1]
return Aug
def backsub(A, B):
n = B.size
X = np.zeros([n, 1], dtype=np.double)
X[n - 1] = B[n - 1] / A[n - 1][n - 1]
for i in range(n - 2, -1, -1):
X[i] = (B[i] - A[i, i + 1:] @ X[i + 1:]) / A[i][i]
return X
def main():
A, B, n = Input()
gaussian = Gaussian(np.hstack((A, B))) # 将A,B以列合并
print(gaussian)
X = backsub(gaussian[:, 0:n], gaussian[:, n:])
print(X)
if __name__ == '__main__':
main()
