Untitled

import numpy as np

def GaussElim(A, b):
    """create a Gaussian elimination matrix for a system
    Args:
        A: N by N array
        b: array of length N
    Returns:
        augmented matrix ready for back substitution
    """

     # making sure input is a square matrix
    [Nrow, Ncol] = A.shape
    N = Nrow
    aug_matrix = np.zeros((N, N + 1))
    aug_matrix[0:N, 0:N] = A
    aug_matrix[:, N] = b

    for column in range(0, N):                                     # designates pivot and starts iterating thru rows
        for row in range(column + 1, N):
            mult = aug_matrix[row, column] / aug_matrix[column, column]   # calculating multiplier to eliminate a variable
            aug_matrix[row, :] -= mult * aug_matrix[column, :]            # executes the elimination
    return aug_matrix


def BackSub(aug_matrix):
    N = aug_matrix.shape[0]                 # how many equations are in system
    x = np.zeros(N)                         # creates empty array w/ # of equations
    for row in range(N - 1, -1, -1):        # starting at bottom row to work our way up. last row has only 1 unknown
        RHS = aug_matrix[row, N]            # already solved for value
        for column in range(row + 1, N):
            RHS -= x[column] * aug_matrix[row, column]         # moves everything to one side of equation to solve for unknown
        x[row] = RHS / aug_matrix[row, row]                    # divides isolated term by coefficient
    return x

def fill_in(n):
    A = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            A[i, j] = 1.0 / (i + j + 1)           # creating the n x n matrix and filling in by 1/i+j+1 for any i/j

    b = np.ones(n)                                # creating vector of size n of all ones
    b[1::2] = -1.0                                # making every other one -1, starting w/ 1

    aug = GaussElim(A, b)            # calling function to initiate Gaussian Elimination
    x = BackSub(aug)                 # calling function to solve through back substitution, solves for x

    Ax = np.dot(A, x)                # dots our x with A, to get the expression Ax (LHS)
    res_vector = Ax - b              # calculates the error of our function bc Ax is supposed to equal b
    residual = np.linalg.norm(res_vector)      # gets the residual in a magnitude (number) form
    return x, residual

for n in [2, 4, 8, 16]:                 # for all of our desired n's
    sol, res = fill_in(n)               # sol is our x, res is the residual
    print(f"n = {n}")
    print(f"sol = {sol}")
    print(f"Residual: {res:.2e}\n")