Programming Assignment: Gaussian Elimination

Objective

Implement Gaussian elimination with back substitution to solve \(A\mathbf{x} = \mathbf{b}\) without using numpy.linalg.solve.

import numpy as np

def gaussian_elimination(A, b):
    n = len(b)
    Ab = np.hstack([A.astype(float), b.reshape(-1, 1).astype(float)])

    # Forward elimination
    for col in range(n):
        # Partial pivoting
        max_row = np.argmax(np.abs(Ab[col:, col])) + col
        Ab[[col, max_row]] = Ab[[max_row, col]]

        if np.isclose(Ab[col, col], 0):
            raise ValueError("Matrix is singular or nearly singular")

        for row in range(col + 1, n):
            factor = Ab[row, col] / Ab[col, col]
            Ab[row] -= factor * Ab[col]

    # Back substitution
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - np.dot(Ab[i, i+1:n], x[i+1:n])) / Ab[i, i]

    return x

Test Case 1: Unique Solution

A = np.array([[2, 1, -1],
              [-3, -1, 2],
              [-2, 1, 2]])
b = np.array([8, -11, -3])

x = gaussian_elimination(A, b)
print("Solution:", x)
print("Verification Ax - b:", np.round(A @ x - b, 10))

Test Case 2: Compare with NumPy

np.random.seed(0)
A_rand = np.random.randint(1, 10, (4, 4)).astype(float)
b_rand = np.random.randint(1, 20, 4).astype(float)

x_ours = gaussian_elimination(A_rand, b_rand)
x_numpy = np.linalg.solve(A_rand, b_rand)

print("Our solution: ", np.round(x_ours, 6))
print("NumPy solution:", np.round(x_numpy, 6))
print("Match:", np.allclose(x_ours, x_numpy))