Résolution de systèmes linéaires
Par des méthodes directes
Elimination de Gauss
Le processus d'élimination de Gauss est probablement le plus simple que l'on puisse envisager. Le principe est d'exprimer une des inconnues en fonction des autres, de substituer la valeur obtenue dans les autres équations, de répéter l'opération jusqu'à l'obtention d'une équation à une seule inconnue, ce qui permet ensuite d'en déduire l'ensemble des autres inconnues. Essayons par exemple de résoudre le système :
x - 7y = -11
2x + 3y = 12
La première équation nous permet d'écrire x = 7y - 11, ce qui donne en réinjectant cette expression dans la deuxième équation :
2(7y-11) + 3y = 12, d'où 17y = 34 et par conséquent y = 2. On remonte ensuite facilement à x = 3.
Voici une description simplifiée de la partie "élimination" de l'algorithme :
Pour i = 1 à n-1,
pour j = i+1 à n
aij
= aij / aii
fin pour
bi = bi
/ aii
pour k = i+1 à n
pour j = i+1 à n
akj = akj - aki.aij
fin pour
bk
= bk - aki.bi
fin pour
fin pour
bn = bn/ann
Il suffit ensuite de résoudre le système triangulaire supérieur ainsi obtenu. Le nombre d'opérations (multiplications et divisions) de cet algorithme est approximativement de (n3/3 - n/3) pour la matrice A et de n2 pour chaque vecteur b.
Si toutes les sous-matrices de a ne sont pas
régulières, on risque de tomber sur un pivot nul.
Il faut alors intégrer une procédure de choix du
plus grand pivot. Le pivotage partiel consiste à
échanger deux équations dans le but d'avoir le
plus grand pivot possible en valeur absolue. Le pivotage complet
consiste à procéder à un
échange de lignes ET de colonnes, mais il requiert un nombre
d'opérations de l'ordre de n3, contre
n2 pour le partiel. Il n'est donc jamais
pratiqué.
Lorsque deux systèmes linéaires assez voisins aboutissent à des résultats très différents, on parle de mauvais conditionnement. Le conditionnement est donc la mesure de la sensibilité d'un système à de faibles modifications de l'un de ses éléments (que ce soit la matrice ou le vecteur du second membre). Supposons que l'on parte de A.x = b. Puisque ||AB|| <,= ||A||.||B||, on a ||b|| <,= ||A||.||x||. Altérons légèrement b en b+db. La solution x+dx du nouveau système est telle que A.(x+dx) = b+db. En soustrayant Ax = b à cette dernière équation, on a dx = A-1.db, donc ||dx|| <,= ||A-1||.||db||. En multipliant ceci par ||b|| <,= ||A||.||x||, on obtient ||b||.||dx|| <,= ||A||.||x||.||A-1||.||db||. En divisant par ||b||.||x||, on obtient :
(||dx|| / ||x||) <,= ||A||.||A-1||.(||db|| / ||b||)
qui montre que "l'erreur" relative sur la solution x est
bornée par le produit de l'erreur relative sur le vecteur b
et de la quantité ||A||.||A-1||.
Une mesure (parmi d'autres) du conditionnement est donnée
par r = |||A|||.|||A-1|||, où |||A|||
est ici la norme spectrale (i.e la plus petite valeur propre) de A. Le
rayon spectral r(A) est la valeur absolue du rapport entre la plus
grande et la plus petite valeur propre d'une matrice A.
La règle suivante peut être retenue : sur un
calculateur à p chiffres significatifs, le nombre de
chiffres significatifs de la solution est E[p-log10 r(A)],
où E[] est la partie entière.
On comprend en tout cas aisément que plus un
système est mal conditionné, plus la solution
sera sensible aux erreurs d'arrondi. Ceci n'implique par contre
évidemment pas qu'un système bien
conditionné soit systématiquement peu sensible
aux erreurs d'arrondi ...
Système surdéterminé
Dans le cas d'un système surdéterminé, c'est-à-dire du type A.x = b avec A matrice [m,n] et m > n, on cherche à déterminer le vecteur x qui minimise la norme de A.x - b. En pré-multipliant A.x = b par At , on obtient At.A.x = At.b, où At.A est une matrice [n,n] symétrique définie positive. Il est facile de montrer que la solution de ce système est également solution du problème aux moindres carrés initial.
Décomposition LU
Si toutes les sous-matrices de A sont régulières, A est décomposable en un produit L.U, où L et U sont des matrices triangulaires inférieure et supérieure, l'une d'entre elles (on prendra U par la suite, selon l'algorithme proposé par Crout) contenant des 1 sur sa diagonale. U correspond à la matrice obtenue par le processus d'élimination de Gauss.
La forme particulière de ces matrices permet d'établir un algorithme où le stockage de U et L se fera directement "dans" A.
Pour i = 2 à n,
a1i = a1i/a11
fin pour
pour k = 2 à n-1
akk = akk
- Σj=1->k-1 akj.
ajk
pour i = k+1 à n
aik
= aik - Σj=1->k-1
aij.ajk
aki
= (aki - Σj=1->k-1
akj.aji) / akk
fin pour
fin pour
ann = ann - Σj=1->n-1
anj.ajn
Lorsque l'on souhaite résoudre un système A.x = b, la décomposition LU permet de résoudre successivement L.y = b puis U.x = y. La résolution du système linéaire pour des seconds membres différents est donc beaucoup plus rapide.
Le nombre d'opérations pour la décomposition LU d'une matrice n x n est de (n3 - n)/3 multiplications/divisions, (2n3 - 3n2 + n)/6 additions/soustractions, et les remontées/descentes coûtent n2 mult./div. et (n2 -n) add./soustr. On peut donc dire que le gros du travail est constitué par la décomposition, qui évolue en n3/3.
Si l'on souhaite inverser une matrice, il est souvent préférable de passer par une décomposition LU pour résoudre n systèmes linéaires (avec un coût total en 4n3/3) plutôt que de calculer réellement l'inverse.
Décomposition de Cholesky
Rappelons tout d'abord que toutes les sous-matrices d'une matrice symétrique définie positive le sont également, et sont donc régulières.
Le théorème de Cholesky permet de dire que toute matrice carrée A symétrique définie positive admet une décomposition unique A = L.Lt où L est une matrice triangulaire inférieure à valeurs diagonales positives.
L'algorithme peut s'écrire de la manière suivante :
a11 = a111/2
Pour i = 2 à n,
ai1 = ai1/a11
fin pour
pour k = 2 à n-1
akk = (akk
- Σj=1->k-1 akj2)1/2
pour i = k+1 à n
aik
= (aik - Σj=1->k-1
aij.akj) / akk
fin pour
fin pour
ann = (ann - Σj=1->n-1
anj2)1/2
Matrices-bande
Si on note l la demi-largeur de bande d'une matrice, on peut montrer que le nombre d'opérations nécessaires pour une décomposition LU (ou LLt) évolue comme l2n. Pour les matrices tridiagonales, qui apparaissent notamment dans la résolution d'équations différentielles par des méthodes implicites, l'algorithme de Thomas qui dérive du procédé de Gauss permet de résoudre le système en 5n-4 multiplications (et seulement 3n-2 opérations par vecteur b supplémentaire). Le pivotage ne peut être utilisé avec l'algorithme de Thomas parce qu'il détruit la tridiagonalité, mais comme les systèmes tridiagonaux sont généralement à diagonale dominante le pivotage n'est pas nécessaire.
Par des méthodes itératives
Méthodes de Jacobi et Gauss-Seidel
Dans de nombreux problèmes numériques la matrice A est dite creuse, ce qui signifie qu'elle contient une majorité de termes nuls. Dans ces cas là les méthodes directes de résolution trouvent fréquemment des adversaires coriaces (en termes de rapidité) dans la famille des méthodes itératives. Celles-ci consistent à partir d'un vecteur initial (plus ou moins bien choisi) et à construire une suite de vecteurs qui converge vers la solution du système linéaire. On décompose pour ce faire A sous la forme A = D - E - F, où D est la matrice diagonale formée des aii, -E la matrice triangulaire inférieure contenant les éléments sous-diagonaux de A, et -F la matrice triangulaire supérieure contenant les éléments sur-diagonaux de A.
A.x = b se réécrit alors D.x = b + (E+F).x, système diagonal, ce qui permet d'obtenir la méthode de Jacobi :
xi = (bi - Σj=1->n,j#i aij.xj) / aii
Petit exemple :
4x1 - x2
=
100
-x1 + 4x2 - x3
= 100
-x2
+ 4x3 = 100
Ce système peut se réécrire sous la forme :
4x1 = x2 + 100
4x2 = x1 + x3
+ 100
4x3 = x2 + 100
En partant d'un vecteur initial (0, 0, 0), on peut donc itérer :
4x1 = 100
x1
= 25
4x2 = 100 =>
x2 = 25
4x3 =
100 x3
= 25
puis :
4x1 = 125
x1
= 31.25
4x2 = 150 =>
x2 = 37.5
4x3 =
125 x3
= 31.25
etc ...
La méthode de Gauss-Seidel consiste à résoudre (D-E).x = b + F.x, qui est un système triangulaire inférieur. Si l'on note x(n) l'itéré de rang n, l'écriture composante par composante montre que cela revient à prendre en compte les termes de x(n+1) venant d'être calculés pour obtenir les composantes restantes. Voilà pourquoi on parle parfois de méthode des "itérations successives", par rapport au nom de "méthode des itérations simultanées" qui désigne parfois celle de Jacobi.
Les méthodes itératives
définies par une relation du type x(n+1)
= M.x(n) + p sont convergentes si le rayon
spectral r(M) = maxj |lj|,
où lj est la j-ième valeur
propre de A est strictement inférieur à 1. On
peut montrer que :
- si A est symétrique définie positive, la
méthode de Gauss-Seidel est convergente.
- si A et 2D - A sont symétriques définies
positives, la méthode de Jacobi est convergente.
- si A est à diagonale dominante, les méthodes de
Gauss-Seidel, Jacobi, et SOR (voir ci-dessous) convergent. En
revanche, si A ne peut être transformée en une
matrice à diagonale dominante, le système peut
converger pour certains vecteurs initiaux mais la convergence n'est pas
garantie quel que soit ce vecteur.
La méthode de Gauss-Seidel converge
généralement (mais pas toujours) plus vite que
celle de Jacobi.
Les erreurs d'arrondi à une itération donnée ne dépendent pas des erreurs d'arrondi à l'itération précédente, alors que pour les méthodes directes on a vu que les erreurs d'arrondi polluent progressivement la solution.
Méthodes de relaxation
En utilisant les méthodes qui viennent
d'être présentées on se rend compte que
la correction apportée au vecteur solution à
chaque itération a tendance à être
sous-estimée. En d'autres termes, le vecteur ne converge pas
"assez vite" vers la solution. De là vient l'idée
de "corriger la correction", à l'aide d'un facteur
multiplicatif
En reprenant la décomposition A = D - E - F, et si on note w
un réel non nul que l'on appellera paramètre de
relaxation, on peut écrire
(D/w - E).x(n+1) = ((1-w).D/w + F).x(n) + b
qui nécessite comme Gauss-Seidel la résolution d'un système triangulaire à chaque itération. On parle de sous- ou de surrelaxation selon que w<1 ou w>1. On peut montrer que cette méthode converge pour 0 < w < 2, et que le paramètre optimal est donné par
wopt = 2 / (1 + sqrt(1-r2(J)))
où r(J) est le rayon spectral de la matrice J = D-1(E+F) de la méthode de Jacobi.
La surrelaxation successive symétrique (SSOR) consiste à faire jouer le même rôle aux matrices E et F, en introduisant un vecteur intermédiaire y(n) :
(D/w - E).y(n) =
((1-w).D/w + F).x(n) + b
(D/w - F).x(n+1) = ((1-w).D/w + E).y(n)
+ b
Méthodes du gradient et du gradient
conjugué
On définit l'application L(x) = ½.xt.A.x
- bt.x. Si A est symétrique
définie positive et si x est solution de A.x = b, alors x
minimise L.
Si xn est le n-ième
itéré, on peut chercher
l'itéré suivant sous la forme
xn+1 = xn + qn+1.wn+1
où qn+1 est un réel est wn+1 un vecteur donné. On cherche le réel qn+1 qui permet de minimiser L(xn+1). On peut montrer qu'il est donné par
qn+1 = (wn+1)t.(b - A.xn) / [(wn+1)t.A.wn+1]
Il faut en outre savoir quelle direction de descente w choisir. On peut choisir la direction de plus grande pente, qui correspond au gradient de L. On constate de suite que grad(L(xn)) = A.xn - b. Le vecteur w est donc égal au résidu. Si A est symétrique définie positive, cette méthode converge en théorie. Mais la convergence peut être très lente, et elle peut être perdue à cause des erreurs d'arrondis.
La solution consiste à corriger la direction de descente en intégrant une composante fonction de la direction à l'itération précédente :
wn+1 = -rn + bnwn
où bn est un réel calculé de manière à minimiser ||x - xn||. On montre que
bn = (rn)t.A.wn / [(wn)t.A.wn]
Le résidu ainsi obtenu à l'itération n est orthogonal aux n résidus précédents, ce qui permet de prouver que si A est une matrice NxN symétrique définie positive, la méthode du gradient conjugué converge en au plus N itérations.
Cette méthode est cependant fréquemment mal conditionnée, la vitesse de convergence dépendant directement du rayon spectral de A. On prémultiplie donc le système par l'inverse d'une matrice C, de façon à ce que le rayon spectral de C-1A soit plus petit que celui de A. C-1A n'étant pas symétrique, il est nécessaire de procéder à une décomposition de Cholesky de C, ce qui permet d'écrire
L-tL-1A.x = L-tL-1.b
puis
L-1AL-tLt.x = L-1.b
ce qui revient à un système A*.x* = b*, en notant L-1AL-t = A*, Lt.x = x*, et L-1.b = b*, et où L-1AL-t est symétrique définie positive. C'est le choix de la matrice de préconditionnement C qui va ... conditionner la rapidité de convergence de l'algorithme.
Il existe également des méthodes proches
du gradient conjugué et qui soient adaptées aux
matrices non symétriques : cf. GMRES et BICG-Stab
Critères de terminaison d'une méthode
itérative
L'éventuelle convergence d'une méthode
se juge à la variation entre deux
itérés successifs. On peut formuler un
critère en termes de variation relative ou de variation
absolue. A moins de connaître l'ordre de grandeur
approximatif de la solution, le premier choix est bien
évidemment préférable.