République Algérienne Démocratique et populaire

Ministère de l’enseignement Supérieur et de la Recherche Scientifique

Université 08 Mai 1945 Guelma

Faculté de Mathématiques et Informatique

Et Science de La Matière

Département de Mathématiques

Polycopie de Cours

Présentés aux Etudiants

Licence Académique en mathématiques

Option : Analyse

Par

Dr. Noura Tabouche

Intitulé

Introduction à l’analyse numérique

Septembre 2016
Table des matières

Introduction 1

1 Résolution numérique de systèmes d’équations 3

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1 Solution théorique . . . . . . . . . . . . . . . . . . . . . . 4
1.1.2 Esprit des méthodes numériques . . . . . . . . . . . . . . . 5
1.2 Méthodes directes . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.1 Méthode d’élimination de Gauss . . . . . . . . . . . . . . . 5
1.2.2 Méthode de factorisation LU . . . . . . . . . . . . . . . . . 11
1.3 Méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.2 Principe des méthodes itératives . . . . . . . . . . . . . . . 14
1.3.3 Le choix de la décomposition . . . . . . . . . . . . . . . . 14
1.3.4 Méthode de Jacobi . . . . . . . . . . . . . . . . . . . . . . 15
1.3.5 Méthode de Gauss-Seidel . . . . . . . . . . . . . . . . . . . 16
1.4 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.4.1 Réponses . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2 Calcul numérique des valeurs propres d’une matrice 24

2.1 Rappels. Dé…nitions . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.1 Rappels élémentaires sur les matrices . . . . . . . . . . . . 25
2.1.2 Décomposition de Schur réelle . . . . . . . . . . . . . . . 27
2.2 Méthode de la puissance . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1 Principe de la méthode et algorithme . . . . . . . . . . . . 28
2.3 Méthode de la puissance inverse . . . . . . . . . . . . . . . . . . . 30
2.4 Méthode de Jacobi pour matrices réelles symétriques . . . . . . . 32
2.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5.1 Réponses . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3 Interpolation numérique 39
3.1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.1 Position du problème . . . . . . . . . . . . . . . . . . . . . 39
3.2 Interpolation polynômiale . . . . . . . . . . . . . . . . . . . . . . 40

1
TABLE DES MATIÈRES 1

3.2.1 Problème . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2.2 Interpolation de LAGRANGE . . . . . . . . . . . . . . . 42
3.2.3 Interpolation dans la base de Newton . . . . . . . . . . . . 45
3.3 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3.1 Réponses . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4 Intégration numérique 51
4.1 Méthode du rectangle . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.2 Méthode du trapèze . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3 Méthode de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.4 Formule des Trapèzes composée . . . . . . . . . . . . . . . . . . . 54
4.5 Formule de Simpson composée . . . . . . . . . . . . . . . . . . . . 57
4.6 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.6.1 Réponses . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5 Résolution numérique des équations di¤érentielles ordinaires 60

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2 Préliminaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.2.1 Problèmes de conditions aux limites . . . . . . . . . . . . . 61
5.2.2 Problèmes de conditions initiales (problème de Cauchy) . . 62
5.2.3 Existence et unicité des solutions . . . . . . . . . . . . . . 63
5.3 Principe d’intégration approchée par des méthodes numériques . . 63
5.4 Méthodes numériques à un pas . . . . . . . . . . . . . . . . . . . 65
5.4.1 Méthode d’EULER . . . . . . . . . . . . . . . . . . . . 65
5.4.2 Méthode de Taylor . . . . . . . . . . . . . . . . . . . . 67
5.4.3 Méthode de Runge- kutta . . . . . . . . . . . . . . . . 70
5.4.4 Méthode de Runge - kutta d’ordre 4 : . . . . . . . . 72
5.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.5.1 Réponses . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

6 Résolution des équations non linéaires 74

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.2 Généralités sur les méthodes de résolution . . . . . . . . . . . . . 74
6.3 Séparation des racines . . . . . . . . . . . . . . . . . . . . . . . . 76
6.3.1 Méthode graphique . . . . . . . . . . . . . . . . . . . . . . 77
6.3.2 Méthodes algébriques . . . . . . . . . . . . . . . . . . . . . 77
6.4 Approximation des racines . . . . . . . . . . . . . . . . . . . . . . 77
6.4.1 Méthode de dichotomie ( ou méthode de la bissection) . . 77
6.4.2 Méthode de Newton-Raphson . . . . . . . . . . . . . . . 80
6.4.3 Méthode de la sécante . . . . . . . . . . . . . . . . . . . . 82
6.4.4 Méthode du point …xe . . . . . . . . . . . . . . . . . . . . 84
6.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.5.1 Réponses . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Introduction

La modélisation des problèmes "réels " observés conduit souvent à un modèle

mathématique mis en équation et dont on cherche à calculer la solution.
Néanmoins, les méthodes analytiques de résolution de problèmes mathéma-
tiques s’adaptent à un nombre limité d’équations algébriques ; notamment, à des
problèmes linéaires. En e¤et, en considérant des hypothèses plus ou moins justi-
…ées il s’ensuit des modèles très simplifés représentant les phénomènes observés
qui sont en réalité très complexes.

La plupart du temps les solutions dégagées à partir des équations simpli…ées

ne représentent le phénomène que dans le domaine où les hypothèses "simpli…ées"
utilisées ont un sens. En plus, lorsque certaines valeurs des paramètres sont
ignorées la solution obtenue n’a pas toujours un grand rapport avec l’observation.
Il est alors naturel de songer à prendre en compte dans les équations tous
les termes pour avoir une approche plus complète du phénomène étudié ; cepen-
dant, la résolution analytique échoue. Il convient de noter que sans moyens de
calculs appropriés, il peut être di¢ cile de résoudre un problème, même si l’on
possède un moyen de résolution ( citons par exemple l’inversion d’une matrice de
grande taille). Ces di¢ cultés augmentent de plus en plus lorsqu’un phénomène
est représenté d’une façon satisfaisante par des équations non linéaires.
Actuellement, le calcul automatique s’est considérablement développé et beau-
coup de méthodes de calcul sont disponibles permettant :
– Le traitement de non-linéarités des équations ainsi que leur résolution ;
– la représentation du phénomène étudié à l’aide de fonctions d’une série de
mesure ;
– la réduction du temps de calcul.

L’objet de l’analyse numérique est de concevoir et fournir des méthodes de

calcul permettant la résolution numérique d’un problème mathématique.

1
TABLE DES MATIÈRES 2

Le présent polycopié est organisé en six chapitres. Pour une revue relative-
ment exhausive, on pourra consulter les ouvrages mis en bibliographie en …n de
ce polycopié.
A la …n de chaque chapitre nous avons proposé des [Link] quelques
corrigés.
Chapitre 1

Résolution numérique de
systèmes d’équations

1.1 Introduction

Soit le système algébrique linéaire de n équations à n inconnues :

X
n
aij xj = bi ; i = 1; :::; n (1.1)
j=1

ou encore 8
>
> a11 x1 + a12 x2 + :::a1n xn = b1
>
>
< a x + a x + :::a x = b
21 1 22 2 2n n 2 0
(1.1 )
>
> ::::
>
>
: a x + a x + :::a x = b
n1 1 n2 2 nn n n

x1 ; x2 ; :::xn sont les inconnues du système.

Le système (1.1) peut s’écrire sous forme matricielle

A:x = b (1.2)

avec 0 1 0 1 0 1
a11 a12 ::: a1n x1 b1
B C B C B C
B a21 a22 ::: a2n C B x2 C B b2 C
A=B
B
C;
C x=B
B
C;
C b=B
B
C:
C
@ ::: ::: ::: ::: A @ ::: A @ ::: A
an1 an2 ::: ann xn bn

3
1. Résolution numérique de systèmes d’équations 4

La notation précédente est équivalente à (l’écriture la plus utilisée) :

0 1
a11 a12 ::: a1n b1
B C
B a21 a22 ::: a2n b2 C
B C
B ::: ::: ::: ::: ::: C
@ A
an1 an2 ::: ann bn

appellée matrice augmentée du système.

Si A est inversible, le système (1:2) a une solution unique qui s’écrit :

x = A 1b (1.3)

Si A est singulière ( i.e. det A = 0) le système n’admet pas de solution.

1.1.1 Solution théorique

On peut mettre (1:3) sous la forme
com (A) 1 com (A)
x= b; car A = :
det A det A
où com (A) est la comatrice de A donnée par com (A) = (cof A)t telle que
cof A = [Aij ](n n) ( matrice des cofacteurs de A ) avec Aij = ( 1)i+j det(aij ); aij
étant une sous-matrice déduite de la matrice A par supression de la ieme ligne et
la jeme colonne.
Ceci constitue une méthode très ine¢ cace car il y a plus de calculs que
nécessaires et plus d’espace mémoire pour moins de précisions

Méthode de Cramer

Si det A 6= 0 alors le système (1:2) a exactement une solution donnée par

det Aj
xj = ; j = 1; 2; :::; n:
det A
Où Aj la matrice obtenue en remplaçant la jième colonne de A par b:
On devra donc calculer (n + 1) déterminants, pour chaque déterminant n:n!
opérations soit (n + 1):n:n! opérations et e¤ectuer n divisions soit au total de
l’ordre de n3 :(n 1)! opérations.

Ces méthodes sont donc lourdes du point de vue programmation du fait

qu’elles demandent un temps de calcul considérable, une place considérable en
1. Résolution numérique de systèmes d’équations 5

mémoire, qu’elles soient imprécises en raison de l’in‡uence d’erreurs d’arrondis

et du fait qu’elles soient di¢ ciles à programmer, car le programme de calcul d’un
déterminant est di¢ cile à écrire.

1.1.2 Esprit des méthodes numériques

Le choix des méthodes de résolution dépend :
- de la taille du système
- du type et de la structure de la matrice A ( forme et symétrie) intervenant
dans le problème.
D’une façon générale on distingue deux types de méthodes :
Les méthodes directes : conduisent à un système équivalent au sys-
tème initial en un nombre …ni d’étapes mais plus facile à résoudre, généralement
on se ramène à un système diagonal ou [Link] mis à part les erreurs
d’arrondis, ces méthodes conduisent à la solution exacte. Cependant, elles ne
sont applicables que lorsque n < 50 environ.
Les méthodes indirecte ou itératives : consistent à fournir la so-
lution en un nombre in…ni d’opérations, par approximations successives et sans
modi…er le système initial. Ces méthodes s’appliquent pour résoudre des systèmes
de grandes tailles.

1.2 Méthodes directes

Les méthodes directes ont pour principe de transformer la matrice augmentée
(le système) en une autre matrice (système équivalent), le nouveau système ayant
la même solution que l’original.

1.2.1 Méthode d’élimination de Gauss

La méthode d’élimination de Gauss a pour but de transformer le système
Ax = b par un système équivalent (c’est-à-dire ayant la même solution) de la
forme U x = bb; où U est une matrice triangulaire supérieure et bb est un second
membre convenablement modi…é.
(k)
La méthode comporte (n-1) étapes On notera aij l’élément aij à l’étape k:
1. Résolution numérique de systèmes d’équations 6

Étape 1 :
- On transforme A en une matrice dont les termes sous diagonaux de la
première colonne sont nuls.
a21
- Pour éliminer le terme a21 on multiplie la ligne l1 par
a11
(1) a21
l2 l2 + l1
a11
on obtient ainsi
(1) a21
a2j = a2j a1j
a11
(1) a21
b2 = b2 b1
a11
a31
- Pour éliminer le terme a31 on multiplie l1 par
a11
(1) a31
l3 l3 + l1
a11
on obtient ainsi
(1) a31
a3j = a3j a1j
a11
:
(1) a31
b3 = b3 b1
a11
- D’une manière générale pour éliminer tout les termes ai1 on utilise la transfor-
mation
(1) ai1
li li + l1
a11

(1) ai1
aij = aij aij
a11
(1) ai1
bi = bi b1
a11
- A la …n de la première étape le système d’équations aura la forme suivante :
0 1
a11 a12 a13 :::: a1n b1
B (1) C
B 0 a(1) 22
(1)
a23
(1)
a2n b2 C
B . .. .. . . . .. C
B . . .. C
@ . . . . A
(1) (1) (1) (1)
0 an2 an3 ::: ann bn

Étape 2 :
Nous éliminons ensuite les termes sous-diagonaux de la seconde colonne.
1. Résolution numérique de systèmes d’équations 7

Comme à la première étape pour éliminer le terme a32 on doit utiliser la

transformation élémentaire
!
(1)
(2) (1) a32 (1)
l3 l3 + (1)
l2
a22
on obtient ainsi
8 !
(1)
>
> (2) (1) a32 (1)
>
> a32 = a32 a22 = 0
>
> (1)
a22 !
>
>
>
> (1) 0 ! 1
>
> (2) (1) a32 (1) (1)
>
> a = a a a32
>
>
33 33
a22
(1) 23 (2) (1)
B a3j = a3j
(1)
a2j C
< B (1) C
a22 !
::::::::::::::::::::::::::: ! ,B
B (1)
C
C
>
> @ (2) a32 A
>
> a32
(1)
b3 = b3
(1) (1)
b1
>
> (2)
a3n = a3n
(1) (1)
a2n (1)
>
> a22
a22
>
> !
>
>
>
> a
(1)
> (2)
> b3 = b3 (1) 32 (1)
b2
: (1)
a22

En continuant de la même façon que l’étape 1 on obtient la …n de la seconde

étape.
Le système d’équations aura la forme suivante :
2 3
a11 a12 a13 :::: a1n b1
6 (1) (1) (1) (1) 7
6 0 a22 a23 ::: a2n b2 7
6 (2) (2) (2) 7
6 0 0 a33 ::: a3n b3 7
6 7
6 .. .. .. . . .. .. 7
4 . . . . . . 5
(2) (2) (2)
0 0 an3 ::: ann bn

Étape k :
A la keme étape; nous éliminons les termes sous-diagonaux de la kieme colonne

(k 1)
!
(k) (k 1) ak+1;k (k 1)
lk+1 lk+1 (k 1)
lk ;
akk
1. Résolution numérique de systèmes d’équations 8

on obtient ainsi
8 !
> (k 1)
>
> (k) (k 1) a k+1;k (k 1)
>
> a = ak+1;k akk = 0
> k+1;k
> (k 1)
akk !
>
> ! (k 1)
ak+1;k
>
> (k 1) (k) (k 1) (k 1)
>
< a(k) (k 1) a k+1;k (k 1) ak+1j = ak+1;j (k 1)
ak;j
k+1;k+1 = ak+1;k+1 (k 1)
ak;k+1 akk !
a , (k 1)
>
> kk ak+1;k
>
> ::::::::::::::::::::::::::
(k) (k 1)
bk+1 = bk+1 b1
>
> ! (k 1)
>
> (k 1) akk
>
> (k) (k 1) ak+1;k (k 1)
>
: ak+1;n+1 = ak+1;j
> (k 1)
akk
ak;n+1

A la …n de la dernière étape on obtient alors le système triangulaire supérieur

suivant :

0 (1) (1) (1) 10 1 0 (1)

1
a11 a12 a1n x1 b1
B (2)
a22
(2)
a2n CB C B (2) C
B 0 CB x2 C B b2 C
B .. CB C B C
B .. . .. CB .. C B .. C
B . . CB . C=B . C (1.4)
B CB .. C B .. C
B .. .. .. CB C B C
@ . . . A@ . A @ . A
(n 1) (n 1)
0 ann xn bn

(k)
Remarque Les opérations précédentes supposent que les termes akk appelés
pivots sont non nuls.

Exemple 1.2.1 Soit le système d’équations représenté par la matrice augmentée

suivante : 2 3
2 1 2 10
6 7
4 6 4 0 26 5
8 5 1 35
(
6
l2 l2 2
l1
8
l3 l3 2
l1
2 3
2 1 2 10
6 7
4 0 1 6 4 5
0 1 7 5
1
l3 l3 l2
1
1. Résolution numérique de systèmes d’équations 9

2 3
2 1 2 10
6 7
4 0 1 6 4 5
0 0 1 1
la solution du système est

x3 = 1; x2 = 2; x1 = 3

Remarque
3
- La méthode de Gauss nécessite 2n3 opérations pour un système de taille n:
- Dans la méthode de Gauss si un pivot est nul on permutera l’équation
correspondante avec l’une de ses suivantes ( permutation de 2 lignes ), donc
pour i = k + 1; n, de manière à mettre un pivot non nul en position k:
Si une telle permutation est impossible la matrice est alors singulière.
(k)
aik
- Si le pivot de valeur proche de zéro c’est-à-dire pivot petit, les quantités (k)
akk
deviennent plus grandes que 1 et ampli…erons les erreurs de troncature ; dans
ce cas il y peut y avoir perte de précision, et on opte pour la modi…cation de
l’algorithme de Gauss, d’où les méthodes suivantes :

Pivot partiel On intervertit les lignes à chaque étape de façon à placer en

pivot le terme de coe¢ cient le plus élevé de la ligne.
C’est la méthode du pivot partiel, à la kieme étape, le pivot est l’élément :
(k) (k)
ai;k = max ap;k
p=k;::::;n

Pivot total On intervertit les lignes et les colonnes de façon à placer en pivot
le terme de coe¢ cient le plus élevé de la matrice : C’est la méthode du pivot
total, à la k-ieme étape, le pivot est l’élément
(k)
ai;j = max a(k)
p;q :
p;q=k;::::;n
0 1 0 1 0 1
1 3 3 x1 0
B C B C B C
Exemple 1.2.2 A = @ 1 1 0 A, x = @ x2 A b = @ 1 A,
3 2 6 x3 11
k = 1 : max(1; 1; 3) = 3 on permute, par exemple, les lignes 1 et 3 on obtient
0 .. 1 0 .. 1
3 2 6 . 11 3 2 6 . 11
B .. C B . C
B 1 1 0 C B
. 1 A!@ 0 3
1
2 .. 38 C
@ A
1 3 3 .. 0 7
0 3 1 .. 3 . 11
.
1. Résolution numérique de systèmes d’équations 10

k = 2 : max( 13 ; 73 ) = 37 .On permute les lignes 2 et 3.

0 .. 1 0 .. 1
3 2 6 . 11 3 2 6 . 11
B C B C
B 0 7 1 ... 11 C ! B 0 7 1 ..
.
11 C
@ 3 3 A @ 3 3 A
0 3 1 .
2 .. 3 8
0 0 715 .. 15
. 7

Il s’ensuit que :
8
>
> x3 = 1
>
> 8
>
> >
< < x3 = 1
Ax = b , x2 = 73 311 x3 , x2 = 2
>
> >
:
>
> x1 = 3
>
>
: 1
x1 = 3 (11 2x2 6x3 )
8
>
< x1 + 3x2 + 3x3 = 2
Exemple 1.2.3 Soit le système 2x1 + 2x2 + 5x3 = 7
>
:
3x1 + 2x2 + 6x3 = 12
Posons : 0 1
.
1 3 3 .. 2
. B .. C
A..b = B @ 2 2 5 . 7 C
A
..
3 2 6 . 12
(1)
k = 1 : max ai;j ; i = 1; 2; 3 j = 1; 2; 3 = 6: La ligne du pivot total sera
alors L3 : En permutant les lignes 1 et 3 on obtient :
0 .. 1
3 2 6 . 12
B C
B 2 2 5 ... 7 C
@ A
..
1 3 3 . 2
La colonne du pivot total est la colonne 3, et on permute les colonnes 1 et 3 :
0 .. 1 0 .. 1
6 2 3 . 12 6 2 3 . 12
B . C B C
B 5 2 2 .. 7 C ! B 0 1 1 .. C
@ A @ 3 2
. 3 A
.. .
3 3 1 . 2 0 2 21 .. 8
(2)
k = 2 : max ai;j ; i = 2; 3 j = 2; 3 = 2:
La ligne du pivot total est alors L3 : Et donc on permute les lignes 2 et 3 :
0 .. 1 0 .. 1
6 2 3 . 12 6 2 3 . 12
B . C B . C
B 0 2 1 .
. 8 C B
! @ 0 2 1 .. 8 C
@ 2 A 2 A
. .
0 13 21 .. 3 0 0 125 .. 35
1. Résolution numérique de systèmes d’équations 11

Du fait de la permutation précédentes des colonnes 1 et 3 on obtient le système

équivalent suivant :

8
>
> 6x3 + 2x2 + 3x1 = 12
>
> 8
>
> > x3 = 1
< <
1
2x2 2 x1 = 8 , x2 = 3
>
> >
:
>
> x1 = 4
>
>
: 5
12 1
x = 35

Remarque Notons bien à chaque permutation de colonnes les inconnues

changent de places.

1.2.2 Méthode de factorisation LU

La méthode L U consiste à transformer la matrice A en un produit de deux
matrices triangulaires inférieures et l’autre supérieure.

A=LU

où L est une matrice triangulaire inférieure.

U est une matrice triangulaire supérieure.
Une fois la décomposition e¤ectuée il faut résoudre deux systèmes à matrices
triangulaires
A x = b , LU x=b
|{z}
Y

On résoudra d’abord le système

Ly=b

une fois y trouvé on résoudra le second système

U x=y

Remarque La méthode impose en général que la diagonale de U soit des

unites .
1. Résolution numérique de systèmes d’équations 12

Exemple 1.2.4 Résoudre le système

2 3
2 1 2 10
6 7
4 6 4 0 26 5
8 5 1 35

0 1 0 10 1
2 1 2 L11 0 0 1 U12 U13
B C B CB C
@ 6 4 0 A = @ L21 L22 0 A @ 0 1 U23 A =
8 5 1 L31 L32 L33 0 0 1
0 1
L11 L11 U12 L11 U13
B C
@ L21 L22 + L21 U12 L21 U13 + L22 U23 A
L31 L32 + L31 U12 L33 + L31 U13 + L32 U23

8 8 8
>
< L11 = 2 >
< L11 U12 = 1 >
< L11 U13 = 2
) (1) L21 = 6 (2) L22 + L21 U12 = 4 (3) L21 U13 + L22 U23 = 0
>
: >
: >
:
L31 = 8 L32 + L31 U12 = 5 L33 + L31 U13 + L32 U23 = 1
8
>
< L11 = 2
L21 = 6
>
:
L31 = 8

8
> 1
>
> L11 U12 = 1 ) U12 =
>
> 2
>
>
<
(2) , L22 + L21 U12 = 4 ) L22 = 1
>
>
>
>
>
>
>
: L +L U =5)L =1
32 31 12 32

8
>
> L11 U13 = 2 ) U13 = 1
>
>
>
>
<
(3) , L21 U13 + L22 U23 = 0 ) U23 = 6
>
>
>
>
>
>
:
L33 + L31 U13 + L32 U23 = 1 ) L33 = 1
1. Résolution numérique de systèmes d’équations 13

0 1 0 1
1
2 0 0 1
1
B C B 2 C
L=@ 6 1 0 A; U =B
@ 0 1 6
C
A
8 1 1 0 0 1
8 8
>
< 2y1 = 10 >
< y1 = 5
Ly = b ) 6y1 + y2 = 26 ) y2 = 4
>
: >
:
8y1 + y2 y3 = 35 y3 = 1

8
> 1
>
> x1 + x2 + x3 = 5
>
> 2
>
>
<
Ux = y ) x2 6x3 = 4
>
>
>
>
>
>
>
: x =1 3

8
> 1
>
> x1 + x2 + x3 = 5
>
> 2 8
>
> >
< < x1 = 3
> x2 6x3 = 4 ) > x2 = 2
>
> :
>
> x3 = 1
>
>
: x3 = 1

1.3 Méthodes itératives

1.3.1 Introduction
Les méthodes itératives consistent à construire une suite de vecteurs
(k) (k) (k)
x = ( x1 ; x2 ; :::; xn )t convergente vers le vecteur solution exacte x = (
(k)
(0) (0) (0)
x1 ; x2 ; :::; xn )t du système (2.2) pour tout vecteur initial x(0) = ( x1 ; x2 ; :::; xn )t
.
1. Résolution numérique de systèmes d’équations 14

1.3.2 Principe des méthodes itératives

Soit le système d’équations suivant :

A x = b:

Ecrivons la matrice A sous la forme suivante

A=M N:

La matrice M doit être inversible. alors le système peut s’écrire :

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

La méthode itérative associée à l’égalité précédente consiste à partir d’un vecteur

initial x(0) à générer la suite x(1) ; x(2) ; :::; x(k+1) de la manière suivante :

x(1) = M 1
N x(0) + M 1
b

x(2) = M 1
N x(1) + M 1
b

:::::::::::::::::::::::::::

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

Cette suite d’égalités peut être représentée par la relation itérative suivante :

x(k+1) = T x(k) + V;

1 1
avec T = M N ;V =M b

1.3.3 Le choix de la décomposition

Les di¤érentes méthodes itératives ne di¤èrent que par le choix de la décom-

position de la matrice A
Soit A la matrice 0 1
a11 a12 ::: a1n
B C
B a21 a22 ::: a2n C
A=B
B
C
C
@ ::: ::: ::: ::: A
an1 an2 ::: ann
1. Résolution numérique de systèmes d’équations 15

Dé…nissons les matrices suivantes

0 1
( a11 0 ::: 0
B C
dii = aii : 8i B 0 a22 ::: 0 C
() D = B
B
C
C
dij = 0 : 8i 6= j @ ::: ::: ::: ::: A
0 0 ::: ann
0 1
( 0 0 ::: 0
B C
lij = aij : 8i > j B a21 0 ::: 0 C
() L = B
B :::
C
C
0 : 8i j @ ::: ::: ::: A
an1 an2 ::: 0
0 1
( 0 a12 ::: a1n
B C
uij = aij : 8i < j B 0 0 ::: a2n C
() U = B
B ::: :::
C
C
0 : 8i j @ ::: ::: A
0 0 ::: 0
A partir de ces matrices on a A = D L U; selon les di¤érentes possibilités
d’écrire cette égalité sous la forme A = M N:

1.3.4 Méthode de Jacobi

La matrice A étant décomposée en

A=M N =D (L + U ) ::::::(D = M )

Le système A x = b peut s’écrire donc sous la forme

D x = (L + U ) x + b:

La méthode de Jacobi s’écrit donc

1
x=D (L + U ) x + D 1 b:

La matrice HJ = D 1 (L + U ) est appelée "matrice de Jacobi" associée à A:

On obtient encore

x(k+1) = D 1
(L + U ) x(k) + D 1 b:
1. Résolution numérique de systèmes d’équations 16

Dans la méthode de Jacobi alors, pour une donnée x(0) ; on calcule x(k+1) selon
la formule
(k+1) bi X aij
xi = xk ; avec aii 6= 0; i = 1; :::; n: (1.5)
aii j=1 aii j
j6=i

La méthode de Jacobi suppose des pivots aii non nuls sinon un simple échange
de ligne su¢ t à résoudre le problème(pour avoir D inversible).

1.3.5 Méthode de Gauss-Seidel

La matrice A étant décomposée en :

A=M N = (D L) U

Le système A:x = b peut s’écrire donc sous la forme

(D L) x = U x + b

La méthode de Gauss-Seidel s’écrit donc

1 1
x = (D L) U x + (D L) b:

La matrice HGS = (D L) 1 U est appelée la matrice de Gauss-Seidel asso-

ciée à A
Soit :
x(k+1) = (D L) 1 U x(k) + (D L) 1 b:
L’inverse de la matrice D L peut être compliqué à calculer, alors on décompose
de la manière suivante :

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

x(k+1) = D 1 L x(k+1) + D 1 U x(k) + D 1 b:

Dans la méthode de Gauss-Seidel alors, pour une donnée x(0) ; on calcule
x(k+1) selon la formule

(k+1) bi X
i 1
aij Xn
aij k
xi = xk+1
j xj ; i = 1; :::; n: (1.6)
aii j=1
aii j=i+1
a ii
1. Résolution numérique de systèmes d’équations 17

Critère d’arrêt

1. En pratique les itérations s’arrêtent quand

x(k+1) x(k) "; " tolérance …xée

2. Pour kT k < 1 si on prends x(0) = V;

on aura x(1) = T x(0) + V ) x(1) x(0) = T V d’où

(k) x(k+1) x(k) kT kk+1 kV k

x x
1 kT k 1 kT k
et en posant
kT kk+1 kV k
"
1 kT k
ce qui nous permet d’estimer le nombre su¢ sant d’itérations k pour approximer
x avec la précision donnée ":

Etude de convergence

Dé…nition 1.3.1 :Une matrice A est dite à diagonale strictement dominante

si :
X
n
jaii j > jaij j ; 8 i: 2 f1; ::; ng
j=1;j6=i
ou
X
n
jajj j > jaij j ; 8 j: 2 f1; ::; ng
i=1;i6=j

Cette dé…nition signi…e que le terme diagonal aii de la matrice A est net-
tement dominant puisque sa valeur absolue est plus grande que la somme des
valeurs absolues de tous les autres termes de la ligne

Proposition 1.3.1 Si la matrice A est à diagonale strictement dominante, les

méthodes de Jacobi et de Gauss-Seidel convergent.

La méthode de Gauss-Seidel est plus rapide que celle de Jacobi.

Exemple 1.3.1 Considérons le système linéaire

0 10 1 0 1
4 2 1 x 4
B CB C B C
@ 1 2 0 A@ y A = @ 2 A
2 1 4 z 9
1. Résolution numérique de systèmes d’équations 18

mis sous la forme 8

y z
>
> x=1 ;
>
> 2 4
>
>
<
x
y =1+
>
>
2
>
>
>
>
: y
z = 94 x
2 4
:
Soit x(0) = (0; 0; 0) le vecteur initial.
- En calculant les itérés avec la méthode de Jacobi en considérant la suite
8
>
> xn+1 = 1 12 yn 41 zn ;
>
>
>
>
<
yn+1 = 1 + 12 xn
>
>
>
>
>
>
:
zn+1 = 94 21 xn 14 yn :
on trouve
0 1 0 1 0 1 0 1
1 1 5
1 16 8 128
B C B C B C B C
B C B C B C B C
B C B C B C B C
x(1) B
=B 1 C
C; x(2) =B
B
3
2
C;
C x(3) =B
B
1
32
C;
C x(4) =B
B
15
16
C:
C
B C B C B C B C
@ A @ A @ A @ A
9 3 61 265
4 2 32 128

La suite x(k) converge vers (0; 1; 2) la solution du système au bout de quatre

itérations.
-En calculant les itérés avec la méthode Gauss-Seidel en considérant la suite
8
>
> xn+1 = 1 21 yn 14 zn ;
>
>
>
>
<
yn+1 = 1 + 21 xn+1
>
>
>
>
>
>
:
zn+1 = 94 12 xn+1 14 yn+1 :
on trouve
0 1 0 1 0 1
3 9
1 32 1024
B C B C B C
B C B C B C
B C B C B C
x(1) = B
B
3
2
C;
C x(2) =B
B
61
64
C;
C x(3) =B
B
2047
2048
C;
C
B C B C B C
@ A @ A @ A
11 527 16349
8 256 8192
1. Résolution numérique de systèmes d’équations 19

La suite x(k) converge vers (0; 1; 2) la solution du système au bout de trois

itérations.

1.4 Exercices
Exercice 1 Soit le système d’équations linéaires suivant
8
>
< 3x1 2x2 + x3 = 2
(S) 2x1 + x2 + x3 = 0
>
:
4x1 3x2 + 2x3 = 4
1) Ecrire la forme matricielle Ax = b de ce système et le résoudre à l’aide de
la méthode de Gauss.
2) Calculer le déterminant de la matrice A:
3) Factoriser la matrice A du système en produit LU , puis résoudre ce sys-
tème.
Exercice 2 Résoudre le système d’équations linéaires suivant
8
>
< 3x1 + x2 x3 = 2
x1 + 5x2 + 2x3 = 17
>
:
2x1 x2 6x3 = 18
par la méthode de Jacobi, ensuite par la méthode de Gauss-Seidel avec x(0) =
(0; 0; 0)T :
Exercice 3 Soit le systèmes d’équations linéaires suivant
8
>
< 2x1 x2 + x3 = 1
x1 + 2x2 x3 = 2
>
:
x1 x2 + 2x3 = 3

1) Ecrire ce système sous la forme matricielle Ax = b et montrer que ce

système admet une solution unique.
2) Peut-on décomposer la matrice A sous la forme A = L LT ; où L est une
matrice triangulaire inférieure ?
3) Si la réponse est oui, résoudre ce système par la méthode de Cholesky, en
déduire la valeur du déterminant de A:
1. Résolution numérique de systèmes d’équations 20

Exercice 4 Soit le systèmes d’équations linéaires suivant

8
>
< 4x1 + x2 + 2x3 = 4
(S) 3x1 + 5x2 + x3 = 7
>
:
x1 + x2 + 3x3 = 3
Donner le nombre su¢ sant d’itérations pour que les deux méthodes ( Jacobi
et Gauss-Seidel ) donnent une solution approchée de (S) à 10 2 près.

1.4.1 Réponses
Exercice 1
1) Le système (S) sous sa forme matricielle
2 32 3 2 3
3 2 1 x1 2
6 76 7 6 7
42 1 15 4x2 5 = 405
4 3 2 x3 4
2 3 2 3 2 3
3 2 1 x1 2
6 7 6 7 6 7
avec la matrice A = 42 1 15 ; x = 4x2 5 et b = 405 :
4 3 2 x3 4
det A = 5 alors le système (S) a une solution admet une solution unique.
2 . 3 2 . 3
3 2 1 .. 2 3 2 1 .. 2
6 . 7 2 4 6 7
62 1 1 .. 07 L L et L L 60 7 1 ... 47
4 5 2
3
1 3
3 4
1
3 3 35
.. 1 2 .. 4
4 3 2 . 4 0 3 3
. 3
2 . 3
3 2 1 .. 2
1 6 .. 7 8 2 2
L3 + L2 6
40 7 1
. 4 7 =) x = ; x =
5 3 2 et x1 =
7 3 3 3 5 5 5
.. 24
0 0 15 21
. 21

La solution du système est x = ( 25 ; 52 ; 85 )T :

2 3 2 3
3 2 1 3 2 1
6 7 6 1 7
2) det A = det 40 73 13 5 = det 40 37 3 5
= 5:
1 2 15
0 3 3
0 0 21
1. Résolution numérique de systèmes d’équations 21

3)
2 32 3
L11 0 0 1 U12 U13
6 76 7
A = L U = 4L21 L22 0 5 40 1 U23 5
L31 L32 L33 0 0 1
2 3 2 3
L11 L11 U12 L11 U13 3 2 1
6 7 6 7
= 4L21 L21 U12 + L22 L21 U13 + L22 U23 5 = 42 1 15
L31 L31 U12 + L32 L31 U13 + L32 U23 + L33 4 3 2
ce qui donne L11 = 3; L11 U12 = 2 =) U12 = 23 ; L21 = 2
et L21 U12 + L22 = 1 =) L22 = 73 ;
L31 = 4 et L31 U12 + L32 = 3 =) L32 = 13 :
L11 U13 = 1 =) U13 = 31 et L21 U13 + L22 U23 = 1 =) U23 = 17 :
L31 U13 + L32 U23 + L33 = 2 =) L33 = 57 :
2 32 2 1
3
3 0 0 1 3 3
6 76 7
A = LU = 42 37 0 5 40 1 17 5
1 5
4 3 7
0 0 1
8
>
< 3y1 = 2
Ly = b =) 2y1 + 37 y2 = 0
>
:
4y1 13 y2 + 75 = 4
ce qui donne y1 = 23 ; y2 = 74 ; y3 = 85 .
Ensuite, on résoud U X = Y
2 2 1
32 3 2 2
3
1 3 3
x1 3
6 17 6 7 6 47
40 1 7 5 4x2 5 = 4 75
8
0 0 1 x3 5

d’où la solution x = (x1 ; x2 ; x3 )T = ( 52 ; 25 ; 85 )T :

Exercice 2
Indication on pose A = D 2 (L + U ), 3 Ax 0 2 = b =) x3 1 D21 (L + U )x +3D
1
1
b
1
3
0 0 0 0 1 0 1 1
1 6 1 7 B6 7 6 7C
avec Hj = D (L + U ) = 4 0 5 0 5 @4 1 0 05 + 40 0 25A ;
1
0 0 6
2 1 0 0 0 0
(Hj ) < 1:
On dé…nit la suite x(k) par
(
x(k) = D 1 (L + U ) x(k 1) + D 1 b
x(0) = (0; 0; 0)
1. Résolution numérique de systèmes d’équations 22

Il s’ensuit que
0 1
1 B X
3
C
81 i 3;
(k)
xi = Bb i aij xj
(k 1) C
aii @ A;
j=1
i6=j

où 8 (k) (k 1) (k 1)
>
< x1 = 13 (2 x2 + x3 )
(k) (k 1) (k 1)
x = 51 (17 x1 2x3 )
> (k)2
: 1 (k 1) (k 1)
x3 = 6 ( 18 2x1 + x2 ):
A l’aide d’un calcul simple et en faisant varier k = 0; ::::; 7 ( 7 itérations)
(7) (7) (7)
on déduit que la (x1 ; x2 ; x3 )T converge vers la solution exacte (1; 2; 3)T :
Pour la méthode de Gauss-Seidel on pose

A = (D L) U)
Ax = b =) (D L) x = U x + b
1 1
x = (D L) U x + (D L) b

On dé…nit la suite x(k) par

(
x(k) = (D L) 1 U x(k 1) + (D L) 1
b
; (HGS ) < 1:
x(0) = (0; 0; 0)

Il s’ensuit que
!
(k) 1 X
i 1
(k)
X (k 1)
81 i 3; xi = bi aij xj aij xj
aii j=1 j=1+1

où 8 (k) (k) (k 1)
>
< x1 = 31 (2 x2 + x3 )
(k) (k) (k 1)
x2 = 15 (17 x1 2x3 )
>
: (k) 1 (k) (k)
x3 = 6 ( 18 2x1 + x2 ):
Au bout de 6 itérations en faisant k = 0; ; 6 on trouve x(6) = (0; 9914; 2; 0057; 2; 9962)T
qui converge vers la solution exacte (1; 2; 3)T :
Exercice 4 2 3 21 3
1 1
0 4 2 4
0 0
6 17 6 1 7
HJ = D 1 (L + U ) = 4 35 0 5 5 avec D 1
= 4 0 5
05
1 1
3 3
0 0 0 13
et kHJ k1 = 45 < 1 méthode de Jacobi converge.
1. Résolution numérique de systèmes d’équations 23

Posons
21
32 3
4
0 0 4
1 6 1 7 677
c = D b = 40 5 05 4 5 5
0 0 13 1
x = HJ x + c:

On sait que
(k) kHkk
x x x(1) x(0) :
1 kHk
Puisque x(k) = HJ x(k 1)
+ D 1 b =) x(1) x(0) = HJ c
et
4 k+1 7
(k) kHJ kk+1 kck 5 5 2
x x = 10 =) k = 29:
1 1 kHJ k 1 45
De même on a HGS = (D L) 1 U et x(k) = HGS x(k 1) + c avec
2 3 21 3
4 0 0 4
0 0
6 7 6 7
D L = 43 5 05 et (D L) 1 = 4 203 15 0 5
1 1 1
1 1 3 30 5 3
2 3
0 41 21
6 3 1 7 3
Alors HGS = 40 20 10 5
et kHGS k1 = < 1 GS est CV.
1 2
4
0 30 15
2 3
1
1 647
c = (D L) b = 4 5 5 et kck1 = 1;
2
5
k+1
kHGS k kck 3 k+1 1
x x(k) 1 1 kHGS k
= 4 1 34
10 2
=) k = 20:
Chapitre 2

Calcul numérique des valeurs

propres d’une matrice

Le calcul des valeurs et vecteurs propres est indispensable dans plusieurs ap-
plications. Notamment, dans la théorie de la stabilité des systèmes dynamiques,
en physique, mécanique, chimie....

Dans ce chapitre nous allons aborder le problème du calcul de valeurs et

vecteurs propres d’une matrice qui s’avère être un problème ardu. En e¤et, dé-
terminer les valeurs propres d’une matrice revient à calculer les racines de son
polynôme caractéristique. Or, il ne su¢ t pas de factoriser ce polynôme a …n
d’obtenir les racines, d’après les résultats de Galois, N. Abel et P. Ruffini il
n’existe pas toujours de possibilité d’exprimer les racines d’un polynôme de de-
gré supérieur ou égal à cinq à partir des coe¢ cients du polynôme et d’opérations
élémentaires. En …n de compte, pour le calcul de valeurs propres d’une matrice
il est nécessaire d’avoir recours aux méthodes itératives.

Nous commençons par donner quelques rappels et dé…nitions qui seront utiles
par la suite.

24
2. Calcul numérique des valeurs propres d’une matrice 25

2.1 Rappels. Dé…nitions

2.1.1 Rappels élémentaires sur les matrices
- On note A = (aij ) une matrice avec 1 i m (lignes)
et 1 j n (colonnes).
- On appelle transposé de A et l’on note At = (aji ):
- La matrice adjointe est donnée par A = (aji ):
- Une matrice carrée A d’ordre n est inversible s’il existe une matrice, notée
A 1 telle que A A 1 = A 1 A = In ; où In est la matrice identité.
- Soit A une matrice carrée réelle, alors
A est symétrique si At = A ;
A est orthogonale si A At = At A = In :
- Soit A une matrice carrée complexe, alors
A est hermitienne si A = A ;
A est unitaire si A A = A A = In ;
A est normale si A A = A A:
- On dit qu’une matrice A réelle symétrique est dé…nie positive si pour tout
u 2 Rn ;
u 6= 0; ut Au > 0; la matrice A est dite semi-dé…nie positive si pour tout
u 2 Rn ; ut Au 0:

Normes matricielles

Dé…nition 2.1.1 Une norme matricielle sur l’espace vectoriel Mn des matrices
carrées d’ordre n; est une application de Mn vers R+ qui véri…e pour tous A; B 2
Mn ; 2 R (ou C)

kAk = 0 , A = 0
k Ak = j j kAk
kA + Bk kAk + kBk
kA Bk kAk kBk

Dé…nition 2.1.2 ( Norme matricielle subordonnée) Soit k:k une norme vec-
torielle, on appelle norme matricielle subordonnée (à cette norme vectorielle )
2. Calcul numérique des valeurs propres d’une matrice 26

l’application de Mn vers R+ dé…nie par

kA xk
jkAkj = sup = sup kA xk = sup kA xk :
x2V;x6=0 kxk x2V;x 1 x2V; kxk=1

Les normes matricielles subordonnées suivantes sont couramment utilisées

Xn
kAk1 = max jaij j "maximum sur les lignes"
1 i n
j=1
Xn
kAk1 = max jaij j "maximum sur les colonnes"
1 j n
j=1
p
kAk2 = (At A)

où (At A) est le rayon spectral de la matrice symétrique semi-dé…nie positive

At A:

Dé…nition 2.1.3 On appelle valeur propre de A, toute racine i 2 C du poly-

nôme caractéristique p( ) = det(A In ): La multiplicité algébrique mi 2 N
est l’ordre de multiplicité de i en tant que racine de p ( ie : multiplicité dans le
polyn^me caractéristique). Si on note r n le nombre de racines distinctes du
polynôme p; on a
Y
r
mi
p( ) = ( i) ; où r nombre de racines distinctes de p; (r n):
i=1

-Un vecteur propre vi 2 Cn nf0g; associé à la valeur propre i; véri…e

Avi = i vi :

On a donc vi 2 ker(A i In )nf0g et on appelle multiplicité géométrique de i

l’entier naturel gi = dim(ker(A i In )); 1 i r:

Remarque 2.1.1 Une valeur propre est dite simple si sa multiplicité algébrique
est 1: Sinon, c’est une valeur propre multiple.

Dé…nition 2.1.4 L’ensemble des valeurs propres de A est appelé le spectre de

A et est noté (A) ou spec (A):

(A) = spec (A) = f 2 C = det(A In ) = 0g = f 1 ; 2 ; :::; ng :

Dé…nition 2.1.5 On appelle rayon spectral de la matrice A le maximum du

module de ces valeurs propres. On le note (A):

(A) = max fj j = 2 spec (A)g = max fj 1 j ; j 2 j ; :::; j n jg :

1 i n
2. Calcul numérique des valeurs propres d’une matrice 27

Matrices diagonalisables, matrices semblables

La proposition suivante fait l’objet de rappel de quelques résultats reliants la

multiplicité des valeurs propres à la réduction des matrices.

Proposition 2.1.1 Soit A une matrice carrée d’odre n, on a les propriétés sui-
vantes :
1- Les multiplicités véri…ent les relations :
X
r X
r
pour 1 i r : gi mi ; mi = n et gi n:
i=1 i=1

X
r
2- La matrice A est diagonalisable si et seulement si gi = n:
i=1
3- Si toutes les valeurs propres i 2 C; 1 i n; sont simples (i.e. mi =
gi = 1 et r = n), alors la matrice est diagonalisable (dans C).
4- La matrice A n’est pas diagonalisable si et seulement si il existe au moins
une valeur propre i telle que mi > gi 1:

Dans ce cas, on dit que la matrice respectivement la valeur propre est défec-
tive.

Proposition 2.1.2 Soit P une matrice carrée d’odre n inversible. On dit que
B = P 1 AP et A sont des matrices semblables. Alors A et B ont les mêmes
valeurs propres et v est un vecteur propre de A si et seulement si P 1 v est un
vecteur propre de B:

2.1.2 Décomposition de Schur réelle

La décomposition de Schur réelle permet de traiter tous les cas de matrice

réelle A diagonalisable ou défective. D’où l’objet de la proposition suivante :
2. Calcul numérique des valeurs propres d’une matrice 28

Proposition 2.1.3 (Décomposition de Schur réelle) Soit A une matrice carée

réelle, alors il existe une matrice réelle et orthogonale U et une matrice triangu-
laire supérieure par blocs T tels que U t AU = T . Sur la diagonale de T se trouvent
des blocs 1 1 correspondant aux valeurs propres réelles de A et des blocs 2 2
correspondant aux valleurs propres complexes conjuguées de A:

Exemple 2.1.1
0 1 0 10 10 1
1 0 0 1 0 0 1 0 1 1 0 0 0 1 0 0
B C B CB CB C
B 0 1 1 0 C B 1 0 0 0 CB 0 1 0 0 CB 0 0 1 0 C
A = BB 0 0
C=B
C B
CB
CB
CB
CB
C
C
@ 1 0 A @ 0 1 0 0 A@ 0 0 1 1 A@ 1 0 0 0 A
1 0 0 1 0 0 0 1 0 0 1 1 0 0 0 1
= U T U t:

Ici spec(A)=f 1 = 1; 2 = 1 + i; 3 = 1 ig ; donc r = 3

et m2 = 2 > g1 = 1; m2 = m3 = g2 = g3 = 1: La matrice est défective car
1 est une valeur propre défective.

2.2 Méthode de la puissance

la méthode de la puissance est la méthode la plus simple pour le calcul des
valeurs propres d’une matrice A. Cette méthode permet de fournir une approxi-
mation de la valeur propre du plus grand module (valeur propre dominante) et
un vecteur propre associé. On obtient donc le rayon spectral d’une matrice.

2.2.1 Principe de la méthode et algorithme

La méthode des puissances appelée aussi méthode de la puissance itérée part
du principe suivant : en appliquant un grand nombre de fois la matrice sur
un vecteur arbitraire x(0) normalisé, on obteint des vecteurs successifs qui se
rapproche de la direction du vecteur propre associé à la plus grande valeur propre
en module.
Soit A une matrice carrée diagonalisable, on suppose que ces valeurs propres
sont ordonnées comme suit : j 1 j > j 2 j ::::: j n j et v1 ; v2 ; :::; vn les vecteurs
proppres associés.
2. Calcul numérique des valeurs propres d’une matrice 29

La méthode consiste à calculer les suites des itérées dé…nies par

y (k) (k)
y (k) = Ax(k 1)
; x(k) = (k)
et e = (x(k) )t Ax(k)
ky k2

Analysons les propriétés de la suite x(k) k2N : Une simple récurrence sur
l’indice k donne
(k) Ak x(0)
8k 2 N; x = (2.1)
kAk x(0) k2
et l’on voit la présence des puissances de la matrice A ce qui justi…e le nom
de cette méthode. l’ensemble fv1 ; :::; vn gde vecteurs propres de A forment une
base, alors le vecteur x(0) peut se décomposer de la façon suivante :
X
n
x(0) = i vi
i=1

k
notons que Avi = i vi ) Ak vi = i vi et l’on obtient, en supposant le coe¢ cient
1 6= 0 et pour tout k 2 N

X
n
Ak x(0) = i (A
k
vi ) (2.2)
i=1
(2.1)
X
n
k
= i i vi
i=1
(2.2)
!
k k
k 2 n
= 1 1 v1 + v2 + : : : + n vn
1 1
!
X
n
i i
k
k
= 1 1 v1 + vi :
i=2 1 1

Comme j 1 j > j i j ; 8i 6= 1, on a i
1
< 1 pour 2 i n; la composante
le long de v1 du vecteur x(k) augmente en module avec l’entier k, par contre les
composantes suivant les autres directions vi , i = 2; : : : ; n, diminuent. Les vecteurs
fvi gi=1;::::;n sont tels que kvi k2 = 1; il s’ensuit que

X
n
i i
k X
n
i i
k
2
k X
n
i 2
k
vi =C ; (2.3)
i=2 1 1 i=2 1 1 1 i=2 1 1
2
2. Calcul numérique des valeurs propres d’une matrice 30

en l’on déduit de ( 2.1), de la dernière égalité de (2.2) et de (2.3) que

k
1 1 (v1 + rk )
8k 2 N; x(k) = k
1 1 (v1 + rk ) 2

où la suite de vecteurs (rk )k 2 a pour limite le vecteur nul quand k ! +1: Par
conséquent, le vecteur x(k) devient peu à peu colinéaire avec le vecteur propre
v1 associé à la valeur propre dominante 1 quand k ! +1 d’autant plus rapi-
dement que 12 est petit. La suite 8k 2 N; (x(k) )t Ax(k) = (k) ; converge donc
vers la valeur propre 1 :

Nous présentons maintenant, l’algorithme de la méthode de la puissance

Algorithme 1
(0) (0)
x(0) = x1 ; :::; xn 2 Rn vecteur arbitraire donné, k = 0
répéter
y (k+1) = Ax(k)
(k+1)
x(k+1) = yy(k+1)
k k2
(k+1)
e = (x(k) )t Ax(k)
k =k+1
jusqu’à [test d’arrêt].
(k+1)
Le test d’arrêt est souvent basé sur e e(k) ; c’est-à-dire que l’ité-
ration s’arrête dès que la di¤érence entre deux estimations de la valeur propre
est su¢ samment petite.

2.3 Méthode de la puissance inverse

Cette méthode permet d’approcher la valeur propre de la matrice A la plus
proche d’un nombre donné n’appartenent pas à son spectre.
On applique la méthode de la puissance à la matrice (A In ) 1 dont les
vecteurs propres sont fv1 ; :::; vn g ceux de la matrice A et les valeurs propres sont
( i ) 1 ; 1 i n:

Il su¢ t d’adapter l’algorithme1 pour obtenir :

2. Calcul numérique des valeurs propres d’une matrice 31

Algorithme 2
Méthode de la puissance inverse
x(0) 2 Rn donné, k = 0
répéter
y (k+1) = (A In ) 1 x(k)
(k+1)
x(k+1) = yy(k+1)
k k
(k+1)
e = (x ) Ax(k)
(k) t

k =k+1
jusqu’à [test d’arrêt]

Remarque
On calcule une décomposition LU de la matrice A In a…n de résoudre le
système
(A In ) y (k+1) = x(k) :
0 1
3 0 0
B C
Exemple 2.3.1 Soit A = @ 17 13 7 A ; avec spec(A) = f 1 = 6; 2 = 3; 3 = 1g
16 14 8
. En appliquant la méthode de la puissance avec
0 1 0 1
0:0938417 0:0002303
0 1 B C B C
1 B C B C
B C B C B C
(0)
x = @ 1 A on obtient x =B
(1)
B 0:7194529 C; x
C
(10)
= B 0:7070492
B
C
C
B C B C
1 @ A @ A
0:6881724 0:7071643
(1)
et e = 4:1986301; e(10) = 6:0036631:

Si l’on applique la méthode de la puissance inverse, avec = 0; et en partant

de
0 1 0 1
0:13938466 0:0000151
0 1 B C B C
1 B C B C
B C B C B C
x(0) = @ 1 A on obtient x(1) =B
B 0:6270597 C;
C x(10) =B
B 0:4472257
C
C
B C B C
1 @ A @ A
0:7664063 0:8944211
(1)
et e = 0:5242718; e(10) = 1:0000478:
2. Calcul numérique des valeurs propres d’une matrice 32

2.4 Méthode de Jacobi pour matrices réelles sy-

métriques
Soit A = (aij )1 i;j n une matrice réelle carrée symétrique. La méthode de Ja-
cobi consiste à construire une suite de matrice A = A0 ; A1 ; :::; Ak ; :::; qui converge
vers une matrice diagonale D: Avec

Ak+1 = Jkt Ak Jk ;

où Jk est une rotation de Jacobi (voir la dé…nition ci-après) et Jkt Jk = In :

Si l’on pose J = J0 J1 :::Jk 1 ; on aura une approximation de D ' J t A J;
respectivement une approximation de A ' J D J t : Comme J est encore une
matrice orthogonale, D contient une approximation de toutes les valeurs propres
de A et les colonnes de J forment une base orthonormée de vecteurs propres.
Une rotation de Jacobi est dé…nie grâce à une matrice de Givens : Les matrices
de rotation de Givens sont des matrices orthogonales qui permettent, d’annuler
certains coe¢ cients d’un vecteur ou d’une matrice. Pour un couple d’indices i et
j véri…ant 1 i < j n; et un nombre donnés, on dé…nit la matrice de Givens
comme suit :
0 1
1 0 0 0 0
B .. C
B 0 1 . C
B C
B .. ... C
B . 0 C
B C
B .. C
B 0 c s . C
G(i; j; ) = B
B .. ..
C
C
B . . 0 C
B .. C
B C
B 0 s c . C
B .. C
B . 1 0 C
@ A
0 0 0 0 1
avec c = cos et s = sin : Pour x 2 Rn ; G(i; j; ) x représente une rotation
d’angle dans le plan (ei ; ej ):
Le produit G(i; j; )t A n’a¤ecte que les i e et j e lignes de A.
Le AG(i; j; )A n’a¤ecte que les i e et j e colonnes de A:

La matrice de Givens est construite à partir d’une valeur de qui permet

d’annuler les éléments aij et aji :
2. Calcul numérique des valeurs propres d’une matrice 33

L’angle véri…e
aii ajj
cot g(2 ) = =
2aij
De cette relation on peut déduire les valeurs de c et s grâce à l’équation
s
t= = tg ; t2 + 2 t 1 = 0:
c
e
L’algorithme 5 suivant remplace la matrice A par une matrice A
telle que e
aij = e
aji = 0:

Algorithme 3
Rotation de Jacobi pour annuler les éléments de coordonnées fi; jg
Si jaij j > alors
aii ajj
= 2aij
signe( )
t= p
j j+ 1+ 2

1
c = p1+t2 et s = cte
Ae = G(i; j; )t AG(i; j; ):

L’algorithme 4 qui suit donne la convergence vers une matrice diagonale

contenant les valeurs propres, et une matrice J contenant les vecteurs propres
de A:

Algorithme 4
Algorithme de Jacobi pour valeurs et vecteurs propres d’une matrice symé-
trique
Seuil > 0 donné, J = In ; k = 0
répéter
Choisir fi; jg et en déduire Jk = G(i; j; ) par l’algorithme 3
J = J Jk
k =k+1
jusqu’à max1 k<l n jakl j < :

Remarque 2.4.1 Dans l’algorithme original de Jacobi, on choisit les coordon-

nées fi; jg de façon à avoir

jaij j = max jakl j :

1 k<l n
2. Calcul numérique des valeurs propres d’une matrice 34

Ceci permet d’avoir une bonne vitesse de convergence. Cependant, pour n assez
2
grand, la recherche du maximum parmi les n 2 n éléments en dehors de la diago-
nale prends du temps. En pratique on parcourt A simplement ligne par ligne par
exemple .
0 1
2 1 3
B C
Exemple 2.4.1 Soit A = @ 1 6 3 A et l’on …xe = 10 5 :
3 3 0

On a
spec(A) = f 7:7155084; 0:5472851; 4:2627936g

0 1
0:8112422 0 0:5847103
B C
Pour J0 = @ 0 1 0 A;
0:5847103 0 0:8112422

on annule a13 = a31 = 3:

0 1
16
4:1622777 0:9428887 2:2210
t B C
A1 = J0 AJ0 = @ 0:9428887 6 3:0184368 A :
16
2:2210 3:0184368 2:1622777
0 1
1 0 0
B C
Après 5 itérations, on a J4 = @ 0 0:9999904 0:0043757 A et
0 0:0043757 0:9999904

0 1
4:2627928 0:0029793 0:0000130
B C
A5 = J4t A4 J4 = @ 0:0029793 7:7155077 6:94 10 18 A :
0:0000130 1:2210 16 0:5472851

0 1
0:7766476 0:2219207 0:5895504
B C
En plus J = J1 J2 J3 J4 = @ 0:1060891 0:8764440 0:4696712 A :
0:6209376 0:4273139 0:6571448
2. Calcul numérique des valeurs propres d’une matrice 35

2.5 Exercices
Exercice1 Soit la matrice :
2 3
1 2 0
6 7
A = 42 1 05
0 0 1

Déterminer ses valeurs et vecteurs propres

1) Par la méthode de la puissance itérée
2) Par la méthode de Jacobi.
Exercice2 Soit la matrice
" #
2 1
A=
1 2

Déterminer sa plus grande valeurs propre ainsi que le vecteur propre associé
par la méthode de la puissance itérée.

2.5.1 Réponses
Exercice1
1) Résumons un peu la méthode :
En pratique, on se donne un vecteur initial x(0) (de norme égale à 1 ou
quelconque),
et on calcule les x(k) par la relation x(k+1) = A x(k) :
De plus, on préfère ramener à 1 la valeur de l’une des composantes de ces
vecteurs ( par exemple la 1ere ):
La composante correspondante du vecteur itéré suivant donne le rapport à
la valeur précédente et tend vers la valeur propre de plus grand module.
A présent on veut trouver les valeurs propres suivantes, on applique la mé-
thode de défalation qui consiste a créer la matrice
T
v1 x(1)
A1 = A 1 T
(x(1) ) v1

où le vecteur x(1) est le vecteur propre de la matrice transposée de A;correspondant

à la valeur propre 1 :
2. Calcul numérique des valeurs propres d’une matrice 36

2 3
1
6 7
On se donne un vecteur initial x(0) = 405 à partir duquel on construit la
0
suite de vecteurs :
x(1) = A x(0) ; x(2) = A x(1) ; ::::; reporté dans le tableau ci-dessous
13 14 121 365
1 5 5 13 41 121
::::
(k)
x 1 1 1 1 1 1 1
0 2 45 14
13
40
41
122
121
364
365
:::::
0 0 0 0 0 0 0

Les composantes de chaque vecteur de ce tableau ont été divisées par la première
composante correspondant à chaque vecteur. 2 3
1
(k) 6 7
La suite des vecteurs x tend vers le vecteur v1 = 415 et que la suite des
0
valeurs de converge vers la valeur 1 = 3:
La matrice A considérée étant symétrique, le vecteur
2 3propre de sa transposée,
1
(1) 6 7
correspondant à la valeur propre 1 sera x = v1 = 415 :
0
2 1 1
3
T 2 2
0
v1 (x(1) ) 6 1 1 7
Construisons la matrice A1 qui est telle que A1 = A 1 T =4 2 2
0 5:
(x(1) ) v1
0 0 1
2 3
1
6 7
On se donne à nouveau un vecteur x(0) = 405 et on procède de la même
0
manière qu’avec la matrice A :
1
2
-1 -1 ....
(k)
x 1 1 1 1
0 -1 -1 -1 ....
0 0 0 0
3 2
1
6 7
On voit que cette suite de vecteurs converge vers le vecteur propre : v2 = 4 15
0
et que la valeur propre qui lui correspond est 2 = 1:
2. Calcul numérique des valeurs propres d’une matrice 37

On recommence le processus en dé…nissons une matrice A2 à partir de A; x(2) =

v2 : 2 3
(2) T 0 0 0
v2 x 6 7
A2 = A1 2 T
= 40 0 0 5
(2)
(x ) v2
0 0 1
Il est clair que la valeur propre2de
3 cette matrice est 3 = 1 et que le vecteur
0
6 7
propre correspondant est v3 = 405
1
2) Le plus grand élément non diagonal de la matrice A étant l’élément
a12 ; l’angle sera tel que :
1 2a12
tg (2 ) = = :
cot g(2 ) a11 a22

Comme a11 = a22 ; = 4 ; la matrice de transformation G12 est

2 p2 p
2
3
2 2
0
6 p2 p
2 7
G12 = 42 2
05
0 0 1

d’où 2 3
3 0 0
6 7
A1 = G12 AG12 = 40 1 05
0 0 1
matrice diagonale. Les valeurs propres de A1 ; et par conséquent de A sont les
éléments de la diagonale : 1 = 3 et 2 = 3 = 1:
Quand aux vecteurs propres, puisque le calcul s’est opéré en une seule étape,
ils sont données directement par les colonnes de la matrice G:
On a donc :
2 3 2 3 2 3 2 3
1 1 0 1
6 7 6 7 6 7 6 7
V1 = G12 405 = 415 ; V2 = G12 415 = 4 15
0 0 0 0

et 2 3 2 3
0 0
6 7 6 7
V3 = G12 405 = 405 :
1 1
2. Calcul numérique des valeurs propres d’une matrice 38

Exercice2 " #
1
En partant du vecteur initial x(0) = ; l’algorithme de la puissance
0
x(k) = Ax(k 1) fournit les valeurs suivantes

k 0 1 2 3 4 5 6 7
x(k) 1 2 5 14 41 122 365 1094
0 1 4 13 40 121 364 1093
2:928 2:975 2:991 2:9973
(k)
x
Le rapport de 2 composantes homologues successives ( (ki 1) ) est donné dans la
xi
dernière ligne, dans la dernière colonne,
il converge vers la valeur 3, qui serait la plus grande des valeurs propres.
" #
1
Le vecteur x(k) semble tendre vers une limite proportionnelle au vecteur
1
qui serait le vecteur propre v1 associé à la valeur propre 1 = 3:
Remarque : Comme il est expliqué dans l’algorithme, il vaut mieux norma-
(k)
liser x(k) = xx(k) dés qu’il est calculé avec l’inconvénient que le calcul est très
k k
lourd à la main.
Je vous donne des résultats caculés à l’aide de Scilab

k 0 1 2 :::: 5 6
v (k) 1 0:8944272 0:7808688 :::: 0:7100107 0:7080761
0 0:4472136 0:6246950 :::: 0:7041909 0:7061361

Ainsi la valeur approchée de la valeur propre est obtenue par la formule

' v (6)T Av (6) = 2:9999962:
Cette approximation est plus précise que celle caculée précedemment.
Chapitre 3

Interpolation numérique

3.1 Généralités
3.1.1 Position du problème
Soit une fonction f dé…nie sur un intervalle [a; b] de R et dont la valeur
numérique n’est connue qu’en (n + 1) points de [a; b] ; soit fx0 ; x1 ; ::::; xn g cet
ensemble de points tels que xi 2 [a; b] ; 8i 2 f0; :::::; ng et f (xi ) la valeur de la
fonction f en un point xi :
Principe L’interpolation est une méthode qui permet d’obtenir une valeur
approchée (x) de la fonction f (x) en un point x 2 [a; b] et ceci de telle sorte
que la condition suivante soit véri…ée :

(xi ) = f (xi ) 8i 2 f0; :::::; ng :

La fonction ainsi déterminée est appelée fonction d’interpolation de la fonction

f sur le maillage fx0 ; x1 ; ::::; xn g : En principe la fonction d’interpolation est une
combinaison linéaire de fonctions simples
– polynômes
– fonctions circulaires
– exponentielles
– fractions rationnelles

Remarque 3.1.1 On peut trouver un problème d’interpolation ayant un ca-

ractère plus général : En (n + 1) points xi 2 [a; b] on connaît les valeurs :
0
f (xi ); f (xi ); ::::; f ( i ) (xi ) la fonction d’interpolation est alors telle que :

39
3. Interpolation numérique 40

0 0 ( i)
(xi ) = f (xi ); (xi ) = f (xi ); :::; (xi ) = f ( i ) (xi ) c’est le problème de
l’interpolation d’Hermite.

3.2 Interpolation polynômiale

3.2.1 Problème
Déterminer un polynôme Pn (x) de degré n tel que :

Pn (xi ) = f (xi ); 8 i 2 f0; :::::; ng :

Question : Un tel polynôme existe-t-il ? Si oui, est-il unique ?

Solution théorique du problème
On choisit comme fonction de base :

0 (x) = 1; 1 (x) = x ; ::::; n (x) = xn :

La fonction d’interpolation devient alors le polynôme d’interpolation Pn (x) de

degré n tel que : Pn (x) = a0 +a1 x+: : :+an xn où a0 ; :::::; an sont des constantes
à déterminer. A partir de la condition

Pn (xi ) = f (xi ); i = 0; n;

les coe¢ cients a0 ; :::::; an sont solution du système linéaire :

X
n
ak xki = f (xi ) ; 8i 2 f0; :::::; ng (S)
k=0

On écrit le système (S) sous la forme matricielle suivante

0 10 1 0 1
1 x0 x20 ::: xn0 a0 f (x0 )
B CB C B C
B 1 x1 x21 ::: xn1 C B a1 C B f (x1 ) C
B CB C B C
B 1 x2 x22 ::: xn2 C B a2 C = B f (x2 ) C
B CB C B C
B CB .. C B .. C
@ ::: ::: ::: ::: ::: A@ . A @ . A
1 xn xnn ::: xnn an f (xn )
| {z }| {z } | {z }
A X B

Pour prouver que le système AX = B possède une solution unique X il su¢ t

d’établir que le déterminant de la matrice A est non nul.
3. Interpolation numérique 41

Lemme 3.2.1 Soit x0 ; x1 ; : : : ; xn (n + 1) points distincts alors

1 x0 x20 ::: xn0

1 x1 x21 ::: xn1 " #
Y
n Y1
n Y
n

n = 1 x2 x22 ::: xn2 = (xi xj ) = (xi xj )

i>j j=0 i=j+1
::: ::: ::: ::: ::: j=0
1 xn xnn ::: xnn

n est appelé le déterminant de Vandermonde.

Remarque 3.2.1 Notons que pour n = 1 on a

1 x0 Y 1

1 = = (xi xj ) = (x1 x0 ):
1 x1 i>j
j=0

Par un raisonnement par récurrence on peut montrer le lemme i.e le déter-

minant ( du type VanderMonde) du système (S) est :

Y
n

n = (xi xj )
i>j
j=0

et est par conséquent non nul si les points xi sont tous distincts.
Conclusion n 6= 0 car les points xi ; i = 0; n; sont distincts. Donc
le système (S) admet une solution unique (a0 a1 : : : an )t : Ainsi le polynôme
d’interpolation Pn de degré n de la fonction f existe et est unique si les points
d’interpolation sont tous distincts i.e : 9! Pn tel que Pn (xi ) = f (xi ); 8i = 0; n:

Exemple 3.2.1 Calculons le polynôme passant par les points (0; 1); (1; 2); (2; 9) et (3; 28):
Etant donné ces 4 points le polynôme recherché est tout au plus de degré 3. Ces
coe¢ cients ai sont solution de
0 10 1 0 1
1 0 0 0 a0 1
B CB C B C
B 1 1 1 1 C B a1 C B 2 C
B CB C B C
B1 2 4 8 C Ba C = B 9 C
@ A @ 2A @ A
1 3 9 27 a3 28

dont la solution (obtenue par décomposition LU) est (1 0 0 1)t : Le polynôme

recherché est donc P3 (x) = 1 + x3 :
3. Interpolation numérique 42

Remarque 3.2.2 En général le système linéaire permettant de déterminer les

coe¢ cients ak est mal conditionné. Il n’est donc pas indiqué de calculer un pô-
lynome d’interpolation Pn en résolvant un système linéaire. Cette méthode est
donc rarement utilisée. Il est par contre classique de l’exprimer sous les formes
suivantes :

– polynôme d’interpolation de Lagrange

– polynôme d’interpolation de Newton

3.2.2 Interpolation de LAGRANGE

Dé…nition 3.2.1 Soient x0 ; x1 ; : : : ; xn (n + 1) points distincts. Le polynôme de
Lagrange au points xi est un polynôme de degré n noté Li (x) et véri…ant la
condition suivante :
(
0 ; i 6= j
8i 2 0; n : Li (xj ) =
1 ; i=j

Cela signi…e que le polynôme Li (x) prend la valeur 1 en xi et s’annule en

tous les autres points de collocation.

Théorème 3.2.1 Pour tout i = 0; n on a

Y n
(x xj )
Li (x) =
i6=j
(xi xj )
j=0

Preuve. On a par dé…nition

Li (xj ) = 0; 8i 6= j :

Par coséquent, pour tout i = 0; n les racines de Li (x) sont : x0 ; x1 ; : : : ; xi 1 ; xi+1 ; : : : ; xn :

Le polynôme Li (x) s’écrit donc comme suit :

Y
n
Li (x) = K(x x0 ) : : : (x xi 1 )(x xi+1 ) : : : (x xn ) = k (x xj ) ;
i6=j
j=0

où k est une constante à determiner.

3. Interpolation numérique 43

Y
n
D’autre part Li (xi ) = 1; d’où k (xi xj ) = 1
i6=j
j=0
Y
n

(x xj )
i6=j
j=0
Y
n
(x xj )
1
)k= Y
n : Finalement, on trouve Li (x)0 i n = Y
n := (xi xj )
:
(xi xj ) (xi xj ) i6=j
j=0
i6=j i6=j
j=0 j=0

À titre d’illustration, on représente sur la …gure (3.1) les graphes sur l’inter-
valle [ 1; 1] des polynômes de Lagrange qu’on note li ; i = 0; 4 associés aux
noeuds 1; 21 ; 0; 21 et 1:

FIG.3.1-
3. Interpolation numérique 44

Revenons à la résolution du problème Pn (xi ) = f (xi ) qui consiste à construire

Pn (x) véri…ant les conditions ci-dessus. Ce polynôme est de la forme
X
n
Pn (x) = f (xi ) Li (x) :
i=0

Pn (x) est bien de degré inférieur ou égal à n: et

X
n
Pn (xj ) = f (xi ) Li (xj ) = f (xj ) Lj (xj ) = f (xj ) ; j = 0; 1; :::; n:
i=0

Ainsi, Le polynôme d’interpolation de Lagrange est dé…ni par

X
n Y n
(x xj )
Pn (x) = f (xi ) :
i=0 i6=j
(xi xj )
j=0

Exemple 3.2.2 On cherche à construire le polynôme d’interpolation de La-

grange qui en 1 vaut 8; en 0 vaut 3 et en 1 vaut 6:
On a
(x x1 ) (x x2 ) (x x0 ) (x x2 ) (x x0 ) (x x1 )
P2 (x) = f (x0 ) + f (x1 ) + f (x2 )
(x0 x1 ) (x0 x2 ) (x1 x0 ) (x1 x2 ) (x2 x0 ) (x2 x1 )
x(x 1) (x + 1)(x 1) (x + 1)x
= 8 +3 +6 = 4x2 x + 3:
2 1 2

Erreur d’interpolation

On admettra le théorème suivant sans démonstration

Théorème 3.2.2 Soient x0 ; x1 ; ::::; xn ; (n + 1) points distincts de [a; b]

et f 2 C n+1 ([a; b]) : Alors, pour chaque point x 2 [a; b] ; il existe 2 [a; b] tel
que :

f (n+1) ( ) Y
n
f (x) Pn (x) = en (x) = (x xj ) ; ou Pn (x) est le polynôme d’interpolation de f:
(n + 1)! j=0

Posons Mn+1 = maxx2[a;b] f (n+1) (x) alors

Mn+1 Y
n
jen (x)j (x xj ) ; x 2 [a; b] :
(n + 1)! j=0

en (x) est appelée erreur d’interpolation au point x:

3. Interpolation numérique 45

Exemple 3.2.3 Déterminons le polynôme d’interpolation de Lagrange de la

fonction
p
f (x) = x sur les points x0 = 100; x1 = 121 et x2 = 144.
p
Avec quelle précision peut-on calculer 115 ?
Calculons P2 (x) :
P2 (x) = f (x0 ) L0 (x) + f (x1 ) L1 (x) + f (x2 ) L2 (x)
= 10L0 (x) + 11L1 (x) + 12L2 (x)
où
(x 121) (x 144) 132 265 1 2
L0 (x) = = x+ x
(100 121) (100 144) 7 924 924
(x 100) (x 144) 4800 244 1 2
L1 (x) = = + x x
(121 100) (121 144) 161 483 483
(x 100) (x 121) 275 221 1 2
L2 (x) = = x+ x
(144 100) (144 121) 23 1012 1012
660 727 1
Donc P2 (x) = 161
+ 10626 x 10626 x2
M3
je2 (115)j 3!
j(115 100) (115 121) (115 144)j où
000 3 5
M3 = max f (x) = 10
[100;144] 8
000 5
(avec f (x) = 83 x 2 ) d’ou je2 (115)j 3
8:3!
10 5 (2610) 0; 163 10 2 :

3.2.3 Interpolation dans la base de Newton

Les di¤érences divisées de Newton

Dé…nition 3.2.2 Soit f une fonction dé…nie sur [a; b] et x0 ; x1 ; ::::; xn

(n + 1) points distincts de [a; b] :On dé…nit les di¤érences divisées de f aux
points x0 ; x1 ; ::::; xn par les relations de récurrence :
f [x0 ] = f (x0 ) ; 8i 2 0; n : f [xi ] = f (xi )

f [x0 ] f [x1 ] f (x0 ) f (x1 )

f [x0 ; x1 ] = =
x0 x1 x0 x1

f [x0 ; x1 ] f [x1 ; x2 ]
f [x0 ; x1 ; x2 ] =
x0 x2
:::::::::::::::::::::::::::::::::::::
f [x0 ; x1 ; ::; xn 1 ] f [x1 ; x2 ; ::; xn ]
f [x0 ; x1 ; ::; xn 1 ; xn ] = ;
x0 xn
3. Interpolation numérique 46

f [x0 ; x1 ; ::; xn 1 ; xn ] s’appelle la di¤érence divisée d’ordre n de la fonction f

aux points x0 ; x1 ; ::; xn 1 ; xn :

Construction de la table des di¤érences divisées

Nous nous contentons de construire une table s’arrêtant aux quatrièmes di¤é-
rences divisées, les autres s’obtiendraient de la même manière.

xi f (xi ) f [xi ; xi+1 ] f [xi ; xi+1 ; xi+2 ] f [xi ; xi+1 ; xi+2 ; xi+3 ] f [xi ; xi+1 ; xi+2 ; xi+3 ; xi+4 ]
x0 f [x0 ]
x1 f [x1 ] f [x0 ; x1 ]
x2 f [x2 ] f [x1 ; x2 ] f [x0 ; x1 ; x2 ]
x3 f [x3 ] f [x2 ; x3 ] f [x1 ; x2 ; x3 ] f [x0 ; x1 ; x2 ; x3 ]
x4 f [x4 ] f [x3 ; x4 ] f [x2 ; x3 ; x4 ] f [x1 ; x2 ; x3 ; x4 ] f [x0 ; x1 ; x2 ; x3 ; x4 ]

On verra dans ce qui suit que seules les valeurs encadrées, sur la diagonale
principale de la table, interviennent dans l’expression du polynôme d’interpola-
tion.

Formule de Newton

Théorème 3.2.3 Soit f : [a; b] ! R fonction dé…nie sur [a; b] et x0 ; ::; xn ; n +

1 points distincts de [a; b] : Le polynôme qui interpole f sur fx0 ; x1 ::; xn g est
donné par :

Pn (x) = f [x0 ] + f [x0 ; x1 ] (x x0 ) + :::: + f [x0 ; x1 ; ::; xn 1 ; xn ] (x x0 ) :::: (x xn 1 )

Xn Y
i 1
= f [x0 ] + f [x0 ; x1 ; ::; xi ] (x xj ) (3.1)
i=1 j=0

Cette forme est appelée forme de Newton d’interpolation pour les di¤érences
divisées.

Preuve. Voir [12; Th 5:3]

Exemple 3.2.4 Trouvons le polynôme d’interpolation de Newton pour la fonc-

tion f (x) = sin x en les 3 points xi = 2 i avec i = 0; :::; 2
3. Interpolation numérique 47

xi f (xi ) f [xi ; xi+1 ] f [xi ; xi+1 ; xi+2 ]

0 0
2
2
1
4
2
0 2

P2 (x) = f [x0 ] + f [x0 ; x1 ] (x x0 ) + f [x0 ; x1 ; x2 ] (x x0 ) (x x1 )

2 4
= 0 + (x 0) 2
(x 0)(x )
2
2 4
= x 2
x (x )
2
4
= 2
x (x ):

Remarque
1- La formule (3.1) est la forme la plus utilisée pour calculer le polynôme
d’interpolation, car elle nécessite le moins de calcul pour obtenir numériquement
ses coe¢ cients.
2- Dans la formule (3:1) posons x = xn : On obtient alors :

Pn (xn ) = f [x0 ]+f [x0 ; x1 ] (xn x0 )+::::+f [x0 ; x1 ; ::; xn 1 ; xn ] (xn x0 ) :::: (xn xn 1 ) :

Cette formule est connue sous le nom de " l’identité des di¤érences divisées".
3- Considérons encore la forme (3:1) du polynôme d’interpolation. Le poly-
nôme Pn (x) peut s’écrire aussi sous la forme suivante :

Pn (x) = a0 + a1 x + : : : + an xn : (3.2)

Par identi…cation des coe¢ cients correspondants à (3:1) et (3:2), on obtient :

f [x0 ; x1 ; ::; xn 1 ; xn ] = an :

D’où la di¤érence divisée d’ordre n d’un polynôme de degré n est constante.

Par conséquent, la di¤érence divisée d’ordre n + 1 d’un polynôme de degré n est
nulle.
4- Outre l’avantage d’être aisemment calculable, l’algorithme de Newton peut
être allongé ou réduit. Si l’on veut, par exemple, tenir compte de p points supplé-
mentaires, on a juste à les ajouter à la suite de la table des di¤érences divisées
et de compléter les di¤érences divisées. Si l’on veut au contraire, négliger les
q points, il su¢ t d’arrêter la table des di¤érences divisées aux (n + 1 q )points.
3. Interpolation numérique 48

Exemple 3.2.5 Reprenons l’exemple précédent et calculons le polynôme d’in-

terpolation de Newton de la même fonction en les 4 points xi = 2 i avec i = 0; :::; 3
i.e on a juste ajouté le point x = 32 Il su¢ t de calculer une di¤érence divisée en
plus, i.e ajouter une ligne au tableau :

xi f (xi ) f [xi ; xi+1 ] f [xi ; xi+1 ; xi+2 ] f [xi ; xi+1 ; xi+2 ; xi+3 ]
0 0
2
2
1
4
2
0 2

8
3 2
2
1 0 3 2

P3 (x) = f [x0 ] + f [x0 ; x1 ] (x x0 ) + f [x0 ; x1 ; x2 ] (x x0 ) (x x1 )

+f [x0 ; x1 ; x2 ; x3 ] (x x0 ) (x x1 ) (x x2 )
= P2 (x) + f [x0 ; x1 ; x2 ; x3 ] (x x0 ) (x x1 ) (x x2 )
4 8
= 2
x (x ) + x (x )(x )
3 2 2
8
= x (x2 3 x + 2 2 ):
3 3

3.3 Exercices
Exercice 1
Construire le polynôme d’interpolation de Lagrange de la fonction
f (x) = sin( x) sur les points x0 = 0; x1 = 16 ; x2 = 12 : En déduire une
valeur approchée de sin 5 et calculer l’erreur d’interpolation au point x = 15
Exercice 2
p
Avec quelle précision peut-on calculer 80 à l’aide de la formule de Lagrange
p
pour la fonction f (x) = x si on considère les points x0 = 64; x1 = 81;
et x2 = 100:
Exercice 3
On considère la fonction dé…nie sur l’intervalle [0; 0:4] par la table de valeurs :

x 0 0:1 0:2 0:3 0:4

f (x) 0 0:097554 0:190428 0:278919 0:363304

1) Calculer à l’aide de la méthode d’interpolation de Lagrange, une valeur

approchée de f (0:15):
3. Interpolation numérique 49

Exercice 4
1) Calculer les di¤érences divisées relatives au tableau suivant

i 1 2 3 4
xi 0 1 2 4
yi 1 1 2 5
2) Déterminer le polynôme de Newton relatif au tableau de la question 1.
Exercice 5
p
On considère la fonction f (x) = x
1) Calculer les di¤érences divisées de f liées aux points x0 = 7; x1 = 9;
x2 = 11; x3 = 13; x4 = 15:
2) Déterminer une approximation de f (8) à l’aide du polynôme de Newton.
3) Calculer l’erreur exacte en se basant respectivement sur le polynôme de
Newton p1 (x); p2 (x) et p3 (x):ensuite comparer les résultats avec l’estimation de
l’erreur donnée par la formule de Newton suivante
(en (x) ' f [x0 ; x1 ; ::; xn 1 ; xn ; xn+1 ] (x x0 )(x x1 )::::(x xn )).

3.3.1 Réponses
Exercice 1
P2 (x) = f (x0 )L0 (x) + f (x1 )L1 (x) + f (x2 )L2 (x); avec
1 1
(x )(x )
L0 (x) = 6
1
2
1 =1 8x + 12x2 ;
(0 6
)(0 2
)
1
(x 0)(x )
L1 (x) = 2
= 9x 18x2 ;
( 16 0)( 61 1
2
)
1
(x 0)(x )
L2 (x) = 6
= x + 6x2 :
( 12 0)( 21 1
6
)
D’où, P2 (x) = 72 x 3x2 ; x 2 0; 21 :
La valeur approchée de sin 5 : sin 5 ' P2 ( 15 ) (car sin 5
= f ( 15 )):par
7 3
conséquent sin 5 ' 10 25
= 0; 58:
L’erreur d’interpolation au point x = 15 : on sait que
1 M3 1 1 1
e2 x0 x1 x2
5 3! 5 5 5
1 3 1
où M3 = maxx2[0; 1 ] f (3) (x) = maxx2[0; 1 ] j 3
cos xj = 3
; e2 5 6 500
=
2 2
1; 0335:10 2 :
3. Interpolation numérique 50

Exercice 2
Pour qu’on puisse trouver la précision demandée, on utilise l’estimation de
l’erreur de l’interpolation suivante

Mn+1 Y
n
jf (x) Pn (x)j (x xi ) ;
(n + 1)! i=0

avec Mn+1 = maxx2[a;b] f (n+1) (x) :

Pour n = 2 : a = 64 et b = 100 on a
5
M3 = maxx2[64;100] f (3) (x) = f (3) (64) = 38 8 5 = 3 8 6
avec f (3) = 38 x 2

qui est une fonction décroissante pour tout x > 0:

Par conséquent pour x = 80 on a
p 8 6 160
80 P2 (80) j(100 80) (81 80) (64 80)j = ' 0; 00061 10 3 :
2 86
p
On peut alors calculer une valeur approchée de 80 qui est P2 (80) avec une
précision de 10 3 :
Exercice 4
Indication En faisant un tableau de di¤érences divisées comme vu au cours
le polynôme d’interpolation de Newton s’écrit

P (x) = f [x0 ] + (x x0 )f [x0 ; x1 ] + (x x0 )(x x1 )f [x0 ; x1 ; x2 ]

+(x x0 )(x x1 )(x x2 )f [x0 ; x1 ; x2 ; x3 ]
1 1
= 1 + (x 0)0 + (x 0)(x 1) + (x 0)(x 1)(x 2)( )
2 12
2 3 1 3
= 1 x + x2 x:
3 4 12
Chapitre 4

Intégration numérique

Principe
En général, il peut être di¢ cile, voir impossible de trouver la valeur exacte
d’une intégrale
Zb
f (x)dx :
a

Sauf dans le cas où on connaît explicitement une primitive de la fonction f ,

ou alors l’expression de f permettant par exemple une intégration par parties
ou un changement de variables approprié.
Zb
L’idée consiste à remplacer f par une approximation fe et calculer fe(x)dx
a
Zb Zb Zb
au lieu f (x)dx .i.e f (x)dx ' fe(x)dx cette formule s’appelle formule de
a a a
quadrature ou formule d’intégration numérique.
Si f est de classe C 0 sur [a; b] ; l’erreur de quadrature

Zb Zb
E(f ) = f (x)dx fe(x) dx
a a

satisfait
Zb Zb
E(f ) = f (x) fe(x) dx f (x) fe(x) dx (b a) max f (x) fe(x) :
x2[a;b]
a a

51
4. Intégration numérique 52

L’approximation fe doit être facilement intégrable, ce qui est le cas si, par
exemple fe est un polynôme.
Un choix naturel consiste à prendre fe le polynôme d’interpolation de La-
grange de f sur un ensemble de n + 1 noeuds distincts xi ; i = 0; n on aura
0 1
Zb Zb Xn Z b
Yn
x xj
I = f (x)dx ' Ie = fe(x)dx = @f (xi Li (x) dxA où Li (x) =
i=0
x
j=0 i
xj
a a a
i6=j

Il s’agit d’un cas particulier de la formule de quadrature suivante

X
n
Ie = i f (xi )
i=0

qui est une somme pondérée des valeurs de f aux points xi : on dit que ces points
sont les noeuds de la formule de quadrature et que les nombres i 2 R sont les
coe¢ cients ou encore les poids.

4.1 Méthode du rectangle

Formule de base La formule du rectangle est obtenue en remplaçant
f par une constante égale à la valeur de f en x0 = a; polynôme qui interpole f
en le point (a; f (a)) et donc de degré 0; ce qui donne

Zb Zb Zb
I= f (x)dx ' fe(x)dx = f (a) dx = (b a)f (a)
a a a

C’est à dire I ' (b a)f (a)

Erreur Si f 2 C 1 ([a; b]) alors il existe 2 ]a; b[ tel que
0
f (x) = f (a) + (x a)f ( ) et l’erreur de quadrature est majorée par

Zb Zb
jEj = f (x) fe(x) dx = f (x) f (a) dx
a a

Zb
0 (b a)2 0 (b a)2
= (x a) f ( )dx = f( ) M1 ;
2 2
a

0
avec M1 = maxx2[a;b] f (x) :
4. Intégration numérique 53

4.2 Méthode du trapèze

Formule de base La formule du trapèze est obtenue en remplaçant
f par le segment qui relie (a; f (a)) à (b; f (b)) polynôme qui interpole f en les
points (a; f (a)) et (b; f (b)) et donc de degré 1, on obtient une fonction de classe
C 2 sur [a; b] et P1 (x) le polynôme qui interpole f aux points fx0 = a; x1 = bg :
Alors
Zb Zb Zb
f (b) f (a) b a
I = f (x)dx ' fe(x)dx = (x a)+f (a) dx = (f (a) + f (b)) :
b a 2
a a a

Erreur Si f 2 C 2 ([a; b]) alors il existe 2 ]a; b[ tel que

00
f (x) fe(x) = f 2( ) (x a) (x b) et l’erreur de quadrature est majorée par
Zb Zb 00
f ( )
jEj = f (x) fe(x) dx = (x a)(x b) dx
2
a a
Zb
1 00 1 00 (b a)3
f ( ) (x a)(x b) dx = f ( )
2 2 6
a
3
(b a) 00
max f (x) :
12 [a;b]

4.3 Méthode de Simpson

Formule de base La formule de Simpson est obtenue en remplaçant f
par la parabole qui interpole f en a, en b et en a+b
2
donc un polynôme de degré
2, on a
Zb Zb
I = f (x)dx ' fe(x)dx
a a
Zb a+b a+b
(x 2
)(x b) (x a)(x b) a+b (x a)(x 2
)
= a+b
f (a) + a+b f + f (b)dx
a 2
(a b) ( 2 a)( a+b
2
b) 2 (b a)(b a+b
2
)
a
b a a+b
= f (a) + 4f ( ) + f (b) :
6 2
Erreur Si f 2 C 4 ([a; b]); l’erreur de quadrature est majorée par
(4) (b a)5
jEj max f (x) :
[a;b] 2880
4. Intégration numérique 54

Z1
Exemple 4.3.1 Approximons l’intégrale I = log (1 + x) dx par la méthode de
0
Simpson puis calculons l’erreur d’intégration. On a :
1 0 1
I = f (0) + 4f ( ) + f (1)
6 2
1 3
= 0 + 4 log( ) + log 2
6 2
= 0:38583

L’erreur d’intégration :
M4 6
jEj (b a)5 ; ou M4 = max f (4) (x) = =6
2880 [0;1] (1 + x)4
6
d0 ou jEj = 0; 2 10 2
2880
Z1
- La valeur exacte de I = log (1 + x) dx = 0:38629
0
- La valeur approchée de I par la méthode du trapèze :
b a 1
I' [f (a) + f (b)] = [0 + log 2] = 0:34657
2 2
L’erreur d’intégration :

a)3
(b 00 1
jEj M2 ; M2 = max f (x) = max 2 = 1
12 [0;1] [0;1] (1 + x)

1
) jEj = 0:8 10 1
12
Conclusion : La méthode de Simpson est plus précise que la méthode du trapèze.

4.4 Formule des Trapèzes composée

Soit fx0 ; x1 ; :::; xn g une partition de [a; b] ;
b a
où x0 = a; x1 = a + h; ::::; xn 1 = x0 + (n 1) h; xn = b et h = n
; on sait
que
Zb Zx1 Zx2 Zb n Zxi
X
I= f (x)dx = f (x)dx + f (x)dx + ::: + f (x)dx = f (x)dx:
a a x1 xn i=1 x
1 i 1
4. Intégration numérique 55

En appliquant la formule du trapèze à chaque intervalle [xi 1 ; xi ] ; i 2 1; n on

obtient :
Zxi
xi xi 1 b a
f (x)dx ' [f (xi 1 ) + f (xi )] ; avec xi xi 1 =h=
2 n
xi 1

c’est-à-dire
Zb
hX
n
f (x)dx ' [f (xi 1 ) + f (xi )]
2 i=1
a

ce qui donne

Zb
b a f (x0 ) f (xn )
f (x)dx ' + f (x1 ) + :::: + f (xn 1 ) + :
n 2 2
a

Cette formule est appelée "formule des Trapèzes composées":

Erreur
On sait que l’erreur d’intégration Ei dans l’intervalle [xi 1 ; xi ] est donnée
par :l’intervalle

(xi xi 1 )3 00
Ei = f ( i) ; i 2 [xi 1 ; xi ]
12
X
n
d’autre part : Ei = E i.e :
i=1

h3 X 00 h3 X 00
n n
h3 n M2 00
E = f ( i ) ) jEj f ( i) ; ou M2 = max f (x)
12 i=1 12 i=1 12 [a;b]

(b a)3
) jEj M2 :
12n2

Illustration de la règle du trapèze composée voir la …gure 4.1.

4. Intégration numérique 56

Fig 4.1 : trapèze composée à 4 sous intervalles sur [a; b]

La valeur approchée de l’intégrale I correspond à l’aire colorée.

Z1
Exemple 4.4.1 Approximons l’intégrale I = log (1 + x) dx par la méthode
0
du Trapèze composée avec n = 3 et estimons l’erreur d’intégration associée.

On a :
1
0 f (0) 1 2 f (1)
I ' +f +f +
3 2 3 3 2
1 1 2 1
) I' log 1 + + log 1 + + log 2
3 3 3 2
I ' 0; 38169:

L’erreur d’intégration :

(10)3 00
jEj M2 ; ou M2 = max f (x) = 1:
108 [0;1]
1
) jEj = 0; 92 10 2 :
108
4. Intégration numérique 57

4.5 Formule de Simpson composée

On décompose l’intervalle [a; b] en n sous intervalles de largeur h = b na avec
n 1: En introduisant les noeuds de quadrature xi = a + i h2 ; i = 0; 1; :::; 2n
(i ;e ; chaque sous-intervalle [x2i ; x2i+2 ] a pour largeur h et donc x2i+1 est son
point milieu), la formule de Simpson composée est
Zb n 1 Z2i+2
X
x
X
n 1
e x2i+2 x2i
I = f (x)dx = f (x)dx = (f (x2i ) + f (x2i+1 ) + f (x2i+2 ))
i=0 i=0
6
a x2i

h X
n 1
= (f (x2i ) + f (x2i+1 ) + f (x2i+2 ))
6 i=0
!
h X
n 1 X
n 1
= f (a) + f (b) + f (x2k ) + 4 f (x2i+1 )
6 i=1 i=0
!
h X
n 1 X
n 1
1
= f (a) + f (b) + f (a + ih) + 4 f a + (i + )h :
6 i=1 i=0
2
4
Erreur Si f 2 C ([a; b]) alors l’erreur de quadrature est majorée par
X
n 1
h5 b a 4
E[a;b] = E[x2i ;x2i+2 ] n max f (4) (x) = h max f (4) (x) :
i=0
2880 [a;b] 2880 [a;b]

Illustration de la règle de Simpson composée voir …gure 4.2 :

Fig 4.2 : Règle de Simpson composée à 4 sous-intervalles sur [a; b]

La valeur approchée de l’intégrale I correspond à l’aire colorée.
4. Intégration numérique 58

4.6 Exercices
Exercice 1
Z3
dx
Calculer à 0.01 près l’intégrale I = 1+x
par la méthode des trapèzes géné-
1
ralisée et comparer le résultat avec la valeur exacte.
Exercice 2
Trouver le nombre n de subdivisions nécessaires de l’intervalle d’intégration
[ ; ] ; pour évaluer à 0:5 10 3 près, grâce à la méthode de Simpson,
Z
l’intégrale cos x dx

Exercice 3
Calculer une valeur approchée de l’intégrale dé…nie
Z1
p
I= 1 + x dx
0
6
par la formule de Simpson à 10 près, puis comparer le résultat trouvé avec
la valeur exacte de l’intégrale.

4.6.1 Réponses
Exercice 1
On sait que
(b a)3 b a 2
jET j 2
M2 = h M2 0; 01
12n 12n2
0; 01 12
h2 ; avec M2 = max f " (x)
M2 (b a) [1;3]

1 2
f (x) = x+1 ; f " (x) = (1+x) 3 fonction décroissante son max est atteind pour la
2 0;01 12
valeur 1 et h 2 0;25
= 0; 24 =) h 0; 48:
On prendra h = 0; 4 par suite n = b h a = 0;4 2
= 5:
x0 = 1 ! f0 = 0; 5; x1 = 1; 4 ! f1 = 0; 416; x2 = 1; 8 ! f2 = 0; 357;
x3 = 2; 2 ! f3 = 0; 312; x4 = 2; 6 ! f4 = 0; 277;
et x5 ="3 ! f5 = 0; 25: #
X
n 1
f0 +fn
I = h 2
+ fi = 0; 4 0;5+0;25
2
+ 0; 416 + 0; 357 + 0; 312 + 0; 277 =
i=1
0; 694:
4. Intégration numérique 59

Z3
IT = 0; 694 0; 01. La valeur exacte dx
1+x
= [ln(1 + x)]31 = 0; 693:
1
Exercice 2
Le pas d’intégration est h = b na = 2
n
: D0 autre part l’erreur théorique de la
méthode de Simpson est donnée par
4
b a 4 (4) 2 2
E(h) = h f ( )= cos ; avec 2 [a; b] :
180 180 n

Par conséquent
4
2 2
jE(h)j :
180 n
16 4
Ainsi pour que jE(h)j 0; 5 10 3 , il su¢ t que n véri…e 90 n4
0; 5 10 3

donc n4 0:5 110 3 90 16 4 ainsi n 18; 6:

On prendra par exemple n = 20; car pour la méthode de Simpson le nombre
de subdivisions de [a; b] doit toujours être pair.
Chapitre 5

Résolution numérique des

équations di¤érentielles
ordinaires

5.1 Introduction
Les équations di¤érentielles, souvent obtenues lors de la modélisation des
phénomènes physiques, peuvent être diviser en deux grandes familles :
Les équations di¤érentielles ordinaires et les équations aux dérivées partielles.
Dans ce dernier cas les équations sont caractériséés par le fait que la variable
dépendante (ou les variables dépendantes) est (ou sont) fonction de plusieurs
variables indépendantes.
L’autre catégorie d’équations di¤érentielles sont les équations di¤érentielles
ordinaires qui sont caractérisées par le fait que la variable dépendante (ou les
variables dépendantes) ne dépend (ou ne dépendent) que d’une seule variable
indépendante et que sa dérivée soit alors, une dérivée totale.
Un exemple de ce genre d’équations est l’équation qui gouverne l’évolution,
en fonction du temps t; de la vitesse v d’une masse m dans le champ de gravité
terrestre g en chute amortie par un coe¢ cient d’amortissement k :
dv k
=g v;
dt m
ou bien encore, l’équation de la position angulaire d’un pendule, de longueur

60
5. Résolution numérique des équations di¤érentielles ordinaires 61

l, qui oscille dans un plan vertical :

d2 g
+ sin = 0:
dt2 l
Il existe d’autres modèles basés sur des équations di¤érentielles ordinaires qui
sont extrêmement courants tel que :
En cinétique chimique, dynamique des populations et en météorologie.
Quand ces équations sont linéaires ou qu’elles peuvent raisonnablement être
rendues linéaires, souvent des solutions analytiques peuvent être obtenues pour
ces équations.
Toutefois, quand ces équations sont complexes ou non linéaires, les méthodes
analytiques échouent et une approche numérique s’avère être une solution. L’ob-
jectif de ce chapitre est de présenter quelques méthiodes numériques pour la
résolution des équations di¤érentielles ordinaires.

5.2 Préliminaires
Soit une équation di¤érentielle de forme générale
0 00
F (y; y ; y ; :::; y (n) ; t) = 0; t 2 [t0 ; T ]: (5.1)

On considère deux types de problèmes :

5.2.1 Problèmes de conditions aux limites

On connaît alors :
8 8
>
> y(t0 ) = y0 >
> y(T ) = yT
>
> >
>
>
> >
>
>
> >
>
>
< y 0 (t0 ) = y00 >
< 0 0
y (T ) = yT
.. .. :
>
> . >
> .
>
> >
>
>
> .. >
> ..
>
> . >
> .
>
: (n) (n)
>
: (n)
y (t0 ) = y0 y (n) (T ) = yT

Ces problèmes se traitent comme des problèmes d’équations aux dérivées

partielles.
5. Résolution numérique des équations di¤érentielles ordinaires 62

5.2.2 Problèmes de conditions initiales (problème de Cau-

chy)
On connaît ici : 8
>
> y(t0 ) = y0
>
>
>
>
>
>
>
< 0 0
y (t0 ) = y0
..
>
> .
>
>
>
> ..
>
> .
>
: (n 1)
y (n 1)
(t0 ) = y0
En général, une équation de type (5.1), se ramène à un système d’équations
di¤érentielles du 1er ordre ; en e¤et
0 00 0 00
F (y; y ; y ; :::; y (n) ; t) = 0 () y (n) = '(y; y ; y ; :::; y (n 1)
; t))
Posons
u1 = y
et
0
u 1 = u2 ()
0 00
u 2 = u3 () u3 = y
.. .. ..
. . .
0
un 1 = un () un = y (n 1)

et
0 0
un = y (n) = '(y; y ; ; y (n 1)
; t) = '(u1 ; :::; un ; t);
d’où en convenant des notations suivantes :

u1 (t) u2
u2 (t) ..
.
U (t) = .. ; (t; U ) = ; U; 2 Rn :
. un
un (t) '(U; t)

On a donc ramené (5.1) à un système de la forme :

8
dU
>
< dt = (t; U ) n o
0 (n 1)
où U (t0 ) = y0 ; y0 ; : : : ; y0 (5.2)
>
:
U (t0 ) donné,
5. Résolution numérique des équations di¤érentielles ordinaires 63

5.2.3 Existence et unicité des solutions

Dans toute la suite du chapitre, on ne s’interessera qu’aux équations di¤é-
rentielles du premier ordre de conditions initiales du type :
(
y 0 (t) = dy
dt
= f (t; y(t)) ; y 2 R; t 2 [t0 ; T ]
(5.3)
y(t0 ) = yo
La variable indépendante t représente très souvent le temps.

Remarque 5.2.1 On a vu que toute équation di¤érentielle de conditions ini-

tiales d’ordre n, se ramène à un système d’équations di¤érentielles du premier
ordre et d’autre part les formules établies par une équation du premier ordre res-
tant valables pour un système d’équations di¤érentielles à condition de remplacer
les fonctions scalaires par des fonctions vectorielles.

Théorème(Cauchy-Lipshitz)
On suppose que la fonction f véri…e les conditions suivantes :
- f est continue par rapport à ses deux variables
-f (t; y(t)) est Lipschitzienne en y uniformément par rapport à t; c’est-à-dire :
9 L > 0 tel que jf (t; y) f (t; z)j L jy zj 8t 2 [t0 ; T ]; 8y; z 2 R:
Alors dans ce cas le problème (5.3) admet une solution unique passant par yo
pour t = t0 :

5.3 Principe d’intégration approchée par des mé-

thodes numériques
Soit à résoudre le problème di¤érentiel :
8 dy
>
< dt = f (t; y(t)) ; y 2 R; t 2 [t0 ; T ]
(P )
>
:
y(t0 ) = yo

1) On va subdiviser l’intervalle [t0 ; T ] en n parties égales par l’intermédiaire

de points tn , en général équidistants, qui déterminent le pas d’intégration h =
tn+1 tn :
5. Résolution numérique des équations di¤érentielles ordinaires 64

2) E¤ectuer le (n+1)ieme pas d’intégration consiste donc à passer de la valeur

approchée yn (approximation de la valeur exacte y(tn )) supposée calculée, à la
valeur approchée suivante
yn+1 (approximation de y(tn+1 )), (voir …gure 5.1).

FIG. 5.1

On obtiendra alors la solution approchée sous forme de table numérique.

On distingue deux classes d’algorithmes (méthodes) numériques :
a) Méthodes à un pas ( à pas séparés)
Ils permettent de calculer yn+1 en n’utilisant que la valeur yn :
b) Méthodes à pas multiples ( à pas liés)
Ils permettent d’obtenir yn+1 en utilisant les valeurs yn ; :::; yn k (k …xé) ;
pour initialiser ce type d’algorithme on e¤ectue k pas d’une méthode à pas
séparés.
5. Résolution numérique des équations di¤érentielles ordinaires 65

5.4 Méthodes numériques à un pas

Dé…nition 5.4.1 Les méthodes numériques à un pas sont Les méthodes de
résolution numérique qui peuvent s’écrire sous la forme :

yn+1 = yn + hn (tn ; yn ); 0 n N

où : [t0 ; T ] R ! R est une fonction que l’on supposera continue.

5.4.1 Méthode d’EULER

La méthode d’Euler est la méthode à un pas associée à la fonction

(t; y; h) = f (t; y):

Reprenons le problème (P ), le but est d’obtenir une approximation de la

solution en t = t1 = t0 + h: Avant d’é¤ectuer la première itération, il faut
déterminer dans quelle direction on doit avancer à partir du point (t0 ; y0 ) pour
obtenir le point (t1 ; y1 ), qui est une approximation du point (t1 ; y(t1 )). Nous
0
n’avons pas l’équation de la courbe y(t); mais nous en connaissons la pente y (t)
en t = t0 . En e¤et, l’équation di¤érentielle assure que :
0
y (t0 ) = f (t0 ; y(t0 )) = f (t0 ; y0 ):

On peut donc suivre la droite passant par (t0 ; y0 ) et de pente f (t0 ; y0 ) l’équation
de cette droite, notée d0 (t); est :

d0 (t) = y0 + f (t0 ; y0 )(t t0 )

En t = t1 ; on a :

d0 (t1 ) = y0 + f (t0 ; y0 )(t1 t0 ) = d0 (t) = y0 + hf (t0 ; y0 ) = y1

En d’autres termes, d0 (t1 ) est proche de la solution analytique y(t1 ), c’est à dire

y(t1 ) ' y1 = d0 (t1 ) = y0 + hf (t0 ; y0 ):

Si on souhaite faire une deuxième itération et obtenir une approximation de

y(t2 ), on peut refaire l’analyse précédente à partir du point (t1 ; y1 ). On remarque
cependant que la pente de la solution analytique en t = t1 est :
0
y (t1 ) = f (t1 ; y(t1 )):
5. Résolution numérique des équations di¤érentielles ordinaires 66

On ne connâit pas exactement y(t1 ), mais nous possédons l’approximation y1 de

y(t1 ). On doit alors utiliser l’expression :
0
y (t1 ) = f (t1 ; y(t1 ) ' f (t1 ; y1 )

et construire la droite :

d1 (t) = y1 + f (t1 ; y1 )(t t1 )

qui permettra d’estimer y(t2 ). On constate que l’erreur commise à la première

itération est réintroduite dans les calculs de la deuxième itération. On a alors :

y(t2 ) ' y2 = d1 (t2 ) = y1 + hf (t1 ; y1 ):

Ceci conduit à l’Algorithme d’Euler :

(
yn+1 = yn + hf (tn ; yn ); 0 n N
tn+1 = tn + h

Remarque 5.4.1 L’erreur introduite à la première itération se répercute sur

les calculs de la deuxième itération , ce qui signi…e que les erreurs se propagent
d’une itération à l’autre. Il en résulte que l’erreur :

jy(ti ) yi j

augmente légerement avec i:

Exemple 1 : soit l’équation di¤érentielle :

( 0
y (t) = y(t) + t + 1
y(0) = 1

On a donc t0 = 0 et y(0) = 1, et on pend un pas de temps h = 0:1. De plus, on

a:
f (t; y) = y + t + 1
En utilisant la méthode d’Euler on obtient successivement des approximations
de y(0:1), y(0:2), y(0:3)...., notées y1 ; y2 ; y3 ::::la première itération donne :

y1 = y0 + hf (t0 ; y0 ) = 1 + 0:1f (0; 1) = 1 + 0:1( 1 + 0 + 1) = 1

5. Résolution numérique des équations di¤érentielles ordinaires 67

la deuxième itération fonctionne de manière similaire :

y2 = y1 + hf (t1 ; y1 ) = 1 + 0:1f (0:1; 1) = 1 + 0:1( 1 + 0:1 + 1) = 1:01
on parvient à :
y3 = y2 + hf (t2 ; y2 ) = 1:01 + 0:1f (0:2; 1:01)
= 1:01 + 0:1( 1:01 + 0:2 + 1)
= 1:029
la tableau suivant rassemble les résultats des dix premières itérations.
La solution analytique, de cette équation di¤érentielle est : y(t) = exp( t) + t .

Le tableau suivant rassemble les résultats des 10 premières itérations et permet

de comparer les solutions numériques et analytiques et regarder la croissance de
l’érreur.
ti y(ti ) yi jy(ti ) yi j
0.0 1.000000 1.000000 0.000000
0.1 1.004837 1.000000 0.004837
0.2 1.018731 1.010000 0.008731
0.3 1.040818 1.029000 0.011818
0.4 1.070302 1.056100 0.014220
0.5 1.106531 1.090490 0.016041
0.6 1.148812 1.131441 0.017371
0.7 1.196585 1.178297 0.018288
0.8 1.249329 1.230467 0.018862
0.9 1.306570 1.287420 0.019150
1.0 1.367879 1.348678 0.019201

5.4.2 Méthode de Taylor

Le développement de Taylor autorise une généralisation immédiate de la mé-
thode d’Euler, qui permet d’obtenir des algorithmes dont l’erreur de troncature
locale est d’ordre plus élevé. Nous nous limitons cependant à la méthode de
Taylor de second ordre. On cherche, au temps t = tn , une approximation de la
solution en t = tn+1 . On a :
y(tn+1 ) = y(tn + h)
00
0 y (tn )h2
= y(tn ) + y (tn )h + + 0(h3 )
2
5. Résolution numérique des équations di¤érentielles ordinaires 68

En se servant de l’équation di¤érentielle (P), on trouve

0
f (tn ; y(tn ))h2
y(tn+1 ) = y(tn ) + f (tn ; y(tn ))h + + 0(h3 )
2
Remarque 5.4.2 Dans la relation précédente, on voit apparaître la dérivée de
la fonction f (t; y(t)) par rapport au temps. On sait que :

0 @f (t; y(t)) @f (t; y(t)) 0

f (t; y(t)) = + y (t)
@t @y
c’est à dire
0 @f (t; y(t)) @f (t; y(t)) 0
f (t; y(t)) = + f (t; y(t))
@t @y
on obtient donc

y(tn+1 ) = y(tn ) + hf (tn ; y(tn )) (5.4)

h2 @f (tn ; y(tn )) @f (tn ; y(tn ))
+ + f (tn ; y(tn ))
2 @t @y
+0(h3 ):

En négligeant les termes d’ordre supérieur ou égal à 3, on en arrive à poser

y(tn+1 ) ' y(tn ) + hf (tn ; y(tn )) (5.5)

h2 @f (tn ; y(tn )) @f (tn ; y(tn ))
+ + f (tn ; y(tn ))
2 @t @y

qui sera à la base de la méthode de Taylor. Ceci conduit à l’algorithme de Taylor

d’ordre 2

yn+1 = yn + hf (tn ; yn )
h2 @f (tn ; yn ) @f (tn ; yn )
+ ( + f (tn ; yn ))
2 @t @y
avec tn+1 = tn + h

Dans cet algorithme, on a remplacé la solution analytique y(tn ) par son approxi-
mation yn dans (5:5). On en conclut que les erreurs se propagent d’une itération
à une autre.

Exemple 5.4.1 Soit l’équation di¤érentielle déja résolu par la méthode d’Euler :
0
y (t) = y(t) + t + 1 y(0) = 1
5. Résolution numérique des équations di¤érentielles ordinaires 69

Dans ce cas :
f (t; y) = y(t) + t + 1
et
@f @f
= 1 et = 1
@t @y
On aura :
h2
yn+1 = yn + hf ( y(t) + t + 1) + (1 + ( 1)( y(t) + t + 1))
2
La première itération de la méthode de Taylor d’ordre 2 donne (avec h = 0:1)
(0:1)2
y1 = 1 + 0:1( 1 + 0 + 1) + (1 + ( 1)( 1 + 0 + 1)) = 1:005
2
Une deuxième itération donne
(0:1)2
y2 = 1:005 + 0:1( 1:005 + 0:1 + 1) (1 + ( 1)( 1:005 + 0:1 + 1))
2
= 1:019025

Les résultats sont donnés dans le tableau suivant

ti y(ti ) yi jy(ti ) yi j
0:0 1:000000 1:000000 0:000000
0:1 1:004837 1:005000 0:000163
0:2 1:018731 1:019025 0:000294
0:3 1:040818 1:041218 0:000400
0:4 1:070302 1:070802 0:000482
0:5 1:106531 1:107075 0:000544
0:6 1:148812 1:149404 0:000592
0:7 1:196585 1:197210 0:000625
0:8 1:249329 1:249975 0:000646
0:9 1:306570 1:307228 0:000658
1:0 1:367879 1:368541 0:000662
Remarque 5.4.3 On remarque que l’erreur est plus petite avec la méthode de
Taylor d’ordre 2 qu’avec la méthode d’Euler.

L’algorithme de Taylor d’ordre p est :

p
X hk
yn+1 = yn + f (k 1)
(tn ; yn )
k=1
k!
avec tn+1 = tn + h
5. Résolution numérique des équations di¤érentielles ordinaires 70

Remarque 5.4.4 En pratique, la méthode de Taylor pour p 2 est presque

(k)
unitilisable pour les raisons suivantes : Le calcul des quantités f est très com-
plexe et très couteux en temps machine. En plus la méthode suppose a priori que
f soit très régulière ; les erreurs risquent donc de ne pas pouvoir être contrôlées si
certaines dérivées de f présentent des discontinuités. Un moyen pour surmonter
cette di¢ culté : les méthodes de Range -kutta.

5.4.3 Méthode de Runge- kutta

Les méthodes de Range-kutta sont inspirées des méthodes de Taylor du même
ordre.

Méthode de Range - kutta d’ordre 2 :

On a vu que le développement de la méthode de Taylor passe par la relation

(4)
y(tn+1 ) = y(tn ) + hf (tn ; y(tn )) (5.6)
2
h @f (tn ; y(tn )) @f (tn ; y(tn ))
+ ( + f (tn ; y(tn ))) + 0(h3 ):
2 @t @y
Le but est de remplacer cette dernière relation par une expression équivalente
possédant le même ordre de précision (0(h3 )). Soit la forme :
y(tn+1 ) = y(tn ) + a1 hf (tn ; y(tn )) (5.7)
+a2 hf (tn + a3 h; y(tn ) + a4 h)
où on doit déterminer les paramètres a1 ; a2 ; a3 et a4 de telle sorte les expressions
(5:6) et (5:7) aient toutes deux une erreur en O(h3 ). On ne trouve par ailleurs
aucune dérivée partielle dans cette expression. Pour y arriver, on doit recourir
au développement de Taylor en deux variables autour du point (tn ; y(tn )). On a
ainsi :
@f (tn ; y(tn ))
f (tn + a3 h; y(tn ) + a4 h) = f (tn ; y(tn ) + a3 h
@t
@f (tn ; y(tn ))
+a4 h + 0(h2 )
@y
la relation (5:7) devient alors :
y(tn+1 ) = y(tn ) + (a1 + a2 )hf (tn ; y(tn )) (5.8)
@f (tn ; y(tn )) @f (tn ; y(tn ))
+a2 a3 h2 + a2 a4 h2 + 0(h3 )
@t @y
5. Résolution numérique des équations di¤érentielles ordinaires 71

On voit immédiatement que les expressions (5:2) et (5:4) sont du même ordre.
Pour déterminer les coe¢ cients ai , il su¢ t de comparer ces deux expressions
terme à terme
coe¢ cients respectifs de f (tn ; y(tn )) :

h = (a1 + a2 )h
@f (tn ;y(tn ))
coe¢ cients respectifs de @t
:

h2
= a2 a3 h2
2
@f (tn ;y(tn ))
coe¢ cients respectifs de @y
:

h2
f (tn ; y(tn )) = a2 a4 h2 :
2
On obtient ainsi un système non linéaire de 3 équations comprenant 4 incon-
nues : 8
>
< 1= (a1 + a2 )
1
= a2 a3 (5.9)
> 2
: f (tn ;y(tn ))
2
= a2 a4
le système (5.9) est insoluble car il y a moins d’équation que d’inconnues. D’où
les variantes de la méthode de Runge-kutta d’ordre 2. Voici deux variantes les
plus couramment utilisées.

Méthode d’Euler modi…ée

1
a1 = a2 = ;
2
a3 = 1 et a4 = f (tn ; y(tn ))

Ainsi on obtient l’algorithme de la méthode d’Euler modi…ée

yb = yn + hf (tn ; yn )
h
yn+1 = yn + f (tn ; yn ) + f (tn+1 ; yb)
2
tn+1 = tn + h
5. Résolution numérique des équations di¤érentielles ordinaires 72

Méthode du point milieu

1
a1 = 0; a2 = 1; a3 =
2
f (tn ; y(tn ))
et a4 =
2

En remplaçant ces valeurs des coe¢ cients ai dans l’équation (5.7) on obtient
l’algorithme suivant :

k1 = hf (tn ; yn )
h k1
yn+1 = yn + h(f (tn + ; yn + ))
2 2
tn+1 = tn + h:

Remarque 5.4.5 On constate que la fonction f (t; y) est évaluée au point milieu
de l’intervalle [tn ; tn+1 ] ; d’où le nom de la méthode.

5.4.4 Méthode de Runge - kutta d’ordre 4 :

L’algorithme de la méthode Runge - kutta d’ordre 4 est :

k1 = hf (tn ; yn )
h k1
k2 = hf (tn + ; yn + )
2 2
h k2
k3 = hf (tn + ; yn + )
2 2
k4 = hf (tn + h; yn + k3 )
1
yn+1 = yn + (k1 + 2k2 + 2k3 + k4 )
6
tn+1 = tn + h:

5.5 Exercices
Exercice 1
Soit le problème de Cauchy
( 0
y =y+t
(1)
y(0) = 1
5. Résolution numérique des équations di¤érentielles ordinaires 73

Approcher à 10 3 la solution du problème 1 en t = 1 à l’aide de la méthode

d’Euler en subdivisant l’intervalle [0; 1] en 10 parties égales.
Sachant que la solution exacte de (1) est y(t) = 2 exp(t) t 1:
Exercice 2
Soit le problème de Cauchy
( 0
y = y 2ty
(1)
y(0) = 1
p
dont la solution exacte est y = 2x + 1:
Approcher la solution de (1) en t = 0:2 en e¤ectuant les calculs avec 6
décimales à l’aide des méthodes d’Euler modi…é et de Runge-Kutta d’ordre 2.

5.5.1 Réponses
Exercice 1
Par l’algorithme d’Euler 0 n 9 et h = 0; 1 : on trouve y(1) ' 3; 187 et la
solution exacte de l’équation (1) est donnée par l’équation y(t) = 2 exp(t) t 1.
Ce qui donne y(1) ' 3; 437:L’approximation calculée est donc très grossière.
Exercice 2
p
La solution exacte étant y = 2x + 1: On considère l’intervalle [t0 ; t1 ] avec
t0 = 0; t1 = t0 + h = 0; 2 et h = 0; 2:
Méthode de Runge-Kutta d’ordre 2 : yb = y0 + h f (t0 ; y0 )
avec f (t0 ; y0 ) = y0 2ty00 = 1 donc yb = 1; 2:
f (t1 ; y0 + h f (t0 ; y0 )) = f (0:2; 1:2) = 1; 2 2 1;20;2 = 0; 866667:
y1 = y0 + h2 (f (t0 ; y0 ) + f (t1 ; yb)) = 1 + 0;2
2
(1 + 0; 866667) = 1; 1866667:
Ainsi y(0; 2) ' y1 = 1; 1866667. La valeur exacte étant
p
y(0; 2) = 1; 4 = 1; 183216:
L’erreur commise est jy(0; 2) y1 j = 0; 003451 < 0; 5 10 2 :
D’où y(0; 2) = 1; 19 0; 01; ainsi y1 approche y(0; 2) avec 3 chi¤res signi…catifs
exacts.
Chapitre 6

Résolution des équations non

linéaires

6.1 Introduction
Beaucoup de problèmes réels se ramènent à la résolution d’une équation de
la forme
f (x) = 0: (6.1)

Ce problème intervient notamment dans l’étude générale d’une fonction d’une

variable (réelle), pour laquelle des solutions exactes ne sont pas explicitement
connues. Même dans le cas d’une équation algébrique, nous rappellons qu’il
n’existe pas de méthode de résolution générale à partir du degré cinq.
Alors, on fait appel aux méthodes numériques de résolution qui sont, en géné-
ral, des méthodes itératives fournissant la solution sous forme d’approximations
successives.

6.2 Généralités sur les méthodes de résolution

Les méthodes numériques pour approcher une racine x^ de f consistent à :
1. localiser grossièrement la racine de f en procédant à l’étude du graphe
ou à des évaluations qui sont souvent de type graphique ; on note x0 cette
solution grossière ;
2. Ensuite, on procède à construire, à partir de x0 une suite x1 ; x2 ; x3 ; :::

74
6. Résolution des équations non linéaires 75

telle que
lim xn = x^ où f (^
x) = 0:
n!+1

On dit alors que la méthode est convergente.

Dé…nition 6.2.1 Méthode itérative à deux niveaux

On appelle méthode itérative à deux niveaux un procédé de calcul de la forme

xn+1 = (xn ); n = 0; 1; 2; ::::
dans lequel on part d’une valeur donné x0 pour calculer x1 ; puis à partir
de x1 on calcule x2 etc. La formule est dite formule de récurrence. Le procédé
est appelé convergent si la suite (xn ) est convergente. Il est bien évident qu’une
méthode itérative n’est utile que s’il y a convergence vers les valeurs cherchées.

Remarque 6.2.1 On peut envisager des méthodes itératives multi-niveaux,

comme par exemple les schémas à trois niveaux dans lesquels on part de deux
valeurs données x0 et x1 pour calculer x2 ; puis à l’aide de x1 et x2 on calcule
x3 etc.. ;

Ordre de convergence

Dé…nition 6.2.2 On dit qu’une méthode convergente est d’ordre p 1 s’il existe
une constante telle que : j^
x xk+1 j x xk jp ou encore
j^

jxk+1 x^j
lim = :
k!+1 jxk x^jp

Si p = 1 et 0 < < 1 la convergence est dite linéaire, si p = 2 la

convergence est dite quadratique..
Si p > 1 ou p = 1 et = 0 alors la convergence est dite super-linéaire et si
= 1 elle est dite sous-linéaire.

Remarque 6.2.2 Pour obtenir l’ordre d’un processus itératif de la forme

xn+1 = (xn )

on utilise un développement de TAYLOR au voisinage de x

xn+1 x = (xn ) x = (xn ) (x);

6. Résolution des équations non linéaires 76

(en e¤et f (x) = 0 peut être ramenée à la forme x (x) il su¢ t de poser
(x) = x + f (x)) or :
(xn x)2 " (xk x)k (k)
(xn ) = (x) + (xn x) (x) + (x) + : : : + (x) : : :
2! k!
donc
X
1
(x xn )i (i)
x xn+1 = (x)
i=1
i!
ainsi si :
(i)
(x) = 0; 8i = 1; 2; :::; (p 1)
le processus itératif est d’ordre p:

Remarque 6.2.3 La convergence d’un algorithme itératif de la forme

xn+1 = (xn )

est d’autant plus rapide que son ordre est grand.

Critère d’arrêt
Il est important de disposer d’un critère pratique pour interrompre le proces-
sus itératif d’approximation de x^ au bout d’un certain nombre …ni d’opérations.
Nous présentons quelques tests d’arrêt :

jxn xn 1 j < c "

jxn xn 1 j < " (1 + jxn j)

0
jxn xn 1 j < c " jxn j
jf (xn )j < "
0
Avec " est l’erreur maximale tolérée, c et c sont des constantes dé…nies par
l’utilisateur.

6.3 Séparation des racines

Dé…nition 6.3.1 Une racine x^ est dite séparable si l’on peut trouver des voi-
sinages de x qui ne contiennent pas d’autres racines que x^ . La racine x^ est
séparée si on peut trouver des nombres a et b tels que : a < x < b:

On utilise les méthodes suivantes pour séparer les racines d’une équation.
6. Résolution des équations non linéaires 77

6.3.1 Méthode graphique

On localise grossièrement les racines de l’équation f (x) = 0 en considérant
par exemple :
le graphe de y = f (x) et on cherche son intersection avec l’axe Ox:
les graphes de y1 = f1 (x) et y2 = f2 (x) où f (x) = f1 (x) f2 (x); puis,
on cherche les points d’intersection des graphes de f1 et f2 ; dont les abscisses
sont exactement les racines de l’équation f (x) = 0.

6.3.2 Méthodes algébriques

Soit une fonction continue f :[a; b] ! R. On utilise le corollaire du théorème
de Rolle :
Entre deux racines successives x0 et x1 de f 0 (x) il y a une seule racine de
f (x) = 0
si f (x0 ):f (x1 ) < 0.
Notons que cela nécessite la connaissance des racines de f 0 (x) = 0:

Utilisons le résultat(précédent :
f (x) = 0 a au moins une racine dans [a; b] ;
- Si f (a):f (b) < 0
f (x) = 0 a une seule dans [a; b] si f 0 (x) 6= 0 dans ]a; b[ :
(
f (x) = 0 n’a pas de racine dans [a; b] ;
- Si f (a):f (b) > 0
ou f (x) = 0 a un nombre pair de racines dans [a; b]

6.4 Approximation des racines

6.4.1 Méthode de dichotomie ( ou méthode de la bissec-
tion)
Principe

La méthode de dichotomie consiste à construire une suite de sous intervalles

emboités [an ; bn ]
contenant une racine de l’équation f (x) = 0; outil de base le théorème des
valeurs intermédiaires.
6. Résolution des équations non linéaires 78

La fonction f est continue sur l’intervalle [a; b] ; véri…e f (a):f (b) < 0,
alors, d’après le théorème des valeurs intermédiaires, il existe au moins un point
b 2 ]a; b[ tel que f (b
x x) = 0.
Le priincipe est le suivant : posons a1 = a; b1 = b; on note la valeur initiale
x1 = a1 +b
2
1
milieu de l’intervalle [a1 ; b1 ], et on évalue f en ce point.
b et le problème est résolu:
a) Si f (x1 ) = 0; le point x1 est la solution x
b 2 ]a1 ; x1 [ et on posera a2 = a1 et b2 = x1
b) Si f (a1 ):f (x1 ) < 0; alors on a x
b 2 ]x1; b1 [ et on posera a2 = x1 et b2 = b1
c) Si f (b1 ):f (x1 ) < 0; alors on a x
Dans les deux derniers cas, on réitère le processus en partant de la valeur :
x2 = a2 +b
2
2
et ainsi de suite..:
Ce procédé permet de construire de manière récurrente la suite
a1 + b 1 a2 + b 2 an + b n
x1 = ; x2 = ; ::::; xn = ; ::::::
2 2 2
b de l’équation f (x) = 0:
qui converge vers la racine x

Critère d’arrêt de la méthode

Avec la méthode de la dichotomie, les itérations s’achèvent à la m-ème étape

quand la longueur de l’intervalle contenant la solution est inférieure à une tolé-
rance " …xée.
Posons L0 = b a la longueur de l’intervalle à la première itération, claire-
L0
ment, après m itérations la longueur de l’intervalle est égale à Lm = m
2

L0 L0
jxm bj
x Lm ") m
" ) log log (") ) log (L0 ) log (2m ) log (")
2 2m
) log (2m ) log (L0 ) log (") ) m log (2) log (L0 ) log (")

L0
log
log (L0 ) log (") "
)m )m
log (2) log (2)
Notons que cette inégalité est générale elle ne dépend pas du choix de la
fonction f:
6. Résolution des équations non linéaires 79

Exemple 6.4.1 Si on désire calculer une approximation xk de x b avec n déci-

b a n
males exactes, il su¢ t de poser 2k 0:5 10 :Ce qui revient à ce qu’on aille
dans les itérations jusqu’à ce que k véri…e l’inégalité :

log 0:5b 10a n

k :
log 2

Exemple 6.4.2 Soit f (x) = 39 43x + 39x2 5x3 : On cherche a estimer

x 2 [1; 5] tel que f (x) = 0: (Voir …gure (6.1) ).

Figure 6.1
6. Résolution des équations non linéaires 80

Etude de convergence

Proposition 6.4.1 Soit f une fonction réelle continue sur [a; b], véri…ant

f (a):f (b) < 0

et soit x^ 2 ]a; b[ l’unique solution de l’équation f (x) = 0. Alors, la suite (xn )n2N
construite par la méthode de dichotomie converge vers x^ et on a l’estimation
b a
8n 2 N; j^
x xn j
2n

Ordre de convergence

La convergence de la méthode de dichotomie est linéaire, en e¤et

jxk+1 x^j 1
lim = < 1:
k!+1 jxk x^j 2
La méthode de dichotomie est simple, malheureusement, assez lente. Cepen-
dant, cette lenteur est compensée par le fait que l’on est assuré d’obtenir toujours
une suite qui converge vers une racine de f (x) = 0 dans [a; b]. On dit que c’est
une méthode globalement convergente.
On peut l’utiliser pour obtenir une approximation grossière ( mais raison-
nable) de la racine recherchée servant d’initialisation à une méthode d’ordre plus
élévé dont la convergence n’est que locale, comme la méthode de Newton-
Raphson (voir dans la section qui suit).

6.4.2 Méthode de Newton-Raphson

La méthode de Newton-Raphson, encore appelée méthode de la tangente ou
encore méthode de Newton. Consiste à construire une suite (xn )n2N ; telle que
l’itéré xn+1 est obtenue comme l’intersection de la tangente à la courbe de f (x)
au point (xn ; f (xn )) avec l’axe des abscisses.
Cette tangente à pour équation :

y = f 0 (xn ) (xn+1 xn ) + f (xn )

En posant y = 0 on obtient la formule de Newton

f (xn )
xn+1 = xn
f 0 (xn )
6. Résolution des équations non linéaires 81

En e¤et, soit x^ une racine recherchée et x0 une estimation à priori de x^;

supposons également que la fonction f est de classe C 2 au voisinage de x^; le
développement de Taylor d’ordre 1 donne
"
f ( )
f (^ 0
x) = f (x0 ) + f (x0 ) (^
x x0 ) + (^
x x0 )2 avec 2 ]^
x; x0 [
2
comme f (^
x) = 0 on obtient :
"
f ( )
0
0 = f (x0 ) + f (x0 ) (^
x x0 ) + (^
x x0 )2 (*)
2
En négligeant le reste :
f" ( )
x x0 )2
(^
2
L’équation (*) peut être approximée par :
f (x0 )
0 ' f (x0 ) + f 0 (x0 ) (^
x x0 ) ) x^ ' x0 = x1
f 0 (x0 )
En répétant le processus on obtient la formule de récurrence :
f (xn )
xn+1 = xn :
f 0 (xn )

Figure 6.2 : Construction des premiers itérés de la mérhode de Newton-Raphson

6. Résolution des équations non linéaires 82

On remarque que la méthode de Newton nécessite à chaque itération l’évalua-

tion des deux fonctions f et f 0 au point courant xn : Cet e¤ort est compensé par
le fait que cette méthode est deconvergence quadratique si la racine recherchée
est simple. En e¤et, on a :
f (x)
(x) = x
f 0 (x)
et
f (x):f " (x)
= = 0 en x
f 0 (x)
donc l’algorithme est d’ordre 2.

Remarque 6.4.1 Si la racine x est de multiplicité m, nous avons l’algorithme

suivant :
f (xn )
8n 2 N; xn+1 = xn m 0 :
f (xn )

La plus grande di¢ culté dans l’utilisation de la méthode de Newton-Raphson

réside dans le caractère local de sa convergence. Si l’initialisation x0 est trop éloi-
gnée de la racine, la méthode peut très bien diverger. Pour cette raison, il est
courant dans les applications de l’associer à une méthode d’encadrement, comme
la méthode de dichotomie, cette dernière permettant d’approcher, bien que len-
tement, la racine recherchée de manière à fournir une bonne initialisation pour
la méthode de Newton-Raphson.
En conclusion, la méthode de Newton-Raphson est la méthode de choix
en termes de vitesse de convergence, puisqu’elle converge de manière quadratique
pour un x0 assez proche de x: Par ailleurs, elle nécessite pour cela que la dérivée
de la fonction f puisse être évaluée en tout point donné. Si cela n’est pas le cas,
on utilisera la méthode de la sécante (voir la prochaine section), dont la vitesse
de convergence est moindre mais ne requiert pas que la dérivée de f existe.

6.4.3 Méthode de la sécante

La méthode de la sécante peut être considérée comme une quasi-méthode de
Newton-Raphson, dans laquelle on aurait remplacé la valeur de f 0 (xk ) par
le taux d’accroissement de f sur un petit intervalle. C’est l’une des méthodes
que l’on peut employer lorsque la dérivée de f est compliquée, voir impossible,
à calculer ou encore coûteuse à évaluer.
6. Résolution des équations non linéaires 83

Plus précisemment,

f (xn ) f (xn 1 )
f 0 (xn ) =
xn xn 1
et l’équation de la sécante aux points xn 1 et xn est :

(f (xn ) f (xn 1 ))
y = f (xn ) + (x xn )
(xn xn 1 )

Ainsi au lieu d’utiliser la tangente au point xn nous allons utiliser la sécante

joignant le point (xn 1 ; f (xn 1 )) au point (xn ; f (xn )) pour en déduire xn+1 .
Ce dernier est obtenu comme intersection de la sécante passant par les points
d’abscisse xn 1 et xn et de l’axe des abscisses, i.e. on cherche x solution du
système .
(
y = f (xn ) + (x xn ) (f (x(xnn) f (xn 1 ))
xn 1 )
y=0
ce qui donne
(xn xn 1 )
x = xn f (xn ):
(f (xn ) f (xn 1 ))
D’où l’algorithme correspondant est l’algorithme de la sécante
8
>
< x0 ; x1 donnés
Pour n = 0; :::
>
:
poser xn+1 = xn (f (x(xnn) xf n(xn1 ) 1 )) f (xn )

Critère d’arrêt

Le processus s’achève dès que l’on ait jxn+1 xn j "; " tolérance …xée.

Convergence

Théorème 6.4.1 Sioit f est une fonction de classe C 2 (I); I étant un voisinage
de x^, l’unique racine de f dans I: Alors, ilexiste un > 0 tel que si x0 ; x1 sont
choisis pdans l’intervalle ]^
x b à l’ordre
; x^ + [ alors la suite (xk ) converge vers x
(1+ 5)
p = 2 ' 1; 618 (le nombre d’or).

Nous donnos une brève preuve du théorème : nous aurons besoin, d’abord,
du lemme suivant :
6. Résolution des équations non linéaires 84

Lemme 6.4.1 Soit ak une suite de réels positifs telle que

ak+1 ak ak 1

Alors il existe c > 0 tel que

k
ak c ap0 :
p
1+ 5
où p = 2
est le nombre d’or.

Preuve Posons rk = ln(ak ): Nous avons

rk+1 rk + rk 1

Il est facile de prouver que rk est majorée par le terme général de la suite
fk+1 = fk + fk 1 ; f0 = r0 ; f1 = 1 (suite de Fibonacci), suite dont le terme
général vaut
1
fk = a pk + b k :
p
Le résultat en découle et d’où la preuve du théorème.

6.4.4 Méthode du point …xe

Principe

Cette méthode que nous allons étudier utilise le fait que le problème f (x) = 0
peut se ramener au problème équivalent g(x) x = 0 pour lequel on a le résultat
suivant

Théorème 6.4.2 Soit [a; b] un intervalle non vide de R et g une application

continue de [a; b] dans lui même. Alors, il existe un point x 2 [a; b] ;appelé point
…xe de la fonction g; véri…ant g (x ) = x :

Preuve Posons f (x) = g(x) x: On a alors f (a) = g(a) a 0 et

f (b) = g(b) b 0; puisque g(x) 2 [a; b] pour tout x 2 [a; b] : Par cosé-
quent, la fonction f; continue sur [a; b] ; est telle que f (a) f (b) 0: Le théorème
des valeurs intermédiaires assure alors l’existence d’un point x 2 [a; b] tel que
0 = f (x ) = g (x ) x :
6. Résolution des équations non linéaires 85

Dé…nition (Méthode de point …xe)

Supposons que x 2 R soit un point …xe de g: La méthode de point …xe
consiste en la construction d’une suite (xk )k2N dé…nie par récurrence comme suit
(
x0 donné,
xk+1 = g(xk ) 8k 2 N:

Naturellement une telle suite n’est pas forcemment convergente. Par contre,
si elle converge, c’est-à-dire si la suite xk a une limte que nous notons x ; et si g
est continue, alors cette limite est nécessairement un point …xe de g puisque

x = lim xk+1 = lim g(xk ) = g( lim xk ) = g(x ):

k!1 k!1 k!1

Etude de convergence

Théorème 6.4.3 Soit une fonction g : [a; b] ! R:On se donne x0 2 [a; b] et on

considère la suite xk+1 = g(xk ) pour k 0: Si les deux conditions suivantes sont
satisfaites :
1. condition de stabilité : g (x) 2 [a; b] pour tout x 2 [a; b]
2. condition de contraction striste : il existe K 2 [0; 1[ tel que

jg(x) g(y)j < K jx yj pour tout x; y 2 [a; b] ;

alors
- g est continue
- g a un et un seul point …xe x 2 [a; b]
- la suite xk+1 = g(xk ) converge vers x pour tout choix de x0 2 [a; b] :

Preuve
Continuité La condition de contraction stricte entraîne que g est continue
puisque, si on prend une suite (yk )k2N 2 [a; b] qui converge vers un élément
x de [a; b] ; alors nous avons jg(x) g(yn )j K jx yn j
et par suite limk!1 g(yk ) = g(x):
Existence Prouvons l’existence d’un point …xe de g:
La fonction f (x) = g(x) x est continue dans [a; b] et, grâce à la condition
de stabilité, on a f (a) = g(a) a 0 et f (b) = g(b) b 0: En appliquant le
théorème des valeurs intermédiaires, on en déduit que f a au moins une racine
dans [a; b] ; i.e.g a au moins un point …xe dans [a; b] :
6. Résolution des équations non linéaires 86

Unicité L’unicité du point …xe découle de la condition de contraction stricte.

En e¤et, si on avait deux points …xes distincts x1 et x2 ; alors

jx1 x2 j = jg(x1 ) g(x2 )j K jx1 x2 j < jx1 x2 j

ce qui est impossible.d’où x1 = x2 :

Convergence Prouvons que la suite xk converge vers l’unique point …xe x
quand k tend vers +1 pour toute donnée initiale x0 2 [a; b] : On a

0 jxk+1 x j = jg(xk ) g(x )j K jxk xj

où K < 1 est la constante de contraction. En itérant k + 1 fois cette relation on

obtient
jxk+1 x j K k+1 jx0 x j ;

i.e, pour tout k 0

jxk+1 x j
K k+1 :
jx0 x j
En passant à la limite quand k tend vers +1 on obtient jxk+1 x j tend vers 0:

Nous avons un résultat plus précis si g est C 1 :

:

Théorème 6.4.4 Soit g une fonction continue sur l’intervalle [a; b] qui véri…e

8x 2 [a; b] : g (x) 2 [a; b] , a g (x) b; 8x 2 [a; b]

Alors la fonction g admet un point …xe dans l’intervalle [a; b] :

Si de plus la fonction g est de classe C 1 ([a; b]) et qu’il existe une constante
2 ]0; 1[ tel que :
0
g (x) 8x 2 [a; b]

Alors g admet un point …xe unique dans [a; b] et la suite xk converge vers x pour
tout x0 2 [a; b] : De plus, nous avons

xk+1 x 0
lim = g (x ):
k!+1 xk x
6. Résolution des équations non linéaires 87

Preuve
La dernière hypothèse implique que g est contractante stricte, c’est une consé-
quence du théorème des accroissement …nis. En e¤et, pour tout x; y de [a; b] ; il
existe un c 2 [a; b] tel que
0
g(x) g(y) = g (c) (x y);

il en découle que
jg(x) g(y)j < jx yj :
Les hypothèses du théorème 6.4.3 sont satisfaites, nous sommes assurés de l’exis-
tence et de l’unicité du point …xe x de g ainsi que la convergence de la suite
xk :
Pour établir le résultat de vitesse de convergence, En vertu du théorème
des accroissement …nis il existe, pour tout entier naturel k, un réel ck tel que
xk < ck < x et on a
0
xk+1 x = g(xk ) g(x ) = g (ck ) (xk x ):

La suite xk converge vers x ainsi ck converge vers x et puis que g est C 1 ([a; b]);
il s’ensuit de la dernière inégalité
xk+1 x 0 0
lim = lim g (ck ) = g (x ):
k!+1 xk x k!+1

Critère d’arrêt

Soit g : [a; b] ! R satisfaisant aux hypothèses du théorème 6.4.4.

On désire avoir une estimation du nombre d’opérations nécessaires pour at-
teindre une précision donnée ":
Soit xk la suite des approximations dé…nie par

xk+1 = g (xk ) ; k = 0; 1; :::

Pour k; n 2 N; n > k :

jxk xn j = jg(xk 1 ) g(xn 1 )j jxk 1 xn 1 j

et par récurrence

2 k n
jxk xn j jxk 1 xn 1 j jxk 2 xn 2 j ::: jx0 xn k j jb aj :
6. Résolution des équations non linéaires 88

En passant à la limite sur n, on obtient :

n
jxk xj (b a):
Cette relation fournit donc une estimation de l’erreur commise à l’étape n du
calcul et un moyen de trouver le nombre d’opérations nécessaires pour atteindre
une précision donnée ".

n n "
jb aj ") ) n log ( ) log (") log (b a)
jb aj
log (") log (b a)
)n :
log ( )
Remarque 6.4.2 Pour évaluer l’erreur de la n-ième approximation xn on peut
utiliser
la formule : jx xn j jxn xn 1 j ".

6.5 Exercices
Exercice 1
A l’aide de la méthode de dichotomie, résoudre l’équation x3 + 2x 7=0
avec deux décimales exactes.
Exercice 2
Calculer la racine de l’équation x3 2x2 5 = 0
par la méthode de Newton avec une précision de 5 C.S.E.
Une fois la racine obtenue, calculer jen j et enen 1 :
Obtenir expérimentalement le taux de convergence de la méthode.

Exercice 3
Trouver la racine positive comprise dans l’intervalle ]1; 1:7[ de l’équation :
x4 2x 4=0
avec une précision de deux décimales exactes, en utilisant la méthode de la
sécante.
Exercice 4
Utiliser la méthode du point …xe pour trouver la racine de l’équation
x3 x 1=0
comprise dans l’intervalle ]1; 2[ ; avec deux décimales exactes.
6. Résolution des équations non linéaires 89

6.5.1 Réponses
Exercice 1
L’intervalle de séparation des racines exactes peut être trouver en construi-
sant les courbes représentatives des fonctions y1 = x3 et y2 = 2x + 7
(f (x) = y1 y2 ) et chercher ensuite le point d’intersection de ces deux courbes.
Vous trouvez que le seul point d’intersection se trouve dans ]1; 2[ vous pouvez
prendre a = 1 et b = 2:
f (1) = 4 < 0; f (2) = 5 > 0 d’après le théorème des valeurs intermédiaires
il y a existence de la solution de l’équation donnée dans ]1; 2[, on notera cette
b
solution x::
Calculons le nombre d’itérations pour avoir la précision " = 0:5 10 2 :
1
log b " a log
0:5 10 2
On sait que k log 2
= log 2
' 7; 6, on prendra k = 8:
1+2 3
Soit x1 = 2 = 2 = 1:5 avec f (x1 ) = 0:625 < 0; ainsi la racine cherchée se
trouve dans ]1:5; 2[ :
On prend x2 = 1:5+2 2
b cherchée
= 1:75 avec f (x2 ) = 1:313 > 0 d’où la racine x:
est dans ]1:5; 1:75[ :
On prend x3 = 1:5+1:75 2
= 1:625 avec f (x3 ) = 0:296 > 0 d’où la racine x: b
cherchée est dans ]1:5; 1:625[ :
On prend x4 = 1:5+1:625 2
= 1:5625 avec f (x4 ) = 0:176 < 0 d’où la racine x: b
cherchée est dans ]1:5625; 1:625[ :
On prend x5 = 1:562+1:625 2
= 1:594 avec f (x5 ) > 0 d’où la racine x:b cherchée
est dans ]1:5625; 1:594[ :
On prend x6 = 1:562+1:594 2
= 1:578 avec f (x6 ) > 0 d’où la racine x:b cherchée
est dans ]1:5625; 1:578[ :
On prend x7 = 1:562+1:578 2
= 1:57 avec f (x7 ) > 0 d’où la racine x: b cherchée
est dans ]1:5625; 1:57[ :
On prend x8 = 1:562+1:57 2
= 1:566 avec f (x8 ) < 0 d’où la racine x: b cherchée
est dans ]1:566; 1:57[ :
On voit que la racine cherchée calculée avec deux décimales exactes est
x = 1:57 [Link]
Exercice 2
On a f (0) = 5 et f (10) = 795 par le théorème des V.I la racine x: b de
l’équation donnée se trouve dans ]0; 10[ :
Réduisant l’intervalle obtenu, étant donné que f (3) = 4 > 0;
d’après le T.V.I 9 x: b 2 ]0; 3[ tel que f (x:)b = 0:
6. Résolution des équations non linéaires 90

Il faut noter que dans cet intervalle toutes les conditions de la méthode de
Newton sont réunies.
Nous pouvons adopter comme approximation initiale x0 = 3:
La suite est générée par la formule de Newton xn+1 = xn ff0(x n)
(xn )
:
Ainsi on obtient x1 = 2:73333:
La précision est de 5 C.S.E alors le dernier C.S.E est du rang -4 et donc
" = 0:5 10 4 ; puis on travaille avec 5 décimales on obtient
x1 = 2:73334; x2 = 2:69163; x3 = 2:69064; x4 = 2:69064;
on s’arrête là car jx4 x3 j = 7:43988 10 6 < ":Donc x: b = x4 :

en en
en e2n en 1 e2n 1

x0 0:30936 0:957 10 2 = =
1 2
x1 0:4269 10 0:18224:10 0:13799 0:446
x2 0:98 10 3 0:96 10 6 0:22996 10 1
0:5377
x3 0 0 0 0
x4 0 0 = =
b
x:j
On a en ! 0:98 10 3 ' 0 et e2en ! 0:5377, on a bien lim jxjxk+1 x:j 2 = C < +1:
n 1 k b
On conclut que la convergence est quadratique.
Bibliographie

[1] O. Axelsson [1994], Iterative Solution [Link] University Press,

Cambridge.
[2] J. Barager [1977], Introduction à l’analyse numérique, Hermann.
[3] N. Ben Salah, [Link] [2005], Introduction aux méthodes numériques ap-
pliquées et à leur [Link] de publication universitaire, Tunis.
[4] C. Brezinski, [Link]-Zaglia [2006], Méthodes numériques itératives. Al-
gèbre linéaire et non linéaire. Ellipses, Paris.
[5] G. Charet [1975], Cours d’analyse numérique, Sedes- Informatique.
[6] P.G. Ciarlet [1982], Introduction à l’analyse numérique et à l’optimisation,
Masson, Paris.
[7] Crouzeix, Mignot, Analyse numérique des équations di¤érentielles, Masson,
Paris.
[8] J.P. Demailly [1991], Analyse numérique et équations di¤érentielles. Presses
Universitaires de Grenoble.
[9] E. Durand [1972], Solutions numériques des équations algébriques, tomes 1
et 2, Masson, Paris.
[10] G. Faccanoni, Analyse numérique, Receuil d’exercices corrigés et aide mé-
moire, http// [Link]/[Link]
[11] [Link], N.H. Allal, Exercices corrigés en Analyse numérique élémentaire,
O¢ ce des publications Universitaires.[2003].
[12] A. Fortin [1995], Analyse numérique pour ingénieurs, éditions de l’école
polytechnique de Montréal.
[13] G. Hacques [1972], Algorithmes numériques, cours- tome 1- Collection U.
Colin.

91
BIBLIOGRAPHIE 92

[14] R. Herbin, Cours d’Analyse numérique, https ://[Link]-

[Link]/~herbin/[Link]
[15] F. Jedrzejzwski [2005], Introduction aux méthodes numériques, deuxième
édition, Springer-Verlag France Paris.
[16] G. Koep‡er, Notes de cours, Méthodes numériques, [Link]-
[Link]/~gk/MN_S6.
[17] M. Lakrib [2003], Cours d’analyse numérique, O¢ ce des publications uni-
versitaires, Alger.
[18] P. Lascaux, [Link]éodor [2005], Analyse numérique matricielle appliquée à
l’art de l’ingénieur. 2.Méthodes itératives. Dunod, Paris.
[19] G. Legendre [2016], Introduction à l’ana-
lyse numérique et au calcul scienti…que,
https ://[Link]/ legendre/enseignement/cours_ananum_dauphine.pdf
[20] Legras [1971], Initiation à l’analyse numérique, Dunod.
[21] J-L. Matret [2001], Aide mémoire Analyse Numérique,
http ://[Link]/ jlm/cours/analnum
[22] P. Pelletier, Techniques numériques appliquées au calcul scienti…que, Mas-
son.
[23] A. Quarteroni, [Link], [Link] [2000], Méthodes numériques pour le calcul
scienti…que. Programme en MATlAB. Springer-Verlag France Paris.
[24] P. Spiteri [2000], Introduction à l’analyse numérique, Polycopié INP Tou-
louse, ENSEEIHT.
[25] M. R. Temam [1970], Analyse numérique, Collection SUP le mathématicien.
Presses Universitaires de FranceUniversité Paris-Sud, Centre d’Orsay.