Forward-backward substitution

ref: https://ubcmath.github.io/MATH307/systems/lu.html#forward-and-backward-substitution

Given $A=LU$ by LU factorization, we can adapt this technique to solve $Ax=b$ by $Ly=b$ and $Ux=y$

Python code example for LU factorization

np.linalg

it is equivalent with codes:

y = np.linalg.solve(L, b)
x = np.linalg.solve(U, y)

custom implementation

import numpy as np
 
def fw_sub(L, b):
    n = L.shape[0]
    y = np.zeros((n, 1))
    
    for i in range(n):
        sum_val = 0
        for j in range(i):
            sum_val += L[i, j] * y[j, 0]
        y[i, 0] = (b[i, 0] - sum_val) / L[i, i]
    
    return y
 
def bw_sub(U, y):
    n = U.shape[0]
    x = np.zeros((n, 1))
    
    for i in range(n-1, -1, -1):
        sum_val = 0
        for j in range(i+1, n):
            sum_val += U[i, j] * x[j, 0]
        x[i, 0] = (y[i, 0] - sum_val) / U[i, i]
    
    return x
 
L = np.array([
	[1, 0, 0],
	[2, 1, 0],
	[1, 1, 1]
])
U = np.array([
	[2, 4, 2],
	[0, 1, 1],
	[0, 0, 1]
])
 
# solve Ax = b
A = L @ U
b = np.array([[-1], [1], [2]])
print(f"L={L}")
print(f"U={U}")
print(f"b={b}")
# A=[[2 4 2]
# [4 9 5]
# [2 5 4]]
print(f"A={A}")
 
# solve Ly = b
y = fw_sub(L, b)
# solve Ux = y
x = bw_sub(U, y)
 
# y=[[-1.]
#  [ 3.]
#  [ 0.]]
print(f"y={y}")
# x=[[-6.5]
#  [ 3. ]
#  [ 0. ]]
print(f"x={x}")
# Ax=[[-1.]
#  [ 1.]
#  [ 2.]]
# Ax should be eq to b
print(f"Ax={A @ x}")