Chapitre 3

Systèmes d’équations linéaires

Mathématiques-Informatique, Sciences de l’Homme et de la

Société

Université Moulay Ismail - Meknès

2017/2018

Math-Info-ScHS (ENSAM) CMN 1 / 105

1 Introduction

2 Élimination de Gauss

3 Décomposition LU

4 Effets de l’arithmétique flottante

5 Conditionnement d’une matrice

6 Méthodes itératives pour la solution d’un système linéaire

Motivations
Les itérations de Jacobi
Les itérations de Gauss Seidel
Les itérations de relaxation
Critères de convergence
Formulation matricielle des itérations
Condition suffisante de convergence
Math-Info-ScHS (ENSAM) CMN 2 / 105
 Introduction Introduction

Après avoir résolu (2ème année) une équation de la forme

f (x) = 0 où x ∈ R et f (x) ∈ R (i.e. f : R →: R)

On complique maintenant les choses en considérant les systèmes

d’équations linéaires de n équations et n inconnues:

Ax = b

où A est une matrice carrée inversible et b un vecteur (F : Rn →: Rn

avec F (x) = Ax − b et on cherche F (x) = 0)

Math-Info-ScHS (ENSAM) CMN 3 / 105

 Introduction Introduction

Beaucoup d’applications dans le domaine de l’ingénierie exigent de

résoudre des systèmes d’équations qui comportent souvent un grand
nombre d’inconnues. On rencontre ces applications, par exemple,
dans les domaines suivants:
la mécanique des fluides: écoulement de l’air autour d’un avion,
écoulement dans une turbine hydraulique, etc...
la mécanique des solides: déformation d’un matériau soumis à
des contraintes extérieures, analyse des structures complexes,
etc....
Optimisation de la distribution des fréquences sur des antennes,
filtrage de signaux, etc....

Math-Info-ScHS (ENSAM) CMN 4 / 105

 Introduction Introduction

Comportement hygro-mécanique de lames de parquet en services

(période de 42 jours, humidité relative augmentée de 20%).

Déformation d’une lame de parquet par variation des conditions d’humidité

relative de l’air ambiant.
Ce calcul exige de résoudre un premier système linéaire de 35890 équations
pour le mouvement de l’humidité puis un système de 107670 équations pour
la déformation. Ces deux résolutions sont faites au total 129 fois pour décrire
le comportement à divers moment dans le temps.
Math-Info-ScHS (ENSAM) CMN 5 / 105
 Élimination de Gauss Généralités

Un système d’équations linéaires est un système de n équations

linéaires à n inconnues


 a11 x1 + a12 x2 + a13 x3 + . . . + a1n xn = b1
 a21 x1 + a22 x2 + a23 x3 + . . . + a2n xn = b2

.. .


 . = ..
an1 x1 + an2 x2 + an3 x3 + . . . + ann xn = bn


où aij et bj sont connus pour i, j = 1, ..., n.

Exemple: Trouver x1 , x2 , x3 solution de


 2x1 + 3x2 − x3 = 5
4x1 + 4x2 − 3x3 = 3
−2x1 + 3x2 − x3 = 1


Math-Info-ScHS (ENSAM) CMN 6 / 105

 Élimination de Gauss Généralités

Écriture matricielle:
Ax = b
où A dénote la matrice carrée d’ordre n des coefficients du système
linéaire et b le vecteur du second membre.
     
a11 ... a1i ... a1n b1 x1
 .. . . .. . . ..   ..   .. 
 .
 . . . . 

 . 
 
 . 
 
A =  ai1 ... aii ... ain 
  b =  bi 
  x = xi


 .. . . .. . . .. . . 
 ..   .. 
   
 . . . . . 
an1 ... ani ... ann bn xn

Exemple: Trouver x ∈ R3 solution de

    
2 3 −1 x1 5
 4 4 −3   x2  =  3 
−2 3 −1 x3 1

Math-Info-ScHS (ENSAM) CMN 7 / 105

 Élimination de Gauss Généralités

Rappel d’algèbre:
Trois énoncés équivalents:
Le système Ax = b admet une solution unique si et seulement si
la matrice A est de rang maximal (autant d’inconnues que
d’équations et pas d’équation qui se "répète").
Le système Ax = b admet une solution unique si et seulement si
detA 6= 0
Le système Ax = b admet une solution unique si et seulement si
A est inversible.
Dans ce cas la solution est fournie par
x = A−1 b
detAi
La formule de Cramer: xi =
detA
où Ai est la matrice A dont la colonne i est remplacée par b.

Math-Info-ScHS (ENSAM) CMN 8 / 105

 Élimination de Gauss Généralités

Remarque
La plupart du temps on traite des cas où la matrice est inversible
(matrice non-singulière). Ce qui assure l’existence d’une solution.
Sauf dans de rare cas, la formule de Cramer peut être considérée
comme purement théorique (sans applications).
Sauf dans des cas particuliers il est plus rapide et simple de
résoudre le système linéaire que d’inverser la matrice (même "à la
main").

Comment résoudre efficacement des systèmes linéaires?

Math-Info-ScHS (ENSAM) CMN 9 / 105

 Élimination de Gauss Méthode des substitutions successives

Méthode des substitutions successives

La méthode la plus simpliste (et classique) consiste a éliminer des
inconnues dans les équations par substitution successive:
    
2 3 x1 8
=
3 4 x2 11

De la première équation on a

8 − 3x2
x1 =
2
Et en substituant dans la deuxième (éliminant x1 ) on a

8 − 3x2
3( ) + 4x2 = 11 ⇒ x2 = 2
2

Math-Info-ScHS (ENSAM) CMN 10 / 105

 Élimination de Gauss Matrices diagonales

Pour une matrice quelconque cette approche est simple pour un

humain mais difficile à traduire en algorithme. On doit abandonner
cette approche du moins sous cette forme.
Cette approche est plus simple à décrire pour certaines structures de
matrices.

Matrice diagonale
Facile lorsque le système linéaire a une matrice diagonale:
 
a11 0 0 ... 0    
x1 b1
 0 a22 0 ... 0 
..  =  ..  ⇔ a x = b
Ax =  i = 1, ..., n
 
.. ..  .   .  ii i i
 . . 0 
xn bn
0 ... 0 ... ann

bi n
Q
alors xi = . Solution unique ssi aii 6= 0 ∀i (detA = aii 6= 0)
aii i=1

Math-Info-ScHS (ENSAM) CMN 11 / 105

 Élimination de Gauss Matrices triangulaires

Définition 3.1
Une matrice est dite triangulaire inférieure (ou supérieure) si toutes
ces entrées aij sont nulles pour i < j (i > j resp.).
Triangulaire inférieure = “nulle au dessus de la diagonale“.
Triangulaire supérieure = “nulle au dessous de la diagonale“.
   
a11 0 0 ... 0 a11 a12 a13 ... a1n

 a21 a22 0 ... 0 


 0 a22 a23 ... a2n 


 a31 a32 a33 ... 0 


 0 0 a33 ... a3n 

 .. .. .. .. ..   .. .. .. .. 
 . . . . .   . . . . 
an1 an2 an3 ... ann 0 0 0 ... ann

Math-Info-ScHS (ENSAM) CMN 12 / 105

 Élimination de Gauss Algorithmes de résolution pour les matrices triangulaires

Matrice triangulaire
Les matrices triangulaires: résoudre sucessivement en choisissant
l’équation avec le moins d’inconnues.
Matrice triangulaire inférieure: descente triangulaire (on
commence à 1 on fini à n ”en bas”).
i−1
X
bi − aik xk
b1 k =1
x1 = xi = i = 2, ..., n
a11 aii

Matrice triangulaire supérieure: remontée triangulaire (on

commence à n on fini à 1 ”en haut”).
n
X
bi − aik xk
bn k =i+1
xn = xi = i = n − 1, n − 2, ..., 1
ann aii
Math-Info-ScHS (ENSAM) CMN 13 / 105
 Élimination de Gauss Unicité dans le cas des matrices triangulaires

Unicité dans le cas des matrices triangulaires

La remontée et la descente triangulaire n’ont de sens qu’a la condition
que les termes diagonaux soient tous non nuls.
Que se passe-t’il si cela n’est pas le cas?
Puisque le déterminant d’une matrice triangulaire est le produit des
termes de la diagonale, alors
n
Y
detA = aii 6= 0 ⇔ aii 6= 0 ∀i
i=1

Donc l’unicité de la solution impose que les termes diagonaux soient

tous non nuls. Dans le cas contraire le système n’a pas une solution
unique et A n’est pas inversible.

Math-Info-ScHS (ENSAM) CMN 14 / 105

 Élimination de Gauss Exemple de résolution

Exemple de remontée triangulaire

Soit à résoudre le système
    
3 1 3 x1 0
 0 2 2   x2  =  30 
0 0 1 x3 15

3
X
b2 − a2k xk
b3 k =3 30 − 2x3
x3 = = 15 x2 = = =0
a33 a22 2
3
X
b1 − a1k xk
k =2 0 − 1x2 − 3x3
x1 = = = −x3 = −15
a11 3

Math-Info-ScHS (ENSAM) CMN 15 / 105

 Élimination de Gauss Comment exploiter les algorithmes des matrices triangulaires?

L’algorithme de remontée (ou de descente) est efficace. On sait que

l’on peut transformer un système quelconque en système triangulaire.
On voudrait
une méthode de passage d’un système de quelconque vers un
système triangulaire tout en conservant la solution du système
original.
utiliser la remontée ou descente triangulaire.
Pour cela on va revoir les opérations élémentaires sur les systèmes
linéaires. On produira ensuite deux méthodes directes: la méthode
de Gauss et la méthode de la décomposition LU.

Definition 3.2
Une méthode de résolution d’un système linéaire est dite directe si la
solution du système peut être obtenue en un nombre fini et
prédéterminé d’opérations.

Math-Info-ScHS (ENSAM) CMN 16 / 105

 Élimination de Gauss Comment exploiter les algorithmes des matrices triangulaires?

Remarque sur les méthodes directes

Une méthode directes est essentiellement une méthode basée
sur la remontée (descente). Pour une telle méthode on peut faire
un estimé a priori du nombre d’opérations en virgule flottante et
du temps de calcul (une remontée ou une descente exige ≈ n2
opérations élémentaires). La plupart du temps la limitation dans
l’utilisation de méthodes directes vient de l’espace mémoire requis
pour la résolution.
À l’opposée de ces méthodes on trouve les méthodes itératives.
Ces méthodes ne garantissent pas la convergence et ne peuvent
garantir un nombre fixe d’opérations pour l’obtention d’une
solution. Cependant la plupart des méthodes itératives exigent
moins d’espace mémoire.

Math-Info-ScHS (ENSAM) CMN 17 / 105

 Élimination de Gauss Transformation des matrices

Transformation = produit par une matrice

On veut garder la linéarité du système alors on transforme avec
les opérations de bases (+, −, ×, ÷) avec des scalaires sur et
entre les lignes du système.
Une transformation de ce type est représentable par une matrice
multipliant la matrice A et le vecteur b du système:

WAx = Wb

Quels conditions sur W ? Si on a la solution (unique) de

WAx = Wb alors on veut que cette solution soit la solution de
Ax = b. Alors il faut que la matrice W soit inversible!
On s’assurera de la validité de notre transformation en s’assurant
que son déterminant est non nul (ou quel est inversible).

Math-Info-ScHS (ENSAM) CMN 18 / 105

 Élimination de Gauss Transformation des matrices

Exemple
Soit     
10 6 2 x1 44
Ax =  4 4 1   x2  =  21 
3 2 1 x3 14
Prenons
     
1 −1 −1 3 0 0 9
W = 0 1 −1  ⇒ WA =  1 2 0  Wb =  7 
0 0 1 3 2 1 14

la solution x = (3 2 1)T correspond à la solution du système original.

     
1 −1 −1 3 0 0 9
W̃ =  0 1 −1  ⇒ W̃ A =  1 2 0  W̃ b =  7 
0 0 0 0 0 0 0

Le système n’est plus valide: pas unicité donc pas correspondance

entre les deux problèmes!
Math-Info-ScHS (ENSAM) CMN 19 / 105
 Élimination de Gauss Opérations élémentaires sur les lignes

Pour réduire un système linéaire à la forme triangulaire on utilise 3

opérations élémentaires sur les lignes
Multiplication d’une ligne Li par une scalaire λ non nul:

λLi −→ Li

Addition de deux lignes:

Li + λLj −→ Li

permuter deux lignes:

Li ←→ Lj
On sait qu’en multipliant (à gauche) le système original par une
matrice inversible la solution sera inchangée. On va montrer que les 3
opérations proposées correspondent a des produits par une matrice
inversible permettant ainsi de conserver la solution originale.

Math-Info-ScHS (ENSAM) CMN 20 / 105

 Élimination de Gauss Opérations sur les lignes: matrice M

A l’opération multiplication de la ligne i par λ non nul: λLi −→ Li

correspond le produit de Ax = b avec la matrice élémentaire M:

mjj = 1 j = 1, ..., i − 1, 1 + 1, ..., n mii = λ mjk = 0 sinon

On sait que detM = λ donc la matrice est inversible si λ 6= 0

   
1 0 ... ... 0 1 0 ... ... 0
 ..   .. 
 0 . 0 ... 0   0 . 0 ... 0 
−1
M= 0 M = 0 ←i
   
1
 ... λ ... 0 
  ... λ ... 0 
 0 ... 1 ...   0 ... 1 ... 
0 ... ... 1 0 ... ... 1

Math-Info-ScHS (ENSAM) CMN 21 / 105

 Élimination de Gauss Opérations sur les lignes: matrice P

A l’opération intervertir deux lignes du système: Li ←→ Lj correspond

le produit de Ax = b avec la matrice élémentaire P (matrice de
permutation) obtenue en permutant les lignes i et j de la matrice
identité  
1 0 ... ... 0
 0 ... 0 . . . 0 
 
P = j →  0 ... 0
 
1 0 
i→  0
 
1 0 ... 
0 ... ... 1
Puisque PP = I la matrice P est inversible et son inverse P −1 = P. On
a aussi que detP = −1

Math-Info-ScHS (ENSAM) CMN 22 / 105

 Élimination de Gauss Opérations sur les lignes: matrice T

A l’opération addition de deux lignes Li + λLj −→ Li correspond le

produit de Ax = b avec la matrice élémentaire T

tkk = 1 k = 1, ..., n tij = λ tkl = 0 sinon

On sait que detT = 1 alors T est inversible (il suffit d’enlever ce que
l’on a ajouter)
   
1 0 ... ... 0 1 0 ... ... 0
 ..   .. 
 0 . 0 ... 0   0 . 0 ... 0 
−1
T = j→  0 T = 0
   
... 1 ... 0  ... 1 ... 0 
i→  0
   
λ 1 ...   0 −λ 1 ... 
0 ... ... 1 0 ... ... 1

Math-Info-ScHS (ENSAM) CMN 23 / 105

 Élimination de Gauss Elimination de Gauss

La méthode d’élimination de Gauss consiste à utiliser les opérations

élémentaires afin de réduire le système sous la forme triangulaire
supérieure.
On distingue deux cas:
Sans pivotement: on ne permet pas la permutation de deux lignes
(la transformation P est interdite); dans certain cas on ne pourra
pas obtenir de matrice triangulaire sup.!
Avec pivotement: on permet la permutation des lignes; dans ce
cas on obtiendra toujours une matrice triangulaire sup. si la
matrice est inversible.
Dans la pratique on ne construit pas les matrices (de types M,T et P)
nécessaires à la transformation.
Cependant la méthode de Gauss est valide parce qu’elle consiste à
multiplier le système de départ par une suite de matrices inversibles.
On verra plus tard comment ces matrices nous seront utiles.
Math-Info-ScHS (ENSAM) CMN 24 / 105
 Élimination de Gauss Elimination de Gauss

Dans l’application de la méthode de Gauss on procède par élimination

des termes sous la diagonale en utilisant la notion de matrice
augmentée.

Définition
La matrice augmentée du système Ax = b est la matrice de dimension
n × n + 1 obtenu en ajoutant le membre de droite b à la matrice A:
 
a11 ... a1i ... a1n b1
 .. . . .. . . .. 
 .
 . . . . 

 ai1 ... aii ... ain bi 
 
 .. . . .. . . .. 
 . . . . . 
an1 ... ani ... ann bn

Cette notation est très utile puisque les opérations élémentaires

doivent être faites sur la matrice A et sur b simultanément.

Math-Info-ScHS (ENSAM) CMN 25 / 105

 Élimination de Gauss Élimination de Gauss sans permutation

Exemple 3.8
On applique la méthode de Gauss sur la matrice augmentée
 
2 1 2 10
A =  6 4 0 26 
8 5 1 35

Math-Info-ScHS (ENSAM) CMN 26 / 105

 Élimination de Gauss Élimination de Gauss avec permutation

Avec une permutation: Exemple 3.9

On applique la méthode de Gauss sur la matrice augmentée
 
1 1 2 1 2
 2 2 5 3 4 
A=  1 3 3 3 −2 


1 1 4 5 −2

Math-Info-ScHS (ENSAM) CMN 27 / 105

 Décomposition LU

Et si je veux changer b?
Trois options:
1 Je recommence le processus avec la nouvelle valeur de b car je
dois transformer A et b simultanément dans la méthode de Gauss
2 j’attend d’avoir TOUT les b que je veux et je fais Gauss sur la
matrice augmentée de tout les b
3 J’essaie de garder l’information ayant servi à la transformation de
la matrice lors de la première résolution pour résoudre avec les
autres b
Revenons à nos exemples pour voir ce qu’on peut faire en tenant
compte des matrices T ,M et P cette fois

Math-Info-ScHS (ENSAM) CMN 28 / 105

 Décomposition LU

Exemple 3.8 revisité

On applique la méthode de Gauss sur la matrice augmentée
 
2 1 2 10
A =  6 4 0 26 
8 5 1 35

Mais on construit les 3 matrices T et on a

A = (T1−1 T2−1 T3−1 )U = LU

Math-Info-ScHS (ENSAM) CMN 29 / 105

 Décomposition LU

Avec une permutation: Exemple 3.9 revisité

On applique la méthode de Gauss sur la matrice augmentée
 
1 1 2 1 2
 2 2 5 3 4 
A=  1 3 3 3 −2 


1 1 4 5 −2

Dans ce cas on a des T et un P et on a

PA = (P4 T1−1 T2−1 T3−1 P4 T5−1 )U = LU ⇔ A = PLU

Math-Info-ScHS (ENSAM) CMN 30 / 105

 Décomposition LU

La décomposition LU permet de résoudre avec différents b

Cette décomposition donne toute l’information nécessaire pour
résoudre à “répétition”. Puisque L et U sont construit on a

Ax = b ⇔ (LU)x = b ⇔ L(Ux) = b

Si on pose y = Ux alors la résolution de Ax = b peut se faire en deux

temps:
1 résoudre Ly = b par une descente triangulaire
2 résoudre Ux = y par une remontée triangulaire
Une fois le problème résolu pour un b, en ayant garder L je peux
changer le b et résoudre à nouveau sans refaire de transformation du
système.

Math-Info-ScHS (ENSAM) CMN 31 / 105

 Décomposition LU

Remarques
la méthode de Gauss est une décomposition LU où on “jete“ L.
Plus la peine d’en parler (sauf pour les exercices :-)
les permutations causent des problèmes (empêchent d’avoir un
LU ”propre”) on en reparlera.
Pour le moment on va faire comme si A était permuter avant qu’on
arrive...

Unicité de LU
La décomposition LU se fait grâce aux transformations de Ax = b mais
on pourrait ajouter des transformations (superflues). Par exemple on
pourrait finir en s’assurant que les pivots (termes diagonaux non nuls)
soit tous égaux à 1 en appliquant des transformations M. Donc la
décomposition LU n’est pas unique.

Math-Info-ScHS (ENSAM) CMN 32 / 105

 Décomposition LU Décomposition de Crout

Pour éviter les ambiguités possible (l’opérateur == est mal défini) on

voudrait assurer l’unicité de la décomposition. Pour se faire on va
imposer une conditions supplémentaire à notre décomposition.
Deux méthodes “classiques“
Décomposition de type Matlab: Lii = 1 ∀i (pas de division dans
la descente)
Décomposition de Crout: Uii = 1 ∀i (pas de division dans la
remontée)

Décomposition de Crout
Algorithme composé de trois parties
1 factorisation LU avec Uii = 1 ∀i
2 Résolution de Ly = b par descente triangulaire
3 Résolution de Ux = y par remontée triangulaire

Math-Info-ScHS (ENSAM) CMN 33 / 105

 Décomposition LU Décomposition de Crout

Comment on fait la factorisation

On identifie termes à termes les composantes de L et U. À noter que
c’est un système linéaire de n2 équations et inconnues qui a une
solution unique quand A est inversibles.
Une particularité de l’approche est que l’on résout alternativement
pour une colonne de L puis pour une ligne de U en commençant par la
colonne 1 donc on construit:

Li1 , U1i , ..., Lin−1 , Un−1i , Lin

Cette stratégie nous permettra de réduire l’espace mémoire lors de la

factorisation.

Math-Info-ScHS (ENSAM) CMN 34 / 105

 Décomposition LU Décomposition de Crout

Factorisation
A1i
Li1 = Ai1 i = 1, ..., n, U1i = L11 i = 2, ..., n
i = 2, ..., n − 1
i−1
P
Lii = Aii − Lik Uki
k =1
j = i + 1, ..., n
i−1
P
Lji = Aji − Ljk Uki
k =1
i−1
P
Aij − Lik Ukj
k =1
Uij = Lii
n−1
P
Lnn = Ann − Lnk Ukn
k =1
descente sur Ly = b
i−1
P
bi − Lik yk
b1 k =1
y1 = L11 yi = Lii i = 2, ..., n
remonté sur Ux = y
n
P
xn = yn xi = yi − Uik xk i = n − 1, ..., 1
k =i+1

Math-Info-ScHS (ENSAM) CMN 35 / 105

 Décomposition LU Décomposition de Crout

Remarque
L’algo ne fonctionne que si Lii 6= 0 ∀i.
Si on a pas fait de permutations alors detA = detLdetU mais
detU = 1 et

detA 6= 0 ⇔ detL 6= 0 ⇔ Lii 6= 0 ∀i

Si on a fait des permutations detA = detPdetLdetU Puisque

detP 6= 0 alors

detA 6= 0 ⇔ detL 6= 0 ⇔ Lii 6= 0 ∀i

Donc si A est inversible l’algorithme va bien.

Math-Info-ScHS (ENSAM) CMN 36 / 105

 Décomposition LU Décomposition de Crout

Remarque
D’un point de vue pratique une fois la factorisation LU faite on a plus
de raison de conserver la matrice A. En fait notre procédure de
factorisation est faite pour que l’on puisse ”écraser“ A au fur et à
mesure de la factorisation. Ce qui permet de minimiser l’espace de
stockage puisqu’on a pas simultanément A, L et U

Math-Info-ScHS (ENSAM) CMN 37 / 105

 Décomposition LU Décomposition de Crout

Definition 3.4
La notation compacte de LU consiste a remplacer la matrice A par la
matrice

Les coefficients de la diagonale de U ne sont pas gardé puisqu’ils sont

tous égaux à 1 mais on ne les néglige pas pour autant.
Cette notation revient à stocker:

L + (U − I) = L + U − I

à la place de A

Math-Info-ScHS (ENSAM) CMN 38 / 105

 Décomposition LU Décomposition de Crout

Exemple 3.11
Utilisation de la décomposition LU avec notation compacte pour
   
3 −1 2 12
A= 1 2 3  b =  11 
2 −2 −1 2

Math-Info-ScHS (ENSAM) CMN 39 / 105

 Décomposition LU Décomposition de Crout avec permutations

Permutations
On gardera un vecteur de permutation (position ou numérotation) qui
indiquera les permutations effectuées sur la matrice. A la fin de la
factorisation on obtient L, U et le vecteur de permutation.
Si nécessaire on reconstruit la matrice P en appliquant les
permutations sur une matrice identié. Mais c’est normalement inutile.
Comment résoudre une fois la factorisation terminée?

Ax = b et PA = LU ⇔ PAx = LUx = Pb ⇔ Ly = Pb et Ux = y

On ne construit pas P mais on permutent b avec le vecteur de

permutation puis on résoud par descente et remontée.
Voir exemple 3.12
ATTENTION P est en général le produit de plusieurs permutations

P = Pp Pp−1 ..P1 ⇔ P −1 = P1 P2 ...Pp = P T

Math-Info-ScHS (ENSAM) CMN 40 / 105

 Décomposition LU Factorisation de Cholesky

Dans le cas d’une matrice symétrique il est possible déconomiser de

l’espace en imposant une condition supplémentaire à la décomposition
LU. On impose que
U = LT ⇒ A = LLT
A étant symétrique on ne stocke que sa partie inférieure et sa
diagonale. Alors on peut faire la même chose sur la décomposition et
ne stocke que L.
La construction est basée sur Crout.

Math-Info-ScHS (ENSAM) CMN 41 / 105

 Décomposition LU Factorisation de Cholesky

Algorithme de Cholesky
1
L11 = (A11 ) 2
Ai1
Li1 = L11 i = 2, ..., n
k = 2, ..., n
−1
kP 1
Lkk = (Akk − L2kj ) 2
j=1
i = k + 1, ..., n
−1
kP
Aik − Lij Lkj
j=1
Lik = Lkk

Remarques
Tout comme pour la décomposition LU on peut introduire une
notation compacte et ”écraser“ A lors de la création de L
La factorisation n’est pas unique: A = (−L)(−LT ) = L̃L̃T on rend
la factorisation unique en imposant des pivots positifs: Lii > 0.

Math-Info-ScHS (ENSAM) CMN 42 / 105

 Décomposition LU Factorisation de Cholesky

Exemple 3.13
Utilisation de la factorisation de Cholesky avec notation compacte pour
 
4 6 2
A =  6 10 5 
2 5 14

Math-Info-ScHS (ENSAM) CMN 43 / 105

 Décomposition LU Existence de la factorisation de Cholesky

Existence de la factorisation de Cholesky

Soit  
−1 2
A=
2 6
√
On ne peut pas faire Cholesky car L11 = A11 n’est pas réels.
En fait on aura pas existence de la factorisation si une des racines (les
pivots) n’est pas réelles

Définition 3.5
Une matrice A est définie positive ssi

Ax · x > 0 ∀x 6= 0

Math-Info-ScHS (ENSAM) CMN 44 / 105

 Décomposition LU Existence de la factorisation de Cholesky

Condition pour l’existence de la factorisation de Cholesky

On peut montrer que

A symétrique définie positive ⇔ A a une factorisation de Cholesky

La vérification de la positivité est difficile, la meilleur approche:

appliquer la factorisation...
Dans le cas ou on n’a pas la positivité de la matrice on applique
simplement la factorisation LU

Math-Info-ScHS (ENSAM) CMN 45 / 105

 Décomposition LU Exemples pour l’existence de Cholesky

Exemple d’existence de la factorisation de Cholesky

Soit  
4 6 2
A =  6 10 5 
2 5 14
On calcule les déterminants des sous-matrices principales.

Exemple d’existence de la factorisation de Cholesky

Soit  
−1 2
A=
2 6
Alors detA11 = −1 < 0 et la factorisation de Cholesky n’existe pas.

Math-Info-ScHS (ENSAM) CMN 46 / 105

 Décomposition LU Retour sur l’inverse de A

Comment inverser A?
La méthode usuelle (i.e. à la main) consiste a augmenter la matrice
avec la matrice identité et à faire une élimination jusqu’à ce qu’on ait
une ”identité à gauche”

En utilisant la décomposition LU cela revient à résoudre n fois le

système Ax = b avec b = ei (les vecteurs de la base euclidienne) pour
i=1,...,n.
Clairement il ne faut pas utiliser cette approche pour résoudre un
problème: on ferait n résolution et un produit matrice vecteur.

Math-Info-ScHS (ENSAM) CMN 47 / 105

 Décomposition LU Combien ça coûte?

Nombre d’opérations ×/÷

Si on veut absolumment calculer l’inverse. Il faudra faire une
factorisation LU suivi de n descente et de n remontée. Utilisant les
ordres de grandeurs sur le nombre d’opérations ×/÷ on a

n Ax = b A−1 b

10 333 1 333
100 333 333 1 333 333
1000 333 333 333 1 333 333 333

Pour un système 1000 × 1000 on fait 109 opérations ×/÷ de plus pour
résoudre le système.
1000 est un petit chiffre dans le monde réel...

Math-Info-ScHS (ENSAM) CMN 48 / 105

 Décomposition LU Exemple

Exemple 3.15
On se donne une matrice A obtenu par factorisation LU avec
permutation:
     
0 2 1 3 0 0 1 0 1/3
A= 1 0 0  L= 0 2 0  U =  0 1 1/2 
3 0 1 1 0 −1/3 0 0 1

Avec le vecteur de permutations

 
3
O= 1 
2

On va calculer l’inverse. Il s’agit SURTOUT d’un exemple de résolution

avec un vecteur de permutation

Math-Info-ScHS (ENSAM) CMN 49 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Exemple 3.16
Soit  
1.012 −2.132 3.104
A =  −2.132 4.096 −7.013 
3.104 −7.013 0.014
On veut appliquer la factorisation LU mais avec une mantisse fixe de 4
chiffres.
Cette restriction engendre des “perturbations” par rapport à la
factorisation “humaine” (mantisse infinie). Quel est l’impact de cette
perturbation lors de la résolution?

Math-Info-ScHS (ENSAM) CMN 50 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Exemple 3.16 suite

Avec virgule flottante a 4 chiffres
 
1.012 −2.107 3.067
L + U − I =  −2.132 −0.3960 1.197 
3.104 −0.4730 −8.940

Avec virgule flottante humaine:

 
1.012 −2.106719... 3.067193...
L + U − I =  −2.132 −0.395525... 1.197757... 
3.104 −0.473743... −8.939140

Quel sont les effets sur la solution de Ax = b pour b donné ?

Math-Info-ScHS (ENSAM) CMN 51 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Effet des opérations arithmétiques en virgules flottantes

On retourne aux considérations du chapitre 1. Pour le résoudre un
système Ax = b il est à prévoir que des erreurs de troncatures vont
apparaitre.
On verra que pour certaines matrices les effets sur la solution seront
négligeables. Pour certaines autres matrices les variations dans les
entrées de la matrice induiront d’importantes erreurs dans la solution
obtenues. De telles matrices sont dites mal conditionnées

Math-Info-ScHS (ENSAM) CMN 52 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Exemple 3.17
Soit une matrice A et une perturbation à de A
     
1. 2 10 1. 2
A= b= Ã =
1.1 2 10.4 1.05 2

La solution exacte avec A est (4, 3) et pour à la solution est (8, 1).
Une “petite” variation de A produit une “grande” variation de la solution.

Puisque l’arithmetique flottante (mantisse finie) mènent à de petites

perturbations dans la décomposition LU. Il est possible que
d’importants effets apparaissent sur la solution du système.
On peut produire d’importantes erreurs sur la solution numériques via
la décomposition LU et le procèdé de descente et remontée.

Math-Info-ScHS (ENSAM) CMN 53 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Recettes simples pour éviter les dégats

Dépendant de la matrice A on pourra avoir ou non de large erreurs.
On verra par la suite comment détecter de telle matrice. Pour le
moment voici deux corrections simples à appliquer lors de la résolution
pour éviter deux des principales causes d’erreurs:
division par des nombres trop petits (pivots trop petits)
la soustraction/addition avec des nombres ayant des ordres de
grandeurs trop éloignées

Math-Info-ScHS (ENSAM) CMN 54 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Exemple 3.18
Soit    
0.0003 3.0000 2.0001
A= b=
1.0000 1.0000 1.0000
La solution est x = [1/3, 2/3]T , en faisant la décomposition LU avec 4
chiffres on a
0.3000 × 10−3 0.1000 × 105 0.6667 × 104
   
L+U−I = x=
0.1000 × 101 −0.9999 × 104 0.6667 × 100

La solution par LU est “loin” de la solution réelle. Comment éviter cette

erreur?
On fait une division “presque par zéro“ à cause de A11 c’est lui qui
modifie U12 et L22 . Alors on va éviter cette division en faisant une
permutation.

Math-Info-ScHS (ENSAM) CMN 55 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Encore des permutations...

Soit    
1.0000 1.0000 1.0000
A= b=
0.0003 3.0000 2.0001
en faisant la décomposition LU avec 4 chiffres on a

0.1000 × 101 0.1000 × 101 0.3333 × 100

   
L+U −I = x=
0.3000 × 10−3 0.3000 × 101 0.6667 × 100

La solution par LU est “proche” de la solution réelle.

Recette 1
On peut éviter une partie des erreurs en faisant une permutation
systèmatique des lignes lors de la décomposition: on s’assure d’avoir
TOUJOURS le pivot de plus grand en valeur absolue sur la diagonale.

Math-Info-ScHS (ENSAM) CMN 56 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Exemple 3.19

1 6 9
A= 2 1 2 
3 6 9

Math-Info-ScHS (ENSAM) CMN 57 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Exemple 3.20
   
2.0000 100000 100000
A= b=
1.0000 1.0000 2.0000
La solution est x = [1.00002, 0.99998]T , en faisant la décomposition
LU avec 4 chiffres on a
0.2000 × 101 0.5000 × 105 0.0000 × 100
   
L+U −I = x=
0.1000 × 101 −0.5000 × 105 0.1000 × 101

La solution par LU est “loin” de la solution réelle. Comment éviter cette

erreur?
En calculant L22 on fait une soustraction de chiffres d’ordre de
grandeur très différents (opérations à éviter Chap 1) à cause de A12 .
Alors on va éviter cette opérations en multipliant la ligne par une
constante: utilisation d’une transformation M.

Math-Info-ScHS (ENSAM) CMN 58 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Matrice M
On divise la première ligne de A par 100000 ce qui donne

0.2000 × 10−4 0.1000 × 101 0.1000 × 101

   
A= b=
0.1000 × 101 0.1000 × 101 0.2000 × 101

en faisant la décomposition LU avec 4 chiffres AVEC permutation on a

0.1000 × 101 0.1000 × 101 0.1000 × 101

   
L+U −I = x =
0.2000 × 10−4 0.1000 × 101 0.1000 × 101

La solution par LU est “proche” de la solution réelle.

Math-Info-ScHS (ENSAM) CMN 59 / 105

 Effets de l’arithmétique flottante Propagation d’erreurs

Definition 3.7
Une mise à l’échelle d’une matrice A consiste à diviser chaque ligne
de A par le plus grand termes en valeur absolue de la ligne
correspondante.

Recette 2
On peut éviter une partie des erreurs en faisant une mise à l’échelle
systèmatique d’une matrice avant la décomposition.

Attention au déterminant
L’introduction de matrice M à un effet important sur le déterminant. Si
on multiplie une matrice par un scalaire on récupère son déterminant
en divisant le déterminant de L par ce scalaire.

Math-Info-ScHS (ENSAM) CMN 60 / 105

 Conditionnement d’une matrice

On veut étudier le comportement de la solution calculée par rapport à

la solution exacte pour un système linéaire. On va d’abord introduite
une mesure de la distance dans l’espace Rn .

Definition 3.8
Une norme vectorielle est une application de Rn dans R notée kk
vérifiant
1 Positivité stricte:
kxk> 0 ∀x ∈ Rn , x 6= 0
2 ∀α ∈ R
kαxk= |α|kxk ∀x ∈ Rn
3 Inégalité du triangle

kx + y k≤ kxk+ky k ∀x, y ∈ Rn

Math-Info-ScHS (ENSAM) CMN 61 / 105

 Conditionnement d’une matrice

Definition 3.9
La norme euclidienne d’un vecteur x notée kxk2 est définie par
 1
2
kxk2 = x12 + x22 + ... + xn2

évidemment la norme euclidienne est une norme vectoriel.

Démonstration p 135-136

Definition 3.10
Il existe d’autre norme sur Rn .
n
P
La norme l1 : kxk1 = |xi |
i=1
n
La norme l∞ : kxk∞ = max |xi |
i=1

Math-Info-ScHS (ENSAM) CMN 62 / 105

 Conditionnement d’une matrice

Definition 3.11
On peut aussi introduire une norme matricielle. On l’a définie comme
une application A :→ kAk vérifiant
1 Positivité stricte:
kAk> 0 ∀A, A 6= 0
2 ∀α ∈ R
kαAk= |α|kAk ∀A
3 Inégalité du triangle

kA + Bk≤ kAk+kBk ∀A, Bde dimensions compatibles

kABk≤ kAkkBk ∀A, Bde dimensions compatibles

Une norme matricielle DOIT vérifier 1-4.

Math-Info-ScHS (ENSAM) CMN 63 / 105

 Conditionnement d’une matrice

Norme induite
On peut construire une norme matricielle à partir d’une norme
vectorielle:
kAk= sup kAxk
kxk=1

Une telle norme est dite norme matricielle induite par la norme
vectorielle. On dit aussi que cette norme est subordonnée à la
norme vectorielle

Toute les normes ne sont pas induite

La norme de Frobenius définie par
 1
n X
n 2
X
kAkF =  A2ij 
i=1 j=1

est une norme matricielle, s’apparentant à la norme euclidienne mais

elle n’est induite par aucune norme.
Math-Info-ScHS (ENSAM) CMN 64 / 105
 Conditionnement d’une matrice

Théorème 3.5
Les normes induites par les normes l1 et l∞ sont
n
P
kAk1 = sup kAxk1 = max |Aij |
kxk1 =1 1≤j≤n i=1
n
P
kAk∞ = sup kAxk∞ = max |Aij |
kxk∞ =1 1≤i≤n j=1

La norme kAk1 consiste à trouver la colonne dont la somme des

éléments (en valeur absolue) est la plus grande.
La norme kAk∞ consiste à trouver la ligne dont la somme des
éléments (en valeur absolue) est la plus grande.

Math-Info-ScHS (ENSAM) CMN 65 / 105

 Conditionnement d’une matrice

Norme induite par la norme euclidienne

La norme induite par la norme euclidienne est définie par
 1
2
kAk2 = sup kAxk2 = ρ(AAT )
kxk2 =1

où ρ est la plus grande valeur propre (en valeur absolue) de la matrice

AAT .
Pourquoi s’intéresser au norme matricielle? On va manipuler des
vecteurs résultants de produit matrice vecteur. Ce qui nous amène à
une définition fort importante

Math-Info-ScHS (ENSAM) CMN 66 / 105

 Conditionnement d’une matrice

Définition 3.12
On dit que la norme matricielle et la norme vectorielle sont
compatibles si
kAxk≤ kAkkxk
pour toutes matrice A et vecteurs x de dimensions cohérentes.
Attention au mélange de normes: à gauche norme vectorielle, à droite
matricielle et vectorielle (respectivement).

Math-Info-ScHS (ENSAM) CMN 67 / 105

 Conditionnement d’une matrice

Normes compatibles
On peut montrer que
kAxk1 ≤ kAk1 kxk1 la norme 1 des matrices et vecteurs est
compatible
kAxk∞ ≤ kAk∞ kxk∞ la norme ∞ des matrices et vecteurs est
compatible
kAxk2 ≤ kAkF kxk2 la norme Frobenius des matrices est
compatible avec la norme euclidienne.
La dernière relation montre qu’il n’y a pas de lien entre norme induite
et norme compatible.

Math-Info-ScHS (ENSAM) CMN 68 / 105

 Conditionnement d’une matrice

Soit le système Ax = b et x ∗ une solution obtenue par arithmétique

flottante. On veut étudier l’erreur par rapport à la solution exacte:
e = x − x∗
Si tout va bien kek= kx − x ∗ k est petit.
Posons r = b − Ax ∗ r est le résidu de notre résolution.
r = b − Ax ∗ = Ax − Ax ∗ = A(x − x ∗ ) = Ae
On cherche à avoir une borne inférieure et supérieure sur kek. On va
se servir de la compatibilité des normes (on choisi une norme
compatible)
A−1 r = e ⇒ kek= kA−1 r k≤ kA−1 kkr k
On a aussi
kr k= kAek≤ kAkkek
Ce qui nous donne
kr k
≤ kek≤ kA−1 kkr k
kAk
Math-Info-ScHS (ENSAM) CMN 69 / 105
 Conditionnement d’une matrice

kr k
≤ kek≤ kA−1 kkr k
kAk
On veut faire apparaitre x dans cette relation (on voudrait l’“erreur
relative“) alors on revient au système:

1 kAk
kbk = kAxk≤ kAkkxk⇒ ≤
kxk kbk

1 1
kxk= kA−1 bk≤ kA−1 kkbk⇒ ≤
kA−1 kkbk kxk
En combinant ces deux inégalités on a

1 1 kAk
≤ ≤
kA−1 kkbk kxk kbk

Math-Info-ScHS (ENSAM) CMN 70 / 105

 Conditionnement d’une matrice

On multiplie la dernière relation par kek

kek kek kAk

≤ ≤ kek
kA−1 kkbk kxk kbk

et on a
kr k kek kA−1 kkAkkr k
≤ ≤
kAkkA−1 kkbk kxk kbk

Definition 3.13
Le conditionnement d’une matrice A noté condA est défini par

condA = kA−1 kkAk

Math-Info-ScHS (ENSAM) CMN 71 / 105

 Conditionnement d’une matrice

Remarque
condA dépend de la norme matricielle utilisée
Puisque 1 ≤ kIk parce que

kAk= kAIk≤ kAkkIk

on a
1 ≤ kIk≤ kAA−1 ≤ kA−1 kkAk
et
1 ≤ condA < ∞

Théorème 3.6
1 kr k kek kr k
≤ ≤ CondA
condA kbk kxk kbk

Math-Info-ScHS (ENSAM) CMN 72 / 105

 Conditionnement d’une matrice

Concernant le conditionnement
1 kr k kek kr k
≤ ≤ CondA
condA kbk kxk kbk

1 le terme du milieu représente l’erreur relative

2 Si condA est près de 1 alors

kek kr k
≈
kxk kbk

est la solution est bonne si kr k est petit.

3 si condA est grand alors l’erreur relative peut être grande.
4 même si kr k est petit il est possible que l’erreur soit grande
(manque le conditionnement petit)
5 Plus condA augmente plus l’algorithme de résolution doit être
précis (pour réduire kr k)
Math-Info-ScHS (ENSAM) CMN 73 / 105
 Conditionnement d’une matrice

Concernant le conditionnement
1 kr k kek kr k
≤ ≤ CondA
condA kbk kxk kbk

1 Même si condA est petit si la résolution est imprécise on pourra

avoir des résultats erronés (manque kr k est petit).
2 On peut montrer que
 
kekkbk kr kkxk
max , ≤ condA
kxkkr k kekkbk

3 condA dépend de la norme matricielle choisie (géneralement la

norme ∞) mais ne dépend pas de l’algorithme de résolution du
système. C’est une mesure de la sensibilité de la matrice au
perturbation lors de manipulations matricielles.

Math-Info-ScHS (ENSAM) CMN 74 / 105

 Conditionnement d’une matrice

Erreurs dans A
On veut résoudre le système Ax = b mais à cause d’erreur de
troncature ou dans la mesure des éléments de A on résout plutôt le
système
(A + E)x = b
Soit x ∗ solution du système perturbé et x la solution du système
original.

x = A−1 b = A−1 ((A + E)x ∗ ) = (I + A−1 E)x ∗ = x ∗ + A−1 Ex ∗

Donc si les normes sont compatibles

kek= kx − x ∗ k= kA−1 Ex ∗ k≤ kA−1 kkEkkx ∗ k

Et on a
kx − x ∗ k kEk
∗
≤ condA
kx k kAk

Math-Info-ScHS (ENSAM) CMN 75 / 105

 Conditionnement d’une matrice

Erreurs de troncature sur A

kx − x ∗ k kEk
∗
≤ condA
kx k kAk

kx−x ∗ k
kx ∗ k approxime l’erreur relative
kEk
kAk ressemble à l’erreur relative sur les entrées de la matrice A
Si condA est petit et kEk est petit alors l’erreur relative sera petite:
bonne approximation de x
Si condA est grand on ne peut rien dire

Math-Info-ScHS (ENSAM) CMN 76 / 105

 Conditionnement d’une matrice

Erreurs dans A
kx − x ∗ k kEk
∗
≤ condA
kx k kAk

Si E représente l’erreur de troncature sur une machine alors

|Eij | ≤ m |Aij |

par définition de la troncature donc

X X
|Eij | ≤ m |Aij | ⇔ kEk∞ ≤ m kAk∞
j j

d’où
kx − x ∗ k
≤ m condA
kx ∗ k
Si le conditionnement est mauvais il faut plus de précision sur les
nombres
Math-Info-ScHS (ENSAM) CMN 77 / 105
 Conditionnement d’une matrice

Quoi faire de condA?

On c’est donné cette mesure qu’en fait-on?
Exemple 3.16 (3.24) Le conditionnement de A = 270.51
correspond à une matrice sensible au perturbations. Ne veut pas
dire que je n’aurai pas la bonne solution mais que j’ai de bonne
chance d’obtenir une solution erronnées sans un bon algorithme
Exemple 3.18 (3.26) condA = 4 ce qui est faible indiquant que si
tout va bien j’aurai une bonne solution. Rappelons que sans
permutations la solution est mauvaise ce qui ne contredit pas le
conditionnement mais indique que mon algorithme est mauvais.
En incluant les permutations on obtient une bonne solution.
   
1. 2 −1 −10. 10.
A= A =
1.1 2 −5.5 5.
cond A = 62 indiquant un mauvais conditionnement.

Math-Info-ScHS (ENSAM) CMN 78 / 105

 Conditionnement d’une matrice

Quoi faire de condA?

Est-ce que je peux diminuer le conditionnement en changeant de
norme? Non, on a ”équivalence des normes“ faisant en sorte que
rien ne changera en prenant d’autre norme.
Alors quoi?
Mauvais conditionnement == appliquer les recettes 1 et 2
OBLIGATOIREMENT et être critique des résultats obtenus et
choisir avec soin l’algorithme de résolution (contrôler la qualité
autrement, de manière ad hoc)
Bon conditionnement == les recettes ne sont pas nécessaires, on
peut croire que tout va bien mais on doit encore être critique des
résultats obtenus (contrôler la qualité autrement, de manière ad
hoc).

Math-Info-ScHS (ENSAM) CMN 79 / 105

 Méthodes itératives pour la solution d’un système linéaire Principe des méthodes itératives

Supposons que nous cherchions à résoudre le système suivant:


 2x1 +4x2 −2x3 −6x4 = −6
1x1 +3x2 −6x4 =0
 3x1 −x2 +x3 −6x4 =8
−x2 +2x3 −6x4 =6

et que nous connaissions un vecteur x ∗ 1 , x ∗ 2 , x ∗ 3 , x ∗ 4 proche de la

solution.  −6 − 4x ∗ 2 + 2x ∗ 3 + 6x ∗ 4
x1 =
2


 x2 = 0 − 1x ∗ 1 − 0x ∗ 3 + 6x ∗ 4


3
8 − 4x ∗ 1 + x ∗ 2 + 6x ∗ 4
 x3 =
 1
6 − 4x ∗ 1 + x ∗ 2 − 2x ∗ 3


x4 =

−6

La politique des « petits » pas...

on cherche à créer une suite de vecteurs {x (k ) , k ∈ IN} qui converge
vers xe solution de Axe = b

Math-Info-ScHS (ENSAM) CMN 80 / 105

Trois méthodes itératives

(k )
Jacobi xJ , k ∈ IN
(k )
Gauss Seidel xGS , k ∈ IN
(k )
Relaxation xR , k ∈ IN

a chaque itération...
le vecteur d’erreurs : e(k ) = x (k ) − xe doit diminuer

comment choisir ?
il en existe plein d’autres...
la question centrale : quand convergent’elles ?
seconde question : à quelle vitesse

parmi celles qui convergent on prend celle qui converge le plus vite !
La méthode de Jacobi
i−1
X n
X
Ax = b ⇐⇒ Aij xj + Aii xi + Aij xj = bi i = 1, n
j=1 j=i+1

Fonction x ← Jacobi(A, b, x)

Tant que (on n’a pas convergé) faire

Pour i = 1, n faire
bi − i−1
P Pn
j=1 Aij xj − j=i+1 Aij xj
yi =
Aii
Fin Pour
Pour i = 1, n faire
xi = yi
Fin Pour
Fait
Exemple : méthode de Jacobi
       
3 1 −1 1 0 1/3
A =  1 2 0  avec b = 1 et x0 = 0 x1 = 1/2
−1 1 4 1 0 1/4

1−1/2+1/4 3
y1 = 3 = 12
1−1/3
y2 = 2 = 13
1+1/3−1/2 5
y3 = 4 = 24
 
3/12
x2 =  1/3 
5/24

norm(A ∗ x0 − b) = 1.7321 
0.2500 0.2941

norm(A ∗ x1 − b) = 0.4488
 
x2 A\b = 0.3333 0.3529
norm(A ∗ x2 − b) = 0.1718 0.2083 0.2353
Les itérations de Gauss Seidel

Fonction x ← Gauss_Seidel(A, b, x)

Tant que (on n’a pas convergé) faire

Pour i = 1, n faire
bi − i−1
P Pn
j=1 Aij xj − j=i+1 Aij xj
xi = ;
Aii
Fin Pour
Fait

plus « simple », mais non parallélisable

Les itérations de relaxation

Fonction x ← Relaxation(A, b, x, ω)

Tant que (on a pas convergé) faire

Pour i = 1, n faire
Pi−1
Aij xj − nj=i+1 Aij xj
P !
bi − j=1
xi = ω + (1 − ω)xi
Aii
Fin Pour
Fait

pour ω = 1 c’est Gauss Seidel

Quand s’arrêter ?
kxnew − xold k < ε
un nombre maximal d’itérations est atteint

Fonction y ← Gauss_Seidel(A, b, x, ε, n_ite_max)

Pour i = 1, n faire
yi = xi xi = yi + 2 ∗ ε
Fin Pour
n_ite = 0
Tant que (kx − yk > ε ET n_ite < n_ite_max) faire
Pour i = 1, n faire
x = xold
xi = yi
Fin Pour y = xnew
Pour i = 1, n faire
bi − i−1
P Pn
j=1 Aij yj − j=i+1 Aij xj
yi =
Aii
Fin Pour
n_ite = n_ite + 1
Fait
Formulation matricielle des itérations

Considérons maintenant la suite de vecteurs x(k ) construit à l’aide des

matrices carrées M et N, ainsi que d’un vecteur b:

Mx(k +1) = Nx(k ) + b

On peut réécrire ces itérations à l’aide de la matrice carrée C = M −1 N

et du vecteur d = M −1 b:

x(k +1) = M −1 N x(k ) + M −1 b

= C x(k ) + d

pour un vecteur x (0) donné.

 Décomposition d’une matrice

Revenons au système Ax = b. Nous pouvons décomposer la matrice

A
 A  =  D  +  U  +  L 
a11 a12 a13 a11 0 0 0 a12 a13 0 0 0
a21 a22 a23  =  0 a22 0  + 0 0 a23  + a21 0 0
a31 a32 a33 0 0 a33 0 0 0 a31 a32 0
diagonale triangulaire sup triangulaire inf

 Ax  =  Dx  +  (L + U)x 
a11 x1 + a12 x2 + a13 x3 a11 x1 +0 +0 0 +a12 x2 +a13 x3
a21 x1 + a22 x2 + a23 x3  =  0 +a22 x2 +0  + a21 x1 +0 +a23 x3 
a31 x1 + a32 x2 + a33 x3 0 +0 +a33 x3 a31 x1 +a32 x2 +0

Ax = Dx + Lx + Ux
Décomposition d’une matrice

Revenons au système Ax = b. Nous pouvons décomposer la matrice

A

De manière générale : A = D + L + U avec :

D = diag(a11 , a22 , ..., ann )
0 0 ... ... 0
 
a21 0 ... ... 0
L = a31 a32 ... ... 0
 
 ... 0 0
an1 an2 ... an,n−1 0
0 a12 a13 ... a1n
 
0 0 ... ... a2n 
0 ... ...
 
U=

... an−2,n−1 an−2,n 

... 0 an−1,n 
0 0 ... 0 0
Les itérations de Jacobi s’écrivent alors

Le point fixe
Ax =b
(D + L + U)x =b
Dx + (L + U)x =b
Dx = −(L + U)x + b

Les itérations de Jacobi

Dx(k +1) = −(U + L)x(k ) + b

x(k +1) = −D −1 (U + L) x(k ) + D −1

| {z b}
| {z }
Cj dj
Les itérations de Jacobi s’écrivent alors

Fonction x ← Jacobi(A, b, x)

Tant que (on a pas convergé) faire

Pour i = 1, n faire
i−1 n Fonction x ← Jacobi(A, b, x)
X X
bi − Aij xj − Aij xj
j=1 j=i+1
Tant que (on a pas convergé)
yi = faire
Aii x = D\ b − (L + U)x


Fin Pour
Fait
Pour i = 1, n faire
xi = yi
Fin Pour
Fait

x(k +1) = Cj x(k ) + dj

|{z} |{z}
−D −1 (U+L) D −1 b=D\b
Gauss Seidel et la relaxation

les itérations de Gauss Seidel s’écrivent

(D + L)x(k +1) = −Ux(k ) + b

x(k +1) = (D + L)−1 (−U) x(k ) + (D + L)−1 b
| {z } | {z }
Cg dg

et les itérations de relaxation s’écrivent

(D + ωL)x(k +1) = (1 − ω)D − ωU x(k ) + ωb

x(k +1) = ((D + ωL))−1 (1 − ω)D − ωU x(k ) +

| {z }
Cr
((D + ωL))−1 (ωb)
| {z }
dr

Pour w = 1, on retrouve bien l’équivalence des formulations de Gauss

Seidel et de la relaxation
Condition suffisante de convergence
x(k +1) = Cx(k ) + d
Si xe est solution du problème, on a xe = C xe + d et donc
x(k +1) − xe = C(x(k ) − xe)
= C 2 (x(k −1) − xe)
...
= C k +1 (x(0) − xe)
Si maintenant on raisonne en norme:
erreur à l’étape k +1
z }| {
kx(k +1) − xek = kC k +1 (x(0) − xe)k
≤ kC k +1 kkx(0) − xek ≤ kCkk +1 kx(0) − xek
| {z } | {z }
→0? erreur initiale

ce faisant, nous avons démontré le résultat suivant :

Theorem (convergence d’une suite vectorielle)
S’il existe une norme matricielle telle que kCk < 1 alors la suite x(k )
converge vers xe = (I − C)−1 d
Normes matricielles de type vectorielle

Norme vectorielle (norme de Frobenius)

XX
kCk2F = Cij2
i j

Ce type de norme n’est pas toujours pertinent pour ce que nous

cherchons à faire. En effet, certaines normes vectorielles ne sont pas
sous multiplicatives (la propriété kABk ≤ kAkkBk n’est pas toujours
vérifié) ce qui est fâcheux pour notre analyse.
Par exemple pour la norme vectorielle kAk∆ = maxi,j |Aij |:
 
1 1
A=B= kABk∆ = 2  kAk∆ kBk∆ = 1
1 1
Normes matricielles d’opérateur
Pour conserver cette propriété, il faut utiliser une Norme d’opérateur
définie par:
kCxk
kCk = sup
x6=0 kxk

Par construction on a bien:

kABxk ≤ kAkkBxk ≤ kAkkBkkxk ⇒ kABk ≤ kAkkBk

Exemples (théorèmes):
X
kCk1 = max |Cij |
j
i
X
kCk∞ = max |Cij |
i
j
kCk2 =
√
max µi où µi est valeur propre de C >C.
i √
pour C symétrique : max |λi | λi = µi est val. propre de C.
i
démonstrations: de manière générale nous allons montrer les deux
inégalités suivantes
φ(C) ≤ kCk ≤ φ(C)
On commence par trouver une borne sur cette norme:
kCxk
kCk1 = sup
x6=0 kxk
P P
i j Cij xj
= sup
x6=0 P Pkxk
i j |Cij ||xj |
≤ sup
x6=0 kxk
X X
|xj |
≤ sup |Cij |
x6=0 j kxk
i
X X
|xj |
≤ sup max |Cij |
x6=0 j j kxk
i
X
X |xj |
≤ sup max |Cij |
x6=0 j kxk
i j
X |xj |
≤ max |Cij | car kxk
≤1
j
i

en choisissant le vecteur ej = (0, . . . , 0, 1, 0, . . . , 0)> , on atteint la borne

Le cas de la norme 2
on passe par le système des vecteurs propres: C > Cvi = µi vi
X X
x= ξi vi ⇔ x = Pξ ; kxk2 = ξ > P > Pξ = ξi2
i i
kCxk2
kCk22 = sup 2
x6=0 kxk
x> C > Cx
= sup
x6=0 kxk2
x> C > C
P

i ξi vi
= sup
x6=0 kxk2

x>
P
i ξi µi vi
= sup
x6=0 kxk2
P 2
ξ µi
= sup i i 2 ≤ max µi
x6=0 kxk i

Là encore, en choisissant le vecteur propre associé, on atteint la borne

:
kCvi k2
= vi C > Cvi = µi
kvi k2
Illustration 2d...
kAxk
kAk = sup = sup kAxk
x6=0 kxk kxk=1

Definition
rayon spectral

ρ(A) = max |λi |

1≤i≤n

Propriétés des normes (Op) Propriétés des normes

P
kAk > 0 si A 6= 0 kAk1 = maxj |Aij |
Pi
kAk = 0 ⇔ A = 0 kAk∞ = maxi j |Aij |
kkAk = |k |kAk A sym. : kAk2 = ρ(A)
kA + Bk ≤ kAk + kBk A sym. : ρ(A) ≤ kAk
p
kABk ≤ kAkkBk kAk2 ≤ |Ak1 |Ak∞
Résumons nous (normes matricielles?)
1 on cherche à résoudre Ax = b
2 on utilise une méthode itérative : x(k +1) = Cx(k ) + d
3 l’erreur... e(k ) = C k e(0)
4 ... peut être contrôlée : ke(k ) k ≤ kCkk ke(0) k
5 si l’on trouve une norme telle que kCk < 1 la suite converge
6 définir une norme telle que kC k k ≤ kCkk
kCxk
7 définition : kCk = supx6=0 kxk
8 théorèmes (bienX
pratiques...):
kCk1 = max |Cij |
j
i
X
kCk∞ = max |Cij |
i
j
√
kCk2 = max µi
i
9 on teste les différentes normes pour avoir kCk < 1

A quelles conditions sur A, la méthode itérative converge?

Le cas des matrices à diagonale dominante
Definition
Une matrice carrée A est a diagonale dominante si
n
X
. ∀i = 1, n |Aii | > |Aij |
j=1,j6=i

Theorem
Si A est a diagonale dominante alors les méthodes de Jacobi et de
Gauss Seidel convergent

Démonstration :
1 X
Jacobi : kCj k∞ = kD −1 (M + N)k∞ = max |Aij | < 1
i |Aii |
j6=i
Gauss Seidel: on montre que si Gauss Seidel diverge, alors A
n’est pas à diagonale dominante :
1 Pi−1 Pn
1 < maxi |Aii | j=1 |Aij | + j=i+1 |Aij |. Les détails de la preuve ne
sont pas évidents...
Exemple : la matrice est-elle à diagonale dominante?

  
  2 −1 0 0 −9 3 1 −1
3 −1 1 −1 2 −1 0  1 2 0 0
A1 = 1 −2 0 ; A2 = 
 0 −1 2 −1 ; A3 = −1 1
 
8 −5
1 1 4
0 0 −1 2 −2 1 4 6

Remarques :
toute matrice à diagonale dominante est régulière
on retrouve ce type de matrice dans les méthodes à éléments
finis .
Le cas des matrices positives

Definition (rappel)
Une matrice carrée A est définie positive si elle est symétrique et si si
. ∀x 6= 0 ∈ IRn x> Ax > 0

Theorem
Si A est strictement symétrique définie positive alors la méthode de
Gauss Seidel converge et la méthode de la relaxation pour 0 < ω < 2
étude des perturbations (ou des erreurs)
     
10 7 8 7 32 1
7 5 6 5  23 1
A=
8 ; b= admet comme solution x =  
6 10 9  33  1
7 5 9 10 31 1

   
32, 1 9, 2
22, 9 −12, 6
b + δb = 
33, 1
 admet comme solution x + δx =  
 4, 5 
30, 9 1, 1

Pour la norme infinie:

kδb k∞ kδx k∞
= 0, 003 = 13, 6
kbk∞ kxk∞

C’est la nature de la matrice A qui est en cause:

comment mesurer cet effet?
étude des perturbations (ou des erreurs)
     
10 7 8 7 32 1
7 5 6 5  23 1
A=
8 ; b= admet comme solution x =  
6 10 9  33  1
7 5 9 10 31 1

   
32, 1 9, 2
22, 9 −12, 6
b + δb = 
33, 1
 admet comme solution x + δx =  
 4, 5 
30, 9 1, 1

Pour la norme infinie:

kδb k∞ kδx k∞
= 0, 003 = 13, 6
kbk∞ kxk∞

C’est la nature de la matrice A qui est en cause:

comment mesurer cet effet?
Conditionnement d’un système linéaire
Ax = b A(x + δx ) = b + δb

1 kAk
kbk = kAxk ≤ kAk kxk ⇒ ≤
kxk kbk
kδx k = kA−1 δb k ≤ kA−1 k kδb k ⇒ kδx k ≤ kA−1 k kδb k

kδx k kδb k
≤ kA−1 kkAk
kxk kbk

cond(A) : conditionnement de A
|λ1 |
cond(A) = kA−1 kkAk =
|λn |

cond(A) proche de 1 : stable

cond(A) grand devant 1 : instable
Conclusion

3 méthodes itératives (Jacobi, Gauss Seidel, relaxation)

chaque itération ≤ O(n2 ) : produit matrice vecteur

encore mieux si A est sparse
peut s’appliquer à des matrices par blocs

convergent à certaines conditions sur A

pour comprendre : il faut connaitre la notion de norme matricielle

l’importance des valeurs propres

il existe des méthodes plus sophistiques

a chaque itération on se rapproche de la cible...