Ce document pédagogique, conçu pour les Travaux Pratiques d'Analyse Numérique à l'Université de Skikda, est destiné aux étudiants de 2ème année LMD Sciences et Technologie. Il vise à explorer les différentes approches de résolution des systèmes d'équations linéaires, tant directes qu'itératives, en mettant l'accent sur leur implémentation pratique.
Il couvre les notions suivantes :
- La méthode directe de Gauss (sans pivotation)
- La méthode directe de décomposition LU (Crout)
- La méthode directe de Cholesky
- Les méthodes itératives de Jacobi et Gauss-Seidel
Chaque section présente le cadre théorique et est illustrée par des exemples concrets et des scripts Matlab.
TP N°2 : Résolution des systèmes d'équations linéaires en Analyse Numérique
Ce Travaux Pratiques (TP) est destiné aux étudiants de 2ème année LMD Sciences et Technologie en Analyse Numérique. Il couvre les principales méthodes numériques pour la résolution de systèmes d'équations linéaires, tant directes qu'itératives, et leur implémentation en Matlab.
1) La méthode directe de Gauss (sans pivotage)
La méthode de Gauss, ou élimination de Gauss, est une technique directe pour résoudre des systèmes d'équations linéaires. Elle consiste à transformer le système en une forme triangulaire supérieure, puis à résoudre le système résultant par substitution rétrograde. Cette version est sans pivotage, ce qui signifie qu'elle peut échouer si un élément diagonal devient nul ou très petit, ou si le système est singulier.
Implémentation Matlab de la méthode de Gauss
function x = gauss(A,b)
n = length(b);
% ------------------ Triangularisation -------------------
for k = 1:n-1
for i = k+1:n
m = A(i,k)/A(k,k);
A(i,k+1:n) = A(i,k+1:n) - m*A(k,k+1:n);
b(i) = b(i) - m*b(k);
end
end
% ---------------------- Solution ------------------------
x = zeros(n,1); % Initialisation du vecteur solution
x(n) = b(n)/A(n,n);
for i = n-1:-1:1
somme = sum(A(i,i+1:n).*x(i+1:n));
x(i) = (b(i)-somme)/A(i,i);
end
x = x';
end
Description des entrées et sorties de la fonction `gauss`
| Entrées | Description | Sorties | Description |
|---|---|---|---|
| A | La matrice des coefficients | x | Le vecteur solution |
| b | Le vecteur du second membre |
Exemples d'utilisation en Matlab
Exemple 1
>> A = [2 -1 1 ; 3 2 -5 ; 1 3 -2];
>> b = [0 ; 1 ; 4];
>> X = gauss(A,b)
X =
0.4643
1.6786
0.7500
Exemple 2
>> A = [4 -1 1 ; 4 -8 1 ; -2 1 5];
>> b = [7 ; -21 ; 15];
>> X = gauss(A,b)
X =
2
4
3
Exemple 3
>> A = [1 3 3 ; 2 2 5 ; 3 2 6];
>> b = [-2 ; 7 ; 12];
>> X = gauss(A,b)
X =
4
-3
1
Exemple 4 : Système non résolvable
>> A = [-1 3 -4 ; 1 -5 4 ; 2 1 8];
>> b = [2;-1;5];
>> X = gauss(A,b)
X =
NaN
NaN
Inf
Le quatrième système ne peut pas être résolu car son déterminant est nul (la matrice A est singulière).
2) La méthode directe de décomposition LU (Crout)
La décomposition LU (Lower-Upper) est une méthode directe qui factorise une matrice A en le produit d'une matrice triangulaire inférieure L et d'une matrice triangulaire supérieure U (A = LU). La méthode de Crout est une variante spécifique où les éléments diagonaux de la matrice U sont tous égaux à 1. Une fois la décomposition effectuée, la résolution du système Ax=b est simplifiée en deux systèmes triangulaires Ly=b et Ux=y.
Implémentation Matlab de la décomposition LU (Crout)
function [L,U,x] = decompositionLU(A,b)
n = length(b);
L = zeros(n,n);
U = zeros(n,n);
x = zeros(n,1);
% --------- Factorisation LU ------------------------------------
for i = 1:n , U(i,i) = 1; end
for k = 1:n
for i = k:n
L(i,k) = A(i,k)-L(i,1:k-1)*U(1:k-1,k);
end
for j = k+1:n
U(k,j) = (A(k,j)-L(k,1:k-1)*U(1:k-1,j))/L(k,k);
end
end
% ---------- Résoudre Ly = b ------------------------------------
y = zeros(1,n);
for i = 1:n
somme = sum(L(i,1:i-1).*y(1:i-1));
y(i) = (b(i)-somme)/L(i,i);
end
% ---------- Résoudre Ux = y ------------------------------------
for i = n:-1:1
somme = U(i,i+1:n)*x(i+1:n);
x(i) = y(i)-somme;
end
end
Description des entrées et sorties de la fonction `decompositionLU`
| Entrées | Description | Sorties | Description |
|---|---|---|---|
| A | La matrice des coefficients | L | La matrice triangulaire inférieure L |
| b | Le vecteur du second membre | U | La matrice triangulaire supérieure U |
| x | Le vecteur solution |
Exemple d'utilisation en Matlab
>> A = [2 -1 1 ; 3 2 -5 ; 1 3 -2];
>> b = [0 ; 1 ; 4];
>> [L,U,X] = decompositionLU(A,b)
L =
2.0000 0 0
3.0000 3.5000 0
1.0000 3.5000 4.0000
U =
1.0000 -0.5000 0.5000
0 1.0000 -1.8571
0 0 1.0000
X =
0.4643
1.6786
0.7500
On peut vérifier si L*U = A en utilisant la commande :
>> L*U
ans =
2 -1 1
3 2 -5
1 3 -2
3) La méthode directe de Cholesky
La méthode de Cholesky est une méthode de décomposition directe efficace pour résoudre des systèmes d'équations linéaires. Elle est applicable uniquement pour les systèmes dont la matrice A est symétrique (A = Aᵀ) et définie positive. Elle factorise A en un produit d'une matrice triangulaire inférieure L et de sa transposée Lᵀ (A = LLᵀ). Cette décomposition est souvent plus rapide et numériquement plus stable que la décomposition LU pour ce type de matrices.
Génération de matrices de test symétriques définies positives en Matlab
Pour obtenir une telle matrice, on peut utiliser les fonctions Matlab suivantes :
A = gallery('moler',n)
Ou bien
A = gallery('lehmer',n)
Où n représente le nombre de lignes/colonnes de la matrice.
Implémentation Matlab de la méthode de Cholesky
function [L,x] = cholesky(A,b)
n = length(b);
L = zeros(n,n);
x = zeros(n,1);
% ----------------------- Factorisation LL' ----------
for i = 1:n
L(i,i) = sqrt(A(i,i)-sum(L(i,1:i-1).^2));
for j = i+1:n
L(j,i) = (A(i,j)-sum(L(i,1:i-1).*L(j,1:i-1)))/L(i,i);
end
end
% ------------------------- Résoudre Ly = b ----------
y = zeros(1,n);
for i = 1:n
somme = sum(L(i,1:i-1).*y(1:i-1));
y(i) = (b(i)-somme)/L(i,i);
end
% ------------------------- Résoudre L'x = y ---------
for i = n:-1:1
somme = sum(L(i+1:n,i).*x(i+1:n));
x(i) = (y(i)-somme)/L(i,i);
end
end
Description des entrées et sorties de la fonction `cholesky`
| Entrées | Description | Sorties | Description |
|---|---|---|---|
| A | La matrice des coefficients | L | La matrice triangulaire inférieure L (avec LL'=A) |
| b | Le vecteur du second membre | x | Le vecteur solution |
Exemples d'utilisation en Matlab
Exemple 1
>> A = gallery('moler',4)
A =
1 -1 -1 -1
-1 2 0 0
-1 0 3 1
-1 0 1 4
>> b = [5 ; -1 ; 4 ; 3];
>> [L,X] = cholesky(A,b)
L =
1 0 0 0
-1 1 0 0
-1 -1 1 0
-1 -1 -1 1
X =
135
67
38
25
On peut vérifier si la réponse est correcte en utilisant la division matricielle inverse de Matlab :
>> A\b
ans =
135
67
38
25
On peut aussi vérifier si L*L' = A en écrivant :
>> L*L'
ans =
1 -1 -1 -1
-1 2 0 0
-1 0 3 1
-1 0 1 4
>> L*L'-A
ans =
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
Exemple 2
>> A = gallery('lehmer',5) ;
>> b = rand(5,1) ; % des valeurs aléatoires
>> [L,X] = cholesky(A,b)
4) Les méthodes itératives de Jacobi et Gauss-Seidel
Les méthodes itératives sont utilisées pour résoudre des systèmes d'équations linéaires en partant d'une approximation initiale et en l'affinant progressivement à chaque itération jusqu'à atteindre une solution convergente. Les méthodes de Jacobi et de Gauss-Seidel sont des exemples classiques de ces techniques.
Pour le système suivant, nous allons utiliser les méthodes de Jacobi et Gauss-Seidel (en considérant que le vecteur initial est (0,0,0)T) :
4x₁ - x₂ + x₃ = 7
4x₁ - 8x₂ + x₃ = -21
-2x₁ + x₂ + 5x₃ = 15
Configuration du système en Matlab
>> A = [4,-1,1 ; 4,-8,1 ; -2,1,5]
A =
4 -1 1
4 -8 1
-2 1 5
>> b = [7;-21;15]
b =
7
-21
15
>> X0 = [0;0;0]
X0 =
0
0
0
Implémentation Matlab de la méthode de Jacobi
function [X,niter] = jacobi(A,b,X0,nmax,tol)
n = length(b);
X = X0;
for niter = 1:nmax % Calculer l'itération suivante
for i = 1:n
j = [1:i-1,i+1:n];
somme = A(i,j)*X0(j);
X(i) = (b(i)-somme)/A(i,i);
end
% Tester la convergence
if norm(X-X0) < tol
return
end
% L'ancien X0 devient le nouveau X
X0 = X;
end
% En cas de divergence
disp('Pas de convergence')
end
Description des entrées et sorties de la fonction `jacobi`
| Entrées | Description | Sorties | Description |
|---|---|---|---|
| A | La matrice A | X | Le vecteur solution. |
| b | Le vecteur b | niter | Le nombre d’itérations effectuées. |
| X0 | Le vecteur initial | ||
| nmax | Le nombre maximal d’itérations | ||
| tol | L’erreur tolérée pour le résultat (tolérance) |
Implémentation Matlab de la méthode de Gauss-Seidel
function [X,niter] = gseidel(A,b,X0,nmax,tol)
n = length(b);
X = X0;
for niter = 1:nmax % Calculer l'itération suivante
for i = 1:n
somme1 = A(i,1:i-1)*X(1:i-1);
somme2 = A(i,i+1:n)*X0(i+1:n);
X(i)=(b(i)-somme1-somme2)/A(i,i);
end
% Tester la convergence
if norm(X-X0) < tol
return
end
% L'ancien X0 devient le nouveau X
X0 = X;
end
% En cas de divergence
disp('Pas de convergence')
end
Comparaison des méthodes de Jacobi et Gauss-Seidel
>> [X,N] = jacobi(A,b,X0,50,0.0001)
X =
2.0000
4.0000
3.0000
N =
11
>> [X,N] = gseidel(A,b,X0,50,0.0001)
X =
2.0000
4.0000
3.0000
N =
7
Pour observer les itérations, on peut modifier le programme de la manière suivante :
Ajouter X' juste après la ligne X0 = X; et avant le end de la boucle for niter.
Puis, on exécute les commandes Matlab suivantes avec format compact pour un affichage plus concis :
>> format compact
>> [X,N] = jacobi(A,b,X0,50,0.0001)
ans =
1.7500
2.6250
3.0000
ans =
1.6563
3.8750
3.1750
ans =
1.9250
3.8500
2.8875
ans =
1.9906
3.9484
3.0000
ans =
1.9871
3.9953
3.0066
ans =
1.9972
3.9944
2.9958
ans =
1.9996
3.9981
3.0000
ans =
1.9995
3.9998
3.0002
ans =
1.9999
3.9998
2.9998
ans =
2.0000
3.9999
3.0000
X =
2.0000
4.0000
3.0000
N =
11
>> [X,N] = gseidel(A,b,X0,50,0.0001)
ans =
1.7500
3.5000
3.0000
ans =
1.8750
3.9375
2.9625
ans =
1.9937
3.9922
2.9991
ans =
1.9983
3.9990
2.9995
ans =
1.9999
3.9999
3.0000
ans =
2.0000
4.0000
3.0000
X =
2.0000
4.0000
3.0000
N =
7
On constate que la méthode de Gauss-Seidel converge plus rapidement (en moins d'itérations) que celle de Jacobi pour cet exemple.
Foire Aux Questions (FAQ)
Quels sont les principaux types de méthodes pour résoudre les systèmes linéaires ?
Il existe deux grandes catégories : les méthodes directes et les méthodes itératives. Les méthodes directes (comme Gauss, LU, Cholesky) calculent la solution en un nombre fini d'étapes. Les méthodes itératives (comme Jacobi, Gauss-Seidel) partent d'une approximation et l'améliorent progressivement jusqu'à convergence, nécessitant des conditions spécifiques pour garantir cette convergence.
Quand utiliser une méthode directe ou une méthode itérative ?
Les méthodes directes sont généralement préférables pour les petits systèmes ou les matrices denses, car elles sont précises et robustes. Cependant, elles peuvent être coûteuses en mémoire pour de très grands systèmes. Les méthodes itératives sont plus adaptées aux grands systèmes ou aux matrices creuses, car elles peuvent être plus économes en mémoire et en temps de calcul, surtout si une bonne approximation initiale est disponible. Elles sont aussi essentielles lorsque la matrice est trop grande pour être stockée entièrement.
Qu'est-ce qu'une matrice symétrique définie positive et pourquoi est-elle importante pour la décomposition de Cholesky ?
Une matrice symétrique est une matrice qui est égale à sa transposée (A = Aᵀ). Elle est définie positive si le produit xᵀAx est strictement supérieur à zéro pour tout vecteur non nul x. Cette propriété est cruciale pour la décomposition de Cholesky car elle garantit que tous les pivots sont positifs et que l'extraction de la racine carrée des éléments diagonaux est toujours possible et donne des nombres réels. La décomposition de Cholesky est spécifiquement conçue pour ces matrices, offrant une méthode de factorisation numériquement stable et efficace (A = LLᵀ).