Numpy 数组在高斯消除程序中无法正确更新

Question

我正在尝试编写一个函数

gaussian_elim

，它接受一个 n x n numpy 数组 A 和一个 n x 1 numpy 数组 b 并对增广矩阵 [A|b] 执行高斯消除。它应返回 M = [U|c] 形式的 n x (n+1) 矩阵，其中 U 是 n x n 上三角矩阵。然而，当我在一个简单的 2x2 矩阵上测试我的代码时，似乎消除步骤没有正确执行。我插入了打印语句来说明矩阵如何未正确更新。

def gaussian_elim(A,b):

    """
    A: n x n numpy array
    b: n x 1 numpy array 

    Applies Gaussian Elimination to the system Ax = b.
    Returns a matrix of the form M = [U|c], where U is upper triangular.
    """
    
    n = len(b) 
    b = b.reshape(-1, 1)        # Ensure b is a column vector of shape (n, 1)
    M = np.hstack((A,b))        #Form the n x (n+1) augmented matrix M := [A|b]
    
    #For each pivot:
    for j in range(n-1):  #j = 0,1,...,n-2
    
        #For each row under the pivot:
        for i in range(j+1,n):               #i = j + 1, j + 2,..., n-1

            if (M[j,j] == 0):
                print("Error! Zero pivot encountered!")
                return 

            #The multiplier for the the i-th row 
            m = M[i,j] / M[j,j]
 
            print("M[i,:] = ", M[i,:])
            print("M[j,:] = ", M[j,:])
            print("m = ", m)
            print("M[i,:] - m*M[j,:] = ", M[i,:] - m*M[j,:])

            #Eliminate entry M[i,j] (the first nonzero entry of the i-th row)
            M[i,:] = M[i,:] - m*M[j,:]
            
            print("M[i,:] = ", M[i,:])   #Make sure that i-th row of M is correct (it's not!)
       
    return M

使用 2x2 矩阵进行测试

A = np.array([[3,-2],[1,5]])
b = np.array([1,1])
gaussian_elim(A,b)

产生以下输出：

M[i,:] =  [1 5 1]
M[j,:] =  [ 3 -2  1]
m =  0.3333333333333333
M[i,:] - m*M[j,:] =  [0.         5.66666667 0.66666667]    <-- this is correct!
M[i,:] =  [0 5 0]                                          <--why is this not equal to the above line???

array([[ 3, -2,  1],
       [ 0,  5,  0]])

我期望的输出是

array([[ 3, -2,  1],[0. 5.66666667 0.66666667]])

。为什么第二行没有正确更新？

Answer 1

因为您使用 numpy 整数数组，所以所有结果都将四舍五入为整数。您需要将

A

和

b

定义为

A = np.array([[3,-2],[1,5]], dtype=np.float64)
b = np.array([1,1], dtype=np.float64)

这样做可以在矩阵中使用浮点值。