Département d'Ingénierie

Cours d'Analyse Numériques

MASTER 1

Dr. DOSSO Mouhamadou

Maître de Conférences
E-mail : [Link]@[Link]
mouhamadoudoss@[Link]
Table des matières

1 Quelques Notions fondammentales 1

1.1 Norme de vecteurs et de matrices . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Conditionnement d'une matrice . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Préconditionnement d'un système linéaire . . . . . . . . . . . . . . . . . . 4
1.3.1 Préconditionnement à gauche . . . . . . . . . . . . . . . . . . . . . 4
1.3.2 Préconditionnement à droite . . . . . . . . . . . . . . . . . . . . . . 4
1.3.3 Préconditionnement bilatéral (ou à gauche et droite) . . . . . . . . 4
2 Polynômes orthogonaux 7
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Dénition, existence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3 Relation avec les polynômes habituels . . . . . . . . . . . . . . . . . . . . 9
2.4 Relation de Récurrence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.5 Equation Dierentielle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.6 Fonction génératrices, Formule de Rodrigués et Identité de Darboux Chris-
tofel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.6.1 Fonction génératrice . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.6.2 Formule de Rodrigués . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.6.3 Identité de Daboux-Christofel . . . . . . . . . . . . . . . . . . . . . 13
2.7 Polynômes particuliers connus . . . . . . . . . . . . . . . . . . . . . . . . 13
2.7.1 Polynômes de Legendre . . . . . . . . . . . . . . . . . . . . . . . . 13
2.7.2 Polynômes d'Hermite . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.7.3 Polynômes de LaGeurre . . . . . . . . . . . . . . . . . . . . . . . . 14
2.7.4 Polynômes de Tschebychef . . . . . . . . . . . . . . . . . . . . . . . 15
2.8 Autres polynômes classiques . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.8.1 Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.8.2 LaGeurre Généralisé . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.9 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3 Coubes de Bézier et Polynômes de Bernstein 19
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2 Quelques notions de base sur les courbes de Bézier : Rappels sur la notion
de barycentres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3 Polynôme de Bernstein . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.4 Courbe de Bézier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
 TABLE DES MATIÈRES

3.4.1 Dénitions et propriétés des courbes de Bézier . . . . . . . . . . . . 22

[Link] Dérivée d'une courbe de Bézier . . . . . . . . . . . . . . 23
[Link] Tangentes aux extrémités . . . . . . . . . . . . . . . . . . 23
3.4.2 Construction par barycentres successifs . . . . . . . . . . . . . . . 24
3.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4 Méthodes itératives de base 29
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2 Convergence de matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.3 Méthodes de relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.3.1 Généralités sur la méthode . . . . . . . . . . . . . . . . . . . . . . 33
4.3.2 Cas particuliers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.3.3 Tests d'arrêt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.3.4 Les méthodes de Jacobi et de Gauss-Seidel . . . . . . . . . . . . . 36
[Link] La méthode de Jacobi . . . . . . . . . . . . . . . . . . . . 37
[Link] La méthode de Gauss-Seidel . . . . . . . . . . . . . . . . . 38
4.3.5 Convergence des méthodes de Jacobi et de Gauss-Seidel . . . . . . 38
[Link] Résultats préliminaires . . . . . . . . . . . . . . . . . . . . 39
[Link] Convergence des méthodes . . . . . . . . . . . . . . . . . . 39
4.3.6 Resultats particuliers de convergence . . . . . . . . . . . . . . . . . 40
4.4 La méthode de sur-relaxation . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.4.1 Convergence de la méthode de sur-relaxation . . . . . . . . . . . . 42
4.4.2 La méthode de sur-relaxation symétrique . . . . . . . . . . . . . . 43
4.5 Les méthodes de Richardson . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.6 Pseudo-codes et codes Matlab . . . . . . . . . . . . . . . . . . . . . . . . 44
4.6.1 La méthode de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.6.2 La méthode de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . 45
4.6.3 La méthode de sur-relaxation . . . . . . . . . . . . . . . . . . . . . 45
4.7 Exercice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5 Calcul de valeurs propres 47
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.2 La méthode de la puissance . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.2.1 Accélération de la convergence . . . . . . . . . . . . . . . . . . . . . 51
5.2.2 Variantes de la méthode de la puissance . . . . . . . . . . . . . . . 52
[Link].1 La méthode des traces . . . . . . . . . . . . . . . 52
[Link].2 La méthode du produit scalaire . . . . . . . . . . 53
5.2.3 Recherche des autres valeurs propres . . . . . . . . . . . . . . . . . 54
5.2.4 La méthodes de la puissance inverse . . . . . . . . . . . . . . . . . 56
5.3 Méthodes de décomposition . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.3.1 La méthode de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3.2 La méthode Greenstadt . . . . . . . . . . . . . . . . . . . . . . . . 60
5.3.3 L'algoritme LR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.3.4 l'Agorithme QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.4 Pseudo-codes et Codes Matlab . . . . . . . . . . . . . . . . . . . . . . . . 65
 iii

5.4.1 L'algorithme de Jacobi, choix cyclique . . . . . . . . . . . . . . . . 65

5.4.2 L'algorithme de Jacobi, choix classique . . . . . . . . . . . . . . . . 65
5.4.3 L'algorithme QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
TABLE DES MATIÈRES
Chapitre 1
Quelques Notions fondammentales

Sommaire
1.1 Norme de vecteurs et de matrices . . . . . . . . . . . . . . . . 1
1.2 Conditionnement d'une matrice . . . . . . . . . . . . . . . . . . 3
1.3 Préconditionnement d'un système linéaire . . . . . . . . . . . 4
1.3.1 Préconditionnement à gauche . . . . . . . . . . . . . . . . . . . 4
1.3.2 Préconditionnement à droite . . . . . . . . . . . . . . . . . . . . 4
1.3.3 Préconditionnement bilatéral (ou à gauche et droite) . . . . . . 4

1.1 Norme de vecteurs et de matrices

Soit E une espace vectoriel sur C. On appelle norme toute application de E dans R
qui, à tout x ∈ E, associe le nombre réel ∥x∥ (appelé norme de x) qui vérie les trois
conditions suivantes :
1. ∥x∥ ≥ 0 et ∥x∥ = 0 si et seulement si x = 0 ∈ E.
2. ∥λx∥ = |λ|.∥x∥, ∀λ ∈ C,
3. ∥x + y∥ ≤ ∥x∥ + ∥y∥, ∀x, y ∈ E.
Comme c'est le cas pour la valeur absolue, une norme sert à mésurer la distance ∥x−y∥
entre deux éléments. Sur un espace vectoriel, il se peut que l'on puisse dénir plusieurs
normes. cependant, si l'espace est de dimension nie, toutes les normes sont équivalentes
c'est-à-dire qu'il existera des inégalités entre elles. Par conséquent si deux éléments sont
voisins, au sens d'une certaine norme (et en particulier si une suite converge pour une
certaine norme), il en sera de même pourune autre norme.
x1
Considérons le cas E = C . Soit
n
 x2 
 ...  ∈ E
x=  un vecteur. Les quantités suivantes
xn
 CHAPITRE 1. QUELQUES NOTIONS FONDAMMENTALES

s'appellent normes de Hölder d'indice k

1/k
∥x∥k = |x1 |k + · · · + |xn |k , k = 1, 2, ...,

∥x∥∞ = max |xi |.

1≤i≤n

Parmi ces normes, les plus utilisées sont ∥x∥ , ∥x∥ et ∥x∥ . On les désigne souvent
respectivement sous les termes de normes l , l et l . Cette dernière représente la longueur
1 ∞ 2

du vecteur x au sens euclidien du terme; elle s'appelle norme euclidienne.

1 ∞ 2

Une norme de matrice peut être dénie à partir d'une norme pour les vecteurs (mais
cela n'est pas obligatoire). Soit A ∈ C une matrce. La quantité
n×m

∥Ax∥
∥A∥ = sup
x̸=0 ∥x∥

est une norme pour la matrice A. Puisque l'on peut, dans cette dénition, remplacer x
par αx, où α est un scalaire, on peut toujours choisir x de norme 1 et l'on a également
∥A∥ = sup ∥Ax∥.
∥x∥=1

Ces normes de matrices vérient les trois (3) propriétés des normes mais, en plus, il existe
deux autres propiétés qui serons très utiles par la suite
∥Ax∥ ≤ ∥A∥.∥x∥

et
∥AB∥ ≤ ∥A∥.∥B∥.
On appelle multiplicative toute norme vériant cette dernière inégalité.
Les normes de matrices les plus utilisées sont celles qui sont reliées à une norme de Hölder
pour les vecteurs c'est-à-dire
∥Ax∥k
∥A∥k = sup
x̸=0 ∥x∥k
où ∥.∥ est une norme de vecteur de Hölder d'indice k.
Les normes de matrice semblent diciles à calculer en pratique puisqu'elle font intervenir
k

une borne supérieure. Cependant, on connait leurs expressions dans trois cas pour les
normes de Hölders :
la norme une
X n
∥A∥ = max
1 |a |, ij
1≤j≤n
i=1

la norme innie
Xn
∥A∥∞ = max |aij |,
1≤i≤n
j=1
p
∥A∥2 = ρ(AT A), la norme 2 ou la norme spectrale
 3

où ρ(A A) désigne le rayon spectral de la matrice A A ∈ C , c'est-à-dire sa plus

T T m×m

grande valeurs propres. A A étant symétrique dénie positive, ses valeurs propres sont
T

réelles et positives. cette dernière norme s'appelle aussi norme spectrale.

Maintenant, on considère que les matrices carrées (m = n).
On montre que quelques soit la norme,
ρ(A) ≤ ∥A∥.
D'autre part, si A est symétrique, ∥A∥ = ρ(A).
On utilise aussi parfois la norme de Frobenius car elle est facile à calculer. Elle est dénie
2

par ! 1/2
n
X
∥A∥F = |aij |2 .
i,j=1

Pour cette norme, on a également

∥AB∥F ≤ ∥A∥F .∥B∥F .
De plus 2 T
∥A∥ = tr(A A).
où tr désigne la trace d'une matrice c'est-à-dire la somme des éléments de sa diagonale.
F

1.2 Conditionnement d'une matrice

Soit A ∈ C . Pour toute norme de Hölder, on a
n×n

∥AA−1 ∥ = ∥I∥ ≤ ∥A∥.∥A−1 ∥.

Or, la norme de la matrice identité est égale à 1 d'après la dénition. Il s'en suit que l'on
a l'inégalité −1
κ(A) = ∥A∥.∥A ∥ ≥ 1.
Le nombre κ(A) s'appelle le conditionnement de A.
On a les propriétés suivantes
1. κ(λA) = κ(A), ∀λ ̸= 0,
2. κ(A ) = κ(A) puisque A et A jouent les rôles symétriques dans la dénition de
−1 −1

κ(A),

3. 0 < ∥A1 ∥ ≤ ∥Ax∥

−1 ∥x∥
≤ ∥A∥, ∀x.
La notion de conditionnement d'une matrice est absolument fondamentale pour les mé-
thodes de résolution des systèmes d'équations linéaires. On dit que la matrice A est bien
conditionnée si κ(A) est voisin de 1. Cependant, si κ(A) est grand par rapport à 1,
on dit que A est mal conditionnée. Naturellement les adjectif voisin et grand sont
subjectifs. Si la précision de l'ordinateur avec lequel on travaille est de 10 , un condi-
−7

tionnement de 10 sera considéré comme grand. Par contre, si la précision est de 10 ,

5 −16

un tel conditionnement sera petit.

 CHAPITRE 1. QUELQUES NOTIONS FONDAMMENTALES

1.3 Préconditionnement d'un système linéaire

Soit à résoudre le système d'équations linéaires
Ax = b. (1.1)
Si A est mal conditionnée, une petite perturbation des données du système linéaire (c'est-
à-dire sur la matrice A et/ou le second membre b) pourra engendrer une grande variation
dans la solution exacte. De plus, dans ce cas, les erreurs provenant de l'arithmétique à
précision nie des ordinateurs pourront aecter les résultats d'une erreur importante. Si
le conditionnement est mauvais, un certain nombre de méthodes itératives convergeront
lentement car plus le conditionnement est grand et plus leur vitesse de convergence est
faible. C'est le cas, par exemple :
• de la méthode de la plus profonde descente
• et de la méthode du grandient conjugué.
Des méthodes convergentes pourront même ne plus converger du tout. Un manière de
remédier à ce problème est de préconditionner le système.
Le préconditionnement consiste à remplacer le système Ax = b par
CAx = Cb.
La matrice C est choisie de sorte que le conditionnement de CA soit plus petit que celui
de A.
1.3.1 Préconditionnement à gauche
On dit que M est une matrice de préconditionnement à gauche si le système (1.1)
original est transformé sous la forme
M −1 Ax = M −1 b.

1.3.2 Préconditionnement à droite

De manière analogue, on dira que M est une matrice de préconditionnement à droite
si le système (1.1) est transformé sous la forme
AM −1 y = b.
Dans ce cas, la solution s'obtient en posant :
x = M −1 y.

1.3.3 Préconditionnement bilatéral (ou à gauche et droite)

En principe, on peut combiner ces deux type de préconditionnement. Ils s'agira d'avoir
deux matrices inversibles M et M et d'écrire le système (1.1) sous la forme
G D

MG−1 AMD−1 y = MG−1 b.

 5

En générale, on utilise un préconditionnement à gauche. Mais il y a un cas très important

où ce type de préconditionnement est nécessaire.
CHAPITRE 1. QUELQUES NOTIONS FONDAMMENTALES
Chapitre 2
Polynômes orthogonaux

Sommaire
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . ........ 7
2.2 Dénition, existence . . . . . . . . . . . . . . . . ........ 8
2.3 Relation avec les polynômes habituels . . . . . . ........ 9
2.4 Relation de Récurrence . . . . . . . . . . . . . . . . . . . . . . 10
2.5 Equation Dierentielle . . . . . . . . . . . . . . . . . . . . . . . 11
2.6 Fonction génératrices, Formule de Rodrigués et Identité de
Darboux Christofel . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.6.1 Fonction génératrice . . . . . . . . . . . . . . . . . . . . . . . . 12
2.6.2 Formule de Rodrigués . . . . . . . . . . . . . . . . . . . . . . . 12
2.6.3 Identité de Daboux-Christofel . . . . . . . . . . . . . . . . . . 13
2.7 Polynômes particuliers connus . . . . . . . . . . . . . . . . . . 13
2.7.1 Polynômes de Legendre . . . . . . . . . . . . . . . . . . . . . . 13
2.7.2 Polynômes d'Hermite . . . . . . . . . . . . . . . . . . . . . . . 14
2.7.3 Polynômes de LaGeurre . . . . . . . . . . . . . . . . . . . . . . 14
2.7.4 Polynômes de Tschebychef . . . . . . . . . . . . . . . . . . . . 15
2.8 Autres polynômes classiques . . . . . . . . . . . . . . . . . . . 16
2.8.1 Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.8.2 LaGeurre Généralisé . . . . . . . . . . . . . . . . . . . . . . . . 16
2.9 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.1 Introduction
Les polynômes orthogonaux apparaissent très couramment en physique mathéma-
tique; en particulier, dans la résolution d'équation aux dérivées partielles (Laplace, Schrö-
dinger) par la méthode de séparation de variables. Certaines d'entre elles sont aussi très
utilisées en analyse numérique. Cependant, l'orthogonalité impose que ces polynömes or-
thogonaux aient en commun un certain nombre de propriétés, en particulier celles de
 CHAPITRE 2. POLYNÔMES ORTHOGONAUX

vérier une relation de récurrence à trois termes et d'obéir à une équation diérentielle
linéaire du second ordre.
2.2 Dénition, existence
Soient [a, b] un intervalle (appelé I dans la suite qui peut être ni ou inni) et w une
fonction strictement positive et intégrable sur cet intervalle.
on appelle polynômes ortogonaux sur I , par rapport à la "fonction de poids" w, une
suite de polynômes G (x), G (x), ...,G (x), ...(où G est le dégré k) tel que :
0 1 n k

w(x)G (x)G (x)dx = 0, si k ̸= n. (2.1)

Z b
(G , G ) =
k n k n
a

Le nombre
(2.2)
Z
(f, g) = f (x)g(x)w(x)dx

est souvent appelé "produit scalaire" de f par g. Lorsqu'il est nul, on dit que les
I

fonctions f et g sont "orthogonaux".

Remarque 2.1
 on a
(2.3)
Z
2
∥f ∥ = (f, f ) = (f (f (x))2 w(x)dx

 (Inégalité de Cauchy-Schwartz) :
I

|(f, g)| ≤ ∥f ∥.∥g∥(2.4)

 Les équations (2.1) et (2.2) ne dénissent pas entièrement les G ; il est donc
necessaire d'imposer les condions supplémentaires :
n

Z b Z b
2
(Gn , Gn ) = w(x) (Gn (x)) dx = 1 ; xw(x)Gn (x)Gn+1 (x)dx > 0
a a
(2.5)
On a alors
(2.6)
Z b
(Gk , Gn ) = w(x)Gk (x)Gn (x)dx = δk,n .

Les polynômes sont ainsi orthonormés.

a

Pour prouver l'existence de ces polynômes, nous allons les construire en utilisant le procedé
d'orthogonalisation de Gram-Schmidt, appliqué à l'ensemble des puissances successives de
x.
◦ G est une constante positive, telle que
0

(G0 , G0 ) = 1,
 9

g1 est une fonction linéaire défnie comme :

g1 = x − (x, G0 )G0

Ainsi, g est manifestement orthogonale à G .

G est dénie comme g normé
1 0
1 1

G1 = g1 /∥g1 ∥.

◦ g2 est une fonction quadratique dénie par

g2 = x2 − (x2 , G0 )G0 − (x2 , G1 )G1 ,

et
G2 = g2 /∥g2 ∥.
En formant le produit scalaire de g successivement par G et G , nous vérions
que cette fonstions est bien orthogonale à G et G .
2 0 1

◦ Nous écrivons à l'ordre n

0 1

gn = xn + an,n−1 Gn−1 + an,n−2 Gn−2 + · · · + an,0 G0

Xn
n
=x + an,n−k Gn−k
k=1

et nous choisissons les constantes a de telles façons que g soit orthogonal aux
G pour j = 0, ..., n − 1. Ceci impose que
n,j n
j

an,j = −(xn , Gj ).

Il reste à normaliser ce polynôme :

gn
Gn = .
∥gn ∥

Nous construisons ainsi la suite des polynômes G , à partir des nombres (x, G ), (x , G ),...,(x , G ),...
k 0
2
0
k
k

2.3 Relation avec les polynômes habituels

Théorème 2.1 Tout polynôme de dégré n est représentable comme combinaison linéaire
de polynômes G , avec k ≤ n.
k

Preuve
Nous avons construit les G comme combinaisons linéaires des x , p ≤ k ; de plus, le
k

coecients de x dans G est non nul. Il en résulte que ces relations linéaires peuvent
k
k

s'inverser et que x est une fonction linéaire des G , k ≤ p. Un polynôme quelconque Φ,

k
k

lui-même combinaison linéaire des x , peu donc s'explimer comme :

k
p
 CHAPITRE 2. POLYNÔMES ORTHOGONAUX

Φ = a0 G0 + a1 G1 + · · · + an Gn .
La détermination des a est facile :
eectuons, en eet, le produit scalaire de l'égalité précédente par G .
i
(p ≤ n)
On a :
k

ap = (Φ, Gp ) , p = 0, 1, 2, ...., n.
Nous avons au contraire (Φ, G ) = 0, quelque soit Φ(x) de dégré n, puisque G est
orthogonal à tous G d'indice inférieur. □
n+1 n+1

Nous en déduisons le théorème suivant

k

Théorème 2.2 Le polynôme G est orthogonal à tout polynôme de degré inférieur.

k

Par ailleur, Φ(x) sera identiquement nul si et seulement si :

(Φ, Gp ) = aq = 0, p = 0, 1, 2, ..., n.

Remarque 2.2  Le produit xcalaire (xG , G ) est le produit d'un polynôme de

degré n, xG , par G ; si b dénit le coecient de G dans le developpement de
n−1 n

xG , le produit scalaire considéré se reduit à :

n−1 n n n
n−1

(xGn−1 , Gn ) = bn (Gn , Gn ) ̸= 0.

 Les zeros de G appartiennent à I et sont tous simples.

n

Les zeros des polynômes successifs adoptent une disposition particulière, résumée par le
théorème ci-dessous
Théorème 2.3 Les zéros de G séparent ceux de G .
k k+1

Preuve
En eet, ce résultat est vrai pour G et G ; il se démontre par récurrence dans le cas
général. □
0 1

On peut dire aussi que la suite des G forme une suite de Sturn.
n

2.4 Relation de Récurrence

Nous allons vérier l'existence d'une relation de résurrence, de la fonction :
Gn+1 (x) = (αn x + βn ) Gn (x) + γn Gn−1 (x) (2.7)
où les α , β , γ sont des constantes.
Cette propriété est une conséquence de l'orthogonalité. Soit c le coecient de x dans
n n n
k

G (terme de plus haut degré).

k

Posons
k

ak = ck+1 /ck .
 11

Considérons le polynôme
F = Gn+1 − an xGn .
Il est de degré n au plus (à cause de la dénition de a ) et peut donc s'exprimer comme
une combinaison linéaire des G , i ≤ n.
n
i

n
X
F = bj Gj .
j=0

Considérons
(F, Gp ) = (Gn+1 , Gp ) − an (xGn , Gp ) = −an (xGn , Gp ),
à conditions que p < n + 1.
Or, l'examen de l'intégrale qui dénit (xG , G ) montre que : (xG , G ) = (G , xG ).
Comme xG est de degré p + 1 au plus, il est orthogonal à G tant que p < n − 1.
n p n p n p

On vient de montrer que F était orthogonal à G , 0 ≤ p ≤ n − 2.

p n

Le développement de F s'écrit alors :

p

F = Gn+1 − an xGn = bn Gn + bn−1 Gn−1

si bien que
Gn+1 = (an x + bn )Gn + bn−1 Gn−1 . (2.8)
Cette relation est bien de la forme annoncée à condition de choisir α ,
= an β n = b n
et γ = b .
n

En posant G = 0, on dénit complétement les G .

n n−1
−1 n

2.5 Equation Dierentielle

On peut montrer que chaque G obéit à une équation diérentielle du second ordre.
Il est plus facile de décomposer une réciproque partielle. Soit l'équation diérentielle :
k

A(x)y” + B(x)y ′ + Cn y = 0

où C est une constante (avec C ̸= C si n ̸= p), A(x) et B(x) deux fonctions régulières
de x. Supposons que l'équation admet une solution polynômiale y = G (x) de g de degré
n n p

n pour chaque valeur de n.

n

On peut alors trouver un intervalle [a, b] et une fonction de poids w tels que la suite des
G soit orthogonale par rapport à ces éléments. En eet, G et G satisfont séparément
à l'équation proposée :
n n p

′′
AG + BG + C G = 0 ′
(2.9)
(2.10)
n n n n

′′
AG + BG + C G = 0.
p
′
p p p (2.11)
 CHAPITRE 2. POLYNÔMES ORTHOGONAUX

En multipliant, la premier relation par wG et la seconde par wG et retranchant membre

à membre, il vient
p n

Aw(Gn Gp ” − Gp Gn ”) + Bw(Gn G′p − Gp G′n ) + w(Cp − Cn )Gn Gp = 0

équation qui s'écrit encore


Aw(Gn G′p − Gp G′n ) + [wB − (wA)′ ] (Gn G′p − Gp G′n ) + w(Cp − Cn )Gn Gp = 0

(2.12)
Considérons l'équation diérenteille
uB − (uA)′ = 0

et notons w la solution qui obéit aux conditions aux limites

Aw |a = Aw |b = 0.

En intégrant terme à terme l'équation (2.12) avec cette dénition de w, nous obtenons
Z b
(Cp − Cn ) wGn Gp dx = 0.
a

Ce qui signie que les polynômes G sont orthogonaux sur l'intervalle I par rapport à la
fonction de poids w.
n

2.6 Fonction génératrices, Formule de Rodrigués et Iden-

tité de Darboux Christofel
2.6.1 Fonction génératrice
On peut en général trouver une fonction g(u, x) telle que l'on puisse écrire
(2.13)
∞
X
g(u, x) = Ck Gk (x)uk
k=0

les G étant des constantes et u une variable réelle auxiliaire. Si la fonction g a une forme
analytique "simple" par rapport à u et à x, on l'appelle génératrice de la suite des G .
k
k

N.B.: 2.1 La fonction gégératrice peut servir à dénir les G ; les propriétés de ces poly-
nômes se démontrent alors par manipulations de g.
k

2.6.2 Formule de Rodrigués

Les G s'expriment en fonction de la dérivée d'ordre n d'une fonction U (x) (forme
de Rodriguées)
n n

(2.14)
n
1 d
G (x) =
n U (x). n n
w(x) dx
 13

La quelle est solution d'un problème diérentielle

dn+1 1 dn Un
 
=0
dxn+1 w(x) dxn
avec les conditions aux limites
Un = Un′ = · · · = Un(n−1) = 0 en x=a et x = b.

2.6.3 Identité de Daboux-Christofel

Une démonstraction laborieuse permet d'établir l'identité de Daboux-Christofel :
(2.15)
n
X 1
G (x)G (y) =k k [G (x)G (y) − G (x)G (y)]
n+1 n n n+1
k=0
(x − y)a n

souvent utilisée pour simplier des expressions impliquant les G (la constante a a deja
été dénie ci-dessus).
n n

2.7 Polynômes particuliers connus

2.7.1 Polynômes de Legendre
Les polynômes de Legendre sont orthogonaux sur [−1, 1] par rapport à la fonction de
poids w = 1.
A partir de P = 1, on peut orthogonaliser les puissances successives de x pour ontenir
0

P0 = 1 ; P1 (x) = x
P2 (x) = (3x2 − 1)/2 ;
P3 (x) = (5x3 − 3x)/2,
P4 (x) = (35x4 − 30x2 + 3)/8 ;
P5 (x) = (63x5 − 70x3 + 15x)/8.
Le polynôme P a le degré de parité de k. En intégrant par partie dans la relation d'or-
thogonalité (P , P ) = 0, on démontre une relation d'Olinde Rodrigués :
k
k l

1 dn 2
Pn (x) = (x − 1)
2n n! dxn
avec P (1) = 1.
À partir de ces dénitions, les P ne sont pas normalisés à 1 mais :
n
n

(Pn , Pn ) = 2/(2n + 1).

Il admettent les fonctions génératrices
+∞
1 X
g ′ (u, x) = √ = Pk (x)uk
1 − 2ux + u2
k=0
 CHAPITRE 2. POLYNÔMES ORTHOGONAUX

d'où l'on tire la relation de recurrence

(n + 1)Pn+1 − (2n + 1)xPn + nPn−1 = 0

Les polynômes de Legendre obéissent à l'équation diérenteille

(x2 − 1)Pn ” + 2nPn′ − n(n + 1)Pn = 0.

2.7.2 Polynômes d'Hermite

Les polynômes d'Hermide sont orthogonaux sur l'intervalle : [−∞, +∞], par rapport
à la fonction de poids : w(x) = exp(−x ) Ils satisfont à la relation de récurrence
2

Hn+1 − 2xHn + 2nHn−1 = 0

et a l'équation diérentielle
Hn ” − 2xHn′ + 2nHn = 0.

Ils peuvent se déduire de la fonction génératrice :

+∞
X
g(u, x) = exp(2ux − u2 ) = Hk (x)uk /k!
k=0

On connait aussi une fonction de Rodrigués

dn
Hn (x) = (−1)n exp(x2 ) exp(−x2 ).
dxn
Les polynômes d'Hermite admettent une représentation générale assez simple
n/2
X 1 1 1
Hn (x) = (− )p ( ) (2x)n−2p
p=0
2 p! (n − 2p)!

et les premiers représentations de l'espace s'écrivent

H0 (x) = 1 ; H1 (x) = 2x ; H2 (x) = 4x2 − 2 ;
H3 (x) = 8x3 − 12x ; H4 (x) = 16x4 − 48x2 + 12
H5 (x) = 32x5 − 160x3 + 120x.

On peut constater que H à la parité de k.

k

2.7.3 Polynômes de LaGeurre

Nous pouvons dénir ces polynômes à partir d'une formule d'Olinde Rodrigués
1 dn
Ln (x) = × ex n (e−x xn ),
n dx
 15

qui nous permet d'écrire les premiers membres de la suite :

1
L0 (x) = 1 ; L1 (x) = −x + 1 ; Lx (x) = (x2 − 4x + 2)
2
1
L3 (x) = (−x3 + 9x2 − 18x + 6).
6
Ils admettent la fonction génératrice
∞
1 xt X
g(x, t) = × exp(− )= Lk (x)tk
1−t 1−t k=0

et la relation de récurrence
(n + 1)Ln+1 + (x − 2n − 1)Ln + nLn−1 = 0

Ils satisfont à l'équation diérentielle

xLn ” + (1 − x)L′n + nLn = 0.

Les polynômes sont orthogonaux sur l'intervalle [0, ∞] par rapport à la fonction de poids
w(x) = e ; Ils sont normalisés à un : (L , L ) = 1.
−x
n n

2.7.4 Polynômes de Tschebychef

Ces polynômes admettent la dénition simple
Tn (x) = cos[n × arccos(x)]

On obtient une relation de récurrence entre T n−1 Tn et T en calculant cos[(n ± 1)x] :

n+1

Tn+1 (x) + Tn−1 (x) = 2xTn (x).

Nous en déduisons les expressions des premiers polynômes

T0 = 1 ; T1 (x) = x ; T2 (x) = 2x2 − 1 ; T3 (x) = 4x3 − 3x ;
T4 (x) = 8x4 − 8x2 + 1.

Les T obéissent à l'équation diérentielle

n

(1 − x2 )Tn ” − xTn′ + n2 Tn = 0.

Ils sont orthogonaux sur l'intervalle [−1, 1] par rapport à la fonction de poids
1
w(x) = √
1 − x2
Leur fonction génératrice s'écrit
+∞
1 − ux X
g(x, u) = 2
= uk Tk (x).
1 − 2ux + u k=0
 CHAPITRE 2. POLYNÔMES ORTHOGONAUX

2.8 Autres polynômes classiques

Voici les dénitions de deux suites polynômiales moins fréquemment rencontrées
2.8.1 Jacobi
I = [−1, 1] w(x) = (1 − x)α (1 − x)β , α, β > −1
Les polynômes de Jacobi recouvrent plusieurs cas particuliers intéressants :
 Polynômes de Legendre, lorsque α = β = 0
 Polynômes de Tschebychef de première espece lorsque α = β = −1/2 ;
 Polynômes de Tschebychef de deuxième espece lorsque α = β = 1/2 ;
 Polynômes de Tschebychef de troisième espece lorsque α = −β = −1/2 ;
 Polynômes de Tschebychef de quatrième espece lorsque α = −β = 1/2 ;
 Polynômes de Gegenbauer lorsque α = β = λ − 1/2 ;
2.8.2 LaGeurre Généralisé
;
I = [0, +∞] w(x) = xα e−x ; α > −1 ;
Les polynômes de Laguerre habituels sont obtenus en faisant α = 0 dans la dénition des
polynômes généralisés.
2.9 Exercices
Soit la liste des exercices du chapitre 1.
Exercice 2.1 Les polynômes de Tschebychef sont dénis pour x ∈ [−1, 1] par
Tn (x) = cos(n arccos x).

1. Former T , T , T et T comme fonctions de x ; il est cependant commode d<'introduire

la variable intermédiaire θ = arccos x ; préciser le domaine de dnition de cette
0 1 2 3

variable.
2. Démontrer la relation de récurrence
Tn+1 (x) = 2xTn (x) − Tn−1 (x).

Utiliser cette relation pour construire T et T . Quelle est la parité de T ?

3. Déterminer les zéros de T . Montrer qu'entre deux zéros successifs de ce polynôme,
4 5 n

on rencontre un zéros de T .
n

4. Majorer T sur [−1, 1].

n−1

5. Vérier la relation
n

Z 1
Tk (x)Tl (x)
√ dx = 0 k ̸= l.
1 − x2
Que vaut cette intégrale quand ?
−1

k=l
 17

6. Explimer x , x , x , x , x , et x en fonction des polynômes T , T , T , T , T et

0 1 2 3 4 5

T.
0 1 2 3 4

7. On note p (x) le polynôme de degré n qui représente le développement de Taylor

5

tronqué à l'ordre n de e au voisinage de zero et q le développement de e sur la

n
x x

base des T :
n
i

n n
x
X xk x
X
e ≃ ≡ pn (x) ; e ≃ ak Tk (x) ≡ qn (x).
k=0
k! k=0

Former p (x) et q (x).

8. Quelle erreur commet-on en remplaçant e par p ou par q sur le segment [−1, 1] ?
5 5
x

9.
5 5

Exercice 2.2 Déterminer, enutilisant le procédé d'orthogonalisation se Schmidt appliqué

aux puissances successives de x, les quatre premiers polynômes P qui vérient les relations
les relations :
k
Z 1
2δ k,l
P0 = 2; Pk (x)Pl (x)dx = .
2l + 1
Les polynômes ainsi obtenus sont appelés polynômes de Legendre.
−1

Exercice 2.3 hh
CHAPITRE 2. POLYNÔMES ORTHOGONAUX
Chapitre 3
Coubes de Bézier et Polynômes de
Bernstein

Sommaire
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2 Quelques notions de base sur les courbes de Bézier : Rappels
sur la notion de barycentres . . . . . . . . . . . . . . . . . . . . 20
3.3 Polynôme de Bernstein . . . . . . . . . . . . . . . . . . . . . . 20
3.4 Courbe de Bézier . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4.1 Dénitions et propriétés des courbes de Bézier . . . . . . . . . . 22
3.4.2 Construction par barycentres successifs . . . . . . . . . . . . . 24
3.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.1 Introduction
Dans les années 60, les ingénieurs Pierre BÉZIER et Paul DE CASTELIAU travaillant
respectivement chez Renault et Citröen, ont réfréchi au moyen de dénir de manière la
plus concise possible la forme d'une carroserie.
Le principe a été énoncé par BÉZIER mais l'algorithme de construction, lui, a été
énoncé par son collègue de la marque aux chevrons qui n'a d'ailleurs été dévoilé que
bien plus tard à cause de la loi du secret industriel qui a primé sur le développement
scientique.
Pierre BÉZIER (diplômé de l'ENSAM et de SUPELEC), à l'origine des prémières
machines à commandes numériques de la CAO à été mise à l'écart par sa direction. Il se
consacra alors presque exclusivement aux mathématiques et à la modélisation des surfaces
et obtient même un doctorat en 1977.
PAUL DE CASTELJAU était lui aussi un mathématicien d'origine, ainsi élève de la
Rue d'ULM a été un temps employé par l'industrie automobile.
Aujourd'hui, les courbes de Bézier (courbe paramétrique aux extrémités imposées avec
des points de contrôle qui dénissent les tangentes à cette courbe à des instants donnés)
sont très utilisées en informatique.
 CHAPITRE 3. COUBES DE BÉZIER ET POLYNÔMES DE BERNSTEIN

3.2 Quelques notions de base sur les courbes de Bézier :

Rappels sur la notion de barycentres
Dénition 3.1 (Barycentre des points pondérés) Soient A et B des points du plan.
Soient α et β des nombres réels tel que α+β ̸= 0. Le barycentre des points pondérés (A, α),
(B, β) est l'unique point G du plan déni par l'égalité :
−→ −−→ →−
αGA + β GB = 0 .

L'ensembre de points du plan tels que :

−−→ −−→ → −
αM A + β M B = 0 .

est : −−→ β −→
AM = AB.
α+β
Cette dernière égalité dénit un unique point dans le plan appartenant à la droite
M
(AB) lorsque A ̸= B .
Dénition 3.2 (Barycentre dans le cas général) Soit n un nombre entier naturel
P 2. Soient A , A , ..., A des points du plan et λ , λ , ..., λ des nombres réels tels que
≥
λ ̸= 0. Le barycentre des points pondérés (A , λ ), (A , λ ), ..., (A , λ ) est l'unique
1 2 n 1 2 n
n

point G du plan déni par l'égalité :

i=1 i 1 1 2 2 n n

n
X −−→ → −
λi GAi = 0
i=1

Ainsi, on montre facilement, en utilisant la relation de Chasles, que :

n
−−→ 1 X −−−→
A1 G = Pn λi A1 Ai
i=1 λi i=2

Ce qui permet d'obtenir à la fois l'existence et l'unicité.

Dénition 3.3 (Isobarycentre.) On dit que G est l'isobarycentre des points A , A , ..., A
lorsque G est barycentre de (A , 1), (A , 1), ..., (A , 1).
1 2 n
1 2 n

3.3 Polynôme de Bernstein

Dénition 3.4 (Polynômes de Bernstein.) Les polynômes de Bernstein sont dénis
par les formules suivantes :
pour
  
n
ti (1 − t)n−i , n ∈ {0, 1, ..., n}.


i

Bi,n (t) =
pour


0, i ̸∈ {0, 1, ..., n}

 21

Rappelons que le coecient  

n n!
=
i i!(n − i)!
est appelé coecient binomial.
Exercice 3.1 Montrer que les B i,n , i = 0, ..., n forment une base de P , puis que
n

Bi,n (t) = (1 − t)Bi,n−1 (t) + tBi−1,n−1 (t), n ≥ 1, i ∈ Z. (3.1)

Théorème 3.1 Soit n un nombre entier naturel. Pour tout t ∈ [0, 1] :
n
X
Bi,n (t) = 1.
i=0

Preuve
On a : n n  
X X n
Bi,n (t) = ti (1 − t)n−i = (t + 1 − t)n = 1
i
i=0 i=0

d'après la formule du binôme de Newton. □

Propriété 3.1 On peut montrer par récurrence que :
B (0) = 1 et ∀i ≥ 1, B
0,n i,n (0) = 0,

ainsi que :
et ∀i ≤ n − 1, B (1) = 0.
Bn,n (1) = 1 i,n

Théorème 3.2 Pour tout entier naturel n ≥ 1, la dérivée du polynôme de Berstien vaut :
d
dt
B (t) = n [B i,n (t) − B (t)]
i−1,n−1 (3.2)
i,n−1

Preuve

d
Bi,n (t) =Cni .[Link]−1 .(1 − t)n−i + Cni .ti .(n − i).(1 − t)n−i−1 .(−1)
dt
n!.i n!(n − i) i
= .ti−1 .(1 − t)n−i − .t .(1 − t)n−i−1
i!(n − i)! i!(n − i)!
n(n − 1)! n(n − 1)!
= .ti−1 .(1 − t)n−1−(i−1) − .ti .(1 − t)n−i−1
(i − 1)!.(n − 1 − (i − 1))! i!(n − i − 1)!
=[Link]−1,n−1 (t) − [Link],n−1 (t)

□
En particulier, nous avons le corrolaire suivant
 CHAPITRE 3. COUBES DE BÉZIER ET POLYNÔMES DE BERNSTEIN

Corollaire 3.1 Pour tout entier naturel n :

B (0) = n et B (1) = n.
′
0,n
′
n,n

Pour tout n ≥ 1 :
B (0) = n et B ′
1,n (1) = −n. ′
n−1,n

Pour tout n ≥ 4 et pour tout i vériant i ≤ n − 2 :

′ ′
Bi,n (0) = Bi,n (1) = 0.

3.4 Courbe de Bézier

3.4.1 Dénitions et propriétés des courbes de Bézier
Dénition 3.5 (Courbe de Bézier.) Soient P , P , ..., P des points du plan P .
1. Soient B , B , ..., B les n + 1 polynômes de Bernstein de degré n. On appelle
0 1 n

courbe de Bézier piloté par les points P , P , ..., P , la courbe Γ décrite par les points
0,n 1,n n,n

M (t), barycentre des points pondérés :

0 1 n

(P0 , B0,n (t)), (P1 , B1,n (t)), ..., (Pn , Bn,n (t))
avec t ∈ [0, 1]. En d'autres termes :
( n
)
X −−−−→ → −
Γ= M (t) ∈ P, ∃t ∈ [0, 1], Bi,n (t)M (t)Pi = 0 .
i=0

2. Les points P , P , ..., P sont les points de contrôle de la courbe de Bézier.

0 1 n

Le théorème ci-dessous nous donne des propriétés des courbes de Bézier :

Théorème 3.3

1. Soit O un point quelconque de P . Alors, pour tout t ∈ [0, 1] :

n n
X −−−−→ → − X −−→ −−−−→
Bi,n (t)M (t)Pi = 0 ⇐⇒ Bi,n (t)OPi = OM (t).
i=0 i=0

2. On suppose que n ≥ 2 et que les points P , P , ..., P ne sont pas tous confondus.
Avec les notations de la dénition 3.5, on a :
0 1 n

(a) M (0) = P et M (1) = P ;

0 n

(b) Soit k le plus petit des nombres entiers i tels que P ̸= P , alors −P−P→ dirige la
tangente à Γ en P ;
0 i 0 k
0

(c) Soit p le plus grand des nombres entiers i tels que P ̸= P , alors −P−−P→, dirige
la tangente à Γ en P .
n i p n
n
 23

[Link] Dérivée d'une courbe de Bézier

A partir de la dérivée d'un polynôme de Bernstein, nous obtenons celle de la courbe
de Bézier :
Proposition 3.1
n−1
d −−−−→ X −−−−→
OM (t) = n. Pi Pi+1 Bi,n−1 (t)
dt i=0

Preuve
n
d −−−−→ X −−→
OM (t) = (Bi,n (t))′ .OPi
dt i=0
n n−1
X −−→ X −−→
=n. Bi−1,n−1 (t).OPi − n Bi,n− (t).OPi
i=1 i=0

en posant j = i − 1, nous obtenons

n−1 n−1
d −−→ X −−−−→ X −−→
OM = n. Bj,n−1 (t).OPj+1 − n. Bi,n−1 (t).OPi
dt j=0 i=0

Donc nous pouvons écrire

n−1 n−1
d −−→ X −−−→ X −−→
OM =n. Bi,n−1 (t).OPi+1 − n. Bi,n−1 (t).OPi
dt i=0 i=0
n−1  n−1
X −−−→ −−→ X −−−−→
=n. OPi+1 − OPi .Bi,n−1 (t) = n. Pi Pi+1 Bi,n−1 (t)
i=0 i=0

D'où la dérivée d'une courbe de Bézier de degré quelconque. □

[Link] Tangentes aux extrémités
Il est possible de voir graphiquement que la tangente au début de la courbe correspond
à la droite reliant P à P . La tangente a la n de la courbe est égale à la droite reliant
P à` P .
0 1

Nous allons le voir en utilisant la formule de la dérivée d'une courbe de Bézier que
n n−1

nous venons d'étudier.

Pour t = 0 : n−1
d −−→ X  −−−−→
OM (t) = n. Pi Pi+1 .Bi,n−1 (0)
dt
Le polynôme de bernstein :
i=0

(n − 1)!
Bi,n−1 (0) .ti .(1 − t)n−i−1 | t = 0
i!(n − i − 1)!
 CHAPITRE 3. COUBES DE BÉZIER ET POLYNÔMES DE BERNSTEIN

vaut 0 sauf quand i = 0 où il vaut 1. La tangente en t = 0 est donc égale à :

d −−→ −−→
OM (0) = n.P0 P1
dt
Ce qui correspond bien a` la droite reliant P à P .
Nous pouvons appliquer la même démarche pour t = 1 et nous trouvons :
0 1

d −−→ −−−−→
OM (1) = [Link]−1 Pn
dt
Les valeurs des tangentes aux extrémités sont logiques ce qui facilite l'utilisation de
courbes de Bézier. Cette propriété est tr`es intéressante pour joindre plusieurs courbes
de Bézier sans rupture de pente.
3.4.2 Construction par barycentres successifs
Théorème 3.4 Soient P , P , ..., P des points du plan. On suppose n ≥ 1. Pour tout
t ∈ [0, 1], on note P (t) le barycentre de
0 1 n

(P0 , B0,n−1 (t)), ..., (Pn−1 , Bn−1,n−1 (t))

et on note Q(t) le barycentre de

(P1 , B0,n−1 (t)), ..., (Pn , Bn−1,n−1 (t)).

Les points P (t) et Q(t) parcourent les courbes de Bézier pilotées respectivemet par P , P , ..., P
et P , P , ..., P quand t parcourt [0, 1].
0 1 n−1

1. Alors, ∀t ∈ [0, 1], le point M (t), barycentre de (P (t), 1 − t) et (Q(t), t) est aussi le
1 2 n

barycentre de
(P , B (t)), ..., (P , B (t)) ;
0 0,n n n,n

autrement dit : le barycentre M (t) de (P (t), 1 − t), (Q(t), t) parcourt la courbe de

Bézier Γ pilotée par les points P , P , ..., P quand t parcourt [0, 1].
2. On a :
0 1 n

d −−−−→ −−−−−→
OM (t) = nP (t)Q(t).
dt
Donc, pour tout t ∈ [0, 1] tel que P (t) ̸= Q(t), la droite (P (t), Q(t)) est tangente à
la courbe de Bézier Γ au point M (t).
Preuve

1. Comme M (t) est le barycentre de (P (t), 1 − t), (Q(t), t), on a :

−−−−−−→ −−−−−−→ → −
(1 − t)M (t)P (t) + tM (t)Q(t) = 0

On en déduit −−−−→ −−−→ −−−→

OM (t) = (1 − t)OP (t) + tOQ(t)
 25

Dans cette dernière égalité, on utilise le fait que les P (t) et Q(t) sont des points
courants de courbes de Bézier. On obtient :
n−1 n−1
−−−−→ X −−→ X −−−−→
OM (t) =(1 − t) Bi,n−1 (t)OPi + Bj,n−1 (t)OPj+1
i=0 j=0
n−1
X (n − 1)! −−→
=(1 − t) tk (1 − t)n−1−k OPk
k=0
k!(n − 1 − k)!
n
X (n − 1)! −−→
+t tk−1 (1 − t)n−1−(k−1) OPk
k=1
(k − 1)!(n − 1 − (k − 1))!
n−1   n  
n−1 k n−k −−→ X n − 1 k −−→
X
= t (1 − t) OPk + t (1 − t)n−k OPk
k k−1
k=0 k=1
n−1     !
−−→ X n − 1 n − 1 −−→ −−→
=(1 − t)n OP0 + + tk (1 − t)n−k OPk + tn OPn .
k k−1
k=1

Comme , on en déduit le théorème; la dernière expres-

     
n n−1 n−1
= +
k−1
sion donne
k k

n−1  
!
−−−−→ −−→ X n −−→ −−→
OM (t) =(1 − t)n OP0 + tk (1 − t)n−k OPk + tn OPn
k
k=1
n
X −−→
= Bk,n (t)OPk .
k=0

2. En écrivant , on obtient
n  
−−−−→ X n k −−→
OM (t) = t (1 − t)n−k OPk
k
k=0

n  
d −−−−→ X n
−−→
OM (t) = ktk−1 (1 − t)n−k − (n − k)tk (1 − t)n−1−k OPk
dt k
k=0
n   n−1  
n k−1 n−k −−→ X n k −−→
X
= k t (1 − t) OPk − (n − k) t (1 − t)n−1−k OPk
k k
k=1 k=0
n   n−1  
n − 1 k−1 n−1−(k−1) −−→ X n−1 k −−→
X
= n t (1 − t) OPk − n t (1 − t)n−1−k OPk
k−1 k
k=1 k=0
−−−→ −−−→
=nOQ(t) − nOP (t)
−−−−−→
=nP (t)Q(t).

□
 CHAPITRE 3. COUBES DE BÉZIER ET POLYNÔMES DE BERNSTEIN

3.5 Exercices
Exercice 3.1 (Courbe de Bezier)

1. On pose t = t = 0, t = 1, t = 2, t = 3, etc.. Calculer B pour k ≤ 3 et

0 ≤ i ≤ 3 − k.
0 1 2 3 4 i,k

2. On pose t = 0, t = 1, t = 2, t = t = 3, t = 4, t = 5, t = 6 et t = 7.
Calculer B pour k ≤ 2 et 0 ≤ i ≤ 5, ainsi que B et B .
0 1 2 3 4 5 6 7 7

Montrer que (sauf pour t entier si k = 0 et t = 3 et k = 1) B (6−t) = B (t).

i,k 0,3 1,3
i,k 6−k−i,k

Exercice 3.2 (Courbes de Hermite et de Bezier)

a) Soient les points P = (1, 1), P = (2, 3), P = (4, 3) et P = (3, 1). Donner
l'équation de la courbe de Bezier approximant ces points.
0 1 2 3

b) Donner l'équation de la courbe interpolant deux points extrémités P et P et de

vecteurs T = (1, 2) et T = (−1, 2) tangents en ces deux points.
0 3

c) A-t-on la même courbe. Quelle condition aurait-il fallu?

0 3

Exercice 3.3 (Polynômes de Bernstein, 10 points) Les polynômes de Bernstein sont

dénis de la matière suivantes :
(3.3)
 
n k n−k
P (x) =
n,k x (1 − x) , 0 ≤ k ≤ n.
k

1. Montrer que la famille (P ) forme une base de l'espace vectoriel des polynômes
de dégré inférieur ou égal à n.
n,k k≤n

2. Montrer que P ≥ 0 sur [0, 1].

n,k

3. Montrer que X P (x) = 1.

n

n,k
k=0

4. Montrer que .
Xn
kPn,k (x) = nx
k=0

5. Montrer que k(k − 1)P (x) = n(n − 1)x .

n
X
2
n,k

6. En utilisant le fait que (nx − k) = k(k − 1) − (2nx − 1)k + n x montrer que

k=0
2 2 2

n
X
(nx − k)2 Pn,k (x) = nx(1 − x).
k=0

Le polynôme de Bernstein d'ordre n associé à f : [0, 1] → R est égal à

(3.4)
n  
X k
Bn (f )(x) = f Pn,k (x).
k=0
n
 27

7.(a) Calculer le polynôme de Bernstein associé à la fonction f dénie par f (x) = x . 2

(b) Existe t-il n tel que B (f ) = f ?

8. Montrer que pour tout n ≥ 1 et pour tout δ > 0, on a
n

   
X k X k
f (x) − Bn (f )(x) = f (x) − f ( ) + f (x) − f ( ) Pn,k (k).
n n
(3.5)
|x−k/n|≤δ |x−k/n|>δ

où les sommes sont faites sur les k variant de 0 à n.

CHAPITRE 3. COUBES DE BÉZIER ET POLYNÔMES DE BERNSTEIN
Chapitre 4
Méthodes itératives de base

Sommaire
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2 Convergence de matrices . . . . . . . . . . . . . . . . . . . . . 30
4.3 Méthodes de relaxation . . . . . . . . . . . . . . . . . . . . . . 32
4.3.1 Généralités sur la méthode . . . . . . . . . . . . . . . . . . . . 33
4.3.2 Cas particuliers . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.3.3 Tests d'arrêt . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.3.4 Les méthodes de Jacobi et de Gauss-Seidel . . . . . . . . . . . 36
4.3.5 Convergence des méthodes de Jacobi et de Gauss-Seidel . . . . 38
4.3.6 Resultats particuliers de convergence . . . . . . . . . . . . . . 40
4.4 La méthode de sur-relaxation . . . . . . . . . . . . . . . . . . . 41
4.4.1 Convergence de la méthode de sur-relaxation . . . . . . . . . . 42
4.4.2 La méthode de sur-relaxation symétrique . . . . . . . . . . . . 43
4.5 Les méthodes de Richardson . . . . . . . . . . . . . . . . . . . 44
4.6 Pseudo-codes et codes Matlab . . . . . . . . . . . . . . . . . . 44
4.6.1 La méthode de Jacobi . . . . . . . . . . . . . . . . . . . . . . . 44
4.6.2 La méthode de Gauss-Seidel . . . . . . . . . . . . . . . . . . . 45
4.6.3 La méthode de sur-relaxation . . . . . . . . . . . . . . . . . . . 45
4.7 Exercice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.1 Introduction
Soit à résoudre le système d'équations linéaires réelles
AX = b,

où A ∈ R et b ∈ R (avec n un entier supérieur ou égal à 2).

n×n n

Une méthode itérative consiste à construire une suite de vecteurs qui, sous certaines
conditions converge vers la solution du système.
 CHAPITRE 4. MÉTHODES ITÉRATIVES DE BASE

4.2 Convergence de matrices

La dénition principale est
Dénition 4.1 Soit (A ) une suite de matrices. On dit que (A ) converge vers A si
k k

lim ∥Ak − A∥ = 0.
k→∞

On dit qu'une matrice carrée A est convergente si la suite de ses puissances converge
vers 0 lorsque k tend vers l'inni.
Pour chercher la condition nécessaire et susante pour que A soit convergente, on consi-
dère la décomposition en Jordan de A.
Il existe une matrice S régulière telle que
SAS −1 = J
 
J1

où J est de la forme de Jordan J2

... avec J une
 
J =  (m = 1, ..., r)
 
 m
 
Jr

. . . . . . 
 
λm 1

matrice n × n de la forme . . . 1  et λ une valeur propre de A


Jm = 

m m m


de multiplicité η telle que λ ̸= λ si p ̸= m. Evidenment, n + · · · + n = n où n est la

λ m

dimension de la matrice A.
n m p 1 r

Par récurrence, on obtient J = SA S , c'est-à-dire que A = S J S. Il est donc

k k −1 k −1 k

équivalent de chercher la condition nécessaire et susante pour que J tende vers 0 quand k

k tend vers l'inni. On a  

J1k
J2k
...
 
Jk = 
 

 
Jrk
Donc A tend vers 0 si et seulement si, pour m = 1, ..., r, J tend vers 0 quand k tend
k k

vers l'inni. Calculons J . On a

k
m
m

λ2m 2λm 1
 
λ2m 2λm
... ... ...

 1 


... ...
 
2
 
Jm =  

...

 1 

 
 2λm 
λ2m
 31

et par récurrence, en notant simplement par d les éléments de J , on trouve que

ij
k
n


 0 si j < i,



k
dij = Cj−i λk−j+i
m si i ≤ j ≤ min(nm , k + i),




0 si k + i < j ≤ nm ,


où C désigne le coecient du binôme (c'est-à-dire C = p!/[q!(p − q)!] avec q ≤ p). Par

p p

conséquent, lorsque k tend vers l'ini, J tend vers 0 si et seulement si |λ| < 1 pour
q
k
q

m = 1, ..., r.
m

Puisque ρ(A) = max |λ | (rayon spectral), nous avons le Théorème suivant qui
regroupe deux (2) propriétés importantes
m m

Théorème 4.1

1. Une condition nécessaire et susante pour que lim A k

=0 est que ρ(A) < 1.
2. La série
k→∞

X ∞
Ak = I + A + A2 + · · ·
k=0

converge vers (I − A) si et seulement si ρ(A) < 1.

−1

Preuve

1. La démonstraction précédente nous renseigne, de plus, sur la vitesse de convergence

vers zéro. En eet
k
∥Jm ∥ ∼ Cnkm −1 |λm |k−(nm −1) (k −→ ∞).

On déduit de ce résultat que

∥J k ∥ ∼ Cp−1
k
[ρ(A)]k−(p−1) ,

où p est la plus grande des dimensions des blocs J pour lesquels λ = ρ(A). De
plus, l'inégalité ∥A ∥ ≤ ∥S∥.∥J ∥.∥S ∥, nous permet d'en déduire.
m
k k −1

2. Supposons que ρ(A) < 1. Soit λ une valeur propre de A. Alors 1 − λ est valeur
propre de I − A. Puisque ρ(A) < 1, on a |λ| < 1 et donc 1 − λ ̸= 0. Par conséquent,
la matrice I − A est régulière.
Posons X k
Ai = I + A + A2 + · · · + Ak .


Sk =
i=0

Alors
(I − A)Sk = I − Ak+1
 CHAPITRE 4. MÉTHODES ITÉRATIVES DE BASE

c'est-à-dire
Sk = (I − A)−1 (I − Ak+1 ).
On a donc
Sk − (I − A)−1 = −(I − A)−1 Ak+1
ce qui, en passant aux normes conduit à
∥Sk − (I − A)−1 ∥ ≤ ∥(I − A)−1 ∥.∥Ak+1 ∥.

D'après le point 1., la suite (A ) converge vers zéro; ce qui démontre que la
k

suite (S ) converge vers (I − A) .

k≥0
−1

Réciproquement, si la suite (S ) admet une lmite, alors la suite (A ) converge

k k≥0
k

vers zéro. Donc, d'après le point 1., ρ(A) < 1.

k k≥0 k≥0

4.3 Méthodes de relaxation

Une méthode de relaxation pour résoudre le système de n équations à n inconnues
Ax = b

consiste à eectuer une décomposition de la matrice A sous la forme

A=M −N

où M est une matrice régulière. Le système peut alors s'écrire

M x = N x + b.

Construisons la suite de vecteurs x , x , x , ... de la façon suivante :

• x arbitraire,
0 1 2

= N x + b, soit encore x = M N x + M b, pour k = 0, 1, 2, ....

0
−1 −1
• Mx
Une telle méthode itérative s'appelle une méthode de relaxation. On la qualie souvent,
k+1 k k+1 k

de méthode itérative stationnaire car la matrice M N est indépendante de k. On peut

−1

mettre les itérations précédentes sous une autre forme. Ainsi, on a

M −1 N = I − M −1 A,

d'où
(4.1)
xk+1 = xk + M −1 (b − Axk ) = xk + M −1 rk
avec r = b − Ax . La matrice M apparaît donc comme une approximation de A ,
−1 −1

c'est-à-dire comme un préconditionneur.

k k

Considérons l'erreur
e k = x − xk .
 33

On a
ek = M −1 N ek−1 ,
d'où par récurrence,
ek = (M −1 N )k e0 .
D'autre part r k = Aek . Donc
rk =(M − N )M −1 N ek−1
=N (I − M −1 N )ek−1
=N M −1 (M − N )ek−1
=(N M −1 )rk−1 .

Il s'en suit par récurrence que
k

rk = N M −1 r0 .
On rappelle que les matrices M N et N M sont semblables et que si u est vecteur
−1 −1

propre de la prémière matrice, alors v = M u est le vecteur propre de la seconde associé à

la même valeur propre.
Remarque 4.1 Il est possible de considérer de façon plus générale, une méthode itérative
de la forme
xk+1 = Bk xk + ck
Si x = x, il est nécéssaire d'avoir x = x. C'est le cas si c = H b et I − B = Hk A
où H est une matrice inversible. Les itérations s'écrivent alors
k k+1 k k k
k

xk+1 =xk + (Bk − I)xk + ck

=xk − Hk Axk + Hk b
=xk + Hk rk .

Hk est donc également un préconditionneur. cette méthode revient, en fait, à considérer

une décomposition A = M − N qui change à chaque itération avec M = H . Une telle −1

méthode est qualiée de non stationnaire.

k k k k

4.3.1 Généralités sur la méthode

De l'égalité e k = (M −1 N )k e0 , on a
∥ek ∥ ≤ ∥M −1 N ∥k .∥e0 ∥

d'où ce premier théorème de convergence

Théorème 4.2

1. Si ∥M −1
N∥ < 1 alors (x ) converge vers x quelque soit x .
k k≥0 0
 CHAPITRE 4. MÉTHODES ITÉRATIVES DE BASE

2. Une condition nécéssaire et susante pour que la suite (x ) converge vers x

quelque soit x est que ρ(M N ) < 1.
k k≥0
−1
0

Plus ρ(M N ) < 1 est petit et plus la convergence de la méthode est rapide. Le choix
−1

de la décomposition de la matrice A du système est donc guidé par les considérations

suivantes :
1. Il faut que ρ(M N ) < 1,
−1

2. la résolution de M y = c, où c est un vecteur quelconque, doit être simple et exiger

le moins d'opérations possibles,
3. la décomposition retenue doit être la meilleure possible, c'est-à-dire que ρ(M N )
−1

doit être le plus petit possible.

4.3.2 Cas particuliers
On a
M −1 N = (A + N )−1 N = (I + G)−1 G avec G = A N. −1

Soit u un vecteur propre de G correspondant à la valeur propre τ . Alors 1 + τ ̸= 0 et l'on

a
Gu = τ u et (I + G) Gu = −1 τ
u.
1+τ
Par conséquant u est aussi vecteur propre de M N et il correspond à la valeur propre
−1

µ = τ /(1 + τ ).
Réciproquement, si µ est valeur propre de (I + G) G et si ν est le vecteur propre
−1

correspondant, alors Gν = µ(I + G)ν .

D'après cette relation, il ne peut pas être égal à 1 et donc
µ
Gν = ν = τ ν.
1−µ

Puisque µ = τ /(1 + τ ) est une fonction strictement croissante de τ , il découle

ρ(A−1 N )
ρ(M −1 N ) = .
1 + ρ(A−1 N )

Il existe des résultats généraux de convergence dans un certain nombre de cas particuliers.
Commençons par les dénitions suivantes
Dénition 4.2

1. Une matrice A est dite non négative si

aij ≥ 0, ∀i, j

On écrira A ≥ 0 où 0 est la matrice nulle. On dénirait de même une matrice

positive par l'inégalité stricte. (A ≥ B signie que A − B ≥ 0)
 35

2. Soit une décomposition M − N de la matrice A (i.e. A = M − N ) où M et

N sont des matrices de la même dimension que la matrice A. On dit que cette
décomposition est régulière si M est régulière, si M ≥ 0 et si N ≥ 0. On dit −1

que la décomposistion est faiblement régulière si M est régulière, si M ≥ 0 et si −1

M N ≥ 0.
−1

Nous avons les résultas de convergence suivants, dans des cas particuliers
Théorème 4.3

1. Soit A = M − N une décomposition régulière de A. Alors,

A est régulière et A ≥ 0 si et seulement si
−1

ρ(M −1 N ) =
ρ(A−1 N )
1 + ρ(A−1 N )
< 1. (4.2)
Donc, si A est régulière avec A ≥ 0, alors (x ) converge vers x quel que soit x .
−1

2. Soit A = M − N une déécomposition faiblement régulière de A. Alors,

k 0

A est régulière et A ≥ 0 si et seulement si ρ(M N ) < 1.

−1 −1

3. Soient A = M − N = M − N deux décompositions régulières de A. Si A ≥ 0 −1

et si N ≥ N ≥ 0, alors
1 1 2 2
2 1

0 ≤ ρ(M1−1 N1 ) ≤ ρ(M2−1 N2 ) < 1.

De plus, si A > 0, si N ≥ N ≥ 0 et si N ̸= N ̸= 0, alors les inégalités

−1

précédentes sont strictes.

2 1 2 1

4. Soient A = M − N = M − N deux décompositions régulières de A. Si A ≥ 0 −1

et si M ≥ N , alors
1 1 2 2
−1 −1
1 2

0 ≤ ρ(M1−1 N1 ) ≤ ρ(M2−1 N2 ) < 1.

Si A −1
>0 et si M −1
1 > M2−1 , alors les inégalités précédentes sont strictes.
Preuve
En Exercice. □

4.3.3 Tests d'arrêt

On a
xk+1 − xk =(M −1 N − I)xk + M −1 b
=(M −1 N − I)(xk − x)

d'où
∥x − xk ∥ ≤ ∥(I − M −1 N )−1 ∥.∥xk+1 − xk ∥.
 CHAPITRE 4. MÉTHODES ITÉRATIVES DE BASE

Mais puisque ∥M −1
N∥ < 1, on a
(I − M −1 N )−1 =I + (M −1 N ) + (M −1 N )2 + · · ·
1
∥(I − M −1 N )−1 ∥ ≤1 + ∥M −1 N ∥ + ∥M −1 N ∥2 + · · · ≤
1 − ∥M −1 N ∥
et l'on obtient nalement
1
∥x − xk ∥ ≤ ∥xk+1 − xk ∥
1 − ∥M −1 N ∥
∥M −1 N ∥
∥x − xk+1 ∥ ≤ ∥xk+1 − xk ∥.
1 − ∥M −1 N ∥

On arrêtera les itérations lorsque la norme ∥x −x ∥ sera inférieure à un seuil xé.Cependant,

si ∥M N ∥ est voisin de 1, la norme de l'erreur pourra être grande. Il sera alors nécessaire
k+1 k
−1

de tenir compte du facteur 1 − ∥M N ∥ qui intervient dans ces inégalités. Puisque

−1

∥xk+1 − xk ∥/∥xk − xk−1 ≤ ∥M −1 N ∥,

ce rapport pourra servir d'approximation pour la norme de M N . −1

4.3.4 Les méthodes de Jacobi et de Gauss-Seidel

Sour D une matrice diagonale constituée des éléments diagonaux de A
 
...
a11
D= .
 
ann

Soit −E la matrice formée par la partie triangulaire stritement inférieure de A

 
0

.. .0. . . . .
 a21 
 
−E = 
.. ... ...
 

 
 
an1 · · · · · · an,n−1 0

et soit −F la matrice de la partie triangulaire strictement supérieure de A

... ...
0 a12 · · · · · · a1n
..
 

... ... ..
 
−F = 
 

 
0 ann

On a donc A = D − E − F .On supposera que la matrice D n'est pas singulière, c'est-à-dire

que ∀i, a ̸= 0.
ii
 37

An de construire une méthode de relaxation, nous allons maintenant regrouper deux de
ces matrices et mettre A sous la forme A = M − N .
On peut également dénir les méthodes de relaxation par blocs en partitionnant la
matrice A en blocs. On a alors D = diag(A , ..., A ) où les blocs diagonaux A sont
carrés et réguliers. Les matrices par blocs E et F sont dénies de manière analogue à ce
11 nn i

qui vient d'être fait.

[Link] La méthode de Jacobi
Cette méthode, qui date de 1845, est due à Jacobi [7], Elle consiste à prendre M = D
et N = E + F . D'où
M −1 N = D−1 (E + F )
c'est-à-dire  a12 a1n 
0 − ··· ··· −


−
a21
a11
0
... ..
a11 


.. ... ... ... ..

 
 a22 
−1
 
..
M N =
... ...

 
 a n−1,n 
 − 
 an−1,n−1 
 an1 an,n−1 
− ··· ··· − 0
ann ann
ce qui nous donne les itérations suivantes
xk+1 = D−1 (E + F )xk + D−1 b (4.3)
qui correspond, pour les composantes x des vecteurs x , à
k,i k

(4.4)
n
!
1 X
xk+1,i = bi − aij xk,j , i = 1, 2, 3, ..., n.
aii j=1,j̸=i

Or nous avons, pour i = 1, 2, ..., n,

n
X
rk,i =bi − aij xk,j
j=1
Xn
=bi − aij xk,j − aii xk,i
j=1,j̸=i

d'où
rk,i
xk+1,i =xk,i +
aii
n
X aij
rk+1,i = − rk,j .
j=1,j̸=i
ajj
 CHAPITRE 4. MÉTHODES ITÉRATIVES DE BASE

On a aussi M N = I − D A.
−1 −1

La méthode de Jacobi est vraiment la plus simple à laquelle on puisse penser. Cependant,
sa convergence peut être très lente. Si on l'applique, par exemple, à une matrice tridia-
gonale de dimension 100 avec des 2 sur la diagonale des −1 sur les diagonales inférieure
et supérieure et si le second membre est nul, alors, en partant d'un vecteur x aléatoire
entre 0 et 1, il faut de l'ordre de 28.000 itérations pour avoir une erreur de 10 .
0
−6

[Link] La méthode de Gauss-Seidel

Cette méthode est due à Gauss (1826) et à Seidel (1874). Elle consiste à prendre
M = D − E et N = F . D'où

M −1 N = (D − E)−1 F

ce qui nous donne l'itération suivante

xk+1 = (D − E)−1 F xk + (D − E)−1 b (4.5)
qui correspond, pour les composantes, à
(4.6)
i−1 n
!
1 X X
xk+1,i = bi − aij xk+1,j − aij xk,j , i = 1, ..., n.
aii j=1 j=i+1

Remarque 4.2 On voit que la méthode de Jacobi demande de conserver en mémoire

deux vecteurs alors que la méthode de Gauss-Seidel n'en demande qu'un seul. En eet,
dans la méthode de Jacobi, le calcul d'une composante du vecteur x fait appel à toutes
les composantes de x alors que, dans la méthode de Gauss-Seidel, le calcul de la i-ième
k+1

composante de x utilise seulement les composantes i + 1 à n du vecteur x .

k
k+1 k

4.3.5 Convergence des méthodes de Jacobi et de Gauss-Seidel

Dans la suite, nous noterons B = (b ) la matrice M N de la méthode de Jacobi (
−1

il ne faut pas confondre ces éléments b avec les composantes b du second membre b).
J ij

Nous avons donc

ij j

BJ = D−1 (E + F ).
Comme nous l'avons vu, c'est une matrie à diagonale nulle, remarque qui aura son im-
portance par la suite. Nous la décomposons en une somme B = L + U où
 la matrice L est srictement triangulaire inférieure (c'est-à-dire que sa diagonale
J

principale ne contient que des termes nuls);

 la matrice U est strictement tringulaire supérieur.
On a
L=D E et U = D F.
−1 −1

Soit L la matrice M N de la méthode de Gauss-Seidel, c'est-à-dire

1
−1

L1 = (D − E)−1 F = (I − F )−1 U.
 39

[Link] Résultats préliminaires

Dénition 4.3 A est réductible s'il existe une matrice de permutation P telle que
 
T A11 A12
P AP =
0 A22

où les sous-matrices A et A sont carrées. Une matrice qui ne vérie pas cette propriété
est appélée irréductible.
11 22

Soit le théorème suivant de Perron-Frobenius. (rappelons que ρ(A) désigne le rayon

spectrale de A, c'est-à-dire ρ(A) = max |λ |, où les λ sont les valeurs propres de A)
i i i

Théorème 4.4 (de Perron-Frobenius) Si A ≥ 0 et est irréductible, alors

1. A a une valeur propre positive ou nulle égale à ρ(A),
2. il existe un vecteur propre positive correspondant à ρ(A) dont toutes les compo-
santes sont strictement positives,
3. ρ(A) croît quand un élément quelconque de A croît;
4. ρ(A) est une valeur propre simple.
Preuve
Pour la démonstraction de théorème, voir [12]. □
[Link] Convergence des méthodes
Le résultat suivant est connu sous le nom du théorème de Stein-Rosenberg (1948)
Théorème 4.5 (de Stein-Rosenberg) Supposons que B ≥ 0. Alors un et seulement
un des cas suivants (qui s'excluent donc mutuellement) est vérié
J

1. ρ(B ) = ρ(L ) = 0
2. 0 < ρ(L ) < ρ(B ) < 1
J 1

3. ρ(B ) = ρ(L ) = 1,
1 J

4. 1 < ρ(B ) < ρ(L ).

J 1

J 1

Preuve
Pour la démonstraction, voir [4] □
Ainsi, sous les hypothèses du Théorème de Stein-Rosenberg, les matrices B et L
sont simultanément convergentes ou non.
J 1

Dans le cas de :
1. la convergence, la méthode de Gauss-Seidel converge plus vite que celle de Jacobi;
2. la divergence, la méthode de Gauss-Seidel diverge plus vite que celle de Jacobi.
 CHAPITRE 4. MÉTHODES ITÉRATIVES DE BASE

Corollaire 4.1 Si la matrice non négative B est telle que 0 < ρ(B ) < 1, alors
J J

0 < ρ(L1 ) < ρ(BJ ) < 1. (4.7)

Une notion importante est donnée dans la Dénition suivante
Dénition 4.4 On dit que A est à dominance diagonale si

(4.8)
n
X
|aii | ≥ |aij |, i = 1, ..., n.
j=1,j̸=i

Si l'inégalité est stricte, on parle de dominance diagonale stricte.

S'il existe une matrice diagonale D > 0 telle que D AD soit dominance diagonale−1

(strite), on dit que A est à dominance diagonale (stricte) généralisée.

Pour les matrices à dominance diagonale stricte, on a
Théorème 4.6 Si la matrice A est non singulière et à dominance diagonale stricte, alors
les méthodes de Jacobi et de Gauss-Seidel sont convergentes.
Preuve
Voir [4]. □

4.3.6 Resultats particuliers de convergence

Soit cette première dénition
Dénition 4.5 On dit que la matrice A est monotone si et seulement si elle est
régulière et A ≥ 0.−1

Par conséquent, si b ≥ 0 alors x = A b ≥ 0. −1

Cette dénition est la version en dimension nie du principe du maximum.

On a d'abord le résultat suivant
Théorème 4.7 Si a ≤ 0, ∀i ̸= j , alors les deux propositions suivantes sont équiva-
lentes
ij

1. A est régulière et A ≥ 0, −1

2. ∀i, a > 0, B = I − D A ≥ 0 est irréductible et ρ(B ) < 1.

ii J
−1
J

Soit cette deuxième Dénition

Dénition 4.6

1. Une matrice régulière A telle que a ≤ 0, ∀i ̸= j et A −1

≥ 0 s'appelle une M -
matrice.
ij
 41

2. Soit M (A) la matrice dont les éléments m sont donnés par m = |a | et m =

−|a | pour i ̸= j . On dit que A est une H -matrice si et seulement si M (A)
ij ii ii ij

est une M -matrice.

ij

3. Une matrice A symétrique dénie positive et telle que a ≤ 0, ∀i ̸= j, s'appelle une

matrice de Stieltjes.
ij

Nous voyons qu'une M -matrice est une H -matrice.

Une conséquence du Théorème précédent est donnée ci-dessous
Corollaire 4.2 Si A est une matrice de Stieltjes, alors c'est une M -matrice. De plus, A
est irréductible si et seulement si A > 0.
−1

Nous énumérons un ensembles de résultats dans le théorème suivants :

Théorème 4.8 1. Soit A une M -matrice et soit C une matrice obtenue en rempla-
çant certains éléments hors diagonaux de A par 0. Alors C est une M -matrice.
2. A est une H -matrice non singulière si et seulement si A est à dominance diagonale
stricte généralisée.
3. Si A est une matrice de Stieltjes et si A = M − N est une décomposition régulière
de A avec N symétrique, alors
ρ(N )ρ(A−1 )
ρ(M −1 N ) ≤ < 1.
1 + ρ(N )ρ(A−1 )

4. Si A est une M -matrice et si M est obtenue en remplaçant certains éléments hors

diagonaux de A par 0, alors A = M − N est une décomposistion régulière de A et
ρ(M N ) < 1.
−1

5. Si A est une M -matrice ou H -matrice, alors les méthodes de Jacobi et de Gauss-

Seidel sont convergentes.
4.4 La méthode de sur-relaxation
Cette méthode est due à David Young (en 1952). Elles consiste à prendre
1
M = (D − ωE)
ω
1
N = [(1 − ω)D + ωF ]
ω
avec ω ∈ R . Nous pouvons facilement vérier que A = M − N .
∗

On pose
Lω = M −1 N = (D − ωE)−1 [(1 − ω)D + ωF ].
 CHAPITRE 4. MÉTHODES ITÉRATIVES DE BASE

La méthode consiste à calculer la i-ième composante du vecteur intermédiaire xe puis,

immédiatement après, la i-ième composante du vecteur x par les relations
k+1
k+1

i−1 n
! 
1 X X 
x
ek+1,i = bi − aij xk+1,j − aij xk,j 

aii

j=1 j=i+1 i = 1, 2, ..., n.



xk+1,i − xk,i ] = (1 − ω)xk,i + ωe
xk+1,i = xk,i + ω[e xk+1,i .


Le paramère ω s'appelle le parmètre de rélaxation (ou facteur de relaxation).

 Lorsque ω < 1, on parle de sous-relaxation
 tandis que pour ω > 1 on parle de sur-relaxation
 Si ω = 1, on retrouve la méthopde de Gauss-Seidel.
On appelle souvent cette méthode SOR (c'est-à-dire successive Over Relaxation ).
Posons
L = D E et U = D F.
−1 −1

L'itération prend la forme

xk+1 = (I − ωL)−1 [(1 − ω)I + ωU ]xk + ω(I − ωL)−1 D−1 b.

On a aussi
Lω = (I − ωL)−1 [(1 − ω)I + ωU ].

4.4.1 Convergence de la méthode de sur-relaxation

Nous régroupons dans ce théorème ci-dessous quelques résultats sur la convergence de
la méthode de sur-relaxation
Théorème 4.9

1. (Théorème de Reuch-Ostrowski) Soit B une matrice avec une diagonale nulle. Met-
tons B sous la forme B = L + U , L étant une matrice triangulaire strictement
inférieure et U une matrice triangualire supérieure.
Posons −1
L = (I − ωL) [(1 − ω)I + ωU ].
ω

Alors quels que soit ω réel, ρ(L ) ≥ |ω − 1|, l'égalité ne se produisant que si toutes
les valeurs propres de L sont de module |ω − 1|.
ω

2. Si la méthode de sur-relaxation converge, alors 0 < ω < 2.

ω

3. Si A est hermitienne dénie positive, alors ρ(L ) < 1 si et seulement si 0 < ω < 2.
ω

Preuve
Pour les preuve, voir dans [4]. □
Remarque 4.3
 43

1. La matrice B de la méthode de Jacobi est justement de la forme de la matrice B

du théorème de Reich-Ostrowski ci-desssus.
J

2. On montre que si l'on suppose A hermitienne et D dénie positive, alors la condi-

tion A dénie positive est nécéssaire.
Finissons la section par un dernier résultat
Théorème 4.10 Si toutes les valeurs propres de B = D (E + F ) sont réelles, alors
−1

les méthodes de Jacobi, Gauss-Seidel et SOR (avec 0 < ω < 2) convergent ou divergent
J

simultanément.
4.4.2 La méthode de sur-relaxation symétrique
Lorsque la matrice A est symétrique, il existe une version symétrique de la méthode
de sur-relaxation, c'est la méthode de sur-relaxation symétrique, également connue sous
l'acronyme SSOR. Elle consiste a prendre
1
M= (D − ωE)−1 (D − ωF )
ω(2 − ω)

avec N = M − A et E = F .
Les itérations de cette méthode peuvent également s'écrire
xk+k/2 =M1−1 N1 xk + M1−1 b

xk+1 =M2−1 N2 xk+1/2 + M2−1 b

avec
M1 = (D − ωE)/ω , M2 = (D − ωF )/ω.

En éliminant x k+1/2 , il vient

xk+1 = M2−1 N2 M1−1 N1 xk + M2−1 (N2 M1−1 + I)b.

On a les résultas suivants :

Théorème 4.11

1. Si A est symétrique dénie positive et si D est dénie positive, alors la méthode

SSOR converge pour 0 < ω < 2.
2. Si A est une H -matrice, alors la méthode SSOR converge si 0 < ω < 2/(1+ρ(|B |).
J
 CHAPITRE 4. MÉTHODES ITÉRATIVES DE BASE

4.5 Les méthodes de Richardson

Elles regroupent plusieurs méthodes de la forme
xk+1 = xk + λk rk

avec x donné. Le vecteur r peut être calculé soit directement par la formule r = b−Ax ,
soit de façon itérative par
0 k k k

r = r − λ Ar .
k+1 k k k

Ces méthodes se distinguent par le choix du paramètre λ . Le choix le plus simple consiste
à prendre λ = λ.
k

On parle alors de méthode de richardson stationnaire. Puisque

k

 
1 1
xk+1 = I − A xk + b,
λ λ
il correspond à la décomposistion M = I/λ et N = I/λ − A.
Cette méthode converge donc si et seulement si ρ(I − λA) < 1. D'où le théorème suivant
Théorème 4.12 Si A est symétrique dénie positive, alors la méthode de Richardson
stationnaire converge si et seulement si 0 < λ < 2/ρ(A).
Il existe dans ce cas une valeur optimale de λ. Elle est donnée par le théorème ci-dessous
Théorème 4.13 Si A est symétrique dénie positive, alors la valeur optimale de λ pour
la méthode de Richardson stationnaire est λ = 2/(λ + λ ), où λ et λ = ρ(A) sont
respectivement la plus petite et la plus grande des valeurs propres de A.
opt m M m M

Preuve
pour la preuve, voir [4]. □
Ce théorème montre que plus κ (A) est voisin de 1 plus la convergence est rapide.
2

4.6 Pseudo-codes et codes Matlab

4.6.1 La méthode de Jacobi
On donne :
en entrée, la matrice A, le second mambre b, le vecteur initial x , la dimension n du sys-
tème, la tolérance toll (qui représente ε) pour le test d'arrêt, et le nombre maximum d'ité-
0

rations k . En sortie, l'algorithme fournit l'approximation de la solution x = (x , ..., x ) T

et le nombre k d'itérations eectuées. l'algorithme se termine soit quand on a atteint le

max 1 n

nombre maximum d'itération k où quand la solution a été obtenue avec la précision

désirée.
max

L'algorithme nécéssite un vecteur auxilliaire x = (x , ..., x ) ) pour garder en mémoire

T

les composantes x du nouvel itéré x qui sont calculées à partir des composantes x de
1 n

l'ancien itéré x.
i j
 45

[x, k] =Jacobi(A, b, x0 , toll, kmax )

k=0
x=0
test = toll + 1
while test ≥ toll and k < kmax do
k =k+1
f or i = 1, ..., n
xi = 0
f or j = 1, ..., n
xi = xi + aij xj
end f or j
f or j = i + 1, ..., n
xi = xi + aij xj
end f or j
xi = (bi − xi )/aii
end for i
test = ∥x − x∥
x=x
end while

4.6.2 La méthode de Gauss-Seidel

Les entrées et les sorties sont les mêmés que le pseudo-code de la méthode de Jacobi.
[x, k] = Gauss− Seidel(A, b, x0, toll, kmax )
k=0
xx 0
r = b − Ax
test = ∥r∥
toll = toll.∥b∥
while test ≥ toll and k < kmax do
k =k+1
f or i = 1, ..., n
xi = 0
f or j = 1, ..., n
xi = xi + aij xj
end for j
f or j = i + 1, ..., n
xi = xi + aij xj
end for j
xi = (bi − xi )/aii
end for i
r = b − Ax
test = ∥r∥
end while

4.6.3 La méthode de sur-relaxation

Dans les arguments d'entrées, par rapport aux algorithmes de Jacobi et de Gauss-
Seidel, on aura besoin en plus du paramètre ω ∈]0, 2]. Cet algorithme, bien qu'il coïncide
théoriquement avec celui de Gauss-Seidel quand ω = 1, nécéssite, comme la méthode de
Jacobi, l'utilisation d'un vecteur auxiliaire x = (x , ..., x ) dans le cas général. Le test
1 n
T
 CHAPITRE 4. MÉTHODES ITÉRATIVES DE BASE

d'arrêt est eectué sur le résidu rk = b − Axk .

[x, k] = SOR(A, b, x0, n, toll, kmax , ω)
k=0
r = b − Ax
test = ∥r∥
toll = toll.∥A∥
while test ≥ toll and k < kmax do
k =k+1
f or i = 1, ..., n
xi = 0
f or j = 1, ..., i − 1
xi = xi + aij xj
end for j
f or j = i + 1, ..., n
xi = xi + aij xj
end for j
xj = (bi − xi )/aii
xi = xi + ω(xi − xi )
end r i
x=x
r = b − Ax
test = ∥r∥
end while

4.7 Exercice
Exercice 4.1 On donne la matrice
 
2 −1 0 0
 −1 2 −1 0 
A= 
 0 −1 2 −1 
0 0 −1 2

1. Former la décomposition LU de A où L est une matrice triangulaire inférieure

dont les éléments diagonaux sont égaux à l'unité et U est une matrice triangulaire
supérieur. Calculer D = det(A).
2. Utiliser la factorisation précédente pour résoudre le système linéaire Ax = b, avec
b = [1, 1, 1, 1] .
T

3. Former la décomposition de Cholesky de A, qui s'écrit A = BB , avec B une ma- T

trice triangulaire inférieure dont les éléments diagonaux peuvent être quelconques.
Á l'aide de ce résultat, calculer à nouveau le déterminant D.
Chapitre 5
Calcul de valeurs propres

Sommaire
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.2 La méthode de la puissance . . . . . . . . . . . . . . . . . . . . 48
5.2.1 Accélération de la convergence . . . . . . . . . . . . . . . . . . 51
5.2.2 Variantes de la méthode de la puissance . . . . . . . . . . . . . 52
5.2.3 Recherche des autres valeurs propres . . . . . . . . . . . . . . . 54
5.2.4 La méthodes de la puissance inverse . . . . . . . . . . . . . . . 56
5.3 Méthodes de décomposition . . . . . . . . . . . . . . . . . . . . 57
5.3.1 La méthode de Jacobi . . . . . . . . . . . . . . . . . . . . . . . 59
5.3.2 La méthode Greenstadt . . . . . . . . . . . . . . . . . . . . . . 60
5.3.3 L'algoritme LR . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.3.4 l'Agorithme QR . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.4 Pseudo-codes et Codes Matlab . . . . . . . . . . . . . . . . . . 65
5.4.1 L'algorithme de Jacobi, choix cyclique . . . . . . . . . . . . . . 65
5.4.2 L'algorithme de Jacobi, choix classique . . . . . . . . . . . . . 65
5.4.3 L'algorithme QR . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

5.1 Introduction
Le calcul des valeurs propres de matrices internient dans de nombreux domaines des
mathématiques appliquées. Citons, par exemple, la résistence des matériaux et le calcul
des structures, l'analyse des phénomènes vibratoires, les chaînes de Markov et les moteurs
de recherche sur le Web, les modèles économiques, l'analyse des données, la physique et
la chimie quantique, etc. C'est donc un domaine particulièrement important à étudier.
Il existe deux grandes classes de méthodes itératives pour le calcul numérique des
valeurs propres d'une matrice :
1. Les méthodes qui permettent de ne calculer qu'une seule valeur propre à la fois (
en général celle de plus grand module).
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

2. les méthodes qui permettent de calculer simultanément toutes les valeurs propres.
Ces dernières ce divisent elles-mêmes en deux catégories :
(a) Les méthodes de détermination de polynôme caractéristique. Ce sont les moins
intéressantes (sauf dans certains cas particuliers) car il faut ensuite calculer les
racines du polynôme caractéristique,
(b) les méthodes de décomposition qui ont un caractère essentiellement itératif et
qui sont basée sur l'utilisation de transformations semblables.
Enn on verra l'intérêt que présente une classe de matrices spéciales : les matrices de la
forme de Hessenberg. Notons qu'l existe beaucoup d'autres méthode que celles décrites
dans ce Chapitre.

5.2 La méthode de la puissance

La méthode de puissance est généralement utilisée pour calculer la valeur propre de
plus grand module d'une matrice ainsi que le vecteur propre correspondant. Dans certains
cas, et après avoir obtenu cette valeur propre et le vecteur propre correspondant, il est
possible de calculer la seconde valeur propre et ainsi de suite.
Cette méthode consiste à eectuer les itérations suivantes :
uk+1 = Auk , k = 0, 1, ...

à partir d'un vecteur donné.

On peut également considérer la suite
Sk = (uk , Auk )/(uk , uk ), k = 0, 1, ....

On parle alors de la méthode du quotient de Rayleigh.

Sous certaines hypothèses, nous allons démontrer que (S ) converge vers la valeur propre
de plus grand module de la matrice A.
k

Cas 1 :

Supposons que les valeurs popres de A forment une base et que

|λ1 | > |λ2 | ≥ ... ≥ |λm |.

Soit u ∈ R un vecteur quelconque tel que (u , x ) ̸= 0 où x est le vecteur propre

n

correspondant à λ . Écrivons u dans la base formées par les vecteurs propres x , ..., x et
0 0 1 1
1 0 1 n
 49

eectuons les itérations. On obtient

n
X
u0 = ai x i
i=1
n
X n
X
u1 =Au0 = ai Axi = ai λ i x i

..
i=0 i=1

n
X
k
uk =Auk−1 = A u0 = ai λ k x i
i=1
n
X
=a1 λk1 + λk1 ai (λi /λ1 )k xi .
i=2

Soit y un vecteur propre tel que (y, x ) ̸= 0. Calculons le rapport S = (y, uk+1 )/(y, uk ) .
On obtient
1 k

a1 (y, x1 ) + ni=2 ai (λi /λ1 )k+1 (y, xi )

P
Sk = λ1 .
a1 (y, x1 ) + ni=2 ai (λi /λ1 )k (y, xi )
P

Or |λ /λ | < 1 pour i = 2, ..., n puisque |λ | > |λ |, i = 2, 3, ..., n , d'où lim k7→∞ (λi /λ1 )
k
=
0 pour tout i = 2, 3, ..., n et par conséquent
i 1 1 i

lim Sk = λ1 .
k7→∞

On voit que la vitesse de convergence est réglée par le rapport λ /λ . En eet 2 1

Sk = λ1 + O[(λ2 /λ1 )k ].

D'autre part
Sk uk+1 = Sk Auk
soit encore uk+1 kk
Sk =A
(y, uk+1 ) (y, uk )
ce qui montre que, sous les mêmes hypothèses
uk
lim = x1 .
k7→∞ (y, uk )

Cette convergence se démontre aisément à partir de l'expression de u et de celle de (y, u ).

Nous allons maintenant étudier plusieurs cas où les hypothèses précédentes ne sont
k k

plus satisfaites.
Cas 2 :
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

Supposnons toujours que x , x , ..., x forment une base et mais que

1 2 n

λ = λ = ··· = λ
1 2 et |λ | > |λ | ≥ · · · ≥ |λ |.
r r r+1 n

De la même façon, on démontrerait que

(y, uk+1 )
= λ1 + O (λr+1 /λ1 )k
 
Sk =
(y, uk )

Cas 3 :
Si l'on a maintenant le cas où λ 1 = −λ2 avec |λ | > |λ | ≥ · · · ≥ |λ |, on voit que
2 2 n

n
X
u2k =(a1 x1 + a2 x2 )λ2k
1 + qi λ2k
1 xi
i=3
n
X
u2k+1 = (a1 x1 − a2 x2 )λ2k+1
1 + ai λi2k+1 xi
i=3

d'où (y, u2k+2 ) (y, u2k+3 )

lim = lim = λ21 .
k7→∞ (y, u2k ) k7→∞ (y, u2k+1 )

Cas 4 :
Considère maintenant le cas où λ est complexe. On a alors λ = λ1 et, comme
précédemment,
1 2

n
X
uk = a1 λk1 x1 + a2 λk2 x2 + ai λki xi
i=3

où a = a puisque u esr réel et x = x . Posons a = ρe , Pa iα

2 = ρe
−iα
, λ = re ,iθ

et y = Re(x ). On trouve u = 2ρr cos(kθ + α)y + . Formons le

2 1 k 2 1 1 1
−iθ k n k
λ = re i=3 ai λi xi
rapport
2 1 1 k 1

(y, uk+1 ) (y, uk+2 )

(y, uk+2 ) (y, uk+3 )
Sk′ =
(y, uk ) (y, uk+1 )
(y, uk+1 ) (y, uk+2 )
On trouve que Sk′ = λ1 λ2 + O (λ3 /λ1 )k = r2 + O (λ3 /λ1 )k
   
|λ2 | > |λ3| ≥ · · · ≥ |λn | si ,
ce qui nous fournit . Calculons :
r θ
On a cos[(k + 1)θ + α] + cos[(k − 1)θ + α] = 2 cos(kθ + α) cos θ d'où
r(y, uk−1 ) + r−1 (y, uk+1 )
= cos θ + O (λ3 /λ1 )k .
 
2(y, uk )

Cas 5 :
 51

Considérons enn le cas où les vecteurs propres de A ne forment plus une base. Pour
simplier nous supposons que A peut se mettre sous la forme de Jordan suivante
 
λ1 1
 0 λ1 
 
...

 λ3 

 
 
λn

avec |λ | > |λ | ≥ · · · ≥ |λ |.
Puisque les vecteurs propres ne forment plus une base, appelons x , ..., x les vecteurs de
1 2 n

la base par rapport à laquelle la matrice A se transforme sous la forme de Jordan

1 n

Ax1 =λ1 x1 + x2
Ax2 =λ1 x2
Ax3 =λ3 x3
..
Axn =λn xn

On trouve
n
X
u1 =Au0 = a1 λ1 x1 + a1 x2 + a2 λ1 x2 + ai λ i x i

..
i=3

n
X
k
uk =A u0 = a1 λk1 x1 + a1 kλ1k−1 x2 + a2 λk1 x2 + a1 λk1 x1
i=3
n
X
=(a1 x1 + a2 x2 )λk1 + a1 kλ1k−1 x2 + ai λki xi
i=3

d'où (y, uk+1 )

Sk = = λ1 [1 + O(1/k)].
(y, uk )
En utilisant la méthode basée sur le rapport de déterminants précédent, on aurait obtenu
Sk′ = λ21 + O[(λ3 /λ1 )k ].

5.2.1 Accélération de la convergence

La convergence de la méthode puissance peut être extrêmement lente s'il existe une
valeur propre de module voisin de |λ |. On peut alors accélérer sa convergence en utilisant
un procédé d'accélération de la convergence fort utile en analyse numérique : le procédé
1

△ d'Aitken ( qui est étudé en détail dans un chapitre suivant).

2
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

Le procédé △ d'Aitken consiste à transformer une suite de nombre (S ) en une nou-

2

velle suite (T ) dénie par

k
k

(△Sk )2
Tk = Sk − , k = 0, 2, ...,
△2 Sk
avec
△Sk = Sk+1 − Sk
2
△ Sk =△Sk+1 − △Sk .
Dans certains cas, (T ) converge plus vite que (S ) vers la même limite S, c'est-à-dire
k k

Tk − S
lim = 0.
k7→∞ Sk − S

D'après ce qui précède, la suite (S ) générée par la méthode de la puissance dans le cas 1
est de la forme
k

S = S + (b + ε )(S
k − S) k k−1

avec lim = 0. Alors

k7→∞

△Sk =(b + εk+1 )(Sk − S) − (b + εk )(Sk−1 − S)

=(b − εk+1 − 1)(b + εk )(Sk−1 − S)
2
△ Sk =[(b + εk+2 − 1)(b + εk+1 ) − (b + εk+1 − 1)](b + εk )(Sk−1 − S),
d'où, si b ̸= 1 ( ce qui est vérié dans le cas 1 de la méthode de la puissance), lim T = S
et (T ) converge vers S plus vite que S ).
k7→∞ k

Dans le cas où |λ | > |λ | et où les vecteurs propres forment une base, nous avons vu que
k k

la suite (S ) obtenue par la méthode de la puissance vériait S = λ + O[(λ /λ ) ]. Si

1 2
k

on applique le procédé △ d'Aitken à cette suite alors on obtient T = λ + O[(λ /λ ) ].

k k 1 2 1
2 k

Donc si |λ | > |λ3|, il y a accélération de la convergence et le gain dépend de la proximité

k 1 3 1

des trois valeurs de plus grands module.

2

Il existe une généralisation du procédé △ d'Aitken qui s'appelle l'ε-algorithme et qui est
2

due à D. Shanks et P. Wynn. Cet algorithme permet dans le cas de la méthode de la

puissance, d'accélérer d'avantage la convergence.
5.2.2 Variantes de la méthode de la puissance
Il existe diverses variantes de la méthode de la puissances (proposées par de nombreux
auteurs)
[Link].1 La méthode des traces On rappelle que la trace d'une matrice est la
somme de ses éléments diagonaux et que cette trace est aussi égale à la somme de ses
valeurs propres. On a, si |λ | > |λ | pour i = 2, ..., n,
1 i

tra(Am ) =λm
1 + · · · + λn
m

" n
#1/m  
m 1/m
X
m 1 m
[tra(A )] =λ1 1 + (λi /λ1 ) = λ1 + O (λ2 /λ1 )
i=2
m
 53

d'où
lim [tr(Am )]1/m = λ1 .
m7→∞

Du npoint de vu pratique, il est plus facile de calculer la suite A , A , A , A ,..., par 2 4 8 16

élévation au carré des matrices succéssivement obtenues, d'où la méthode

h i1/2k 1
 
2k 2k
tr(A ) = λ1 + O k (λ2 /λ1 ) .
2

[Link].2 La méthode du produit scalaire Elle consiste à générer les deux suites
de vecteurs
donnév =
uk+1 =Ak , u0A v , k+1
T
k v0 , donné
Si les valeurs propres λ de A sont distinctes, on peut écrire
i

et
n
X n
X
u0 = ai x i v0 = bi yi
i=1 i=1

où x et y sont les vecteurs propres de A et A respectivement. On a

i i
T

(uk , vk ) =(Ak u0 , (AT )k v0 ) = (v0 , A2k u0 )

n n
!
X X
= bi yi , ai λ2k
i xi .
i=1 i=1

Or (y , x ) = 0 si λ ̸= λ , d'où
i j i j

n
X
(uk , vk ) = ai bi λ2k
i (xi , yi ).
i=1

De même
(uk , vk−1 ) =(Ak u0 , (AT )k−1 v0 ) = (v0 , A2k−1 u0 )
n
X
= ai bi λi2k−1 (xi , yi ).
i=1

D'où, si |λ | > |λ |. 1 2
(uk , vk )
= λ1 + O (λ2 /λ1 )2k .
 
(uk , vk−1 )
Cette methode est partiellement bien adaptée au cas des matrices sumétriques car si
u = v alors, pour tout k , k = v .
0 0 k k
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

5.2.3 Recherche des autres valeurs propres

Supposons, qu'à l'aide de l'une des méthodes précédentes, nous ayons déterminé la
valeur propre de plus grand module λ et éventuellement le vecteur propre correspondant
x . Nous voulons maintenant déterminer la valeur propre de module immédiatement in-
1

férieur à |λ |.
1

Il existe trois façon de procéder :

1

 Modication de la méthode itérative an qu'elle converge vers λ au lieu de conver-

ger vers λ . c'est la λ-diérence.
2

 Modication de la matrice an qu'elle possède la valeur propre 0 au lieu de λ et

1

que ses autres valeurs propres λ , ...., λ restent inchangées. C'est la défraction.
1

 Utilisations du procédé △ d'Aitken.

2 n
2

1. La λ-diérence :
On a X n
uk+1 − λ1 uk = a2 (λ2 − λ1 )λk2 x2 + ai (λi − λ1 )λki xi
i=3

d'où (y, uk+2 − λ1 uk+1 )

= λ2 + O (λ3 /λ2 )k .
 
(y, uk+1 − λ1 uk )
Ce rapport converge donc vers λ si |λ | > |λ | ≥ · · · ≥ |λ |, si (u − λ u , x ) ̸= 0
et (ynx ) ̸= 0.
2 2 3 n 1 1 0 2
2

2. La déation :
Soit x le vecteur propre de A correspondant à λ et y le vecteur propre de A T

correspondant à λ . On sait que (x , y ) = 0 si λ ̸= λ . Considérons la matrice A

i i i

dénie par
i i j i j 1

x1 y1T
A1 = A − λ1 .
(x1 y1 )
Cette matrice admet les valeurs propres 0, λ , ..., λ . En eet
2 n

x1 y1T x1
A1 x1 =Ax1 − λ1
(x1 , y1 )
(x1 , y1 )
=λ1 x1 − λ1 x1 =0
(x1 , y1 )

ce qui prouve que 0 est valeur propre de A avec x comme vecteur propre. De
même, pour i = 2, ..., n, on a
1 1

x1 y1T xi
A1 xi =Axi − λ1
(x1 , y1 )

=λi xi − λ1 x1
(xi , y1 )
(x1 , y1 )
= λ i xi , si λ 1 ̸= λi .
 55

L'inconvénient de cette méthode de déaction est qu'elle nécessitye la connaissance

de y . Si A est symétrique alors x = y et cette méthode est avantageuse. cepen-
dant, en pratique, cette méthode de déaction est numériquement instable. Son
1 1 1

usage n'est donc pas recommandé. On utilise à la place une méthode de déaction
faisant appel à une transformation orthogonale (ou tout simplement semblable)
comme on va le voir maintenant.
Soit H une matrice non singulière quelconque telle que
Hx = λe , avec λ ̸= 0.
1 1 (5.1)
Ax = λ x peut s'écrire HA(H H)x = λ Hx en multipliant à gauche par H .
−1

D'ou
1 1 1 1 1 1

−1
HAH e = λ e . 1 1 1

Par conséquent la première colonne de HAH est égale à λ e −1

1 1
 
λ1 bT1
A1 = HAH −1 =
0 B1
où B est une matrice carrée de dimension n − 1 qui admet λ , ..., λ comme va-
leurs propres et b ∈ R . On peut donc calculer λ en appliquant la méthode de la
1 2 n
n

puissance à la matrice B .
1 2
1

Pour obtenir x , le vecteur propre de A correspondant à λ . Notons y le vecteur

propre A correspondant à λ (il s'obtient aussi grâce à la méthode de la puissance
2 2 2

à la matrice B ) et écrivons
1 2
1  
α
y2 =
z2
où α ∈ R et z est un vecteur de B correspondant à la valeur propre λ . On a
2 1 2
    
λ1 bT1 α α
= λ2
0 B1 z2 z2
d'où
(λ1 − λ2 )α + (b1 , z2 ) = 0.
On a
HAH −1 y2 =λ2 y2
AH −1 y2 =λ2 H −1 y2
et par conséquent x 2 = H −1 y2 .
Construction de la matrice H
Elle s'eectue de la façon suivante (méthode de Householder [5] ). On peut dé-
montrer qu'il existe une matrice H orthogonale élémentaire qui satisfait 5.1 et qui
s'écrire sous la forme
H = I − 2uuT , u ∈ Rn , (u, u) = 1.
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

On a
∥Hx1 ∥2 = ∥x1 ∥2
d'où
|λ| = ∥x1 ∥2 .
Posons µ = 2(u, x ). On a
1

Hx1 = (I − 2uuT )x1 = λe1

et donc
µu = x1 − λe1 .
Multiplions scalairement 2x , on trouve µ = 2λ(λ − x ), où x est la première
2

composante du vecteur x . Ces équations déterminent donc u ∈ R , λ et µ ∈ R

1 1,1 1,1
n

tels que Hx = λe avec H = I − 2uu .

1
T

Pour des raisons de stabilité numérique on choisit λ = −(signe de x )∥x ∥ . On

1 1

peut aussi éliminer la racine dans le calcul de µ. Pour plus de détails, voir [5].
1,1 1 2

3. Le procédé △ d'Aitken :
2

Dans ce qui précède, nous avions utilisé le procédé △ d'Aitken pour accélérer la
2

convergence de la suite (S ) obtenue par la méthode de la puissance. Cette suite

(S ) était donnée par
k
k
Sk = yk+1 /yk , k = 0, 1, ...
avec y = (y, u ) et u = Au . Nous allons maintenant appliquer le procédé △ 2

d'Aitken à la suite (y ). Nous obtenons la suite

k k k+1 k
k

zk = yk − (△yk )2 /△2 yk , k = 0, 1, ...

On peut montrer que, sous certaine hypothèses, la suite

Zk = zk+1 /zk , k = 0, 1, ...

converge vers λ .2

5.2.4 La méthodes de la puissance inverse

Dans certaines applications, on a besoin de calculer la valeur propres de plus petit
module d'une matrice. Puique les valeurs propres de A sont les inverses de celles de A,
−1

cela revient à calculer la valeur propre de plus grand module de A . Pour ce faite, on −1

utilise la méthode de la puissance inverse qui consiste à faire les itérations suivantes
uk+1 = A−1 uk , k = 0, 1, ...

c'est-à-dire que u est solution du système linéaire Au = u . Du point de vue pra-

tique, avant de commencer les itérations, on eectue une décomposition de A en un produit
k+1 k+1 k

A = LU avec L triangulmaire inférieure à diagonale unité et U triangulaire supérieure

 57

(voir méthode de Gauss). On obtient ensuite u en résolvant deux systèmes linéaires

triangulaires
k+1

uk+1 = U −1 (L−1 uk ).
On peut aussi adapter la méthode de la puissance au calcul de la valeur propre la plus
voisine d'un nombre donné σ. Ilsut, pour cela, d'appliquer la méthode
uk+1 = (A − σI)−1 uk = U −1 (L−1 uk )

avec A − σI = LU .
Cette méthode est dicilement applicable aux matrices de grande dimension.
5.3 Méthodes de décomposition
Nous allons maintenant étudier les méthodes de décomposistion qui fournissent simul-
tanément toutes les valeurs propres d'une matrice A.
Ces méthodes sont basées sur l'idée suivante :
La décomposistion de la matrice A en un produit de deux matrices A = BC . Si B est
inversible, on peut alors écrire que C = B A d'où
−1

CB = B −1 AB

ce qui montre que les matrices A = BC et CB sont semblables et donc qu'elles ont les
mêmes valeur propres. On peut de nouveau décomposer la matrice CB en un produit de
deux matrices à l'aide de la même décomposition et ainsi de suite.
On pose donc A = A et, dans une méthode de décomposistion, on construit la suite de
matrice (A ) par
0
k

A0 =B0 C0
A1 =C0 B0 = B1 C1
A2 =C1 B1 = B2 C2
..
Ak =Ck−1 Bk−1 = Bk Ck
Ak+1 =Ck Bk = Bk+1 Ck+1
..
où toutes les matrices B sont de même type (par exemple, toutes triangulaires inférieurs
à diagonale unité ou prthogonales) et toutes les matrices C aussi (par exemple, toutes
k

triangulaires supérieures).
k

Toutes les matrices A sont semblables à A. On va donc choisir la décomposistion de façon

que
i

1. (A ) converge vers une matrice que nous appellerons A ,

k ∞

2. les valeurs propres de A soient faciles à calculer.

∞
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

Propriétés 5.1 Posons −1

Pk = Pk−1 Pk ,
d'où −1 −1
lim Bk = lim Pk−1 Pk = P ∞ P∞ = I.
k→∞ k→∞

Puisque A k = Bk Ck , on aura
A∞ = lim Bk Ck = lim Ck .
k→∞ k→∞

La décomposistion devra par conséquent être choisie de sorte que les valeurs propres
des matrices C soient simples à calculer : c'est le cas lorsque les matrices C sont des
matrices triangulaires supérieures, les valeurs propres sont alors sur la diagonale. Or nous
k k

connaissons justement deux décompositions de A en un produit BC où C est triangulaire

supérieure (voir [5]) :
 La décomposistion LU qui est utilisée dans la méthode de Gauss pour la résolution
des systèmes d'équations linéaires (L est triangulaire inférieure avec des 1 sur la
diagonale et U est triangulaire sypérieure).
 la décomposistion QU de Householder où Q es une matrice orthogonale et U est
triangulaire supérieure.
Les méthodes que nous allons maintenant étudier seront basées sur ce genre de dé-
composistions. La décomposition LU conduit à l'algorithme LR dû à Rutishouser (1955);
la décomposistion QU conduit à l'algorithme QR, trouvé indépendemment par Francis et
Kublanovskaya en 1961, et à la méthode de Jacobi (1846) dans le cas où A est symétrique.
Remarque 5.1 1. On a vu que A = B A B . La connaissance des matrices C −1

n'est donc pas nécéssaire pour la mise en oeuvre d'une méthode de décomposition.
k+1 k k k k

On verra que la méthode de Jacobi est construite à partir de cette remarque et à

partir de la propriété de symétrie de A.
2. Si l'on pose Q = C · · · C C , on a
k k 1 0

−1
Pk Qk =Pk−1 Bk Ck Qk−1 = Pk−1 Ak Qk−1 = Pk−1 (Pk−1 APk−1 )Qk−1
=APk−1 Qk−1
=A2 Pk−2 Qk−2 = · · · = Ak+1 .

On voit donc que P Q est la décomposition de A .

k k
k+1

3. Soit x le vecteur propre de A correspondant à la valeur propre λ . On a Ax = λ x .

Posons z = P x où P est une matrice inversible quelconque. Alors z est un
i i i i i
−1

vecteur propre de la matrice P AP correspondant à la même valeur propre λ . En

i i i
−1

eet x = P z d'où AP z = λ P z . Multiplions par P , on obtient P AP z =

i
−1 −1

λ P P z = λ z . Donc une fois calculée les valeurs propres de A , on calculera

i i i i i i
−1

ses vecteurs propres z puis les vecteurs propres x de la matrice A par x = P z .

i i i i ∞
i i i ∞ i
 59

5.3.1 La méthode de Jacobi

La méthode de Jacobi (voir [?]) s'applique à une matrice A réelle symétrique. Comme
dans la remarque précédente, on peut mettre en oeuvre une méthode de décomposition en
n'utilisant que les matrices B et leurs inverses. Dans la méthode de Jacobi, on prendra,
comme matrices B , les matrices
k
k

...
1
 
 
 

 1 


..
cos θ · · · · · · · · · − sin θ
.. ← p
 
 
 
.. 1
... ..
 
 

.. ..
 
 
Bk = 
 
1

 
sin θ · · · · · · · · · cos θ ← q
 
 

...
1
 
 
 
 
1
↑ ↑
p q
Il est facile de vérier que B est une matrice orthogonale, c'est-à-dire que B = B et −1 T

que toutes les matrices A obtenues sont symétriques (et semblable à A).
k k k

Lorsque l'on passe de A à A , seuls sont modiés les éléments des lignes et des colonnes
k

p et q . Notons a les éléments de A . Les nouveaux éléments sont donnés par

k k+1
(k)
ij k

(k+1) (k) (k) (k) (k) (k+1)

aip =aip cos θ + aiq sin θ = api cos θ + aqi sin θ = api , i ̸= p, q
(k+1) (k) (k) (k) (k) (k+1)
aiq = − aip sin θ + aiq cos θ = −api sin θ + aqi cos θ = aqi , i ̸= p, q
a(k+1)
pp =a(k) 2 (k) (k) 2
pp cos θ + 2apq sin θ cos θ + aqq sin θ (5.2)
a(k+1)
qq =a(k) 2 (k) (k) 2
pp sin θ − 2apq sin θ cos θ + aqq cos θ

a(k+1)
pq =(a(k) (k) (k) 2 2 (k+1)
qq − app ) sin θ cos θ + 2apq (cos θ − sin θ) = +aqp .
Ces relations montre que, si A est symétrique, alors A l'est aussi. L'angle θ est choisi
an d'annuler a et a , d'où
k k+1
(k+1) (k+1)
pq qp

(k) (k)
(a − a ) sin θ cos θ = a (cos − sin θ).
pp qq (5.3)
(k)
pq
2 2

Divisons par cos θ, il vient

2

(a(k) (k) (k) 2

pp − aqq )tgθ = apq (1 − tg θ).

Or, puisque tg2θ(1 − tg θ) = 2tgθ, on a, si a

2 (k)
pp ̸= aqq
(k)
,
(k)
2apq
tg 2θ = (k) (k)
,
app − aqq
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

où |θ| < π/4. Si a = a on prendra θ = π/4.

(k) (k)

Si l'on pose µ = cotg 2θ et t = tg θ, alors, puisq ue µ = (1 − t )/2t, il est susant,

pp qq
2

pour obtenir un angle de rotation θ tel que |θ| < π/4, de choisir la plus petite racine de
l'équation t + 2µt − 1 = 0, c'est-à-dire
2

signe(µ)
t = tg θ = p ,
|µ| + 1 + µ2

et l'on a
c = cos θ = √
1
1 + t2
, et s = sin θ = cos θ tan θ = ct.

En utilisant les relations (5.2) et (5.3 et le fait qur a = 0, on obtient les nouvelles
(k+1)

formules suivantes qui semblent être numériquement plus stables

pq

(k+1) (k) (k) (k) (k+1)

aip =aip + s(aiq − τ aip ) = api , i ̸= p, q
(k+1) (k) (k) (k) (k+1)
aiq =aiq − s(aip − τ aiq ) = aqi , i ̸= p, q
a(k+1)
pp =a(k) (k)
pp + tapq

a(k+1)
qq =a(k) (k)
qq − tapq ,

avec τ = s/(1 + c) = tg(θ/2).

Les indices p et q varient à chaque itération. Donc les termes a et a qui avaient été (k) (k)

annulés lors de la (k−1)-ième itération ne seront plus obligatoirement nuls après la k-ièmes
pq qp

itération. Cependant en itérant le procédé pour tous les couples (p, q) avec p ̸= q on arrive
peu à peu à annuler tous les éléments extradiagonaux. La matrice A convergera donc
vers une matrice diagonale semblable à A et les valeurs propres seront sur sa diagonale.
k

5.3.2 La méthode Greenstadt

Cette méthode proposée par Greenstadt en 1955, est une généralisation de la méthode
de Jacobi au cas où la matrice A n'est pas symétrique. Dans la méthode de Jacobi, l'angle
θ intervenant dans les matrices B était choisi à chaque itération de sorte que a = 0.
(k)

Comme les matrices étaient symétriques, on avait également a = 0. si la matrice

k pq
(k+1)

A n'est pas symétrique, on va annuler a , q > p ; c'est idée de base de la méthode de

qp
(k+1)

Greenstadt. On aboutira ainsi à une matrice A triangulaire supérieure. Comme matrices

qp
∞
 61

Bk on prendra les matrices orthogonales

...
1
 
 
 

 1 


..
c · · · · · · · · · −s
.. ← p
 
 
 
.. 1
... ..
 
 

.. ..
 
 
Bk = 
 
1

 
s ··· ··· ··· c ← q
 
 

...
1
 
 
 
 
1
↑ ↑
p q
avec c + s = 1 et c ̸= 0. De nouveau, seuls sont modiés les éléments des lignes et des
2 2

colonnes p et q. On obtient
(k+1) (k) (k)
aip =caip + saiq , i ̸= p, q
(k+1) (k) (k)
aiq = − saip + caiq , qquadi ̸= p, q
(k+1) (k) (k)
apj =capj + saqj , qquadj ̸= p, q
(k+1) (k) (k)
aqj = − sapj + caqj , qquadj ̸= p, q
a(k+1)
pp =c2 a(k) (k)
pp + cs(apq + a(k)
qp ) + s2 a(k)
qq

a(k+1)
qq =c2 a(k) (k)
qq − cs(apq + a(k)
qp ) + s2 a(k)
pp

a(k+1)
pq =c2 a(k) (k)
pq + cs(aqq − a(k)
pp ) − s2 a(k)
qp

a(k+1)
qp =c2 a(k) (k)
qp + cs(aqq − a(k)
pp ) − s2 a(k)
pq .

On veut annuler a , d'où

(k+1)
qp

c2 a(k) (k) (k) 2 (k)

qp + cs(aqq − app ) − s apq = 0
c2 + s2 = 1.
Il nous faut résoudre ce système de deux équations non linéaires à deux inconnues. Posons
x = s/c et 2△ = a − a . On trouve que a x − 2△ x − a = 0 d'où
(k) (k) (k) 2 (k)
qp qq pp pq qp qp

c =(1 + x2 )−1/2
s =(1 + x2 )−1/2 x
où x est l'une des deux racines
(k) (k)
△qp ± (△2qp + aqp apq )1/2
x= (k)
apq
à condition que a (k)
pq ̸= 0 .
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

5.3.3 L'algoritme LR
Cet algorithme, formulé par Heinz Rutishauser en 1955 (voir [11]), s'applique à des
matrices quelconques. Il est basé su la méthode de Gauss pour résoudre les systèmes
d'équations linéaires (voir [5]), c'est-à-dire, pour des matrices B , on predra des matrices
triangulaires inférieures avec des 1 sur la diagonale ( elles sont notées L dans la termi-
k

nologie habituelle) et pour C des matrices triangulaires supérieures, notées maintenant

k

R . On aura donc A = L R , puis A = R L . Les éléments r et l des matrices

k
(k) (k)

R et L sont donnés par les formules

k k k k k+1 k k ij ij
k k

(k) (k) Pi−1 (k) (k)


rij = aij − p=1 lip rpj , i = 1, ..., j 





(k)

ljj = 1

j = 1, ..., n.


(k) 1 h (k) Pj−1 (k) (k)
i 

lij = aij − p=1 lip rpj , i = j + 1, ..., n 

(k) 
rjj


Si lim A = A existe, alors c'est une matrice triangulaire supérieure. En eet

nous avons vu précédemment que
k→∞ k ∞

A∞ = lim Rk .
k→∞

Les valeurs propres de A se trouveront donc sur la diagonale de A . ∞

Remarque 5.2 Si A est symétrique, on utilisera, à la place de la méthode de Gauss,

la méthode de Cholesky (voir [5]). Les matrices L et R sont alors trasposées l'une de
l'autre; on aura donc
k k

Ak =Lk LTk
Ak+1 =LTk Lk .

Si la méthode converge, on voit que l'on aboutira, par raison de symétrie, à une matrice
diagonale.
Le théorème ci-dessous regroupe quelques propriétés de convergence de l'algorithme
LR.

Théorème 5.1 1. Si A est diagonalisable, si |λ | > |λ | > · · · > |λ | et si tous les

mineurs fondamentauxde la matrice des vecteurs propres de A et de son inverse
1 2 n

sont diérents de zéro, alors l'algorithme LR converge et de plus on a

 
λ1

...
 λ2 
A∞ =  .
 
 
λn

2. Si A est hermitienne dénie positive, alors l'algorithme LR converge.

 63

5.3.4 l'Agorithme QR
Cet algorithme a été trouvé indépendamment par J. F. G. Francis [6] et Vera Kubla-
novskaya en 1961 [9]. Il consite à prendre pour matrices B des matrices orthogonales
(notées Q dans la terminologie habituelle) et pour matrices C des matrices triangulaires
k

supérieures (notées R ).
k k

Pour ce faire, on utilise la décomposition de Householder pour résoudre les sytèmes d'équa-
k

tions linéaires. Nous avons déjà utilisé cette décomposition pour mettre une matrice sous
forme de Hessenberg supérieur [5]. la décomposition de Householder est eectuée à l'aide
de n−1 produits par des matrices orthogonales élémentaires (au lieu de n−2 pour obtenir
la forme Hessenberg).
Soit à décomposer A en A = Q R . On pose
k k k k

A(1) = Ak
et l'on eectue une succession de transformations orthogonales
A(r+1) = Hr A(r) , r = 1, ..., n − 1,
où A (r)
,r > 1, est de la forme
 (2) (2) (2) (2)

a11 ··· ··· a1,r−1 a1r ··· a1n
(2) (2) (2) (2) 
a22 · · ·
...
a2,r−1
.. a2r
.. ··· a2n 
..


 
 
 
(2) (2) (2)
A(r) = ar−1,r−1 ar−1,r · · · ar−1,n 
 
(2) (2) 

.. ..


 arr · · · arn  
 
 
(2) (2)
anr · · · ann
c'est à dire, en l'écrvant sous une forme partitionnée en 4 blocs
!
(r) (r)
A11 A12
A(r) = (r)
0 A22

On passe de A à A à l'aide d'une matrice orthohonale élémeantaire H

(r) (r+1)
= I−
2u u , avec u = 0 pour i = 1, ..., r − 1.
r
(r) (r)T (r)

Désigons par ue ∈ R le vecteur de composantes u , ..., u et notons He

i
(r) n−r+1 (r) (r)
= I−
e . On a
r n r
(r) (r)T
2e
u u
 
I 0
Hr =
0 H
er
! !
  (r) (r) (r) (r)
I 0 A11 A12 A11 A12
A(r+1) = (r) = .
0 H
er 0 A22 0 e r A(r)
H 22

Les éléments de la première colonne de He A ne sont pas tous nuls. La matrice He sera
(r)

donc choisie à annuler tous les éléments de première collone de He A sauf un seul, le
r 22 r
(r)
r 22
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

premier, ce qui correspond à annuler tous les éléments de la rième colonne de A , à partir (r)

de la ligne (r + 1) jusqu'à la ligne n.

La matrice A ainsi obtenue est triangilaire supérieure et l'on a
n)

A(n) = Hn−1 ...H1 Ak

d'où la décomosition QR cherchée en prenant R = A et Q = H ...H puisque (n)

Q =Q =Q .
k k n−1 1
−1 T

Comme dans la méthode de Householder pour triangulariser une matrice ( voir [?]),
k k k

la décomposition A = Q R se résume aux régles suivantes.

On calcule, pour r = 1, ..., n − 1,
k k k

!1/2
signe de a
n
X (r)
a(r+1)
rr =−( (r)
rr ) |air |2
i=r
(r)
v =a(r+1)
rr (a(r+1)
rr − a(r)
rr )
(r)
wi =0, i = 1, ..., r − 1
wr(r) =a(r) (r+1)
rr − arr
(r) (r)
wi =air , i = r + 1, ..., n
(k)
air =0, i = r + 1, ..., n

puis
n
X
(r) (r) (r)
βj = wi aij , j = r + 1, ..., n
i=r
(r) (r)
γj =βj /v (r) , j = r + 1, ..., n
(r+1) (r) (r) (r)
aij =aij − γj wi , i = r, ..., n; j = r + 1, ..., n,

où H = I − 2w w /v .
r
(r) (r)T (r)

Nous regoupons quelques propriétés inmportantes de l'algorithme QR dans le théoreme

suivante (voir [4] pour les détails) :
Théorème 5.2 1. Si A est une matrice quelconque, non singulière, et si ses valeurs
propres sont de modules tous diérents, alors l'algorithme Qr converge et la matrice
A est triangulaire supérieure.
∞

2. La forme de Hessenberg supérieure est invariante par l'algotithme QR.

3. Si on applique l'algorithme QR à une matrice de Hessenberg supérieure A telle que
a ̸= 0 pour i = 1, ..., n − 1, alors l'algorithme QR converge et la matrice A est
triangulaire par blocs et semblable à A. Les valeurs propres de chaque bloc diagonal
i+1,i ∞

sont toutes égales entre elles en module.

 65

5.4 Pseudo-codes et Codes Matlab

Dans cette section, nous nous limitons aux pseudo-codes de la méthode de Jacobi et
de l'algorithme QR. Deux (2) pseudo-codes de la méthode de jacobi sont donnés, l'un avec
le chois cyclique pour les indices de l'élément à annuler, l'autre avec le choix classique de
Jacobi.
5.4.1 L'algorithme de Jacobi, choix cyclique
Dans le pseudo-code suivant, on considère la méthode de Jacobi pour le calcul des
valeurs et vecteurs propres d'une matrice A réelle symétrique. Le choix de p et q est
eectué dans l'ordre cyclique. chaque cycle complet d'itérations est constitué de n(n−1)/2
rotations. Les éléments de A sont dénis, et en sortie ses valeurs propres se trouvent sur
sa diagonale. En rentrée, on doit donner la valeur de toll (qui correspond à ε) pour le
test d'arrêt N (A ) ≤ ε et nombre maximum k de cycles complets de rotations. Dans
chacun de ces cycles ( qui correspondent aux instructions qui se trouvent dans la double
k max

boucle en p et en q), la rotation n'est eectuée que si la valeur de l'élément que l'on
droit annuler satisfait à la condition |a | > tolla où tolla = toll/(n − n). Les valeurs
p(k) 2

propres se trouvent en sortie dans la matrice V . On fournit égélement le nombre k de

pq

cycle complets eectués ainsi que le nombre n de rotatiuons. r

[A, V, k, nr ] = Jacobi-c(A, n, toll, kmax )

Initialisation for i = p + 1, ..., q − 1
V = In tmp = api
nr = 0 api = api + s(aiq − τ api )
k=0 p aiq = aiq − s(tmp + τ aiq )
tolla = toll/(n2 − n) end for i
test = 2 p<q |apq |2 for i = q + 1, ..., n
P

while test > toll and k < kmax do tmp = api

k =k+1 api = api + s(aqi − τ api )
for p = 1, ..., n − 1 aqi = aqi − s(tmp + τ aqi )
for q = p + 1, ..., n end for i
if |apq | > tolla and test > toll then for i = 1, ..., n
µ = (apq − ap qq /(2apq ) tmp = vip
t = 1/(|µ| + 1 + µ2 ) vip = vip + s(viq − τ aip )
if µ < 0 then viq = viq − s(tmp + τ viq )
t = −t end for i
end if p for j = 1, ..., n − 1
c = 1/ 1 − t2 for i = j + 1, ..., n
s = tc aij = aji
τ = s/(1 + c) end for i
app = app + tapq end for j
aqq = aqq − tapq nr = nr + 1
test = test − 2a2 pq end if
apq = 0 end for q
for i = 1, ..., p − 1 end for p
tmp = aip end while
aip = aip + s(aiq − τ aip )
aiq = aiq − s(tmp + τ aiq )
end for i

5.4.2 L'algorithme de Jacobi, choix classique

Dans ce pseudo-code, le choix de p et q est celui utilisé par Jacobi. Donc, à chaque
itération, p et q correspondent à l'élément extradiagonal de plus grand module, c'est à
dire |a | = max |a |.
(k) (k)

Les éléments de A sont détruits, et en sortie ses valeurs propres se trouvent sur sa diago-
pq i<j ij

nale. En entrée, on doit donner la valeur de toll(qui correspond à ε ) pour le test d'arrêt
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

N (Ak ) ≤ ε et le nombre maximum n de rotations. Les vecteurs propres se trouvent

en sortie dans la matrice V . On fournit également le nombre n de rotation.
rmax
r

[A, V, k, nr ] = Jacobi-j(A, n, toll, n

rmax )

Initialisation for i = q + 1, ..., n

V = In tmp = api
nr = 0 P api = api + s(aqi − τ api )
test = 2 p<q |apq |2 aqi = aqi − s(tmp + τ aqi )
while test > toll and nr < nrmax do end for i
Déterminer p et q tels que |apq | = maxi<j |aij |. for i = 1, ..., n
µ = (app − aqq
p)/(2apq ) tmp = vip
t = 1/(|µ| + 1 + µ2 ) vip = vip + s(viq − τ vip )
if µ < 0 then viq = viq − s(tmp + τ viq )
t = −t end for i
end if√ for j = 1, ..., n − 1
c = 1/ 1 + t2 for i = j + 1, ..., n
s = tc aij = aji
τ = s/(1 + c) end for i
app = app + tapq end for j
aqq = aq − tapq nr = nr + 1
test = test − 2a2pq end while
apq = 0
fori = 1, ..., p − 1
tmp = aip
aip = aip + s(aiq − τ aip )
aiq = aiq − s(tmp + τ aiq )
end for i
for i = p + 1, ..., q − 1
tmp = api
api = api + s(aiq − τ api )
aiq = aiq − s(tmp + τ aiq )
end for i

5.4.3 L'algorithme QR
Ce pseudo-code eectue l'algorithme QR (sans shift) sur une matrice A déja mise sous
forme de Hessenberg supérieure, voir [?].
La matrice A est détruite et, à la n des itérations, elle est remplacée par une matrice
triangulaire surpérieure semblable à la matrice de départ. Par conséquent, ses valeurs
propres se trouvent sur sa diagonale.
Puisque l'algorithme QR doit converger vers une matrice triangulaire supérieure, on
arrête les itérations lorsque les termes de la première sous-diagonale sont susamment
 67

petits ou lorsque le nombre maximum d'itérationd k est atteint. En sortie, k donne le

nombre d'itérations eectuées.
max

[A, k] = QR(A, n, toll, kmax )

Initialisations wr = arr − tmp
k=0 wr+1 = ar+1,r
test = maxi=1,...,n−1 |ai+1,i | arr = tmp
while test ≥ toll and
k < kmax do ar+1,r = 0
k =k+1 for j = r + 1, ..., n
Q = In γj = (wr arj + wr+1 ar+1,j )/v
for r = 1, ..., n − 1 arj = arj + γj wr+1
if arr < 0 then ar+1,j = ar+1,j − γj wr+1
s = −1 end for j
else Qr = In − wwT /v
s=1 Q = QQr
end if q end for r
tmp = −s a2rr + a2r+1,r A = AQ
v = tmp(tmp − arr ) test = maxi=1,...,n−1 |ai+1,i |
for i = 1, ..., n end while
wi = 0
end for i

5.5 Exercices
Exercice 5.1 Soient L et M deux matrices carrées triangulaires inférieures d'ordre N .
Montrer que le produit LM est également triangulaire inférieur. Déduire de ce résultat la
propriété : si A est une matrice régulière d'ordre N qui possède une décomposition LU
(avec l = 1, i = 1, 2, ..., N ), alors cette décomposition est unique.
ii

Exercice 5.2 (Calcul de valeurs propres,10 points) Étant donné un vecteur propre
réel non nul u et une matrice réelle symétrique d'ordre n A, on sait que le quotient de
Rayleigh Q(u) s'écrit
uT Au
Q(u) =
uT u
1. Si v un vecteur de A, associé à la valeur propre λ , calculer Q(v ) (1 ≤ i ≤ n).
2. Soient θ (i = 1, 2, ..., n) les coecients du développement de u sur la base des v .
i i i

(a) Explimer Q(u) en fonction des θ et des λ .

i i

(b) On suppose que les valeurs propres λ ont été numérotées en ordre décrois-
i i

sant :λ > λ > · · · > λ . Montrer que, quel que soit u ̸= 0, on a :

i
1 2 n

λ1 ≥ Q(u) ≥ λn .
 CHAPITRE 5. CALCUL DE VALEURS PROPRES

3. On donne la matrice B et le vecteur u :

   
2 1 cosθ
B= ; u=
1 2 sinθ

Former Q(u) = Q(θ) et en déduire les valeurs propres de B.

4. Montrer comment on peut utiliser le quotient de Rayleigh pour accélérer le calcul
de la valeur propre dominante par la méthode de la puissance n-ième.
Bibliographie

[1] A. Angot. Complements de Mathématique à l'usage des ingénieurs. Masson, PAris,

1972.
[2] M. Mabramowitz, I.A. Stegun. Handbook of mathematical functions. Dven, New
Yprk, 1972.
[3] E. Belrizky. Outils mathématiques à l'usage des sceintiques et Ingénieur, ch. 9, EDP,
Grenoble Sciences, 2007.
[4] Claude Brezinski, Michela Redivo-Zaglia. Méthodes numériques itératives. Mathé-
matiques a l'Universite. Cours et exercices. Ellipses. 2006.
[5] C. Brezinski, M. Redivo-Zaglia, Méthodes Numeriques Directes de l'Algèbre Matri-
cielle, Ellipses, Paris, 2004.
[6] J.G. F. Francis, The QR transformation : a unitary analogue to the LR transforma-
tion, I, Comput. J., 4(1961/1962)265-271; II, Comput. J., 4(1961/1962)332-345.
[7] C. G. J. Jacobi, Uber eine neue Auösungart der bei der Methode der kleinsten
Quadrate vorkommenden lineären Gleichungen, Astron. Nachr., 22(1845)297-306.
[8] Jean-Philippe GRIVET. Methodes numériques appliquées pour le scientique et l'in-
génieur. Collection Grenoble Sciences. EDP Sciences, 2009.
[9] V. N. Kublanovskaya. On some algorithms for the solution of the complete eigenvalue
problem (en russe), Z. Vych. Mat. Fiz., 1(1961)555-570; traduction anglaise dans
USSR Comput. Math. Math. Phys., 1(1962)637-657.
[10] Jean-Louis Merrien. Analyse numériques avec Matlab. Exercices & Problèmes, DU-
NOD, Paris, 2007.
[11] H. Rutishauser, Une méthode pour la détermination des valeurs propres d'une ma-
trice, C. R. Acad. Sci. Paris, 240(1955), 34-36.
[12] R. S. Varga. Matrix Itérative Analysis , 2nd ed. , Springer -Verlag, Berlin, 2000.
[13] Ahmed ZIDNA. Contribution a la modélisation des carreaux troues. PHD, Université
de Metz, 1990.