path: root/tutorials/module_3/1_3_LU_Decomposition.md

# LU Decomposition
Imagine you’re designing the heat shield of a spacecraft. Engineers provide you with a structural model of the heat shield, represented by a coefficient matrix $A$. To evaluate how the design performs under different operating conditions, you need to test many different load cases, each represented by a right-hand side vector $\mathbf{b}$.

Now suppose your model has a 50x50 matrix $A$, and you want to simulate 100 different load cases. If you rely on Gaussian elimination alone, you would need to **repeat the elimination process 100 times** - once for each $\mathbf{b}$. This quickly becomes computationally expensive, especially for larger systems.

This is where **LU decomposition** comes in. Instead of redoing elimination for every load case, we can factorize the matrix $A$ **just once** into two simpler pieces:
$$
A = LU
$$
expanding it:
$$
\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{bmatrix} 
=
\begin{bmatrix} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \\ \end{bmatrix} 
\begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \\ \end{bmatrix} 
$$
Once we have this factorization, solving $A\mathbf{x} = \mathbf{b}$ becomes a two-step process:

1. Solve $L\mathbf{d} = \mathbf{b}$ by **forward substitution**.
$$
d_m \;=\; \frac{\,b_m \;-\; \sum_{n=1}^{m-1} a_{mn}\, d_n\,}{a_{mm}}
$$
2. Solve $U\mathbf{x} = \mathbf{d}$ by **backward substitution**.
 $$
 x_m=\frac{\,d_m-\sum_{n=m+1}^{N} u_{mn}\,x_n\,}{u_{mm}}
 $$

The big advantage is that the factorization is reused. For 100 different load vectors, we only do the expensive decomposition once, and then solve two triangular systems (fast operations) for each new $\mathbf{b}$.
<img 
    style="display: block; 
           margin-left: auto;
           margin-right: auto;
           width: 40%;"
    src="https://upload.wikimedia.org/wikipedia/commons/thumb/b/bd/Tadeusz_Banachiewicz_%28NAC%29.jpg/1024px-Tadeusz_Banachiewicz_%28NAC%29.jpg" 
    alt="Tadeusz Banchiewicz">
First formalized in 1938 by Tadeusz Banachiewicz, based on Gauss’s elimination ideas.

Example


<img 
    style="display: block; 
           margin-left: auto;
           margin-right: auto;
           width: 75%;"
    src="lu_decomp.png" 
    alt="LU Decomposition">
## Problem 1
```python
"""
LU Decomposition Solver (with partial pivoting)

Solves A x = b by computing P, L, U such that P A = L U,
then performing:
  1) y from L y = P b  (forward substitution)
  2) x from U x = y    (back substitution)

- Works for a single RHS (b shape: (n,)) or multiple RHS (b shape: (n, m))
- Raises a ValueError if a near-zero pivot is encountered.
"""

import numpy as np

def lu_decompose(A, pivot_tol=1e-12):
    """
    Compute LU factorization with partial pivoting.
    Returns P, L, U such that P @ A = L @ U.

    Parameters
    ----------
    A : (n, n) ndarray
    pivot_tol : float
        If the absolute value of the pivot falls below this, the matrix is
        treated as singular.

    Returns
    -------
    P, L, U : ndarrays of shape (n, n)
    """
    A = np.array(A, dtype=float, copy=True)
    n = A.shape[0]
    if A.shape[0] != A.shape[1]:
        raise ValueError("A must be square")

    U = A.copy()
    L = np.eye(n)
    P = np.eye(n)

    for k in range(n - 1):
        # Partial pivoting: find index of largest |U[i,k]| for i >= k
        pivot_idx = np.argmax(np.abs(U[k:, k])) + k
        pivot_val = U[pivot_idx, k]

        if abs(pivot_val) < pivot_tol:
            raise ValueError(f"Numerically singular matrix: pivot ~ 0 at column {k}")

        # Swap rows in U and P (and the part of L constructed so far)
        if pivot_idx != k:
            U[[k, pivot_idx], :] = U[[pivot_idx, k], :]
            P[[k, pivot_idx], :] = P[[pivot_idx, k], :]
            if k > 0:
                L[[k, pivot_idx], :k] = L[[pivot_idx, k], :k]

        # Gaussian elimination to build L and zero out below the pivot in U
        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] -= L[i, k] * U[k, k:]
            U[i, k] = 0.0

    # Final pivot check on U[-1, -1]
    if abs(U[-1, -1]) < pivot_tol:
        raise ValueError(f"Numerically singular matrix: pivot ~ 0 at last diagonal")

    return P, L, U

def forward_substitution(L, b):
    """
    Solve L y = b for y, where L is unit lower-triangular.
    b can be (n,) or (n, m).
    """
    L = np.asarray(L)
    b = np.asarray(b, dtype=float)
    n = L.shape[0]

    if b.ndim == 1:
        b = b.reshape(n, 1)

    y = np.zeros_like(b, dtype=float)
    for i in range(n):
        # L has ones on the diagonal (Doolittle), so we don't divide by L[i,i]
        y[i] = b[i] - L[i, :i] @ y[:i]
    return y if y.shape[1] > 1 else y.ravel()

def back_substitution(U, y):
    """
    Solve U x = y for x, where U is upper-triangular.
    y can be (n,) or (n, m).
    """
    U = np.asarray(U)
    y = np.asarray(y, dtype=float)
    n = U.shape[0]

    if y.ndim == 1:
        y = y.reshape(n, 1)

    x = np.zeros_like(y, dtype=float)
    for i in range(n - 1, -1, -1):
        if abs(U[i, i]) < 1e-15:
            raise ValueError(f"Zero (or tiny) pivot encountered at U[{i},{i}]")
        x[i] = (y[i] - U[i, i + 1:] @ x[i + 1:]) / U[i, i]
    return x if x.shape[1] > 1 else x.ravel()

def lu_solve(A, b, pivot_tol=1e-12):
    """
    Solve A x = b using LU with partial pivoting.
    Returns x and the factors (P, L, U).
    """
    P, L, U = lu_decompose(A, pivot_tol=pivot_tol)
    Pb = P @ np.asarray(b, dtype=float)
    y = forward_substitution(L, Pb)
    x = back_substitution(U, y)
    return x, (P, L, U)

def main():
    # Example system:
    #  3x + 2y -  z =  1
    #  2x - 2y + 4z = -2
    # -x + 0.5y -  z =  0
    A = np.array([
        [ 3.0,  2.0, -1.0],
        [ 2.0, -2.0,  4.0],
        [-1.0,  0.5, -1.0]
    ])
    b = np.array([1.0, -2.0, 0.0])

    x, (P, L, U) = lu_solve(A, b)
    print("Solution x:\n", x)
    print("\nP:\n", P)
    print("\nL:\n", L)
    print("\nU:\n", U)

    # Residual check
    r = A @ x - b
    print("\nResidual A@x - b:\n", r)
    print("Residual norm:", np.linalg.norm(r, ord=np.inf))

if __name__ == "__main__":
    main()
```

## Problem 2
Ohm's and Kirchhoff's laws have been used to find a set of equations to model a circuit. Using the built-in Numpy solvers. Plot the function of current at different voltages.


```python
import numpy as np
import matplotlib.pyplot as plt

# Unknown ordering
VR1, VR2, VR3, VR4, VR5, I1, I2, IM1, IM2, IM3 = range(10)

def build_system_matrices(R1, R2, R3, R4, R5):
    A = np.zeros((10,10), float)
    b_const = np.zeros(10, float)

    # KVL/KCL + Ohm's law equations
    A[0, VR1] = -1; A[0, VR2] = -1            # -VR1 - VR2 = -Vs + 5
    A[1, VR2] =  1; A[1, VR3] = -1; A[1, VR4] = -1; b_const[1] = -2  # VR2 - VR3 - VR4 = -2
    A[2, VR4] =  1; A[2, VR5] = -1; b_const[2] = 2                   # VR4 - VR5 = 2

    A[3, VR1] =  1; A[3, I1 ] = -R1           # VR1 = I1 R1
    A[4, VR2] =  1; A[4, IM1] = -R2           # VR2 = IM1 R2
    A[5, VR3] =  1; A[5, I2 ] = -R3           # VR3 = I2 R3
    A[6, VR4] =  1; A[6, IM2] = -R4           # VR4 = IM2 R4
    A[7, VR5] =  1; A[7, IM3] = -R5           # VR5 = IM3 R5

    A[8, I1 ] =  1; A[8, IM1] = -1; A[8, I2 ] = -1  # I1 = IM1 + I2
    A[9, I2 ] =  1; A[9, IM2] = -1; A[9, IM3] = -1  # I2 = IM2 + IM3
    return A, b_const

def b_vector(Vs, b_const):
    b = b_const.copy()
    b[0] = -Vs + 5.0
    return b

# Example component values
R1, R2, R3, R4, R5 = 10.0, 5.0, 8.0, 6.0, 4.0

A, b_const = build_system_matrices(R1, R2, R3, R4, R5)

# Sanity: ensure A is non-singular
condA = np.linalg.cond(A)
condA

Vs_vals = np.arange(0.0, 48.0 + 0.1, 0.1)

# Option 1 (straightforward): loop with numpy.linalg.solve (LAPACK under the hood)
IM1_vals = np.empty_like(Vs_vals)
IM2_vals = np.empty_like(Vs_vals)
IM3_vals = np.empty_like(Vs_vals)

for k, Vs in enumerate(Vs_vals):
    b = b_vector(Vs, b_const)
    x = np.linalg.solve(A, b)  # Uses LU factorization internally
    IM1_vals[k], IM2_vals[k], IM3_vals[k] = x[IM1], x[IM2], x[IM3]

# Plot
plt.figure()
plt.plot(Vs_vals, IM1_vals, label="IM1 (solve)")
plt.plot(Vs_vals, IM2_vals, label="IM2 (solve)")
plt.plot(Vs_vals, IM3_vals, label="IM3 (solve)")
plt.xlabel("Supply Voltage Vs (V)")
plt.ylabel("Motor Current (A)")
plt.title(f"Motor Currents vs Supply Voltage  (cond(A) ≈ {condA:.1f}, diff≈{max_diff:.2e})")
plt.legend()
plt.show()

```