Master E . S . E . R .

Méthodes Numériques

2023-2024
Table des matières

1 Errurs et conditionnement 3
1.1 Erreur absolue et erreur relative . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1 Erreur absolue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.2 Erreur relative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Différentes erreurs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1 Sources des erreurs . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.2 Erreur totale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Sensibilité aux erreurs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.1 Problème bien posé . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Problème bien conditionnée . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2 Méthodes et algorithmes 9
2.1 Méthodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1 Différents types de méthodes . . . . . . . . . . . . . . . . . . . . . . 9
2.1.2 Méthodes exactes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.3 Méthodes approchées . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.4 Méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Algorithme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.1 Complexité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.2 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3 Méthode exacte pour la résolution de systemes linéaires 14

3.1 Problème posé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Formule de Cramer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3 Résolution des systèmes triangulaires supérieurs . . . . . . . . . . . . . . . 15
3.3.1 Algorithme et code : . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3.2 Coût de l’opération . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4 Factorisation LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4.1 Ecriture matricielle des transformations élémentaires : exemple . . . 16
3.4.2 Cas général . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.4.3 Condition de factorisation . . . . . . . . . . . . . . . . . . . . . . . 19
3.4.4 Implantation (version élémentaire) . . . . . . . . . . . . . . . . . . 19
3.5 Stratégie des pivots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.5.1 Factorisation PA=LU . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.5.2 Stratégie du pivot partiel . . . . . . . . . . . . . . . . . . . . . . . . 20
3.5.3 Implantation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.6 Matrices symétriques définies positives . . . . . . . . . . . . . . . . . . . . 21
3.6.1 Théorème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

1
TABLE DES MATIÈRES 2

3.6.2 La méthode de Cholevski . . . . . . . . . . . . . . . . . . . . . . . . 22

3.6.3 Implantation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.7 Etude du conditionnement de la résolution d’un système de CRAMER . . 22
3.7.1 Normes matricielles . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.7.2 Erreur sur le second membre . . . . . . . . . . . . . . . . . . . . . . 24
3.7.3 Erreur sur toutes les données . . . . . . . . . . . . . . . . . . . . . 24
3.7.4 Exemple de système mal conditionné . . . . . . . . . . . . . . . . . 25

4 Méthodes itératives pour la résolution de f(x)=0 27

4.1 Séparation des racines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2 Méthode de la dichotomie . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

5 Méthodes itératives pour la résolution des systèmes de CRAMER 29

Chapitre 1

Errurs et conditionnement

1.1 Erreur absolue et erreur relative

1.1.1 Erreur absolue
Définition 1.Si x est une valeur approchée de la valeur exacte x, la valeur x − x
s’apelle l’erreur absolue commise sur x. Elle peut être de signe quelconque.

Ne connaissant pas la valeur exacte, on ne connait pas l’erreur absolue. On cherche

donc à avoir un majorant de la la valeur absolue, quantité que l’on appelle incertitude
absolue
Dire que x est une valeur approchée de x avec une incertitude absolue de ² revient à dire
que l"on a :
|x − x| ≤ ² ⇐⇒ x−²≤x≤x+² (1)
Régle : L’incertitude absolue d’une somme ou d’une différence est égale à la somme des
incertitudes absolues de chacun des termes

1.1.2 Erreur relative

Définition 2. Si x est une valeur approchée de la valeur x non nulle, la valeur x−x x
s’appelle erreur relative commise sur x. Ne pouvant connaitre la valeur exacte x, on rem-
place souvant l’erreur relative par la valeur x−x
x
. Ces deux valeurs peuvent être de signe
quelconque.

Pour la même raison que précédemment on cherche à avoir un majorant de la valeur

absolue de l’erreur relative , quantité que l’on appelle incertitude relative.
Dire que x est une valeur approchée du réel x avec une incertitude relative ² revient à
dire que l’on a :
x−x
| |≤² ⇐⇒ x(1 − ²) ≤ x ≤ x(1 + ²) (2)
x
Régle pratique(obtenue en négligeant les termes du second ordre) : L’incertitude
relative d’un produit ou d’un quotient est égale à la somme des incertitudes relatives de
chacun des termes

3
CHAPITRE 1. ERRURS ET CONDITIONNEMENT 4

Fig. 1.1 – Variation de l’erreur totale

1.2 Différentes erreurs

1.2.1 Sources des erreurs
Lorsqu’on éffectue des calculs scientifiques plusieurs sources d’erreurs peuvent appa-
raitre.
Erreurs dues aux données : les données sont souvent des résultats expérimentaux ou
des valeurs approchées ; elles ne sont pas connues de manière exacte.
Erreurs due à la méthode : les méthodes employées sont souvent des méthodes ap-
prochées ; il y a alors une erreur systématique.
Erreurs dues aux arrondis : lors des calculs en arithmétique flottante, des arrondis
sont éffectués.

Exemple : cherchons à évaluer une valeur approchée de sin π7

Pour cela, en utilise le developpement en série entière de sin x (ce qui n’est pas une trés
bonne méthode...) :
∞
X x2n+1
sin x = (−1)n
n=0 (2n + 1)!
en remplaçant x par π7 .
Une erreur sur les données est commise en remplaçant π7 par une valeur approchée qui
est un flottant ; une erreur de méthode apparait en remplaçant la série par une somme
partielle et des erreurs d’arrondis sont commises dans les calculs des termes de la somme
à une précision donnée.

1.2.2 Erreur totale

C’est la somme des erreurs de sources différentes ; c’est à dire des erreurs de données,
des erreurs de méthode et des erreurs d’arrondi. Il est donc important que ces erreurs
soient du même ordre de grandeur, même si pour cela il faut faire des compromis.

Exemple : on veut calculer la dérivée d’une fonction f en un point x0 ; pour cela, on

emploie la formule :
f (x0 + h) − f (x0 − h)
f 0 (x0 ) =
2h
On montre qu’il existe une constante K > 0 telle que l’erreur de méthode est majorée par
Kh2 . En conséquence, plus h est petit, meilleure sera l’approximation, mais la dévision
par un nombre h proche de 0 va occasionner une erreur d’arrondi importante ; il y a donc
une valeur optimale h à évaluer ou, au moins, à estimer.
CHAPITRE 1. ERRURS ET CONDITIONNEMENT 5

1.3 Sensibilité aux erreurs

1.3.1 Problème bien posé
Tout problème de calcul numérique consiste en fait à déterminer une solution x dé-
pendant de données a :
x = F (a) (3)
Défintion 3.Un problème est dit bien posé au sens de HADAMARD1 s’il possède une
solution unique, continue par rapport aux données. Il est dit mal posé dans le cas contraire.

Il n’est pas trés facile de donner un contenu mathématique à ce qui précéde car cela
nécessite de définir des normes sur l’ensemble des données et sur l’ensemble des résultats.
Or ce sont souvent des ensembles de fonctions où il est délicat de définr de telles normes.

Exemples :
1. Résoudre le système linéaire Ax = b, où A ∈ MN (R) avec det(A) 6= 0 et b ∈ RN
sont des données et où x ∈ RN l’inconnue ; c’est un système de CRAMER2 . La
solution x unique s’exprime à l’aide des formules de CRAMER comme des fractions
rationnelles en les coefficients de la matrice A et les composantes du vecteur b : c’est
donc un problème bien posé.
2. Résoudre l’équation différentielle y 0 (x) = y(x) dans R avec la condition initiale
y(0) = α. La solution y(x) = αex est unique et dépend continûment de la donnée
α : cest donc un problème bien posé.
On peut remarquer sur ces deux exemples simples qu’il est ne serait pas immédiat de
formaliser la conclusion.

Contre-exemple :
Déterminer le nombre de solutions réelles de l’équation α − 12x(x − 1)2 = 0 avec α réel.
Selon la valeur de α on peut distinguer plusieurs cas :
– α < 0 : un zéro
– α = 0 : deux zéros
– 0 < α < 169
: trois zéros
Ne pouvant définir de notion de continuité pour une fonction à valeur dans N , le problème
est mal posé.

1.4 Problème bien conditionnée

Un problème x = F (a) étant bien posé il est souhaitable de savoir mesurer sa sensibilité
aux perturbations sur les données.

1
Jacques HADAMARD, mathématicien Français (versaille 1865- Paris 1963). Théorie des fonctions.
Théorie des nombres. Analyse fonctionnelle. Enseignant à l’Ecole Polytechnique et au Collège de France.
2
Gabriel CRAMER, mathématicien Suisse (Genève 1704 - Bagnol-Sur-Cèze 1752). Courbes planes.
Systèmes d’équations linéaires.
CHAPITRE 1. ERRURS ET CONDITIONNEMENT 6

Définition 4. Un problème bien posé au sens de HADAMARD est conditionné de

manière absolue s’il existe une constante C telle que pour tout a au voisinage de a, on
ait :
kF (a) − F (a)k ≤ Cka − ak (4)
La plus petite des constantes C vérifiant l’inégalité précédente est appelée conditionne-
ment absolue du problème.

De la même manière que l’on parle d’erreur absolue et d’erreur relative, on utilise aussi
souvent une autre notion de conditionnement :

Définition 5. Un problème bien posé au sens de HADAMARD est conditionné de

manière relative s’il existe une constante C telle que pour tout a au voisinage de a, on
ait :
F (a) − F (a) a−a
k k ≤ Ck k (5)
F (a) a
La plus petite des constantes C vérifiant l’inégalité précédente est appelée conditionne-
ment relative du problème.

On peut remarquer, comme pour les problèmes bien posés, la difficulté à définir cela
de manière correcte et la dépendance à l’égard des normes utilisées. Il ne sera pas toujours
facile de déterminer si un problème est bien posé, bien ou mal conditionné.
Selon les problèmes considérés et quand cela sera possible, on utilisera l’un ou l’autre de
ces deux conditionnements. Ce qui est important c’est l’ordre de grandeur davantage que
leur valeur précise qui de toute manière dépend des normes.
Si la costante C est petite la solution x sera donc peu sensible aux perturbations sur
les données : on dit que le problème est bien conditionné ; un algorithme de calcul de la
solution sera stable.
Si elle est grande les perturbations sur les données seront amplifiées : on dit que le pro-
blème est mal conditionné ou «raid» (sfiff). Il est à craidre qu’un algorithme de calcul de
la solution soit instable : le chois de cet algorithme sera alors important.
Le caractère «petit» ou «grand» dépend bien sûr des problèmes traités et du type de condi-
tionnement : on peut raisonnablement parler de problème mal conditionné si C ≥ 10 et
de problème bien conditionné si C < 2.

Exercice. On considère le système linéaire suivant :

à ! à ! à !
13 21 x1 34
. =
21 34 x2 55

1. Le résoudre d’une maniàre exacte.

2. Modifiér légèrement le second membre, puis calculer d’une manière exacte la solu-
tion.
3. Modifiér légèrement la matrice, puis calculer à nouveau la solution.
4. Que pensez-vous de ces résultats.
CHAPITRE 1. ERRURS ET CONDITIONNEMENT 7

Exemple de problème mal conditionné

On considère le problème de CAUCHY 3 suivant :
(
y 0 = 10(y − x2 )
y(0) = C
La solution exacte est :
x 1 1
y(x) = x2 + + + (C − )e10x
5 50 50
La détermination d’une valeur approchée de y(5), par exemple, est délicate. En effet,
elle dépend très fortement de la condition initiale : le conditionnement absolu du problàme
est égal à e50 ≈ 5, 18 × 1021 .

Exercice. Déterminer la valeur exacte de y(5), puis une valeur approchée obtenue par
une méthode numérique.

Quand il y a plusieurs données il est quelquefois plus simple de les traiter séparément
et de considérer des conditionnements relatifs par rapport à chacune de ces données.

Exemple :
Si x = F (a) ∈ R avec a = (a1 , a2 , ..., an ) ∈ Rn et F de classe C 1 , on peut écrire en
négligeant les termes du second ordre :
x − x = F (a) − F (a)
n
X ∂F
x−x≈ (a) (ai − ai )
i=1 ∂ai
n
x−x X ai ∂F (ai − ai )
≈ (a)
x i=1 F (a) ∂ai ai
n
|x − x| X ai ∂F (ai − ai )
≤ | (a)| | |
|x| i=1 F (a) ∂ai ai
En posant Ci = | Fa(a)
i ∂F
∂ai
(a)|, que l’on peut appeler conditionnement par rapport à la
donnée ai , on obtient l’inégalité :
|x − x| (ai − ai )
≤ Ci | | (6)
|x| ai
qui permet de mesurer l’effet des perturbations des données a = (a1 , a2 , ..., an ) sur x.
Dans le cas particulier où x = F (a1 , a2 ) = a1 + a2 , on obtient :
∆x a1 ∆a1 a2 ∆a2
| |≤| || |+| || |
x a1 + a2 a1 a1 + a2 a2
Si a1 + a2 est proche de zéro le calcul de la somme sera mal conditionné : on retrouve le
problàme de lélimination catastrophique ; dans les autres cas, il sera bien conditionné.
3
Augustin Louis CAUCHY, mathématicien français (Paris 1789 - Sceaux 1857). Fondements de l’Ana-
lyse où il introduit la rigueur connue actuellement. Equations différentielles. Fonctions de variable com-
plexe. Elasticité. Professeur au Collège de France et à l’Ecole Polytechnique. Il dut s’exiler à la révolution
de 1830 pendant plusieurs années en raison de ses opinions l"gitimistes, opinions qu’il gardera toute sa
vie.
Chapitre 2

Méthodes et algorithmes

2.1 Méthodes
2.1.1 Différents types de méthodes
On considère le problème de calcul numérique bien posé au sens de HADAMARD
x = F (a)
consistant à déterminer la solution x en fonction des données a.
On distingue plusieurs types de méthodes selon le résultat obtenu (valeur exacte de la
solution ou valeur approchée) et le moyen de l’obtenir :
méthodes directes : le résultat est donné par une formule et calculé en un nombre fini
d’opérations ; on distingue les méthodes exactes où le résultat est la valeur exacte de
la solution et les méthodes approchées où c’est une valeur approchée de la solution ;
méthodes itératives : la valeur exacte de la solution est la limite x = limn→∞ xn d’une
suite définie par récurrence et on prend une valeur xn avec n bien choisi comme
valeur approchée de la solution.

2.1.2 Méthodes exactes

En fait la valeur exacte de la solution n’est jamais obtenue : elle serait obtenue si les
calculs étaient fait en arithmétique exacte ce qui n’est pas le cas dans le cadre du calcul
numérique car ils sont, bien sûr, éffectués en arithmétique flottante.
L’erreur sur la valeur approchée obtenue n’est dûe qu’aux erreurs d’arrondi et éventuel-
lement aux erreurs sur les données ; par contre, il n’y a pas d’erreur de méthode.
On trouve parmi ces méthodes :
– le calcul de la valeur P (x) où P est un polynôme et x un nombre,
– la résolution d’un système linéaire par la méthode de CRAMER ou par la méthode
de GAUSS et ses variantes,
– le calcul du polynôme d’interpolation d’une fonction f en n points ;
– etc.

2.1.3 Méthodes approchées

La valeur xh choisie comme solution approchée de x est donnée par une formule ex-
plicite du type
xh = Ψ(a, h) (7)

8
CHAPITRE 2. MÉTHODES ET ALGORITHMES 9

où h est un paramètre réel «petit».

Dans le cadre de ces méthodes, on cherche à avoir une majoration de l’erreur de méthode
kx − xh k en fonction du paramètre h. Plus précisement :

Définition 6. On dit que l’erreur de la méthode approchée est d’ordre p s’il existe
une constante K > 0 telle que

kx − Ψ(a, h)k < Khp (8)

L’erreur sur la valeur approchée obtenue est dûe à l’erreur de méthode, aux erreurs
d’arrondi et éventullement aux erreurs sur les données.
On trouve parmi ces méthodes :
– les formules de dérivation numérique et d’intégration numérique,
– les méthodes de résolution des équations différentielles,
– les méthodes des différences finies pour la résolution des équations aux dérivées
partielles,
– etc.

2.1.4 Méthodes itératives

La valeur xn choisie comme solution approchée de x est donnée par la formule de
récurrence du type :
xn = Φ(xn−1)
ou plus généralement
xn = Φ(xn−1 , . . . , xn−p ) (9)
Dans tous les cas le calcul de xn nécessite donc la connaissance des valeurs x0 , x1 , . . . , xn−p .
Dans le cadre de ces méthodes on cherche à avoir une majoration de l’erreur de méthode
kx − xn k en fonction de l’erreur kx − xn−1 k. Plus précisement :

Définition 7. On dit que la méthode itérative est convergente d’ordre p s’il existe
une constante K > 0 telle que, pour n assez grand :

kx − xn k < Kkx − xn−1 kp (10)

Si p = 1 la constante K doit de plus verifiér K < 1.

Remarques
– dans le cas p = 1, on parle de convergence linéaire et dans le cas p = 2 de convergence
quadratique.
– si une méthode est d’ordre p on peut estimer qu’on multiplie approximativement le
nombre de décimales significatives par p à chaque itération.

Comme pour les méthodes approchées, l’erreur sur la valeur approchée obtenue est
dûe à l’erreur de méthode, aux erreurs d’arrondi et éventuellement aux erreurs sur les
données. La question importante dans ces méthodes est : quand s’arrêter ? Pour quelle
valeur de de n considère-t-on que xn est une «bonne» valeur approchée de la valeur exacte
x?
CHAPITRE 2. MÉTHODES ET ALGORITHMES 10

Le test attendu serait une évaluation de l’erreur absolue kx − xn k ou de l’erreur relative

kx−xn k
kxn k
mais, x étant inconnu, un tel test est impossible. Dans la pratique, on utilise
principalement deux tests :
– un test d’arrêt sur l’écart entre deux valeurs approchées successives :

kxn−1 − xn k ≤ ² (11)

– un test d’arrêt sur l’écart relatif entre deux valeurs approchées successives :
kxn−1 − xn k
≤² (12)
kxn k
Il faut bien réaliser que ces deux tests d’arrêt n’apportent en général aucune garantie sur
l’erreur absolue ou sur l’erreur relative. Selon la nature des problèmes il est quelquefois
possible de les améliorer.
On trouve parmi ces méthiodes :
– la résolution des équations non linéaires f (x) = 0 par la méthode de dichotomie, la
méthode de NEWTON1 ou la méthode des sécantes.
– la résolution des systèmes linéaires par la méthode de JACOBI, la méthode de
GAUSS-SEIDEL ou la méthode de relaxation,
– etc.

2.2 Algorithme
Une fois la méthode est choisie, il s’agit de la mettre en oeuvre.

Définition 8. On appelle algorithme une suite finie et ordonnée de règles à appliquer

à des données pour arriver à un résultat.
Une méthode se traduit donc par un algorithme ; mais il n’y a pas obligatoirement unicité
de cette mise en oeuvre. Il faut donc sans doute choisir. De même il faut quelquefois
choisir entre plusieurs méthodes. Parmi les critères de chois important on trouve :
la complexité mesurée en nombre d’opérations élémentaires nécessaires à l’execution de
l’algorithme,
la stabilité qui caractérise l’insensibilité de l’algorithme à la propagation des erreurs.
La complexité évoquée ici est une complexité temporelle ; on peut s’intéresser aussi à
la complexité spatiale qui caractérise la place mémoire nécessaire à la mise en oeuvre
informatique de l’algorithme.

2.2.1 Complexité
On mesure la complexité temporelle d’un algorithme par l’ordre de grandeur du
nombre d’opérations élémentaires nécessaire à son exécution. On appelle opérations élé-
mentaires les quatre opérations de base (addition, soustraction, multiplication et division)
en arithmétique flottante. On fait ainsi l’hypothèse que leur temps d’exécution est iden-
tique. . .ce qui est faux en général. De même, en général, on ne compte pas le nombre
1
Isaak NEWTON, astronome, mathématicien et physicien anglais (Woolsthorpe 1642 - Londres 1727).
Fondateur avec LEIBNITZ (avec lequel il aura une longue polémique), du calcul différentiel et intégral.
Classification des cubiques. Mécanique. Théorie de la gravitation. Mouvement des planétes. Réalisation
du premiet téléscope. Théorie des couleurs. Alchimie. A laissé une oeuvre exceptionnelle.
CHAPITRE 2. MÉTHODES ET ALGORITHMES 11

d’instructions informatiques (instructions d’entrées/sorties, test, appel de fonction, appel

de procédure. . .). Cela donne néanmoins un ordre de grandeur qui permet éventuelle-
ment de comparer des algorithmes différents pour un même problème ; cela fournit aussi,
pour un algorithme donné, un ordre de grandeur de l’augmentation du temps d’exécution
lorsque la taille du problàme augmente. Cet ordre de grandeur s’exprime avec la notation
O de LANDAU.2

Exemple : Algorithme de HORNER

Pour calculer P (x) = an xn + an−1 xn−1 + . . . + a1 x + a0 on dispose de deux algorithmes :
Algorithme classique Algorithme de HORNER3
valeur=a0 valeur=a0
do i=1 :n do i=n-1 :-1 :0
puissance=1 valeur=valeur*x+a(i)
do j=1 :i end
puissance=puissance*x
end
valeur=valeur+a(i)*puissance
end
Si on décompte les opérations élémentaires, on obtient :
– pour l’algorithme classique : n(n+1)
2
multiplications et n additions, soit une complexité
2
en O(n ),
– pour l’algorithme de HORNER : n multiplications et n additions, soit une complexité
en O(n)
Pour de grandes valeurs de n (n > 10) l’avantage de l’algorithme de HORNER est évident ;
on peut y ajouter sa simplicité.
Prenons un deuxième exemple avec la résolution d’un système linéaire de dimension n de
la forme Ax = b où A est une matrice inversible. La méthode de CRAMER conduit à une
complexité qu’on peut qualifier de déraisonnable : les formules de CRAMER nécessitent
(n + 1)(n! − 1) additions, (n2 − 1)n! multiplications et n divisions , soit une complexité
de l’ordre de O(n2 n!).
On verra que la méthode GAUSS 4 nécessite environ 23 n3 opérations élémentaires, soit
une complexité de l’ordre de O(n3 ). La aussi, le choix de l’algorithme est vite fait. . .

Exercice

1. Décompter le nombre d’opérations nécessaire avec la formule de CRAMER.

2. En déduire l’ordre de grandeur du temps d’exécution d’un système de dimension n =
20 sur un ordinateur effectuant 106 opérations par seconde, puis sur un ordinateur
effectuant 109 opérations/seconde.
2
Edmund LANDAU, mathématicien allemand (Berlin 1877 - Berlin 1938) Théorie des nombre.
4
Carl friedrich GAUSS, astronome, mathématicien et physicien allemand (Brunswick 1777 - Göttingen
1855) Géométrie des courbes et des surfaces. Géométrie non euclidienne. Mécanique céleste. Magnétisme
terrestre. Optique. Enfant prodige, surnommé le prince des mathématiciens, laissé une oeuvre trés origi-
nale et a ouvert la voie dans de nombreux domaines en y introduisant la rigueur et de nouvelles méthode.
CHAPITRE 2. MÉTHODES ET ALGORITHMES 12

2.2.2 Stabilité
Un algorithme est stable (on dit aussi numériquement stable) s’il est peu sensible à la
propagation des erreurs.
Dans certains domaines -résolution numérique des équations différentielles par exemple-
une définition précise de la stabilité sera donnée.
Il est évident qu’un problème très mal conditionné conduira à des algorithmes instables.
Mais un problàme bien conditionné ne se traduira pas obligatoirement par un algorithme
stable. Il convient d’être prudent.
Chapitre 3

Méthode exacte pour la résolution de

systemes linéaires

3.1 Problème posé

On veut résoudre le système de Cramer suivant :
−
→
A→
−
x = b (13)
A = (aij ) : matrice carrée de dimension (N × N ) inversible
→
−
x = (x1 , x2 , . . . , xN )T vecteur inconnu
→
−
b = (b1 , b2 , . . . , bN )T vecteur second membre

3.2 Formule de Cramer

det(A(i) )
xi = (14)
det(A)
−
→
A(i) : matrice obtenue en substituant dans A la colonne i par b
X
det(A) = a1σ(1) × . . . × a1σ(N ) × ε(σ)
σ∈PN

PN : ensemble des permutations de {1, 2, . . . , N }

ε(N ) = +1, −1 : signature de la permutation

Calculons le nombre d’opérations élémentaires nécéssires pour résoudre le

système :

– Complexité d’un déterminant : C(N ) = N ! − 1 + (N − 1)N ! ∼ N

| {z } | ∗
{zN}!
| {z }
additions multiplications oprations

– complexité globale : C(N ) = (N + 1)C(N ) = N (N + 1)! . C’est de l’ordre O(N 2 N !)

1 √
– Pour N = 100, N ! ∼ N N + 2 e−N π (formule de Stirling), C(N ) ∼ 9.4 × 10161
opérations élémentaires. Sur un ordinateur réalisant 109 opérations flottantes par
seconde (1G flop), la résolution du système durera 3.0 × 10145 années ! ! !

13
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES14

3.3 Résolution des systèmes triangulaires supérieurs

3.3.1 Algorithme et code :
     
t1,1 ... ... ... t1,N x1 b1
     
 0 t2,2 ... ... t2,N   x2   b2 
−
→      
T→
−
x = b ⇐⇒ 
 ... ... ∗
  ... =
  ... 

     
 ... ...   ...   ... 
0 0 ... 0 tN,N xN bN

On a det(T ) 6= 0 d’où ti,i 6= 0 (i=1,2,...,N)

Par substitution, on calcule les composantes successives de −

→
x en commençant par la
Nième équation
 b
 xN = t N
N,NPN
 (bi − j=j+1 ti,j xj )
xi = ti,i
i = N − 1, N − 2, . . . , 1 (15)

Commande MATLAB :
»x=b/T
Code
function x=Tsup(T,b)
% résout le système linéaire T*x=b
N=length(b)
x=zero(N,1)
x(N)=b(N)/T(N,N)
for i=N-1 :-1 :1
x(i)=(b(i)-T(i,i+1 :N)*x(i+1 :N))/T(i,i)
end

3.3.2 Coût de l’opération

– N divisions
– 2 + 4 + 6 + ... + 2(N − 1) = 2 ∗ N2 ∗ (N − 1) multiplications et additions
Soit une complexité de l’ordre O(N 2 )

Question : donner le code pour résoudre un système triangulaire inférieur.

3.4 Factorisation LU
On met A sous forme de produit de deux matrices régulières L et U :

A = LU

– L est une matrice triangulaire inférieure (Lower)

– U est une matrice triangulaire supérieure (Upper)
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES15

−
→
La résolution du système de Cramer A− →
x = b se réduit au cas trés simple de résolution
des systèmes triangulaires. En effet :
( −
→
−
→ → − → L−
→
y = b
A→
−
x = b ⇐⇒ L U −
| {zx} = b ⇐⇒
→
− U−
→
x =→−
y (16)
y

3.4.1 Ecriture matricielle des transformations élémentaires : exemple

     
2x1 + x2 + x3 = 1 2 1 1 x1 1
6x1 + 2x2 + x3 = −1 ⇐⇒ 
 6 2 1    
 .  x2  =  −1 
−2x1 + 2x2 + x3 = 7 −2 2 1 x3 7
| {z } | {z } | {z }
A −
→
x −
→
b

Avant de transformer le système changeons son écriture :A(1) −

→
x = b(1) avec A(1) = A
→
− (1)
et b = b.

Etape 1 :
     
A(2) (1, :) = A(1) (1, :) 2 1 1 x1 1
(2) (1) 6 (1)      
A (2, :) = A (2, :) − 2 A (1, :) =⇒  0 −1 −2  .  x2  =  −4 
A(2) (3, :) = A(1) (3, :) − −2
2
A(1) (1, :) 0 3 2 x3 8
| {z } | {z }
A(2) −
→(2)
b

Montrons que A(2) = E (1) A(2) . En Effet :

     
A(2) (1, :) A(1) (1, :) 0
 (2)   (1)   
 A (2, :)  =  A (2, :)  −  3 ∗ A(1) (1, :) 
A(2) (3, :) A(1) (3, :) −A(1) (1, :)
| {z } | {z }
A(2) A(1)

or      
0 0 0 0 0 0 0
   (1)  
0 0 
(1)  (1) (1)  (1)
 3 ∗ A (1, :)  =  3a11 3a12 3a13  =  3 ∗A
(1) (1) (1) (1)
−A (1, :) −a11 −a12 −a13 −1 0 0
d’où  
0 0 0
A(2) = (I − 
 3 0 0 
) ∗ A
(1)

−1 0 0
   
0 1
→
−1   −
→  
Introduisons l =  3  et e 1 =  0  premier vecteur de la base canonique de
−1 0
R3 . Nous avons  
0 0 0
→
− 1− 
→ T
l e1 = 3 0 0 

−1 0 0
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES16

ce qui permet d’écrire A(2) sous forme

−
→1 →T (1)
A(2) = (I − l −
e 1 ) A = E (1) A(1)
| {z }
E (1)

Etape 2
     
A(3) (1, :) = A(2) (1, :) 2 1 1 x1 1
A(3) (2, :) = A(2) (2, :) =⇒  0 −1 −2  .  x2  =  −4 
    

3
A(3) (3, :) = A(2) (3, :) − −1
A(2)
(2, :) 0 0 −4 x 3 −4
| {z } | {z }
A(3) −
→(3)
b

De la même manière que la matrice A(2) on peut mettre la matrice A(3) sous forme
A(3) = E (2) A(2) . En effet
   
0 0 0 0 0 0
  
A(3) = A(2) −  0 0 0  = A(2)
−  0 0 0 
∗A
(2)
(2) (2) (2)
−3a21 −3a22 −3a23 0 −3 0
   
0 0
−
→2 
En posant l =  0  et e 2 =  1 
 −
→ 
, on a
−3 0
−
→2 →T (2)
A(3) = (I − l −
e 2 ) A = E (2) A(2)
| {z }
E (2)

Notons U = A(3) la matrice triangulaire supérieure. Nous avons

U = E (2) E (1) A

Les matrices E (1) et E (2) sont inversibles. L’inverse de E (1) est

−
→1 →(T )
(E (1) )−1 = L(1) = I + l −
e1

En effet
→
− 1 →T ) −
→1 − T ) −
→1 →T −→1 →T
(I − l −e 1 )(I + l →
e 1 ) = I − l (−
e 1 l )−
e1 =I
De la même manière on note (E (2) )−1 = L(2) .
On obtient la matrice triangulaire inférieure en écrivant

L(1) L(2) U = LU = A

En explicitant :
    
1 0 0 1 0 0 1 0 0
   
L= 3 1 0  0 1 0 = 3 1 0 

−1 0 1 0 −3 1 −1 −3 0
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES17

3.4.2 Cas général

On pose :
A(1) = A
Etape 1

Si le pivot a11 6= 0, on calcule les termes de la matrice A(2) :


 (2) (1)

 a1j = a1j j = 1, . . . , N

 (2)

 ai1 = 0 i = 2, . . . , N
(2) (1) (1)

 aij = aij − li1 a1j i, j = 2, . . . , N

 (1)

 a
 li1 = i1(1)
a11
 
1 0 0 ... 0
 
 −l21 1 0 ... 0 
 
A(2) = E (1) A(1) avec E (1) = 
 ... ... ... 


 −lN −1,1 0 ... 1 0 

−lN 1 0 ... 0 1
" #
(1)
a11 , A(1) (1, 2 : N )
En notation MATLAB : A(2) =
0 , A(2) (2 : N, 2 : N
Etape 2
(2)
Si le pivot a22 6= 0, on renouvelle l’opération pour A(2) (2 : N, 2 : N ) ∈ MN −1 (R) au
moyen de la matrice E (2) (2 : N, 2 : N ) construite de façon analogue à E (1) . On complète
E (2) pour en faire une matrice de MN (R) telle que le produit à gauche laisse invariante
la première ligne et la première colonne.
" #
1 , 0 →
− 2 →T
E (2)
= =I− l −e2
0 , E (2) (2 : N, 2 : N
−
→2
où −
→
e 2 est le 2ème vecteur de la base canonique de RN et l = (0, 0, l32 , l42 , . . . , lN 2 )T
(2)
ai2
avec li2 = (2)
a22
 
(1) (1)
a
 11
a12 A(1) (1, 3 : N ) 
A(3) = E (2) A(2) = E (2) E (1) A(1) 
= 0 (2)
a22 A(2) (2, 3 : N ) 

0 0 A(3) (3 : N, 3 : N )
(k)
avec A(3) (3 : N, 3 : N ) ∈ MN −2 (R) : Si au dela, les akk 6= 0, on pourra réitérer l’opération
pour obtenir finalement

A(N ) = E (N −1) . . . E (2) E (1) A(1) = U

La factorisation A = LU cherchée se déduit alors du lemme suivant :

−
→k → T
Lemme1 : Soit E (k) la mlatrice définie pour 1 ≤ k ≤ N − 1 par E (k) = I − l − e k,
→
−k
avec l = (l1k , . . . , lN k )T telle que lik = 0 si i ≤ k, et −
→
e k le k-ième vecteur de la base
N
canonique de R . Alors,
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES18

→
− k− T
– (i) la matrice E (k) est inversible L(k) = I + l → e k,
– (ii) le produit L = L(1) ...L(N −1) s’assemble sans calcul de telle façon que

 0
 sii < j
L(i, j) = 1 sii = j


lij sii > j

→
− k− T →
− k− T −
→k →T −→k − T → −
→k P
En effet, (I − l →
e k )(I + l →
e k ) = I − l (−
e k l )→
e k et −
e Tk l = N
i=k+1 δki lik = 0

−
→1 →T −
→N −1 → −
→1 − T →
− N −1 −
D’autre part , (I + l −
e 1 )...(I + l −
e TN −1 ) = I + l →
e 1 +...+ l →
e TN −1 puisque
→
−i
les produits −
→
e Ti l sont nuls dés que j ≥ i.

Partant de E (N −1) ...E (1) A = U , on obtient donc :

A = L(1) ...L(N −1) U = LU

C’est la factorisation cherchée, avec L triangulaire inférieure d’éléments diagonaux égaux

(j)
aij
à 1, et dont les éléments sous-diagonaux sont lij = (j) , 1 ≤ j < i ≤ N , et U triangulaire
ajj
supérieure inversible puisque sa diagonale est constituée des pivots qui sont tous différents
de 0.

3.4.3 Condition de factorisation

On admettra le théorème suivant :
Théorème 1 : Soit A = (aij )1≤i,j≤N ∈ MN (R) ; si pour tout k, 1 ≤ k ≤ N , la sous
matrice principale A = (aij )1≤i,j≤k de A est telle que det(Ak ) 6= 0, il existe une matrice
triangulaire inférieure L, dont les éléments diagonaux sont égaux à 1, et une matrice
triangulaire supérieure inversible U telles que A = LU . Cette factorisation est unique.
Cette condition nécéssaire et suffisante n’est pas trés satisfaisante car elle nécéssite le
calcul de déterminants.
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES19

3.4.4 Implantation (version élémentaire)

Code
function [L,U]=Gauss(A) ;
% Gauss factorise A=LU par la méthode d’élimination.
% L’appel se fait selon
% »[L,U]=Gauss(A)
[N, N ] = size(A) ;
for k=1 :N-1
if A(k,k)==0
error(’Matrice non factorisable....’)
end
indices=k+1 :N
% Colonne k de la matrice L :
A(indices,k)=A(indices,k)/A(k,k) ;
%Mise à jour de la matrice A(k) :
A(indices,indices)=A(indices,indices)-A(indices,k)*A(k,indices) ;
end
L=tril(A,-1)+eye(N) ;
U=triu(A) ;

Coût de l’opération
Pour chaque valeur de k :
– (N − k) divisions
– (N − k)2 multiplications et (N − k)2 soustractions pour la mise à jour de A

N −1
X N (N − 1) N2
(N − k) = ∼
k=1 2 2
N
X −1
2 1 2
2 (N − k)2 = N (N − )(N − 1) ∼ N 3
k=1 3 2 3
Le coût est d’environ : C(N ) ∼ 23 N 3
Retour à l’exemple  
2 1 1

Observons sur la matrice A =  6 2 1  les différentes étapes de la fonction ci-
−2 2 1
dessus ; par souci de lisibilité, nous imprimons en caractères gras les valeurs calculées de
la matrice L et leur positionnement dans  la matrice A:
2 1 1
 
– k=1 : A est remplacée par A =  3 −1 −2 
-1 3 2
 
2 1 1
– k=2 : A est remplacée par A =   3 −1 −2 


-1 -3 −4
Les instructions MATLAB
L=tril(A,-1)+eye(N)
U=triu(A)
permettent bien de retrouver L et U.
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES20

3.5 Stratégie des pivots

3.5.1 Factorisation PA=LU
(k)
Si aprés la k-ième étape le pivot akk = 0, comme la matrice A(k) (k : N, k : N ) est
(1) (k−1)
régulière (det(A) = a11 ...ak−1k−1 det(Ak (k : N, k : N )) 6= 0), il existe i > k tel que
(k)
aik 6= 0. On permute alors les lignes k et i avant de poursuivre l’algorithme. La matrice
de permutation s’écrit Pk = I − − →
u→−
u T , avec →
−
u =−
→ek−− →e i.
Pour tout j < k, on aura

−
→j →T −
→gj T
Pk E (j) Pk = I − Pk l −
ej =I− l −
→
e j = Eg
(j)

→
−gj −
→j
Le vecteur l s’obtient du vecteur l en permutant les termes lkj et lij . Puis on a :

Pk E (k−1) ...E (1) = (Pk E (k−1) Pk )(Pk ...Pk )(Pk E (1) Pk )Pk = Eg g
(k−1) ...E (1) P
k

En introduisant des matrices Pk à chaque étape, avec éventuellement Pk = I lorsqu’aucune

permutation n’est nécéssaire, on obtiendra :
g−2) g
(E (N −1) E (N ...E (1) ) (PN −1 ...P1 ) A = U,
| {z }
P

puis
P A = LU. (17)

3.5.2 Stratégie du pivot partiel

La factorisation P A = LU est nécessaire lorsqu’on rencontre un pivot nul. Il se peut
que l’on rencontre un pivot qui sans être nul est trés petit. Dans ce cas , l’effet des erreurs
d’arrondi peut-être catastrophique.
(k) (k)
En effet, comme akk intervient en dénominateur dans les calculs, |akk | a interêt à être
choisi assez grand pour des raisons de stabilité numérique. A l’étape k, plutôt que de
(k)
rechercher le premier coefficient non nul parmi les aik , i > k, on recherche l’indice im tel
que
(k) (k)
aim k = maxk≤i≤N {|aik |}
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES21

3.5.3 Implantation
Code
function [L,U,permute]=Gausspiv(A) ;
% Gauss calcule la factorisation LU de A avec la stratégie du pivot partiel.
% L’appel se fait selon :
% » [L,U,permute]=Gauss(A)
% Elle est telle que A(Permute, :)=LU.
% La solution du système Ax=b est donnée par :
% »x=(b(Permute)/L)/U
[N,N]=size(A) ;
Permute=1 :N
for k=1 :N-1

% Recherche du pivot, permutation des lignes correspondantes

[maxi,ligne]=max(abs(A(k :N, :))) ;
ligne=ligne+k-1 ;
A([k,ligne], :)=A([ligne,k], :) ;
Permute([k,ligne])=Permute([ligne,k]) ;

if A(k,k)==0.
error(’Matrice non inversible...’)
end

% colonne k de la matrice L :
A(k+1 :N,k)=A(k+1 :N,k)/A(k,k) ;

%Mise à jour de la matrice A(k)

A(k+1 :N,k+1 :N)=A(k+1 :N,k+1 :N)-A(k+1 :N,k)*A(k,k+1 :N) ;

end
L=tril(A,-1)+eye(N) ;
U=triu(A) ;

3.6 Matrices symétriques définies positives

On suppose que A est une matrice symétrique (A = AT ) et définie positive (−
→
x T A−
→
x >
−
→ N −
→ →
− −
→ T →
− −
→ →
−
0 ∀ x ∈ R , x 6= 0 et x A x = 0 ⇔ x = 0 ) donc réguliére.

3.6.1 Théorème
Théoràme 2.Toute matrice d’ordre N, symétrique et définie positive est factorisable
sous la forme LLT où L est une matrice triangulaire inférieure inversible. De plus, la
décomposition est unique, quand on choisit les éléments diagonaux de L piositifs.
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES22

3.6.2 La méthode de Cholevski

Le calcul des termes de L se fait par identification. on obtient :
( √
l11 = q a11 li1 = ai1 /l11 i = 2...N
Pj−1 2 Pj−1
ljj = ajj − k=1 ljk lij = (aij − k=1 lik ljk )/ljj j = 2 . . . N, i = j + 1, . . . N (18)

Le calcul s’effectuant colonne par colonne de la 1ère à la N-ième. Le nombre d’opérations

3 2
élémentaires est de l’ordre de N6 additions et multiplications, N2 divisions et N racines
carrées. dans le cas des matrices symétriqus définies positives il y a avantage à employer
la méthodede cholevski, de préférence à lméthode de Gauss.

3.6.3 Implantation
Code
function [L]=Cholevski(A) ;
% Cholevski calcule la factorisation LLT de A qui doit être symétrique et définie positive.
% L’appel se fait selon :
% » [L]=Cholevski(A)
% La solution du système Ax=b est donnée par :
0
% »x = (b/L)/L
N=size(A,1) ;
L=zeros(size(A)) ;
L(1,1)=sqrt(A(1,1)) ;
for k=2 :N

% Clcul des termes sous la diagonale

for m=1 :k-1
0
L(k, m) = (A(k, m) − L(m, 1 : m − 1) ∗ L(k, 1 : m − 1) )/L(m, m) ;
end ;
% Clcul des termes diagonaux
0
L(k, k) = sqrt(A(k, k) − L(k, 1 : k − 1) ∗ L(k, 1 : k − 1) ;
end ;

3.7 Etude du conditionnement de la résolution d’un

système de CRAMER
L’étude du conditionnement (4) et (5) nécessite l’utilisation d’une norme. Nous allons
rappeler ici les normes matricielles induites par les normes vectorielles et non induites.

3.7.1 Normes matricielles

Définnition 9.On appelle norme matricielle une application | k . k | : MN (K) −→ R
telle que :
– (i) | k . k | est une norme,
– (ii) ∀A, B ∈ MN (K), | k AB k | ≤ | k A k |.| k B k |

Normes matricielles induites.

CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES23

Soit A ∈ MN (K), et k . k une norme sur RN . On définit la norme matricielle induite

par cette norme "vectorielle" en posant :
−→ kA→
−xk
| k A k | = supk−
→ kA x k = sup −
→ →
− { }
x k=1 x 6= 0 k− →
xk
−
→
– Cette quantité est bien définie ; il existe au moins un − →x 0 6= 0 tel que | k A k
|.k−
→
x 0 k = kA−→
x 0k
– Par ailleurs, on verifie pour pour tout −→x ∈ RN , on a :
kAB −→x k ≤ | k A k |.kB −
→
x k ≤ | k A k |.| k B k |.k−
→
xk
de sorte que | k . k | définie une norme matricielle.
Normes induites par les normes vectorielles
Quelle est la norme matricielle | k . k |∞ induite par la norme vectorielle k.k∞ ?
Pour tout i, 1 ≤ i ≤ N , on vérifie que :
N
X N
X N
X
| aij xj | ≤ |aij | |xj | ≤ ( |aij |) max |xj |
1≤j≤N
j=1 j=1 j=1
N
X
kA→
−
x k∞ ≤ max ( |aij |) kxk∞
1≤i≤N
j=1
N
X
| k A k |∞ ≤ max |( |aij |)|
1≤i≤N
j=1
PN PN
Soit i0 un indice tel que max1≤i≤N ( j=1 |aij |) = j=1 |ai0 j | ; on peut choisir →
−
y ∈ RN tel
que : (
−1 si ai0 j < 0
yj =
+1 sinon
P P
de sorte que N
j=1 ai0 j yj =
N
j=1 |ai0 j |, et l’égalité est réalisée.
De la même manière on peut montrer que

X N X N
k−
→
x k1 = |xi | −→ | k A k |1 = max { |aij |} (19)
1≤j≤N
i=1 i=1
N
X
k−
→ 1 1
x k2 = ( x2i ) 2 −→ | k A k |2 = ρ(AT A) 2 (20)
i=1
X N
k−
→
x k∞ = max {|xi |} −→ | k A k |∞ = max { |aij |} (21)
1≤i≤N 1≤i≤N
j=1

ρ(A) étant le rayon spectral de la matrice A. Il est égal à : ρ(A) = max1≤i≤N |λi |, où les
λi sont les valeurs propres.

Normes matricielles non induites

Il existe d’autres normes matricielles utilisées en analyse numérique. Nous citons ici la
plus connue la norme de Frobenius, ou norme de la trace, définie par :
N
X 1
| k A k |F = ( |aij |2 ) 2 (22)
i,j=1
1
Exercice. ontrer que | k A k |F = (tr(AT A)) 2
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES24

3.7.2 Erreur sur le second membre

Supposons que le calcul de la solution est fait par un algorithme quelconque implanté
−
→ −
→
sur ordinateur. Notons →−x = − →
x v + δx la solution obtenue où δx représente l’ erreur
d’arrondi. On va supposer que →
−x est la solution exacte dans RN d’un problème perturbé
−
→ →
− −
→
A(−→
x + δx) = b + δb . On a ainsi choisi l’ensemble des perturbations admissible, qui ne
portent que sur le second membre ; d’autres chois sont possibles. Il vient :
−
→ →
−
δx = A−1 δb
On choisit une norme quelconque sur RN , notée k.k, et on note | k . k | la norme matricielle
induite. On a alors :
−
→ −
→
kδxk ≤ | k A−1 k |k δbk
−
→
puis en utilisant k b k ≤ | k A k |.k−
→x k, on obtient :
−
→ −
→
kδxk | k A−1 k |.k δbk
≤ −
→
| k A k |.k→
−xk kbk
−
→ −
→
kδxk −1 k δbk
≤ | k A k |.| k A k | −
→ (23)
k−
→xk kbk
En comparant avec (5), on déduit le conditionnemnt γ(A) = | k A−1 k |.| k A k | de la
matrice A, relativement à la norme choisie, pour la résolution des systèmes linéaires.

3.7.3 Erreur sur toutes les données

−→ −
→ − →
(A + ∆A)(−
→
x + ∆x) = b + ∆b
En supposant que ∆A est suffisamment petit pour que A + ∆A reste inversible. Comme
A est supposée inversible, on peut écrire :
A + ∆A = A(I + A−1 ∆A) (24)
1
Supposons en outre que, pour une norme matricielle induite, | k ∆A k | < |kA−1 k|
. On en
déduit que | k ∆A k |.| k A−1 k | < 1.

Lemme 2. Si B ∈ MN (R) est telle que | k B k | < 1 pour une norme matricielle
induite, alors (I − B) est inversible, et vérifie :
∞
X
−1
(I − B) = B∞1
k=1

1
|kI −B k|<
1−|kB k|
D’aprés ce lemme, la matrice (I +A−1 ∆A) est alors inversible ; et donc la matrice A+∆A.
En l’appliquant à (24), on obtient :
1
| k (A + ∆A)−1 k | ≤ | k A−1 k |.| k (I + A−1 ∆A)−1 k | ≤ | k A−1 k |
1 − | k A−1 ∆A k |
1
Généralisation du résultat bien connu sur les séries numériques géométriques
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES25

−→ −
→
En remarquant que ∆x = (A + ∆A)−1 (−∆A−
→
x + ∆b)on obtient :

−→ | k A−1 k | −→ −
→
k∆xk ≤ −1
(| k ∆A k |.k∆xk + k∆bk)
1 − | k A ∆A k |
−→ −
→
k∆xk | k A−1 k | k∆bk
≤ (| k ∆A k | + −
→ | k A k |)
k−→
xk 1 − | k A−1 ∆A k | kbk
−→ −→
k∆xk | k ∆A k | k∆bk
≤ γ(A)( + → − )(1 + O(| k ∆A k |)
k→−
xk |kAk| kbk

3.7.4 Exemple de système mal conditionné

à ! à ! à !
1000 999 x1 1999
. =
999 998 x2 1997
à !
−1 −998 999
Calculons la matrice inverse : A = . Nous avons | k A k |∞ = 1999,
999 −1000
| k A−1 k |∞ = 1999 et γ(A) = 19992 . C’est une matrice trés mal conditionnée ; une petite
−
→
perturbation des données (A ou/et b ) donne un solution trés éloignée de la solution
initiale. Ã !
−
→ 1
En effet, la solution du système ci-dessus est u = , ajoutons au second membre
1
à !
→
− →
− −2 −
→ −1
b une perturbation δb = 10 v avec v = . On obtient le nouveau système
1
à ! à !
−
→ →
− −
→ −
→ 1 19, 97
A(−→x + δx) = b + δb dont la solution est −
→
x + δx = →−
u +10−2 A−1 →−v = +
1 −19, 99
−
→ −
→
. Remarquons le trés grand écart entre les valeurs de − →x et de δx. On a kδx k∞
→ = 19, 99
−
k x k∞
−
→ −
→ →
−
pour k− δb k∞ 10−2 kδxk∞ δb k∞
= 19992 k→
→ = 1999 soit k− →
x −
k b k∞ k ∞ k b k∞

Interpretation
Le système précédent est sensible aux perturbations. On sait qu’on obtient un système
équivalent en multipliant chaque ligne par un scalaire non nul. On peut donc le remplacer
par le système obtenu en remplaçant chaque équation (Eq : i) par a1i1 (Eq : i)
à ! à ! à !
1 0.999 x1 1.999
. =
1 0.9989998 x2 1.998998

On voit ici que les lignes A(1, :) et A(2, :) sont presque dépendantes puisque

A(1, :) − A(2, :) = [0, 2.010−6 ]

A vrai dire, les colonnes de la matrice A sont aussi presque dépendantes et on peut vérifier
que
A(:, 1) − A(:, 2)/A(1, 2) = [0, 2.010−6 ]
Lorsque les lignes ou les colonnes d’une matrice ne sont pas indépendantes, cette ma-
trice est non inversible. On peut donc penser, que d’une certaine façon, c’est parce que la
CHAPITRE 3. MÉTHODE EXACTE POUR LA RÉSOLUTION DE SYSTEMES LINÉAIRES26

matrice A est proche d’une matrice singulière que le système est très sensible aux pertur-
bations.

Retour sur la factorisation A = LU

On a vu précédemment que la factorisation A = LU s’éffectue de préférence sous la

forme P A = LU , en adoptant la stratégie du pivot partiel.
On peut mieux comprendre à présent que si à l’étape k de la factorisation, le pivot akkk
est très petit, les coefficients lik , k ≤ i ≤ N seront en général très grands . Les lignes de
la matrice A(k+1) seront défines par l’expression MATLAB

A(i, k + 1 : N ) − l(i, k) ∗ A(k, k + 1 : N )

Si les lik sont très grands, le premier terme sera négligeable, et ces lignes seront presque
dépendantes, puisqu’à peu près proportionnelles à A(k, k + 1 : N ) ! En terme de méthode
−
→(k+1)
d’élimination, le système A(k+1) −→
x = b sera très sensible aux petrturbations, puisque
proche d’un système non inversible.

3.8 Exercice
3.8.1 Exercice
On veut résoudre le système linéaire Ax = b, où
     
1 1 1 1 x1
A= 2 2 5 

 , b= 2 

 et x =  x2 


4 6 8 5 x3

par la méthode d’élimination de Gauuss.

1. Vérifier que l’algorithme de gauss ne peut pas être éxécuté jusqu’au bout.
2. On considère la matrice de permutation P suivante :
 
1 0 0
P = 0 0 1 


0 1 0

Ecrire le système équivalent à Ax = b (c.-à-d. ayant la même solution x) qui a P A

comme matrice associée.
3. Appliquer l’algorithme de Gauss à la matrice P A, et calculer la factorisation LU de
P A.
4. Calculer x à partir de la factorisation LU en résolvant les deux systèmes triangu-
laires.

3.8.2 Exercice
Pour résoudre le système A2 x = b, connaissant la matrice A, vaut-il mieux calculer
puis factoriser A2 , ou procéder d’une autre façon.
Chapitre 4

Méthodes itératives pour la résolution

de f(x)=0

On considère l’équation f (x) = 0 où f est continue et on suppose l’existence et l’unicité

d’une solution x∗ dans l’intervalle [a, b]. Différentes méthodes peuvent être utilisées pour
trouver une valeur approchée de x.

4.1 Séparation des racines

Définition 10. Une racine x∗ est separée dans l’intervalle [a, b] si cet intervalle
contient uniquement la racine x∗ .

Pour localiser les racines on peut utiliser l’une des méthodes suivantes :
méthode graphique : on localise les zéros en considèrant les points ou le graphe y =
f (x) coupe l’axe des abscisses. Il peut être commode d’écrire f (x) = G(x) − H(x) ;
les racines sont alors à l’intersection des graphes de y1 = G(x) et y2 = h(x).
Méthode algébrique : on utlise un corollaire du théorème de Rolle, sous la forme des
conditions suffisantes suivantes :

f (a) ∗ f (b) < 0 f (x) a au moins une racine dans ]a, b[

f (x) a une seule racine si f 0 (x) 6= 0
f (a) ∗ f (b) > 0 f (x) n’a pas de racine ou un nombre pair de racines dans ]a, b[
méthode numérique : on éffectue une tabulation de la fonction sur [a, b] pour détecter
les intervalles tels que f (xi ) ∗ f (xi+1 ) < 0. On aura ensuite une estimation de la
racine en faisant un resserrement de l’intervalle [xi , xi+1 ].

4.2 Méthode de la dichotomie

4.2.1 Déscription de la méthode
On suppose que f (a)f (b) < 0 et que f est monotone dans [a, b]

La méthode de dichotomie consiste à créer une suite {xk } qui converge vers la racine
x∗ . On procède de la manière suivante :

27
CHAPITRE 4. MÉTHODES ITÉRATIVES POUR LA RÉSOLUTION DE F(X)=0 28

Fig. 4.1 – Résolution d’une équation non linéaire

a+b
x0 = 2
<0 x∗ ∈]a, x0 [, on définit a1 = a et b1 = x0
si f (a)f (x0 ) = 0, x∗ = x0
>0 x∗ ∈]x0 , b[, on définit a1 = x0 et b1 = b
Dans chacun des cas x∗ ∈]a1 , b1 [ tel que |a1 −b1 | = 12 |b−a|. On définit de même x1 = a1 +b1
2
etc . . .En poursuivant la dichotomie on forme la suite :
a0 + b0 a1 + b1 an + bn
x0 = , x1 = , . . . , xn = ,...
2 2 2

4.2.2 Convergence de la suite

Nous avons :
b−a b−a
|an − bn | = n
et |xn − x∗ | ≤ n+1
2 2
c’est une suite qui converge vers la racine. La convergence est lineaire (voir (10)) puisque
|xn − x∗ | ≤ 12 |xn−1 − x∗ |

4.2.3 Arrêt de l’algorithme

|an −bn |
Si l’on choisit x∗ = an +b
2
n
, l’incertitude sur la racine sera ² = 2
. Pour une incerti-
tude ² fixée on arrête l’algorithme avec le test :

|an − bn |
≤² (25)
2

4.3 Méthode de la sécante

4.3.1 Déscription de la méthode
La valeur xn+1 est obtenue de la manière suivante : xn+1 est l’abscisse du point d’in-
tersection de la sécante joignant les deux points (xn , f (xn )) et (xn−1 , f (xn−1 )) avec l’axe
Ox(FIG.4.1). Quand les valeurs de départ x0 et x1 sont proches de la valeur cherchée, la
suite (xn ) converge vers x∗ . De manière plus précise, on a :
CHAPITRE 4. MÉTHODES ITÉRATIVES POUR LA RÉSOLUTION DE F(X)=0 29

Théorème 3. On suppose f de classe C 2 sur l’intervalle ]a, b[ ; on suppose aussi qu’il

existe une valeur x∗ telle que f (x∗ ) = 0 et f 0 (x∗ ) 6= 0. Alors il existe un réel δ > 0 tel que
la suite définie par :
(
x0 et x1 quelconques appartenant [x∗ − δ, x∗ + δ],
xn+1 = xn − f (xn ) f (xxnn)−f
−xn−1
(xn−1 )
pour n ≥ 0, (26)

soit bien définie et converge vers x∗ .

4.3.2 Convergence de la suite

Pour n assez grand, on a la majoration suivante :
√
∗ ∗ p 1+ 5
|x − xn | ≤ |x − xn | avec p=
2
√
La convergence est donc d’ordre p = 1+2 5 ≈ 1, 618.
Dans la pratique, il est rare de savoir si l’on se trouve dans le voisinage [x∗ − δ, x∗ + δ]
dont l’existence est assurée par le théorème précédent. La suite (xn ) peut alors diverger. . .

4.3.3 Arrêt de l’algorithme

Théoriquement, la solution n’est atteinte qu’aprés une infinité d’itérations (et si la
méthode converge). En pratique, on arrête les opérations par l’un des tests suivants :
1. si f (xn ) est quasiment nulle :
f (xn ) ≤ ²1 (27)

2. si l’amélioration de xn d’une itération à la suivante ne justifie pas l’effort de calcul

supplémentaire, soit si :
xn − xn−1
| | ≤ ²2 ou |xn − xn−1 | ≤ ²3 (28)
xn
3. si la convergence n’est pas obtenue avant un nombre d’itérations pédéfini, soit si :
n ≥ nmax (29)

4.4 Méthode du point fixe

On suppose que l’équation de départ f (x) = 0 est mise sous la forme équivalente :
x = F (x) (30)
Exemple

f (x) = x2 + 3ex − 12 peut s’écrire :

2 x
x = F1 (x) = x
√ + 3e − 12 + x
x = F2 (x) = 12 − 3e x
2
x = F3 (x) = ln 12−x
3
CHAPITRE 4. MÉTHODES ITÉRATIVES POUR LA RÉSOLUTION DE F(X)=0 30

4.4.1 Déscription de la méthode

La méthode consiste à utiliser un estimé x0 , de la solution exacte x∗ qui vérifie f (x∗ ) =
0 et donc x∗ = F (x∗ ), cette valeur estimée étant substituée à x dans le terme de droite
de l’équation (30). On obtient une nouvelle approximation x1 de x∗ . Par substitutions
successives on construit la suite (FIG.4.1) :




x0


 x1 = F (x0 )

... (31)



 xn = F (xn−1 )


 ...

4.4.2 Convergence de la suite

Théoràme 4. Si |F 0 (x)| ≤ m < 1 lorsque x ∈ [a, b], la suite (31) est convergent et
tend vers x∗ solution de x = F (x) donc de f (x) = 0. De plus la racine est unique.

En effet, on obtient par le théorème des accroissements finis :

F (xn1 ) − F (x∗ ) = xn − x∗ = F 0 (c)(xn−1 − x∗ ) o c ∈]xn−1 , x∗ [

Désignons par m = sup |F 0 (x)| lorsque x ∈]a, b[, m < 1, nous avons :


 |xn − x∗ | ≤ m|xn−1 − x∗ |

 ∗ ∗
|xn−1 − x | ≤ m|xn−2 − x |


 ...

|x1 − x∗ | ≤ m|x0 − x∗ |
D’où
|xn − x∗ | ≤ mn |x0 − x∗ |
Etant donné ² > 0, il suffit de choisir :
ln ² − ln(b − a)
n> (32)
ln m
pour obtenir |xn − x∗ | ≤ ². D’où la convergence de l’algorithme.

4.4.3 Ordre de l’algorithme

Développons en série de Taylor F (xn ) autour de x∗ :
(xn − x∗ ) 0 ∗ (xn − x∗ )2 00 ∗
F (xn ) = F (x∗ ) + F (x ) + F (x ) + . . .
1! 2!
Soit :
(xn − x∗ ) 0 ∗ (xn − x∗ )2 00 ∗
xn+1 − x∗ = F (x ) + F (x ) + . . .
1! 2!
D’où en désignant par ²n = xn − x∗ , on obtient :

0 ²2n 00 ∗
∗
²n+1 = ²n F (x ) + F (x ) + . . . (33)
2!
Donc si :
CHAPITRE 4. MÉTHODES ITÉRATIVES POUR LA RÉSOLUTION DE F(X)=0 31

– F 0 (x∗ ) 6= 0 et ²n assez petit pour être negligé :

²n+1 ≈ ²n F 0 (x∗ ) (34)

on parle alors de convergence linéaire.

– Mais si F 0 (x∗ ) = 0 et F 00 (x∗ ) 6= 0 alors :

²2n 00 ∗
²n+1 ≈ F (x ) (35)
2
On parle de convergence quadratique, et ainsi de suite pour les ordres plus élevés.

4.5 Méthode de Newton Raphson

4.5.1 Déscription de la méthode
Le passage de xn à xn+1 se fait de la manière suivante : xn+1 est l’abscisse du point
d’intersection de la tangente au point (xn , f (xn )) avec l’axe Ox (FIG.4.1). Quand la valeur
de départ x0 est proche de la valeur cherchée, la suite (xn ) converge vers x∗ . De manière
plus précise, on a :

Théorème 4. soit l’intervalle [a, b] tel que :


 f (a)f (b) < 0

f 0 (x) 6= 0 ∀x ∈ [a, b]


f 00 (x) 6= 0 ∀x ∈ [a, b]

Alors f (x) = 0 possède une seule racine dans [a, b] et ∀x0 ∈ [a, b] la suite :

f (xn )
xn+1 = xn − (36)
f 0 (xn )
converge quadratiquement.

Avec les hypothèses ci-dessus sur la fonction f , la méthode de Newton Raphson est
un cas particulier de la méthode du point fixe en choisissant :

f (x)
F (x) = x −
f 0 (x)

Etant donnée que F 0 (x∗ ) = 0 , d’aprés (35), la convergence de la méthode de Newton

Raphson est quadratique.

L’inconvénient de la méthode réside dans le fait qu’il faut savoir calculer la valeur de
la dérivée aux points xn , ce qui n’est pas toujours possible ou facile ; si l’on utilise une
valeur approchée de la dérivée -ce que font certains logiciels- la convergence n’est pas
quadratique en général.

Dans la pratique, il est rare de savoir si l’on se trouve dans le voisinage [x∗ − δ, x∗ + δ]
dont l’existence est assurée par le théorème précédent. La suite (xn ) peut alors diverger. . .
CHAPITRE 4. MÉTHODES ITÉRATIVES POUR LA RÉSOLUTION DE F(X)=0 32

Le choix du point de départ est crutial. Pour assurer la convergence on choisira un x0

telle que la condition de Fourier soit vérifiée à savoir :f 00 (x0 )f (x0 ) > 0

Pour les tests d’arrêt on peut utiliser les mêmes tests que ci-haut : (27), (28) ou (29).
Exemple. Comme application de Newton Raphson calculer la racine pième d’un nombre
a.

4.6 Exercices
4.6.1 Exercice sur MATLAB
On considère l’équation : 3 sin x − 1 = 0
1. Montrer qu’il existe une seule racine dans l’intervalle [0, π2 ] que l’on appelle x∗ .
2. Ecrire les algorithmes correspondant aux méthodes présentées ci-dessus et les tra-
duire en MATLAB .
3. Calculer une valeur approchée de x∗ à 10−8 près à l’aide de ces méthodes en com-
parant la rapidité de convergence.

4.6.2 Exercice
√ √
On considère le problème de calculer 2. Cela revient à trouver le zéro positif x∗ = 2
de la fonction
f (x) = x2 − 2
c’est-à-dire à résoudre une équation non linéaire.
1. On applique la méthode de dichotomie à la fonction f sur l’intervalle [a, b], avec
a = 1 et b = 2. Ecrire l’algorithme correspondant, et prouver que si x(m) est la
valeur approchée de x∗ trouvée à la m-ième itération de la méthode, on a
1
|x(m) − x∗ | ≤
2m+1
a+b
en commençant les itérations par x(0) = 2
.
√
2. Combien d’iterations sont nécessaires pour trouver une valeur approchée de 2 qui
soit exacte jusqu’au dixième chiffre aprés la virgule ?
3. Pour résoudre le problème, on peut utiliser aussi la méthode du point fixe

x(k+1) = F (x(k) )

Vérifier que x∗ est un point fixe de la fonction

1 1
F (x) = − x2 + x +
4 2
4. A l’aide du théorème des accroissements finis, prouver que dans l’intervalle [a, b] la
fonction F satisfait l’inégalité

|x(k+1) − x∗ | ≤ C|x(k) − α|

et donner une valeur pour la constante C.

CHAPITRE 4. MÉTHODES ITÉRATIVES POUR LA RÉSOLUTION DE F(X)=0 33

5. A partir de l’inégalité ci-dessus, vérifier que

|x(k) − x∗ | ≤ C k |x(0) − x∗ |, ∀k ≥ 0

Quel est le comportement de la suite {x(k) } lorsque k → ∞ ? Combien d’itérations

de
√ la méthode de point fixe sont nécessaires pour trouver une valeur approchée de
2 qui soit exacte jusqu’au dixième chiffre aprés la virgule ?

4.6.3 Exercice sur MATLAB

On considère la fonction √
f (x) = ex + 3 x − 2
dans l’intervalle [0, 1].
1. Tracer la courbe de y = f (x) et relever une valeur approchée de son zéro x∗ .
2. On veut calculer le zéro x∗ de la fonction f par la méthode de point fixe convenable.
En particulier on se donne deux méthodes de point fixe x = ϕi (x), i = 1, 2, où les
fonctions ϕ1 et ϕ2 sont définies respectivement comme
√ (2 − ex )2
ϕ1 (x) = ln(2 − 3 x) et ϕ2 (x) =
9
Pourquoi doit-on choisir la fonction ϕ2 plutôt que ϕ1 ? Justifier votre réponse.
3. Calculer la valeur de x∗ en partant de x0 = 0. Afficher le nombre de boucles néces-
saires pour une tolérance |x(k + 1) − x(k)| ≤ 10−10
Chapitre 5

Méthodes itératives pour la résolution

des systèmes de CRAMER

Les méthodes présentées ici sont la généralisation au cas n-dimensionnel des méthodes
de résolution de f (x) = 0 étudiées dans le cas mono-dimensionnel du chapitre précédent.
A partir d’une valeur estimée initiale : →
− −
→(∗) =
(0) (0)
x (0) = (x1 , x2 , . . . , x(0) T
n ) de la solution exacte : x
(∗) (∗) −
→ −
→
(x1 , x2 , . . . , x(∗) T
n ) du système A x = b on génère une suite de vecteurs :

−
→ −
→(k) →(k−1)
x (k+1) = f (− x ,...,−
→
x (k−i) ) (37)
−
→(k)
Définition 11. Si la fonction vectorielle d’itération f est indépendante de l’itération
k, on dit que l’on a une itération stationnaire, dans le cas contraire elle est dite non
stationnaire.
Définition 12. Si le calcul de −
→
x (k) demande la connaissance de i vecteurs estimés qui le
précèdent , on a une formule itérative à i pas ou multipas.

Dans la suite, on ne présentera que des méthodes itératives stationnaires à un pas.

5.1 Principe
Ecrivons la matrice A sous forme A = M − N (M doit être régulière) alors le système
(37) peut s’écrire :
−
→
M− →x = N− →
x + b
La méthode itérative consiste , à partir d’un vecteur initial →
−x (0) , à générer la suite
→
−x (1) , →
−
x (2) , . . . , −
→
x (k+1) de la manière suivante :



−
→x (1) = M −1 N −
→
x (0) + M −1 b


 −
→x (2) = M −1 N −
→
x (1) + M −1 b
 ...



 −
→(k+1) −1 −→(k) −1
x =M Nx +M b
Soit, sous forme condensée :

−
→(k+1) = T −
→x (k) + −
→
 x
 v k = 0, 1, 2, 3, ... (38)
−1
T =M N

 − −
→
→
v = M −1 b
La matrice d’itération T et le vecteur −
→
v sont indépendants de k.

34
CHAPITRE 5. MÉTHODES ITÉRATIVES POUR LA RÉSOLUTION DES SYSTÈMES DE CRAMER

5.2 Convergence
Théorème 5. La suite définie par → −x (0) et −
→
x (k+1) = T −
→x (k) + →
−
v converge vers
→
−∗ −1 −
→ −
→
x = (I − T ) v quelque soit x (0)
si et seulement si ρ(T ) < 1.

La démonstration de ce théorème se fait en deux temps :

– montrons d’abord que −→
x (k) − −
→
x ∗ = T k (−
→
x (0) − −
→
x ∗ ). En effet, nous avons :



−
→
x (k) − →
−
x ∗ = T (−
→
x (k−1) − −
→
x ∗)


 −
→(k−1) − →∗ →
− (k−2) → −∗
x − x = T( x − x )

 ...


 −
→(k)
x −→
−
x ∗ = T k (→
−
x (0) − −
→
x ∗)

La convergence de ce processus est assuré si, quelque soit −

→x (0) , on a limk→∞ T k = 0.
Où 0 désigne ici la matrice nulle.
– les normes matricielles (19),(20) et (21) sont équivalentes. Il est aisé de montrer que
|kAk| ≥ ρ(A). En effet, si λ est valeur propre de A associée au vecteur propre − →
x,
alors :
|λ|k−→x k = kλ→−
x k = kA− →x k ≤ k−
→x k |kAk| ⇒ |kAk| ≥| λ |

Comme on a |kT k| ≥ ρ(T ), une condition suffisante de convergence est telle que
|kAk| < 1 =⇒ ρ(T ) < 1

5.3 Décomposition de A
La matrice M sera choisie de façon qu’elle soit facilement inversible et vérifie ρ(M −1 N ) <
1 (ou ce qui est suffisant |kM −1 N k| < 1).
Définissons les matrices A = D − L − U telles que :

dij = δij aij

(
−aij si i > j
lij = (39)
0 sinon
(
−aij si i < j
uij =
0 sinon
Le chois de M correspond aux méthodes suivantes :

Méthode Matrice M Matrice N Matrice T

Jacobi M =D N =L+U TJ = D−1 (L + U )
Gauss-Seidel M =D−L N =U TGS = (D − L)−1 U
Relaxation M=D ω
−L 1
N = ω−1 D+U Tω = M −1 N

5.4 Méthode de Jacobi

−
→ →
−
x (k+1) = D−1 (L + U )−
→
x (k) + D−1 b
CHAPITRE 5. MÉTHODES ITÉRATIVES POUR LA RÉSOLUTION DES SYSTÈMES DE CRAMER

La composante i du vecteur −
→
x (k+1) s’écrit :
X n
(k+1) 1 (k)
xi = (bi − aij xj ) (40)
aii j6=i

Tests d’arrêt
Pour arrêter les itérations on peut utiliser l’un des tests suivants :
1. Lorsqu’on connait l’ordre de grandeur de la solution, on peut arrêter les itérations
dés que :
k−
→
x (k+1) − −
→
x (k) k < ²
Attention ! Ce test n’implique pas que k→ −
x (k) − −
→x ∗ k < ².
2. L’erreur relative est inférieure à une précision donnée :
n
ke(k) k X (k)
≤² avec eki = bi − aij xj
kbk j6=i

3. Un autre critère courant concerne l’amélioration relative sur −

→
x :
k−
→
x (k+1) − −→x (k) k
<²
k−
→
x (k) k

Convergence
Pour la convergence la condition suffisante |kTJ k| < 1 se traduit par :
n
X
|aij | < |aii | i = 1, 2, 3, . . .
j6=i

La matrice A doit être à diagonale dominante.

5.5 Méthode de Gauss-Seidel

−
→ −
→
x (k+1) = D−1 L−→x (k+1) + D−1 U −
→
x (k) + D−1 b
La composante i du vecteur −
→x (k+1) s’écrit :
j=i−1
X Xn
(k+1) 1 (k+1) (k)
xi = (bi − aij xj − aij xj ) (41)
aii j=1 j=i+1

La convergence est généralement deux fois plus rapide que celle de la méthode de Ja-
cobi (mais il y a de nombreux contres-exemples). On a encore une condition suffisante de
convergence : on montre de la même manière qu’il suffit que la matrice A soit à diagonale
dominante.

Remarque
1. On remarque, d’une part, qu’il suffit de stocker en mémoire un seul vecteur que
l’on «modifie sur place » ; on peut alors omettre les indices supérieures et écrire
l’algorithme si contre :
CHAPITRE 5. MÉTHODES ITÉRATIVES POUR LA RÉSOLUTION DES SYSTÈMES DE CRAMER

−
→x est choisie arbitrairemant
−
→ → −
while k b −A
→
−
xk <²
k b k
for i = 1 : n
P
xi = xi + a1ii (bi − nj=1 aij xj ) ;
end
end
2. On voit ainsi clairement apparaitre le terme correctif appliqué à xi ; en augmentant
l’importance de cette correction on obitent la méthode de relaxation.

5.6 Méthode de relaxation

Soit ω ∈]0, 2[ un facteur de relaxation ; l’algorithme s’écrit :
→
−x est choisie arbitrairemant
−
→ − →
while k b −A
−
→
xk <²
k b k
for i = 1 : n
P
xi = xi + aωii (bi − nj=1 aij xj ) ;
end
end
qui se traduit matriciellement :
( (
x(0) arbitraire T = (D − ωL)−1 ((1 − ω)D + ωU )
où −
→
pour k ≥ 0, −
→
x (k+1) = T −
→
x (k) + −
→
v −
→
v = (D − ωL)−1 ω b
Reamrque : pour ω = 1 on retrouve la méthode de Gauss-Seidel.