1 INTRODUCTION 1

Chapitre 8

RÉSOLUTION DE SYSTÈMES LINÉAIRES

1 Introduction
La résolution de systèmes d’équations linéaires est peut-être le problème le plus connu et le
mieux étudié de l’analyse numérique. Ces études s’appuient sur l’ensemble des connaissances en
algèbre linéaire. Malgré ces atouts, la résolution pratique d’un système linéaire peut poser des
difficultés, soit pour certains choix des coefficients (systèmes «mal conditionnés»), soit lorsque
le nombre d’inconnues devient grand (plus de quelques milliers), comme c’est le cas lors de la
résolution d’équations aux dérivées partielles ou dans certains domaines de l’économétrie. Il peut
aussi arriver que le système soit «surdéterminé» (plus d’équations que d’inconnues).
Quelles sont les méthodes que propose l’arsenal mathématique ? Tout d’abord une méthode que
je qualifierai de formelle : la solution du problème

Ax = b (1)

(où A est une matrice n×n, b et x des vecteurs à n coordonnées) s’écrit x = A−1 b, puis la méthode
de Cramer, qui se révèle inutilisable en pratique dès qu’il y a plus de trois inconnues, (le calcul d’un
déterminant n × n nécessite n! opérations, ce qui est prohibitif). Je rencontre ensuite un groupe
de méthodes qui procèdent par éliminations successives des inconnues : elles sont associées aux
noms de Gauss et de Gauss-Jordan et sont d’application aisée. D’autres algorithmes «factorisent»
la matrice des coefficients (A =⇒ FG) avant de pratiquer des éliminations successives qui sont
alors très rapides à cause de la forme particulière des facteurs bf F, G : ce sont les méthodes de
Crout, Doolittle ou Cholesky. En dépit des apparences, ces deux ensembles sont mathématiquement
équivalents.
Dans un esprit bien différent, je pourrai faire appel à des procédés itératifs de résolution, dont
la convergence n’est pas toujours assurée : ce sont les méthodes associées aux noms de Jacobi et
Gauss-Seidel.
La méthode des gradients conjugués (non abordée ici) est exacte, mais peut-être interrompue
en cours de calcul.
Comme vous le verrez, le calcul de l’inverse A−1 d’une matrice A est plus long et moins stable
que le calcul direct de la solution du système linéaire (1). Une remarque analogue s’applique en
algèbre. Pour résoudre 10x = 3, je peux faire 10−1 = 0, 1; 0, 1 × 3 = 0, 3 ou x = 3/10 = 0.3.
On évite donc de calculer l’inverse de la matrice des coefficients, à moins que le problème posé
ne le demande explicitement. Si tel est le cas, je calculerai A−1 en résolvant n systèmes linéaires
particuliers, dont les seconds membres sont respectivement égaux à chacun des vecteurs de base.
Dans la mise en oeuvre des méthodes d’élimination ou de factorisation, l’ordinateur ne commet
que des erreurs d’arrondi. Dans le cas d’une méthode itérative, je devrais tenir compte aussi d’une
erreur de troncation, puisque l’algorithme s’arrête au bout d’un nombre fini d’opérations.
Pour l’analyse numérique, la difficulté de résolution d’un système linéaire tient aux effets conju-
gués de sa taille et de son «conditionnement». Voici une interprétation graphique de cet anglicisme.
Soit le système linéaire :
ax + by = m,
cx + dy = n.
1 INTRODUCTION 2

√ √
Je divise la première équation par un a2 + b2 et la seconde par c2 + d2 , pour obtenir

n1x x + n1y y = d1 , n2x x + n2y y = d2 .

Les coefficients satisfaisant maintenant aux relations n2ix + n2iy = 1, i = 1, 2, les vecteurs ni de
coordonnées {nix , niy } sont unitaires.
Le système linéaire ainsi écrit admet l’interprétation géométrique suivante. Chaque équation est
celle d’une droite Di du plan (x, y). La droite Di est perpendiculaire au vecteur ni ; la distance de
la droite à l’origine est di . La solution cherchée est représentée par le point d’intersection des deux
droites. Dans les cours de mathématiques, on distingue trois cas : D1 confondue avec D2 (système
indéterminé), D1 parallèle à D2 (système impossible), et le cas normal où D1 coupe D2 .
Numériquement, il faut distinguer ici deux possibilités. Si D1 et D2 font entre elles un angle
notable, le point d’intersection est bien défini et la solution du système sera facile à trouver. Si
D1 et D2 se coupent sous un angle très faible, le système est dit «mal conditionné». Dans ce cas,
la moindre variation de l’un des di (les seconds membres) ou de l’un des nα,i (les coefficients,
α = x, y) entraînera un déplacement très important du point d’intersection. Autrement dit, la
solution est extrêmement sensible aux perturbations des coefficients ou des seconds membres. Un
déterminant très petit est une indication qu’un système est mal conditionné ; ce n’est pas une
condition nécessaire surtout lorsque le nombre de variables devient grand.
Je peux donner un aspect plus quantitatif aux considérations précédentes. J’admets qu’il existe
une norme vectorielle k x k et une norme matricielle k A k, cohérente avec la précédente (c’est-à-
dire que, pour tout x, k A k · k x k≥k Ax k). En plus du système (1), je considère un système
«perturbé» ; le second membre a subi une petite modification ∆b, ce qui a entraîné une petite
variation ∆x de la solution
A(x + ∆x) = b + ∆b
d’où ∆x = A−1 ∆b et :
k ∆x k≤k A−1 kk ∆b k .
J’en tire une majoration de la variation relative :
k ∆x k k ∆b k k ∆b k
=k A kk A−1 k = cond(A) .
kxk kbk kbk

Le nombre cond(A) = k A kk A−1 k est une mesure de la sensibilité de l’erreur relative (sur la
solution) aux variations du second membre.
Je montre maintenant que cond(A) représente aussi la sensibilité aux variations des coefficients
des premiers membres. Je considère un nouveau système perturbé

(A + δA)(x + δx) = b.

En retranchant membre à membre (1), je trouve

δx = −A−1 δA(x + δx),

ce qui implique l’inégalité

k δx k≤k A−1 kk δA k (x + δx) k,
soit encore
k δx k k δA k
= cond(A) .
k x + δx k kAk
2 MÉTHODE DE GAUSS 3

2 Méthode de Gauss
2.1 Exemple
Je vais traiter un exemple qui, mieux qu’un long discours, fera saisir l’esprit de la méthode
d’élimination de Gauss et son équivalence avec une factorisation de la matrice des coefficients. Il
s’agit du système
    
2 1 0 4 x1 2
 −4 −2 3 −7     
Ax =    x2  =  −9  = x. (2)
 4 1 −2 8   x3   2 
0 −3 −12 −1 x4 2

Je vais procéder par éliminations successives en utilisant trois théorèmes bien connus d’algèbre
linéaire. La solution du système linéaire reste inchangée si :

* on multiplie une équation (une ligne de A et la coordonnée

correspondante de b) par une constante ;

* on ajoute une équation à une autre (une ligne de A à une autre

et la coordonnée correspondante de b à une autre) ;

* on permute deux équations (deux lignes de A et les coordonnées

correspondantes de b).

Je peux traiter A et b en parallèle ou construire une matrice A0 "augmentée" dont la n+1-ième

colonne est identique à b.
Je commence par éliminer x1 . Pour cela :
je multiplie la ligne (1) par 2 et je l’ajoute à la ligne (2),
je multiplie la ligne (1) par -2 et je l’ajoute à la ligne (3).
Il vient :    
2 1 0 4 2
 0 0 3 1   
A(1) =   ; b(1) =  −5 
 0 −1 −2 0   −2 
0 −3 −12 −1 2
On remarque que le multiplicateur relatif à la ligne i est −ai,1 /a11 ; a11 est appelé le "pivot" ; il a
intérêt à n’être pas nul ! Beaucoup d’auteurs emploient la convention opposée : ils multiplient la
ligne (1) par +ai,1 /a11 et la retranchent de la ligne (i) Le résultat est donné par la matrice A(1) et
(1)
le vecteur b(1) . Je dois maintenant éliminer x2 ; problème : le pivot correspondant (a22 ) est nul !
Mon élan parait devoir être brisé près du départ. Mais, heureusement, je peux permuter les lignes
(2) et (3) ; le nouveau pivot est maintenant -1 ; j’ajoute -3 fois (2) à (4) pour trouver :
   
2 1 0 4 2
 0 −1 −2 0   
A(3) =   ; b(3) =  −2 
 0 0 3 1   5 
0 0 −6 −1 8
2 MÉTHODE DE GAUSS 4

Finalement, en utilisant le pivot 3, j’ajoute 2 fois (3) à (4) :

   
2 1 0 4 2
 0 −1 −2 0   −2 
A =
(4)
 0 0
 ; b(4) = 
 

3 1 −5 
0 0 0 1 −2
La matrice A(4) , équivalente à A, a la forme « triangulaire supérieure ». Elle rassemble les coeffi-
cients d’un système triangulaire facile à résoudre par substitution à partir de la dernière ligne :

(4) ⇒ x4 = −2,
(3) ⇒ 3x3 − 2 = −5 ⇒ x3 = −1,
(2) ⇒ −x2 − 2(−1) = −2 ⇒ x2 = 4,
(1) ⇒ 2x1 + 4 + 4(−2) = 2 ⇒ x1 = 3.
Q (4)
Je gagne en prime la valeur de det A = (−1)p i aii , où p est le nombre de permutations de
lignes, ici 1, d’où det A = 6. Si l’un des éléments diagonaux de A(4) (l’un des pivots de l’élimination
de Gauss) est nul, le déterminant est nul ; réciproquement, si det A 6= 0, je suis assurés de trouver
des pivots non nuls (quitte à permuter des lignes).

2.2 Gauss-Jordan
Pour ceux qui auraient pris goût à l’élimination, je propose une variante. Plutôt que de substi-
tuer dans le système triangulaire, je poursuis l’élimination au-dessus de la diagonale principale :
c’est la méthode de Gauss-Jordan. Pour changer, je travaille sur la matrice augmentée A0 , réunion
de A et de b. J’ajoute (2) à (1) :
 
2 0 −2 4 0
(5)  0 −1 −2 0 −2 
A0 =   0 0

3 1 −5 
0 0 0 1 −2

puis 2/3 de (3) à (1) et à (2) :

 
2 0 0 14/3 −20/3
(6)  0 −1 0 2/3 −16/3 
A0 =
 0

0 3 1 −5 
0 0 0 1 −2

Pour finir, j’ajoute -14/3 (4) à (1), -2/3 (4) à (2) et -(4) à (3) ; il vient :
 
2 0 0 0 6
0 (8)
 0 −1 0 0 −2 
A =
 0 0

3 0 −3 
0 0 0 1 −2

écriture matricielle d’un système dont la résolution est immédiate.

2 MÉTHODE DE GAUSS 5

2.3 Représentation matricielle de l’élimination

Multiplier la coordonnée 1 d’un vecteur par le coefficient m, puis ajouter le résultat à la coor-
donnée 2 revient à multiplier ce même vecteur par une matrice particulière :
    
1 0 0 0 b1 b1
 m 1 0 0   b2   mb1 + b2 
    
 0 0 1 0   b3  =  b3 
0 0 0 1 b4 b4
Plus généralement, pour ajouter m fois la ligne i à la ligne j, il suffit de disposer le nombre m à
la place de l’élément j, i de la matrice unité (i < j). Ce qui est vrai pour un vecteur l’est également
pour une matrice. La prémultiplication de A carrée quelconque par la matrice précédente revient
à multiplier successivement chaque colonne de A, donc à ajouter la première ligne de A à la
deuxième. Encore un exemple ; je veux ajouter -3 fois la 3-ième ligne à la 4-ième de la matrice
générale : -3 est en ligne 4, colonne 3 :
  
1 0 0 0 a11 a12 a13 a14
 0 1 0 0   a21 a22 a23 a24 
  
 0 0 1 0   a31 a32 a33 a34  =
0 0 −3 1 a41 a42 a43 a44
 
a11 a12 a13 a14
 a21 a22 a23 a24 
 
 a31 a32 a33 a34 
−3a31 + a41 −3a32 + a42 −3a33 + a43 −3a34 + a44

Les matrices qui viennent d’être définies (dites matrices de Frobenius) jouissent de propriétés
intéressantes. Si f (i, j, m) désigne la matrice de Frobenius contenant m en ligne i et colonne j, alors
f −1 = f (i, j, −m) : soustraire m fois la ligne i de la ligne j est bien l’opération inverse d’ajouter m
fois (i) à (j). f (i, j, m)f (i0 , j, m0 ) = f (i0 , j, m0 )f (i, j, m) contient m et m0 en colonne j, aux lignes i
et i’, ce qui représente bien le fait qu’ajouter un multiple de la ligne j aux lignes i et i0 sont des
actions indépendantes.
Dans l’exemple précédent, j’ai obtenu la matrice A(2) à partir de la matrice A(1) = A, en
utilisant deux multiplicateurs (2 et -2), ce qui peut s’écrire :

A(2) = M1 A(1) ; b(2) = M1 b(1) ,

avec  
1 0 0 0
 2 1 0 0 
M1 = 
 −2

0 1 0 
0 0 0 1
Je pose, en vue d’une utilisation future, L1 ≡ M−1
1 . Revenant à l’exemple, j’ai ensuite permuté les
lignes (2) et (3) pour faire apparaître un pivot non nul, ce que je représente par une prémultipli-
cation par la matrice de permutation P :
 
1 0 0 0
 0 0 1 0 
P=  0 1 0 0 


0 0 0 1
2 MÉTHODE DE GAUSS 6

Cette matrice est symétrique et orthogonale : P−1 = PT = P et, par conséquent, P2 = I, ce qui
traduit le fait qu’échanger deux fois les mêmes lignes revient à ne rien faire. J’ai ensuite éliminé
x2 avec l’aide implicite de la matrice M2 .
 
1 0 0 0
 0 1 0 0 
M2 =  0 0 1 0 


0 −3 0 1

Je définis au passage son inverse, noté L2 ≡ M−1 2 .

Le système linéaire s’écrit maintenant : A(3) x = b(3) , avec A(3) = M2 PA(2) et b(3) = M2 Pb(1) .
Pour éliminer x3 , j’ai utilisé M3 :
 
1 0 0 0
 0 1 0 0 
M3 =   0 0 1 0 


0 0 2 1

La dernière forme du système s’écrit : A(4) x = b(4) . Je pose encore A(3) = L3 A(4) et b(3) = L3 b(4) .
En résumé, en ce qui concerne le tableau des coefficients, j’ai

A(4) = M3 A(3) = M3 M2 PA(2) = M3 M2 PM1 A(1)

Le vecteur b a subi les mêmes transformations, résumées par :

b(4) = M3 M2 PM1 b(1) .

Ces formules résument bien la «triangularisation» du système, mais c’est la relation inverse qui
est vraiment intéressante. Pour le voir, j’introduis les notations U = A(4) et c = b(4) (le système
linéaire est donc équivalent à : Ux = c). Je calcule ensuite

A = A(1) = L1 PL2 L3 A(4) = L1 PL2 L3 U.

Si je pose L̂ ≡ L1 PL2 L3 , je mets la matrice de départ sous la forme A = L̂U : c’est un produit
de facteurs, dont l’un est triangulaire supérieur (U) ; malheureusement, à cause du facteur P, la
matrice L̂ n’a pas de forme simple. Pour obtenir un résultat plus élégant, je commence par définir
L01 ≡ PL1 PT ; cette nouvelle matrice s’obtient en permutant les lignes (2) et (3) de L1 (prémulti-
plication par P) et en permutant les colonnes (2) et (3) de la matrice obtenue (postmultiplication
par PT ; le résultat s’écrit :  
1 0 0 0
 2 1 0 0 
L01 = 
 −2 0 1 0 


0 0 0 1
Comme vous le constatez, L01 est à nouveau triangulaire inférieure : la perturbation apportée par
P a été «absorbée» dans L01 . Je définis ensuite L = L01 L2 L3 :
   
1 0 0 0 1 0 0 0 1 0 0 0
 2 1 0 0  0 1 0 0  0 1 0 0 
L= 
 −2 0 1 0   0 0 1 0   0 0 1 0 
 

0 0 0 1 0 3 0 1 0 0 −2 1
2 MÉTHODE DE GAUSS 7

Le calcul explicite donne  

1 0 0 0
 2 1 0 0 
L=
 −2

0 1 0 
0 3 −2 1
Il suffit de recopier à leurs places dans L les éléments sous la diagonale des facteurs L01 , L2 et L3 ,
comme pour les matrices élémentaires de Frobenius. Je forme ensuite le produit LU :

LU = L01 L2 L3 A(4) = L01 L2 A(3) = L01 PA(2) = PL1 A(2) = PA(1) = PA.

La matrice PA a été décomposée en un produit de deux facteurs, l’un (L, lower) triangulaire
inférieur, avec des éléments diagonaux égaux à un, l’autre (U, upper) triangulaire supérieur. Vous
pouvez aussi effectuer numériquement le produit LU pour vous convaincre du résultat (matrice de
départ sauf permutation des lignes (2) et (3)).
Pour trouver la solution de (1), il suffit maintenant de résoudre successivement deux systèmes
triangulaires. D’abord :
L0 y = Pb,
dont la solution est le vecteur y, puis :
Ux = y
dont la solution est le vecteur x cherché. Vous pourrez vérifier que y = b(4) .
D’autre part, j’ai toujours sous la main le déterminant de A :
Y
(−1)p det(A) = det(PA) = det(LU) = det(L) det(U) = uii .
i

2.4 Généralisation
L’exemple précédent se généralise sans difficulté au cas où A est une matrice carré, régulière,
n × n. On peut démontrer l’existence et l’unicité de la factorisation. Je me contente de démontrer,
par l’absurde, cette dernière propriété. Supposons qu’il existe deux décompositions d’une même
matrice régulière :

A = L1 U1 = L2 U2 ,
d’où :
X = L−1 −1
2 L1 = U2 U1

Or, l’inverse d’une matrice triangulaire inférieure (supérieure) est elle-même triangulaire inférieure
(supérieure). Le produit de deux matrices triangulaires inférieures est une matrice de même forme ;
de même, le produit de deux matrices triangulaires supérieures est encore triangulaire supérieur.
X, qui doit réunir ces deux propriétés, ne peut être que la matrice unité, ce qui entraîne L1 = L2
et U1 = U2 .

2.5 Mise en oeuvre pratique

Si la matrice A est régulière, je suis assuré de toujours pouvoir trouver un pivot non nul :
l’élimination de Gauss est toujours possible en théorie. On montre qu’elle est encore utile pour des
matrices de rang < n ou pour des matrices rectangulaires. Je n’aborderais pas ces points ici.
2 MÉTHODE DE GAUSS 8

Les calculs réels sont en fait entachés d’erreurs d’arrondi. Ceci a pour conséquence que le
résultat final dépendra fortement du choix des pivots (ou encore de l’ordre des variables éliminées).
On peut voir en gros ce qui se passe. Je suppose que les variables x1 , x2 , . . . , xk−1 ont été éliminées.
Pour éliminer l’inconnue xk , je multiplie la ligne k par aik /akk . Si le pivot akk est petit, toutes les
erreurs qui affectent les coefficients aik seront amplifiées. En conséquence, on adopte toujours la
stratégie de «recherche du pivot partiel» ; pour éliminer la variable xk et quelque soit la valeur de
|akk |, on recherche l’élément de plus grand module dans la colonne k et c’est lui que l’on prendra
comme pivot effectif, après avoir fait la permutation de lignes nécessaire.
Certains auteurs recommandent la stratégie de la recherche complète du pivot maximal («pivot
total»). Pour éliminer xk , on cherche pour i ≥ k et pour j ≥ k l’élément aij de plus grand
module. On l’amène en position de pivot par des permutations de lignes et de colonnes. Permuter
les colonnes n’est pas neutre, c’est équivalent à permuter les inconnues. Il faut donc garder trace
des ces permutations pour restaurer l’ordre initial à la fin. La programmation est un peu plus
compliquée. Il semble par ailleurs que l’analyse précise des erreurs d’arrondi ne démontre pas une
supériorité nette de la recherche globale sur la recherche partielle.
Les erreurs d’arrondi sont plus faibles si les éléments des matrices successives restent à peu près
du même ordre de grandeur. Pour cela, on peut, avant de résoudre le système, multiplier certaines
lignes par des facteurs tels que le plus grand coefficient de chaque ligne soit compris entre 0,5 et 1
(équilibrage).
Le décompte détaillé des opérations arithmétiques nécessaires fait apparaître le résultat suivant :

factorisation n(n − 1/2)(n − 1)/3 additions et multiplications

n(n − 1)/2 divisions
substitution n(n − 1)/2 additions et multiplications
n divisions

Vous retiendrez les valeurs asymptotiques : pour l’élimination : O(n3 /3) opérations et pour la
résolution du système linéaire : O(n2 ) opérations.
L’algorithme de Gauss est assez économe en mémoire, si l’on accepte de détruire A. En effet,
à chaque étape, je range les éléments de L à la place que devraient occuper les zéros de A après
élimination (dans le triangle inférieur). Les 1 de la diagonale de L sont sous-entendus. Le triangle
supérieur de A est progressivement occupé par les éléments de U.

2.6 Calcul de l’inverse de A

Pour certaines applications, il faut calculer A−1 = X explicitement, soit encore trouver X telle
que AX = I. La colonne j de X est solution du système linéaire :

Axj = ej .
Je vais avoir à résoudre n systèmes linéaires dont les premiers membres sont identiques, mais
les seconds membres différents. Il est commode d’utiliser la décomposition PA = LU effectuée une
seule fois, puis de résoudre les systèmes triangulaires Lyj = Pej puis Uxj = yj . On sait que les
coordonnées de yj sont nulles jusqu’à un certain rang, ce qui permet de gagner un peu de temps.
Le nombre d’opérations nécessaires est O(n3 ), pas plus que pour calculer A2 mais trois fois plus
que pour résoudre un système linéaire de même n.
3 FACTORISATION DIRECTE 9

3 Factorisation directe
Je sais que la matrice régulière A admet une représentation sous forme d’un produit de matrices
triangulaires LU. Je peux alors me poser la question suivante : est-il possible de trouver les éléments
de L et U sans utiliser l’algorithme de Gauss ? La réponse est affirmative et c’est même assez
simple. Cependant, comme la décomposition LU est unique, ce nouvel algorithme est certainement
équivalent à celui de Gauss.
Soient `i,j les éléments de L, avec `i,j = 0 si i < j ; de même, les éléments de U sont tels que
ui,j = 0 si i > j ; enfin `i,i = 1. Je procéderai par identification.
X
a11 = `1k uk1 = `11 u11 = u11 ,
X
a12 = `1k uk2 = u12 ,
......................
X
a1n = `1k ukn = u1n .

En utilisant la première ligne de A, j’ai trouvé la première ligne de U. L’astuce pour continuer
(due à Crout, d’où le nom de l’algorithme), consiste à identifier maintenant la première colonne de
A.
X
a21 = `2k uk1 = `21 u11 + `22 u21 ⇒ `21 = a21 /u11
X
a31 = `3k uk1 = `31 u11 + `32 u21 + `33 u31 ⇒ `31 = a31 /u11
................

Je détermine ainsi la première colonne de L. Vous pouvez imaginer que la deuxième ligne de U
s’obtient à partir de la deuxième ligne de A, puis que la deuxième colonne de L découle de la
deuxième colonne de A, etc... Les formules générales s’écrivent :
i−1
X
uij = aij − `ik ukj ; j = i, i + 1, ...n;
k=1
" j−1
#
1 X
`i,j = aij − `ik ukj ; i = j + 1, j + 2, ...n (3)
ujj
k=1
(4)

avec toujours `ii = 1. Il y a une petite difficulté : cet algorithme est très sensible aux erreurs
d’arrondis, et ne fonctionnera pas sans recherche du pivot maximal dans chaque colonne. D’autre
part, en général, je forme la décomposition LU pour résoudre un système linéaire ; ici, je n’ai
permuté que les lignes de A et non celles de b (qui n’apparaît pas) ; il faut donc garder trace des
permutations pour pouvoir retrouver les inconnues par la suite.
On rencontre fréquemment deux cas particuliers. La premier est celui d’une matrice définie
positive, le second celui d’une matrice tridiagonale. Je vais les examiner brièvement.
3 FACTORISATION DIRECTE 10

3.1 Matrice tridiagonale

Soit A une matrice tridiagonale ; Je pose ai,i−1 = xi , ai,i = yi et ai,i+1 = zi . J’admets en plus
que A possède une décomposition LU sans permutation de lignes. Par identification, on vérifie que
L et U sont bidiagonales :
 
y1 z1 0 ... ... 0
 x2 y2 z2 0 ... 0 
 
 0 ... ... ... . . . 0 =
 
 0 . . . 0 xn−1 yn−1 zn−1 
0 ... ... 0 xn yn
  
1 0 0 ... 0 d1 u1 0 ... 0
 l2 1 0 ... 0   0 
   0 d2 u2 ... 
 0 ... 1 . . . 0  0 ... ... ... 0 
  
 0 . . . ln−1 1 0   0 . . . 0 dn−1 un−1 
0 ... 0 ln 1 0 ... 0 0 dn

Vous prouverez sans peine les relations suivantes :

y1 = d1 ; ; z1 = u1 ; xi = li di−1 ;

yi = li ui−1 + di ; ; zi = ui .
Vous pourrez finalement établir les formules très simples de résolution d’un système linéaire dont
las coefficients seraient les xi , yi , zi . Dans le cas où l’on doit pratiquer une recherche du pivot
partiel, le pivot de l’étape k se trouve nécessairement dans l’une des lignes k ou k + 1 (matrice
tridiagonale). Il en résulte une décomposition PA = LU, où L est toujours bidiagonale inférieure,
alors que U est tridiagonale supérieure.

3.2 Matrice symétrique définie positive

Si A est symétrique, sa décomposition s’écrit : A = LLT . Je suppose en plus que A est
définie positive. Je rappelle deux définitions équivalentes : xT Ax > 0 quelque soit le vecteur
x ; toutes les sous-matrices principales ont des déterminants positifs). On montre alors que les
pivots de l’élimination de Gauss sont strictement positifs et donc que la factorisation LLT se
fait sans permutation (décomposition régulière de Cholesky). On détermine les éléments de L
(`ij = 0sij > i) par identification :

`211 = a11 > 0,

`11 `21 = a21 ,
`21 + `222 = a22
2

.......................
" j−1
#
1 X
`i,j = aij − `ik `jk ,
`jj
k=1
" j−1
#
1 X
`i,j = aij − `ik `jk
`jj
k=1
4 MÉTHODES ITÉRATIVES DE RÉSOLUTION DES SYSTÈMES LINÉAIRES 11

3.3 Variantes
J’ai montré que la matrice régulière A admet la factorisation LU en un produit de matrices
triangulaires, inférieure pour L à diagonale unitaire, supérieure pour U. Cette répartition des rôles
n’est pas la seule possible. On peut concevoir la décomposition A = L0 U0 où U0 est maintenant
triangulaire à diagonale unitaire, L0 étant simplement triangulaire ; un algorithme du à Doolittle
et très voisin des précédents permet d’obtenir L0 et U0 .
Plutôt que de favoriser l’un des facteurs grâce à la possession d’une diagonale unitaire, on peut
imposer que chaque facteur triangulaire présente des éléments diagonaux égaux à 1. On a alors :

A = LDU , `ii = uii = 1

où D est une matrice diagonale. Dans le cas particulier où A est symétrique, il vient :

A = LDLT

4 Méthodes itératives de résolution des systèmes linéaires

Les très grands systèmes linéaires posent un problème pratique : il peut être impossible de
caser en mémoire centrale l’ensemble des coefficients. D’autre part, on rencontre, surtout lorsque le
système linéaire provient de la discrétisation d’une équation aux dérivées partielles, une structure
particulière des coefficients : les éléments diagonaux sont plus grands que tous les autres. Dans
ces cas là, on peut utiliser avec profit des méthodes itératives de résolution, beaucoup moins
gourmandes en mémoire et plus rapides.

4.1 Méthode de Jacobi

Je cherche la solution de (1), avec aii 6= 0, 1 ≤ i ≤ n, en connaissant une approximation de
la solution : x[0] . Je vais utiliser une méthode d’approximation successive en traitant de façon
équivalente chaque coordonnée de x :

[k] 1 X [k−1]
xi = (bi − aij xj ) ; 1 ≤ i ≤ n. (5)
aii
j6=i

[k]
Le calcul de xi n’utilise que la ligne i de A : cette matrice peut donc être rangée sur disque et
lue à la demande.
Comme pour toute méthode d’itération, je dois me préoccuper d’un critère d’arrêt. On rencontre
couramment l’une des deux définitions classiques suivantes. Je définis le «résiduel» à l’itération k :

r[k] = b − Ax[k] .

Je peux choisir la condition de convergence

k r[k] k≤ ² k b k

, ² étant un seuil fixé à l’avance ou la condition :

k x[k] − x[k−1] k≤ ² k x[k−1] k .

4 MÉTHODES ITÉRATIVES DE RÉSOLUTION DES SYSTÈMES LINÉAIRES 12

L’algorithme (5) admet une représentation matricielle commode pour les analyses théoriques. Je
décompose A en une somme de trois matrices (rien à voir avec les L et U des paragraphes précé-
dents) :
A=D−L−U
où D est diagonale, E strictement triangulaire inférieure et F strictement triangulaire supérieure
(j’ai donc aii = dii ), `ii = uii = 0). J’introduis le vecteur δx[k] ≡ x[k+1] − x[k] . En utilisant les
relations de récurrence et la définition du vecteur résiduel, vous montrerez facilement que :
δx[k] = D−1 r[k]
d’où :
x[k+1] = D−1 (L + U)x[k] + D−1 b.
Exemple. Étant donné le système linéaire ci-dessous et le vecteur initial [0, 0]T :
3x + y = 2
−x + 2y = −2

ou
x = (1/3)(2 − y)
1
y = (x − 2)
2
je trouve la suite de solutions approchées suivante :
x 2/3 1 8/9 5/6 23/27 31/36 139/162
y −1 −2/3 −1/2 −5/9 −7/12 −31/54 −41/72
qui converge effectivement mais lentement vers la solution exacte [6/7, −4/7]T .

La méthode de Jacobi converge lentement, mais je peux facilement modifier cet algorithme pour
le rendre plus rapide.
4.3. Méthode de Gauss-Seidel
[k+1] [k] [k] [k] [k] [k]
L’algorithme de Jacobi me prescrit de calculer xi en utilisant les valeurs x1 , x2 , . . . , xi−1 , xi+1 , ...xn ;
[k+1] [k+1] [k+1]
il semble que je pourrais gagner du temps ou de la précision, puisque je connais déjà x1 , x2 , . . . , xi−1 ,
qui sont plus précises que les valeurs des mêmes variables à l’itération précédente (k). C’est ce que
fait l’algorithme de Gauss-Seidel. Les relations de récurrence s’écrivent :
 
i−1
X Xn
k+1 1  [k+1] [k] 
xi = bi − aij xj − aij xj (6)
aii j=1 j=i+1

Ces équations sont aussi plus faciles à programmer que celles qui découlent de la méthode de
[k]
Jacobi. En effet, il n’est pas nécessaire de conserver l’ensemble des valeurs xi pendant le calcul
[k+1]
des xi : chaque nouvelle valeur remplace la précédente.
Je vous laisse le soin de montrer que l’équivalent matriciel de ces formules est :
x[k+1] = (D − L)−1 Ux[k] + (D − L)−1 b.
Exemple. Appliquant l’algorithme de Gauss-Seidel à l’exemple du paragraphe précédent, je trouve
la suite de valeurs :
5 UN PEU DE PROGRAMMATION 13

x 2/3 8/9 23/27 139/162

y -2/3 -5/9 -31/54 -185/324
Certains de ces nombres sont apparus lors du calcul précédent ; on observe toutefois que la conver-
gence est beaucoup plus rapide.
4.4. Méthode de surrelaxation
Cette variante de l’algorithme de Gauss-Seidel (elle peut utiliser Jacobi aussi, mais sans avan-
tage) procède comme suit. Connaissant un vecteur x[k] , je calcule un vecteur provisoire x̂[k+1] à
l’aide de l’algorithme (6), puis je définis le vecteur x[k+1] définitif comme :

x[k+1] = ωx̂[k+1] + (1 − ω)x[k] . (7)

ω est un paramètre réel dont il faut optimiser la valeur pour chaque classe de problème. Le choix
ω = 1 me ramène à (6). On a toujours ω < 2, et l’expérience montre que la meilleure valeur de ω
est voisine de 1,7.

5 Un peu de programmation
Il est frustrant d’écrire un beau programme d’inversion de matrice, par exemple, et de devoir
recompiler le programme pour chaque valeur de l’ordre de la matrice ; en général, les programmeurs
tournent la difficulté en déclarant une taille de matrice très supérieure aux besoins prévisibles, d’où
un encombrement de la mémoire qui peut devenir gênant. En Pascal, il est possible de déclarer des
tableaux, arguments de sous-programmes, sans spécifier à l’avance leur taille. C’est évidemment
dangereux, car si l’on s’affranchit du strict contrôle de taille effectué par le compilateur, on travaille
sans filet. Je propose ci-dessous deux possibilités (tirées de la notice Borland). Dans ces ébauches
de programme, les instructions d’affichage des résultats ont été omises.

PROGRAM dim_var_1;

{$R-}
PROCEDURE vext(n: INTEGER; VAR y; VAR ymin, ymax, dy: REAL);
TYPE tmin = ARRAY[0..0] OF REAL;
VAR i: INTEGER;
BEGIN
ymin := tmin(y)[0]; ymax := tmin(y)[0];
FOR i := 1 TO n{\-}1 DO BEGIN
IF tmin(y)[i] < ymin THEN ymin := tmin(y)[i];
IF tmin(y)[i] $>$ ymax THEN ymax := tmin(y)[i];
END;
dy := ymax {\-} ymin;
END;
{$R+}

CONST n = 128;
TYPE tab = ARRAY[1..n] OF REAL;
VAR x, max, min, delta: REAL; vec: tab; i: INTEGER;
BEGIN
FOR i := 1 TO n DO
5 UN PEU DE PROGRAMMATION 14

vec[i] := x := 100.0*RANDOM - 50.0;

END;
WRITELN;
vext(n,vec,max,min,delta);
END.

PROGRAM dim_var_2;

{$R-}
PROCEDURE vext(n: INTEGER; VAR y; VAR ymin, ymax, dy: REAL);
TYPE tmin = ARRAY[0..0] OF REAL;
VAR i: INTEGER; ya: tmin ABSOLUTE y;
BEGIN
ymin := ya[0]; ymax := ya[0];
FOR i := 1 TO n{\-}1 DO BEGIN
IF ya[i] $<$ ymin THEN ymin := ya[i];
IF ya[i] $>$ ymax THEN ymax := ya[i];
END;
dy := ymax - ymin;
END;
{$R+}

CONST n = 128;
TYPE tab = ARRAY[1..n] OF REAL;
VAR x, max, min, delta: REAL;
vec: tab; i: INTEGER;
BEGIN
FOR i := 1 TO n DO
vec[i] := 100.0*RANDOM - 50.0;
vext(n,vec,max,min,delta);
END.

Chaque programme remplit un vecteur de nombres aléatoires, puis la procédure vext trouve
les valeurs extrêmes et leur différence. Elle est programmée pour accueillir un argument (y) de
taille quelconque ; en contrepartie, le contrôle de portée des indices est désactivé (commutateur de
compilation {$R-} ; il est réactivé à la fin de la procédure ({$R+}). Le premier exemple utilise le
transtypage pour affecter les coordonnées de y à celle de tmin (ymin := tmin(y)[i]). Le second
utilise l’instruction ABSOLUTE pour parvenir au même résultat (le premier élément de ya coïncide
avec le premier élément de y et tous les autres se correspondent deux à deux). Les deux procédés
peuvent être employés simultanément, comme ci-dessous.
PROGRAM dim_var_3;

{$R-}
PROCEDURE vext(n: INTEGER; VAR y; VAR ymin, ymax, dy: REAL);
TYPE tmin = ARRAY[0..0] OF REAL;
VAR i: INTEGER; ya: tmin ABSOLUTE y;
BEGIN
5 UN PEU DE PROGRAMMATION 15

ymin := ya[0]; ymax := tmin(y)[0];

FOR i := 1 TO n-1 DO BEGIN
IF ya[i] < ymin THEN ymin := tmin(y)[i];
IF ya[i] > ymax THEN ymax := ya[i];
END;
dy := ymax - ymin;
END;
{$R+}

CONST n = 128;
TYPE tab = ARRAY[1..n] OF REAL;
VAR x, max, min, delta: REAL;
vec: tab; i: INTEGER;
BEGIN
FOR i := 1 TO n DO
vec[i] := 100.0*RANDOM - 50.0;
vext(n,vec,max,min,delta);
END.
En Pascal une structure (tableau, enregistrement) ne peut pas excéder 64500 octets ; de plus,
le nombre de structures de grande taille utilisables simultanément est limité. Ces contraintes pro-
viennent de ce que Pascal range les variables dans une zone de mémoire appelée la «pile» («stack»),
dont la taille est de 64500 octets. On peut manipuler simultanément plusieurs grandes matrices en
utilisant le «tas» («heap»), accessible uniquement par l’intermédiaire de pointeurs. Dans l’exemple
suivant, j’emploie des pointeurs sur des tableaux de pointeurs pour effectuer un calcul bidon met-
tant en oeuvre trois matrices 100 x 100.

program grosse_m;
CONST max_l = 100; max_c = 100;
TYPE
p{_reel = ^REAL; p_mat = ^mat;
mat = ARRAY[1..max_l,1..max_c] OF p_reel;
VAR
pa,pb,pc:p_mat;
i,j,k: INTEGER;
s: REAL;
BEGIN
NEW(pa);NEW(pb);NEW(pc);
FOR i := 1 TO max_l DO BEGIN
FOR j := 1 TO max_c DO BEGIN
NEW(pa^[i,j]); pa^[i,j]^ := 1/(i+j-1);
NEW(pb^[i,j]); pb^[i,j]^ := pa^[i,j]^;
END;
WRITE(i:4);
END;
FOR i := 1 TO max_l DO BEGIN
FOR j := 1 TO max_c DO BEGIN
s := 0;
5 UN PEU DE PROGRAMMATION 16

FOR k := 1 TO max_l DO s := s + pa^[i,k]^*pb^[k,j]^;

NEW(pc^[i,j]); pc^[i,j]^ := s;
END;
WRITE(I:4);
END;
FOR i := 1 TO max_l DO BEGIN
FOR j := 1 TO max_l DO BEGIN
WRITE(pc^[i,j]^:10:4); IF ((j+1) mod 8 = 0) THEN WRITELN;
END;
WRITELN; WRITELN;
END;
END.
Cependant, un réel étant représenté par 6 octets, les matrices précédentes occupent chacune 60 000
octets, soit presque la taille maximale.