Chapitre 2
Résolution des systèmes d'équations linéaires
2.1 Introduction
2.2
Méthode de Gauss
2.2.1 Exposé de la méthode
2.2.2 Méthode de Gauss à pivot partiel
2.2.3 Méthode de Gauss à pivot total
2.2.4 Nombre d'opérations effectuées
2.3
Méthode de décomposition QR (Méthode De Householder)
2.4
Factorisation de Cholesky
2.5
Méthodes itératives
2.5.1 Méthode de Jacoby
2.5.2 Méthode de Gauss-Seidel
2.6
Quantification de l'erreur et de précision
2.7
Inversion de matrice
2.1 Introduction
Le but de ce chapitre est de passer en revue quelques
méthodes de résolution des systèmes d'équations linéaires. On supposera que
dans ce cas on dispose d'autant d'équations que d'inconnues. On examinera par
la suite le cas où cette condition n'est pas satisfaite.
Considérons le système d'équations linéaires suivant que l'on cherche à résoudre :
|
(2.1) |
Ce système peut aussi s'écrire :
|
(2.2) |
avec A = [aij ](i,j = 1...n), xT = [x1,...,xn] et bT = [b1,...,bn].
Théorème 23 Le système linéaire précédent admet une solution unique si et seulement si la matrice A est inversible (on dit aussi matrice non singulière). La solution est donnée par x = A-1.b.
Il existe différentes méthodes de résolution de ces systèmes. Nous examinerons les avantages et les inconvénients de chacune d'elles.
2.2 Méthode de Gauss
La méthode de Gauss est basée sur le constat suivant : le
système linéaire reste invariant pour les trois opérations suivantes effectuées
dans n'importe quel ordre et un nombre de fois indéterminé :
- Permutation de lignes de la matrice A (et donc
de b)
- Multiplication d'une ligne par une constante non nulle
- Addition d'une ligne à une autre ligne
2.2.1 Exposé de la méthode
En combinant les opérations précédentes, on peut transformer
la matrice A en une matrice triangulaire supérieure par le processus
itératif suivant :
- On garde la première ligne inchangée
- On suppose a11 ¹
0. On soustrait à la ième ligne de A et de
b (i = 2,...,n) la première ligne multipliée par la
quantité (ai1/a11). Les
matrices A et b prennent alors la forme :
|
(2.3) |
|
On répète l'étape précédente en supposant cette fois-ci a22(1)
¹ 0. On garde alors les deux premières
lignes inchangées et on retranche (ai1(1)/a22(1))×(2ème
ligne) à la ième ligne (i = 3,...,n). En
tout le processus est effectué (n-1)
fois. La matrice A devient alors triangulaire supérieure.
3 Soit R la matrice triangulaire inférieure obtenue et c la
transformée de la matrice b. Nous avons :
|
|
|
|
|
|
|
4 Finalement la solution est obtenue par simple substitution :
|
(2.6) |
Remarque 24 Il est possible de donner une interprétation matricielle de
l'algorithme de Gauss comme nous allons le voir dans théorème 4.2. A titre
d'exmple, les équations 2.4 et qui correspondent à la
transformation de la matrice A en la matrice A(1) équivaut à la
multiplication matricielle
|
(2.7) |
de même
|
(2.8) |
Il est alors facile de voir que les matrices A(k)
et b(k) au pas k de l'algorithme s'écrivent :
|
(2.9) |
de même
|
(2.10) |
pour k allant de 1 à (n-1), en admettant que A(0) = A et b(0) =
b.
Exemple 25 Soit à résoudre le système d'équations
|
Les opérations successives donnent
|
La dernière équation donne x4 = 1 ; on déduit
de la troisième (et en remplaçant x4 par sa valeur) la valeur de x3
= -2. En procédant de la même
manière en obtient x2 = 1 puis x1 = 3.
Théorème 26 Soit A une matrice non singulière. Supposons que l'ensemble des pivots successifs sont non nuls (a11,a22(1),a33(2),...), alors il existe deux matrices L et U respectivement triangulaire inférieure avec des 1 sur la diagonale et triangulaire supérieure telles que
|
La matrice triangulaire supérieure U est celle obtenue par l'application de la méthode de Gauss
|
(2.11) |
On peut vérifier que la matrice L est constituée des
éléments aik(k-1)/akk(k-1) permettant d'annuler les éléments
de colonne k pour i > k. Il suffit pour cela de se
rappeler que la soustraction de ligne ou de colonnes d'une matrice équivaut à
la multiplication à gauche par une matrice de structure adéquate (cf chapitre
3). On a alors
|
(2.12) |
avec
|
(2.13) |
On peut aussi obtenir L à partir de la matrice U
|
(2.14) |
Ce théorème signifie aussi que la résolution du système
|
(2.15) |
peut se résoudre de la manière suivante en remarquant que le système s'écrit aussi
|
(2.16) |
Définissons une variable auxiliaire y telle que
|
(2.17) |
La variable y vérifie aussi
|
(2.18) |
Finalement, la résolution du système de départ s'est simplifiée en la résolution de deux systèmes triangulaires (qui se résolvent par simple substitution).
2.2.2 Méthode de Gauss à pivot partiel
La quantité aik(k)/akk(k)
est appelée pivot. Il peut arriver que
ce pivot soit nul. Alors, la méthode précédente échoue.
De plus, pour des raisons de stabilité de la méthode, on préfère à chaque étape prendre l'élément akk(k) de module le plus grand. Au pas k de l'élimination, on examine les (n-k+1) lignes restantes et on détermine la jème ligne dans l'élément ajj(k) de module le plus grand, on permute alors la kème ligne avec la jème.
2.2.3 Méthode de Gauss à pivot total
Au pas k de l'élimination, on examine l'ensemble de la matrice d'ordre (n-k+1). On détermine l'élément aij(k)
de module le plus grand, on permute alors la kème
ligne avec la ième et la kème
colonne avec la jème. L'ordre des composantes du
vecteur x est alors changé du fait de la permutation des colonnes, il ne
faudra donc pas oublier de les remettre dans l'ordre initial.
Le théorème ci-dessous donne l'existence de la matrice triangulaire supérieure à laquelle on doit aboutir.
Théorème 27 Soit A une matrice non singulière alors il existe une matrice de permutation P (PP = I) et deux matrices L et U respectivement triangulaire inférieure et triangulaire supérieure telles que
|
(2.19) |
Le système se résout alors en résolvant par substitution
amont et aval les systèmes
|
Ce théorème sera démontré et la méthode illustrée sur un exemple dans le TD 1.
2.2.4 Nombre d'opérations effectuées
Nous n'allons compter que le nombre d'opérations complexes
et coûteuses en temps de calcul, c'est-à-dire les multiplications et les divisions.
Considérons d'abord la décomposition de PA sous la forme LU. Au pas k de l'élimination, on effectue
- (n-k) divisions par aik
- (n-k)2 multiplications pour obtenir les
éléments aij(k)
Le nombre d'opérations total est donc
|
Le calcul de l'élément yi du vecteur y à partir de l'équation 2.20 nécessite (i-1) multiplications. Le nombre d'opérations nécessaires au calcul de y est donc
|
(2.23) |
De même le calcul de xi à partir de
l'équation (2.21) nécessite (n-i+1)
multiplications. Donc le calcul de x nécessite
|
(2.24) |
Le nombre d'opérations total est donc
|
Le nombre d'opérations est donc globalement de l'ordre de [2/3]n3. La méthode présentée ci-après nécessite quant à elle un nombre d'opérations de l'ordre de [1/3]n3.
2.3 Méthode de décomposition QR (Méthode De Householder)
On démontre que toute matrice carrée A peut être
décomposée de manière unique sous la
forme
|
(2.26) |
où Q est une matrice unitaire (Q-1 = QT) et R
une matrice triangulaire supérieure telle que (rii ³ 0). La résolution du système est alors
identique à celle de la décomposition LU. Il suffit d'écrire
|
(2.27) |
Définissons une variable auxiliaire y telle que
|
(2.28) |
et
|
(2.29) |
Par conséquent, la variable y vérifie aussi
|
(2.30) |
Ce système est facile à résoudre puisque la matrice R est triangulaire.
Notons que la méthode de décomposition QR est une des plus utilisées car elle est très robuste aux erreurs numériques.
2.4 Factorisation de Cholesky
Dans le cas où la
matrice A est symétrique définie positive, la décomposition LU
prend une forme particulière appelée factorisation de Cholesky. La matrice A
se décompose sous la forme
|
(2.31) |
où B est une matrice triangulaire inférieure telle
que bii > 0. Cette décomposition est unique.
Remarque 28 Il faut noter que les méthodes de décomposition LU, QR et celle de Cholesky sont généralement celles utilisées pour le calcul du déterminant de la matrice et de son inverse. Nous rappelons que le déterminant d'une matrice diagonale ou triangulaire est le produit des éléments diagonaux.
2.5 Méthodes itératives
Le principe des méthodes itératives est de trouver une suite de vecteurs x(k+1)
de la forme
|
(2.32) |
qui converge vers la solution du système
|
(2.33) |
Le problème revient donc à trouver une matrice M et
un vecteur v assurant la convergence de la suite. Différentes méthodes
existent.
2.5.1 Méthode de Jacoby
La ième équation du système 2.33 s'écrit
|
(2.34) |
En exprimant xi en fonction des
autres variables, on obtient
|
(2.35) |
On peut bien sûr faire de même avec l'ensemble des inconnues
xi .
Les quantités M et v s'obtiennent en affectant le pas (k) au membre de droite de 2.35 et le pas (k+1) à son membre de gauche. L'équation 2.35 se réécrit alors
|
(2.36) |
Une écriture matricielle de l'équation précédente donne
|
avec
|
2.5.2 Méthode de Gauss-Seidel
Pour cette méthode, on prend
|
Cette formule est issue d'une équation similaire à 2.36
dans laquelle tous les xj(k) sont remplacés
par xj(k+1) pour (j < i)
|
(2.41) |
Cette méthode permet de diminuer l'espace mémoire nécessaire et accélère la convergence de la suite.
Remarque 29 Il existe d'autres méthodes qui sont spécifiques à la structure de la matrice A que nous n'exposerons pas ici :
- la méthode de relaxation pour les matrices creuses
- les systèmes tridiagonaux
- le système de Toeplitz (Algorithme de Levinson)
2.6 Quantification de l'erreur et de précision
Les solutions obtenues numériquement (à l'aide d'un calculateur) ne sont que
des valeurs approchées de la solution exacte. Soit x la solution exacte
du système A.x = b et [(x)\tilde] la solution
approchée calculée. Nous avons donc les relations suivantes :
|
(2.43) |
où r est appelé le résidu.
Soit z l'erreur commise sur la solution : x = [(x)\tilde]+z. On montre que z vérifie la relation
|
(2.44) |
Ainsi, on peut calculer z à partir du résidu r par la méthode de gauss par exemple.
Examinons maintenant la sensibilité de la solution vis à vis de variations possibles de A et de b. Supposons donc que chacune des quantités A,x et b de l'équation 2.42 subissent respectivement des variations DA,Dx et Db. L'équation 2.42 devient
|
(2.45) |
On démontre alors que l'erreur relative sur x peut
être majorée
|
(2.46) |
Cette formule montre l'importance du nombre condition de A.
De plus ce résultat est important car il permet de connaître le nombre de
chiffres significatifs nécessaires au calcul de la solution. Supposons qu'on
utilise d chiffres significatifs, alors
|
(2.47) |
Supposons de plus que k(A)
» 10a,
alors
|
Ce qui veut dire que x ne sera précis qu'après ( d-a-1) chiffres décimaux.
Exemple 30 Le système avec A et b donnés par
|
admet pour solution
|
Considérons maintenant (avec d = 5)
|
Il admet pour solution
|
De plus
|
donc a =
4.
On calcule
|
La formule 2.46 donne dans ce cas une
bonne estimée comme le montre le calcul
|
On vérifie donc qu'aucun des chiffres de la solution
n'est précis.
2.7 Inversion de matrice
Trouver l'inverse d'une matrice A, est équivalent à la résolution du système
d'équation
|
(2.49) |
qui donne alors le vecteur x en fonction du vecteur y
sous la forme
|
(2.50) |
On peut alors envisager un algorithme qui itère sur
l'ensemble des composantes xi de x de la manière
suivante :
- On répète les deux opérations suivantes pour i
allant de 1 à n.
- On calcule xi en fonction
de yi et des xj (j > i).
A titre d'exemple au premier pas, nous aurons
|
(2.51) |
- On remplace xi par l'expression
obtenue dans les équations des yj (j > i)
On voit que le but de cet algorithme est d'échanger les
rôles de x et de y. Partant des équations du type
|
(2.52) |
on aboutit à des équations du type
|
Les aij
sont alors les éléments de la matrice A-1.
Remarque 31 Il peut arriver qu'un des pivots soit nul. Il suffit d'échanger alors xi avec un autre des yj restant. Il faudra alors se rappeler de la permutation effectuée.
Nous donnons ci-dessous le formalise de l'algorithme dans le cas général où on permute xq avec yp
|
(2.53) |
Le choix du pivot apq permet d'exprimer xq
en fonction des xi (i ¹q)
et de yp
|
(2.54) |
En remplaçant xq par l'équation 2.54 dans les expressions de yi (i ¹p) et de yp,on
obtient
|
On obtient finalement le système ci-dessous où xq
et yp ont été permutés.
|
(2.56) |
avec
|