假设理想状态下的弹道轨迹满足抛物线方程 y = a + bx + cx^2,已知弹道轨迹经过点(1,6)、(3,5)、(7,2),求此弹道轨迹抛物线方程,采用高斯列主元消去法求解,请给出用Python编写的可运行的实现代码,并在每次消去的时候输出一次矩阵,并给出每次消去时选取的列主元值、列主元行标、列主元列标
时间: 2023-05-24 10:04:30 浏览: 68
以下是用Python实现的可运行代码:
import numpy as np
# 输入矩阵A和向量b
A = np.array([[1, 1, 1], [1, 3, 9], [1, 7, 49]])
b = np.array([6, 5, 2])
# 高斯列主元消去法
n = len(b)
for k in range(n-1):
# 选取列主元
p = np.argmax(abs(A[k:, k])) + k
print(f"第{k+1}次选取的列主元值为:{A[p, k]},行标为:{p+1},列标为:{k+1}")
# 交换k行和p行
if p != k:
A[[p, k], :] = A[[k, p], :]
b[[p, k]] = b[[k, p]]
# 消元
for i in range(k+1, n):
# 计算倍数因子
m = A[i,k] / A[k,k]
A[i,k+1:n] = A[i,k+1:n] - m * A[k,k+1:n]
b[i] -= m * b[k]
print(f"第{k+1}次消元后的矩阵为:\n{np.concatenate((A, b.reshape(-1,1)), axis=1)}")
# 回代求解
x = np.zeros(n)
x[n-1] = b[n-1] / A[n-1, n-1]
for k in range(n-2, -1, -1):
x[k] = (b[k] - np.dot(A[k,k+1:n], x[k+1:n])) / A[k,k]
print(f"方程解为:{x}")