Triangular System Solvers

Triangular System Solvers

After LU decomposition, solving $Ax = b$ reduces to two cheap triangular solves:

$A = LU \ \Rightarrow \ LUx = b$

• Step 1: $Ly = b$ (forward substitution)
• Step 2: $Ux = y$ (backward substitution)

Forward Substitution ($Ly = b$)

$L$ is lower triangular, so solve top-to-bottom:

$y_{0} = \frac{b_{0}}{L_{00}}, \qquad y_{i} = \frac{b_{i} - \sum_{j < i} L_{ij}\, y_{j}}{L_{ii}}$

Backward Substitution ($Ux = y$)

$U$ is upper triangular, so solve bottom-to-top:

$x_{n-1} = \frac{y_{n-1}}{U_{n-1,n-1}}, \qquad x_{i} = \frac{y_{i} - \sum_{j > i} U_{ij}\, x_{j}}{U_{ii}}$

Example

$L = \begin{pmatrix} 1&0&0\\2&1&0\\4&3&1 \end{pmatrix}, \quad U = \begin{pmatrix} 2&1&1\\0&1&1\\0&0&2 \end{pmatrix}, \quad b = \begin{pmatrix} 1\\3\\9 \end{pmatrix}$

$Ly = b \;\ \Rightarrow \;\ y = \begin{pmatrix}1.0\\1.0\\2.0\end{pmatrix}, \qquad Ux = y \;\ \Rightarrow \;\ x = \begin{pmatrix}0.0\\0.0\\1.0\end{pmatrix}$

Your Task

Implement forward_sub(L, b) and backward_sub(U, b), then combine them in solve_lu(L, U, b).