Python Program to Implement Elementary Row Operations Assignment Solution

Instructions

Objective

Write a program to implement elementary row operations in python.

Requirements and Specifications

Source Code

import numpy as np from sympy import * def replace_rows(mat, i, j, alpha): mat[i,:] += alpha*mat[j,:] def exchange_rows(mat, i, j): temp = mat[j,:] # save the row j in a temporary variable mat[j,:] = mat[i,:] mat[i,:] = temp def scale_rows(mat, i, alpha): mat[i,:] = mat[i,:]*alpha def gauss_elimination(mat): """ Performs gauss elimination using pivoting to get the row-echelon form of the matrix :param mat: 2-D array representing the matrix :return: row echelon form of the matrix """ # getting number of rows and columns rows = len(mat) cols = len(mat[0]) visited_pivot = 0 # visited pivot keeps track of the pivot rows that have been used # selecting the columns first for col in range(cols): pivot = visited_pivot # set the current pivot to visited pivot # get the first non-zero element of the column while pivot < rows and mat[pivot][col] == 0: pivot += 1 if pivot < rows: # found a non-zero element # exchange the pivot row with the top most row which has not been pivoted exchange_rows(mat, pivot, visited_pivot) pivot = visited_pivot # scale the rows so that the element at the current column of the # pivot row becomes 1 scale_rows(mat, pivot, 1 / mat[pivot][col]) # now for all the rows below the visited pivot row # in order to make all elements at col below the pivot 0 # we perform row-replacement for i in range(pivot + 1, rows): col_val = mat[i][col] replace_rows(mat, i, pivot, -col_val) # increment visited pivot by 1 visited_pivot += 1 # return the matrix return mat def print_matrix(mat): """ Prints the matrix :param mat: 2-D array representing the matrix :return: None """ for i in range(len(mat)): for j in range(len(mat[i])): print(mat[i][j], end=' ') print('\n') def solve_linear_equations(A, B): """ Solve linear equations for A and B such that AX=B """ X = np.linalg.solve(A,B) return X def find_pivot(A): """ Pivots matrix A - finds row with maximum first entry and if nessecary swaps it with the first row. Input Arguments --------------- Augmented Matrix A Returns ------- Pivoted Augmented Matrix A """ B = np.zeros((1,2)) B[0,0]=A.shape[0] B[0,1]=A.shape[1] nrows =B[0,0] #This stores dimensions of the ncols =B[0,1] #matrix in an array pivot_size = np.abs(A[0,0]) #checks for largest first entry and swaps pivot_row = 0; for i in range(int(0),int(nrows)): if np.abs(A[i,0])>pivot_size: pivot_size=np.abs(A[i,0]) pivot_row = i if pivot_row>0: tmp = np.empty(int(ncols)) tmp[:] = A[0,:] A[0,:] = A[pivot_row,:] A[pivot_row,:] = tmp[:] return A def backsub(A): """ backsub(A) solves the upper triangular system Input Argument --------------- Augmented Matrix A Returns ------- vector b, solution to the linear system """ B = np.zeros((1,2)) B[0,0]=A.shape[0] B[0,1]=A.shape[1] n =B[0,0] ncols =B[0,1] n=int(n) x=np.zeros((n,1)) x[n-1]=A[n-1,n]/A[n-1,n-1] for i in range(int(n-1),int(-1),int(-1)): for j in range(int(i+1),int(n)): A[i,n] = A[i,n]-A[i,j]*x[j] x[i] = A[i,n]/A[i,i] return x def elim(A): """ elim(A)uses row operations to introduce zeros below the diagonal in first column of matrix A Input Argument --------------- Augmented Matrix A Returns ------- A with eliminated first column """ A = find_pivot(A) B = np.zeros((1,2)) B[0,0]=A.shape[0] B[0,1]=A.shape[1] nrows =B[0,0] ncols =B[0,1] #row operations if(float(A[0][0])): rpiv = 1./float(A[0][0]) else: rpiv = 1.0 for irow in range(int(1),int(nrows)): s=A[irow,0]*rpiv for jcol in range(int(0),int(ncols)): A[irow,jcol] = A[irow,jcol] - s*A[0,jcol] return A def gaussfe(A): """ gaussfe(A)uses row operations to reduce A to upper triangular form by calling elim and pivot Input Argument --------------- Augmented Matrix A Returns ------- A in upper triangular form """ B = np.zeros((1,2)) B[0,0]=A.shape[0] B[0,1]=A.shape[1] nrows =B[0,0] ncols =B[0,1] for i in range(int(0),int(nrows-1)): A[i:int(nrows),i:int(ncols)]=find_pivot(np.array(A[i:int(nrows),i:int(ncols)])) A[i:int(nrows),i:int(ncols)]=elim(A[i:int(nrows),i:int(ncols)]) return A def solve(A,b): """ Solve augments the nxn matrix A and column vector b then calls upon the functions in this module to solve the linear system Input Argument --------------- nxn matrix A and column vector b Returns ------- x, the solution to the linear system """ B = np.zeros((1,2)) B[0,0]=A.shape[0] B[0,1]=A.shape[1] nrows =B[0,0] ncols =B[0,1] Aug= np.zeros((int(nrows),int(ncols+1))) #these 2 loops augment the matrix with column vector for i in range(int(0),int(nrows)): for j in range(int(0),int(ncols)): Aug[i,j] = A[i,j] for k in range(int(0),int(nrows)): Aug[k,int(ncols)]= b[k] A = Aug A=gaussfe(A) x=backsub(A) x=x.T return x if __name__ == '__main__': print("Original matrix") # A = [[1, 4, 2], [0, 1, -4], [2, 7, 9]] A = np.array([[-1, 1, 2, -8, 16, 30], [4,-4,-8,28,-60,-108], [1,-1,-2,0,-12,-10], [4, -4,-8,24,-60,-100]]) print_matrix(A) print() M = Matrix([[-1, 1, 2, -8, 16, 30], [4,-4,-8,28,-60,-108], [1,-1,-2,0,-12,-10], [4, -4,-8,24,-60,-100]]) M_rref = M.rref() print("The Row echelon form of matrix M and the pivot columns : {}".format(M_rref)) print("The Row-echelon form:") B = gauss_elimination(np.array(M_rref[0])) print_matrix(B) print("Pivoting matrix") print(find_pivot(np.array(A))) print() print("Using gaussian with partial pivot") print(gaussfe(elim(find_pivot(np.array(A))))) print() print("Using gaussian with full pivot") print(gauss_elimination(np.array(A))) print() a = np.array([[-1, 1, 2, -8, 16, 30], [4,-4,-8,28,-60,-108], [1,-1,-2,0,-12,-10], [4, -4,-8,24,-60,-100]],float) #the b matrix constant terms of the equations b = np.array([-89, 328, 49, 316],float) print("Solve system of linear equations") print(solve(a,b)) # Testing 'solve(A,b)' when A is square A = np.array([[3, -1, 4], [17, 2, 1], [1, 12, -77]], float) b = np.array([2, 14, 54], float).T # Solve by using normal matrix multiplication x_exact = np.matmul(np.linalg.inv(A), b) print("Real sol: ", x_exact) # Using solve x = solve(A, b) print("Using solve: ", x)