Solving 3x3 systems of equations using Gaussian Elimination

solve3x3.py

'''Using Gaussian Elimination Method to solve 3x3 systems of equations.
We'll use numpy arrays to store and manipulate the equations.  For a video demonstrating the process,
I found https://youtu.be/Dj84vEb4Zko to be helpful.
'''

import numpy as np

A = np.zeros((3, 3)) #This will store the coefficients of the variables as a 3x3 numpy array
x = np.zeros((3, 1)) #This stores the values of the variables as a 3x1 numpy array
b = np.zeros((3, 1)) #This will store what the equations are equal to as a 3x1 numpy array

varNames = ["x", "y", "z"]

A = np.array([(1, -2, 1), #Use the example problem in the video, you can edit this to whatever you want.
              (2, -1, 2),
              (3, -1, 2)
              ])
b = np.reshape(np.array([8, 10, 11]), (3, 1))

'''Now that we have populated our arrays, we'll use Gaussian Elimination to solve for x (the variables) in Ax=b
Our goal is to get A formatted so it looks like: [1 _ _]
                                                 [0 1 _]
                                                 [0 0 _]
This format is called row echelon form (REF).  We'll accomplish this by eliminating portions of A.
'''

if not A[0][0] == 1:
  multiplier = 1/A[0][0] #We'll multiply all of row 1 by 1/(the number in the top left) to convert it to a one
  multiplierMatrix = np.array([(multiplier, 0, 0),
                               (0, 1, 0),
                               (0, 0, 1),
                              ])

  A = np.dot(multiplierMatrix, A)
  b = np.dot(multiplierMatrix, b)

'''Next, we use elimination to get A[1][0] to be 0'''

if not A[1][0] == 0:
  multiplier = -1 * A[1][0] / A[0][0] #Calculates how much of row 1 to add to row 2
  multiplierMatrix = np.array([(1, 0, 0),
                               (multiplier, 1, 0),
                               (0, 0, 1)
                              ])
  
  A = np.dot(multiplierMatrix, A)
  b = np.dot(multiplierMatrix, b)

'''Now, we use elimination once again to get A[2][0] to be 0'''
if not A[2][0] == 0:
  multiplier = -1 * A[2][0] / A[0][0] #Calculates how much of row 1 to add to row 3
  multiplierMatrix = np.array([(1, 0, 0),
                               (0, 1, 0),
                               (multiplier, 0, 1)
                              ])
  
  A = np.dot(multiplierMatrix, A)
  b = np.dot(multiplierMatrix, b)

'''Now, we use division to get A[1][1] to be 1'''
if not A[1][1] == 1:
  multiplier = 1 / A[1][1] #Calculates what to multiply A[1][1] by to make it equal to 1
  multiplierMatrix = np.array([(1, 0, 0),
                               (0, multiplier, 0),
                               (0, 0, 1)
                              ])
  
  A = np.dot(multiplierMatrix, A)
  b = np.dot(multiplierMatrix, b)

'''Next, we use elimination to get A[2][1] to be 1'''
if not A[2][1] == 0:
  multiplier = -1 * A[2][1] / A[1][1] #Calculates how much of row 2 to add to row 3
  multiplierMatrix = np.array([(1, 0, 0),
                               (0, 1, 0),
                               (0, multiplier, 1)
                              ])
  
  A = np.dot(multiplierMatrix, A)
  b = np.dot(multiplierMatrix, b)

'''Next, we use division to get A[2][2] to be 1'''
if not A[2][2] == 1:
  multiplier = 1 / A[2][2] #Calculates what to multiply A[2][2] by to make it equal to 1
  multiplierMatrix = np.array([(1, 0, 0),
                               (0, 1, 0),
                               (0, 0, multiplier)
                              ])
  
  A = np.dot(multiplierMatrix, A)
  b = np.dot(multiplierMatrix, b)

'''Now, that we have A in REF, we can start finding the values of variables'''
x[2][0] = b[2][0] #We know the last variable is equal to the last constant, because we've isolated it using elimination
x[1][0] = b[1][0] #Similarly, the second variable was also isolated, so it's equal to the second constant
x[0][0] = b[0][0] - A[0][1] * x[1][0] - A[0][2] * x[2][0]
#Finding first variable is a little bit trickier.  We accomplish this by substituting in the values for
#The variables we already found, and then solving for the unknown using subtraction.

'''Now, we solved for each variable, and we are done!'''
print ("=======SOlUTION======")
print (varNames[0] + ":\t" + str(x[0][0]))
print (varNames[1] + ":\t" + str(x[1][0]))
print (varNames[2] + ":\t" + str(x[2][0]))