Résolution de systèmes linéaires

Cours · Python (école d'ingénieur), chapitre 12 · école d'ingénieur, informatique · le début du chapitre, la suite sur Adloun

12.1 Introduction et Motivation

De nombreux problèmes en physique, en ingénierie, et en apprentissage automatique se ramènent à la résolution d'un système d'équations linéaires de la forme $A x = b$ , où $A \in M_{n} (R)$ est une matrice inversible et $b \in R^{n}$ est un vecteur colonne donné. Ce chapitre présente l'implémentation de la méthode du pivot de Gauss à l'aide de la bibliothèque NumPy, ainsi que des techniques plus avancées comme la décomposition LU et l'analyse du conditionnement. Nous y aborderons la manipulation avancée des tableaux bidimensionnels (coupes, copies, fancy indexing), les opérations élémentaires sur les lignes d'une matrice, et l'importance cruciale du choix du pivot partiel pour assurer la stabilité numérique de l'algorithme face aux erreurs d'arrondis machine.

12.2 Calcul Matriciel avec NumPy

Pour modéliser des vecteurs et des matrices, NumPy fournit le type ndarray. Contrairement aux listes Python, les tableaux NumPy ont une taille fixe et sont homogènes (tous les éléments ont le même type, généralement float).

12.2.1 Création et Types de Données


import numpy as np

# Matrice à partir d'une liste de listes
A = np.array([[3.0, 6.0, 2.0], [2.0, -3.0, 5.0]], dtype=float)

# Matrices particulières
I = np.eye(3)        # Matrice identité 3x3
Z = np.zeros((4, 4)) # Matrice nulle 4x4
O = np.ones((2, 3))  # Matrice de 1 2x3

12.2.2 Coupes et Vues vs Copies

Une coupe (slicing) d'un tableau NumPy ne crée pas un nouvel objet en mémoire mais une vue (un pointeur vers la même zone mémoire). Pour dupliquer réellement les données, il faut utiliser la méthode copy() :


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

# Vue : modification affecte A
ligne1 = A[0]
ligne1[0] = 99
print(A)  # [[99, 2], [3, 4]]

# Copie : modification n'affecte pas A
ligne2 = A[1].copy()
ligne2[0] = 88
print(A)  # [[99, 2], [3, 4]] (inchangé)

12.2.3 Fancy Indexing (Indexation Avancée)

NumPy permet d'échanger des lignes ou des colonnes de manière concise :


# Permutation des lignes i et j
A[[i, j]] = A[[j, i]]

# Permutation des colonnes i et j
A[:, [i, j]] = A[:, [j, i]]

12.2.4 Produit Matriciel

Attention, l'opérateur * effectue un produit élément par élément (produit de Hadamard). Le produit matriciel mathématique $A x$ ou $A B$ s'effectue à l'aide de la fonction np.dot() ou de l'opérateur @ :


A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

C = A @ B  # Produit matriciel
# [[1*5+2*7, 1*6+2*8], [3*5+4*7, 3*6+4*8]] = [[19, 22], [43, 50]]

12.3 Méthode du Pivot de Gauss

Pour résoudre $A x = b$ , nous appliquons des opérations élémentaires sur les lignes du système augmenté $[A ∣ b]$ afin de rendre la matrice triangulaire supérieure (descente), puis diagonale (remontée), avant d'isoler les inconnues.

12.3.1 L'Importance de la Stabilité : Choix du Pivot Partiel

Le choix du pivot lors de l'élimination est crucial. Diviser par un pivot trop petit en valeur absolue amplifie dramatiquement les erreurs d'arrondis.

Exemple 12.1Désastre numérique avec un petit pivot

Considérons le système à 4 chiffres significatifs en base 10 :

${0.003 x + 59.14 y = 59.17 5.291 x - 6.130 y = 46.78$

La solution exacte est $x = 10$ , $y = 1$ .

Sans permutation (pivot $0.003$ très petit) : Le calcul donne la constante multiplicatrice $λ = 5.291/0.003 \approx 1764$ . L'opération $L_{2} \leftarrow L_{2} - λ L_{1}$ conduit aux approximations d'arrondis. Lors de la remontée, on trouve $y \approx 1.001$ et $x \approx - 6.667$ (soit 167% d'erreur sur $x$ !).
Avec permutation (pivot partiel) : On cherche le coefficient maximal en valeur absolue sur la colonne 1, qui est $5.291$ . On permute $L_{1} \leftrightarrow L_{2}$ . La constante devient $λ = 0.003/5.291 \approx 5.670 \cdot 1 0^{- 4}$ . La remontée donne $y \approx 1.000$ et $x \approx 10.00$ .

Définition 12.2Stratégie du Pivot Partiel

Avant de traiter la colonne $j$ de la matrice, on recherche la ligne $i \geq j$ pour laquelle $∣ a_{i, j} ∣$ est maximal. On permute alors les lignes $i$ et $j$ de $A$ et $b$ .

12.3.2 Algorithme et Implémentation en Python


def recherche_pivot(A, b, j):
    """
    Recherche le pivot partiel sur la colonne j.
    Permute les lignes j et i_max dans A et b.
    """
    n = A.shape[0]
    i_max = j
    for i in range(j + 1, n):
        if abs(A[i, j]) > abs(A[i_max, j]):
            i_max = i
    if i_max != j:
        A[[j, i_max]] = A[[i_max, j]]
        b[[j, i_max]] = b[[i_max, j]]

def descente(A, b):
    """
    Étape de descente : transforme A en matrice triangulaire supérieure.
    """
    n = A.shape[0]
    for j in range(n - 1):
        recherche_pivot(A, b, j)
        pivot = A[j, j]
        if abs(pivot) < 1e-15:
            raise ValueError("Matrice singulière ou presque singulière.")
        for i in range(j + 1, n):
            facteur = A[i, j] / pivot
            A[i] = A[i] - facteur * A[j]
            b[i] = b[i] - facteur * b[j]

def remontee(A, b):
    """
    Étape de remontée : résout le système triangulaire supérieur.
    """
    n = A.shape[0]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        somme = 0.0
        for j in range(i + 1, n):
            somme += A[i, j] * x[j]
        x[i] = (b[i] - somme) / A[i, i]
    return x

def solve_gauss(A, b):
    """
    Résout le système Ax = b par la méthode de Gauss avec pivot partiel.
    """
    Ac = A.astype(float).copy()
    bc = b.astype(float).copy()
    descente(Ac, bc)
    return remontee(Ac, bc)

12.3.3 Complexité Algorithmique

Proposition 12.3Complexité de la méthode de Gauss

Pour un système de $n$ équations à $n$ inconnues :

La recherche du pivot sur la colonne $j$ effectue $O (n - j)$ comparaisons.
L'élimination lors de la descente requiert $O (n^{3})$ opérations.
La remontée requiert $O (n^{2})$ opérations.

La complexité temporelle globale est donc $O (n^{3})$ , et la complexité spatiale $O (n^{2})$ .

12.4 Décomposition LU

La décomposition $A = LU$ (où $L$ est triangulaire inférieure avec des 1 sur la diagonale et $U$ triangulaire supérieure) est intimement liée à l'élimination de Gauss. Les facteurs d'élimination sont stockés dans $L$ et la matrice triangulée dans $U$ .

12.4.1 Principe et Implémentation


def lu_decomposition(A):
    """
    Décompose A = LU sans pivotage.
    Retourne L et U.
    """
    n = A.shape[0]
    L = np.eye(n)
    U = A.astype(float).copy()

    for j in range(n - 1):
        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

# Résolution de Ax = b via LU : LUx = b
def solve_lu(A, b):
    L, U = lu_decomposition(A)

    # Ly = b (descente)
    n = len(b)
    y = np.zeros(n)
    for i in range(n):
        y[i] = b[i] - np.dot(L[i, :i], y[:i])

    # Ux = y (remontée)
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]
    return x

12.4.2 Avantages de LU

Une fois la décomposition $A = LU$ calculée (coût en $O (n^{3})$ ), on peut résoudre efficacement des systèmes avec plusieurs seconds membres $b_{1}, b_{2}, \dots$ en $O (n^{2})$ opérations chacun.

12.5 Déterminant et Inverse

12.5.1 Déterminant via la décomposition LU

Le déterminant d'une matrice triangulaire est le produit de ses éléments diagonaux. Pour $P A = LU$ (avec pivotage partiel), le déterminant est :

$det (A) = det (P^{- 1}) \cdot det (L) \cdot det (U)$

Comme $det (L) = 1$ et $det (P^{- 1}) = (- 1)^{s}$ où $s$ est le nombre de permutations de lignes effectuées, on a :

$det (A) = (- 1)^{s} i = 0 \prod n - 1 U [i, i]$


def determinant(A):
    """Calcule le déterminant par la méthode de Gauss avec pivot partiel."""
    A_copy = A.astype(float).copy()
    n = A_copy.shape[0]
    det_sign = 1

    for j in range(n - 1):
        # Recherche du pivot
        i_max = j
        for i in range(j + 1, n):
            if abs(A_copy[i, j]) > abs(A_copy[i_max, j]):
                i_max = i
        if i_max != j:
            A_copy[[j, i_max]] = A_copy[[i_max, j]]
            det_sign *= -1

        if abs(A_copy[j, j]) < 1e-15:
            return 0.0

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

    det = det_sign * np.prod([A_copy[i, i] for i in range(n)])
    return det

12.6 Conditionnement d'une Matrice

Le conditionnement d'une matrice mesure sa sensibilité aux erreurs sur les données.

Définition 12.4Conditionnement

Le conditionnement d'une matrice $A$ (pour la norme 2) est :

$κ (A) = ∥ A ∥_{2} \cdot ∥ A^{- 1} ∥_{2} = \frac{σ _{m a x} ( A )}{σ _{m i n} ( A )}$

où $σ_{m a x} (A)$ et $σ_{m i n} (A)$ sont les valeurs singulières extrêmes de $A$ .

Un système est bien conditionné si $κ (A)$ est proche de 1, et mal conditionné si $κ (A)$ est grand.

Exemple 12.5Matrice de Hilbert

La matrice de Hilbert $H_{n}$ d'éléments $H_{i, j} = \frac{1}{i + j - 1}$ est célèbre pour être extrêmement mal conditionnée même pour de petites valeurs de $n$ . Par exemple, pour $n = 4$ :

$H_{4} = 1 \frac{1}{2} \frac{1}{3} \frac{1}{4} \frac{1}{2} \frac{1}{3} \frac{1}{4} \frac{1}{5} \frac{1}{3} \frac{1}{4} \frac{1}{5} \frac{1}{6} \frac{1}{4} \frac{1}{5} \frac{1}{6} \frac{1}{7}$


def hilbert(n):
    """Génère la matrice de Hilbert d'ordre n."""
    H = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            H[i, j] = 1.0 / (i + j + 1)
    return H

# Analyse du conditionnement
for n in [2, 3, 4, 5, 6, 7, 8]:
    H = hilbert(n)
    cond = np.linalg.cond(H)
    print(f"n = {n}, kappa(H) = {cond:.2e}")

12.7 Utilisation de SciPy

Pour le calcul scientifique professionnel, on utilise les fonctions du module scipy.linalg :


from scipy.linalg import solve, lu, det, inv, cond

# Résolution directe
x = solve(A, b)

# Décomposition LU
P, L, U = lu(A)  # PA = LU

# Déterminant
d = det(A)

# Inverse
A_inv = inv(A)

# Conditionnement
c = cond(A)

12.8 Exercices résolus

Niveau 1 : Manipulation de matrices avec NumPy

Exercice 1 - Création et opérations de base

Créer les matrices suivantes et effectuer les opérations demandées :

$A = (1324)$ et $B = (5768)$ .
Calculer $A + B$ , $A - B$ , $A \times B$ (produit matriciel).
Calculer la transposée de $A$ et la trace de $A$ .

Démonstration (Solution)


import numpy as np

# Création des matrices
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

print("A =")
print(A)
print("\nB =")
print(B)

# Opérations élément par élément
print("\nA + B =")
print(A + B)

print("\nA - B =")
print(A - B)

# Produit matriciel
print("\nA @ B (produit matriciel) =")
print(A @ B)

# Transposée
print("\nTransposée de A =")
print(A.T)

# Trace (somme des éléments diagonaux)
print(f"\nTrace de A = {np.trace(A)}")

Résultat dans la console :


A =
[[1 2]
 [3 4]]

B =
[[5 6]
 [7 8]]

A + B =
[[ 6  8]
 [10 12]]

A - B =
[[-4 -4]
 [-4 -4]]

A @ B (produit matriciel) =
[[19 22]
 [43 50]]

Transposée de A =
[[1 3]
 [2 4]]

Trace de A = 5

Exercice 2 - Vues vs copies

Expliquer et illustrer la différence entre une vue et une copie dans NumPy.

La suite du cours est sur Adloun

Le chapitre en entier, formules et figures comprises, avec le tuteur, les fiches, les QCM et les exercices.

Lire la suite du cours

En accès libre, sans inscription.

Exercices corrigés de ce chapitre

Lire la suite du cours →