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