在下面的代码中,我对一般的平方线性系统Ax=b实现了带有部分旋转的高斯消去。我测试了我的代码,它产生了正确的输出。不过,现在我正在尝试做以下事情,但我不太确定如何编码它,寻找一些帮助与此!

我想通过求解Ax=b来测试我的实现,其中A是随机的100x100矩阵,b是随机的100x10向量。

在我的代码中,我把矩阵

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))