Résolution de systèmes linéaires

Résumé du document

Article expliquant quatre méthodes de résolution de systèmes linéaires du type Ax=b. Le pseudo-code est fourni pour une programmation rapide.

Introduction:

Le but est de résoudre dans R des systèmes linéaires de type Ax = b où A est une matrice carrée symétrique définie positive de dimension N × N (connue), x (inconnu) et b (connu) deux vecteurs de dimension N. Quatre méthodes sont étudiées :

- gradient simple à pas optimal ;
- gradient conjugué non préconditionné ;
- gradient conjugué préconditionné diagonal ;
- factorisation LDLt.

Extraits

[...] Influence du point de départ x0 Convergence du gradient simple : méthode auto-corrective car l'énergie diminue à chaque itération. Cette méthode est très sensible au point de départ : plus x0 est proche de la solution, plus la convergence est rapide donc moins on a besoin d'itérations. Cependant, la correction rk est de plus en plus petite au fur et à mesure que l'on se rapproche de la solution ainsi l'efficacité del'algorithme diminue près de la solution. Avec la précision ǫ = les résultats trouvés ne sont pas exacts mais approchés. [...]

[...] Le "coût" de cette méthode est de 2N 2 opérations à chaque itération : - produit scalaire : N multiplications et N 1 additions ; - combinaison linéaire x ρr : N multiplications et N additions ; - produit matrice-vecteur : N 2 multiplications et N additions Algorithme du gradient conjugué Dans la méthode du gradient simple, la convergence ralentit lorsque l'on se rapproche de la solution. Ceci est dû au fait que xk+1 = xk ρk rk où rk = Axk b (résidu d'ordre k). On dit que xk est de moins en moins "corrigé". L'idée du gradient conjugué est de rendre les directions de descente consécutives "suffismment différentes" mais complémentaires. Pour cela, on découple drésidu er direction de descente. [...]

[...] Il existe d'autres préconditionnements que nous n'aborderons pas ici Conclusion Si la matrice est creuse, autant utiliser le gradient conjugué non préconditionné. Dans le cas contraire, il est préférable d'utiliser un préconditionnement (diagonal par exemple) avec le vecteur d'initialisation nul afin d'avoir moins d'itérations et un résultat plus probable. Une matrice mal conditionnée demande beaucoup plus d'itérations qu'une matrice creuse. Il existe bien sûr d'autres conditionnements comme : - SSOR (Symmetric Successive Over-Relaxation) on décompose A en D E E t avec alors C = D : diagonale de A E : partie strictement triangulaire inférieure . [...]

[...] moins d'itérations nécessaires). Plus la précision ǫ est petite et plus la solution trouvée est exacte mais avec un nombre plus grand d'itérations pour le graident simple (faux pour le conjugué) Gradient conjugué préconditionné L'idée est d'améliorer le conditionnement du système en remplaçant le système initial Ax = b par C−1Ax = C b où C est la matrice de conditionnement (donnée et facilement inversible) Algorithme Initialisation : x0 donné (vecteur nul par exemple) r0 = Ax0 b (résidu) résoudre Cq 0 = r0 ω0 = q0 k k Itération k à k + 1 : ρk = (Aωk ) k k+1 = xk ρk ω k x rk+1 = rk ρk Aω k résoudre Cq k+1 = rk+1 Test d'arrêt αk+1 = (rk ,qk ) ) ω k+1 = q k+1 + αk+1 ω k rk+ ǫ2 b 2 k+1 k+1 Effectuons deux remarques importantes sur chaque itération : - stockage nécessaire de Aω k et (rk , q k ) pour un gain de temps ; 4 - on doit résoudre un système de matrice C donc pour que la résolution soit peu onéreuse, cette matrice doit être facilement inversible Préconditionnement Le préconditionnement étudié ici sera le "diagonal" qui est le préconditionnement minimal à utiliser. [...]

[...] La résolution est effectuée ainsi : pour i = 0 à N faire pour j = 0 à i faire Bi = Bi Aij Bj fin pour j fin pour i Résolution de Dz = y avec z = Lt x D étant diagonale, on a Di zi = yi et la résolution est : pour i = 0 à N faire Bi = fin pour i Résolution de Lt x = z Une telle résolution est appelée "remontée" car Lt est une matrice triangulaire supérieure et la résolution est : pour i = N à 0 faire pour j = i à N faire Bi = Bi Aji Bj fin pour j fin pour i Bi Aii À la fin de ces trois étapes, le vecteur B (qui est en fait contient le résultat attendu Conclusion Il suffit de stocker les matrices L et D une fois pour toute et de faire des appels à la descente-remontée si l'on a plusieurs applications numériques à effectuer. Ainsi, un programmeur découpera l'algorithme en deux programmes : un pour la factorisation et le second pour la descente-remontée. [...]