TP 2- Conditionnement et Décomposition PA=LU

Analyse Numérique Matricielle

Polytech SU

import numpy as np
from scipy.linalg import hilbert
import matplotlib.pyplot as plt

def tri_inf(matrice, vecteur):
    n = len(vecteur)
    x = np.zeros(n)

    for i in range(n):
        x[i] = (vecteur[i] - np.dot(matrice[i][:i], x[:i])) / matrice[i][i]

    return x

def tri_sup(matrice, vecteur):
    n = len(vecteur)
    x = np.zeros(n)

    for i in range(n-1, -1, -1):
        x[i] = (vecteur[i] - np.dot(matrice[i][i+1:], x[i+1:])) / matrice[i][i]

    return x

Exercice 1

Soit

$$ A=\begin{pmatrix} \epsilon &1\\1&10 \end{pmatrix}. $$

1. Calculer à la main les deux décompositions $PA=LU$ de $A$.

1. On a $$ P = \begin{pmatrix} 0&1\\1&0 \end{pmatrix} \quad L=\begin{pmatrix} 1&0\\\epsilon &1 \end{pmatrix} \quad U=\begin{pmatrix} 1&10\\0&1- 10\epsilon \end{pmatrix} $$ et $$ P_2 = \begin{pmatrix} 1&0 \\0 & 1 \end{pmatrix} \quad L_2=\begin{pmatrix} 1&0\\ 1/ \epsilon &1 \end{pmatrix} \quad U=\begin{pmatrix} \epsilon&1\\0& 10 -1/\epsilon \end{pmatrix} $$

1. Dans chacun des cas, pour $\epsilon = 10^{-2},10^{-3},10^{-4},10^{-5}$, calculer les conditionnements de U et de L. Les réprésenter sur un graphique. On pourra utiliser la fonction cond de linalg.

epsilons=[10**(-2),10**(-3),10**(-4),10**(-5)]
condL_LU=[]
condL_PALU=[]
condU_LU=[]
condU_PALU=[]
for epsilon in epsilons: 
    L2=np.array([[1,0],[1/epsilon,1]])
    L=np.array([[1,0],[epsilon,1]])
    U=np.array([[1,10],[0,1-10*epsilon]])
    U2=np.array([[epsilon,1],[0,10-1/epsilon]])
    condU_PALU.append(np.linalg.cond(U, p=2))
    condU_LU.append(np.linalg.cond(U2, p=2))
    condL_PALU.append(np.linalg.cond(L, p=2))
    condL_LU.append(np.linalg.cond(L2, p=2))

# Tracer les graphiques
plt.figure(figsize=(10, 6))

plt.loglog(epsilons, condL_LU, label='condL_LU')
plt.loglog(epsilons, condL_PALU, label='condL_PALU')
plt.loglog(epsilons, condU_LU, label='condU_LU')
plt.loglog(epsilons, condU_PALU, label='condU_PALU')

plt.title('Conditions en fonction de epsilon')
plt.xlabel('Epsilon')
plt.ylabel('Valeur de conditionnement')
plt.legend()
plt.grid(True)
plt.show()

1. On prend maintenant $\epsilon = 10^{-3}$. Utiliser la fonction lude linalg. Que constatez-vous?

# 3. Utiliser scipy pour obtenir la décomposition LU
import numpy as np
from scipy.linalg import lu
epsilon=10**(-3)
A = np.array([[epsilon, 1], [1, 10]])

P, L, U = lu(A)

print("P:\n", P)
print("L:\n", L)
print("U:\n", U)

P:
 [[0. 1.]
 [1. 0.]]
L:
 [[1.    0.   ]
 [0.001 1.   ]]
U:
 [[ 1.   10.  ]
 [ 0.    0.99]]

1. On considère maintenant $b = \begin{pmatrix}2 \\2 \end{pmatrix}$. Résoudre $Ax = b$ à l'aide des deux décompositions précédentes. On utilisera les fonctions tri_inf et tri_sup. Comparer avec la vraie solution $x_{true} = \frac{1}{10\epsilon -1} \begin{pmatrix} 18 \\ 2\epsilon -2\end{pmatrix}$ pour $\epsilon = 10^{-2},10^{-3},10^{-4},10^{-5}$. On tracera les erreurs sur un même graphique.

### résolution
errors_LU=[]
errors_PALU=[]

epsilon=0.001

for epsilon in epsilons:
    L2=np.array([[1,0],[1/epsilon,1]])
    U2=np.array([[epsilon,1],[0,10-1/epsilon]])
    A = np.array([[epsilon, 1], [1, 10]])
    P, L, U = lu(A)
    x_true= 1/(10*epsilon -1)*np.array([18 ,2*epsilon-2])
    b=np.array([2,2])
    y = tri_inf(L,P@b)
    x = tri_sup(U,y)
    y2 = tri_inf(L2,b)
    x2 = tri_sup(U2,y2)
    errors_PALU.append(np.linalg.norm(x-x_true))
    errors_LU.append(np.linalg.norm(x2-x_true))

print(errors_PALU) 
plt.figure(figsize=(10, 6))

plt.loglog(epsilons, errors_LU, label='errors_LU')
plt.loglog(epsilons, errors_PALU, label='errors_PALU')

plt.title('Erreurs en fonction de epsilon')
plt.xlabel('Epsilon')
plt.ylabel('Erreurs')
plt.legend()
plt.grid(True)
plt.show()

[3.580361673049448e-15, 3.580361673049448e-15, 0.0, 4.440892098500626e-16]

Exercice 2

1. Implémenter l'algorithme de la décomposition $PA=LU$ en Python. Calculer la décomposition de la matrice A obenue grâce à la fonction create_matrix.
2. Vérifier que l'on trouve bien les mêmes résultats que la fonction scipy.linalg.lu.
3. Implémenter l'algorithme de recherche de solution de l'équation $Ax=b$ à partir de la décomposition LU de $A$. Résolvez l'équation $Ax=b$ avec le vecteur $b$ donné en 4.(b) de l'exercice $2$.
4. En déduire une fonction pour obetnir l'inverse dune matrice $A$.

# Fonctrion create_matrix

def create_matrix(n):
    A = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
    return A

# 1. Implémentation de l'algorithme PA=LU
def lu_decomposition(A):
    n = A.shape[0]
    P = np.eye(n)
    L = np.eye(n)
    U = A.copy()

    for k in range(n - 1):
        # Échange de lignes pour le pivot
        max_row = np.argmax(np.abs(U[k:, k])) + k
        U[[k, max_row], :] = U[[max_row, k], :]
        P[[k, max_row], :] = P[[max_row, k], :]
        L[[k, max_row], :k] = L[[max_row, k], :k]

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

    return P, L, U

# Tester la décomposition sur la matrice de 4.(a) de l'exercice 2
A_test = create_matrix(5)  # À partir du code précédent
P_test, L_test, U_test = lu_decomposition(A_test)
print("P :\n", P_test)
print("L  :\n", L_test)
print("U  :\n", U_test)

P :
 [[1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]]
L  :
 [[ 1.          0.          0.          0.          0.        ]
 [-0.5         1.          0.          0.          0.        ]
 [ 0.         -0.66666667  1.          0.          0.        ]
 [ 0.          0.         -0.75        1.          0.        ]
 [ 0.          0.          0.         -0.8         1.        ]]
U  :
 [[ 2.         -1.          0.          0.          0.        ]
 [ 0.          1.5        -1.          0.          0.        ]
 [ 0.          0.          1.33333333 -1.          0.        ]
 [ 0.          0.          0.          1.25       -1.        ]
 [ 0.          0.          0.          0.          1.2       ]]

array([[ 2., -1.,  0.,  0.,  0.],
       [-1.,  2., -1.,  0.,  0.],
       [ 0., -1.,  2., -1.,  0.],
       [ 0.,  0., -1.,  2., -1.],
       [ 0.,  0.,  0., -1.,  2.]])

#2. Vérification
import numpy as np
from scipy.linalg import lu

P, L, U = lu(A_test)

print("Ptest:\n", P)
print("Ltest:\n", L)
print("Utest:\n", U)

Ptest:
 [[1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]]
Ltest:
 [[ 1.          0.          0.          0.          0.        ]
 [-0.5         1.          0.          0.          0.        ]
 [ 0.         -0.66666667  1.          0.          0.        ]
 [ 0.          0.         -0.75        1.          0.        ]
 [ 0.          0.          0.         -0.8         1.        ]]
Utest:
 [[ 2.         -1.          0.          0.          0.        ]
 [ 0.          1.5        -1.          0.          0.        ]
 [ 0.          0.          1.33333333 -1.          0.        ]
 [ 0.          0.          0.          1.25       -1.        ]
 [ 0.          0.          0.          0.          1.2       ]]

import numpy as np
from scipy.linalg import lu



# 3. Implémentation de la recherche de solution Ax=b à partir de la décomposition LU
def lu_solve(P, L, U, b):
    # Étape 1 : appliquer la permutation
    Pb = P @ b

    # Étape 2 : résoudre L y = Pb (descente)
    y = tri_inf(L, Pb)

    # Étape 3 : résoudre U x = y (remontée)
    x = tri_sup(U, y)

    return x

# Résoudre Ax=b avec b donné en 4.(b) de l'exercice 2
b_test = np.full(5, 2)
x_test = lu_solve(P_test, L_test, U_test, b_test)
print("Solution x:\n", x_test)

Solution x:
 [5. 8. 9. 8. 5.]

# 4. Fonction inverse

def inverse_lu(A):
    n = A.shape[0]
    A_inv = np.zeros((n, n))
    
    # On recalcule complètement PA=LU
    P, L, U = lu_decomposition(A)
    
    # On résout Ax_i = e_i
    for i in range(n):
        e_i = np.zeros(n)
        e_i[i] = 1
        A_inv[:, i] = lu_solve(P, L, U, e_i)
    
    return A_inv


A_test=create_matrix(5)

A_inv = inverse_lu(A_test)
print(A_inv @ A_test)

[[ 1.00000000e+00 -1.66533454e-16  5.55111512e-17  2.77555756e-17
   0.00000000e+00]
 [ 0.00000000e+00  1.00000000e+00  1.11022302e-16  5.55111512e-17
   0.00000000e+00]
 [ 0.00000000e+00 -2.22044605e-16  1.00000000e+00  2.22044605e-16
   0.00000000e+00]
 [ 0.00000000e+00 -2.22044605e-16  0.00000000e+00  1.00000000e+00
   0.00000000e+00]
 [ 0.00000000e+00 -1.11022302e-16  0.00000000e+00  1.11022302e-16
   1.00000000e+00]]

Exercice 3

On souhaite comparer la précision de deux méthodes d'inversion matricielle :

• inverse_lu, basée sur la décomposition LU
• numpy.linalg.inv, l'inversion classique fournie par NumPy

Pour cela, on génère des matrices aléatoires à l’aide de la fonction create_matrix pour différentes tailles :

n ∈ {5, 10, 20, 100}

1. Pour chaque matrice A, on calcule une approximation de son inverse notée Â⁻¹ à l’aide des deux méthodes, puis on calcule l’erreur relative définie par :

    || A Â⁻¹ - Iₙ || / || A ||

où || · || désigne la norme de Frobenius.

1. Pour chaque valeur de n, mesurer également le temps d'exécution des deux méthodes d’inversion. On pourra utiliser a fonction time.perf_counter.

#1. Calcul des erreurs

def relative_error(A, A_inv):
    I = np.eye(A.shape[0])
    return np.linalg.norm(A @ A_inv - I, ord='fro') / np.linalg.norm(A, ord='fro')


sizes = [5, 10, 20, 100,100]
errors_lu = []
errors_np = []

for n in sizes:
    A = create_matrix(n)

    A_inv_lu = inverse_lu(A)
    A_inv_np = np.linalg.inv(A)

    errors_lu.append(relative_error(A, A_inv_lu))
    errors_np.append(relative_error(A, A_inv_np))


    plt.figure(figsize=(8, 5))
plt.plot(sizes, errors_lu, marker='o', label='Inverse via LU')
plt.plot(sizes, errors_np, marker='s', label='numpy.linalg.inv')

plt.xscale('log')
plt.yscale('log')
plt.xlabel("Dimension de la matrice (n)")
plt.ylabel("Erreur relative (log scale)")
plt.title("Comparaison des erreurs d'inversion matricielle")
plt.legend()
plt.grid(True, which="both")
plt.show()

<Figure size 800x500 with 0 Axes>

<Figure size 800x500 with 0 Axes>

<Figure size 800x500 with 0 Axes>

<Figure size 800x500 with 0 Axes>

# 2. Calcul du temps

import time   

sizes = [5, 10, 20, 100]
times_lu = []
times_np = []

for n in sizes:
    A = create_matrix(n)

    # Mesure du temps LU
    start = time.perf_counter()
    _ = inverse_lu(A)
    times_lu.append(time.perf_counter() - start)

    # Mesure du temps numpy.inv
    start = time.perf_counter()
    _ = np.linalg.inv(A)
    times_np.append(time.perf_counter() - start)

# ---- Visualisation ----
plt.figure(figsize=(8, 5))
plt.plot(sizes, times_lu, marker='o', label='Inverse via LU')
plt.plot(sizes, times_np, marker='s', label='numpy.linalg.inv')

plt.xscale('log')
plt.yscale('log')
plt.xlabel("Taille de la matrice (n)")
plt.ylabel("Temps d'exécution (secondes)")
plt.title("Comparaison des temps d'inversion matricielle")
plt.grid(True, which="both")
plt.legend()
plt.show()