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. 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.

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

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.

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 AP=LU

#2. Vérification

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

# 4. Fonction Inverse

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 la fonction time.perf_counter.

# 1. Calcul de l'erreur

# 2. Calcul du temps

import time