Décomposition LU sans pivot

import numpy as np

def lu_decomposition(A):
    """
    Décompose A = LU où L est triangulaire inférieure avec des 1 sur la diagonale.
    A doit être inversible.
    """
    n = A.shape[0]
    L = np.eye(n)
    U = A.astype(float).copy()

    for j in range(n - 1):
        if abs(U[j, j]) < 1e-15:
            raise ValueError(f"Pivot nul en position ({j}, {j})")

        for i in range(j + 1, n):
            factor = U[i, j] / U[j, j]
            L[i, j] = factor
            U[i, j:] -= factor * U[j, j:]

    return L, U

# Test
A = np.array([[4, 3, 2],
              [2, 2, 1],
              [3, 1, 2]])

L, U = lu_decomposition(A)
print("Matrice A :")
print(A)
print("\nL :")
print(L)
print("\nU :")
print(U)
print("\nL @ U :")
print(L @ U)
print(f"\nErreur de décomposition : {np.linalg.norm(A - L @ U):.2e}")

Résultat dans la console :

Matrice A :
[[4 3 2]
 [2 2 1]
 [3 1 2]]

L :
[[1.   0.   0.  ]
 [0.5  1.   0.  ]
 [0.75 1.25 1.  ]]

U :
[[4.   3.   2.  ]
 [0.   0.5  0.  ]
 [0.   0.   0.5 ]]

L @ U :
[[4. 3. 2.]
 [2. 2. 1.]
 [3. 1. 2.]]

Erreur de décomposition : 0.00e+00