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?
Because you use numpy array of integers, all your result will be rounded off to integers. You need to define A and b as
A = np.array([[3,-2],[1,5]], dtype=np.float64)
b = np.array([1,1], dtype=np.float64)
Doing so allows for float values in your matrices.