Résolution des systèmes linéaires

Jean-Marc Blanc

Ce cours introduit l'ensemble des techniques classiques de résolution numérique de systèmes linéaires : méthodes itératives (Jacobi, Gauss-Seidel, SOR), méthodes directes orientées mathématiques (élimination de Gauss, Cholesky…), méthode pour certains cas pathologiques (SVD, QR…)

Partagez votre avis : Commentez

- Titre : Résolution des systèmes linéaires

- Auteur : Jean-Marc Blanc

- Parution : 19 juillet 2008

- Mise à jour : 8 mars 2013

- Licence : Copyright ® 2008-2013 Jean-Marc Blanc. Aucune reproduction, même partielle, ne peut être faite de ce site et de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.

Sommaire

- I. Introduction

- II. Méthodes directes

- II-A. Systèmes avec matrice triangulaire

- II-B. Systèmes avec matrice orthogonale

- II-B-1. Vecteurs orthogonaux en algèbre linéaire

- II-B-2. Matrices orthogonales

- II-C. Méthode d'élimination de Gauss

- II-D. Systèmes mal conditionnés

- II-E. Échange de pivots

- II-F. Factorisation LU (méthode de Crout)

- II-G. Factorisation LU avec échange de pivots

- II-H. Raffinement itératif de la solution

- II-I. Inversion d'une matrice

- II-J. Calcul d'un déterminant

- II-K. Matrices symétriques et matrices hermitiennes

- II-L. Méthode de Cholesky

- II-M. Factorisation QR (méthode de Householder)

- II-N. Factorisation SVD (valeurs singulières)

- III. Problèmes pathologiques

- III-A. Matrice singulière ; système indéterminé

- Méthode SVD

- III-B. Matrice singulière ; système impossible

- Méthode SVD

- III-C. Matrice carrée mal conditionnée

- Méthode SVD

- III-D. Matrices rectangulaires : systèmes sous-déterminés

- Méthode QR

- Méthode SVD

- III-E. Matrices rectangulaires : systèmes surdéterminés

- Méthode QR

- Méthode SVD

- IV. Méthodes itératives

- IV-A. Méthode de Jacobi

- Algorithme

- Erreur

- Convergence

- Résidu

- IV-B. Méthode de Gauss-Seidel

- Algorithme

- Variante

- Erreur

- Convergence

- Résidu

- IV-C. Méthode de surrelaxation (SOR)

- Variante

- Erreur

- Convergence

- Résidu

- Méthode de surrelaxation symétrique (SSOR)

I - Introduction

De nombreux problèmes numériques se ramènent d'une manière plus ou moins approchée à la résolution de systèmes d'équations linéaires. Un système linéaire d'ordre 4 :

peut aussi être écrit sous forme matricielle :

ou en notation abrégée :

La matrice A et le vecteur B sont connus. La résolution du système consiste à calculer le vecteur des inconnues x . La matrice A et les vecteurs B et x peuvent être réels ou complexes.

La méthode de Cramer a été fréquemment enseignée dans les cours de mathématiques élémentaires. La valeur de chaque inconnue est exprimée explicitement comme le quotient de deux déterminants. On calcule tout d'abord le déterminant de la matrice A :

On remplace ensuite successivement chaque colonne de la matrice A par le vecteur B et on calcule les déterminants correspondants. On en déduit les valeurs des inconnues :

Une telle méthode utilisant des déterminants permet de faire une théorie générale des systèmes linéaires et de discuter l'existence et l'unicité de leurs solutions. En revanche, elle est totalement inadaptée pour la résolution numérique des gros systèmes linéaires, car le nombre d'opérations arithmétiques élémentaires et le temps de calcul augmentent comme la factorielle de l'ordre du système. Avec la méthode de Cramer, la résolution d'un système de 100 équations à 100 inconnues sur un super-ordinateur nécessiterait un temps de calcul supérieur à l'âge de l'univers, alors que le même travail peut être effectué en moins d'une seconde en utilisant une « bonne méthode ».

Les méthodes utilisables pour la résolution numérique de gros systèmes linéaires sont subdivisées en deux catégories :

- les méthodes directes où la solution est obtenue par un nombre d'opérations élémentaires fini et connu d'avance (le plus souvent un polynôme du 3ème degré en l'ordre n du système): méthodes de substitution, méthode d'élimination de Gauss, factorisation LR par la méthode de Crout, factorisation QR, par la méthode de Householder, méthode de Cholesky, décomposition selon les valeurs singulières, etc. ;

- les méthodes itératives à l'aide desquelles on tend par approximations successives vers la solution: méthodes de Jacobi, de Gauss-Seidel, de sur-relaxation, des gradients conjugués, etc.

La bonne stratégie consiste à choisir les algorithmes en fonction des particularités éventuelles de la matrice du système à résoudre. Les critères de choix seront, entre autres :

- la place occupée en mémoire ;

- le temps de calcul ;

- l'exactitude des résultats.

Les méthodes générales donnent, dans l'ensemble, des résultats satisfaisants, mais l'effort de mise en œuvre d'algorithmes plus subtils est presque toujours récompensé.

II - Méthodes directes

II-A - Systèmes avec matrice triangulaire

Il peut arriver que la matrice du système soit triangulaire, c'est-à-dire que tous les termes situés d'un côté de la diagonale principale soient nuls. Selon le côté, on parlera de matrice triangulaire inférieure ( L , en anglais: lower ou left ) :

ou de matrice triangulaire supérieure ( U , en anglais: upper , ou R , en anglais: right ) :

Lorsque la matrice d'un système linéaire est triangulaire, les calculs pour sa résolution sont grandement simplifiés. Prenons tout d'abord le cas d'une matrice triangulaire inférieure gauche. Le système s'écrit :

Le calcul de x_1 est immédiat :

La valeur de x_1 étant maintenant connue, on peut l'introduire dans les trois équations restantes et passer à droite les termes connus. :

C'est maintenant le calcul de x_2 qui est immédiat :

La valeur de x_2 étant calculée, on peut l'introduire dans les deux équations restantes et passer à droite les termes connus :

On a donc, de manière analogue :

puis :

et enfin :

Il va sans dire que cet algorithme ne peut fonctionner que si la matrice n'est pas singulière, c'est-à-dire si aucun des termes diagonaux n'est nul. Dans le cas exposé ci-dessus, les coefficients situés sur la diagonale principale sont quelconques. Dans d'autres cas, et en particulier lors de l'utilisation de la méthode de Crout, que nous verrons ultérieurement, on sait à l'avance que tous ces coefficients sont obligatoirement égaux à l'unité. Il est par conséquent inutile de les stocker, de tester s'ils sont nuls et de diviser par leur valeur.

Dans le cas où la matrice est triangulaire supérieure droite, et non plus inférieure gauche, on procède de la même manière, mais en commençant par la dernière équation et la dernière inconnue.

Le fait que la matrice d'un système linéaire soit triangulaire est une circonstance très favorable pour sa résolution : le nombre d'opérations élémentaires et le temps de calcul augmentent comme un polynôme du 2ème degré en l'ordre n du système). Les méthodes de substitution sont donc très rapides.

Malheureusement, les problèmes qui conduisent à des systèmes aussi simples à résoudre sont plutôt rares. Un certain nombre de méthodes « directes » ont pour but de transformer un système dont la matrice est quelconque en un système qui admet les mêmes solutions et dont la matrice est triangulaire.

II-B - Systèmes avec matrice orthogonale

II-B-1 - Vecteurs orthogonaux en algèbre linéaire

En algèbre linéaire, on fait la distinction entre les vecteurs-lignes qui sont des matrices à une seule ligne :

et les vecteurs-colonnes qui sont des matrices à une seule colonne :

La transposée d'un vecteur-ligne est un vecteur-colonne :

tandis que la transposée d'un vecteur-colonne est un vecteur-ligne :

En vertu des règles qui régissent le produit matriciel, il n'est pas possible de multiplier deux vecteurs-lignes ou deux vecteurs-colonnes. En revanche, on peut définir le produit scalaire ou produit intérieur :

et le produit extérieur :

On appelle carré scalaire du vecteur-ligne A la grandeur :

et du vecteur-colonne B la grandeur :

En vertu du théorème de Pythagore, la norme euclidienne d'un vecteur-ligne ou d'un vecteur-colonne est égale à la racine carrée de son carré scalaire :

On dit que deux vecteurs-lignes A et C de normes non nulles sont orthogonaux lorsque leur produit scalaire est nul :

De manière analogue, on dit que deux vecteurs-colonnes B et D de normes non nulles sont orthogonaux lorsque leur produit scalaire est nul :

II-B-2 - Matrices orthogonales

Une matrice réelle Q est dite orthogonale si et seulement si tous ses vecteurs-lignes sont orthogonaux entre eux et de norme unité. Il en est de même pour les vecteurs-colonnes. De cette définition il résulte que le produit d'une matrice orthogonale et de sa transposée donne la matrice unité 1 :

et, par conséquent, que l'inverse d'une matrice orthogonale est identique à sa transposée :

Le fait que la matrice d'un système linéaire soit orthogonale est une circonstance très favorable pour sa résolution puisque celle-ci se résume au produit de la transposée de la matrice et du vecteur de second membre.

Le nombre d'opérations élémentaires et le temps de calcul augmentent comme un polynôme du 2 ème degré en l'ordre n du système.

II-C - Méthode d'élimination de Gauss

La méthode d'élimination de Gauss utilise la propriété suivante :

- A toute équation d'un système linéaire on peut ajouter une combinaison linéaire quelconque des autres équations sans que cela ne change les valeurs des solutions.

En utilisant de manière judicieuse cette propriété, il est possible d'annuler successivement tous les termes situés sous la diagonale de la matrice, et ainsi de la transformer par étapes en une matrice triangulaire supérieure droite.

Soit le système

On soustrait de chaque terme de la deuxième équation le terme correspondant de la première, préalablement multiplié par A_21/A_11 . Le terme qui figure au dénominateur, ici A_11 , est appelé pivot . Il va sans dire que la division ne peut être effectuée que si le pivot n'est pas nul. Les divers termes de la deuxième équation deviennent :

On procède de la même manière avec la troisième équation :

On procède encore de même avec la quatrième équation :

Après ces opérations, le système a pris la forme :

Il s'agit maintenant d'éliminer les coefficients A_32 et A_42 . Pour ce faire, on soustrait de chaque terme de la troisième équation le terme correspondant multiplié par A_32/A_22 . Les divers termes de la troisième équation deviennent :

On procède de la même manière avec la quatrième équation :

Après ces opérations, le système a pris la forme :

Il reste encore à éliminer le coefficient A_43. Pour ce faire, on soustrait de chaque terme de la quatrième équation le terme correspondant multiplié par A_43/A_33. Les divers termes de la quatrième équation deviennent :

Le système a maintenant pris la forme voulue : sa matrice est triangulaire supérieure droite.

Il peut donc être résolu par la méthode de substitution vue plus haut. Le nombre d'opérations élémentaires et le temps de calcul augmentent comme un polynôme du 3ème degré en l'ordre n du système.

II-D - Systèmes mal conditionnés

Lorsqu'on étudie la résolution des systèmes linéaires d'un point de vue purement théorique, on distingue le cas où la matrice est régulière, qui conduit à une solution unique, et le cas où la matrice est singulière, d'où aucune solution ou une infinité de solutions. Lorsqu'il s'agit de résolution numérique, un cas supplémentaire doit être envisagé: la matrice peut être « presque singulière ».

Chacune des équations d'un système linéaire d'ordre n définit un hyperplan à n-1 dimensions dans un espace à n dimensions. La solution est représentée par les coordonnées de l'intersection de ces n hyperplans. Dans le cas d'un système d'ordre 2, il s'agit de l'intersection de 2 droites dans le plan.

Considérons tout d'abord une matrice régulière et bien conditionnée, comme par exemple celle du système :

Les deux droites se coupent de manière franche et on a une solution bien déterminée :

Fig 1.1 - Système bien conditionné

Si la matrice est singulière, cela implique que les deux droites ont la même pente. Pour le système :

elles sont parallèles. Le système n'admet donc aucune solution.

Fig 1.2 - Système impossible

En revanche, pour le système ayant la même matrice mais des seconds membres différents :

elles sont confondues. Le système admet donc une infinité de solutions.

Fig 1.3 - Système indéterminé

Considérons maintenant deux systèmes ayant la même matrice et des seconds membres très peu différents :

Leurs matrices sont presque singulières. Les équations définissent donc des droites qui se coupent sous un angle très aigu.

On voit que de très petites différences dans les seconds-membres conduisent à des solutions très différentes, puisqu'on a pour l'un des systèmes x_1 = -0,500 et x_2 = 1,600 et pour l'autre x_1 = 3,500 et x_2 = 0,000

On voit qu'un tel système est très sensible à de petites variations du second-membre. On dit que sa matrice est mal conditionnée . Le plus souvent, la cause de ce mauvais conditionnement réside dans une maladresse lors de la mise en équations du problème physique. S'il n'est pas possible de l'éviter, parce que le problème est intrinsèquement instable, il faudra utiliser une méthode de résolution plus stable, mais plus coûteuse que la méthode de Gauss.

On définira plus loin un estimateur du conditionnement Rcond qui permet d'exprimer de manière chiffrée si une matrice est plus ou moins mal conditionnée. Cet estimateur vaut 1 si tous les hyperplans correspondant aux équations sont perpendiculaires. Il est nul si la matrice est singulière. S'il est du même ordre de grandeur que la précision de calcul, on doit considérer la matrice comme tellement mal conditionnée que les résultats perdront, dans la plupart des cas, toute signification.

Fig 1.4 - Systèmes mal conditionnés

II-E - Échange de pivots

L'ordre dans lequel on range les équations et celui dans lequel on numérote les inconnues sont en principe indifférents. Deux raisons peuvent nous conduire à permuter des équations ou des inconnues :

- Nous avons vu que le mécanisme d'élimination faisait intervenir une division par chaque terme diagonal, appelée aussi pivot . Cette opération n'est évidemment possible que si le pivot n'est pas nul. Si l'on tombe sur un pivot nul, il est possible de contourner la difficulté en permutant l'équation où il se trouve avec une des suivantes. Si, à la fin de l'élimination de Gauss, il reste un pivot nul lorsque plus aucune permutation n'est possible, cela signifie que la matrice est singulière ;

- On constate d'autre part que, pour la détermination de chaque nouveau coefficient, on calcule la différence de l'ancienne valeur de ce coefficient et d'un terme qui est le produit de deux coefficients divisé par le pivot. Ce terme est entaché d'une erreur d'arrondi résultant de la multiplication et de la division. Pour que l'erreur relative sur la différence soit minimale, il faut que ce terme soit le plus petit possible en comparaison de l'ancienne valeur. Cela signifie que l'ordre optimal des équations et des inconnues est celui qui conduit à diviser par des pivots de plus grande valeur absolue.

Il est en principe possible d'échanger aussi bien les équations que les inconnues. On parle de pivotage partiel lorsqu'on n'échange que les équations et de pivotage complet lorsqu'on échange aussi bien les inconnues que les équations. On peut démontrer que si la matrice du système est régulière, la méthode de Gauss avec pivotage partiel est toujours exécutable et fournit des résultats presque aussi bons qu'avec le pivotage complet. C'est pourquoi le pivotage complet, beaucoup plus lourd, n'est presque jamais utilisé. En revanche, le pivotage partiel est essentiel si la matrice du système est mal conditionnée ou si son ordre est élevé.

II-F - Factorisation LU (méthode de Crout)

Il arrive fréquemment que l'on ait à résoudre plusieurs systèmes linéaires présentant la même matrice des coefficients A, mais des seconds-membres B différents. Si tous ces vecteurs B sont connus a priori, la méthode d'élimination de Gauss permet de résoudre simultanément les divers systèmes. Dans d'autres cas, chaque second-membre est calculé à partir de la solution du système précédent. La factorisation LU permet de n'effectuer qu'une fois la partie la plus coûteuse des calculs. La méthode repose sur le fait que toute matrice régulière peut être mise sous la forme du produit d'une matrice triangulaire inférieure gauche L dont tous les termes diagonaux sont égaux à l'unité et d'une matrice triangulaire supérieure droite U.

Effectuons littéralement le produit matriciel, soit :

Ces relations permettent de déterminer successivement les termes des deux matrices triangulaires :

On remarquera que la méthode de Crout ne nécessite pas l'utilisation de tableaux auxiliaires autres que celui utilisé à l'entrée pour le stockage de la matrice originale A et à la sortie pour le stockage des deux matrices triangulaires L et U.

Sous forme matricielle condensée, le système :

devient :

On définit maintenant le vecteur auxiliaire :

Le système peut alors être résolu par deux applications successives de la méthode de substitution, tout d'abord en avant à :

puis en arrière à :

Il convient encore de remarquer que, sans qu'il n'y paraisse de manière évidente, les méthodes de Crout et de Gauss sont équivalentes: la matrice U est identique à la matrice triangulaire obtenue par la méthode d'élimination de Gauss. C'est pourquoi certaines bibliothèques de sous-programmes intitulent "élimination de Gauss" des sous-programmes qui effectuent en fait une factorisation LU.

II-G - Factorisation LU avec échange de pivots

Pour la factorisation LU, l'échange de pivots se fait de la même manière que pour la méthode d'élimination de Gauss. Les permutations de lignes qui sont effectuées doivent naturellement être mémorisées. Cela peut se faire dans une matrice de permutation P ; une telle matrice contient des termes égaux à 1, à raison d'un seul par ligne et un seul par colonne, tous les autres termes étant nuls. On a alors :

II-H - Raffinement itératif de la solution

Il est évident qu'un système linéaire ne peut pas être résolu avec une précision meilleure que celle de l'ordinateur utilisé. Malheureusement, pour de gros systèmes, il est souvent difficile d'atteindre, voire d'approcher cette précision. Avec les méthodes directes, comme la méthode d'élimination de Gauss, les erreurs d'arrondis s'accumulent, et cela d'autant plus que la matrice est plus mal conditionnée. Le raffinement itératif de la solution est une technique qui permet de retrouver la précision de l'ordinateur utilisé. Soit le système linéaire :

La méthode d'élimination de Gauss ne nous donne pas la solution exacte x mais une solution approchée :

où est l'erreur inconnue. En multipliant la matrice A par la solution approchée, on obtient un vecteur :

d'où par soustraction :

La méthode d'élimination de Gauss, appliquée à ce nouveau système ne nous donne pas la valeur exacte de l'erreur commise précédemment, mais une valeur approchée . Par soustraction, on obtient :

qui est une meilleure approximation de x . Le procédé peut être répété plusieurs fois. L'utilisation de la méthode de Crout, plutôt que de la méthode d'élimination de Gauss, est particulièrement avantageuse lorsqu'on prévoit de devoir procéder au raffinement itératif de la solution. En effet, il suffit de n'exécuter qu'une fois la factorisation, qui est la partie la plus longue de la résolution.

II-I - Inversion d'une matrice

Contrairement à ce qu'on lit parfois, il est rarement utile de calculer numériquement l'inverse d'une matrice. Mais lorsque c'est nécessaire, il est aisé de l'obtenir à l'aide des matrices L et U. Montrons tout d'abord comment on peut calculer M = L-1, on sait que :

où I est la matrice identité de mêmes dimensions que L . On montre facilement que M est aussi triangulaire inférieure gauche.

Effectuons littéralement le produit matriciel, soit :

On en tire les termes de M :

On procède de manière analogue pour V = U-1. On a pour terminer :

II-J - Calcul d'un déterminant

Avant d'entreprendre le calcul d'un déterminant, il convient de se demander à quoi cela va servir. Si c'est seulement pour s'assurer qu'une matrice est régulière, il est beaucoup plus simple et rapide d'essayer la méthode de Crout: si, malgré toutes les permutations, il subsiste un pivot nul, cela signifie que la matrice est singulière.

On peut sérieusement douter de la signification physique du déterminant d'une matrice: en effet, si, par exemple les termes d'une matrice d'ordre n sont exprimés en ohms, son déterminant sera exprimé en ohms à la puissance n . Si malgré cela la valeur du déterminant doit véritablement être trouvée, il faut éviter à tout prix la méthode consistant à sommer des produits de termes et de leurs mineurs: le temps de calcul est proportionnel à la factorielle de la taille de la matrice. On peut procéder beaucoup plus simplement si la matrice a été factorisée, sachant que :

- le déterminant du produit de deux matrices est égal au produit de leurs déterminants ;

- le déterminant d'une matrice triangulaire est égal au produit de ses termes diagonaux.

On aura donc :

d'où :

Il faut prendre garde au fait que, si la matrice est grande, son déterminant peut devenir très grand ou très petit, et provoquer un dépassement de capacité vers le haut ou vers le bas.

II-K - Matrices symétriques et matrices hermitiennes

Dans le cas des matrices réelles symétriques et des matrices complexes symétriques ou hermitiennes, les calculs peuvent être effectués d'une manière plus efficace.

Prenons le cas d'une matrice réelle symétrique A. Elle peut être transformée en le produit d'une matrice triangulaire supérieure droite U dont tous les termes diagonaux sont égaux à l'unité, de sa transposée UT et d'une matrice diagonale D.

Effectuons littéralement le produit matriciel, soit :

Ces relations permettent de déterminer successivement les termes de la matrice diagonale et de la matrice triangulaire :

Sous forme matricielle condensée, le système :

devient :

On définit les vecteurs auxiliaires :

Le système peut alors être résolu par deux applications successives de la méthode de substitution et une division par une matrice diagonale :

Notons que le pivotage doit être complet (permutation simultanée des lignes et des colonnes) si l'on désire conserver la symétrie de la matrice.

Le cas des matrices complexes symétriques et des matrices hermitiennes se traitent de manière analogue.

II-L - Méthode de Cholesky

On rencontre souvent des problèmes qui conduisent à résoudre un système linéaire dont la matrice est symétrique et définie positive (toutes ses valeurs propres sont strictement positives). On démontre qu'une telle matrice peut être mise sous la forme du produit d'une matrice triangulaire supérieure droite S et de sa transposée ST.

Effectuons littéralement le produit matriciel, soit :

Ces relations permettent de déterminer successivement les termes de la matrice triangulaire :

Sous forme matricielle condensée, le système :

devient :

On définit le vecteur auxiliaire :

Le système peut alors être résolu par deux applications successives de la méthode de substitution, tout d'abord à :

puis à :

La méthode de Cholesky est presque deux fois plus rapide que la méthode d'élimination de Gauss. La question est de savoir, de cas en cas, si on est en droit de l'utiliser. Pour s'en assurer, il suffit de vérifier, au fur et à mesure de la factorisation, que l'argument des racines carrées est toujours strictement positif. D'autre part, de nombreux problèmes physiques reviennent à chercher le minimum d'une forme quadratique , c'est-à-dire d'une expression de la forme :

Pour annuler les dérivées partielles par rapport à chacune des variables indépendantes, on doit résoudre un système linéaire dont la matrice est nécessairement symétrique définie positive. Le nombre d'opérations élémentaires et le temps de calcul augmentent comme un polynôme du 3ème degré en l'ordre n du système.

II-M - Factorisation QR (méthode de Householder)

Comme vu précédemment, une matrice Q est orthogonale si :

Il résulte de cette définition que l'inverse d'une matrice orthogonale est identique à sa transposée :

La solution d'un système linéaire de la forme :

s'obtient en multipliant le vecteur des seconds-membres B par la transposée de la matrice Q :

Le nombre d'opérations élémentaires et le temps de calcul augmentent comme un polynôme du 2ème degré en l'ordre n du système.

La méthode de Householder consiste à écrire la matrice A comme le produit d'une matrice orthogonale Q et d'une matrice triangulaire supérieure droite que nous appellerons R pour ne pas la confondre avec celle obtenue par la méthode de Crout.

La factorisation QR obéit à la relation :

La méthode de fractionnement repose sur deux faits :

- la matrice décrivant une rotation est orthogonale ;

- le produit de matrices orthogonales est encore une matrice orthogonale.

On pose tout d'abord :

où représente un angle de rotation. On peut écrire :

Effectuons le produit de droite :

que l'on peut transformer en appliquant :

Il s'ensuit que :

Comme on le voit, on a annulé le terme 21. De manière analogue, il est possible d'annuler successivement tous les termes situés au-dessous de la diagonale principale à l'aide de matrices de rotation Q_31, Q_41, Q_32, Q_42 et Q_43 . On a ainsi transformé la matrice A en une matrice triangulaire supérieure droite R . La matrice orthogonale Q est le produit des matrices de rotation élémentaires :

Sous forme matricielle condensée, le système :

devient :

On définit le vecteur auxiliaire :

Le système peut alors être résolu en résolvant tout d'abord :

puis :

Le nombre d'opérations élémentaires et le temps de calcul sont à peu près le double de ceux qu'on aurait avec la méthode d'élimination de Gauss. C'est pourquoi la factorisation QR n'est utilisée que dans des cas spéciaux.

II-N - Factorisation SVD (valeurs singulières)

Pour une matrice A quelconque, on appelle valeur singulière et vecteurs singuliers correspondants un scalaire non négatif s et deux vecteurs et tels que :

On regroupe les vecteurs singuliers et en deux matrices U et V et les valeurs singulières, par ordre décroissant s en une matrice diagonale S :

Il s'ensuit que :

Pour résoudre le système :

on calcule les matrices U , S et V . Le système devient :

On définit les vecteurs auxiliaires :

Le système peut être résolu en résolvant successivement :

puis :

et enfin :

On appelle pseudo-inverse de Moore-Penrose le produit M = VS-1UT.

La solution du système est donc :

La méthode des valeurs singulières est plus coûteuse en temps et en place en mémoire que celles que nous avons vues précédemment, mais c'est aussi la plus stable, c'est-à-dire la plus précise, même sans utiliser le raffinement itératif ; c'est pourquoi on ne l'utilise que dans les cas exceptionnels où aucune des autres méthodes n'est applicable.

III - Problèmes pathologiques

Cette section montre comment faire face à des situations dues à trois formulations d'une même loi fondamentale :

Loi de Murphy S'il y a plus d'une façon de faire quelque chose, et que l'une d'elles conduit à un désastre, alors il y aura quelqu'un pour le faire de cette façon ;

Loi de Finagle Si quelque chose de mal peut se produire, cela arrivera ;

Loi de la tartine beurrée Une tartine beurrée tombe toujours sur le côté beurré.

Dans cette perspective, nous étudierons divers cas :

- la matrice est carrée et singulière : système indéterminé ;

- la matrice est carrée et singulière : système impossible ;

- la matrice est carrée et mal conditionnée ;

- la matrice est rectangulaire : plus d'inconnues que d'équations ;

- la matrice est rectangulaire : plus d'équations que d'inconnues.

III-A - Matrice singulière ; système indéterminé

On considère le système linéaire A x = B dont la matrice A est singulière.

La factorisation LU ne fonctionne pas car on tombe inévitablement sur un pivot nul. La factorisation QR peut être réalisée, mais conduit à une matrice R dont au moins un terme diagonal est nul ; c'est alors la substitution arrière qui tombe sur une division par zéro.

Méthode SVD

La factorisation peut toujours être réalisée. Si la matrice A est singulière, alors une ou plusieurs des valeurs singulières, stockées sur la diagonale principale de la matrice S sont nulles.

On définit les vecteurs auxiliaires y et z en posant :

La matrice U étant orthogonale, le calcul du vecteur z ne pose aucun problème :

Pour déterminer chaque composante y k du vecteur y , on doit diviser la composante zk du vecteur z par la valeur singulière correspondante Sk , ce qui est évidemment impossible si celle-ci est nulle. On alors recours à l'astuce suivante : on met à zéro la composante yk qui sans cela serait infinie, résultant d'une division par zéro.

Pour k de 1 à n
Si S(k) plus grand que zéro
y(k) = z(k)/S(k)
Sinon
y(k) = 0
Fin Si
Fin Pour

La matrice V étant orthogonale, le calcul du vecteur x ne pose aucun problème :

Il peut sembler choquant de remplacer l'inverse de zéro par zéro lui-même. En fait, parmi l'infinité de solutions qui satisfont le système à résoudre, x est celle dont la norme euclidienne est la plus petite.

III-B - Matrice singulière ; système impossible

On considère le système linéaire Ax = B dont la matrice A est singulière.

La factorisation LU ne fonctionne pas car on tombe inévitablement sur un pivot nul. La factorisation QR peut être réalisée, mais conduit à une matrice R dont au moins un terme diagonal est nul; c'est alors la substitution arrière qui tombe sur une division par zéro.

Méthode SVD

La factorisation peut toujours être réalisée. Si la matrice A est singulière, alors une ou plusieurs des valeurs singulières, stockées sur la diagonale principale de la matrice S sont nulles.

On définit les vecteurs auxiliaires y et y en posant :

La matrice U étant orthogonale, le calcul du vecteur z ne pose aucun problème :

Pour déterminer chaque composante y k du vecteur y , on doit diviser la composante zk du vecteur z par la valeur singulière correspondante Sk , ce qui est évidemment impossible si celle-ci est nulle. On alors recours à l'astuce suivante: on met à zéro la composante yk qui sans cela serait infinie, résultant d'une division par zéro.

Pour k de 1 à n
Si S(k) plus grand que zéro
y(k) = z(k)/S(k)
Sinon
y(k) = 0
Fin Si
Fin Pour

La matrice V étant orthogonale, le calcul du vecteur x ne pose aucun problème :

Il peut sembler choquant de remplacer l'inverse de zéro par zéro lui-même. En fait, aucune solution ne satisfait exactement le système à résoudre, mais le vecteur trouvé x est celui qui la satisfait le moins mal, c'est-à-dire celui pour lequel la norme euclidienne de Ax - B est la plus petite.

III-C - Matrice carrée mal conditionnée

On considère le système linéaire Ax = B dont la matrice A est mal conditionnée.

Méthode SVD

La factorisation peut toujours être réalisée :

On définit les vecteurs auxiliaires y et z en posant :

Si la matrice A est mal conditionnée, alors les dernières valeurs singulières, stockées au bas de la diagonale principale de la matrice S sont très petites par rapport à la plus grande. La matrice U étant orthogonale, le calcul du vecteur z ne pose aucun problème :

Pour déterminer chaque composante yk du vecteur y , on doit diviser la composante z k du vecteur z par la valeur singulière correspondante Sk , ce qui donnerait un résultat démesuré si celle-ci est très petite. On alors recours à l'astuce suivante : on met à zéro la composante yk qui sans cela serait démesurée, résultant d'une division par une valeur singulière très petite.

Pour k de 1 à n
Si S(k) plus grand que epsilon*S(1)
y(k) = z(k)/S(k)
Sinon
y(k) = 0
Fin Si
Fin Pour

La matrice V étant orthogonale, le calcul du vecteur x ne pose aucun problème :

Il peut sembler choquant de remplacer l'inverse d'un nombre très petit par zéro. En fait, c'est cette manière de faire qui donne les résultats les plus raisonnables.

III-D - Matrices rectangulaires : systèmes sous-déterminés

On considère le système linéaire A x = B dont la matrice A comporte m lignes et n colonnes, avec m<n :

Certaines des méthodes vue jusqu'ici ne s'appliquent, du fait de leur principe même, qu'à des matrices carrées (substitution, élimination de Gauss, factorisation LU, méthode de Cholesky). D'autres, en revanche, peuvent aussi être utilisées pour des matrices rectangulaires.

Méthode QR

La factorisation QR donne une matrice orthogonale carrée Q à m lignes et m colonnes, et une matrice « trapézoïdale » R à m lignes et n colonnes :

On définit le vecteur auxiliaire y en posant :

La matrice Q étant orthogonale, le calcul du vecteur y ne pose aucun problème :

En revanche, la matrice R ayant moins de lignes que de colonnes, la solution du système Rx = y est indéterminée. On contourne la difficulté en attribuant une valeur arbitraire (souvent zéro) à une (ou éventuellement plusieurs) inconnue(s), en général la(les) dernière(s), puis on calcule les autres par la méthode de substitution.

Méthode SVD

La factorisation SVD donne :

On définit les vecteurs auxiliaires y et z en posant :

III-E - Matrices rectangulaires : systèmes surdéterminés

Dans les systèmes surdéterminés, on a plus d'équations que d'inconnues. La matrice A a donc un nombre de lignes m supérieur au nombre de colonnes n :

Méthode QR

La factorisation QR donne :

La matrice Q est carrée; elle a m lignes et m colonnes. La matrice R est rectangulaire; elle a m lignes et n colonnes; tous les termes de la dernière ligne sont nuls.

Méthode SVD

La factorisation SVD donne :

La matrice U est rectangulaire; elle a m lignes et n colonnes. La matrice S est diagonale; elle comporte n termes La matrice V est carrée; elle a n lignes et n colonnes.

IV - Méthodes itératives

IV-A - Méthode de Jacobi

La méthode de Jacobi est une méthode itérative de résolution de système linéaire de la forme A x = B où :

Pour cela, on va construire une suite de vecteurs :

qui converge vers x , solution du système d'équations linéaires.

Remarque : chacun des vecteurs successifs est identifié par un numéro placé en exposant et entre parenthèses.

Algorithme

Un vecteur initial x(0) étant donné, l'algorithme suivant permet de déterminer les éléments successifs de la suite.

On décompose la matrice A en trois matrices L , D et U . La matrice L est constituée des termes qui se trouvent au-dessous de la diagonale principale de A(j < i) ; la matrice D contient les termes diagonaux de A(j = i) ; la matrice U est constituée des termes qui se trouvent au-dessus de la diagonale principale de A(j > i) .

Le système à résoudre peut alors s'écrire :

d'où l'on tire la formule de récurrence :

qui permet de calculer x(k+1) lorsque x(k) est connu :

On remarquera que toutes les composantes de x(k) sont utilisées pour le calcul de chaque composante de x(k+1) . Ces deux vecteurs doivent donc être stockés dans deux tableaux distincts.

Erreur

A chaque itération, le vecteur trouvé x(k+1) comporte une certaine erreur :

On pose P = D -1 (L+U) . Il vient alors :

Convergence

L'algorithme converge si ou, ce qui revient au même, .

Théorème : une condition nécessaire et suffisante pour que est que les modules de toutes les valeurs propres de P soient strictement inférieures à 1.

Théorème : la formule de récurrence converge, quel que soit x(0) , si la matrice A est à diagonale dominante, c'est-à-dire si la valeur absolue de chaque terme diagonal est supérieure à la somme des valeurs absolues des termes rectangles placés sur la même ligne.

Résidu

On appelle résidu le vecteur R (k) = B - A x (k) . La précision exigée étant donnée, on arrête les itérations lorsque :

IV-B - Méthode de Gauss-Seidel

La méthode de Gauss-Seidel est une méthode itérative de résolution de système linéaire de la forme A x = B où :

Pour cela, on va construire une suite de vecteurs :

qui converge vers x , solution du système d'équations linéaires.

Remarque : chacun des vecteurs successifs est identifié par un numéro placé en exposant et entre parenthèses.

Algorithme

Un vecteur initial x(0) étant donné, l'algorithme suivant permet de déterminer les éléments successifs de la suite.

On décompose la matrice A en trois matrices L , D et U . La matrice L est constituée des termes qui se trouvent au-dessous de la diagonale principale de A(j < i) ; la matrice D contient les termes diagonaux de A(j = i) ; la matrice U est constituée des termes qui se trouvent au-dessus de la diagonale principale de A(j > i) .

Le système à résoudre peut alors s'écrire :

d'où l'on tire la formule de récurrence :

qui permet de calculer les composantes de x(k+1) lorsque celles de x sont connues :

On remarquera que chaque composante de x(k) n'est utilisée que jusqu'au calcul de la composante correspondante de x(k+1) . Ces deux vecteurs peuvent donc être stockés dans le même tableau.

Variante

Il est aussi possible de calculer les xi(k+1) à partir du dernier. Cela revient à permuter le rôle des matrices L et U . La formule de récurrence devient :

qui permet de calculer les composantes de x(k+1) lorsque celles de x sont connues :

Erreur

A chaque itération, le vecteur trouvé x(k+1) comporte une certaine erreur :

On pose P = D +L -1U . Il vient alors :

Convergence

L'algorithme converge si ou, ce qui revient au même, .

Théorème : une condition nécessaire et suffisante pour que est que les modules de toutes les valeurs propres de P soient strictement inférieures à 1.

Théorème : la formule de récurrence converge, quel que soit x(0), si la matrice A est à diagonale dominante, c'est-à-dire si la valeur absolue de chaque terme diagonal est supérieure à la somme des valeurs absolues des termes rectangles placés sur la même ligne.

Résidu

On appelle résidu le vecteur R (k) = B - A x (k) . La précision exigée étant donnée, on arrête les itérations lorsque :

IV-C - Méthode de surrelaxation (SOR)

En utilisant la méthode de Gauss-Seidel, on observe qu'à chaque itération la correction apportée au vecteur solution a tendance à être sous-estimée. En d'autres termes, le vecteur converge trop lentement vers la solution. D'où l'idée d'augmenter la correction, à l'aide d'un facteur multiplicatif , appelé paramètre de relaxation .

Comme dans la méthode de Gauss-Seidel, on décompose la matrice A en trois matrices L , D et U :

Le système à résoudre peut alors s'écrire :

d'où l'on tire la formule de récurrence :

qui permet de calculer les composantes de x(k+1) lorsque celles de x sont connues :

Remarque : si on pose = 1, on retrouve la méthode de Gauss-Seidel

Variante

Il est aussi possible de calculer les xi(k+1) à partir du dernier. Cela revient à permuter le rôle des matrices L et U . La formule de récurrence devient :

Ce qui permet de calculer les composantes de x(k+1) lorsque celles de x sont connues :

Erreur

A chaque itération, le vecteur trouvé x(k+1) comporte une certaine erreur :

On pose P = D +L -1U . Il vient alors :

Convergence

L'algorithme converge si ou, ce qui revient au même, .

Théorème : une condition nécessaire et suffisante pour que est que les modules de toutes les valeurs propres de P soient strictement inférieures à 1.

Théorème : la formule de récurrence converge, quel que soit x(0), si la matrice A est à diagonale dominante, c'est-à-dire si la valeur absolue de chaque terme diagonal est supérieure à la somme des valeurs absolues des termes rectangles placés sur la même ligne.

Résidu

On appelle résidu le vecteur R (k) = B - A x (k) . La précision exigée étant donnée, on arrête les itérations lorsque :

Méthode de surrelaxation symétrique (SSOR)

La surrelaxation successive symétrique (SSOR) est une variante qui consiste à faire jouer le même rôle aux matrices L et U , en alternant une itération dans laquelle on commence par la première composante du vecteur x et une dans laquelle on commence par la dernière. On a donc une paire de formules de récurrence :

qui permettent de calculer les composantes de x(k+1) et x(k+2) lorsque celles de x sont connues :