[Python] Numpy array does not correctly update in Gaussian elimination program

I am trying to write a function gaussian_elim which takes in an n x n numpy array A and an n x 1 numpy array b and performs Gaussian elimination on the augmented matrix [A|b]. It should return an n x (n+1) matrix of the form M = [U|c], where U is an n x n upper triangular matrix. However, when I test my code on a simple 2x2 matrix, it seems that the elimination step is not being performed properly. I have inserted print statements to illustrate how the matrix is not being updated properly.

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


Testing with a 2x2 matrix

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


yields the following output:

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


The output I expected is array([[ 3, -2, 1],[0. 5.66666667 0.66666667]]). Why did the second row not update properly?

Continue reading...