# Code for testing the steps of the Gauss-Jordan method # for finding the solution 'x' to a linear system A x = b, # along with the inverse matrix A^{-1}. import numpy as np np.set_printoptions(threshold=np.inf) # print the whole augmented array def swaprows(M,i,j): M[i-1,:], M[j-1,:] = M[j-1,:].copy(), M[i-1,:].copy() def augmented(A,b): n = A.shape[0] aug = np.zeros((n,2*n+1)) aug[:,:n] = A aug[:,n:2*n] = np.identity(n) aug[:,2*n] = b.transpose() return aug def GaussJordan(A,b): """The Gauss-Jordan method for solving linear equations. Returns the inverse matrix of A, and the solution x to Ax=b.""" n = A.shape[0] M = augmented(A,b) # forward elimination with pivoting for i in range(n): B = np.abs(M[:,:n]) argmax = np.argmax(B,axis=1) # the index of the maximum absolute value in each column j = argmax[i] if j > i: swaprows(M, i+1, j+1) M[i,:] = M[i,:] / M[i,i] for j in range(i+1,n): M[j,:] = M[j,:] - M[i,:]*M[j,i] print(M,"\n") # back substitution for i in list(range(n))[::-1]: for j in range(i): M[j,:] = M[j,:] - M[j,i]*M[i,:] print(M,"\n") return M[:,n:2*n], M[:,2*n] A = np.matrix("0,2,1;2,-1,1;1,3,2", dtype=np.float64) b = np.matrix("7;3;13", dtype=np.float64) #A = np.matrix("0,2,3;0,1,2;1,0,5", dtype=np.float64) #b = np.matrix("1;1;7", dtype=np.float64) GaussJordan(A,b) ### The final augmented matrix should look like this: #### sol = np.array(np.matrix("1,0,0,-5,-1,3,1;0,1,0,-3,-1,2,2;0,0,1,7,2,-4,3"), dtype=np.float64) #sol = np.array(np.matrix("1,0,0,5,-10,1,2;0,1,0,2,-3,0,-1;0,0,1,-1,2,0,1"), dtype=np.float64) print("-------") print(sol,"\n")