Solving Linear Systems

A system of linear equations — several equations, several unknowns, all linear — comes up constantly: balancing a portfolio across assets, fitting a regression line, solving for equilibrium prices. For a 2×2 system, Cramer's rule gives a direct formula using determinants (from the Linear Algebra course):

def solve_2x2(a, b, c, d, e, f):
    # ax + by = e
    # cx + dy = f
    det = a*d - b*c
    x = (e*d - b*f) / det
    y = (a*f - e*c) / det
    return [x, y]

print(solve_2x2(2, 1, 1, 3, 5, 10))   # [1.0, 3.0]

Cramer's rule doesn't scale well past 2×2 or 3×3 (the determinant calculations get expensive fast) — for larger systems, you need a different technique.

Gaussian elimination scales to any size: systematically eliminate each variable from every equation below it, producing an upper-triangular system where the LAST equation has only one unknown — solve that, then substitute backward through the rest.

def gaussian_elimination_3x3(A, b):
    A = [row[:] for row in A]
    b = b[:]
    n = 3
    for i in range(n):
        for j in range(i + 1, n):
            factor = A[j][i] / A[i][i]
            for k in range(i, n):
                A[j][k] -= factor * A[i][k]
            b[j] -= factor * b[i]
    # A is now upper-triangular -- solve backward
    x = [0] * n
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - sum(A[i][k] * x[k] for k in range(i + 1, n))) / A[i][i]
    return x

This is exactly the algorithm behind numpy.linalg.solve and every production linear-algebra library — real implementations add pivoting (swapping rows to avoid dividing by a near-zero number) for numerical stability, but the core idea is precisely this.