Gaussian Elimination

Gaussian Elimination#

# Math 545: uMass Amherst Spring 2024
%load_ext autoreload
%autoreload 2

import os
import numpy as np
import scipy as sp
import scipy.linalg as LA
import matplotlib.pyplot as plt

# For output
homedir = os.path.expanduser('~')

Special solving functions#

def backSolve(A):
    """
        Assumes that A is an upper triangular augmented matrix.
    """
    A = A.astype(float)
    nrows = A.shape[0]
    ncols = A.shape[1]
    x = np.empty(ncols-1)
    for i in range(x.size - 1, -1, -1):
        # Compute the sum
        sum = 0
        for j in range(i+1, ncols-1):
            sum += A[i, j] * x[j]
        
        x[i] = (A[i, ncols-1] - sum) / A[i, i]
        
    return x

Gaussian Elimination implementation#

def gaussElim(A):
    """
        Perform GE. Does overwrite elements below diagonal
    """
    nrows = A.shape[0]
    ncols = A.shape[1]
    
    for k in range(ncols-1):         # Iterate over all columns
        for i in range(k+1, nrows):  # Walk down the column
            lik = A[i, k] / A[k, k]
            
            # Do row operation
            for j in range(k, ncols): 
                A[i, j] -= lik * A[k, j]
    
    return A

# Example
A = np.array([[2, 1, 3, 1],
              [4, 4, 7, 1],
              [2, 5, 9, 3]])

A = gaussElim(A)
x = backSolve(A)

print(A)
print(x)

[[ 2  1  3  1]
 [ 0  2  1 -1]
 [ 0  0  4  4]]
[-0.5 -1.   1. ]
Contents