CHAPITRE 1

Introduction Générale à l’Analyse Nu-

mérique

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Définition de l’analyse numérique . . . . . . . . . 1
1.2 Analyse d’erreur . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 Erreurs de modélisation . . . . . . . . . . . . . . 2
1.2.2 Erreurs de représentation sur ordinateur . . . . . 3
1.2.3 Erreurs de troncature . . . . . . . . . . . . . . . . 3
1.3 Notion d’erreurs . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Chiffres significatifs . . . . . . . . . . . . . . . . . . . . . 5

1.1 INTRODUCTION

1.1.1 DÉFINITION DE L’ANALYSE NUMÉRIQUE

L’Analyse Numérique comprend deux mots : l’analyse qui fait référence aux mathématiques
et le mot numérique qui fait référence au traitement informatique. En d’autres termes c’est
l’élaboration de méthodes de calcul mathématiques adaptées au traitement par ordinateur. En
général le but de ces méthodes est la résolution de problèmes concrets qui se posent dans dif-
férentes disciplines : Physique, Économie, Mécanique, ....etc.
Un exemple représentatif de ces méthodes est la prévision météorologique : Les données col-
lectées par les satellites et les stations d’observations donnent un aperçu sur l’état actuel du

1 de 45
temps ; La simulation numérique permet à partir de cet état initial de prévoir le temps qui
fera les jours suivants. Or cette simulation est la mise en œuvre sur ordinateur de méthodes de
résolutions numériques des équations mathématiques de la mécanique de fluide.

1.2 ANALYSE D’ERREUR

Nous nous sommes habitués à résoudre d’une façon analytiques un certain nombre de pro-
blèmes mathématiques. C’est le cas notamment des techniques d’intégration et de résolution
d’équations algébriques ou différentielles. Pour des problèmes plus compliqués l’analyse numé-
rique fournit plusieurs techniques de résolutions qui résultent en différents algorithmes.
Un ordinateur ne peut fournir que des solutions approximatives. Ces approximations dépendent
des contraintes physiques (espace mémoire, ...) et du choix des algorithmes retenus par le pro-
grammeur. Ces algorithmes dépendent de certains paramètres qui influent sur la précision du
résultat. Ces algorithmes dépendent de certains paramètres qui influent sur la précision du
résultat. De plus, on utilise en cours de calcul des approximations plus ou moins précises. Par
exemple, on peut remplacer une dérivée par une différence finie de façon à transformer une
équation différentielle en une équation algébrique. Le résultat final et son ordre de précision
dépendent des choix que l’on fait. Une partie importante de l’analyse numérique consiste donc
à contenir les effets des erreurs ainsi introduites, qui proviennent de trois sources principales :
• Erreurs de modélisation ;
• Erreurs de représentation sur ordinateur ;
• Erreurs de troncature.

1.2.1 ERREURS DE MODÉLISATION

La première étape de la résolution d’un problème, et peut-être la plus délicate, consiste
à modéliser le phénomène observé. c’est l’étape de mathématisation du phénomène physique
auquel on s’intéresse. Cette étape consiste à faire ressortir les causes les plus déterminantes
du phénomène observé et à les mettre sous forme d’équations. Le plus souvent, on obtient des
équations différentielles ou aux dérivées partielles. L’effort de modélisation produit en général
des systèmes d’équations complexes qui comprennent un grand nombre de variables inconnues.
Pour réussir à les résoudre, il faut simplifier certaines composantes et négliger les moins impor-
tantes ce qui entraîne une erreur de modélisation. Par exemple, pour le problème du pendule
qui est connu depuis très longtemps, une brève analyse montre qu’il est raisonnable de négliger
la friction de l’air, car les vitesses de mouvement sont faibles. Il s’agit donc d’une erreur de

2 de 45
modélisation, qui paraît acceptable à cette étape de l’analyse.

1.2.2 ERREURS DE REPRÉSENTATION SUR ORDINATEUR

La deuxième catégorie d’erreurs est liée à l’utilisation de l’ordinateur. En effet, la représen-
tation sur ordinateur (généralement binaire) des nombres introduit souvent des erreurs. Même
infimes au départ, ces erreurs peuvent s’accumuler lorsqu’on effectue un très grand nombre
1
d’opérations. Par exemple, la fraction n’a pas de représentation binaire finie, pas plus qu’elle
3
ne possède de représentation décimale finie. On ne pourra donc pas représenter exactement
cette fraction, ce qui introduit une erreur. Ces erreurs se propagent au fil des calculs et peuvent
compromettre la précision des résultats.

1.2.3 ERREURS DE TRONCATURE

Les erreurs de troncature proviennent principalement de l’utilisation du développement de
Taylor, qui permet par exemple de remplacer une équation différentielle par une équation al-
gébrique. Cette erreur est générée lors du choix de l’ordre de précision de la discrétisation,
c’est-à-dire du terme à partir duquel nous négligeons le reste du développement limité.
Le développement de Taylor est le principal outil mathématique du numéricien. Il est primor-
dial d’en maîtriser l’énoncé et ses conséquences. Tout au long de ce cours, nous abordons des
méthodes de résolution qui comportent des erreurs de troncature plus ou moins importantes.
L’ordre d’une méthode dépend du nombre de termes utilisés dans les développements de Taylor
appropriés.

1.3 NOTION D’ERREURS

Definition 1: Erreur absolue

Soit x, un nombre, et x∗ , une approximation de ce nombre. L’erreur absolue est définie
par :
∆x = |x − x∗ | (1.1)

Definition 2: Erreur relative

Soit x, un nombre, et x∗ , une approximation de ce nombre. L’erreur relative est définie
par :
|x − x∗ | ∆x ∆x
Er (x∗ ) = = = ∗ (1.2)
|x| |x| |x |

3 de 45
 De plus, en multipliant par 100 %, on obtient l’erreur relative en pourcentage.

En pratique, il est difficile d’évaluer les erreurs absolue et relative, car on ne connaît géné-
ralement pas la valeur exacte de x et l’on n’a que x∗ . C’est pourquoi on utilise l’approximation
∆x
pour l’erreur relative. Dans le cas de quantités mesurées expérimentalement dont on ne
|x∗ |
connaît que la valeur approximative, on dispose souvent d’une borne supérieure pour l’erreur
absolue qui dépend de la précision des instruments de mesure utilisés. Cette borne est quand
même appelée erreur absolue, alors qu’en fait on a :

|x − x∗ | ≤ ∆x

ce qui peut également s’écrire :

x∗ − ∆x ≤ x ≤ x∗ + ∆x

et que l’on note parfois x = x∗ ± ∆x. On peut interpréter ce résultat en disant que l’on a estimé
la valeur exacte x à partir de x∗ avec une incertitude de ∆x de part et d’autre.
L’erreur absolue donne une mesure quantitative de l’erreur commise et l’erreur relative en
mesure l’importance. Par exemple, si l’on fait usage d’un chronomètre dont la précision est
de l’ordre du dixième de seconde, l’erreur absolue est bornée par 0,1 s. Mais est-ce une erreur
importante ? Dans le contexte d’un marathon d’une durée approximative de 2 h 20 min, l’erreur
relative liée au chronométrage est très faible :

0.1
= 0, 0000119
2 × 60 × 60 + 20 × 60

et ne devrait pas avoir de conséquence sur le classement des coureurs. Par contre, s’il s’agit d’une
course de 100 m d’une durée d’environ 10 s, l’erreur relative est beaucoup plus importante :

0.1
= 0, 01
10.0

soit 1 % du temps de course. Avec une telle erreur, on ne pourra vraisemblablement pas dis-
tinguer le premier du dernier coureur. Cela nous amène à parler de précision et de chiffres
significatifs au sens de la définition suivante.

4 de 45
 1.4 CHIFFRES SIGNIFICATIFS

Definition 3
Si l’erreur absolue vérifie :
∆x ≤ 0, 5 × 10m

alors le chiffre correspondant à la méme puissance de 10 est dit significatif et tous ceux
à sa gauche, correspondant aux puissances de 10 supérieures à m, le sont aussi. On
arrête le compte au dernier chiffre non nul. Il existe une exception à la règle. Si le chiffre
correspondant à la me puissance de 10 est nul ainsi que tous ceux à sa gauche, on dit
qu’il n’y a aucun chiffre significatif. Inversement, si un nombre est donné avec n chiffres
significatifs, on commence à compter à partir du premier chiffre non nul à gauche et
l’erreur absolue est inférieure à 0, 5 fois la puissance de 10 correspondant au dernier
chiffre significatif.

En pratique, on cherchera à déterminer une borne pour ∆x aussi petite que possible
et donc la valeur de m la plus petite possible.

22
On obtient une approximation de π(x = π) au moyen de la quantité ,
7
22
 
x∗ = = 3, 142857... . On en conclut que :
7

Exemple
22
∆x = x − = 0.00126... = 0.126 × 10−2
7

Puisque l’erreur absolue est plus petite que 0.5 × 10−2 , le chiffre des centièmes est
significatif et on a en tout 3 chiffres significatifs (3,14).

Une méthode numérique est dite stable si elle donne de bons résultats quelque soit
la nature de ses données. Les données initiales qui vont être traités par ordinateur
sont souvent tronquées. Il faut donc que les méthodes numériques utilisées soient
insensibles aux petites variations des données initiales ; on dit dans ce cas que la
méthode numérique est bien conditionnée. Une méthode est mal conditionnée si de
petites variations sur les données peuvent produire de grandes perturbations sur les
résultats obtenues. Il faut donc prendre en compte tous ces paramètres au cours de
l’élaboration des méthodes numérique. En Analyse Numérique, pour résoudre un

5 de 45
problème donné, souvent le problème qui se pose n’est pas de trouver une méthode
numérique mais la difficulté réside dans la démonstration que ce problème est bien
conditionné et que la méthode utilisée est stable.

6 de 45
 CHAPITRE 2
Systèmes d’équations linéaires

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Rappel mathématiques . . . . . . . . . . . . . . . . . . . 8
2.2.1 Les matrices . . . . . . . . . . . . . . . . . . . . . 8
2.2.2 Applications linéaires . . . . . . . . . . . . . . . . 11
2.3 Systèmes linéaires . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Résolution des systèmes d’équations linéaires . . . . . . . 14
2.4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . 14
2.4.2 Méthodes directes de résolution de systèmes li-
néaires . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.3 Méthodes itératives de résolution de systèmes li-
néaires . . . . . . . . . . . . . . . . . . . . . . . . 25

2.1 INTRODUCTION

Les progrès de l’informatique et de l’analyse numérique permettent d’aborder des problèmes

de taille prodigieuse. On résout couramment aujourd’hui des systèmes de plusieurs centaines de
milliers d’inconnues. On rencontre ces applications en mécanique des fluides et dans l’analyse
de structures complexes. On peut par exemple calculer l’écoulement de l’air autour d’un avion
ou l’écoulement de l’eau dans une turbine hydraulique complète. On peut également analy-
ser la résistance de la carlingue d’un avion à différentes contraintes extérieures et en vérifier
numériquement la solidité.

7 de 45
 Ces calculs complexes requièrent des méthodes évoluées. La modélisation de ces types des
problèmes génère généralement des systèmes d’équations linéaires ou non linéaires de taille
considérable, qu’on doit résoudre à l’aide de méthodes efficaces qui minimisent le temps de
calcul et l’espace mémoire requis.
Dans ce chapitre, nous allons aborder les principales méthodes de résolution des systèmes
linéaires à savoir les méthodes de Gauss.

Dans la plupart des cas, nous traitons des matrices non singulières ou inversibles,
c’est-à-dire dont la matrice inverse existe. Nous faisons donc de révision systéma-
tique de l’algèbre linéaire élémentaire.

2.2 RAPPEL MATHÉMATIQUES

2.2.1 LES MATRICES

— Une matrice A est un tableau rectangulaire d’éléments de K.

— Elle est dite de taille n × p si le tableau possède n lignes et p colonnes.
— Les nombres du tableau sont appelés les coefficients de A.
— Le coefficient situé à la i-ème ligne et à la j-ème colonne est noté ai,j .

Un tel tableau est représenté de la manière suivante :

 
a
 1,1
a1,2 . . . a1,j . . . a1,p 
 
a
 2,1 a2,2 . . . a2,j . . . a2,p 

 
 
... ... ... ... ... ...  
A= ou A = ai,j 1≤i≤n .
 
 
 ai,1 ai,2 . . . ai,j ... ai,p  1≤j≤p
 
 
 
... ... ... ... ... ...
 
 
an,1 an,2 . . . an,j . . . an,p

Opérations sur les matrices

1. Addition de matrices

8 de 45
 Definition 4
Soient A et B deux matrices ayant la même taille n × p. Leur somme C = A + B est la
matrice de taille n × p définie par

cij = aij + bij .

   
3 −2 0 5 
Si A=  et B= 
1 7 2 −1

alors

Exemple
 
3 3
A+B = .
3 6
 
−2
Par contre si B′ =   alors A + B′ n’est pas définie.
8

2. Produit d’une matrice par un scalaire

Definition 5
 
Le produit d’une matrice A = aij de Mn,p (K) par un scalaire α ∈ K est la matrice
 
αaij formée en multipliant chaque coefficient de A par α. Elle est notée α · A (ou
simplement αA).

Exemple

   
1 2 3 2 4 6
Si A=  et α=2 alors αA =  .
0 1 0 0 2 0

3. Multiplication de matrices : le produit AB de deux matrices A et B est défini si et

seulement si le nombre de colonnes de A est égal au nombre de lignes de B.

Definition 6: Produit de deux matrices

Soient A = (aij ) une matrice n × p et B = (bij ) une matrice p × q. Alors le produit
C = AB est une matrice n × q dont les coefficients cij sont définis par :
p
X
cij = aik bkj = ai1 b1j + ai2 b2j + · · · + aik bkj + · · · + aip bpj .
k=1

9 de 45
Matrices particulières

• Matrice carrée : Si n = p (même nombre de lignes que de colonnes). On note Mn (K) au

lieu de Mn,n (K).
• Matrice diagonale : si aij = 0. On note A = diag(a11 , . . . , ann .
• Matrice Identité : une matrice diagonale avec aii = 1. On note In = diag(1, . . . , ann .
• Matrices triangulaire inférieure : si ses éléments au-dessus de la diagonale sont nuls,
autrement dit :
i < j =⇒ aij = 0.

• Matrice triangulaire supérieure : on dit que A est triangulaire supérieure si ses

éléments en-dessous de la diagonale sont nuls, autrement dit :

i > j =⇒ aij = 0.

• Matrices symétriques : une matrice A de taille n × n est symétrique si A = AT , ou

encore si aij = aji pour tout i, j = 1, . . . , n.

Trace et déterminant

• Trace : dans le cas d’une matrice carrée de taille n × n, les éléments a11 , a22 , . . . , ann
sont appelés les éléments diagonaux, sa trace est la somme de ses termes diagonaux.
On a les propriétés suivantes :

T race(A + B) = T race(A) + T race(B)

T race(AB) = T race(A)T race(B)

• Déterminant : Le déterminant d’une matrice carrée (n, n) est noté det(A) ou encore

a11 a12 . . . a1n

. .. .. ..
det(A) = .. . . .
am1 am2 . . . amn

Il est défini par la formule suivante :

n
(−1)i+j aij det(Aij )
X
det(A) =
i=1

10 de 45
 X
det(A) = ϵσ aσ(1),1 aσ(2),2 aσ(3),3 . . . aσ(n),n (2.1)
σ∈Sn

où Sn est l’ensemble des permutations de 1, . . . , n et ϵσ est la signature de la permutation

σ (ϵσ ∈ {−1, 1}).
Les déterminants possèdent les propriétés suivantes :

1. det(I) = 1,

2. det(AT ) = det(A)

3. det(AB) = det(A)det(B)

4. pour tout scalaire α ∈ K, det(αA) = αn det(A)

5. le déterminant d’une matrice triangulaire (a fortiori diagonale) est égal au produit

de ses termes diagonaux.

Théorème 1
Soient A une matrice carrée, alors

A est inversible, si est seulement si det A ̸= 0.

De plus, si A est inversible, et C est sa comatrice, on peut calculer l’inverse de A par la

formule
1
A−1 = CT
det A

soit A ∈ Mn (K) on dit que A est inversible si

∃B ∈ Mn (K), AB = In ou ∃C ∈ Mn (K), CA = In

La matrice inverse de A est notée A−1 et vérifie AA−1 = A−1 A = In .

Nous avons les résultats suivants sur les inverses :

(AB)−1 = A−1 B −1

(AT )−1 = (A−1 )T

2.2.2 APPLICATIONS LINÉAIRES

Dans toute la suite, K désigne un corps. On peut penser à R ou C.
On considère deux espaces vectoriels E et F définis tous les deux sur K. On définit alors

11 de 45
l’application linéaire :

Definition 7
On appelle application linéaire f : E −→ F , une application possédant les propriétés
suivantes :
• f (x + y) = f (x) + f (y) ∀x ∈ E ∀y ∈ F
• f (λx) = λf (x) ∀x ∈ E ∀λ ∈ K

On pourra écrire toute application linéaire f : E −→ F comme suit :

   
 x1   a1,1 x1 a1,2 x2 . . . a1,j . . . a1,p xp 
 .   
 . 
 .   ... ... ... ... ... ... 
 
   
   
f x 
 j = a x
 i,1 1 ai,2 x2 . . . ai,j xj . . . ai,p xp 

 .. 
   
 
 .   ... ... ... ... ... ... 
   
   
xn an,1 x1 an,2 x2 . . . an,j xj . . . an,p xp

On appelle matrice associée à l’application f la matrice :

 
 a1,1 a1,2 . . . a1,j . . . a1,p 
 
... ... ... ... ... ...
 
   
 
A= a
 i,1 ai,2 . . . ai,j ... ai,p 
 ou A = ai,j 1≤i≤n .
  1≤j≤p
 
... ... ... ... ... ...
 
 
an,1 an,2 . . . an,j . . . an,p

On remarque que :
    
 x1   a1,1 a1,2 . . . a1,j . . . a1,p   x1 
 .    . 
 . 
 .   . . . . . . . . . . . . . . . . . .   .. 
  
    
    
f x 
 j = a
 i,1 ai,2 . . . ai,j ... ai,p 
  xj 
 
 ..    .. 
    

 . 
 . 
... ... ... ... ... ...  
  
    
xn an,1 an,2 . . . an,j . . . an,p xn

f : R3 −→ R3    
Exemple

x 2x + 3y + z 
   
f y  =  −x + 2z 
   
   
   
z x − y + 4z

12 de 45
 La matrice associée à l’application f est :
 
 2 3 1
 
A = −1 0 2
 
 
 
1 −1 4

2.3 SYSTÈMES LINÉAIRES

Definition 8
On appelle système linéaire S la donnée de n équations linéaires :

a11 x1 + a12 x2 + a13 x3 + · · · + a1n xn = b1








a21 x1 + a22 x2 + a23 x3 + · · · + a2n xn = b2



.. (2.2)
.








+ an2 x2 + an3 x3 + · · · + ann xn = bn

an1 x1


où x1 , x2 , . . .xp sont les p inconnues du système, k1 , k2 , . . .kn sont les n termes du second
membre ou constantes et les aij sont les n × p coefficients du système.

La résolution d’un système d’équations linéaires consiste à trouver un vecteur X = (x1 , x2 , ·, xn )T

solution de (S) 2.2. On peut utiliser la notation matricielle, qui est beaucoup plus pratique et
surtout plus compacte. On écrit alors le système précédent sous la forme :

AX = b

où A est la matrice :  
 a11 a12 a13 . . . a1n 
 
 a21 a22 a23 . . . a2n 
 
 .

.. .. . . .. 
 (2.3)
 .. . . . . 
 
 
an1 an2 an3 . . . ann

et b = (b1 , b2 , . . . , bn ) est le membre de droite.

Le problème 2.2 est un système de n équations et n inconnues. En pratique, la valeur de n varie
considérablement et peut s’élever jusqu’à plusieurs centaines de milliers. Notons que le coût de
la résolution croît rapidement avec n.

13 de 45
 2.4 RÉSOLUTION DES SYSTÈMES D’ÉQUATIONS LINÉAIRES

2.4.1 MOTIVATION

la solution de l’équation 2.12 peut s’écrire :

X = A−1 b

. et la discussion peut sembler close. Nous verrons cependant que le système linéaire 2.12 peut
devenir très grand, et donc difficile voire impossible à résoudre à la main. Pour comprendre
le problème, essayons de compter le nombre d’opérations nécessaires pour calculer la solution
en utilisant les formules de Cramer. Nous devons calculer n + 1 déterminants de taille n. On
sait que le déterminant de A est donné par la formule 2.1. L’ensemble Sn contient n! éléments
donc le calcul de det(A) se ramène à n! − 1 additions et n!(n − 1) multiplications soit au total
n!−1+n!(n−1) = nn!−1 opérations à virgule flottante. Comme nous avons n+1 déterminants
à calculer, le nombre d’opérations nécessaires pour résoudre le système à l’aide des formules de
Cramer est de (n + 1)(n × n! − 1) opérations à virgule flottante qui est équivalent quand la
dimension n tend vers l’infini à n(n + 1)!. Essayons d’évaluer cette quantité lorsque n = 100.
√
Pour ceci on peut utiliser la formule de Stirling n! ≃ nn+1/2 e−n 2π qui est très précise pour
évaluer la factorielle :

100 × 101! = 100 × 101 × 100!

√
≃ 100 × 101 × 100100,5 e−100 2π

≃ 9, 4 × 10161

Par exemple, avec ordinateur fonctionnant à 100 megaflops (flops = opérations à virgule flot-
tante par secondes - floating point operations per second), il faudrait environ 3 × 10146 années
pour résoudre notre système !

On en déduit donc qu’il est impossible d’utiliser les formules de Cramer pour des systèmes
de grande taille. En pratique, on utilise les méthodes directes pour des systèmes de dimension
peu élevée (n ≤ 100). Pour de très grands systèmes qui peuvent par exemple apparaître en
économie ou pour l’approximation d’équations aux dérivées partielles.

14 de 45
2.4.2 MÉTHODES DIRECTES DE RÉSOLUTION DE SYSTÈMES LINÉAIRES

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

Principe des méthodes directes

L’idée des méthodes directes est de se ramener à la résolution d’un (ou de deux) système(s)
triangulaire(s). Ces méthodes reviennent souvent à écrire une décomposition de la matrice
A = BC, où B et C ont des propriétés particulières. Ainsi résoudre AX = b équivaut à
résoudre B(Cx) = b ce qui est équivalent à résoudre :

By = b,

puis
Cx = y.

La plus célèbre des méthodes est appelée pivot de Gauss.

1. Méthode de Gauss

La méthode d’élimination de Gauss consiste à éliminer tous les termes sous la diagonale
de la matrice A. La stratégie de résolution est basée sur la question suivante : quelles
opérations sont permises sur les lignes du système 2.12 pour le ramener à un système
triangulaire ? Ou encore : pour ramener un système linéaire quelconque à un système
triangulaire, quels sont les coups permis, c’est-à-dire ceux qui ne changent pas la solution
du système de départ ? Pour transformer un système quelconque en système triangulaire,
il suffit d’utiliser trois opérations élémentaires sur les lignes de la matrice qui sont :

(a) échange de deux lignes.

(b) multiplication d’une ligne par un nombre non nul.

(c) addition d’un multiple d’une ligne à une autre ligne

L’algorithme de GAUSS

L’algorithme de résolution de GAUSS opère en deux phases à savoir la phase detriangularisation

(la transformation de la matrice en une matrice triangulaire supérieure) et la phase de

15 de 45
remontée (détermination de la solution du système linéaire). Pour la présenter, nous
allons prendre l’exemple suivant :

x + 2y + 2z = 2 L1







x + 3y − 2z = −1 L2




 3x + 5y + 8z = 8

L3

(a) La phase de triangularisation

On conserve la ligne L1 , qui sert de pivot pour éliminer l’inconnue x des autres
lignes ; pour cela, on retire L1 à L2 , et 3 fois L1 à L3 . On obtient :

x + 2y + 2z = 2 L1







 y − 4z = −3 L′2 ← L2 − L1



−y + 2z = 2 L′3 ← L3 − 3L1




On conserve alors la ligne L2 qui sert de pivot pour éliminer y de la troisième ligne ;
pour cela, remplace la ligne L3 par L3 + L2 . On trouve :

x + 2y + 2z = 2 L1







 y − 4z = −3 L′2



−2z = −1 L′′3 ← L3 + L2




(b) La phase de remontée

Ce dernier système, triangulaire, est facile à résoudre : la dernière ligne donne z, en

reportant, la deuxième ligne donne y, etc...
  
  1 

 x + 2y + 2z = 2,

 + 2 = 2,
x + 2y 
 x = 3,
2

 
 


 
 

 
1 
y − 4z = −3, =⇒ y = −4 = −3, =⇒ y = −1,





 2 


 
 1 
 1
−2z = −1 z = . z = .

 
 

2 2

16 de 45
Algorithm 1 Algorithme de GAUSS
Élimination
Require: matrice A, vecteur b
Ensure: Une matrice échelonnée A
for i=1 :n-1
for k=i+1 :n
Ak ← Ak − Ak,i /A,i × Ai ;
Endfor
Endfor

recherche de pivot non null

Soit le système linéaire suivant :


x1 + x2 + 2x3 + x4 = 2








2x1 + 2x2 + 5x3 + 3x4 = 4




x1 + 3x2 + 3x3 + 3x4 = −2









+ 4x3 + 5x4 = −2

 x1 + x2


dont la matrice augmentée est :

 
1 1 2 1 2
 
2 2 5 3 4





 (2.4)
1
 3 3 3 −2
 Exemple
 
1 a1 4 5 −2

En faisant les opérations indiquées (le pivot a11 est 1), on élimine les termes
non nuls sous la diagonale de la première colonne et l’on obtient :
 
1 1 2 1 2
 
0 0 1 1 0





 (2.5)
0
 2 1 2 −4

 
0 0 2 4 −4

Ici, la procédure est interrompue par le fait que le nouveau pivot serait 0 et
qu’il n’est pas possible d’éliminer les termes sous ce pivot. Toutefois, on peut
encore, parmi les opérations élémentaires, inter-changer deux lignes. Le seul

17 de 45
 choix possible dans cet exemple est d’intervertir la ligne 2 et la ligne 3. On se
rend immédiatement compte qu’il n’y a plus que des 0 sous le nouveau pivot
et que l’on peut passer à la colonne suivante :
 
1 1 2 1 2
 
0 2 1 2 −4
 



 (2.6)
0 0 1 1 0 
 
 
0 0 2 4 −4

En effectuant cette dernière opération, on obtient le système triangulaire :

 
1 1 2 1 2
 
0 2 1 2 −4
 



 (2.7)
0 0 1 1 0 
 
 
0 0 0 2 −4

La remontée triangulaire donne la solution :

X = [1 − 12 − 2]T

Complexité de l’algorithme

La méthode de Gauss nécessite un nombre total d’opérations élémentaires (additions,

soustractions, multiplications et divisions) à peu prés égal à 1/3n3 , soit une complexité
n3
de O( ) quand n devient grand (n ≥ 20).
3

Notion de pivotage

La méthode précédente s’applique lorsque tous les pivots ne sont pas nuls. Si tel n’est pas
le cas, on recherche dans les équations suivantes un coefficient non nul.
Afin de diminuer les risques d’incidents numériques, le pivot choisi est le plus grand en
valeur absolue. Lorsque ce pivot est trouvé, on effectue les permutations nécessaires pour
faire apparaître ce nouveau pivot à la place du pivot nul. Il existe deux stratégies de
pivotation (recherche de pivot non nul) : la pivotation partielle et la pivotation totale.
— Stratégie de pivot partiel : A la k éme étape, on choisira comme pivot l’élément akik

18 de 45
 vérifiant
|akik | = max |akip |
k≤p≤n

et on permute les lignes i et k.

— Stratégie de pivot total : A la k éme étape, on choisira comme pivot l’élément akij
vérifiant
|akij | = max |akip |
k≤p≤n

et on permute les lignes i et k puis les colonnes j et k.

en dépit que la stratégie de pivot total est toujours plus sure, l’expérience
montre que la stratégie de pivot partiel est suffisamment robuste dans la
plupart des cas. Mentionnons que la stratégie de pivot total est plus coûteuse.
De plus, il existe une classe importante de matrices ou‘ il n’est pas nécessaire
d’utiliser une stratégie de pivot, c’est le cas des matrices symétriques définies-
positives.

Exercice 1
Utiliser l’algorithme de GAUSS (avec pivotation partielle puis pivotation totale) pour
résoudre le système linéaire suivant (le mettre sous la forme Ax=b) :

2x + 4y − z + 5t = −10








 x + 2y + 7t = −13



x+y + 3z + t = 4








+ 2z + 4t = −5

2x + y


2. La décomposition LU

Principe de la méthode
Supposons qu’on doive, dans un programme donné, résoudre plusieurs fois le même sys-
tème linéaire AX = b, avec différents seconds membres b, mais la même matrice A. Si tous
les seconds membres b sont initialement connus, on peut alors effectuer simultanément
sur tous les seconds membres les manipulations intervenant dans l’élimination de Gauss.
Le plus souvent, dans la pratique, il s’agit de résoudre des systèmes avec des b calculés
au cours du processus et donc inconnus initialement. Il n’est alors pas raisonnable de
recommencer la triangulation de Gauss à chaque résolution : il faut conserver ce résultat

19 de 45
qui consiste en fait à factoriser la matrice A sous la forme A = LU ou L est une matrice
triangulaire inférieure et U une matrice triangulaire supérieure (L pour lower, U pour
upper) : on conservera en mémoire L et U . Par la suite, tout système

AX = b =⇒ LU X = b

sera remplacé par la résolutions des deux systèmes équivalents



 LY = b
(2.8)
U X =Y


Il s’agit alors de systèmes triangulaires qui se résolvent par une méthode de descente puis
de remontée, soit un nombre d’opérations égal à 2n2 .

Comment définir une décomposition LU ?

— Condition d’existence de la décomposition

Théorème 2
Soit A une matrice inversible. Une condition nécessaire et suffisante pour que A soit
décomposable en un produit LU est que tous ses mineurs principaux soient différents de
zéro.

— L’unicité de la décomposition

Proposition 1

La décomposition de GAUSS n’est pas unique, mais si l’on spécifie la diagonale de L ou

de U alors on a unicité de la décomposition.

Les deux choix les plus populairesde factorisation LU sont :

✓ Décomposition de Crout : elle consiste à imposer que la matrice U ait des 1
sur sa diagonale.
✓ Décomposition de Doolittle : elle consiste à imposer que la matrice L ait des
1 sur sa diagonale.

20 de 45
 Décomposition de Croûte

Pour obtenir cette décomposition, on doit déterminer les coefficients lij et uij des matrices
L et U de telle sorte que A= LU . En imposant que la diagonale de U soit composée de
1, on doit avoir :
    
 a11 a12 a13 . . . a1n   l11 0 0 ... 0  1 u12 u13 . . . u1n 
    
 a21

a22 a23 . . . a2n  l

l 0 ... 0  0 1 u23 . . . u2n 
 
 =  21 22

 .

.. .. . . ..   . .. .. .. .. 

. .. .. . . .. 
 (2.9)
 .. . . . .   .. . . . .   .. . . . . 
    
    
an1 an2 an3 . . . ann ln1 ln2 ln3 . . . lnn 0 0 0 ... 1

Il suffit de procéder de façon systématique par identification des coefficients du produit

L × U et la matrice A. Pour illustrer cette décomposition nous considérons l’exemple
suivant, une matrice de dimension 4 sur 4, le cas général étant similaire.
    
2 4 2 6  x 9
    
4 9 6 15
 y  23
    


 
  =  
  (2.10)
2 6 9 18 22
 z 
 
  
    
6 15 18 40 t 47

(a) Calculons le produit A = LU

   
2 4 2 6 l11 l11 u12 l11 u13 l11 u13 
   
4 9 6 15 l21 l21 u12 + l22 l21 u13 + l22 u23 l21 u14 + l22 u24
  
 



 = 



2 6 9 18 l l31 u12 + l32 l31 u13 + l32 u23 + l33 l31 u14 + l32 u24 + l33 u34
 31

  
   
6 15 18 40 l41 l41 u12 + l42 l41 u13 + l42 u23 + l43 l41 u14 + l42 u24 + l43 u34 + l44
(2.11)

(b) Identification des coefficients

i. Produit des lignes de L par la première colonne de U

l11 = 2 l21 = 4 l31 = 2 l41 = 6

ii. Produit de la première ligne de L par les colonnes de U


4




 l11 u12 = 4,


u12 = = 2,







 l11
 
2
l11 u13 = 2, =⇒ u13 = = 1,


   l11
6

 

l11 u13 =6
 

u13

 = = 3.
l11

21 de 45
 iii. Produit des lignes de L par la deuxième colonne de U . Les différents produits
donnent :
 
 
l21 u12 + l22 = 9, l22 = 9 − l21 u12 = 1,

 


 


 

 
l31 u12 + l32 = 6, =⇒ l32 = 6 − l31 u12 = 2,

 


 

 
l41 u12 + l42 = 15 l42 = 15 − l41 u12 = 3.

 


iv. Produit de la deuxième ligne de L par les colonnes de U :



6 − l21 u13
= = 2,

l + l22 u23 = 6, u23
21 u13

 
l22

=⇒
  15 − l21 u14
l
21 u14 + l22 u24 = 15 u24 = = 3.
 

l22

v. Produit des lignes de L par la troisième colonne de U

 
l31 u13 + l32 u23 + l33 = 9, l = 9 − l31 u13 − l32 u23 = 3,

 

33
=⇒
 
l
41 u13 + l42 u23 + l43 = 18 l = 18 − l41 u13 − l42 u23 = 6.
 
43

vi. Produit de la troisième ligne de L par la quatrième colonne de U

(
18 − l31 u14 − l32 u24

l31 u14 + l32 u24 + l33 u34 = 18 =⇒ u34 = =2
l33

vii. Produit de la quatrième ligne de L par la quatrième colonne de U

 
l41 u14 + l42 u24 + l43 u34 + l44 = 40 =⇒ l44 = 40 − l41 u14 − l42 u24 − l43 u34 = 1

Exercice 2
Soit le système :     
3 −1 2 x 12
    
1 2 3
 y 
= 11
    
 
  
    
2 −2 1 z 2

Décomposer le en un produit LU , puis résoudre le système.

22 de 45
Algorithm 2 Décomposition de Crout
Entrées : matrice A.
Sortis : matrices L et U satisfaisant A = LU
— Fixant la diagonale de la matrice U à 1 : Uii = 1 i = 1, . . . , n
— Première colonne de L : li1 ← ai1 i = 1, . . . , n
— Première ligne de U : u1i ← a1i /l11 i = 1, . . . , n
j−1
X
— Les valeurs de lij :lij ← aij − ljk ukj j≤i
k=1
i−1
1 X
— Les valeurs de uij :uij ← (aij − lik ukj ) j>i
lii k=1

L’algorithme précédent ne fonctionne que si les pivots lii sont tous non nuls.
Ce n’est pas toujours le cas et il est possible qu’il faille permuter deux lignes
pour éviter cette situation, tout comme pour l’élimination de Gauss. Le co-
efficient lii est encore appelé pivot.

Résoudre ce système avec la décomposition de Crout.

Exemple
    
0 2 1 x 1
    
1 0 0
 y 
= 2
    
 
  
    
3 0 1 z 3

3. Factorisation de Cholesky

Théorème 3
Si A est une matrice symétrique réelle définie positive, il existe au-moins une matrice B
triangulaire inférieure telle que A = LLT . De plus, on peut imposer que les coefficients
diagonaux de B soient strictement positifs et alors la décomposition est unique.

Dans le cas où la matrice est symétrique, on peut diminuer l’espace mémoire nécessaire à
la résolution d’un système linéaire. la factorisation de Cholesky qui consiste à décomposer
A sous la forme LLT où L est encore ici une matrice triangulaire inférieure. La matrice
triangulaire supérieure U de la décomposition LU est ainsi remplacée par la matrice
transposée de L, réduisant ainsi l’espace mémoire nécessaire. Si on souhaite résoudre un
système linéaire, on procède encore ici en deux étapes. On résout d’abord le système
triangulaire inférieur LY = b et ensuite le système triangulaire supérieur LT X = Y . On
refait donc le même raisonnement qui nous a mené à la décomposition LU générale mais

23 de 45
cette fois en considérant une matrice A symétrique :
    
 a11 a12 a13 . . . a1n   l11 0 0 ... 0  l11 l21 l31 . . . ln1 
    
 a21

a22 a23 . . . a2n  l

l 0 ... 0   0 l22 l32 . . . ln2 
 
 =  21 22

 . .. .. . . . (2.12)
 ... .. .. ..   . .. .. . . .. 
  
 .. ..
. ..  .
 
 . .   . . . .  .
  . . . . 

    
an1 an2 an3 . . . ann ln1 ln2 ln3 . . . lnn 0 0 0 . . . lnn

Pour une matrice 3 × 3, en faisant le produit, on trouve facilement :

√
 
2 l11 = a11 ,
 
a11 = l11 ,

 


 

 
a21

 

a21 = l21 l11 , =⇒ 21 l = ,





 l11

 
 a31
a31 = l31 l11 l31 =

 

l11

qui complète le calcul de la première colonne de L. Notons que le pivot l11 doit être non
nul. On poursuit avec le calcul de la deuxième colonne :

 
q
2

2 2 l22 = a22 − l21 ,
a = l21 + l22 ,

 

22
=⇒
  a32 − l31 l21
a

32 = l31 l21 + l32 l22 l32

 =
l22
q
2 2 2 2 2
Enfin : a33 = l31 + l32 + l33 qui entraîne que l33 = a33 − l31 − l32 .

Effectuer la factorisation de Cholesky de la matrice suivant :

 
Exemple

1 1 1 1
 
1 5 5 5


 
 
1 5 14 14
 
 
1 5 14 15

Dans le cas général, on obtient donc l’algorithme :

24 de 45
Algorithm 3 Factorisation de Cholesky
Entrées : matrice symétrique A.
Sortis : matrice L qui satisfait A = LLT
√
— Première Premier pivot l11 : li1 ← a11 i = 1, . . . , n
— Première colonne de L : li1 ← ai1 /l
v11 i = 2, . . . , n
u
u k−1
X
— Terme diagonal (pivot) lii : lii ← takk − a2kj k = 2, . . . , n
u
j=1
k−1
éme 1 X
— k colonne de L : lik ← (aik − lij lkj ) i>k
lkk j=1

2.4.3 MÉTHODES ITÉRATIVES DE RÉSOLUTION DE SYSTÈMES LI-

NÉAIRES
Dans la section précédente, nous avons vu les méthodes appelées méthodes directes qui four-
nissent des solutions exactes (aux erreurs d’arrondi près). Cependant, elles ont l’inconvénient de
nécessiter une assez grande place mémoire car elles nécessitent le stockage de toute la matrice
en mémoire vive. Si la matrice est pleine, c.à.d. si la plupart des coefficients de la matrice sont
non nuls et qu’elle est trop grosse pour la mémoire vive de l’ordinateur dont on dispose, ces
méthodes peuvent être lourds pour de très grandes matrices. En revanche,les méthodes itéra-
tives permettent de ne placer en mémoire que les coefficients non nuls d’une matrice de grande
taille. Cela est particulièrement important avec les matrices creuses, dont une grande partie
des coefficients sont nuls. Les méthodes itératives possèdent donc des avantages suffisamment
importants pour justifier une recherche active dans ce domaine. Cependant, contrairement à
la décomposition LU, le succès n’est pas assuré quelle que soit la matrice A pour laquelle on
souhaite résoudre un système linéaire. La convergence des méthodes itératives n’est réalisée que
dans certaines conditions que nous préciserons. les méthodes itératives, lorsqu’elles convergent,
ne deviennent vraiment avantageuses que pour les systèmes linéaires de très grande taille.

Definition 10: Méthode itérative

Une méthode de résolution de système linéaire est dite itérative si la résolution du système
peut être obtenue en construisant une suite (xk )k∈N censée converger vers x∗ solution du
système.

Definition 11: Méthode itérative convergente

On dit qu’une méthode itérative est convergente si pour tout choix initial x(0) ∈ Rn , on
a : (xk ) −→ x lorsque k −→ ∞

25 de 45
Principe des méthodes itératives :

Le principe de base de de telles méthodes est d’engendrer une suite de vecteurs xk (les itérés)
convergente vers la solution x∗ du système linéaire Ax = b. On veut que cette suite soit simple
à calculer. Une idée naturelle est de travailler avec une matrice P inversible qui soit “proche” de
A, mais plus facile que A à inverser. On appelle matrice de pré-conditionnement cette matrice
P. On écrit alors : A = P − (P − A).
Posons : P − A = N , Nous avons alors :

AX = b ⇐⇒ (P − N )X = b ⇐⇒ P X = N X + b

Sous cette forme, et à partir d’un choix initial X (0) donné nous cherchons à approcher la solution
du problème par la méthode itérative suivante :

X (0) ∈ Rn


(2.13)
P X k+1 = N Xk + b



A chaque itération, nous devons résoudre un système linéaire de matrice P pour calculer X k+1
en fonction de X k .

Montrons que cette méthode, si elle converge, approche bien la solution du système
AX = b :
Soit X ∗ solution de AX = b,
On a P X k+1 = N X k + b pour tout k ∈ N, donc à la limite : X k −→ X ∞ lorsque
k −→ ∞, alors,

P X∞ = N X∞ + b

P X ∞ = (P − A)X ∞ + b

AX ∞ = b

Et comme A est inversible, par unicité de la solution nous avons X ∞ = X ∗

26 de 45
 Definition 12: Erreur d’approximation

L’erreur d’approximation à la kième étape ek s’écrit :

ek = X k − X

où : X k est construit par 2.16 et X = A−1 b

La méthode converge si et seulement si lim ek = 0

k−→∞
Pour tout k ∈ N :

e(k+1) = X k+1 − X

Nous choisirons P inversible dans cette décomposition. De telle sorte que :

P X k+1 = N X k + b

équivalent à
X k+1 = P −1 N X k + P −1 b

de même
AX = b = P X − N X

qui est équivalent à

X = P −1 N X + P −1 b

Nous pouvons en déduire que

e(k+1) = P −1 N e(k)

Par une récurrence immédiate nous obtenons

ek = (P −1 N )k e(0)

Théorème 4: Convergence des méthodes itératives

Soit A, P ∈ Mn (K) deux matrices inversibles. La méthode itérative associée à la décom-

position A = P − N converge si et seulement si le rayon spectral ρ(P −1 N ) < 1 (max (des
valeurs propres de la matrices <1).

Tout le problème consiste dès lors à bien choisir la décomposition P − N . Les trois décom-

27 de 45
positions proposées ci-dessous sont deux méthodes les plus connues. Soit A une matrice d’ordre
n telle que aii ̸= 0. On décompose A sous la forme : A = D − E − F avec :

1. D : la diagonale de A

2. −E : la partie inférieure stricte

3. −F : la partie supérieure stricte.

1. Méthode de Jacobi :
Dans cette méthode, nous prenons : P = D, N = E + F La méthode de Jacobi s’exprime
donc par la suite : 
X (0) ∈ Rn donné


(2.14)
DX k+1 k

= (E + F )X + b k ≥ 1


Il est facile à chaque itération de calculer X k+1 en fonction de X k car la matrice diagonale
D est inversible.
n
xk+1 aij xkj + bi )/aii
X
i = (−
j=1,j̸=i

Sous forme matricielle :

X k+1 = D−1 (E + F )X k + D−1 b

On appelle la matrice J = D−1 (E + F ) la matrice d’itération de Jacobi. La méthode de

Jacobi converge ssi ρ(J) < 1.

2. Méthode de Gauss-Seidel :
Dans cette méthode, nous prenons ; P = D − E, N = F La méthode de Gauss-Seidel
s’exprime donc par la suite :

X (0) ∈ Rn donné


(2.15)
(D − E)X k+1 k

= FX + b k ≥ 1


D − E est une matrice triangulaire inférieure et pour calculer xk+1 en fonction de xk , on

applique l’algorithme de descente suivant : pour i = 1, n

i−1 n
xk+1 aij xk+1 aij xkj + bi )/aii
X X
i = (− j −
j=1,j̸=i j=i+1

Sous forme matricielle :

X k+1 = (D − E)−1 F X k + (D − E)−1 b

28 de 45
 On appelle la matrice G = (D−E)−1 F la matrice d’itération de Gauss-Seidel. La méthode
de Jacobi converge ssi ρ(G) < 1.

3. Méthode de relaxation :
C’est une méthode intermédiaire entre les deux méthodes précédentes. Nous mettons un
peu de diagonale de chaque côté en posant Dans cette méthode, nous prenons ; P =
1 1−ω
D − E, N = D + F , où ω ̸= 0 est un paramètre de relaxation.
ω ω
La méthode de relaxation s’exprime donc par la suite :

X (0) ∈ Rn donné


(2.16)
1
( D − E)X k+1
1−ω
D + F )X k + b k ≥ 1

=(

ω ω

D − E est une matrice triangulaire inférieure et pour calculer xk+1 en fonction de xk , on

applique l’algorithme de descente suivant : pour i = 1, n

i−1 n
xk+1 aij xk+1 aij xkj + bi )/aii
X X
i = (− j −
j=1,j̸=i j=i+1

Sous forme matricielle :

X k+1 = (D − E)−1 F X k + (D − E)−1 b

On appelle la matrice G = (D−E)−1 F la matrice d’itération de Gauss-Seidel. La méthode

de Jacobi converge ssi ρ(G) < 1.

Exercice 3
Résoudre le système : 
9x1 − 2x2 + x3 = 13







 −x1 + 5x2 + −x3 = 9



x1 − 2x2 + 93x3 = −11



à l’aide des méthodes de Jacobi et de Gauss-Seidel à partir de X 0 = (0 0 0)T (faire les

3 premières itérations seulement).

29 de 45
 CHAPITRE 3
Équations non linéaires

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2 Méthode de Newton . . . . . . . . . . . . . . . . . . . . 31
3.2.1 Interprétation géométrique . . . . . . . . . . . . . 31
3.2.2 Analyse de convergence . . . . . . . . . . . . . . . 32
3.3 Méthode de Sécante . . . . . . . . . . . . . . . . . . . . . 33
3.3.1 Analyse de convergence . . . . . . . . . . . . . . . 34
3.4 Méthode de Sécante . . . . . . . . . . . . . . . . . . . . . 35

3.1 INTRODUCTION

L’objectif des méthodes présentées ici est de résoudre pour des fonctions numériques, des
équations de la forme f (x) = 0. Puisque’il n’existe pas de formule générale pour des fonctions
aussi simples que des polynômes, il est peu probable que l’on puisse résoudre analytiquement
des équations non linéaires. Il faudra donc recourir aux méthodes numériques. Dans ce qui
suit, nous verrons plusieurs méthodes pour y arriver. Dans ce cours on va etudier certaines
méthodes itérative dont le principe est le suivant : on élabore une fonction ϕ continue telle que
f (x) = 0 =⇒ ϕ(x) = x pourx ∈ [a, b]. Puis on construit une suite (xn+1 = ϕ(x)) en partant
d’un point initial x0 ∈ [a, b]. Cette suite converge vers x∗ solution de f (x) = 0.

30 de 45
 3.2 MÉTHODE DE NEWTON

La méthode de Newton est l’une des méthodes les plus utilisées pour la résolution des
équations non linéaire pour des fonctions f assez régulières. Cette méthode possède également
une belle interprétation géométrique. L’idée de cette méthode est de remplacer la fontion non
linéaire f (x) par son approximation affinée en se basant sur l’utilisation du développement de
Taylor.
À partir d’une valeur initiale x0 de la solution, on cherche une correction δx telle que : f (x0 +
δx) = 0 Le développement de Taylor de f au voisinage de x0 d’ordre 1 est :

f (x) = f (x0 ) + f ′ (x0 )δx + o lim x −→ x0 (δx)

Donc on a
0 = f (x0 ) + f ′ (x0 )δx

et la correction
f (x0 )
δx = −
f ′ (x0 )
Puisque nous avons négligé les termes d’ordre supérieur à 1 dans le développement de Taylor,
cette correction n’est pas parfaite et l’on pose : x1 = x0 + δx.
On recommence le processus en cherchant à corriger x1 d’une nouvelle quantité δx. On obtient
alors la définition suivante :

(0)
∈ Rdonné

x



f (xn ) (3.1)
n+1
x = xn − ′ n∈N



f (xn )

3.2.1 INTERPRÉTATION GÉOMÉTRIQUE

La figure 3.2.1 permet de donner une interprétation géométrique assez simple de la méthode
de Newton. Sur cette figure, on a représenté la fonction f (x), la valeur initiale x0 et le point
(x0 , f (x0 )) qui est sur la courbe. La droite tangente à la courbe en ce point est de pente f (x0 )
et a pour équation :
y = f (x0 ) + (x − x0 )f ′ (x0 )

31 de 45
qui correspond au développement de Taylor de degré 1 autour de x0 . Cette droite coupe l’axe
des x en y = 0, c’est-à-dire en :
f (x0 )
x1 = x 0 −
f ′ (x0 )
qui devient la nouvelle valeur estimée de la solution. On reprend ensuite le même raisonnement
à partir du point (x1 , f (x1 )) et ainsi de suite.

Figure 3.1 – Interprétation géométrique de la méthode de Newton

3.2.2 ANALYSE DE CONVERGENCE

Théorème 5: Newton convergence global

Soit f : R −→ R telle que

1. f de classe C 2

2. f (x∗ ) = 0

si f ′ (x∗ ) ̸= 0 Alors
la suite : 
(0)
∈ Rdonné

x



f (xn ) (3.2)
n+1
x = xn − ′ n∈N



f (xn )
converge vers x∗ .

32 de 45
 De plus, si f ”(x) garde un sign constante (f ”(x) ̸= 0) on a la majoration de l’erreur :

|xn+1 − x∗ | < C|xn − x∗ |2

La convergence est alors au moins d’ordre 2.

Proposition 2

Si les conditions suivantes sont vérifiées :

— f ′ (x0 ) et f ”(x0 ) existent
— f (x0 ) × f ”(x0 ) > 0
cela et suffisant pour dire que la méthode de Newton converge vers une racine.

√
On veut calculer 2 en résolvant : f (x) = x2 − 2 = 0.
√ √
Dans ce cas on a f ′ (x) = 2x et f ′ ( 2) = 2 2 ̸= 0, On doit donc s’attendre à une
convergence quadratique, ce que l’on peut constater dans le tableau suivant.

Méthode de Newton

Exemple
n xn e
0 2
1 1.5 0.5
2 1,416 6666 0.0833334
3 1,414 2157 0.0024509
4 1,414 2136 0.0000021

Tout comme c’était le cas avec la méthode des points fixes, la convergence de la
méthode de Newton dépend de la valeur initiale x0 . Malgré ses belles propriétés
de convergence, une mauvaise valeur initiale peut provoquer la divergence de cette
méthode.

3.3 MÉTHODE DE SÉCANTE

La méthode de Newton possède de grands avantages, mais elle nécessite le calcul de la

dérivée de f (x) Si la fonction f (x) est complexe, cette dérivée peut être difficile à évaluer et
peut résulter en une expression complexe. On contourne cette difficulté en remplaçant le calcul

33 de 45
de la pente f (xn ) de la droite tangente à la courbe par l’expression suivante :

f (xn ) − f (xn−1 )
f ′ (x) =
xn − xn−1

Cela revient à utiliser la droite sécante passant par les points (xn , f (xn )) et (xn−1 , f (xn−1 ))
plutôt que la droite tangente passant par (xn , f (xn )). Ce choix est représenté à la figure 3.3. Il
en résulte l’algorithme suivant.

∈ Rn donné

x(0) , x(1)



f (xn )(xn − xn−1 ) (3.3)

xn+1 = xn − n∈N



f (xn ) − f (xn−1 )

Figure 3.2 – Interprétation géométrique de la méthode de Sécante

3.3.1 ANALYSE DE CONVERGENCE

Théorème 6: Convergence de méthode de sécante

Soit f : R −→ R une fonction de classe C 2 et (x∗ ) un zéro de f (x∗ ) tel que f ′ (x∗ ) ̸= 0.
Alors il existe un voisinage I de (x∗ ) tel que la suite définie par :

∈ Rdonné

x(0) , x(1)



f (xn )(xn − xn−1 ) (3.4)

xn+1 = xn − n∈N



f (xn ) − f (xn−1 )

34 de 45
 existe et converge vers (x∗ )
De plus, si f ”(x) garde un signe constante (f ”(x) ̸= 0) on a

x∗ − xn+1  1 f ”(x∗ ) p−1

lim =
x−→∞ (x∗ − x )p 2 f ′ (x∗ )
n

1 √
où p = (1 + 5)
2

3.4 MÉTHODE DE SÉCANTE

Les phénomènes non linéaires sont extrêmement courants en pratique. Ils sont sans doute
plus fréquents que les phénomènes linéaires. Dans cette section, nous examinons les systèmes
non linéaires et nous montrons comment les résoudre à l’aide d’une suite de problèmes linéaires,
auxquels on peut appliquer diverses techniques de résolution étudiées dans le chapitre 2. Le
problème consiste à trouver le ou les vecteurs x = [x1 x2 x3 ...xn ]T vérifiant les n équations non
linéaires suivantes : 
 f1 (x1 , x2 , x3 , . . . , xn ) =0







 f (x , x , x , . . . , x ) =0


2 1 2 3 n
.. (3.5)
.









fn (x1 , x2 , x3 , . . . , xn )

=0

où les fi sont des fonctions de n variables que nous supposons différentiables. Contrairement
aux systèmes linéaires, il n’y a pas de condition simple associée aux systèmes non linéaires qui
permette d’assurer l’existence et l’unicité de la solution. Le plus souvent, il existe plusieurs
solutions possibles et seul le contexte indique laquelle est la bonne. Les méthodes de résolution
des systèmes non linéaires sont nombreuses. Notamment, presque toutes les méthodes du cha-
pitre 2 peuvent être généralisées aux systèmes non linéaires. Pour éviter de surcharger notre
exposé, nous ne présentons que la méthode la plus importante et la plus utilisée en pratique,
soit la méthode de Newton. C’est à peu près la même idée que dans R. Nous remplaçons f ′ par
une matrice des dérivées partielles appelé "La Jacobienne" noté J. On aura alors :

X (0) ∈ Rn donné


(3.6)
X n+1 = X n − J −1 (X n )f (X n ) n ∈ N



35 de 45
En multipliant la formule de récurrence par J(X n ) on aura le système linéaire :

X (0) ∈ Rn donné








 J(X )(X n n
− X n+1 ) − f (X n ) = 0 n ∈ N (3.7)




J(X n )(X n − X n+1 ) = f (X n )



On pose Y = (X n − X n+1 ), la résolution du système 3.7 revient donc à résoudre un système

linéaire de la forme AY = b avec A est la Jacobienne (J),

 
∂f1 ∂f1 ∂f1
(X) (X) . . . (X)
 ∂x1 ∂x2 ∂xn


 ∂f2 ∂f2 ∂f2
 
 ∂x (X) (X) . . . (X)


J = 1 ∂x2 ∂xn 
(3.8)
.. .. .. ..

 

 . . . . 

 ∂fn ∂fn ∂fn
 
(X) (X) . . . (X)

∂x1 ∂x2 ∂xn
et b= f (X n )  
 f1 (X) 
 
f2 (X))
 
b= 
 .. 
 (3.9)

 . 

 
fn (X)

On a un résultat de convergence similaire à celui vu pour d = 1.

Exercice 4:
n cherche à trouver l’intersection de la courbe x2 = ex1 et du cercle de rayon 4 centré à
l’origine d’équation x21 + x22 = 16. L’intersection de ces courbes est une solution de :



 ex1 − x2 = 0
 x2

+ x22 − 16 = 0
1

La première étape consiste à calculer la matrice jacobienne de dimension 2. Dans ce cas,

on a :  
x1
e −1 
J =  (3.10)
2x1 2x2

Prenons le vecteur X = [2.8 2.5]T comme approximation initiale de la solution de ce

système non linéaire.

36 de 45
 1. Itération k=1 on a le système linéaire :
    
2.8
 e −1  y11  e2.8 − 2.8
= (3.11)
 
   
2(2.8) 2(2.8) y21 (2.8)2 + (2.8)2 − 16

dont la solution est Y = [0, 77890 − 0, 83604]T . La nouvelle approximation de la

solution est donc :

x11 = x01 − y11 = 2.80, 77890 = 2, 0211

x12 = x02 + y21 + = 2.8 + 0, 83604 = 3, 63604

2. Itération k=2 : On effectue une deuxième itération à partir de X 1 =

[2, 0211 3, 63604]T . Le système devient alors :
    
2,0211 1 2,0211
e −1  y1  e − 3, 63604

   =


 (3.12)
2(2, 0211) 2(3, 63604) y21 2 2
(2, 0211) + (3, 63604) − 16

dont la solution est Y = [0, 5048 − 0, 10106]T . La nouvelle approximation de la

solution est donc :

x21 = x11 − y12 = 2, 0211 − 0, 5048 = 1, 5163

x22 = x12 − y22 = 3, 63604 + 0, 101064 = 3, 7371

et on re-résout a chaque itération un système linéaire pour obtenir la solution du système

non linéaire à l’itération k.

37 de 45
 CHAPITRE 4
Interpolation

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.2 Interpolation de Lagrange . . . . . . . . . . . . . . . . . 39
4.2.1 Base de Lagrange . . . . . . . . . . . . . . . . . . 39
4.2.2 Interpolation de Lagrange . . . . . . . . . . . . . 40
4.3 Polynôme de Newton . . . . . . . . . . . . . . . . . . . . 41

4.1 INTRODUCTION

On se donne un ensemble de points (xi , fi ), les xi sont appelés abscisses ou noeuds d’inter-
polation tandis que les couples ((xi , f (xi )) pour i = 0, 1, 2, ..., n) sont les points de collocation
ou points d’interpolation obtenus suite à une mesure expérimentale (fi représente la tempé-
rature, pression, débit, ....) peut-on construire une approximation de f (x), et ce, pour tout x
pour connaître la valeur de la fonction mesurée en d’autres points dans le domaine ? Il s’agit
d’un problème d’interpolation, dont la solution est relativement simple. Il suffit de construire
un polynôme de degré suffisamment élevé dont la courbe passe par les points de collocation
P (xi ) = f (xi ).
Il y a cependant des éléments fondamentaux qu’il est important d’étudier. En premier lieu, il
convient de rappeler certains résultats cruciaux relatifs aux polynômes, que nous ne démontrons
pas.

38 de 45
 Théorème 7
Un polynôme de degré n dont la forme générale est :

Pn (x) = a0 + a1 x + a2 x2 + a3 x3 + ... + an xn (4.1)

avec an ̸= 0 possède, tenant compte des multiplicités, très exactement n racines qui
peuvent être réelles ou complexes conjuguées. Rappelons que r est une racine de Pn (x) si
Pn (r) = 0

Proposition 3

Par (n + 1) points de collocation d’abscisses distinctes ((xi , f (xi )) pour i = 0, 1, 2, ..., n),
on ne peut faire correspondre qu’un et un seul polynôme de degré n.

Definition 13
L’unique polynôme de degré n passant par les points (xi , f (xi )) pour i = 0, 1, 2, ..., n, est
appelé l’interpolant de f (x) de degré n aux abscisses (noeuds) x0 , x1 , ..., xn .

Rien n’oblige à ce que le coefficient an de l’interpolant soit différent de 0. L’interpo-

lant passant par les n + 1 points d’interpolation peut donc être de degré inférieur
à n. Si on choisit par exemple 10 points sur une droite, l’interpolant sera quand
même de degré 1.

4.2 INTERPOLATION DE LAGRANGE

4.2.1 BASE DE LAGRANGE

Soit x0 , x1 , ..., xn , n + 1 réels donnés distincts. On définit n + 1 polynômes Li pour i = 0 à
n par :
(x − x0 )(x − x1 )(x − x2 )....(x − xn )
Li (x) =
(xi − x0 )(xi − x1 )(xi − x2 )....(xi − xn )
Le numérateur de chacun de ces polynômes est un produit de n termes (x − xk ) et est donc un
polynome de degré n. Le dénominateur est une constante. On a donc
— Li est un polynôme de degré n
— Li (xk ) = 0 si i ̸= k et 0 ≤ k ≤ n
— Li (xi ) = 1

39 de 45
Réciproquement, pour i fixé, il existe un unique polynôme li vérifiant les trois propriétés pré-
cédentes.

Definition 14
Les polynômes Li (x) sont les polynômes de Lagrange de Rn [X] associés aux points
x0 , x1 , ..., xn .

4.2.2 INTERPOLATION DE LAGRANGE

Soit f une fonction donnée définie sur R à valeurs dans R et x0 , x1 , ..., xn , n + 1 réels donnés
distincts. Interpoler la fonction f par un polynôme de degré n aux points x0 , x1 , ..., xn consiste
à trouver un polynôme P de degré ≤ n tel que :

P (xi ) = f (xi ), 0≤i≤n

Si un tel polynôme existe, il s’écrit de manière unique

n
X
P (x) = f (xi )Li (x)
i=0

On doit calculer le polynôme passant par les points (0 , 1), (1 , 2), (2 , 9) et (3 , 28).
Si l’on cherche le polynôme de degré 3 passant par les 4 points on doit construire
quatre fonctions Li (x).
Exemple

(x − x1 )(x − x2 )(x − x3 ) (x − 1)(x − 2)(x − 3) (x − 1)(x − 2)(x − 3)

L0 = = =
(x0 − x1 )(x0 − x2 )(x0 − x3 ) (0 − 1)(0 − 2)(0 − 3) −6
(x − x0 )(x − x2 )(x − x3 ) (x − 0)(x − 2)(x − 3) (x − 0)(x − 2)(x − 3)
L1 = = =
(x1 − x0 )(x1 − x2 )(x1 − x3 ) (1 − 0)(1 − 2)(1 − 3) 2
(x − x0 )(x − x1 )(x − x3 ) (x − 0)(x − 2)(x − 3) (x − 0)(x − 2)(x − 3)
L2 = = =
(x2 − x0 )(x2 − x1 )(x2 − x3 ) (2 − 0)(2 − 1)(2 − 3) −2
(x − x0 )(x − x1 )(x − x2 ) (x − 0)(x − 1)(x − 2) (x − 0)(x − 1)(x − 2)
L3 = = =
(x3 − x0 )(x3 − x1 )(x3 − x2 ) (3 − 0)(3 − 1)(3 − 2) 6

40 de 45
 L’interpolation de Lagrange donne dans ce cas :

(x − 1)(x − 2)(x − 3) (x − 0)(x − 2)(x − 3)

P3 (x) = 1 +2
−6 2
(x − 0)(x − 2)(x − 3) (x − 0)(x − 1)(x − 2)
+9 + 28
−2 6
(x − 1)(x − 2)(x − 3)
P3 (x) = − + (x − 0)(x − 2)(x − 3)
6
(x − 0)(x − 2)(x − 3) (x − 0)(x − 1)(x − 2)
−9 + 14
2 3
P3 (x) = x3 + 1.

qui est l’expression du polynôme de degré 3 passant par les 4 points donnés. Enfin,
le polynôme calculé permet d’obtenir une approximation de la fonction inconnue
f (x) partout dans l’intervalle contenant les points de collocation, c’est-à-dire [0 ,
3]. Ainsi, on a : f (2.5) = P3 (2.5) = (2.5)3 + 1 = 16.625

4.3 POLYNÔME DE NEWTON

Lorsqu’on écrit l’expression générale d’un polynôme, on pense immédiatement à la forme

4.1, qui est la plus utilisée. Il en existe cependant d’autres qui sont plus appropriées au cas de
l’interpolation, par exemple :

Pn (x) = a0

+ a1 (x − x0 )

+ a2 (x − x0 )(x − x1 )

+.

+ an−1 (x − x0 )(x − x1 )(x − x2 )....(x − xn−2 )

+ an (x − x0 )(x − x1 )(x − x2 )....(x − xn−1 )

On remarque que le coefficient de an comporte n monomes de la forme (x − xi ) et qu’en consé-

quence le polynôme est de degré n.
L’aspect intéressant de cette formule apparaît lorsqu’on essaie de déterminer les (n + 1) coef-

41 de 45
ficients ai de telle sorte que Pn (x) passe par les (n + 1) points de collocation (xi , f (xi )) pour
i = 0, 1, 2, ..., n. On doit donc s’assurer que :

P (xi ) = f (xi ), 0≤i≤n

Les coefficients de la forme ?? s’annulent tous en x = x0 , sauf le premier. On peut ainsi montrer
que :
Pn (x0 ) = a0 = f (x0 )

Le premier coefficient est donc :

a0 = f (x0 ) (4.2)

On doit ensuite s’assurer que Pn (x1 ) = f (x1 ), c’est-à-dire :

Pn (x1 ) = a0 + a1 (x1 − x0 ) = f (x0 ) + a1 (x1 − x0 ) = f (x1 )

ce qui permet d’isoler a1 pour obtenir :

f (x1 ) − f (x0 )
a1 =
(x1 − x0 )

Definition 15
On définit les premières différences divisées de la fonction f (x) par :

f (xi+1 ) − f (xi )
f [xi , xi+1 ] =
(xi+1 − xi )

Ainsi, le coefficient a1 peut s’écrire :

a1 = f [x0 , x1 ] (4.3)

Il est facile de démontrer que le polynôme de degré 1 :

P1 (x) = f (x0 ) + f [x0 , x1 ](x − x0 )

obtenu en ne considérant que les deux premiers coefficients de ?? et les expressions

4.2 et 4.3, passe par les points (x0 , f (x0 )) et (x1 , f (x1 )). Il représente donc l’unique
polynôme de collocation de degré 1 passant par ces deux points.

42 de 45
 Le troisième coefficient (a2 ) est à son tour déterminé par :

Pn (x2 ) = a0 +a1 (x2 −x0 )+a2 (x2 −x0 )(x2 −x1 ) = f (x0 )+f [x0 , x1 ](x1 −x0 )+a2 (x2 −x0 )(x−x1 ) = f (x2 )

En isolant a2, on obtient :

1
(f [x1 , x2 ] − f [x0 , x1 ])
(x2 − x0 )
On en arrive donc à une expression qui fait intervenir une différence divisée de différences
divisées.

Definition 16
Les deuxièmes différences divisées de la fonction f(x) sont définies à partir des premières
différences divisées par la relation :

f [xi+1 , xi+2 ] − f [xi , xi+1 ]

f [xi , xi+1 , xi+2 ] =
(xi+2 − xi )

De même, les n-ièmes différences divisées de la fonction f (x) sont définies à partir des
(n − 1)-ièmes différences divisées de la façon suivante :

f [x1 , x2 , ..., xn ] − f [x0 , x1 , x2 , ..., xi−1 ]

f [x0 , x1 , x2 , ..., xn ] =
(xn − x0 )

Notons que les toutes premières différences divisées de f (x) (soit les 0-es différences) sont
tout simplement définies par f (xi ).

Suivant cette notation, on a :

a2 = f [x0 , x1 , x2 ] (4.4)

Il est facile de démontrer que le polynôme :

Pn (x2 ) = a0 + a1 (x1 − x0 ) + a2 (x2 − x0 )(x − x1 )

= f (x0 ) + f [x0 , x1 ](x1 − x0 ) + a2 (x2 − x0 )(x − x1 ) = f (x2 )

passe par les trois premiers points de collocation. On remarque de plus que ce
polynôme de degré 2 s’obtient simplement par l’ajout d’un terme de degré 2 au
polynôme P1 (x) déjà calculé. En raison de cette propriété, cette méthode est dite
récursive.

43 de 45
 On peut soupçonner à ce stade-ci que le coefficient a3 est :

a3 = f [x0 , x1 , x2 , x3 ] (4.5)

Théorème 8
L’unique polynôme de degré n passant par les (n + 1) points de collocation ((xi , f (xi ))
pour i = 0, 1, 2, ...., n) peut s’écrire selon la formule d’interpolation de Newton ?? ou
encore sous la forme récursive :

Pn (x) = Pn−1 (x) + an (x − x0 )(x − x1 )(x − x2 ).....(x − xn−1 )

Les coefficients de ce polynôme sont les différences divisées :

ai = f [x0 , x1 , x2 , x3 , ..., xi ] (4.6)

f [x1 , x2 , ..., xi ] − f [x0 , x1 , x2 , ..., xi−1 ]

ai = f [x0 , x1 , x2 , ..., xi ] =
(xi − x0 )

Il reste maintenant à calculer efficacement la valeur de ce polynôme. La manière la plus

simple consiste à construire une table dite de différences divisées de la façon suivante
Exemple

Calculer le polynôme passant par les points : (0 , 1), (1 , 2), (2 , 9) et (3 , 28) par
l’interpolation de Newton

44 de 45
 Suivant la formule de Newton avec x0 = 0, le polynôme de collocation est :

Pn (x) = 1 + 1(x − 0) + 3(x − 0)(x − 1) + 1(x − 0)(x − 1)(x − 2) = x3 + 1

qui est le même polynôme (en vertu de l’unicité) que celui obtenu par la méthode
de Lagrange.
Si l’on souhaite ajouter un point de collocation et calculer un polynôme de degré
4, il n’est pas nécessaire de tout recommencer. Par exemple, si l’on veut inclure le
point (5 , 54), on peut compléter la table de différences divisées déjà utilisée.

45 de 45