TP 1- Conditionnement et Décomposition LU¶

Analyse Numérique Matricielle¶

Polytech SU¶

In [3]:
import numpy as np
from scipy.linalg import hilbert
import matplotlib.pyplot as plt

Exercice 1¶

On considère la matrice de Hilbert $ H^{[n]} \in M_n(\mathbb{R}) $ définie pour tout $n \geq 1$ par $$ H^{[n]}_{ij}=\frac{1}{i+j-1},\, 1\leq i,j\leq n. $$

  1. Construire la matrice $ H^{[5]} $ en utilisant la bibliothèque scipy et la fonction scipy.linalg.hilbert ou avec des boucles for.
  2. Calculer $ Cond(H^{[n]}) $ pour la norme $ \Vert\cdot\Vert_2 $ et pour$ 1\leq n\leq 10 $ à l'aide de la fonction numpy.linalg.cond. Visualisez le résultat en échelle semilogarithmique en $ y $ (à l'aide de la fonction matplotlib.pyplot.semilogy). La matrice $ H^{[n]} $ est-elle bien ou mal conditionnée selon vous ?
  3. Répéter le point précédent avec $ Id+H^{[n]} $ au lieu de $ H^{[n]} $. Que remarquez vous ?
In [4]:
# 1. Construire la matrice H^{[5]}
H5_scipy = hilbert(5)  # Utilisation de scipy

# Alternative avec boucles for
n = 5
H5_loops = np.zeros((n, n))
for i in range(1, n+1):
    for j in range(1, n+1):
        H5_loops[i-1][j-1] = 1 / (i + j - 1)

# Assurez-vous que les deux méthodes produisent le même résultat
assert np.allclose(H5_scipy, H5_loops), "Les deux méthodes devraient donner le même résultat"
In [5]:
# 2. Calculer Cond(H^{[n]}) pour la norme \Vert\cdot\Vert_2 et visualiser
n_values = range(1, 11)  # n de 1 à 10
condition_numbers = [np.linalg.cond(hilbert(n)) for n in n_values]

plt.figure(figsize=(10, 6))
plt.semilogy(n_values, condition_numbers, marker='o')
plt.xlabel('n')
plt.ylabel('Cond(H^{[n]})')
plt.title('Conditionnement de la matrice de Hilbert en fonction de n')
plt.grid(True)
plt.show()
In [6]:
# 3. Répéter avec Id+H^{[n]}
condition_numbers_id = [np.linalg.cond(np.eye(n) + hilbert(n)) for n in n_values]

plt.figure(figsize=(10, 6))
plt.semilogy(n_values, condition_numbers_id, marker='o', color='red')
plt.xlabel('n')
plt.ylabel('Cond(Id + H^{[n]})')
plt.title('Conditionnement de Id + Matrice de Hilbert en fonction de n')
plt.grid(True)
plt.show()

Exercice 2¶

On considère la matrice de Vandermonde $$A=\begin{pmatrix} x_1^{n-1}&\dots&x_1^2&x_1&1\\ \vdots&\vdots&x_2^2&x_2&1\\ \vdots&\vdots&\vdots&\vdots&\vdots\\ x_n^{n-1}&\dots&x_n^2&x_n&1 \end{pmatrix},$$ où les points $x_1,\dots,x_n$ sont équirépartis sur $[0,1]$. Cette matrice peut être construite en Python à l'aide des fonctions suivantes : y = np.linspace(0, 1, n) A = np.vander(y, increasing=False) On souhaite résoudre le système linéaire $Ax=b$ où $$b=\begin{pmatrix} 1+y_1^2\\\vdots\\1+y_n^2 \end{pmatrix}.$$ La solution exacte est connue et vaut $x=(0,\ldots,0,1,0,1)^T$.

  1. Résoudre le système linéaire en utilisant la fonction numpy.linalg.solve. Evaluer l'erreur relative pour $n=4$.
  2. La variable d'arrondi en Python peut être obtenue via numpy.finfo(float).eps. Comparer l'erreur relative calculée au point précédent avec la borne obtenue pour $n=4$.
  3. Répéter les points $1.$ et $2.$ pour $n=4,6,8,\ldots,20$. Visualiser l'erreur $\epsilon_n$, la borne supérieure $\eta_n$ ainsi que le résidu normalisé $r_n$ en fonction de $n$. Pour cela, tracer deux graphes, un avec une échelle logarithmique à l'aide de la fonction matplotlib.pyplot.loglog et l'autre en échelle semilogarithmique à l'aide de la fonction matplotlib.pyplot.semilogy. Quel type de croissance de l'erreur $\epsilon_n$ observe-t-on ? Est-ce que le résidu $r_n$ est un bon indicateur de l'erreur $\epsilon_n$ ? Et $\eta_n$ ?
  4. (a) Construire la matrice $$A=\begin{pmatrix} 2&-1&&&\\-1&2&\ddots&&\\&\ddots&\ddots&\ddots&\\&&\ddots&\ddots&-1\\ &&&-1&2 \end{pmatrix}.$$ pour $n\geq 1$ en utilisant la bibliothèque numpy. (b) Répéter le point $3.$ pour cette matrice avec $b=(2,2,\ldots,2)^T$ pour $n=5,10,\ldots,100$. Commenter les résultats obtenus.
In [19]:
### 1. Résoudre le système linéaire pour \( n = 4 \) et évaluer l'erreur relative :

import numpy as np

n = 4
y = np.linspace(0, 1, n)
A = np.vander(y, increasing=False)
b = 1 + y**2

x_exact = np.zeros(n)
x_exact[-3:] = [1, 0,1]

x_computed = np.linalg.solve(A, b)

# Calcul de l'erreur relative
error = np.linalg.norm(x_computed - x_exact, ord=2) / np.linalg.norm(x_exact, ord=2)
print(f"Erreur relative pour n=4: {error}")

Erreur relative pour n=4: 2.1268813444517158e-15
In [8]:
### 2. Comparer l'erreur relative avec la borne pour \( n = 4 \):


epsilon = np.finfo(float).eps
eta_4 = epsilon * np.linalg.cond(A, p=2)
print(f"Erreur d'arrondi en Python: {epsilon}")
print(f"Borne pour n=4: {eta_4}")

Erreur d'arrondi en Python: 2.220446049250313e-16
Borne pour n=4: 2.195304793666867e-14
In [23]:
### 3. Répéter pour \( n=4,6,8,\ldots,20 \) et visualiser les erreurs:

ns = list(range(4, 21, 2))
errors = []
etas = []
residuals = []

for n in ns:
    y = np.linspace(0, 1, n)
    A = np.vander(y, increasing=False)
    b = 1 + y**2

    x_exact = np.zeros(n)
    x_exact[-3:] = [1, 0,1]

    x_computed = np.linalg.solve(A, b)

    error = np.linalg.norm(x_computed - x_exact, ord=2) / np.linalg.norm(x_exact, ord=2)
    eta = epsilon * np.linalg.cond(A, p=2)
    residual = np.linalg.norm(b - np.dot(A, x_computed), ord=2) / (np.linalg.norm(np.dot(A, x_computed)))

    errors.append(error)
    etas.append(eta)
    residuals.append(residual)

# Visualisation
import matplotlib.pyplot as plt

plt.figure()
plt.loglog(ns, errors, label="Erreur $\epsilon_n$", marker='o')
plt.loglog(ns, etas, label="Borne supérieure $\eta_n$", marker='x')
plt.loglog(ns, residuals, label="Résidu normalisé $r_n$", marker='.')
plt.xlabel("n")
plt.ylabel("Valeur")
plt.legend()
plt.title("Log-log plot")
plt.grid(True)

plt.figure()
plt.semilogy(ns, errors, label="Erreur $\epsilon_n$", marker='o')
plt.semilogy(ns, etas, label="Borne supérieure $\eta_n$", marker='x')
plt.semilogy(ns, residuals, label="Résidu normalisé $r_n$", marker='.')
plt.xlabel("n")
plt.ylabel("Valeur")
plt.legend()
plt.title("Semilogy plot")
plt.grid(True)

plt.show()
In [22]:
### 4. Matrice et visualisation pour une nouvelle matrice

#(a) Construction de la matrice :


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


# (b) Répéter le point 3 :


ns = list(range(5, 101, 5))
errors_new = []
etas_new = []
residuals_new = []

for n in ns:
    A = create_matrix(n)
    b = np.full(n, 2)
    x_exact = np.array([i * (n - i + 1) for i in range(1, n+1)])
    x_computed = np.linalg.solve(A, b)

    error = np.linalg.norm(x_computed - x_exact, ord=2) / np.linalg.norm(x_exact, ord=2)
    eta = epsilon * np.linalg.cond(A, p=2)
    residual = np.linalg.norm(b - np.dot(A, x_computed), ord=2) / (np.linalg.norm(np.dot(A, x_computed)))

    errors_new.append(error)
    etas_new.append(eta)
    residuals_new.append(residual)

# Visualisation
plt.figure()
plt.loglog(ns, errors_new, label="Erreur $\epsilon_n$", marker='o')
plt.loglog(ns, etas_new, label="Borne supérieure $\eta_n$", marker='x')
plt.loglog(ns, residuals_new, label="Résidu normalisé $r_n$", marker='.')
plt.xlabel("n")
plt.ylabel("Valeur")
plt.legend()
plt.title("Log-log plot pour nouvelle matrice")
plt.grid(True)

plt.figure()
plt.semilogy(ns, errors_new, label="Erreur $\epsilon_n$", marker='o')
plt.semilogy(ns, etas_new, label="Borne supérieure $\eta_n$", marker='x')
plt.semilogy(ns, residuals_new, label="Résidu normalisé $r_n$", marker='.')
plt.xlabel("n")
plt.ylabel("Valeur")
plt.legend()
plt.title("Semilogy plot pour nouvelle matrice")
plt.grid(True)

plt.show()

Exercice 3¶

  1. Créer une fonction pour obtenir la décomposition $QR$ d'une matrice (voir feuille de TD 1).
  2. La tester sur la matrice $$ A = \begin{pmatrix} -4 & 5 & 1 \\ -2 & 7 & -1 \\ 4 & 4 & 2 \end{pmatrix} $$
In [37]:
# 1. Décomposition QR
import numpy as np

def gram_schmidt(A):
    A = np.array(A, dtype=float)
    m, n = A.shape
    Q = np.zeros((m, n))
    R = np.zeros((n, n))

    for j in range(n):
        v = A[:, j].copy()
        for i in range(j):
            R[i, j] = np.dot(Q[:, i], A[:, j])
            v -= R[i, j] * Q[:, i]
        R[j, j] = np.linalg.norm(v)
        Q[:, j] = v / R[j, j]
    
    return Q, R
In [38]:
# 2. Application
A = np.array([
    [-4, 5, 1],
    [-2, 7, -1],
    [ 4, 4, 2]
])

Q, R = gram_schmidt(A)
Q, R
Out[38]:
(array([[-0.66666667,  0.33333333,  0.66666667],
        [-0.33333333,  0.66666667, -0.66666667],
        [ 0.66666667,  0.66666667,  0.33333333]]),
 array([[ 6., -3.,  1.],
        [ 0.,  9.,  1.],
        [ 0.,  0.,  2.]]))
In [39]:
print('Q^T Q =')
print(Q.T @ Q)

print('\nQR =')
print(Q @ R)

Q^T Q =
[[ 1.00000000e+00 -2.46716228e-17  2.00456935e-16]
 [-2.46716228e-17  1.00000000e+00  6.16790569e-18]
 [ 2.00456935e-16  6.16790569e-18  1.00000000e+00]]

QR =
[[-4.  5.  1.]
 [-2.  7. -1.]
 [ 4.  4.  2.]]
In [ ]: