FST-Settat

Université Hassan Premier

Méthodes numériques
pour l'Ingénieur
1ème Année Génie Industriel/Logistique

Pr. Bouchaïb Radi

Faculté des Sciences et Techniques,
Route de Casablanca, Settat.
Maroc.
TABLE DES MATIÈRES 2

Table des matières

1 Interpolation 4
1.1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Techniques de l'interpolation polynomiale . . . . . . . . . . . . 5
1.3 Interpolation dans la base de Lagrange . . . . . . . . . . . . . 5
1.3.1 Base de Lagrange . . . . . . . . . . . . . . . . . . . . . 5
1.3.2 Erreur d'interpolation polynomiale . . . . . . . . . . . 6
1.4 Interpolation dans la base de Newton . . . . . . . . . . . . . . 6
1.4.1 Erreur d'interpolation de Newton . . . . . . . . . . . . 7

2 Norme matricielle et conditionnement 8

2.1 Norme vectorielle . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Norme matricielle . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Conditionnement d'une matrice . . . . . . . . . . . . . . . . . 10

3 Résolution des systèmes linéaires 12

3.1 Méthodes directes . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.1.1 Méthode de Gauss . . . . . . . . . . . . . . . . . . . . 12
3.1.2 Méthode de Cholesky . . . . . . . . . . . . . . . . . . . 14
3.2 Méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2.1 La méthode de Jacobi . . . . . . . . . . . . . . . . . . 15
3.2.2 La méthode de Gauss-Seidel . . . . . . . . . . . . . . . 16

4 Résolution des équations non linéaires 17

4.1 Position du problème . . . . . . . . . . . . . . . . . . . . . . . 17
4.2 Approximation d'une racine séparée . . . . . . . . . . . . . . . 18
4.2.1 Méthode de dichotomie (ou de bissection) . . . . . . . 18
4.2.2 Méthode du point xe (approximations successives) . . 18
4.2.3 Premier critère de convergence . . . . . . . . . . . . . . 19
TABLE DES MATIÈRES 3

4.2.4 Critères d'arrêt des itérations . . . . . . . . . . . . . . 19

4.2.5 Deuxième critère de convergence (critère local) . . . . . 20
4.2.6 Méthode de Newton (ou méthode des tangentes) . . . . 20

5 Intégration numérique 22
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.2 La méthode des trapèzes . . . . . . . . . . . . . . . . . . . . . 23
5.3 La méthode de Simpson . . . . . . . . . . . . . . . . . . . . . 24
5.4 Methode de Gauss-Legendre . . . . . . . . . . . . . . . . . . . 25
5.4.1 Position du problème . . . . . . . . . . . . . . . . . . . 25
5.4.2 Polynômes de Legendre . . . . . . . . . . . . . . . . . . 25

6 Equations diérentielles 26
6.1 Problème de Cauchy . . . . . . . . . . . . . . . . . . . . . . . 26
6.2 Résolution numérique . . . . . . . . . . . . . . . . . . . . . . . 26
6.2.1 Méthode d'Euler . . . . . . . . . . . . . . . . . . . . . 27
6.2.2 Méthode de Runge-Kutta à pas unique . . . . . . . . . 28
6.2.3 Méthode de Prédicteur-Correcteur . . . . . . . . . . . . 30

7 Méthodes numériques de calcul des valeurs propres et vec-

teurs propres 31
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
7.2 Calcul direct du det (A − λI) . . . . . . . . . . . . . . . . . . 32
7.3 Méthode de Krylov . . . . . . . . . . . . . . . . . . . . . . . . 32
7.4 Méthode de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . 34
7.5 Méthode de Givens-Householder . . . . . . . . . . . . . . . . . 36
 4

Chapitre 1
Interpolation

1.1 Généralités

Nous considérons les polynômes comme outils pour résoudre des pro-
blèmes d'interpolation. Les algorithmes considérés seront donc uniquement
des algorithmes d'évaluation de polynômes ou d'expression de polynômes
dans diérentes bases.
La situation type d'un problème d'interpolation est la suivante :
Une fonction f est connue en un certain nombre d'abscisses x1 , x2 ,..., xn et
uniquement en ces points. f peut avoir une expression analytique compliquée,
ou bien les f (xi ) proviennent de mesures physiques.
Étant donné a appartenant à l'intervalle [min xi , max xi ], il s'agit d'évaluer
f (a). (Lorsque a 6∈ [min xi , max xi ], il s'agit d'un problème d'extrapolation).
Une fonction d'interpolation est une fonction p vériant

p(xi ) = f (xi ) i = 1, · · · ,n

on dit aussi interpolant" ?.

La résolution du problème d'interpolation se scinde généralement en deux
problèmes :

Problème 1 : Détermination d'une fonction d'interpolation (dans le but de

remplacer le calcul de f (a) par p(a)).

Problème 2 : Évaluation de l'erreur d'interpolation (lorsque l'on remplace

f (a) par p(a) on commet une erreur e(a) = f (a) − p(a) que l'on doit être en
Pr. B. RADI 5

mesure d'estimer).

1.2 Techniques de l'interpolation polynomiale

Cette partie est exclusivement consacrée au Problème 1 :

 Détermination de p ∈ Pn interpolant les (n+1) couples (xi ,yi ) (vériant

i 6= j ⇒ xi 6= yj ) ?.
 Détermination de p(a).

p est entièrement déterminé par la donnée des (n + 1) couples ; déterminer p

revient donc à calculer ses coecients dans une base de Pn . Les propriétés de
l'algorithme permettant de déterminer les coecients seront donc liées aux
propriétés de la base choisie. Trois bases sont généralement utilisées :

 La base canonique,

 La base de Lagrange,

 La base de Newton.

1.3 Interpolation dans la base de Lagrange

1.3.1 Base de Lagrange

La base de Lagrange est adaptée à un problème d'interpolation particulier
et dépend donc de (n + 1) abscisses distincts xi (i = 0, · · · ,n).
Posons

n
(x − x0 )(x − x1 )...(x − xi−1 )(x − xi+1 )...(x − xn ) Y x − xj
li (x) = =
(xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn ) xi − xj
j=0
j 6= i

Chaque polynôme li est de degré n et on a :


1 si i=j
li (xj ) = δij =
0 si i 6= j

Ces (n + 1) polynômes sont donc linéairement indépendants entre eux et

forment une base de Pn .
Méthodes Numériques : 2011/2012 6

On vérie immédiatement que le polynôme :

n
X
p(x) = yi li (x)
i=0

réalise les conditions p(xi ) = yi (i = 0, · · · ,n). La représentation de Lagrange

est donc particulièrement bien adaptée au problème d'interpolation puisque
la donnée des couples (xi ,yi ) est en fait la donnée des composantes yi du
polynôme d'interpolation sur la base de Lagrange associée aux xi .

Intérêt de la méthode de Lagrange

Elle a un intérêt plus théorique que pratique car :

 coût élevé en nombre d'opérations.

 formulation peu aisée : une modication d'une abscisse oblige à recalcu-

ler une nouvelle base de Lagrange. De même, si on ajoute un nouveau
couple (xn+1 ,yn+1 ) il faut tout recommencer.

 par contre on peut changer les yi sans recalculer la base.

1.3.2 Erreur d'interpolation polynomiale

Soit p un polynôme de degré inférieur ou égal à n interpolant une fonction
f en (n + 1) abscisses xi . La quantité e(x) = f (x) − p(x) est appelée erreur
d'interpolation.
Soit I l'intervalle
I = [min{(xi ),x}, max{(xi ),x}]
i i

Théorème 1 Si f est (n + 1) fois continûment diérentiable sur I , alors :

f (n+1) (ξ)
∃ξ ∈ I : e(x) = (x − x0 )(x − x1 ) · · · (x − xn )
(n + 1)!

1.4 Interpolation dans la base de Newton

pk le polynôme
Soit d'interpolation de la fonction f aux points x0 , x1 ,...,
xk pour k = 0,1, · · · ,n.
Problème : Construire Pn par récurrence.
 On a p0 (x) = f (x0 ) (polynôme de degré 0).
Méthodes Numériques : 2011/2012 7

 Pour k ≥ 1, le polynôme pk − pk−1 est de degré k et s'annule aux points

x0 , x1 , · · · , xk−1 . Il est donc de la forme :
pk (x) − pk−1 (x) = f [x0 ,x1 , · · · ,xk ](x − x0 )(x − x1 ) · · · (x − xk−1 )
où f [x0 ,x1 , · · · ,xk ] désigne le coecient de xk dans pk (x).
En ajoutant ces relations pour k = 1,2, · · · ,n, il en résulte que
n
X
pn (x) = f (x0 ) + f [x0 ,x1 , · · · ,xk ](x − x0 )(x − x1 ) · · · (x − xk−1 ) (1.1)
k=1

Remarque 1 L'ensemble des n + 1 polynômes

{1; (x − x0 ); (x − x0 )(x − x1 ); · · · ; (x − x0 )(x − x1 ) · · · (x − xn−1 )}
est une base de l'espace Pn , appelée base de Newton.
La formule (1.1) ramène le calcul par récurrence de pn à celui des coecients
f [x0 ,x1 , · · · ,xk ].

Dénition 1 La quantité f [x0 ,x1 ,...,xk ] est appelée diérence divisée d'ordre
k de f aux points x0 ,x1 , · · · ,xk .
Lemme 1
f [xi ] = f (xi ), i = 0, · · · ,k
f [x1 , · · · ,xk ] − f [x0 ,x1 , · · · ,xk−1 ]
∀k ≥ 1 f [x0 ,x1 , · · · ,xk ] =
xk − x0

1.4.1 Erreur d'interpolation de Newton

L'erreur d'interpolation de Newton est donnée par le corollaire suivant :

Corollaire 1 L'erreur d'interpolation de Newton est donnée par :

en (x) = f (x) − pn (x) = Π(x).f [x0 ,x1 , · · · ,xn ,x]
Si f ∈ C n+1 [a,b] alors ∃ξ ∈ [a,b] tel que :
f (n+1) (ξ)
f [x0 ,x1 , · · · ,xn ,x] =
(n + 1)!
Remarque 2 D'après sa dénition, la diérence divisée d'ordre k de f aux
points x0 , x1 ,..., xn , f [x0 ,x1 , · · · ,xk ] est indépendante de l'ordre des points
xi .
 8

Chapitre 2
Norme matricielle et
conditionnement

2.1 Norme vectorielle

Il existe plusieurs manières de dénir la norme pour un vecteur. La plus

connue de toutes est certainement la norme euclidienne:
v
u n
uX
k(x1 ,x2 , · · · ,xn )k = t x2i , xi ∈ R.
i=1

Dénition 2 Une application k.k : V → R (ici, V est un espace vectoriel

réel) est dite une norme vectorielle si elle respecte les propriétés suivantes :
1. k v k≥ 0, ∀v ∈ V et k v k= 0 ⇔ v = 0.
2. k λv k=| λ |k v k, ∀λ ∈ R, v ∈ V .
3. k v + u k≤k v k + k u k, ∀v , u ∈ V .
On peut facilement vérier que la dénition de la norme euclidienne satisfait
les propriétés ci-dessus. Comme nous l'avons mentionné, il existe plusieurs
normes pour un même espace vectoriel. Si on se limite aux normes dénies
n
sur R , on trouve les normes p, dites usuelles, d'un vecteur x. Celles-ci se
calculent de la manière suivante :

n
!1/p
X
k x kp = |xpi | .
i=1
Pr. B. RADI 9

Lorsque p = 2, on retrouve la dénition de la norme euclidienne. Il y a

également la norme dite innie qui se dénit comme suit :

kxk∞ = max |xi |

1≤i≤n

Cette dernière norme est l'une des plus utilisée avec la norme 2 et la norme
1. Il existe également des normes un peu plus exotiques (exemple norme
elliptique), mais nous allons travailler avec les normes que nous venons de
donner.

2.2 Norme matricielle

Comme pour les vecteurs de Rn , les matrices possèdent une norme.

Dénition 3 Soient A, B ∈ Mn×m , ensemble des matrices réelles de format

n × m. Alors, l'application k . k : Mn×m −→ R+ est une norme matricielle
si elle satisfait les propriétés suivantes :
1. k A k≥ 0 et k A k= 0 ⇔ A = 0.
2. k βA k=| β | . k A k et ce, pour tout scalaire β .
3. k A + B k≤k A k + k B k.
4. k AB k≤k A k . k B k.
Regardons maintenant quelques normes matricielles. Celles qui sont les
plus utilisées, encore une fois nommées usuelles, découlent des normes vecto-
rielles décrites plus haut.
Dénition 4 Soit A une matrice n×m, étant donnée une norme vectorielle,
on appelle norme matricielle subordonnée, la norme matricielle dénie par :
k Ax k
k A k= max = max k Ax k .
x∈Rn , x6=0 k x k kxk=1

En réalité, cette relation signie que la norme d'une matrice correspond au

maximum de la norme du vecteur résultant de la multiplication de la matrice
A par tous les vecteurs de norme égale à 1. Sous cette forme, il est dicile
de calculer la norme d'une matrice A. Nous allons énoncer quelques résultats
qui vont rendre la tâche plus facile.

Théorème 2 Si A = {aij } ∈ Mn×m , la norme innie satisfait la formule

m
X
k A k∞ = max | aij | .
1≤i≤n
j=1
Pr. B. RADI 10

On a un résultat similaire pour la norme 1.

Théorème 3 Si A = {aij } ∈ Mn×m , la norme 1 satisfait la formule
n
X
k A k1 = max | aij | .
1≤j≤m
i=1

Il est cependant plus dicile d'arriver à une formule pareille pour la

norme 2. Mais avant de regarder cette formule, dénissons ce qu'est le rayon
spectral d'une matrice.

Dénition 5 Le rayon spectral d'une matrice A carrée, noté ρ(A), est déni
de la manière suivante :
ρ(A) = max | λi |
i=1,··· ,n

où les λi sont les valeurs propres de A. Il est à noter que les valeurs propres
peuvent être complexes et | λi | représente le module d'un nombre complexe,
en général (qui est un nombre réel positif ou nul).
Théorème 4 Soit A la matrice de format n × m, alors :
p
k A k2 = ρ(AT A).
Corollaire 2 Si A est une matrice carrée symétrique de format n×n, alors :
k A k2 = ρ(A).

2.3 Conditionnement d'une matrice

Dénition 6 On appelle conditionnement d'une matrice A carrée non-singulière

relatif à une certaine norme, le nombre
Kp (A) = kAkp .kA−1 kp
Si le problème est bien conditionné, le nombre de conditionnement sera
proche de 1. On peut constater cela en calculant le nombre de condition-
nement de la matrice identité, qui est très bien conditionnée. Il est dit mal
conditionné si K(A) est grand (et mal posé si K(A) est inni).

Théorème 5 Si A est une matrice carrée inversible alors

• Perturbation de la matrice
k δx k k δA k
(A + δA)(x + δx) = b ⇒ ≤ K(A) .
k x + δx k kAk
Pr. B. RADI 11

• Perturbation du second membre

k δx k k δb k
A(x + δx) = b + δb ⇒ ≤ K(A) .
kxk kbk

Théorème 6 s
1. K(A) = K(A−1 ).
ρ(AT A)
2. K2 (A) = où (B) = min |λ|.
(AT A) λ∈Sp (B)

ρ(A)
3. Si A est hermitienne, alors K2 (A) = .
(A)
4. Si A est symétrique dénie positive alors K2 (A) = λmax
λmin
= ρ(A)ρ(A−1 ).
5. ∀α ∈ R∗ , K(αA) = K(A).
6. Si U est unitaire, alors K2 (U ) = 1.
7. Si U est unitaire, K2 (U.A) = K2 (A).
 12

Chapitre 3
Résolution des systèmes linéaires

3.1 Méthodes directes

Le problème de résolution des systèmes linéaires se rencontre très fré-

quemment dans beaucoup de problèmes ayant une origine physique. En eet,
pour résoudre des problèmes plus compliqués d'analyse numérique (équations
aux dérivées partielles, par exemple), on est souvent ramené à résoudre un
système linéaire de n équations à n inconnues de la forme :


 a11 x1 + a12 x2 + ... + a1n xn = b1
a21 x1 + a22 x2 + ... + a2n xn = b2




. = .

(S)

 . = .
. = .




an1 x1 + an2 x2 + ... + ann xn = bn


Avec aij (1 ≤ i ≤ n ; 1 ≤ j ≤ n) et bi (1 ≤ i ≤ n) sont des nombres réels

donnés et xi (1 ≤ i ≤ n) sont des inconnues.
En utilisant la notation matricielle, le système linéaire (S) s'écrit d'une
manière plus simple sous la forme Ax = b, avec A la matrice carrée formée
par les éléments aij , x et b les matrices colonnes formées respectivement par
les éléments xi et bi .

3.1.1 Méthode de Gauss

On dénit une matrice A = (aij )1≤i≤n,1≤j≤n+1 , associée au système (S) tel
que ai,n+1 = bi .
Méthodes Numériques : 2011/2012 13

Pour obtenir le système triangulaire équivalent, il sut d'appliquer (n−1)

itérations de Gauss sur les lignes de la matrice A.
On dénit l'itération k de Gauss de la façon suivante :
1ier
Cas
Si akk 6= 0, on laisse les lignes 1 à k inchangées et on construit une nouvelle
matrice à partir de A de la façon suivante:

aik
aij ← aij − akj , i = k + 1, · · · ,n; j = k, · · · ,n + 1. (3.1)
akk

2ème Cas
Si akk = 0, on échange la ligne k avec une ligne l (k + 1 ≤ l ≤ n) telle que
alk 6= 0. Les lignes 1 à k − 1 restent inchangées et on transforme le reste de
la matrice A comme précédemment c-à-d :

akj ←→ alj , j = k, · · · ,n + 1 (3.2)

et
aik
aij ← aij − akj , i = k + 1, · · · ,n; j = k, · · · ,n + 1 (3.3)
akk
3ème Cas
Si akk = 0, et alk = 0 ∀k ≤ l ≤ n, on arrête le processus (le système est
singulier).

Ainsi, après l'itération k, la matrice A est de la forme :

 
a11 a12 . . . . . . ... . . . a1n a1,n+1
 0 a22 . . . . . . ... . . . a2n a2,n+1 
 
 .. . . . . .
.
.
.

 . . . . . 
 .
..

 . . akk
A= . ak,k+1 . . . akn ak,n+1 

 . . 
 .. 0 ak+1,k+1 .
.

 
 .. .
.
.
.
.
.

 . . . . 
0 ... ... 0 an,k+1 . . . ann an,n+1
avecaii 6= 0, ∀ 1 ≤ i ≤ k .
Enn, après (n − 1) itérations la matrice A devient :
Méthodes Numériques : 2011/2012 14

 
a11 . . . . . . . . . a1n a1,n+1
.. . .
 0

. . . 
. . 
 . .. .. . .

 ..
A= . . .
.
.
.


 . .. .. . .
 .. . .

. . . . 
0 ... ... 0 an,n an,n+1
et le système linéaire associé est alors triangulaire supérieur et admet la même
solution que le système initial.

Remarque 3 Pour l'itération k, l'élément akk s'appelle le pivot de l'itéra-

tion.
Ce pivot ne doit pas être nul et pour des raisons de précision des calculs on
a intérêt à choisir un pivot dont la valeur est la plus éloignée possible de 0,
c-à-d la plus grande possible en valeur absolue.
Ainsi, à chaque itération k, on commencera par chercher la ligne l tel que :
|alk | = maxk≤i≤n |aik |.
Si |alk | = 0, on arrête le processus (système singulier), sinon on échange la
ligne k avec la ligne l et on transforme la matrice comme précédemment.

3.1.2 Méthode de Cholesky

On va maintenant étudier la méthode de Cholesky, qui est une méthode
directe adaptée au cas où A est symétrique dénie positive. On rappelle
T
qu'une matrice A ∈ Mn (R) de coecients (ai,j )1≤i,j≤n est symétrique si A =
T
A, où A désigne la transposée de la matrice A, dénie positive si (Ax,x) > 0
pour tout x ∈ Rn tel que x 6= 0. On rappelle aussi que si A est symétrique
dénie positive, alors A est inversible.

Théorème 7 Si A est une symétrique dénie positive, alors elle peut être
décomposée en A = L.LT .

Décomposition de la matrice A
Pn
Développons l'équation A = L.LT , ai,j = k=1 li,k lj,k , or li,k =0 si i < k.
On a donc
min(i,j)
X
ai,j = li,k lk,j , 1 ≤ i, j ≤ n.
k=1
Pr. B. RADI 15

ième
Pour la partie triangulaire supérieure de A, la r ligne s'écrit :
r
X
ar,j = lr,k lj,k j = r, · · · ,n
k=1
Soit
r−1
X
ar,j = lr,k lj,k + lr,r lj,r j = r, · · · ,n.
k=1
Ce qui donne l'algorithme suivant :
" r−1
# 21 
X 
(lr,k )2

lr,r = ar,r −





k=1 

r = 1, · · · ,n
r−1
X 

lik ljj

ar,j − 



k=1 
lj,r = lr,r
j = r + 1, · · · ,n 

3.2 Méthodes itératives

3.2.1 La méthode de Jacobi

On remarque que si les éléments diagonaux de A sont non nuls, le système
linéaire Ax = b est équivalent à :
n
!
1 X
xi = bi − aij xj , i = 1, · · · ,n. (3.4)
aii j=1,j6=i

Pour une donnée initiale x(0) choisie, on calcule x(k+1) par

n
!
(k+1) 1 X (k)
xi = bi − aij xj , i = 1, · · · ,n. (3.5)
aii j=1,j6=i

Cela permet d'identifer le splitting suivant pour A: A = D − E − F

D : diagonale de A
E : triangulaire inférieure avec des 0 sur la diagonale
F : triangulaire supérieure avec des 0 sur la diagonale
La matrice d'itération de la méthode de Jacobi est donnée par :

J = D−1 (E + F ) = I − D−1 A. (3.6)

L'algorithme de Jacobi nécessite le stockage des deux vecteurs x(k) et x(k+1) .

Méthodes Numériques : 2011/2012 16

3.2.2 La méthode de Gauss-Seidel

Pour cette méthode, on a :

i−1 n
!
(k+1) 1 X (k+1)
X (k)
xi = bi − aij xj − aij xj , i = 1, · · · ,n. (3.7)
aii j=1 j=i+1

Il s'écrit aussi
x(k+1) = (D − E)−1 (b + F x(k) ) (3.8)

Dans ce cas, le splitting de A est

M = D − E, N = F (3.9)

et la matrice d'itération associée est

BGS = (D − E)−1 F. (3.10)

(k)
L'algorithme de Gauss-Seidel ne nécessite qu'un vecteur de stockage, x
(k+1)
étant remplacé par x au cours de l'itération. Il est en général plus rapide
que l'algorithme de Jacobi, donc préférable.
 17

Chapitre 4
Résolution des équations non
linéaires

4.1 Position du problème

Soit F : R → R une fonction dénie sur une partie D ⊂ R. On considère

l'équation F (x) = 0 dont on cherche à déterminer les racines éventuelles.
A l'exception de quelques cas particuliers, on ne sait pas résoudre analy-
tiquement ce problème. C'est la raison pour laquelle, très souvent, on fait
appel aux méthodes numériques pour la résolution approchée de l'équation
F (x) = 0.
Les méthodes numériques proposées sont nombreuses, chacune a ses avan-
tages et ses inconvénients. (on se limitera dans ce chapitre à présenter les
méthodes les plus usuelles).
Notons qu'avant toute tentative de résolution numérique d'un problème donné,
on doit d'abord s'assurer de l'existence et de l'unicité d'une solution.
Ainsi, pour résoudre numériquement l'équation F (x) = 0, nous procéderons
en deux étapes :
La première étape consistera à montrer l'existence des solutions réelles et à
séparer chaque solution à approximer dans un intervalle [a,b] ⊂ D.
La deuxième étape consistera à choisir une méthode numérique pour approxi-
mer les racines séparées.
Pr. B. RADI 18

4.2 Approximation d'une racine séparée

Dans ce paragraphe, on note x une racine de l'équation F (x) = 0 et on

suppose que x est séparée dans l'intervalle [a,b].
Nous allons proposer quelques méthodes usuelles permettant de trouver une
valeur approchée de x et pour lesquelles, nous étudierons la convergence et
la vitesse de convergence.
Notons que l'approximation de x à  près consiste à trouver une valeur ap-
prochée x
e telle que |x − x
e| ≤ . ( un nombre positif petit; plus  est petit
plus la précision est meilleure).

4.2.1 Méthode de dichotomie (ou de bissection)

On suppose x ∈ [a,b] racine simple séparée (F (a)F (b) < 0).
La méthode de dichotomie consiste à construire à partir de l'intervalle [a,b],
une suite d'intervalles In = [an ,bn ] contenant x tel que la largeur de In soit
la moitié de celle de In−1 .
Principe de la méthode :(construction de 3 suites (an ),(bn )et(xn ))
On pose a0 = a; b0 = b, x0 = 21 (a0 + b0 ).
et pour n ≥ 1 :
Si F (an−1 ).F (xn−1 ) < 0, on pose an = an−1 , bn = xn−1 , xn = 12 (an + bn )
Si F (an−1 ).F (xn−1 ) > 0, on pose an = xn−1 , bn = bn−1 , xn = 12 (an + bn )
Si F (an−1 ).F (xn−1 ) = 0 alors x = xn−1 : la solution exacte est trouvée, on
arrête le processus.
De telle sorte qu'à chaque itération, on choisit le demi-intervalle contenant
la racine x.
b−a
Après n itération, la largeur de l'intervalle contenant x devient
2n
.

Théorème 8 La suite (xn ) dénie par l'algorithme de dichotomie converge

vers la racine x.
L'erreur commise en approximant x par xn est majorée par : |xn − x| ≤ b−a
2n+1
.

4.2.2 Méthode du point xe (approximations succes-

sives)
Dénition 7 On appelle point xe d'une fonction f , toute solution de l'équa-
tion f (x) = x.
Méthodes Numériques : 2011/2012 19

Théorème 9 (corollaire du TVI)

Si f est continue sur [a,b] et si f ([a,b]) ⊂ [a,b], alors il existe entre a et b au
moins un point xe de f (càd ∃c ∈]a,b[/f (c) = c)
On suppose x ∈ [a,b] F (x) = 0.
racine séparée de l'équation
La méthode du point xe consiste à associer à la fonction F une fonction f
tel que x est un point xe de f (c.à d. tel que F (x) = 0 ⇔ f (x) = x).

Ainsi, au lieu de chercher x comme racine de F (x), on le cherchera comme

point xe de f (x), en utilisant la propriété suivante :
Si la suite récurrente (xn ) dénie par : x0 donné; xn = f (xn−1 ),n ≥ 1 converge
alors sa limite est forcément un point xe de la fonction f.

4.2.3 Premier critère de convergence

Théorème 10 Soit x ∈ [a,b] racine séparée de F (x) et soit f (x) une fonc-
tion continue et dérivable sur [a,b] telle que f (x) = x
S'il existe un réel k, 0 < k < 1 tel que :
(i) f ([a,b]) ⊂ [a,b],
(ii) ∀x ∈ [a,b], |f 0 (x)| ≤ k
Alors la suite récurrente (xn ) dénie par : x0 ∈ [a,b]; xn = f (xn−1 converge
vers x.
Une majoration de l'erreur commise en approximant x par xn est donnée
par :
|x − xn | ≤ k n |x − x0 |

4.2.4 Critères d'arrêt des itérations

La méthode proposée, consiste à déterminer x comme limite d'une suite
numérique convergente. Dans la pratique, on doit déterminer une valeur ap-
prochée de x avec une précision  donnée, en un nombre ni d'itérations (le
plus petit possible).
Puisque x est inconnue, on ne peut pas utiliser directement le test d'arrêt
|x − xn | ≤ .
Méthodes Numériques : 2011/2012 20

4.2.5 Deuxième critère de convergence (critère local)

Théorème 11 Soit x ∈ [a,b] racine séparée de F (x) et soit f (x) une fonc-
tion de classe C 1 telle que f (x) = x
1. Si |f 0 (x)| < 1, alors il existe un voisinage V de x tel que la suite
récurrente (xn ) associée à f converge vers x, ∀x0 ∈ V .
De plus, (xn ) converge d'une façon monotone si 0 ≤ f 0 (x) < 1 et (xn )
converge en oscillant de part et d'autre de x si −1 < f 0 (x) < 0.
V s'appelle le domaine d'attraction de x
2. Si |f 0 (x)| > 1, alors ∀x0 , la suite xn = f (xn−1 ) ne converge pas vers x.
3. Si |f 0 (x)| = 1, on ne peut pas conclure. Le comportement de la suite
(xn ) est instable (peut converger ou diverger) et dépend fortement des
erreurs d'arrondi.

4.2.6 Méthode de Newton (ou méthode des tangentes)

On suppose avoir isolé une racine x de l'équation F (x) = 0 dans un
0
intervalle [a,b] et on suppose que F (x) ne s'annule pas sur [a,b].
La méthode de Newton consiste à choisir comme fonction f associée à F, la
fonction dénie par :
F (x)
f (x) = x −
F 0 (x)
Il est évident que x est un point xe de f . Alors, que peut-on dire de la suite
récurrente (xn ) dénie par :

x0 ∈ [a,b], xn+1 = f (xn )?

Théorème 12 Soit x racine séparée de l'équation F (x) = 0 dans [a,b]. On

suppose que F est de classe C 2 et que F 0 ne s'annule pas sur [a,b].
Alors il existe un voisinage V de x tel que la suite (xn ) dénie par la méthode
de Newton converge vers x, ∀x0 ∈ V .
Théorème 13 Soit x une solution séparée de l'équation F (x) = 0 dans [a,b].
On suppose que F 0 ne s'annule pas sur [a,b] et que f 00 est continue sur [a,b];
Soit M = max |f 00 (x)|/2 sur [a,b].
Si x0 ∈ [a,b] et |x0 − x| < 1/M , alors la suite (xn ), dénie par la méthode de
Newton converge vers la solution x.
Méthodes Numériques : 2011/2012 21

Une majoration de l'erreur commise en approchant x par xn est donnée par :

1 n
|xn − x| ≤ (M |xn − x|2 ).
M
 22

Chapitre 5
Intégration numérique

5.1 Introduction

Étant donnée une fonction f intégrable sur [a,b], intégrer numériquement

Rb
f consiste à obtenir une valeur approchée de
a
f (x)dx au moyen d'une for-
mule de quadrature :

Z b n
X
f (x)dx = Ai f (xi ) + Rn (f )
a | {z }
|i=0 {z } erreur de quadrature
formule de quadrature

où xi ∈ [a,b] sont les noeuds de la formule et les Ai sont les poids de la

formule.
Une manière classique d'établir des formules de quadrature est la sui-
vante :
Choisissons (n + 1) abscissesxi ∈ [a,b], distinctes. Remplaçons f par p (po-
lynôme d'interpolation aux abscisses xi ). On a

n
X
f (x) = li (x)f (xi ) + En (f )
i=0

où li (x) est le polynôme de Lagrange et En (f ) est l'erreur d'interpolation.

On obtient une formule de quadrature par intégration ; l'erreur de qua-

Pr. B. RADI 23

drature s'obtient par utilisation de la formule de la moyenne.

Z b n
X
f (x)dx = Ai f (xi ) + Rn (f )
a i=0

Z b Z b
Ai = li (x)dx Rn (f ) = En (f )
a a

5.2 La méthode des trapèzes

Prenons xi = a+ih avec h = (b−a)/n. Sur chaque sous-intervalle [xi ,xi+1 ],

utilisons une formule de quadrature de Newton-Côtes à deux points. On a :
Z xi+1
h
f (x)dx = [f (xi ) + f (xi+1 )] + Ri (f )
xi 2

avec, en supposant f deux fois continuement diérentiable sur [a,b] :

Z 1
h3
Ri (f ) = f ”(ξi ) t(t − 1)dt ξi ∈ [xi ,xi+1 ]
2 0
h3
= − f ”(ξi )
12
On obtient donc :

Z b n−1 Z
X xi+1
f (x)dx = f (x)dx avec x0 = a
a i=0 xi
n−1 n−1
h X h3 X
= [f (xi ) + f (xi+1 )] − f ”(ξi )
2 i=0
12 i=0

d'où nalement la méthode des trapèzes :

Z b
f (a) f (b)
f (x)dx = h[ + f (a + h) + ... + f (a + (n − 1)h) + ] + R(f )
a 2 2
3
avec : R(f ) = − (b−a)
12n2
f ”(ξ) ξ ∈ [a,b]
On voit immédiatement que : lim R(f ) = 0.
n→∞
Pr. B. RADI 24

5.3 La méthode de Simpson

La méthode de Simpson consiste à utiliser une formule de Newton-Côtes

à trois points sur chaque sous-intervalle. Cette formule de Newton-Côtes fait
intervenir le milieu de chaque sous-intervalle. Il est plus simple de numéroter
de deux en deux les extrémités des sous-intervalles et de réserver un indice
impair aux milieux de ces intervalles. En d'autres termes on prend :
xi = a + ih pour i = 0, · · · ,n et h = (b − a)/h avec n pair
Z b Z x2 Z x4 Z b
f (x)dx = f (x)dx + f (x)dx + ... + f (x)dx
a a x2 xn−2

On a :
Z x2i+2
h
f (x)dx = [f (x2i ) + 4f (x2i+1 ) + f (x2i+2 )] + Ri (f )
x2i 3
en supposant f quatre fois continuement diérentiable sur [a,b] :
Z 2
h5 (4)
Ri (f ) = f (ξi ) t2 (t − 1)(t − 2)dt ξ ∈ [x2i ,x2i+2 ]
4! 0
5
h
= − f (4) (ξi )
90
On obtient donc :
n
Z b 2
−1 Z x2i+2
X
f (x)dx = f (x)dx
a i=0 x2i
n n
2
−1 2
−1
hX h5 X
= [f (x2i ) + 4f (x2i+1 ) + f (x2i+2 )] − f (4) (ξi )
3 i=0 90 i=0

d'où nalement la méthode de Simpson :

Z b
h
f (x)dx = [f (a) + 4f (a + h) + 2f (a + 2h) + ... + 2f (a + (n − 2)h)
a 3
+4f (a + (n − 1)h) + f (b)] + R(f )
5
avec : R(f ) = − (b−a)
180n4
f (4) (ξ) ξ ∈ [a,b].

On voit immédiatement que : lim R(f ) = 0.

n→∞
Pr. B. RADI 25

5.4 Methode de Gauss-Legendre

5.4.1 Position du problème

Soit [a,b] un intervalle de R et soit g une fonction réelle, on veut calculer

Z b X
g(x)dx = βk g(yk ) + R.
a

La méthode de Gauss-Legendre consiste à

Rb R1
1.
a
g(x)dx = K −1 f (t)dt avec x = 21 (b − a)t + 21 (a + b),
2. Choisir (xi )i=0,··· ,n et (αi )i=0,··· ,n tels que

Z 1 n
X
f (t)dt = αk f (xk ) + R
−1 k=0

soit exacte pour tout polynôme de P2n+1 .

5.4.2 Polynômes de Legendre

Dénition 8 Les polynômes de Legendre sont dénis par
1 dn
Ln (x) = ((−1 + x2 )n )
2n n! dxn

0 n 6= m
Proposition 1 1. In,m
R1
= −1 Ln (x)Lm (x)dx = 2
2n+1
n=m
2. {L0 ,L1 ,...,Ln } est une base de Pn .
3. Ln admet n racines réelles distinctes.
4. Soit ψ(u,x) = ∞ i=0 Li (x)u s'appelle la fonction génératrice des (Ln )
i
P

−1
ψ(u,x) = (1 − 2xu + u2 ) 2 .

5. Ln (x) = 2n−1
n
xLn−1 (x) − n−1
n
Ln−2 (x).
 26

Chapitre 6
Equations diérentielles

6.1 Problème de Cauchy

Soient f : [a,b]×R −→ R et η ∈ R. On s'intéresse au problème de Cauchy

suivant :
Existe-t-il une fonction y : [a,b] −→ R, dérivable sur [a,b] et vériant :

y(a) = η
(P)
y 0 (x) = f (x,y(x)) ∀x ∈ [a,b]

Dénition 9 On dit que la fonction f : [a,b] × R −→ R est Lipschitzienne

en y dans [a,b] × R s'il existe A > 0 tel que:
|f (x,y) − f (x,z)| ≤ A|y − z| ∀x ∈ [a,b], ∀y,z ∈ R

A est appelée constante de Lipschitz.

Théorème 14 Soit le problème de Cauchy (P ). Si f : [a,b] × R → R vérie
les hypothèses :
i) f continue
ii) |f (x,y) − f (x,z)| ≤ K|y − z| ∀x ∈ [a,b], ∀y,z ∈ R
Le problème (P ) admet une solution unique.

6.2 Résolution numérique

On discrétise l'intervalle [a,b] en une suite (xi )0≤i≤N

x0 = a < x1 < x2 < · · · < xN = b (xi = a + ih, h = b−a N
)
Méthodes Numériques : 2011/2012 27

On cherche N nombres y1 ,y2 , · · · ,yN où yi est une valeur approchée de y(xi ).

Puis on reliera ces points par interpolation pour dénir une fonction yh sur
[a,b].
Enn, on pourra estimer l'erreur de discrétisation en = yn − y(xn ) qui dé-
pendra de h.
Une méthode de calcul numérique fournissant des valeurs approchées yn de
y(xn ) telles que |en | ≤ K.hp est dite d'ordre p.
On va voir dans ce qui suit deux méthodes numériques permettant de calculer
la solution approchée yi+1 au point ti+1 en utilisant celle au point ti .

6.2.1 Méthode d'Euler


 yi+1 = yi + hf (ti ,yi )
y0 = y(a) = η
1≤i≤N


Convergence
Théorème 15 Si f vérie les hypothèses:
i) f ∈ C 1 ([a,b] × R)
ii) f est Lipschitzienne en y:
|f (x,y) − f (x,z)| ≤ K|y − z| ∀x ∈ [a,b], ∀y,z ∈ R, K > 0
La méthode d'Euler converge.
Plus précisément, si on pose M = max{|y 0 (t)|, t ∈ [a,b]}, on a la majoration:
eK(b−a) − 1 M
|en | = |yn − y(xn )| ≤ . .h
K 2
et
lim max |en | = 0
n→∞ n=1,··· ,N
.
Remarque 4 1. Le résultat du théorème précédent s'exprime sous la forme
|en | ≤ A.h, A > 0 une constante ; c'est à dire que la méthode d'Euler
est du premier ordre.
2. Une méthode d'ordre 1 ne converge pas assez vite pour donner des ré-
sultats pratiques intéressants. On va voir un autre exemple de méthodes
numériques d'ordre plus élevé.
Méthodes Numériques : 2011/2012 28

6.2.2 Méthode de Runge-Kutta à pas unique

Les méthodes de Runge-Kutta sont bien utilisées dans la pratique, car
elles présentent plusieurs avantages (facilité de programmation, stabilité de
(0)
la solution, modication simple du pas et la connaissance de y sut pour
intégrer l'équation diérentielle). Les inconvénients de cette méthode se ré-
sument au temps de calcul lent et à la diculté de l'estimation de l'erreur
locale.

Méthode Runge-Kutta d'ordre 2

Cette méthode est obtenue en prenant les diérences centrées au 1er
ordre :
ti1 = ti + θ1 h, 0 ≤ θ1 ≤ 1
yi1 = yi + hA10 f (ti ,yi )
yi+1 = yi + h[A20 f (ti ,yi ) + A21 f (ti1 ,yi1 )]
Les coecients θ1 , A10 , A20 , A21 sont à déterminer de sorte que la méthode
soit d'ordre 2 (voir erreur).
Notons fi = f (ti ,yi ),

yi+1 = yi + h [A20 fi + A21 f (ti + θ1 h,yi + hA10 fi )]

yi + h A20 fi + A21 fi + A21 θ1 hfti0 + A21 A10 fi hfyi
0
 
= + ...
h2 00 h3
y(ti+1 ) = y(ti ) + hy 0 (ti ) + y (ti ) + y (3) (ζi )
2 6
ei+1 = yi+1 − y(ti+1 )

On obtient donc le système


 A20 + A21 = 1
A21 .A10 = 12
A21 .θ1 = 12


système non linéaire de 3 équations à 4 inconnues. Généralement, on xe θ1 :

1
A21 =
2θ1
2θ1 − 1
A20 =
2θ1
A10 = θ1
Méthodes Numériques : 2011/2012 29

Méthode Runge-Kutta d'ordre 4

Dans la plus part des cas et pour des fonctions susamment régulières,
la méthode de Runge Kutta utilisée est celle d'ordre 4 (R.K.4) :

tij = ti + θj h
j−1
X
yij = yi + h Ajk f (tik ,yik )
k=0
A10 = θ1

avec j = 1,2,3,4 et yi3 = yi+1 .

On s'arrange à choisir les coecients θj de sorte que la méthode soit
d'ordre 4. On trouve alors :

1 1
θ 1 = , θ2 = , θ3 = θ 4 = 1
2 2

1
A10 =
2
1
A20 = 0, A21 =
2
A30 = 0, A31 = 0, A32 = 1
1 1 1 1
A40 = , A41 = , A42 = , A43 =
6 3 3 6
Donc, pour calculer yi+1 , on procède comme suit :

4t
yi+1 = yi + (f (yi ,ti ) + 2f (ŷi+1/2 ,ti+1/2 ) + 2f (ỹi+1/2 ,ti+1/2 ) + f (ŷi+1 ,ti+1 ))
6
4t
ŷi+1/2 = yi + f (yi ,ti )
2
4t
ỹi+1/2 = yi + f (ŷi+1/2 ,ti+1/2 )
2
ŷi+1 = yi + 4tf (ỹi+1/2 ,ti+1/2 )

Le calcul de yi+1 nécessite alors 4 évaluations de la fonction f , et par suite

pour les fonctions compliquées le temps de calcul devient important.
Méthodes Numériques : 2011/2012 30

6.2.3 Méthode de Prédicteur-Correcteur

La combinaison des deux méthodes (formules d'Adams ouvertes et for-
mules d'Adams fermées) constitue la méthode de prédicteur-correcteur. Il
(0)
sut de trouver une valeur estimée yi+1 proche de la valeur nale yi+1 . Le
(0)
prédicteur yi+1 est obtenu par une formule ouverte de même ordre. Cette
(0)
valeur yi+1 est introduite dans la méthode à formule fermée. Ceci permet
d'accélérer la convergence vers la solution. Ainsi, on prote des avantages
liés à la méthode d'Adams : temps de calcul réduit et erreur estimée de tron-
cature réduite considérablement.
Pour la méthode d'Adams au 4ème ordre, on a :

Prédicteur : formules d'Adams ouvertes au 4ème ordre :

(0) 4t
yi+1 = yi + (55fi − 59fi−1 + 37fi−2 − 9fi−3 ) (6.1)
24
(0)
au pas précédent, on a yi qui est la valeur du prédicteur (modica-
teur) :
(0) (0) 251 (0)
ŷi+1 = yi+1 − (yi − yi ) (6.2)
720

Correcteur : obtenu à partir de la formule d'Adams fermée (au 4ème ordre

par exemple) :

(
(r+1) (r)
yi+1 = yi + 4t24
.(9.fi+1 + 19fi − 5fi−1 + fi−2 )
(0) (0) (6.3)
fi+1 = f (ŷi+1 ,ti+1 )

Les valeurs y1 , y2 et y3
peuvent être calculées à partir des formules de Runge-
(0)
Kutta. Étant donné que y3 n'existe pas, la valeur y4 est calculée à partir des
(0)
équations (1.5) et (1.6) et y4 va initialiser les calculs d'itérations (formule
(0) (0)
(1.7)). À partir de y5 , on estime y5 par l'équation (5), puis ŷ5 par la
formule (1.6). Cette valeur est injectée dans la formule (1.7) pour amorcer
les itérations.
 31

Chapitre 7
Méthodes numériques de calcul
des valeurs propres et vecteurs
propres

7.1 Introduction

Dans les chapitres précédents, de nombreuses méthodes numériques sup-

posent la connaissance des valeurs propres ou du rayon spectral d'une ma-
trice. En outre, le calcul des valeurs propres et vecteurs propres est indispen-
sable dans toutes les branches de la science, en particulier pour chercher la
solution des systémes des équations diérentielles linéaires en théorie de sta-
bilité, pour les questions de convergence de processus itératifs en physique et
en chimie (mécanique, circuits, cinétique chimique, équations de Schrödinger)
et en automatique.

Théorème 16 (Cercles de Gerschogrin-Hadamard) Soit A une matrice

carrée d'ordre n, alors les valeurs propres de A appartiennent à la réunion
n
des n disques dénis par |z − aii | ≤ |aij |.
X

j=1

Corollaire 3 Une borne supérieure pour la plus grande valeur propre est
donc : ( " n
# " n
#)
X X
|λ| ≤ max |aij | , max |aij |
1≤i≤n 1≤j≤n
j=1 i=1
Pr. B. RADI 32

Exemple 1 Soit la matrice A dénie par

 
4 −1 1
A= 1 1 1 
−2 0 −6

or det (A − λI) = −λ3 − λ2 + 23λ − 26 = P (λ) donc P (−6) = 16 > 0,

P (0) = −26 > 0, P (2) = 8, P (4) = −14.
A admet trois valeurs propres réelles λ1 ∈ ]−6,0[, λ2 ∈ ]0,2[, et λ3 ∈ ]2,4[.
On pourra par exemple utiliser la méthode de Newton (voir chapitre ??) pour
approcher ces valeurs propres, on trouve λ1 ≈ −5.7685 ; λ2 ≈ 1.2992 et
λ3 ≈ 3.4694.

7.2 Calcul direct du det (A − λI)

Nous avons vu antérieurement, des méthodes permettant de calculer le dé-
terminant d'une matrice après triangularisation. Nous pouvons donc utiliser
ces méthodes pour calculer le polynôme caractéristique d'une matrice A. On
pourra aussi utiliser l'interpolation de Lagrange, puisque le polynôme carac-
téristique d'une matrice A d'ordre n est de degré n, on se donne (n + 1) va-
leurs distinctes µ1 ,µ2 , · · · ,µn+1 quelconques, et on calcule vi = det (A − µi I) ;
i = 1, · · · , (n + 1), alors le polynôme d'interpolation passant par les points
{(µi ,vi )}i=1,··· ,(n+1) sera identique au polynôme caractéristique de A. On peut
alors chercher ses racines par l'une des méthodes itératives pour l'approxi-
mation des valeurs propres de A.

7.3 Méthode de Krylov

La méthode consiste à calculer les coecients du polynôme caractéristique

dont on approche les racines à l'aide des méthodes classiques (dichotomie,
Newton, sécante ,...).
n
!
X
Soit P (λ) = (−1)n λn − ak λn−k le polynôme caractéristique de la
k=1
matrice A. D'après le théorème de Caley-Hamilton, on a P (A) = 0, donc
X n
n
A = ak An−k .
k=1
Pr. B. RADI 33

Prenons un vecteur x0 quelconque non nul, on a

n
X
An x 0 = ak An−k x0 . (7.1)
k=1

Notons a le vecteur de composantes (ai )ni=1 ; x1 = Ax0 , x2 = A2 x0 , ..., xn−1 =

An−1 x0 .
Notons B la matrice dont les colonnes sont les vecteurs xn−1 , xn−2 ,..., x0 .
L'équation (7.1) peut s'écrire :

An x0 = Ba (7.2)

Pour x0 donné arbitrairement, on cherche a solution de (7.2), on obtient ainsi,

si le système est inversible, les coecients du polynôme caractéristique.

Exemple 2 On considère la matrice

 
4 1 −1
A = −6 −1 2 
2 1 1

On a det (A − λI) = −λ3 + 4λ2 − 5λ + 2 = (λ − 1)2 (2 − λ).

On pose
     
1 8 2 −3 14 3 −7
x0 =  0  , A2 = −14 −3 6  et A3 = −26 −5 14  ;
0 4 2 1 6 3 1

le système (7.2) donne :

    
8 4 1 a1 14
−14 −6 0 a2  =  −26 
4 2 0 a3 6

−1 −3 
 

a1
0 14
 
4

 2 2 
⇒  a2  = 0 1
 7 
  −26  =  −5  .
a3
 2  6 2
1 0 −2
Pr. B. RADI 34

7.4 Méthode de Jacobi

La méthode de Jacobi est une méthode itérative applicable à une matrice

A symétrique. Elle consiste à faire opérer le groupe des rotations planes sur
A, c-à-d à multiplier A par des transformations orthogonales an de la mettre
sous forme diagonale, les éléments diagonaux étant les valeurs propres de la
matrice A.
Le principe de la méthode est le suivant : on considère la matrice H dont
les éléments sont égaux à ceux de la matrice identité sauf pour les quatre
valeurs suivantes : hpp = cos α, hpq = sin α, hqq = cos α et hqp = − sin α, avec
p < q. La matrice H est une matrice orthogonale c.-à-d. H T H = In .

On commence par calculer la matrice A1 = H1−1 AH = H T AH , en remar-

quant que seules les lignes et les colonnes p et q sont modiées, pour j = p
ou q, on a :
 (1) (1)

 apj = ajp = ajp cos α − ajq sin α

 (1) (1)
 aqj = ajq = ajp cos α + ajq sin α



(1)
app = app cos2 α + aqq sin2 α − 2 sin α cos α

(1)
= aqq cos2 α + app sin2 α + 2apq sin α cos α




 aqq
 (1)
 (1)
apq = aqp = qpq (cos2 α − sin2 α) + (app − aqq ) sin α cos α
(1) (1)
On peut donc choisir α de sorte que apq = aqp = 0, tel que

2apq π
tan 2α = ≡θ avec |α| ≤
app − aqq 4
ou encore
s  
1 1
cos α = 1+ √
2 1 + θ2
s  
1 1
sin α = sgn(θ) 1− √
2 1 + θ2
Si app = aqq , on choisira

1
cos α = √
2
sgn(apq )
sin α = √
2
Pr. B. RADI 35

On a alors
(a(1) 2 (1) 2 2 2 2
pp ) + (aqq ) = app + aqq + 2apq
(1)
et comme aii = aii pour i 6= p ou q
n
X n
X
(1)
(aii )2 = a2ii + 2a2pq
i=1 i=1

en passant de A à A1 la somme des carrés des éléments diagonaux augmente

2
de la quantité 2apq . En itérant ce processus, on obtient

Ak+1 = (H1 H2 . . . Hk )T A(H1 H2 . . . Hk ).

La suite des matrices Ak converge vers une matrice diagonale dont les élé-
ments diagonaux sont les valeurs propres de la matrice initiale A. La suite
des matrices Pk = H 1 H 2 . . . H k converge vers la matrice dont les colonnes
sont constituées des vecteurs propres. Au cours des itérations un terme peut
redevenir nul, mais on démontre que

X (k)
lim (aij )2 = 0.
k→∞
i6=j

On arrête l'itération quand

n
X (k)
(aii )2
i=1
1− n < .
X (k+1)
(aii )2
i=1

En pratique, on a le choix à chaque pas d'itération du couple (p,q). On dénit

diérentes stratégies :

• Dans la méthode de Jacobi classique, on choisit (p,q) tels que

(k)
|a(k)
pq | = sup |aij |.
i6=j

• Dans la méthode de Jacobi cyclique, on eectue un balayage systématique

en prenant pour (p,q) les couples (1,2), (1,3), ..., (1,n) puis (2,3),...,
(2,n), etc, jusqu'à (n − 1,n).
Pr. B. RADI 36

• Dans la méthode de Jacobi cyclique avec seuil, on eectue comme précé-

demment un balayage sur les éléments triangulaires supérieurs, chaque
élément aij étant pris comme élément à annuler apq , mais on ne retient
le couple (p,q) que si |aij | est supérieur à certain seuil qui peut être
réajusté à chaque itération.

Remarque 5 La méthode de Jacobi est stable, mais sa convergence est lente,

ce qui en fait une méthode très peu utilisée.

7.5 Méthode de Givens-Householder

La méthode Givens-Householder est la réunion de deux algorithmes. La

méthode de Householder met la matrice initiale A sous la forme tridiagonale
symétrique. La méthode de Givens calcule les valeurs propres d'une matrice
tridiagonale symétrique. On suppose que A soit mise sous la forme

 
b1 c1 0 ··· 0
.. .
c 1 b 2 . .
c2
 
. 
.. ..
 
B=
 0 c2 . . 0 
. . .. ..
 .. ..

. . cn−1 
0 ··· 0 cn−1 bn

et on note Bk la sous-matrice

 
b1 c1 0 ··· 0
.. .
c1 b2 . .
c2
 
. 
.. ..
 
Bk = 
 0 c2 . . 0 
. . .. ..
 .. ..

. . ck−1 
0 ··· 0 ck−1 bk

Les polynômes caractéristiques des matrices Bk vérient les relations de ré-

currence

p0 (λ) = 1
p1 (λ) = b1 − λ
pk (λ) = (bk − λ)pk−1 (λ) − c2k−1 pk−2 (λ) pour k = 2, · · · ,n.
Pr. B. RADI 37

Ils vérient les propriétés suivantes

lim pk (λ) = +∞
λ→−∞

Si pk (λ0 ) = 0, alors pk−1 (λ0 )pk+1 (λ0 ) < 0 pour k = 1, . . . ,n − 1. Le polynôme

pk a k racines réelles distinctes qui séparent les (k + 1) racines du polynômes
pk+1 (i.e. x < y < z avec pk+1 (x) = pk+1 (z) = 0 et pk (y) = 0).
Soit a un réel quelconque, si on pose


sgn(pk (a)) si pk (a) 6= 0
sgn(pk (a)) =
sgn(pk−1 (a)) si pk (a) = 0

alors on démontre que le nombre N (k,a) de changements de signes entre

éléments de l'ensemble ordonné {+,sgn(p1 (a)), . . . ,sgn(pk (a))} est égal au
nombre de racines du polynôme pk qui sont strictement inférieures à a.

Algorithme de Givens.
Pour déterminer une valeur propre de la matrice B, on se donne un in-
tervalle arbitraire [a0 ,b0 ] contenant λi . Soit c0 le milieu de l'intervalle [a0 ,b0 ]

Si N (n,c0 ) ≥ i alors λi ∈ [a0 ,c0 [

Si N (n,c0 ) < i alors λi ∈ [c0 ,b0 ]

On restreint alors l'intervalle de recherche à [a1 ,b1 ] dans lequel on peut trouver
λi . On détermine ainsi une suite d'intervalles emboîtés [ak ,bk ] contenant λi
et de longueur (b0 − a0 )/2k .