Сравните результат исключения по Гауссу с выводом 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".