Ce document ne remplace pas le cours dispensé en classe

Notes de cours
Analyse numérique

E. ZAOUI

c École Nationale Supérieure des Mines de Rabat
2017
 Ce document ne remplace pas le cours dispensé en classe
Table des matières

1 Interpolation polynômiale 4
1.1 Objectif . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Existence et unicité du polynôme d’interpolation . . . . . . . . . . . . . . 4
1.3 Méthodes de construction du polynôme d’interpolation . . . . . . . . . . 5
1.3.1 Interpolation de Lagrange . . . . . . . . . . . . . . . . . . . . . . 5
1.3.2 Polynôme de Newton . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3.3 Algorithme de Neville . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Erreur d’interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Interpolation d’Hermite . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.6 Approximation au sens des moindres carrés discrète . . . . . . . . . . . . 9
1.7 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 Intégration numérique 13
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Intégration numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.1 Formules de Newton-Côtes . . . . . . . . . . . . . . . . . . . . . . 14
2.2.2 Méthode de Romberg . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.3 Quadratures de Gauss - Legendre . . . . . . . . . . . . . . . . . . 18
2.3 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3 Équations différentielles 22
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 Méthodes numériques de résolution . . . . . . . . . . . . . . . . . . . . . 22
3.2.1 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.2 Méthode de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.3 Méthodes de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . 24
3.3 Systèmes d’équations différentielles . . . . . . . . . . . . . . . . . . . . . 25
3.4 Équations d’ordre supérieur . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4 Résolution numérique des systèmes d’équations linéaires 29

4.1 Introduction et rappel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2 Conditionnement de la matrice . . . . . . . . . . . . . . . . . . . . . . . 31
4.3 Méthodes directes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.3.1 Méthode de remontée . . . . . . . . . . . . . . . . . . . . . . . . 32
4.3.2 Élimination de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . 32

2
 4.3.3 Décomposition LU . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.3.4 Méthode de Choleski . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.3.5 Système tridiagonal - Méthode de Thomas . . . . . . . . . . . . . 34
4.4 Méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.4.1 Méthode de Jacobi - algorithme . . . . . . . . . . . . . . . . . . . 35

Ce document ne remplace pas le cours dispensé en classe

4.4.2 Méthode de Gauss-Seidel - algorithme . . . . . . . . . . . . . . . 36
4.5 Système surdéterminé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.6 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

5 Système d’équations non linéaires 39

5.1 Résolution d’une équation non linéaire . . . . . . . . . . . . . . . . . . . 39
5.1.1 Approche géométrique . . . . . . . . . . . . . . . . . . . . . . . . 39
5.1.2 Méthodes de points fixes . . . . . . . . . . . . . . . . . . . . . . . 41
5.2 Résolution d’un système d’équations non linéaires . . . . . . . . . . . . . 43
5.2.1 Méthode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.2.2 Méthode de points fixes . . . . . . . . . . . . . . . . . . . . . . . 45
5.3 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3
 Ce document ne remplace pas le cours dispensé en classe
Chapitre 1

Interpolation polynômiale

1.1 Objectif
Dans ce chapitre, nous considérons le problème d’interpolation suivant :
Soit f une fonction connue seulement en (n + 1) (n ∈ N∗ ) points de la forme (xi , f (xi ))
pour i = 0, 1, 2, . . . , n.
On cherche un interpolant (une fonction) p qui vérifie :

p(xi ) = f (xi ), i = 0, . . . , n.

Les xi sont les noeuds et sont deux à deux distincts.

Les points (xi , f (xi )) sont appelés points d’interpolation et peuvent provenir de données
expérimentales ou d’une table.

Restrictions :
L’interpolant p devrait être :
1) explicite et facilement calculable,
2) facile à dériver et intégrer,
3) appartenir à une classe de fonctions suffisamment riche pour assurer l’approxi-
mation précise et efficace de toute fonction f soit possible.

Remarque 1 -
Les fonctions les plus faciles à évaluer sont les fonctions polynômes.

1.2 Existence et unicité du polynôme d’interpola-

tion
Théorème 1 Un polynôme de degré n dont la forme générale est :

pn (x) = a0 + a1 x + . . . + an−1 xn−1 + an xn (an 6= 0)

possède au plus n racines. (r est une racine de pn si pn (r) = 0)

4
Théorème 2 Par les (n + 1) points d’interpolation (xi , f (xi ))i=0,...,n , on ne peut faire
correspondre qu’un unique polynôme de degré n vérifiant pn (xi ) = f (xi ) pour i = 0, ..., n.

Exercice 1 -
Pour démontrer l’existence et l’unicité du polynôme d’interpolation, résoudre le système

Ce document ne remplace pas le cours dispensé en classe

linéaire de n + 1 équations suivant :

i=n
nX
ai xji = f (xi ), 0≤i≤n
i=0

1.3 Méthodes de construction du polynôme d’inter-

polation
1.3.1 Interpolation de Lagrange
L’interpolation de Lagrange est une façon simple de construire le polynôme d’inter-
polation :

i=n
X
pn (x) = f (xi )Li (x),
i=0

où
n
Y (x − xj )
Li (x) = .
j=0, j6=i
(xi − xj )

Exemple 1.3.1 Donner le polynôme de degré 1 qui passe par les points (0, 1) et (1, 2).

Exemple 1.3.2 Trouvez l’interpolant polynomial pour les 3 points d’interpolation : (0, 1),
(1, 2) et (2, 9).

Exemple 1.3.3 Trouvez l’interpolant polynomial pour les 4 points d’interpolation : (0, 1),
(1, 2), (2, 9), (3, 28) ?

Remarque 2 -
La méthode de Lagrange ne nous permet en aucun cas d’obtenir le polynôme pn+1 de
degré (n + 1)(en ajoutant un point d’interpolation) à partir du polynôme pn déjà calculé.
Autrement dit, cette méthode n’est pas récursive.

Exercice 2 -
Montrer que les polynômes (Li )0≤i≤n forment une famille libre et déduire que c’est une
base de l’espace vectoriel Pn des fonctions polynômes sur R à coefficients réels de degré
inférieur ou égal à n.

5
1.3.2 Polynôme de Newton
Cette fois-ci, on écrit le polynôme pn (x) dans la base de Newton :

pn (x) = b0 + b1 (x − x0 ) + b2 (x − x0 )(x − x1 ) + . . . + bn (x − x0 )(x − x1 ) . . . (x − xn−1 )

Ce document ne remplace pas le cours dispensé en classe

Xn i−1
Y
= b0 + bi (x − xj ).
i=1 j=0

Il est claire que pn (x) est un polynôme de degré n (le coefficient de bn comporte n
monômes de la forme (x − xi )).

On peut déterminer les coefficients bi par récurrence :

pn (x0 ) = b0 = f (x0 ) = p0 (x0 ).

De même, on trouve b1 :
f (x1 )−f (x0 )
pn (x1 ) = b0 + b1 (x1 − x0 ) = f (x1 ) ⇒ b1 = x1 −x0
.

Le même raisonnement donne b2 :

f (x2 )−f (x1 ) f (x1 )−f (x0 )
x2 −x1
− x1 −x0
b2 = .
x2 − x0
Exemple 1.3.4 Trouvez l’interpolant polynomial pour les 3 points d’interpolation : (0, 1),
(1, 2) et (2, 9).

Remarque 3 Pour déterminer facilement les coefficients bi , on utilise la méthode des

Différences divisées.

Définition 1 Les 0èmes différences divisées sont données par :

f [xi ] = f (xi ).

Définition 2 Les premières différences divisées de la fonction f (x) sont définies par :
f (xi+1 ) − f (xi )
f [xi , xi+1 ] = .
xi+1 − xi
Définition 3 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
D’une façon générale, les k-èmes différences divisées de la fonction f (x) sont définies à
partir des (k-1)-èmes différences divisées par la relation :
f [xi+1 , . . . , xi+k ] − f [xi , . . . , xi+k−1 ]
f [xi , xi+1 , . . . , xi+k ] = .
xi+k − xi

6
Théorème 3 -
Dans la base de Newton, l’unique polynôme de degré n qui passe par les (n + 1) points
d’interpolation (xj , f (xj )) pour j = 0, . . . , n) vérifie la relation :

pi (x) = pi−1 (x) + bi (x − x0 )(x − x1 ) . . . (x − xi−1 ), i = 1, . . . , n.

Ce document ne remplace pas le cours dispensé en classe

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

bi = f [x0 , x1 , . . . , xi ] 0 ≤ i ≤ n.

Algorithme pratique (tableau des différences divisées) :

Étape 0 Étape 1 Étape 2 ··· Étape n

x0 f [x0 ]
x1 f [x1 ] f [x0 , x1 ]
x2 f [x2 ] f [x1 , x2 ] f [x0 , x1 , x2 ]
.. .. .. .. ..
. . . . .
xn f [xn ] f [xn−1 , xn ] f [xn−2 , xn−1 , xn ] ··· f [x0 , . . . , xn ]

Remarque 4 -
Une fois les coefficients f [x0 , . . . , xi ] sont connus, on peut utiliser par exemple le schéma
d’Hörner pour évaluer le polynôme d’interpolation :

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

Ainsi, on réduit le nombre d’opérations nécessaires à l’évaluation du polynôme. En outre, cette

forme est moins sensible aux erreurs d’arrondi.

Exemple 1.3.5 Trouvez la forme de Newton du polynôme d’interpolation pour les données

xk 1 2 4 5
f (xk ) 4 16 106 208

1.3.3 Algorithme de Neville

L’algorithme découle du théorème suivant :

Théorème 4 Soient f une fonction définie sur l’ensemble {x0 , x1 , . . . , xk } et xi , xj deux points
distincts de cet ensemble. Alors, le polynôme
(x − xj )P0,1,...,j−1,j+1,...,k (x) − (x − xi )P0,1,...,i−1,i+1,...,k (x)
P (x) =
xi − xj

est le polynôme qui interpole f en {x0 , x1 , . . . , xk }.

P0,1,...,j−1,j+1,...,k est le polynôme qui interpole f en {x0 , x1 , . . . , xj−1 , xj+1 , . . . , xk }.
P0,1,...,i−1,i+1,...,k est le polynôme qui interpole f en {x0 , x1 , . . . , xi−1 , xi+1 , . . . , xk }.

7
Algorithme de Neville
Notons par Qij (x) pour 0 ≤ j ≤ i, le polynôme qui interpole f en {xi−j , . . . , xi }

Qij (x) = Pi−j,...,i (x).

Ce document ne remplace pas le cours dispensé en classe

L’algorithme de Neville est :

Qi0 (x) = f (xi ), i = 0, . . . , n

(x − xi−j )Qi,j−1 (x) − (x − xi )Qi−1,j−1 (x)

Qij (x) =
xi − xi−j
i = 0, . . . , n; j = 0, . . . , i

1.4 Erreur d’interpolation

L’interpolation permet de faire l’approximation de la fonction f en tout point x. Or, cette
opération entraı̂ne une erreur d’interpolation que l’on peut exprimer par :

En (x) = f (x) − pn (x)

Théorème 5 Soient x0 < x1 < . . . < xn des points d’interpolations d’un intervalle [a, b]
et f une fonction (n + 1) fois dérivable sur [a, b]. Alors, pour tout x de [a, b], il existe ξx ∈
] min (x, x0 ), max (x, xn )[ tel que :

f (n+1) (ξx )
En (x) = (x − x0 )(x − x1 ) . . . (x − xn ) (∗)
(n + 1)!

La relation (∗) est l’expression analytique de l’erreur d’interpolation.

Remarque 5 -

1) En (xi ) = 0, ∀i ∈ {0, . . . , n} ;
2) ξx est aussi inconnu et varie avec x ;
3) Comme l’erreur en point x fait intervenir des coefficients de la forme x − xi , il y a
intérêt de prendre les xi qui sont situés le plus près possible du point x. Il n’est pas
nécessaire d’utiliser tous les points disponibles.
4) La fonction (x − x0 ) . . . (x − xn ) est un polynôme de degré (n + 1) qui peut osciller avec
de fortes amplitudes, d’où le risque de grandes erreurs.

Théorème 6 -
Soient x0 < x1 < . . . < xn (n + 1) points et f une fonction n fois continûment dérivable La
n-ème différence divisée vérifie :

f (n) (µ)
f [x0 , . . . , xn ] =
n!
où µ ∈]x0 , xn [.

Le théorème précédent permet de construire un algorithme d’interpolation adaptative :

8
Problème :
Étant donné une tolérance tol (), on cherche une approximation pn (x) telle que

|f (x) − pn (x)| < .

Ce document ne remplace pas le cours dispensé en classe

Il suffit de calculer successivement p0 (x), p1 (x), . . . et d’estimer à chaque étape l’erreur de
l’approximation pi (x) à l’aide de

f [x0 , x1 , . . . , xi+1 ](x − x0 )(x − x1 ) · · · (x − xi ).

Exemple 1.4.1 Trouvez le polynôme d’interpolation de degré le plus bas possible de telle sorte
que l’erreur d’approximation de f (0.4) soit inférieure à 0.15

x 0 0.25 0.5 0.75 1.0

f (x) 2 5 4 4 7

1.5 Interpolation d’Hermite

Charles Hermite (1822-1901) a généralisé l’interpolation de Lagrange en faisant coı̈ncider
non seulement f et pn aux points xi , mais aussi leurs dérivées d’ordre ki aux points xi .

Théorème 7 Soient x0 , x1 , . . . , xn (n+1) points et f une fonction admettant des dérivées

jusqu’à l’ordre ki aux points xi . Posons m = n + k0 + . . . + kn . Alors, il existe un unique
polynôme pm de degré ≤ m appelé polynôme d’Hermite tel que :

p(j)
m (xi ) = f
(j)
(xi ) 0 ≤ i ≤ n et 0 ≤ j ≤ ki .

Dans le cas où les ki = 1, pm prend la forme :

i=n
X
pm (x) = ri (x)f (xi ) + Si (x)f 0 (xi )
i=0
avec :
ri (x) = [1 − 2(x − xi )L0i (x)]L2i (x)
et
Si (x) = (x − xi )L2i (x)

1.6 Approximation au sens des moindres carrés discrète

Soit f une fonction continue sur l’intervalle [a, b]. On suppose qu’on donnait exactement ou
d’une façon approchée les valeurs yi = f (xi ) que prend la fonction f en des points xi (0 ≤ i ≤ n)
de [a, b].

Définition 4 Soit p∗ ∈ Pn . On dit que p∗ est le polynôme de meilleure approximation au sens

des moindres carrés si :
n
X n
X
∗ 2
ω(xi )(yi − p (xi )) = min ω(xi )(yi − p(xi ))2
p∈Pn
i=0 i=0

où ω est fonction poids non nulle sur tous les xi .

9
 Par la suite on prend ω(x) = 1. Soit l ≤ n, pour trouver le polynôme p∗ de degré l qui
approxime au sens des moindres carrées la fonction f , on doit déterminer les coefficients qui
minimisent l’écart :
n
X
E(a) = (yi − pl (xi ))2

Ce document ne remplace pas le cours dispensé en classe

i=0
n
X n
X n
X
= yi2 − 2 yi pl (xi ) + (pl (xi ))2
i=0 i=0 i=0
n n l n l
aj xji + aj xji )2
X X X XX
= yi2 − 2 yi (
i=0 i=0 j=0 i=0 j=0
n l n n X l
yi xji ) aj xji )2 .
X X X X
= yi2 −2 aj ( + (
i=0 j=0 i=0 i=0 j=0

Pour que la quantité E(a) soit minimale, il faut que

n l n
∂E j
xj+k
X X X
=− yi x i + ak i = 0, j = 0, . . . , l.
∂aj
i=0 k=0 i=0

Théorème 8 -
Pl xi i(0 ≤ i ≤ n) sont deux à deux disjoints, alors il existe un polynôme unique
Si les points
q̄(x) = i=0 āi x de degré ≤ n qui rend minimale la quantité :

n X
X l
E(a) = ( ai xi − yi )2 , (a0 , . . . , al ) ∈ Rl+1 .
i=0 i=0

Remarque 6 Les coefficients āi (0 ≤ i ≤ l) forment l’unique solution du système linéaire

suivant :
    
S0 S1 . . . Sl ā0 v0
 S1 S2 . . . Sl+1  ā1   v1 
A= . =
    
.. .. .. .. ..
 ..
 
. . .  .   . 
Sl Sl+1 . . . S2l āl vl

où
n
X n
X
Sk = xki , vk = yi xki .
i=0 i=0

1.7 Exercices
Exercice 1
Le tableau suivant contient les valeurs d’une fonction f (x) aux 4 noeuds x0 , x1 , x2 et x3 .

x x0 = 1 x1 = 2 x3 = 4 x4 = 5
f (x) 0 4 3 2

10
 a) Trouvez la forme de Newton de l’interpolant polynômial de degré 2 pour f aux noeuds
x0 , x1 , x2 .
b) Trouvez la forme de Lagrange de l’interpolant polynômial de degré 2 pour f aux noeuds
x1 , x2 , x3 .
c) Déduire la forme de l’interpolant polynômial de degré 3 pour f aux noeuds x0 , x1 , x2 , x3 .

Ce document ne remplace pas le cours dispensé en classe

d) Soit p l’interpolant polynômial de degré 2 pour f aux noeuds x0 , x1 , x2 . Estimez l’erreur
d’approximation de f (3) par la valeur p(3).
e) Si on apprenait que la mesure originale f (x2 ) = 3 était fausse et que la vraie valeur était
f (x2 ) = 4, alors quel serait le nouveau interpolant polynômial de degré 2 aux noeuds
x0 , x1 , x2 .

Exercice 2
Une fonction f a été évaluée en 5 points et les valeurs suivantes ont été obtenues :

xi 2.0 2.5 3.0 3.5 4.0

f (xi ) 1 3 4 1 0

a) À l’aide d’un interpolant polynômiale de de degré 3, trouvez la meilleure estimation

possible de f (2.7).
b) Estimez l’erreur de l’approximation obtenue en (a).

Exercice 3
On souhaite concevoir un virage d’une voie de chemin de fer entre les points (0, 0) et (1, 1).
Le virage est décrit par une courbe de la forme y = f (x) qui satisfait :

f (0) = 0, f (1) = 1.

De plus, pour assurer une transition en douceur, la pente de la courbe doit satisfaire :
1
f 0 (0) = 0 et f 0 (1) = .
3
On représente la courbe à l’aide d’un polynôme dans l’intervalle [0, 1].

a) Quel est le degré minimale que ce polynôme devra avoir pour remplir toutes les condi-
tions ?
b) Calculer ce polynôme.

Exercice 4
Soient les valeurs du tableau suivant que l’on a obtenues en mesurant la vitesse d’un véhicule
(en km/h) toutes les 5 secondes :

temps (secondes) 0 5 10 15 20 25 30 35 40 45
vitesse (km/h) 55 60 58 54 55 60 54 57 52 49

11
 — En utilisant la fonction polynew.sci, calculer le polynôme d’interpolation de degré 9
passant par les données du tableau et tracer sur le même graphique le polynôme obtenu
et les données expérimentales.
— Utiliser le polynôme obtenu à la question (a) pour estimer la valeur de la vitesse au
temps t = 42.5s.

Ce document ne remplace pas le cours dispensé en classe

— Reprendre la question (b) avec le polynôme d’interpolation obtenu en prenant seulement
les 3 dernières données expérimentales.
— Quelles conclusions peut-on tirer des résultats obtenus en (b) et (c) ?

Exercice 5 (Convergence de l’interpolation polynomiale)

Nous voulons étudier la convergence des polynômes d’interpolation en utilisant l’exemple
de la fonction f (x) = sin(x). Le but de cette expérience numérique est de construire une suite
de polynômes d’interpolation définis à l’aide de points d’interpolation équidistants définis sur
l’intervalle [0 , 3π]. Il s’agit donc de construire un polynôme d’interpolation de degré n à l’aide
des points d’interpolation (xi , sin(xi )) pour i = 0, 1, . . . , n où x0 = 0 et xn = 3π.
— En utilisant la fonction polynew.sci, calculer les polynômes d’interpolation Pn (x) de
degré n = 4, 5, 6 et 7 et faire le graphe de la fonction sin(x) sur l’intervalle [0 , 3π] et
celui des polynômes d’interpolation obtenus.
— Évaluer l’erreur En = maxx∈[0,3π] | sin(x) − Pn (x)| pour n = 4, 5, 6 et 7 et visualiser le
graphe de En en fonction de n.
— Que peut-on conclure des résultats obtenus en (a) et (b) ?

Exercice 5 (Le phénomène de Runge)

Considérons la fonction
1
u(x) = .
1 + 25 x2
— Tracer le graphe de la fonction u(x) sur l’intervalle [−1 , 1].
— En utilisant la fonction polynew.sci, calculer les polynômes d’interpolation de degré
n = 6, 8 et 10 passant par les n + 1 points d’interpolation (xi , u(xi )) choisis de telle
sorte que les abscisses xi soient équidistantes dans l’intervalle [−1 , 1]. Tracer sur le
même graphique que celui de la question (a) le graphe des polynômes d’interpolation
obtenus.
— Reprendre la question (b) en utilisant des points d’interpolation (xi , u(xi )) aux abscisses
de Chebyshev :  
(2i + 1)π
xi = cos pour i = 0, 1, 2, . . . , n.
2(n + 1)
Ces abscisses ne sont pas équidistantes mais sont concentrées aux extrémités de l’inter-
valle [−1 , 1].
— Que peut-on conclure des résultats obtenus en (b) et en (c) ?

12
 Ce document ne remplace pas le cours dispensé en classe
Chapitre 2

Intégration numérique

2.1 Introduction
Le contenu de ce chapitre prolonge celui du chapitre précédent. Dans l’interpolation, on
cherchait à évaluer une fonction f connue seulement en quelques points. Dans le présent cha-
pitre, le problème consiste à trouver des approximations des intégrales de la forme
Z b
f (x)dx.
a

Rappelons qu’en approximant une fonction f au point x par un polynôme de degré n, on

peut procurer une certaine erreur En (x) :

f (x) = pn (x) + En (x). (∗)

13
2.2 Intégration numérique
Soient f une fonction continue sur [a, b] et a = x0 < x1 < . . . < xN = b (N = kn + 1, k ∈
N, n ∈ N∗ ) une subdivision de [a, b].
Rb
Dans les méthodes d’intégration numérique, l’intégrale a f (x)dx est remplacée par une somme

Ce document ne remplace pas le cours dispensé en classe

finie. Ces méthodes se répartissent en deux grandes catégories :

1) Les formules de Newton-Côtes ([a, b] est borné) dans lesquelles f est remplacée par
des polynômes d’interpolation ;
2) Les méthodes de Gauss ([a, b] est particulier) pour lesquelles les points xi sont
imposés.

2.2.1 Formules de Newton-Côtes

Les formules de Newton-Côtes consistent à remplacer la fonction f par des polynômes :
Z b=xN n−1
X Z xk(i+1) n−1
X Z xk(i+1)
f (x)dx = f (x)dx ≈ Pki,ki+1,...,k(i+1) (x)dx.
a=x0 i=0 xki i=0 xki


où Pki,ki+1,...,k(i+1) est le polynôme interpolant f aux points xki , . . . , xk(i+1) .

Remarque 7 -
• Si n = 1 la formule est dite simple.
• Si n ≥ 2 la formule est dite composée.

Méthode des trapèzes (k = 1)

Sur chaque intervalle [xi , xi+1 ], on remplace la fonction f par le polynôme Pi,i+1 , on obtient :
Z xi+1 Z xi+1 Z xi+1
f (x)dx = Pi,i+1 (x)dx + Ei1 (x)dx
xi x xi
Z ixi+1
= [f (xi ) + f [xi , xi+1 ](x − xi )]dx
xi
Z xi+1
f ”(ξ(x))
+ (x − xi )(x − xi+1 )dx
xi 2!
xi+1 − xi
= (f (xi ) + f (xi+1 ))
2
Z xi+1
f ”(ξi (x))
+ (x − xi )(x − xi+1 )dx.
xi 2!

Posons s = (x − x0 )/h. Nous avons :

(x − xi ) = (s − i)h et dx = hds.

Le terme d’erreur s’écrit :

Z xi+1 Z 1
f ”(ξi (x)) f ”(ξi (s))
(x − xi )(x − xi+1 )dx = s(s − 1)h3 ds.
xi 2! 0 2!

Faisant appel au second théorème de la moyenne suivant :

14
Théorème 9 Soient f1 , une fonction continue sur [a, b] et f2 , une fonction intégrable qui ne
change pas de signe sur [a, b]. Il existe alors un η ∈ [a, b] tel que :
Rb Rb
a f1 (x)f2 (x)dx = f1 (η) a f2 (x)dx. 

Puisque la fonction s(s − 1) ne change pas de signe sur [0, 1], alors :

Ce document ne remplace pas le cours dispensé en classe

Z 1 Z 1
f ”(ξi (s)) f ”(η) 3 f ”(η) 3
s(s − 1)h3 ds = h s(s − 1)ds = − h .
0 2! 2! 0 12
La méthode des trapèzes simple se résume donc à l’égalité :
Z xi+1
h f ”(ηi ) 3
f (x)dx = (f (xi ) + f (xi+1 )) − h , η ∈ [xi , xi+1 ].
xi 2 12

Méthode des trapèzes composée

Dans chaque sous-intervalle [xi , xi+1 ], on peut utiliser la méthode des trapèzes. On obtient
alors la formule des trapèzes composée :
Z b n−1
X Z xi+1 n−1
X h
f (x)dx = f (x)dx ≈ (f (xi ) + f (xi+1 ))
a xi 2
i=0 i=0

Remarque 8 On peut montrer que l’erreur totale commise dans la formule des trapèzes
composée est :
(b − a)
− f ”(η)h2 .
12

Définition 5 Les formules d’intégrations numériques sont appelées aussi formules de quadra-
ture.

Définition 6 Le degré de précision d’une formule de quadrature est N si cette formule ap-
prochée est exacte pour tout polynôme de degré ≤ N et inexacte pour au moins un polynôme
de degré N + 1.

Formule de Simpson 1/3

Reprenant le raisonnement utilisé avec la méthode des trapèzes, mais cette fois avec k = 2 dans
l’équation (??)). Ainsi, le polynôme utilisé est p2 dont la courbe passe par les points (x0 , f (x0 )),
(x1 , f (x1 )) et (x2 , f (x2 )) et qui s’écrit :

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

L’approximation qu’on cherche est :

Z x2 Z x2
f (x)dx ≈ p2 (x)dx.
x0 x0

Et puisque les noeuds xi sont également distancés on a :

Z x2
h
f (x)dx ≈ (f (x0 ) + 4f (x1 ) + f (x2 ))
x0 3

qui est la formule de Simpson 1/3 simple.

15
Remarque 9 L’analyse d’erreur est plus délicate dans ce cas. En fait,
Z x2 Z x2 (3)
h f (ξ(x))
f (x)dx − (f (x0 ) + 4f (x1 ) + f (x2 )) = (x − x0 )(x − x1 )(x − x2 )dx
x0 3 x0 3!
et on s’aperçoit que le théorème de la moyenne ne s’applique pas ici.

Ce document ne remplace pas le cours dispensé en classe

Nous admettons que l’erreur est :
Z x2
h f (4) (η) 5
f (x)dx − (f (x0 ) + 4f (x1 ) + f (x2 )) = − h , η ∈ [x0 , x2 ]
x0 3 90
Ce qui est un ordre de plus que ce quoi on s’attendait !

Méthode de Simpson 1/3 composée

Pour
R b améliorer la précision d’une quadrature simple, il faut la composer. Ainsi, pour estimer
a f (x)dx et puisque la méthode de Simpson 1/3 simple utilise 3 noeuds, il est préférable de
diviser l’intervalle [a, b] en 2n sous-intervalles et d’utiliser la méthode de Simpson 1/3 sur chaque
intervalle élémentaire.
On trouve alors :
Z b n−1
X Z x2i+2 n−1
Xh
f (x)dx = f (x)dx ≈ (f (x2i ) + 4f (x2i+1 ) + f (x2i+2 ))
a x2i 3
i=0 i=0

Dans la méthode de Simpson 1/3 composé, on utilise n fois la méthode de Simpson 1/3
simple. Ce qui fait, qu’on commet n fois l’erreur et par conséquent l’erreur totale est :
! !
f (4) (η) 5 (b − a) f (4) (η) 5 (b − a) (4)
n − h = − h =− f (η)h4
90 2h 90 180

Remarque 10 *) La méthode de Simpson 1/3 est d’ordre 4.

*) Le degré de la précision est 3.
Formule de Simpson 3/8
Dans l’équation (??), on donne à k la valeur 3, on obtient la formule de Simpson 3/8 simple
qui s’écrit :
Z x3
3h 3f (4)(η) 5
f (x)dx = (f (x0 ) + 3f (x1 ) + 3f (x2 ) + f (x3 )) − h
x0 8 80

pour un certain η ∈ [x0 , x3 ]. Pour composer cette méthode, on subdivise [a, b] en 3n sous-
intervalles de longueur b−a
3n .
L’erreur totale est :
(b − a)f (4) (η) 4
− h
80
Formule de Boole
Dans l’équation (??), on donne à k la valeur 4, on trouve la formule de Boole simple qui s’écrit :
Z x4
2h 8f (6) (η) 7
f (x)dx = [7f (x0 ) + 32f (x1 ) + 12f (x2 ) + 32f (x3 ) + 7f (x4 )] − h .
x0 45 945

Où η ∈ [x0 , x4 ]. On compose cette méthode en subdivisant l’intervalle [a, b] en 4n sous-

intervalles de longueur b−a
4n et on applique la formule de Boole simple sur chaque quadruplet

16
de sous-intervalles.
L’erreur totale est :
2(b − a)f (6) (η) 6
− h
945

Ce document ne remplace pas le cours dispensé en classe

Remarque 11 La méthode de Boole conduit à une approximation d’ordre 6 au lieu d’ordre 5.
Elle a de plus un degré de precision de 5.

2.2.2 Méthode de Romberg

La méthode de Romberg est une méthode d’intégration qui permet d’atteindre des résultats
très précis. Elle est basée sur une utilisation de la méthode des trapèzes et de la technique
d’extrapolation de Richardson.

Extrapolation de Richardson
Cette technique permet d’augmenter la précision d’une méthode d’approximation par une
technique d’extrapolation :
Notons par Qapp (h) une approximation dépendant d’un paramètre h de la quantité exacte à
approximer Qexa (inconnue). Généralement plus h est petit plus l’approximation est précise.
On suppose de plus que cette approximation est d’ordre n :

Qexa = Qapp (h) + O(hn )


(2.1)
Qexa = Qapp (h) + cn hn + cn+1 hn+1 + cn+2 hn+2 + . . .

Les cn dépendent de la méthode numérique utilisée.

La méthode de Richardson consiste à obtenir une nouvelle approximation d’ordre au moins

n + 1, à partir de l’approximation (2.1) d’ordre n. Pour ce faire, on remplace h par h/2 dans
(2.1) :

2n Qapp (h/2) − Qapp (h)

Qexa = + O(hn+1 ). (2.2)
2n − 1
Pour avoir plus de précision, la méthode de Romberg applique la technique de Richardson
sur la méthode des trapèzes. On peut montrer que le terme d’erreur de cette dernière s’écrit :

−(b − a)f ”(η)h2 /2 = c2 h2 + c4 h4 + c6 h6 + ....

Les ci sont des constantes.

Procédure :
Notons par T1,i le résultat obtenu par la méthode des trapèzes composée avec 2i−1 inter-
valle(s). Les T1,i sont donc des approximations d’ordre 2. Pour passer de T1,i à T1,i+1 , on doit
doubler le nombre de sous-intervalles, ce qui revient à diviser la valeur de h par 2.

L’extrapolation de Richardson avec n = 2 donne des approximations d’ordre 4 :

22 T1,i+1 − T1,i
T2,i = .
22 − 1

17
De même, on obtient des approximations respectivement d’ordres 6, 8, 10 et 12 :

24 T2,i+1 − T2,i
T3,i =
24 − 1

Ce document ne remplace pas le cours dispensé en classe

6
2 T3,i+1 − T3,i
T4,i =
26 − 1
8
2 T4,i+1 − T4,i
T5,i =
28 − 1
10
2 T5,i+1 − T5,i
T6,i = .
210 − 1

2.2.3 Quadratures de Gauss - Legendre

La méthode de Gauss appliquée à une fonction f conduit à une approximation de la forme :
Z b n
X
f (x)dx ' ωi f (xi ). (2.3)
a i=1

Les points xi et les poids ωi d’intégrations sont choisis de façon à ce que la quadrature(2.3)
soit exacte dans le cas des polynômes de degré le plus élevé possible.

(b−a)t+(a+b)
Avec ce changement de variable x = 2 , nous avons :

b 1
b−a
Z Z
f (x)dx = g(t)dt
a 2 −1

avec g(t) = f ( (b−a)t+(a+b)

2 ).
Il est donc possible de travailler sur [−1, 1] et considérer :
Z 1 n
X
g(t)dt ' ωi g(ti )
−1 i=1

Remarque 12 Comme il y a 2n inconnus à déterminer en ( 2.3), il est raisonnable de penser

que l’on peut atteindre le degré (2n − 1).

Quadrature de Gauss à 1 point

On cherche à déterminer ω1 et t1 tel que :
Z 1
g(t)dt = ω1 g(t1 )
−1

soit exacte dans les cas de deux polynômes 1 et t.

Nous avons donc : Z 1
1dt = ω1 = 2
−1

et Z 1
tdt = 0 = ω1 t1 = 2t1 ⇒ t1 = 0
−1

18
Ainsi, on retrouve la formule du point milieu :
Z 1
g(t)dt ' 2g(0)
−1

La quadrature de Gauss à 1 point a le même degré de précision (1) que la méthode de trapèze

Ce document ne remplace pas le cours dispensé en classe

qui utilise 2 points.

Quadrature de Gauss à 2 points

On doit maintenant déterminer les 4 coefficients inconnues de l’expression :
Z 1
g(t)dt ' ω1 g(t1 ) + ω2 g(t2 ).
−1

qui est exacte pour g(t) = 1, g(t) = t, g(t) = t2 et g(t) = t3 .

on se retrouve avec un système non linéaire mais facile à résoudre.
on trouve : r
1
ω1 = ω2 = 1; t1 = −t2 = − .
3
Remarque 13 Le degré de précision de la quadrature de Gauss à 2 points est plus élevé que
celui de la méthode de trapèze même si les deux méthodes utilisent le même nombre (2) de
points d’intégration.
Quadrature de Gauss à n points
En général, on peut montrer que les points d’intégration de Gauss-Legendre sont les racines
des polynômes (orthogonaux) de Legendre définis par :
L0 (x) = 1 et L1 (x) = x
et par la formule de récurrence :
(n + 1)Ln+1 (x) = (2n + 1)xLn (x) − nLn−1 (x).
Les ωi sont donnés par :
1 n
(t − tj )
Z Y
ωi = dt.
−1 j=1,j6=i (ti − tj )

Théorème 10 Le terme d’erreur de la quadrature de Gauss-Legendre est donné par :

22n+1 (n!)4
f (2n) (ξ) où ξ ∈ [−1, 1].
(2n + 1)((2n)!)3

n ti ωi Degré de précision
1 0 2 1
2 −0.5773502629 1 3
+0.5773502629 1
3 −0.774596669 0.555555556 5
0.0 0.888888889
0.774596669 0.555555556
4 −0.861136312 0.347854845 7
−0.339981044 0.652145155
0.339981044 0.652145155
0.861136312 0.347854845

19
2.3 Exercices
Exercice 1
Rb
Soit f une fonction convexe. L’approximation de a f (x)dx fournie par la méthode des

Ce document ne remplace pas le cours dispensé en classe

trapèzes composée est-elle plus grande ou plus petite que la valeur exacte ?

Exercice 2
Rb
Pour approximer a f (x)dx, on peut remplacer f par la constante f ( a+b
2 ).
a) Montrer que cela donne la formule d’intégration numérique simple suivante :
a+b
M1 (f ) = (b − a)f ( ).
2
b) Donner la forme composée de la formule trouvée en (a).
c) Quel est le degré de précision de cette quadrature.

Exercice 3
Soient g(t) une fonction continue et définie sur l’intervalle [−1, 1] et t1 , t2 et t3 trois points
d’intégration tels que Rt1 = −1, t2 = α et t3 = 1, où α est un nombre réel vérifiant |α| < 1. Pour
1
approcher l’intégrale −1 g(t)dt, on utilise la formule de quadrature suivante :

I(g) = ω1 g(t1 ) + ω2 g(t2 ) + ω3 g(t3 ).

a) Trouver les poids d’intégration ω1 , ω2 et ω3 en fonction de α pour que la formule de

quadrature soit de degré de précision
R 1 2.
b) Trouver ensuite α tel que I(g) = −1 g(t)dt pour tout polynôme de degré 3.
c) Réécrire cette formule sur l’intervalle [a, b] pour intégrer une fonction continue f (x) sur
cet intervalle. Identifier la formule d’intégration numérique de Newton-Côtes qui est
obtenue.

Exercice 4
On note tn le polynôme de Chebychev de degré n vérifiant :

t0 (x) = 1, t1 (x) = x, tn+2 (x) = 2xtn+1 − tn (x) n≥0

et aussi : Z 1
tn (x)tk (x) π
√ dx = δn,k
−1 1−x 2 2
où δn,k est le symbol de Kronecker.
R 1 n tm (x)
a) Calculer −1 x√1−x 2
dx pour n < m.
b) On note x0 , x1 et x3 les racines de t3 . Déterminer trois réels A0 , A1 et A2 tels que pour
tout polynôme P de degré ≤ 2 on ait
Z 1
P (x)
√ dx = A0 P (x0 ) + A1 P (x1 ) + A2 P (x2 ).
−1 1 − x2

20
 c) Montrer que l’égalité est encore vérifiée si P est de degré ≤ 5.
R1 4
d) Utiliser la question (c) pour calculer la valeur de 0 √ x dx.
x(1−x)

Exercice 5

Ce document ne remplace pas le cours dispensé en classe

1
a) Quelle sera l’erreur si l’on utilise la quadrature composée de Simpson 3 avec 4 sous-
intervalles pour estimer Z 1
x2 (1 − x) dx?
0
b) La circonférence d’un demi-cercle
√ de rayon 1 centré à l’origine est décrit par le graphe
2
de la fonction f (x) = 1 − x . Il est donc clair que la moitié de l’aire d’un cercle de
rayon 1 satisfait Z 1p
π
= 1 − x2 dx.
2 −1
Nous avons utilisé la quadrature composée du trapèze et la quadrature composée de
Simpson 13 avec h = 1/10 et 1/20 pour estimer cette intégrale. Les résultats ont été
placés dans le tableau suivant.
h trapèze Simpson 13
1/10 1.55225916 1.56350408
1/20 1.56423440 1.56822353
Estimez numériquement l’ordre de convergence de chaque suite d’approximations. Est-
ce que l’ordre de convergence observé est le même que celui prédit par la théorie ?
Donnez une explication de ce phénomène.

Exercice 6
Il est facile de vérifier que
Z 1
ln(1 + x) dx = 2 ln 2 − 1.
0

a) Utilisez la quadrature composée de Simpson 3/8 avec h = 1/4.

b) Comment doit être h pour que l’erreur soit inférieur à 10−6 ?

21
 Ce document ne remplace pas le cours dispensé en classe
Chapitre 3

Équations différentielles

3.1 Introduction
On s’intéresse à résoudre numériquement l’équation différentielle :
 0
y (t) = f (t, y(t))
(3.1)
y(t0 ) = y0

1) La variable t représente très souvent le temps.

2) La variable y dépend de t.
3) La fonction f à deux variables est supposée suffisamment différentiable.
4) La condition y(t0 ) = y0 est la Condition Initiale (CI).
5) On cherche à obtenir une approximation de la solution analytique (inconnue) de y(t).

Définition 7 L’équation différentielle ( 3.1) est dite d’ordre un, car seule la dérivée d’ordre 1
de la variable y(t) est présente.

Remarque 14 *) Avec les outils numériques de résolution, on ne peut obtenir qu’un

nombre fini d’approximation de la solution analytique aux instants ti (i = 1, . . . , N )
distancées d’une valeur hi = ti+1 − ti . On peut avoir hi = h constante pour tout i . On
appelle h le pas de temps.
*) On note par y(ti ) la solution exacte de ( 3.1) en t = ti
*) On note par yi la solution approchée en t = ti obtenue à l’aide d’une méthode numérique.

3.2 Méthodes numériques de résolution

3.2.1 Méthode d’Euler
Méthode d’Euler explicite - algorithme

1) Étant donné un pas de temps h, une CI (t0 , y0 ) et un nombre maximal d’itérations N .

2) Pour 0 ≤ n ≤ N :


yn+1 = yn + hf (tn , yn )
(3.2)
tn+1 = tn + h;

22
 3) Arrêt.

Définition 8 Une méthode de résolution d’équations différentielles est dite à un pas si elle est
de la forme :
yn+1 = yn + hφ(tn , yn ) (3.3)

Ce document ne remplace pas le cours dispensé en classe

où φ est une fonction quelconque.

Remarque 15 *) Dans une méthode à un pas on n’utilise que yn pour obtenir yn+1 .
*) Dans les méthodes à pas multiples les solutions yn−1 , yn−2 , yn−3 , ... sont également
exigées.
*) La méthode d’Euler est une méthode à un pas avec :

φ(t, y) = f (t, y).

Définition 9 On définit l’erreur de troncature locale au point t = tn par :

y(tn+1 ) − y(tn )
τn+1 (h) = − φ(tn , y(tn )). (3.4)
h
Exercice :
À l’aide du développement de Taylor autour du point t = tn , vérifier que dans le cas de la
méthode d’Euler, l’erreur de troncature locale est d’ordre 1.

3.2.2 Méthode de Taylor

Développement de Taylor en 2 dimensions
Soit U (x, y) une fonction réelle assez régulière. Le développement de Taylor autour du point
(x, y) est de la forme :
∂ ∂
U (x + h, y + k) = U (x, y) + h ∂x U (x, y) + k ∂y U (x, y)
2
h ∂ 2 2
k ∂ 2 ∂2
+ 2! ∂x2 U (x, y) + 2! ∂y2 U (x, y) + hk ∂x∂y U (x, y)
+ ...
 (n−1)
1 ∂ ∂
+ (n−1)! h ∂x + k ∂y U (x, y) + Rn .

On obtient facilement l’approximation :

h2 ∂f (tn ,y(tn )) ∂f (tn ,y(tn ))
y(tn+1 ) ≈ y(tn ) + f (tn , y(tn ))h + 2 ( ∂t + ∂y f (tn , y(tn )))

qui sera la base de la méthode de Taylor.

Méthode de Taylor d’ordre 2-algorithme

1) Étant donné un pas de temps h, une CI (t0 , y0 ) et un nombre maximal d’itérations N .
2) Pour 0 ≤ n ≤ N :

(
h2 ∂f (tn ,yn ) ∂f (tn ,yn )
yn+1 = yn + hf (tn , yn ) + 2 ( ∂t + ∂y f (tn , yn ))
(3.5)
tn+1 = tn + h;

3) Arrêt.

23
Remarque 16 On peut obtenir des méthodes de Taylor encore plus précises en poursuivant le
développement de Taylor (??) jusqu’à des termes d’ordre élevé. Toutefois, on doit calculer des
dérivées de types :
∂2f ∂2f ∂2f
, , ,....
∂t2 ∂y 2 ∂t∂y

Ce document ne remplace pas le cours dispensé en classe

Ceci rend la méthode difficile à utiliser.

3.2.3 Méthodes de Runge-Kutta

Les méthodes de Runge-Kutta ont été développées pour avoir des méthodes d’ordre plus
élevée tout en évitant les difficultés des méthodes de Taylor.

Méthodes de Runge-Kutta d’ordre 2

Ici, on cherche à remplacer la relation :
h2 ∂f (tn ,y(tn )) ∂f (tn ,y(tn ))
y(tn+1 ) = y(tn ) + hf (tn , y(tn )) + 2 ( ∂t + ∂y f (tn , y(tn ))) + O(h3 ). (3.6)

par une expression équivalente ayant le même ordre de précision O(h3 ). On propose la forme :

y(tn+1 ) = y(tn ) + a1 hf (tn , y(tn )) + a2 hf (tn + a3 h, y(tn ) + a4 h). (3.7)

Les paramètres a1 , a2 , a3 et a4 sont à déterminer de sorte que les expressions (3.6) et (3.7)
aient toutes les deux une erreur en O(h3 ). Comme :

f (tn + a3 h, y(tn ) + a4 h) = f (tn , y(tn )) + a3 h ∂f (tn∂t

,y(tn ))
+ a4 h ∂f (tn∂y
,y(tn ))
+ O(h2 ).

Alors la relation (3.7) devient :

y(tn+1 ) = y(tn ) + (a1 + a2 )hf (tn , y(tn ))

+a2 a3 h2 ∂f (tn∂t
,y(tn ))
+ a2 a4 h2 ∂f (tn∂y
,y(tn ))
(3.8)
+O(h3 ).

On constate que les expressions (3.6) et (3.8) sont du même ordre (O(h3 )). On détermine les
ai en comparant ces deux expressions terme à terme. On trouve :
1 f (tn , y(tn ))
a1 + a2 = 1, a2 a3 = et a2 a4 = .
2 2
Nous avons moins d’équations que d’inconnues, donc pas de solution unique. Autrement
dite, il y a plusieurs variantes de la méthode de Runge-Kutta. Voici les choix les plus utilisés.

Méthode d’Euler modifiée-algorithme

a1 = a2 = 21 , a3 = 1 et a4 = f (tn , y(tn )).

1) Étant donné un pas de temps h, une CI (t0 , y0 ) et un nombre maximal d’itérations N .

2) Pour 0 ≤ n ≤ N :

ŷ = yn + hf (tn , yn )
yn+1 = yn + h2 (f (tn , yn ) + f (tn+1 , ŷ)) (3.9)
tn+1 = tn + h;

24
 3) Arrêt.

Méthode du point milieu - algorithme

a1 = 0, a2 = 1, a3 = 1/2 et a4 = f (tn , y(tn )).

Ce document ne remplace pas le cours dispensé en classe

1) Étant donné un pas de temps h, une CI (t0 , y0 ) et un nombre maximal d’itérations N .
2) Pour 0 ≤ n ≤ N :

k1 = hf (tn , yn )
k1
yn+1 = yn + h(f (tn + h2 , yn + 2 ))
(3.10)
tn+1 = tn + h;

3) Arrêt.

Méthodes de Runge-Kutta d’ordre 4

1) Étant donné un pas de temps h, une CI (t0 , y0 ) et un nombre maximal d’itérations N .

2) Pour 0 ≤ n ≤ N :

k1 = hf (tn , yn )
k2 = hf (tn + h2 , yn + k21 )
k3 = hf (tn + h2 , yn + k22 )
(3.11)
k4 = hf (tn + h, yn + k3 )
yn+1 = yn + 16 (k1 + 2k2 + 2k3 + k4 )
tn+1 = tn + h;

3) Arrêt.

Remarque 17 La méthode de Runge-Kutta d’ordre 4 est très fréquemment utilisée en raison

de sa grande précision.

3.3 Systèmes d’équations différentielles

La forme générale d’un système de m équations différentielles avec conditions initiales
s’écrit :
 0
 y10 (t) = f1 (t, y1 (t), y2 (t), . . . , ym (t)) (y1 (t0 ) = y1,0 )

 y20 (t) = f2 (t, y1 (t), y2 (t), . . . , ym (t)) (y2 (t0 ) = y2,0 )



y3 (t) = f3 (t, y1 (t), y2 (t), . . . , ym (t)) (y3 (t0 ) = y3,0 ) (3.12)
 ..
 .. ..


 . = . .
 0
ym (t) = fm (t, y1 (t), y2 (t), . . . , ym (t)) (ym (t0 ) = ym,0 )

Méthode de Runge-Kutta d’ordre 4 - algorithme

1) Étant donné un pas de temps h, des conditions initiales (t0 , y1,0 , . . . ym,0 ) et un nombre
maximal d’itérations N .
2) Pour 0 ≤ n ≤ N :
Pour i = 1, 2, 3, . . . , m :

25
 ki,1 = hfi (tn , y1,n , y2,n , . . . , ym,n )
k k2,1 km,1
ki,2 = hfi (tn + h2 , y1,n + 1,1 2 , y2,n + 2 , . . . , ym,n + 2 )
k k km,2
ki,3 = hfi (tn + h2 , y1,n + 1,2 2,2
2 , y2,n + 2 , . . . , ym,n + 2 ) (3.13)

Ce document ne remplace pas le cours dispensé en classe

ki,4 = hfi (tn + h, y1,n + k1,3 , y2,n + k2,3 , . . . , ym,n + km,3 )
yi,n+1 = yi,n + 61 (ki,1 + 2ki,2 + 2ki,3 + ki,4 )
tn+1 = tn + h;
3) Arrêt.
On a alors :
f1 (t, y1 (t), y2 (t)) = y2 (t)
f2 (t, y1 (t), y2 (t)) = 2y2 (t) − y1 (t).
et la condition initiale (t0 , y1,0 , y2,0 ) = (0, 2, 1). Si on prend h = 0.1, on trouve :

y1,1 = y1,1 + 16 (0.1 + 2(0.1) + 2(0.09975) + 0.09945)

= 2.099825
y2,1 = y2,0 + 16 (+2(−0.005) + 2(−0.0055) + (−0.011075))
= 0.994651667.

3.4 Équations d’ordre supérieur

La forme générale d’une équation différentielle d’ordre m avec conditions initiales est :

y (m) (t) = f (t, y(t), y (1) , y (2) , . . . , y (m−1) ) (3.14)

(1) (2) (m−1)
y(t0 ) = c1 , y (t0 ) = c2 , y (t0 ) = c3 , . . . , y (t0 ) = cm (3.15)

y (i) (t) désigne la ième dérivée de y(t).

Théorème 11 L’équation différentielle d’ordre m ( 3.14) est équivalente au système de m

équations d’ordre 1 suivant :
 0
 y1 (t) = y2 (t) y1 (t0 ) = c1
0 (t)

 y
 20

 = y 3 (t) y (t
2 0 ) = c2
 y3 (t)

= y4 (t) y3 (t0 ) = c3
.. ..


 . .

 y 0 (t) = y (t) y (t ) = cm−1
m m−1 0
 m−1


0
ym (t) = f (t, y1 (t), y2 (t), . . . , ym (t)) ym (t0 ) = cm

Remarque 18 Pour voir l’équivalence, poser yi (t) = y (i−1) (t), i = 1, . . . , m.

3.5 Exercices
Exercice 1
Considérons l’équation différentielle

ty 0 (t) + 2y(t) = 4t2 ,

y(1) = 3.

26
 a) À l’aide de la méthode d’Euler explicite, faites deux pas de temps de grandeur ∆t = 0.5
pour estimer y(2).
b) À l’aide de la méthode du point milieu, faites un pas de temps de grandeur ∆t = 1.0
pour estimer y(2).
c) Ci-dessous, vous trouverez les erreurs des approximations de y(2) pour différentes va-

Ce document ne remplace pas le cours dispensé en classe

leurs du pas de temps ∆t et obtenus à l’aide des méthodes d’Euler explicite et du point
milieu.
∆t erreur avec Euler erreur avec pt milieu
1.000 3.500000 1.833333
0.500 0.833333 0.309524
0.250 0.369048 0.056297
0.125 0.175000 0.011965
Pour la même quantité de travail, quelle méthode est plus précise ? Justifiez.

Exercice 2
La méthode de Ralston pour la résolution numérique d’une équation différentielle y 0 (t) =
f (t, y) est

k1 = f (ti , yi ),
3 3

k2 = f ti + h, yi + hk1 ,
4 4
h
yi+1 = yi + (k1 + 2k2 ).
3
À l’aide du développement de Taylor montrez que l’erreur de troncature est O(h2 ). Indice : En
dérivant l’équation différentielle, on obtient une identité très utile y 00 (t) = ∂f /∂t+∂f /∂y ·y 0 (t).

Exercice 3
Écrivez les équations suivantes comme un système d’équations d’ordre 1.
(a) L’équation de Van der Pol :
y 00 = y 0 (1 − y 2 ) − y.
(b) L’équation de Blasius ;
y 000 = −yy 00 .

Exercice 4
Les équations non linéaires de Lotka-Volterra modélisent la dynamique de systèmes biolo-
giques dans lesquels un prédateur et sa proie interagissent. Les équations sont
dx
(t) = ax(t) − bx(t)y(t),
dt
dy
(t) = cx(t)y(t) − dy(t),
dt
où x(t) est l’effectif des proies, y(t) est l’effectif des prédateurs, et a, b, c, d sont des paramètres
caractérisant les intéractions entre les deux espèces.

27
 Faites une itération de la méthode du point milieu avec h = 0.5, a = 1, b = c = 0.025,
d = 0.1 et avec comme valeur initiales (x(0), y(0)) = (50, 10).

Ce document ne remplace pas le cours dispensé en classe

28
 Ce document ne remplace pas le cours dispensé en classe
Chapitre 4

Résolution numérique des systèmes

d’équations linéaires

4.1 Introduction et rappel

La résolution des systèmes d’équations linéaires est sans doute le problème numérique que
l’on rencontre le plus fréquemment. Le système général, qui comporte n équations et n inconnues
s’écrit :

a11 x1 +a12 x2 . . . +a1n xn = b1

a21 x1 +a22 x2 . . . +a2n xn = b2
.. .. .. .. ..
. . . . .
an1 x1 +an2 x2 . . . +ann xn = bn

Avec la notation matricielle, le système précédent s’écrit :

Ax = b, (4.1)

où A = (aij )i,j=1,...,n et b = [b1 , . . . , bn ]T et x = [x1 , . . . , xn ]T .

La matrice A et le vecteur b sont connus.

Remarque 19 La notation zT désigne le transposé de l’objet z, vecteur ou matrice.

Définition 10 Une norme vectorielle est une application de Rn dans R qui associe à un vecteur
x un scalaire noté ||x|| et qui vérifie les trois propriétés suivantes :
1) ||x|| > 0 sauf si x = ~0
2) Si α est un scalaire, alors :
||αx|| = |α|||x||
3) L’inégalité triangulaire entre deux vecteurs x et y quelconques :

||x + y|| ≤ ||x|| + ||y||.

Exemple 4.1.1 Les normes les plus connues sont :

a) La norme euclidienne : q
||x||2 = x21 + x22 + . . . + x2n

29
 b) La norme l1 :
n
X
||x||1 = |xi |
i=1
c) La norme l∞ :

Ce document ne remplace pas le cours dispensé en classe

||x||∞ = max |xi |.
1≤i≤n

Puisque nous nous intéressons aux systèmes linéaires, il importe de pouvoir définir des
normes relatives aux matrices.

Définition 11 Une norme matricielle est une application qui associe à une matrice A un
scalaire noté ||A|| et qui vérifie les quatre propriétés suivantes :
1) ||A|| > 0 sauf si A = 0
2) Si α est un scalaire, alors :
||αA|| = |α|||A||
3) L’inégalité triangulaire entre deux matrices quelconques :

||A + B|| ≤ ||A|| + ||B||

4)
||AB|| ≤ ||A||||B||.

Remarque 20 On peut construire une norme matricielle à partir d’une norme vectorielle
quelconque en posant :
||A|| = sup ||Ax||.
||x||=1

Une norme matricielle ainsi construite est dite norme induite par la norme vectorielle. On peut
aussi dire que la norme matricielle ainsi construite est subordonnée à la norme vectorielle.

Définition 12 (norme de Frobenius) On définit la norme de Frobenius par :

v
u n 2
uX
||A||F = t aij .
i,j=1

Cette norme n’est induite par aucune norme vectorielle.

Définition 13 Une norme vectorielle et une norme matricielle sont dites compatibles si la
condition :
||Ax|| ≤ ||A||||x||
est valide quels que soient la matrice A et le vecteur x.

Définition 14 Une matrice est dite triangulaire inférieure (ou supérieure) si tous les aij (ou
tous les aji ) sont nuls pour i < j.

Définition 15 Une matrice A est dite symétrique si :

AT = A.

Définition 16 Une matrice symétrique est dite définie positive si :

Ax · x > 0 ∀x 6= ~0.

30
Définition 17 Soient v un vecteur non nul de Rn et A une matrice carrée d’ordre n. v est
un vecteur propre de A associé à la valeur propre λ ∈ C, si la relation : Av = λv est vérifiée.

Définition 18 Le rayon spectral d’une matrice A est défini par :

ρ(A) = max |λi |

Ce document ne remplace pas le cours dispensé en classe

1≤i≤n

où |λi | est le module de la valeur propre λi de A.

Remarque 21 La norme matricielle induite par la norme euclidienne vérifie :

q
||A||2 = sup ||Ax||2 = (ρ(AT A)).
||x||2 =1

Définition 19 Une matrice A est dite à diagonale strictement dominante si :

Xn
|aii | > |aij |, ∀i
j=1,j6=i

Par la suite, on suppose que la matrice A est inversible.

Définition 20 Une méthode de résolution d’un système linéaire est dite directe si elle permet
d’obtenir la solution du système en un nombre fini et prédéterminé d’opérations.

4.2 Conditionnement de la matrice

Définition 21 Le conditionnement d’une matrice A (noté condA) est défini par :
condA = ||A||||A−1 ||
où A−1 est la matrice inverse de A.

Remarque 22 Le conditionnement dépend de la norme matricielle utilisée. On utilise le plus

souvent la norme ||A||∞ .
Considérons un système  perturbé  : le second membre a subi une petite modification
∆b, ce qui a entrainé une petite variation ∆x de la solution
A(x + ∆x) = b + ∆b
d’où ∆x = A−1 ∆b et :
||∆x|| ≤ ||A−1 || ||∆b||.
De l’inégalité précédente, nous tirons la majoration suivante :
||∆x|| ||∆b|| ||∆b||
≤ ||A−1 || ||A|| = condA . (4.2)
||x|| ||b|| ||b||
Dans (4.2), condA représente la sensibilité de l’erreur relative aux variations du second membre.

En considérant un nouveau système perturbé :

(A + ∆A)(x + ∆x) = b.
On trouve :
||∆x|| ||∆A||
≤ condA . (4.3)
||x|| ||A||
Ce qui montre que condA représente aussi la sensibilité de l’erreur relative aux variations
des coefficients des premiers membres.

31
4.3 Méthodes directes
4.3.1 Méthode de remontée
Lorsque A est une matrice triangulaire supérieure (ou inférieure), la résolution est immédiate :

Ce document ne remplace pas le cours dispensé en classe



 a11 x1 + a12 x2 + . . . + a1n xn = b1
...


 an−1,n−1 xn−1 + an−1,n xn = bn−1
an,n xn = bn


On calcule successivement xn à partir de la dernière équation, puis xn−1 à partir de l’avant

dernière équation et ainsi de suite :


 xn = bn /ann
xn−1 = (bn−1 − (an−1,n xn ))/an−1,n−1


 . ..
x1 = (b1 − (a12 x2 + . . . + a1n xn ))/a11


Remarque 23 La méthode de remontée nécessite n(n − 1)/2 additions, n(n − 1)/2 multipli-
cations et n divisions.

4.3.2 Élimination de Gauss

Le procédé de l’élimination de Gauss consiste à transformer le système (4.1) en un nouveau
système triangulaire.

Exemple 4.3.1 Utiliser l’élimination de Gauss pour résoudre le système linéaire suivant :

E11 : 3x1 + 2x2 + x3 = 1

E21 : 1x1 + 3x2 + 2x3 = 2
E31 : 2x1 + 4x2 + 6x3 = 3.

Solution 4.3.1 Solution donnée en classe.

Plus généralement, lorsque les divisions par zéro n’apparaissent pas, on remplace à chaque étape
(k)
la matrice A = A(1) par une matrice A(k) = (aij ) dont les k-ièmes premiers vecteurs colonnes
correspondent au début d’une matrice triangulaire. À la (k + 1)-ième étape, on conserve les k
premières lignes et les (k − 1) premières colonnes de A(k) :
(k)
(k+1) (k) aik (k)
aij = aij − (k) akj , i = k + 1, . . . , n et j = k + 1, . . . , n
akk
(k+1)
aik = 0, i = k + 1, . . . , n
(k)
(k+1) (k) aik (k)
bi = bi − (k) bk
akk

Après la n-ème étape on obtient le système triangulaire :

 (1) (1) (1) (1)
   (1) 
a11 a12 a13 · · · a1n x1 b1
(2) (2) (2)    (2)
 0 a22 a23 · · · a2n   x2   b2
  

  ..  =  .
 . ..  .
..

 .
. .   .   ..

 . 
(n) xn (n)
0 0 0 · · · ann bn

32
équivalent au système initial (4.1).

(k)
Remarque 24 Les akk sont les pivots.
On note généralement la matrice triangulaire supérieure A(n) par U .

Ce document ne remplace pas le cours dispensé en classe

4.3.3 Décomposition LU
Dans la méthode de Gauss, l’opération d’addition à une ligne de A d’une autre ligne mul-
tipliée par un scalaire peut se décrire comme un produit matriciel.
Il est facile de vérifier que :

A(k+1) = Mk A(k) .

Où Mk est la matrice donnée par :

1 ··· 0 ··· 0
 
0
 .. . . .. .. .. 
 . . . . . 
 
 0 1 0 0 
Mk =  
 0
 −mk+1,k 1 ··· 0  
 . . .. .. . . .. 
 .. .. . . . . 
0 · · · −mn,k 0 ··· 1

(k) (k)
Avec mik = aik /akk ; i = k + 1, · · · , n.
Et que :
Mn−1 Mn−2 · · · M1 A = A(n) = U.

Si on définit la matrice L par :

−1
L = (Mn−1 Mn−2 · · · M1 )−1 = M1−1 M2−1 · · · Mn−1

on obtient la factorisation suivante :

A = LU.

Une fois les matrices L, U obtenues, la résolution du système Ax = b se fait en deux étapes :

Ly = b
U x = y.

Théorème 12 Soit A ∈ Rn×n , La factorisation LU de A avec Lii = 1 pour i = 1, . . . n existe

et unique si et seulement si les n sous matrices Ak = (aij )ij=1,...,k , k = 1, . . . , n de A sont
inversibles.

Théorème 13 Soient 2 ≤ n et A ∈ Rn×n , il existe une matrice de permutation P ∈ Rn×n , une

matrice triangulaire inférieure L ∈ Rn×n avec 1 sur la diagonale et une matrice triangulaire
supérieure U ∈ Rn×n tel que :
P A = LU.

33
Algorithme de la décomposition LU (sans permutation) :
Pour k = 1, . . . , n
k−1
X
ukj = akj − lkr urj , j = k, . . . , n

Ce document ne remplace pas le cours dispensé en classe

r=1
k−1
1 X
lik = (aik − lir urj ), i = k + 1, . . . , n
ukk
r=1

avec toujours Lii = 1.

Exemple 4.3.2 Utiliser la méthode de décomposition LU pour résoudre le système linéaire

suivant :

3x1 + 2x2 + x3 = 1
x1 + 3x2 + 2x3 = 2
2x1 + 4x2 + 6x3 = 3.

4.3.4 Méthode de Choleski

Définition 22 On dit que la matrice A est définie positive si :

∀x ∈ Rn x 6= 0 ⇒ (Ax, x) > 0.

Théorème 14 Une matrice A est symétrique définie positive si et seulement si il existe une
matrice L triangulaire inférieure inversible telle que A = LLT .

Algorithme de Choleski

Pour k = 1, . . . , n
v
u
u k−1
X
lkk = akk −
t l2 kj
j=1
k−1
X
lik = (aik − lij lkj )/lkk i = k + 1, . . . , n .
j=1

4.3.5 Système tridiagonal - Méthode de Thomas

De nombreuses applications mènent à un système d’équations linéaires tridiagonal à diago-
nale dominante. Dans ces conditions, la matrice des coefficients admet une décomposition LU
sans permutation de lignes :

    
a1 c1 0 ... ... 0 1 0 0 ... ... 0 v1 u1 0 ... ... 0

 b2 a2 c2 0 ... 0  
  l2 1 0 0 ... 0 
 0 v2 u2 0 ... 0 


 0 ... ... ... ... 0 =
  0 ... ... ... ... 0 
 0 ... ... ... ... 0 

 0 . . . . . . bn−1 an−1 cn−1   0 . . . . . . ln−1 1 0  0 ... . . . 0 vn−1 un−1 
0 ... ... 0 bn an 0 ... ... 0 ln 1 0 ... ... 0 0 vn

34
 avec
v1 = a1 , u1 = c1 ,
ui = ci , 1 ≤ i ≤ n − 1;
li = bi /vi−1 ,

Ce document ne remplace pas le cours dispensé en classe

vi = ai − li ui−1 , 2 ≤ i ≤ n.

4.4 Méthodes itératives

Dans les méthodes itératives, le système Ax = b est mis sous la forme :
M~x = N~x + b. (A = M − N )
Lorsque la matrice M est inversible on obtient le problème de point fixe suivant :
~x = M −1 N~x + M −1 b.
Par conséquent, les méthodes itératives sont des méthodes de points fixes. La détermination
de la solution repose sur l’itération de l’équation :
~xk+1 = M −1 N~xk + M −1 b. (4.4)
L’algorithme est initialisé par x~0 et s’arrête quand la quantité ||~xk+1 − ~xk || devient assez petit
par rapport à une tolérance  choisie par l’utilisateur.

Théorème 15 Les conditions suivantes sont équivalentes :

1. La matrice A est convergente.

2. Pour toute norme matricielle :
lim ||An || = 0.
n→∞
3. Pour tout vecteur ~x :
lim An ~x = 0
n→∞
.
4. ρ(A) < 1.

Théorème 16 Soit x~k une suite construite à partir de ( 4.4). La suite x~k converge vers ~x
indépendamment du choix de x~0 si et seulement si ρ(M −1 N ) < 1.
Selon les choix des matrices M et N , on a différentes méthodes itératives. Par la suite, on
note D la matrice formée des éléments diagonaux de A, E la matrice formée par les −aij , (i > j)
et F la matrice formée par les −aij , (i < j), de sorte que A = D − (E + F ).

4.4.1 Méthode de Jacobi - algorithme

Cette méthode consiste à prendre :
M =D et N = E + F.
On construit l’algorithme :
* Pour x(0) donné
* Pour i = 1, . . . , n
n
(k+1) 1 X (k)
xi = (bi − aij xj ).
aii
j=1,j6=i

35
4.4.2 Méthode de Gauss-Seidel - algorithme
En prenant :
M =D−E et N = F,
on obtient l’algorithme de la méthode de Gauss-Seidel suivant :

Ce document ne remplace pas le cours dispensé en classe

* Pour x(0) donné
* Pour i = 1, . . . , n
i−1 n
(k+1) 1 X (k+1)
X (k)
xi = (bi − aij xj − aij xj ).
aii
j=1 j=i+1

4.5 Système surdéterminé

Dans cette section, nous considérons un système d’équations linéaires Ax = b avec A ∈
Rm,n , b ∈ Rm où (m > n) et rang(A) = n. Dans ce cas, on dit que le système est surdéterminé.

La méthode des moindres carrés consiste à rechercher, parmi les x ∈ Rn , celui (ou ceux)
qui minimisent la quantité :
r(x) = ||Ax − b||22 (4.5)
appelée fonction résidu. Tout vecteur x∗ vérifiant :

r(x∗ ) = infn r(x)

x∈R

est appelé solution au sens de moindres carrés du système Ax = b.

Théorème 17 Soient A ∈ Rm,n , rangA = n et b ∈ Rm . Le problème de moindres carrés

infx∈Rn ||Ax − b||22

a une unique solution :

x∗ = (AT A)−1 AT b
.

Définition 23 L’équation
AT Ax = AT b
est appelée l’équation normale.

Exemple 4.5.1 Régression linéaire

On suppose que des mesures physiques ont été effectuées :

(xi , yi ), 1 ≤ i ≤ m, m > 2,

où xi ∈ R est le paramètre de la mesure et yi ∈ R le résultat obtenu. Le modèle linéaire consiste

à supposer yi = αxi + β. En général, il est impossible de trouver α et β pour les quels il y a
égalité ∀i. On cherche donc une solution au sens des moindres carrés :
n
X
inf (αxi + β − yi )2 .
α, β
i=1

36
4.6 Exercices
Exercice 1
En utilisant l’élimination de Gauss, résoudre le système A~x = ~b , où :

Ce document ne remplace pas le cours dispensé en classe

   
1 2 1 0
A=  2 2. 3  ~
et b =  3 
−1 −3 0 2

Exercice 2
On veut résoudre le système linéaire A~x = ~b où :
   
2. 0. 0. 1. 3
 0. 2. 1. 0.  3 
~ 

A=  0. 2. 4. 1.  et b = 
 
7 
2. 0. 3. 6. 11
par des méthodes directe et itérative.
1) Utiliser l’élimination de Gauss pour triangulariser la matrice A.
2) À partir de la question précédente, déduire les matrices L et U telles que A = LU .
3) Résoudre le système d’équations linéaires en question.
4) À partir de ~x(0) = [1.2, 0.8, 1.2, 0.8]T , effectuer deux itérations de la méthode de Gauss-
Seidel pour estimer la solution.
5) La méthode de Gauss-Seidel converge t-elle pour ce système ? Justifier.

Exercice 3
Plusieurs applications mènent à des matrices symétriques. Dans ce cas, une factorisation
de la forme LLt , où L est une matrice triangulaire inférieure, est plus simple à exécuter qu’une
factorisation LU .
a) Expliquer pourquoi on s’intéresse à ce cas particulier de la factorisation LU .
b) Résoudre le système linéaire

9x + 3y + 3z = 15;

3x + 5y + z = 9;

3x + y + 5z = 9,


à l’aide d’une factorisation LLt .

Exercice 4
En utilisant l’élimination de Gauss, résoudre le système A~x = ~b , où :
   
1 1 2 1 2
 2 2. 5 3   et ~b =  4 
 
A=  1 3 3 3   −2 
1 1 4 5 2

37
Exercice 5
Soient    
4 1 0 1 4
 1 5 1 0   7 
A= 
 0
b =  
1 6 1 16 

Ce document ne remplace pas le cours dispensé en classe

 
1 0 1 4 14
et la solution du système Ax = b, x = [0, 1, 2, 3]T .

a) À partir de x(0) = [0, 0, 0, 0]T , faites trois itérations de la méthode de Jacobi pour estimer
la solution.
b) Répéter la question précédente avec la méthode de Gauss-Seidel.

Exercice 6
Considérons les deux systèmes d’équations linéaires suivants :
   
2 1 x1 3
S1 :
2 1.001 x2 0
   
2 1 x1 3
S2 :
2 1.002 x2 0

a) Vérifier que [1501.5, −3000]T est solution de S1.

b) Vérifier que [751.5, −1500]T est solution de S2.
c) Commenter et expliquer le résultat.

38
 Ce document ne remplace pas le cours dispensé en classe
Chapitre 5

Système d’équations non linéaires

Dans ce chapitre, on se propose d’approcher numériquement les racines d’une fonction réelle
et aussi présenter quelques méthodes d’approximations des solutions des systèmes non linéaires.

5.1 Résolution d’une équation non linéaire

On se propose d’approcher les racines de l’équation ;

f (x) = 0 (5.1)

Où f : [a, b] ⊂ R → R est une fonction suffisamment régulière sur l’intervalle [a, b].

Définition 24 Soit α ∈ [a, b]. On dit que α est une racine de f , si f (α) = 0.

Définition 25 On dit qu’une suite xn est convergente vers α avec l’ordre p si :

|xn+1 − α|
∃C > 0 : ≤ C; ∀n ≥ n0
|xn − α|p

où n0 ∈ N. Dans ce cas la méthode est dite d’ordre p.

Remarque 25 .
Si p = 1, pour que la suite xk converge vers α il faut que C < 1.

5.1.1 Approche géométrique

Méthode de la dichotomie

Cette méthode repose sur la propriété suivante :

Proposition 1 Soit f : [a, b] ⊂ R → R une fonction continue telle que f (a)f (b) < 0, alors
∃α ∈ [a, b] tel que f (α) = 0.
Posons a0 = a, b0 = b et x0 = a0 +b 2 .
0

Si f (a)f (b) < 0, alors la racine α se trouve dans l’un des deux intervalles [a0 , x0 ] ou [x0 , b0 ].
Pour savoir lequel, on utilise de nouveau la propriété (1) : Si f (a0 )f (x0 ) < 0, on pose a1 = a0 ,
b1 = x0 et x1 = a1 +b2 .
1

39
Sinon, on pose a1 = x0 , b1 = b0 et x1 = a1 +b
2 .
1

ak +bk
En itérant ce procédé, on obtient une suite de valeurs (xk = 2 )k qui vérifie :

b−a
|xn − α| < . (5.2)
2n+1

Ce document ne remplace pas le cours dispensé en classe

Remarque 26 .

1 La formule ( 5.2) permet d’estimer à l’avance le nombre d’itérations nécessaires pour

approcher α avec une précision donnée.
2 La convergence de la méthode de dichotomie n’est pas rapide, mais elle est sûre à partir
du moment où on a un intervalle avec changement de signe.
Méthode de la sécante

Soient x0 et x1 deux point de [a, b] . La droite qui passe par M0 = (x0 , f (x0 )) et M1 =
(x1 , f (x1 )) coupe l’axe des x en un point d’abscisse x2 . On recommence avec les points M1 et
M2 = (x2 , f (x2 )) pour obtenir l’abscisse x3 .
En itérant ce procédé, on obtient ainsi une suite xn définie par :

f (xn )
xn+1 = xn − (xn − xn−1 )
f (xn ) − f (xn−1 )

Remarque 27 .

* On choisit les valeurs initiales le plus près possible de la racine. Il n’est pas nécessaire
qu’il y ait un changement de signe dans l’intervalle [x0 , x1 ].
* La convergence de la méthode est linéaire.
Méthode de Régula-Falsu

Cette méthode vient pour améliorer la méthode de la sécante en s’inspirant de la méthode

de dichotomie. Plus précisément, posons :

a0 = a, et b0 = b.

On définit x1 comme l’intersection de l’axe des x et la courbe joignant Ma = (a, f (a)) et

Mb = (b, f (b)) :
a0 f (b0 ) − b0 f (a0 )
x1 = .
f (b0 ) − f (a0 )
Si f (x1 ) est de même signe que f (b0 ), on pose :

a1 = a0 , b1 = x1 .

Si f (x1 ) est de même signe que f (a0 ), on pose :

a1 = x1 , b1 = b0 .

Et en itérant le processus, on définit ainsi une suite d’abscisse :

an f (bn ) − bn f (an )
xn+1 =
f (bn ) − f (an )

40
Remarque 28 Généralement, Cette méthode converge rapidement.
Méthode de Newton

Soit x0 un point de [a, b]. La tangente de la courbe y = f (x) au point M0 = (x0 , f (x0 ))

Ce document ne remplace pas le cours dispensé en classe

coupe l’axe des x en un point d’abscisse x1 . De même, la tangente de la courbe y = f (x) coupe
l’axe des x en un point d’abscisse x2 .

En itérant ce procédé, on obtient une suite d’abscisse xn définie par :

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

5.1.2 Méthodes de points fixes

Définition 26 Un point fixe d’une fonction F (x) est une valeur qui vérifie :
F (x) = x

Proposition 2 Soit F : I ⊂ R → I une fonction lipschitzienne de rapport K, avec 0 < K < 1.

∀x1 , x2 ∈ I, |F (x1 ) − F (x2 )| ≤ K|x1 − x2 |. (5.3)
Alors, la suite définie par xn+1 = F (xn ) converge vers le point fixe α = F (α) et ceci ∀x0 .

Remarque 29 .

1 L’équation ( 5.1) peut toujours s’écrire sous la forme :

F (x) = x.
On posera par exemple F (x) = x + f (x). Cette transformation n’est pas unique.
2 La condition ( 5.3) peut être remplacée par
|F 0 (x)| < 1, ∀x ∈ I.
3 On peut facilement déterminer des points fixes à partir de l’algorithme :

x0 donné
xn+1 = F (xn )

Exemple 5.1.1 Résoudre l’équation x2 − 2x − 3 = (x + 1)(x − 3) = 0. On a au moins les

transformations suivantes :
√
F1 (x) = 2x + 3 = x;
3
F2 (x) = = x;
x−2
x2 − 3
F3 (x) = =x
2
Appliquons l’algorithme des points fixes à chacune des fonctions Fi (x), i = 1, . . . , 3 en
partant de x0 = 4, on obtient :
x10 = F1 (x9 ) = 3.0000157
x10 = F2 (x9 ) = −1.000338
x4 = F3 (x3 ) = 18252.43.

41
Remarque 30 À partir de l’exemple précédent, on peut conclure qu’un algorithme des points
fixes, selon le choix de la fonction F , converge vers l’une ou l’autre racine et peut même diver-
ger complètement :
L’algorithme xn+1 = F1 (xn ) semble converger vers la racine 3.
L’algorithme xn+1 = F2 (xn ) semble converger vers la racine −1.

Ce document ne remplace pas le cours dispensé en classe

L’algorithme xn+1 = F3 (xn ) semble diverger.

Convergence de la méthode des points fixes

Théorème 18 Soit xk+1 = F (xk ), pour k = 0, 1, . . . et x0 donné. Si

1- F : [a, b] → [a, b];
2- F ∈ C 1 ([a, b]) ;
3- ∃K < 1 : |F 0 (x)| ≤ K, ∀x ∈ [a, b].
Alors, F a un point fixe α ∈ [a, b] et la suite {xn }n converge vers α dans [a, b] pour tout
x0 ∈ [a, b]. De plus, nous avons :
xn+1 − α
lim = F 0 (α).
n→∞ xn − α

Théorème 19 Soit α un point fixe d’une fonction F de classe C 1 sur un voisinage J de α. Si

|F 0 (α)| < 1 alors il existe δ > 0 tel que la suite (xn )n converge vers α pour tout x0 vérifiant :
|x0 − α| < δ.

Soit α tel que f (α) = 0 et F (α) = α.

Notons en = xn − α, l’erreur à l’étape n.
xn tend vers α si et seulement si en tend vers 0.

Nous avons :
en+1 = xn+1 − α = F (xn ) − F (α).
On constate aussi :
xn = α + (xn − α) = α + en .
Le développement de Taylor de F autour de α donne :

en+1 = F (α + en ) − F (α)
e2n 00
en+1 = (F (α) + en F 0 (α) + F (α) + . . .) − F (α)
2
Si F 0 (α) 6= 0 et si on néglige les termes d’ordre supérieur ou égal à 2, on obtient :

en+1 = F 0 (α)en . (5.4)

Remarque 31 .

* L’erreur à l’étape (n + 1)est proportionnelle à l’erreur à l’étape n.

* L’erreur ne peut diminuer que si :

|F 0 (α)| < 1. (5.5)

42
 * Le signe de F 0 (α) a une influence sur la convergence.
* La convergence de la méthode des points fixes dépend aussi du choix de la valeur ini-
tiale x0 . Un mauvais choix de x0 peut résulter en un algorithme divergent même si la
condition ( 5.5) est respectée.

Ce document ne remplace pas le cours dispensé en classe

Définition 27 Le bassin d’attraction de la racine α pour la méthode de points fixes xn+1 =
F (xn ) est l’ensemble de valeurs initiales x0 pour lesquelles xn converge vers α quand n tend
vers l’infini.

Définition 28 Un point fixe α de la fonction F est dite attractif si :

|F 0 (α)| < 1

et répulsif si :
|F 0 (α)| > 1.
Le cas |F 0 (α)| = 1 est indéterminé.

5.2 Résolution d’un système d’équations non linéaires

Dans cette section, nous nous intéressons aux systèmes non linéaires. Le problème consiste à
trouver le vecteur x∗ = (x∗1 , x∗2 , . . . , x∗n ) ∈ Rn qui vérifie les n équations non linéaires suivantes :

f1 (x1 , x2 , . . . , xn ) = 0 

f2 (x1 , x2 , . . . , xn ) = 0 



.. ≡

f (x) = 0
. 
fn−1 (x1 , x2 , . . . , xn ) = 0 



fn (x1 , x2 , . . . , xn ) = 0


où f : Rn → Rn et les fi sont des fonctions réelles de n variables que nous supposons
différentiables.

Remarque 32 Les méthodes de résolution des systèmes non linéaires sont nombreuses. Nous
ne présentons ici que les deux méthodes les plus importantes et les plus utilisées, soient la
méthode de Newton et la méthode de points fixes.

5.2.1 Méthode de Newton

Considérons le système suivant (n = 2) :

f1 (x1 , x2 ) = 0
f2 (x1 , x2 ) = 0

Soit x(0) = (x01 , x02 ), une approximation initiale de la solution de système x∗ = (x∗1 , x∗2 ).
(0) (0)
Le but est de déterminer une correction (δx1 , δx2 ) = (x∗1 − x01 , x∗2 − x02 ) de telle sorte que :
(0) (0)
f1 (x01 + δx1 , x02 + δx2 ) = 0
(0) (0)
f2 (x01 + δx1 , x02 + δx2 ) = 0

43
 (0) (0)
Le développement de Taylor autour du point (x01 , x02 ) permet de déterminer (δx1 , δx2 ) :
∂f1 (0) ∂f1 (0)
0 = f1 (x01 , x02 ) + 0 0
∂x1 (x1 , x2 )δx1 + 0 0
∂x2 (x1 , x2 )δx2 + ...
∂f2 (0) ∂f2 (0) (5.6)
0 = f2 (x01 , x02 ) + 0 0
∂x1 (x1 , x2 )δx1 + 0 0
∂x2 (x1 , x2 )δx2 + ...

Ce document ne remplace pas le cours dispensé en classe

En négligeant les termes d’ordre supérieur, (5.6) s’écrit :
∂f1 0 0 (0) ∂f1 0 0 (0)
(x , x )δx + (x , x )δx = −f1 (x01 , x02 )
∂x1 1 2 1 ∂x2 1 2 2
∂f2 0 0 (0) ∂f2 0 0 (0)
(x , x )δx + (x , x )δx = −f2 (x01 , x02 )
∂x1 1 2 1 ∂x2 1 2 2
ou encore sous forme matricielle :
" #" #
∂f1 0 , x0 ) ∂f1 (x0 , x0 ) (0)
f1 (x01 , x02 )
 
∂x1 (x1 2 ∂x2 1 2 δx1
∂f2 0 , x0 ) ∂f2 (x0 , x0 ) =−
∂x1 (x1 2 ∂x2 1 2 δx
(0)
2
f2 (x01 , x02 )
Ce système linéaire peut s’écrire aussi sous une forme plus compacte :

Jf (x(0) )δx(0) = −f (x(0) )

(0) (0) ∂fi
avec f = (f1 , f2 )T , x(0) = (x01 , x02 ) , δx(0) = (δx1 , δx2 )T et (Jf (x(0) ))ij = ( ∂xj
)(x(0) ), i, j =
1, 2.
Jf (x(0) ) est la matrice jacobienne évaluée au point x(0) = (x01 , x02 )T .
Son déterminant qu’on appelle le jacobien doit être non nul. En posant :
 0
x1 + δx01

x(1) =
x02 + δx02

et en itérant ce procédé, on définit une suite x(k) qui peut converger vers x∗ .

De manière plus générale, on peut formuler la méthode de Newton comme suit :

Pour x(0) donné, et pour k = 0, 1, 2, . . .

On résout : Jf (x(k) )δx(k) = −F(x(k) ) (5.7)
et on pose : x(k+1) = x(k) + δx(k) .

Théorème 20 Soit f : Rn → Rn une fonction différentiable sur un convexe ouvert D ∈ Rn

contenant x∗ . Supposons que Jf−1 (x∗ ) existe et il existe trois constantes positives R, L et C
telles que :

||Jf−1 (x∗ )|| ≤ C

||Jf (x) − Jf (y)|| ≤ L||x − y||, ∀x, y ∈ B(x∗ , R)

(les normes matricielle et vectorielle sont consistantes). Alors, il existe r > 0 tel que ∀x(0) ∈
B(x∗ , r), la suite définie en ( 5.7) converge vers x∗ et on a :

||x(k+1) − x∗ || ≤ CL||x(k) − x∗ ||2

Remarque 33 .
*) La convergence de la méthode de Newton dépend du choix de la valeur initiale x(0) .
*) La convergence est généralement quadratique.
*) À chaque itération, on décompose la matrice jacobienne de taille n × n.

44
Exemple 5.2.1 Considérons le système de deux équations suivants :

e x1 − x 2 = 0
x21 + x22 − 16 = 0

Ce document ne remplace pas le cours dispensé en classe

On commence par trouver la matrice jacobienne qui est dans ce cas :
 x 
e 1 −1
2x1 2x2

Graphiquement, on peut voir qu’il y’a deux solutions à ce système. La première solution se
situe près du point (−4, 0) et la deuxième près de (2.8, 2.8). Posons donc x(0) = (2.8, 2.8).
En résolvant le système :
 2.8 " (0))
#
e2.8 − 2.8
 
e −1 δx1
= −
2(2.8) 2(2.8) (0))
δx2 (2.8)2 + (2.8)2 − 16

dont la solution est δx(0) = (−0.77890, 0.83604)T on trouve le terme x(1) = (2.0211, 3.63604)T .

De la même manière, on trouve x(2) et ainsi de suite.

5.2.2 Méthode de points fixes

Soient F : Rn → Rn , x∗ ∈ Rn un point fixe de F (F(x∗ ) = x∗ ) et (x(k) )k une suite définie
par :

x(0) donné
(5.8)
x(k+1) = F(x(k) ).

Définition 29 Une application F : D ⊂ Rn → Rn est dite contractante sur un D0 ⊂ D s’il

existe 0 < α < 1 tel que :

||F(x) − F(y)|| ≤ α||x − y||, ∀x, y ∈ D0

Théorème 21 Soit F : D ⊂ Rn → Rn une contraction sur l’ensemble fermé D0 ⊂ D et

F(x) ∈ D0 pour tout x ∈ D0 . Alors F a un seul point fixe dans D0 .

Proposition 3 Supposons que F : D ⊂ Rn → Rn a un point fixe x∗ ∈ D et que F est de classe

C 1 dans un voisinage de x∗ . Si ρ(JF (x∗ )) < 1, alors il existe un voisinage S ⊂ D de x∗ tel que
∀x(0) ∈ S la suite ( 5.8) converge vers x∗ .

45
5.3 Exercices
Exercice 1
On cherche à résoudre l’équation
2

Ce document ne remplace pas le cours dispensé en classe

2
e−x = .
3(x + 1)
a) Montrez qu’il existe une solution dans l’intervalle [0, 3].
b) Faites une itération de chacune des méthodes suivantes pour trouver une des solutions
du problème dans l’intervalle [0, 3].
(i) méthode de la bissection
(ii) méthode de la sécante
(iii) méthode de Newton
2
c) Soit f (x) = e−x − 3(x+1)
2
. Sachant que f (x) < 0, f 0 (x) > 0 et f 00 (x) < 0 pour x ≥ 2,
montrez que la méthode de Newton appliquée à f ne peut pas converger si x(0) ≥ 2.

Exercice 2
Pour chaque valeur du paramètre s, on a une fonction

gs (x) = s3 + x(1 − s) + (x − s2 )3 .

a) Montrez que pour chaque valeur de s, x = s2 est un point fixe de gs .

b) La méthode des points fixes converge vers le point fixe x = s2 pour quelles valeurs de
s ? Si la convergence est linéaire, donnez le taux de convergence. Si la convergence est
plus que linéaire, donnez l’ordre de convergence.

Exercice 3
Considérez l’application qui dépend du paramètre k :
 
2ky − x + 2
G(x, y) = .
x/k − y + 1

a) Montrez que (x∗ , y ∗ ) = (k + 2, k+1

k ) est un point fixe de l’application.
b) Soit une approximation initiale (x(0) , y (0) ) suffisamment proche du point fixe (x∗ , y ∗ ).
Est-ce que la méthode des points fixes convergera et si oui, pour quelles valeurs du
paramètre k ?

Exercice 4 Soit le système d’équations non linéaire

u2 +2v+v 2
 
d(u,v)
G(~x) = G(u, v) =  
2u2 −2u+v 2
d(u,v)

où  
u
~x = et d(u, v) = 5u2 + 5v 2 − 8u + 4v + 4.
v

a) Vérifiez que p~ = [0, 0]T et ~q = [1, 0]T sont des points fixes de G.
b) À partir de ~x(0) = [0.5, 0]T , faites une itération de la méthode des points fixes. Vers quel
point fixe s’approche-t-on ?

46
c) Un long calcul montre que les jacobiennes de G évaluées aux deux points fixes sont
   
? 1/2 0 −2
JG (~
p) = et JG (~q) = .
−1/2 0 2 0

Ce document ne remplace pas le cours dispensé en classe

Faites un calcul détaillé de la composante inconnue de JG (~ p).
d) Si ~x(0) est suffisament proche de p~, est-ce que la méthode des points fixes convergera
vers p~ ? Répondez à la même question pour le point fixe ~q ? Justifiez vos réponses.

47