Méthodes itératives pour les systèmes linéaires

Objectifs d'apprentissage

À la fin de cette leçon, vous serez en mesure de :

• Comprendre pourquoi on utilise des méthodes itératives plutôt que directes
• Implémenter les méthodes de Jacobi et Gauss-Seidel
• Comprendre intuitivement pourquoi ces méthodes convergent (ou pas !)
• Utiliser la relaxation (SOR) pour accélérer la convergence
• Savoir quand choisir quelle méthode

Prérequis

• Systèmes linéaires et élimination de Gauss
• Notion de convergence d'une suite

Pourquoi des méthodes itératives ?

Le problème des grands systèmes

Jusqu'ici, nous avons résolu avec des méthodes directes comme l'élimination de Gauss ou la factorisation LU. Ces méthodes donnent la solution exacte (aux erreurs d'arrondi près) en un nombre fini d'opérations.

Mais imaginez que vous devez résoudre un système avec n = 1 000 000 d'inconnues (ce qui arrive couramment en simulation physique, traitement d'images, apprentissage automatique...).

La bonne nouvelle : les matrices creuses

Dans la plupart des applications réelles, les matrices sont creuses (sparse) : la grande majorité des éléments sont nuls.

Exemple : En traitement d'images, chaque pixel n'interagit qu'avec ses voisins immédiats. Pour une image 1000×1000, on a pixels, mais chaque ligne de la matrice n'a que ~5 éléments non nuls au lieu de !

L'idée des méthodes itératives

Au lieu de calculer la solution exacte en une fois, on va construire une suite qui converge vers la solution.

Avantages :

• On ne stocke jamais la matrice complète, seulement les éléments non nuls
• On peut s'arrêter dès qu'on a la précision souhaitée
• Souvent, quelques dizaines d'itérations suffisent !

L'intuition derrière les méthodes itératives

L'idée simple : résoudre une équation à la fois

Prenons un système simple à 2 équations :

Idée naïve : isolons chaque variable dans "son" équation :

Maintenant, partons d'une estimation initiale (par exemple ) et itérons :

Itération 1 : et

Itération 2 : et

Itération 3 : ...

La suite converge vers qui est bien la solution exacte !

Visualisation interactive

Observez comment les deux méthodes (Jacobi et Gauss-Seidel) convergent vers la solution :

Remarquez que Gauss-Seidel converge plus vite avec le même nombre d'itérations. Nous allons comprendre pourquoi dans les sections suivantes.

Principe général : forme itérative et convergence

Avant d'entrer dans les détails des différentes méthodes, formalisons le cadre général.

La forme itérative

Toutes les méthodes itératives pour résoudre peuvent s'écrire sous une forme unifiée :

où :

• est la matrice d'itération (dépend de la méthode choisie)
• est un vecteur constant (dépend de et de la méthode)
• est l'approximation à l'itération

Chaque méthode (Jacobi, Gauss-Seidel, SOR) correspond à un choix particulier de et .

Le rayon spectral

La question fondamentale est : la suite converge-t-elle vers la solution ?

La réponse dépend entièrement du rayon spectral de la matrice d'itération .

Théorème de convergence

Interprétation intuitive : le rayon spectral mesure le "facteur d'amplification maximal" de la matrice. Si , chaque itération réduit l'erreur (en moyenne). Si , l'erreur est amplifiée et la méthode diverge.

De plus, plus est petit, plus la convergence est rapide. C'est pourquoi on cherche des méthodes dont la matrice d'itération a un petit rayon spectral.

En pratique

Calculer les valeurs propres d'une grande matrice est coûteux — parfois plus que de résoudre le système lui-même ! On utilise donc des critères suffisants plus faciles à vérifier, comme la dominance diagonale stricte (que nous verrons plus loin).

La décomposition D - L - U

Pour construire les matrices d'itération des différentes méthodes, on décompose la matrice en trois parties :

où :

• D (Diagonal) : la partie diagonale
• L (Lower) : la partie triangulaire inférieure stricte, avec signe opposé
• U (Upper) : la partie triangulaire supérieure stricte, avec signe opposé

Exemple concret :

Vérifions : ✓

Méthode de Jacobi

Le principe

La méthode de Jacobi utilise l'idée simple qu'on a vue : on isole chaque variable et on met à jour toutes les composantes simultanément avec les valeurs de l'itération précédente.

Pour la i-ème équation , on isole :

En français : la nouvelle valeur de est calculée à partir des anciennes valeurs de toutes les autres variables.

Forme matricielle

En notation matricielle :

La matrice d'itération est .

Méthode de Gauss-Seidel

L'amélioration clé

Dans Jacobi, on utilise les anciennes valeurs pour tout. Mais si on vient de calculer , pourquoi ne pas l'utiliser immédiatement pour calculer ?

C'est l'idée de Gauss-Seidel : utiliser les valeurs les plus récentes disponibles.

• Pour : on utilise les nouvelles valeurs (en vert)
• Pour : on utilise les anciennes valeurs (en bleu)

Exemple pas à pas

Reprenons notre système :

Jacobi (pour comparaison) :

Gauss-Seidel :

• ← On utilise le x qu'on vient de calculer !

Après une seule itération, Gauss-Seidel a déjà alors que Jacobi a . La solution exacte est , donc Gauss-Seidel est plus proche !

Forme matricielle

Matrice d'itération :

Comparaison des vitesses de convergence

Le graphique montre que Gauss-Seidel converge environ 2× plus vite que Jacobi. Mais peut-on faire encore mieux ? Oui, avec la relaxation !

Méthode SOR (Successive Over-Relaxation)

L'idée de la relaxation

Gauss-Seidel calcule une correction à appliquer à chaque composante. Et si on amplifiait cette correction ?

Soit la valeur calculée par Gauss-Seidel. Au lieu de simplement prendre cette valeur, on fait :

Le paramètre (oméga) contrôle l'ampleur de la correction :

• : c'est exactement Gauss-Seidel
• : sur-relaxation — on va "plus loin" que Gauss-Seidel (accélération)
• : sous-relaxation — on y va doucement (stabilisation)

Visualisation interactive

Explorez l'effet du paramètre sur la convergence :

Trouver le ω optimal

Le choix de est crucial. Un mauvais choix peut rendre SOR plus lent que Gauss-Seidel !

L'objectif est de minimiser le rayon spectral de la matrice d'itération. Pour certaines matrices (notamment celles issues de la discrétisation d'équations aux dérivées partielles), on peut calculer le optimal :

où est le rayon spectral de la matrice de Jacobi (introduit en début de leçon).

En pratique, on fait souvent quelques essais autour de pour trouver une bonne valeur.

Conditions de convergence

Nous avons vu en introduction que la convergence dépend du rayon spectral : . Mais comment vérifier cette condition en pratique sans calculer les valeurs propres ?

Critère pratique : dominance diagonale

Pourquoi ça marche ? On peut démontrer que si est strictement diagonalement dominante, alors et . La dominance diagonale implique que les éléments hors diagonale ont une influence "limitée", ce qui se traduit par un rayon spectral inférieur à 1.

Exemple : convergence garantie

Vérifions la dominance diagonale :

• Ligne 1 : ✓
• Ligne 2 : ✓
• Ligne 3 : ✓

Cette matrice est strictement diagonalement dominante → → convergence garantie !

Contre-exemple : divergence

• Ligne 1 : ✗
• Ligne 2 : ✗

La diagonale ne domine pas. Que vaut le rayon spectral ? La matrice d'itération de Jacobi est :

Les valeurs propres de sont , donc : la méthode diverge.

Vérifions-le concrètement avec le système où :

Partons de :

• Itération 1 :
• Itération 2 :
• Itération 3 :
• Itération 4 :

Les valeurs oscillent avec une amplitude croissante : La suite diverge ! (La solution exacte est .)

Comparaison des méthodes

Quand utiliser quelle méthode ?

• Jacobi : quand vous avez beaucoup de processeurs (GPU, calcul distribué)
• Gauss-Seidel : choix par défaut pour un code séquentiel
• SOR : quand la convergence est lente et que vous pouvez estimer ω

Critères d'arrêt

Comment savoir quand s'arrêter ? Deux approches courantes :

1. Différence entre itérations

Avantage : facile à calculer Inconvénient : peut être trompeur si la convergence est lente

2. Résidu

Avantage : mesure directement à quel point on satisfait l'équation Inconvénient : peut être trompeur pour les matrices mal conditionnées (voir leçon précédente !)

Résumé

Points clés à retenir :

1. Principe général : toute méthode itérative s'écrit
2. Convergence : garantie si et seulement si (rayon spectral)
3. Jacobi : mise à jour simultanée, parallélisable
4. Gauss-Seidel : utilise les nouvelles valeurs immédiatement, plus rapide
5. SOR : accélère Gauss-Seidel avec un paramètre
6. Critère pratique : la dominance diagonale stricte garantit
7. Applications : grands systèmes creux (simulations, traitement d'images...)

Implémentations Python

Méthode de Jacobi

import numpy as np

def jacobi(A, b, x0=None, max_iter=100, tol=1e-8):
  """
  Méthode de Jacobi pour résoudre Ax = b.

  Paramètres:
  -----------
  A : matrice (n×n) du système
  b : vecteur second membre
  x0 : estimation initiale (zéro par défaut)
  max_iter : nombre maximum d'itérations
  tol : tolérance pour le critère d'arrêt

  Retourne:
  ---------
  x : solution approchée
  """
  n = len(b)
  x = x0.copy() if x0 is not None else np.zeros(n)

  # Extraire la diagonale et le reste
  D = np.diag(A)
  R = A - np.diag(D)  # R = -(L + U) avec les signes d'origine

  for k in range(max_iter):
      # Formule de Jacobi : x_new = D^{-1}(b - R·x)
      x_new = (b - R @ x) / D

      # Test de convergence : ||x_new - x||_∞ < tol
      if np.linalg.norm(x_new - x, np.inf) < tol:
          print(f"Convergence en {k+1} itérations")
          return x_new

      x = x_new

  print(f"Non convergé après {max_iter} itérations")
  return x

Méthode de Gauss-Seidel

def gauss_seidel(A, b, x0=None, max_iter=100, tol=1e-8):
  """
  Méthode de Gauss-Seidel pour résoudre Ax = b.

  Plus rapide que Jacobi car utilise immédiatement
  les nouvelles valeurs calculées.
  """
  n = len(b)
  x = x0.copy() if x0 is not None else np.zeros(n)

  for k in range(max_iter):
      x_old = x.copy()  # Pour le test de convergence

      for i in range(n):
          # Somme des termes avec les NOUVELLES valeurs (j < i)
          sum_new = np.dot(A[i, :i], x[:i])

          # Somme des termes avec les ANCIENNES valeurs (j > i)
          sum_old = np.dot(A[i, i+1:], x[i+1:])

          # Mise à jour immédiate de x[i]
          x[i] = (b[i] - sum_new - sum_old) / A[i, i]

      # Test de convergence
      if np.linalg.norm(x - x_old, np.inf) < tol:
          print(f"Convergence en {k+1} itérations")
          return x

  print(f"Non convergé après {max_iter} itérations")
  return x

Méthode SOR

def sor(A, b, omega=1.5, x0=None, max_iter=100, tol=1e-8):
  """
  Méthode SOR (Successive Over-Relaxation).

  Paramètres:
  -----------
  omega : paramètre de relaxation (0 < omega < 2)
          - omega = 1 : Gauss-Seidel
          - omega > 1 : sur-relaxation (accélération)
          - omega < 1 : sous-relaxation (stabilisation)
  """
  n = len(b)
  x = x0.copy() if x0 is not None else np.zeros(n)

  for k in range(max_iter):
      x_old = x.copy()

      for i in range(n):
          # Calculer la valeur Gauss-Seidel
          sum_new = np.dot(A[i, :i], x[:i])
          sum_old = np.dot(A[i, i+1:], x[i+1:])
          x_gs = (b[i] - sum_new - sum_old) / A[i, i]

          # Appliquer la relaxation
          x[i] = (1 - omega) * x_old[i] + omega * x_gs

      if np.linalg.norm(x - x_old, np.inf) < tol:
          print(f"SOR (omega={omega}): Convergence en {k+1} itérations")
          return x

  print(f"Non convergé après {max_iter} itérations")
  return x

Exemple d'utilisation

# Système test
A = np.array([[4, -1, 0],
            [-1, 4, -1],
            [0, -1, 4]], dtype=float)
b = np.array([15, 10, 10], dtype=float)

print("=== Jacobi ===")
x_j = jacobi(A, b)

print("\n=== Gauss-Seidel ===")
x_gs = gauss_seidel(A, b)

print("\n=== SOR (omega=1.2) ===")
x_sor = sor(A, b, omega=1.2)

print(f"\nSolution : {x_sor}")
print(f"Vérification A·x = {A @ x_sor}")

Pour aller plus loin

La prochaine leçon abordera les systèmes d'équations non linéaires : que faire quand les équations ne sont plus linéaires en les inconnues ? On verra comment généraliser la méthode de Newton au cas multidimensionnel.