大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。
Jetbrains全家桶1年46，售后保障稳定
在下面的代码中,我对一般的平方线性系统Ax=b实现了带有部分旋转的高斯消去。我测试了我的代码,它产生了正确的输出。不过,现在我正在尝试做以下事情,但我不太确定如何编码它,寻找一些帮助与此!
我想通过求解Ax=b来测试我的实现,其中A是随机的100×100矩阵,b是随机的100×10向量。
在我的代码中,我把矩阵
A = np.array([[3.,2.,-4.],[2.,3.,3.],[5.,-3.,1.]])
b = np.array([[3.],[15.],[14.]])
得到以下正确的输出:
[3. 1. 2.]
[3. 1. 2.]
但是现在我如何改变它来生成随机矩阵呢?
下面是我的代码:
import numpy as np
def GEPP(A, b, doPricing = True):
”’
Gaussian elimination with partial pivoting.
input: A is an n x n numpy matrix
b is an n x 1 numpy array
output: x is the solution of Ax=b
with the entries permuted in
accordance with the pivoting
done by the algorithm
post-condition: A and b have been modified.
”’
n = len(A)
if b.size != n:
raise ValueError(“Invalid argument: incompatible sizes between”+
“A & b.”, b.size, n)
# k represents the current pivot row. Since GE traverses the matrix in the
# upper right triangle, we also use k for indicating the k-th diagonal
# column index.
# Elimination
for k in range(n-1):
if doPricing:
# Pivot
maxindex = abs(A[k:,k]).argmax() + k
if A[maxindex, k] == 0:
raise ValueError(“Matrix is singular.”)
# Swap
if maxindex != k:
A[[k,maxindex]] = A[[maxindex, k]]
b[[k,maxindex]] = b[[maxindex, k]]
else:
if A[k, k] == 0:
raise ValueError(“Pivot element is zero. Try setting doPricing to True.”)
#Eliminate
for row in range(k+1, n):
multiplier = A[row,k]/A[k,k]
A[row, k:] = A[row, k:] – multiplier*A[k, k:]
b[row] = b[row] – multiplier*b[k]
# Back Substitution
x = np.zeros(n)
for k in range(n-1, -1, -1):
x[k] = (b[k] – np.dot(A[k,k+1:],x[k+1:]))/A[k,k]
return x
if __name__ == “__main__”:
A = np.array([[3.,2.,-4.],[2.,3.,3.],[5.,-3.,1.]])
b = np.array([[3.],[15.],[14.]])
print (GEPP(np.copy(A), np.copy(b), doPricing = False))
print (GEPP(A,b))
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/215876.html原文链接：https://javaforall.net
