比较高斯消去的结果和numpy.linalg.solve的输出

Question

在下面的代码中，我为一般的方形线性系统 Ax=b 实现了具有部分旋转的高斯消去法。我已经测试了我的代码，它产生了正确的输出。我用它来求解 Ax=b，其中 A 是一个随机的 100x100 矩阵，b 是一个随机的 100x1 向量。**

但是现在我正在寻求一些帮助，将我的解决方案与使用 numpy.linalg.solve 获得的解决方案进行比较。如何将此比较添加到我的代码中？

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.round(np.random.rand(100, 100)*10)
    b =  np.round(np.random.rand(100)*10)
    print (GEPP(np.copy(A), np.copy(b), doPricing = False))

Answer 1

要比较这两个解决方案，请使用 np.allclose 来验证两个数组是否一致，逐个元素，在一些合理的范围内。人们不应该期望浮点运算中完全相等。

my_solution = GEPP(np.copy(A), np.copy(b), doPricing=False)
numpy_solution = np.linalg.solve(A, b)
if np.allclose(my_solution, numpy_solution):
    print("Yay")
else:
    print("NOOOOO")

这会打印 "Yay".

比较高斯消去的结果和numpy.linalg.solve的输出

Compare the result of Gaussian elimination with the output of numpy.linalg.solve

python

numpy

numerical-methods