Méthodes numériques

PREMIERE PARTIE

Méthodes Numériques

1
Méthodes numériques

Programme Méthodes numériques

Chapitre 1 : Résolution des équations non linéaires f(x)=0

1. Introduction générale,
2. Méthode de bissection,
3. Méthode des approximations successives (point fixe),
4. Méthode de Newton-Raphson.

Chapitre 2 : Interpolation polynomiale

1. Introduction générale,
2. Polynôme de Lagrange,
3. Polynômes de Newton.

Chapitre 3 : Intégration numérique

1. Introduction générale,
2. Méthode du trapèze,
3. Méthode de Simpson,
4. Formules de quadrature.

Chapitre 4 : Résolution des équations différentielles ordinaires (problème de la

condition initiale ou de Cauchy).
1. Introduction générale,
2. Méthode d’Euler,
3. Méthode d’Euler améliorée,
4. Méthode de Runge-Kutta d’ordre quatre.

2
Méthodes numériques

Chapitre I : Résolution numérique des équations non

linéaires à une seule variable

Dans ce chapitre on va étudier trois méthodes pour la résolution des équations

non linéaires à une variable aussi dites équations à variable transcendante. Comme
exemple de ces équations, on peut citer 𝐟(𝐱) = 𝐬𝐢𝐧(𝒙) + 𝒙 = 𝟎, 𝐟(𝐱) = 𝐥𝐧(𝒙) − 𝟐𝒙 + 𝟑 = 𝟎. Ces
équations ne possèdent pas une ou des racines exactes qui peuvent être calculées
directement, c’est pourquoi on fait recours aux méthodes numériques pour trouver
les solutions approchées de ces équations. Les racines calculées sont autant précises
que l’on veut surtout lorsqu’on dispose des moyens de calcul.
Ces méthodes numériques permettent seulement le calcul d’une seule racine
sur un intervalle bien choisi. Donc si l’équation possède plus d’une racine, il est
nécessaire de les localiser dans des intervalles choisis soigneusement et de faire le
calcul pour chaque racine à part.

1 Localisation des racines d’une équation f(x)=0.

Soit une équation f(x)=0 dont on cherche la solution sur un intervalle [a,b], on
commence par un tracé grossier de la fonction sur l’intervalle donnée puis on isole
chaque racine dans un sous intervalle le plus étroit possible. La Fig. 1 montre le tracé
d’une fonction f qui coupe l’axe des x en trois points , c’est-à-dire que l’équation
f(x)=0 possède trois racines, on note les racines exates par 𝒙
̅𝟏 , 𝒙 ̅𝟑 .
̅𝟐 et 𝒙

1 x3
x1 x2
f(x)
a b a' b' a'' b''
0

-1

-2
-2 -1 0 x 1 2 3

Fig. 1.1. Localisation des racines

3
Méthodes numériques

On remarque que la fonction est continue sur chaque sous intervalle, aussi chaque
sous intervalle :

• Englobe une seule racine tel que 𝒙

̅𝟏 ∈ [𝒂 , 𝒃 ], 𝒙
̅𝟐 ∈ [𝒂′, 𝒃′]et 𝒙
̅𝟑 ∈ [𝒂′′, 𝒃′′].
• Vérifie la condition 𝒇(𝒂)𝒇(𝒃) < 𝟎, 𝒇(𝒂′)𝒇(𝒃′) < 𝟎 𝐞𝐭 𝒇(𝒂′′)𝒇(𝒃′′) < 𝟎.

La forme de l’équation f(x)=0 peut etre compliquée, dans ce cas s’il est possible on
peut la décomposer en deux parties simples g(x)=h(x).

Par exemple 𝒇(𝒙) = 𝐥𝐧(𝒙) − 𝒙𝟐 + 𝟐 = 𝟎, qui est assez compliqué pour le tracé, peut
être décomposée en :

𝒈(𝒙) = 𝒉(𝒙) avec 𝒈(𝒙) = 𝐥𝐧(𝒙) 𝐞𝐭 𝒉(𝒙) = 𝒙𝟐 − 𝟐

Le tracé de g et h est très simple, la solution de f(x)=0 se situe à l’intersection de g

et h. Ensuite, on projete les points des intersection sur l’axe des x et on localise les
racines. On peut facilement vérifier les intervalles trouvées en calculant 𝒇(𝒂𝒊 )𝒇(𝒃𝒊 ) < 𝟎.

Exemple :
Prenons l’équation : 𝐥𝐧(𝒙) − 𝒙𝟐 + 𝟐 = 𝟎 localisant ses racines, à ce point on ne connait
pas le nombre de racines de cette équation. On trace la courbe de la fonction :
𝒇(𝒙) = 𝐥𝐧(𝒙) − 𝒙𝟐 + 𝟐 = 𝟎, les intersections de la courbe avec l’axe des x représentent
les racines de cette fonction.
2 2
f2(x)=x2-2

1 1

x1 x2 x1 f1(x)=ln(x) x2
f(x)
f(x)

0 0

-1 -1

-2 -2
0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0
x x

Fig. 1.2. Tracé de la fonction f Fig. 1.3. Tracés des fonction f1 et f2

On note que cette équations a deux racines (Fig. 2) qui appartiennent par exemple
aux intervalles [0.1,0.5] et [1,2]. On peut aussi rééerire la fonction f(x)=0 sous une

forme plus simple, par exemple : 𝐥𝐧(𝒙) − 𝒙𝟐 + 𝟐 = 𝟎 ⇔ 𝐥𝐧(𝒙) = 𝒙𝟐 − 𝟐 ou bien

4
Méthodes numériques

𝑓1 (𝑥) = 𝑓2 (𝑥) ; avec 𝑓1 (𝑥) = 𝐥𝐧(𝒙) et 𝑓2 (𝑥) = 𝒙𝟐 − 𝟐 . Ensuite on trace ces deux fonctions
(Fig. 3), qui sont faciles sur les même axes, leurs intersections représentent les
racines de f(x)=0.
2 Méthode de bissection où de dichothomie
C’est la méthode la plus simple et qui nécessite le plus de calculs, elle est basée
sur le cas particulier (Théorème de Bolzano) du théorème des valeurs intermédiaires
qui dit que :

1. Si f(x) est continue sur l’intervalle [a,b],

2. Si f(a) et f(b) ne sont pas de même signe, il existe au moins un réel c compris
entre a et b tel que f(c) = 0 (car 0 est compris entre f(a) et f(b)).
2.1 Principe de la méthode
Une fois les racines sont localisées chacune dans un intervalle, pour
simplifier l’écriture soit par exemple [a,b] :
𝒂+𝒃
1. On divise l’intervalle en deux parties égales tel que 𝒙𝟎 = 𝟐
.

2. On obtient deux sous intervalles [a, 𝒙𝟎 ] et [𝒙𝟎 ,b] , la racine 𝑥̅ doit

obligatoirement appartenir à l’un d’eux. Pour vérifier, on calcule
𝒇(𝒂)𝒇(𝒙𝟎 ) et 𝒇(𝒙𝟎 )𝒇(𝒃) le produit qui est négatif c’est celui qui
correspond à l’intervalle qui contient la solution.

Notons le nouvel intervalle par [𝒂𝟏 , 𝒃𝟏 ]tel que (Fig. 4):

̅ ∈ [𝒂, 𝒙𝟎 ]
𝒂 𝒔𝒊 𝒙 ̅ ∈ [𝒂, 𝒙𝟎 ]
𝒙𝟎 𝒔𝒊 𝒙
𝒂𝟏 = { 𝑒𝑡 𝒃𝟏 = {
̅ ∈ [𝒙𝟎 , 𝒃]
𝒙𝟎 𝒔𝒊 𝒙 ̅ ∈ [𝒙𝟎 , 𝒃]
𝒃 𝒔𝒊 𝒙

b1- a1
a1 b2- a2 b1
2
a2 b2
b3- a3
a3 b3
1
f(x) a4 b4
x1 x x2
0
x3

-1

-2
0 x 1 2

Fig. 1.4. Illustration de la méthode de bissection

5
Méthodes numériques

En répétant (itérant) la même méthode pour l’intervalle obtenu on aura les valeurs :

𝒂𝟏 +𝒃𝟏 𝒂𝟐 +𝒃𝟐 𝒂𝒏 +𝒃𝒏

𝒙𝟏 = 𝟐
, 𝒙𝟐 = 𝟐
,………… , 𝒙𝒏 = 𝟐

La suite {𝒙𝒏 }𝒏=𝟎,∞ converge vers la solution 𝒙

̅ de f(x)=0 lorsque n.

2.2 Nombre de divisions pour avoir une précision ε donné.

Puisque chaque fois on divise l’intervalle en deux parties égales, on a :

𝒃−𝒂
𝒃𝟏 − 𝒂𝟏 = ;
𝟐
𝒃𝟏 − 𝒂𝟏 𝟏 𝒃 − 𝒂 𝒃 − 𝒂
𝒃𝟐 − 𝒂𝟐 = = = ;
𝟐 𝟐 𝟐 𝟐𝟐
𝒃𝟐 − 𝒂𝟐 𝟏 𝒃 − 𝒂 𝒃 − 𝒂
𝒃𝟑 − 𝒂𝟑 = = = ;
𝟐 𝟐 𝟐𝟐 𝟐𝟑

………………………………………….………………

𝒃𝒏−𝟏 − 𝒂𝒏−𝟏 𝟏 𝒃 − 𝒂 𝒃 − 𝒂
𝒃𝒏 − 𝒂𝒏 = = = .
𝟐 𝟐 𝟐𝒏−𝟏 𝟐𝒏

𝑎𝑛 𝑥𝑛 𝑥̅ 𝑏𝑛
𝟏 𝒃−𝒂 𝒃−𝒂
Puisque 𝒙
̅ ∈ [𝒂𝒏 , 𝒃𝒏 ] = [𝒂𝒏 , 𝒙𝒏 ] ∪ [𝒙𝒏 , 𝒃𝒏 ] on a |𝒙𝒏 − 𝒙
̅| ≤
𝟐 𝟐𝒏
= 𝟐𝒏+𝟏

Il faut que la différence |𝒙𝒏 − 𝒙

̅| qui est l’erreur du calcul soit inférieure à une
précision donnée ε, c’est-à-dire :

|𝒙𝒏 − 𝒙
̅| ≤ 𝜺
𝒃−𝒂
Alors, il suffit que 𝟐𝒏+𝟏
≤𝜺

𝒃−𝒂
𝐥𝐧( )
Cela donne 𝒏≥ 𝐥𝐧(𝟐)
𝟐𝜺

Le nombre de divisions ne dépend que de la longueur de l’intervalle est de la

précision.

Cette méthode est inconditionnellement convergente, son problème c’est

qu’elle est lente c’est pourquoi elle est utilisée pour démarrer d’autres méthodes plus
élaborées.

6
Méthodes numériques

Exemple : Calculons la première racine de l’équation 𝐥𝐧(𝒙) − 𝒙𝟐 + 𝟐 = 𝟎 qui appartient

à [0.1, 0.5] avec une précision de 0.01.

Calculons le nombre de divisions à faire :

𝒃−𝒂 𝟎.𝟓−𝟎.𝟏
𝐥 𝐧( ) 𝐥 𝐧( )
𝒏≥ 𝐥 𝐧(𝟐)
𝟐𝜺
= 𝟐∗𝟎.𝟎𝟏
𝐥 𝐧(𝟐)
= 𝟒. 𝟑𝟐 on prend n=5 puisque n est entier et supérieur à 4.32.

𝒇(𝒂𝟏 ) = 𝒇(𝟎. 𝟏) = −𝟎. 𝟑𝟏𝟑 et 𝒇(𝒃𝟏 ) = 𝒇(𝟎. 𝟓) = 𝟏. 𝟎𝟓𝟕

𝒂𝟏 + 𝒃𝟏 𝟎. 𝟏 + 𝟎. 𝟓
𝒙𝟏 = = = 𝟎. 𝟑𝟎; 𝒇(𝟎. 𝟑) = 𝟎. 𝟕𝟎𝟔 > 𝟎 donc 𝒂𝟐 = 𝟎. 𝟏 et 𝒃𝟐 = 𝟎. 𝟑
𝟐 𝟐
𝒂𝟐 + 𝒃𝟐 𝟎. 𝟏 + 𝟎. 𝟑
𝒙𝟐 = = = 𝟎. 𝟐𝟎; 𝒇(𝟎. 𝟐) = 𝟎. 𝟑𝟓𝟏 > 𝟎 donc 𝒂𝟑 = 𝟎. 𝟏 et 𝒃𝟑 = 𝟎. 𝟐
𝟐 𝟐
𝒂𝟑 + 𝒃𝟑 𝟎. 𝟏 + 𝟎. 𝟐
𝒙𝟑 = = = 𝟎. 𝟏𝟓; 𝒇(𝟎. 𝟏𝟓) = 𝟎. 𝟎𝟖𝟎 > 𝟎 donc 𝒂𝟒 = 𝟎. 𝟏 et 𝒃𝟒 = 𝟎. 𝟏𝟓
𝟐 𝟐
𝒂𝟒 + 𝒃𝟒 𝟎. 𝟏 + 𝟎. 𝟏𝟓
𝒙𝟒 = = = 𝟎. 𝟏𝟐𝟓; 𝒇(𝟎. 𝟏𝟐𝟓) = −𝟎. 𝟎𝟗𝟓 < 𝟎 donc 𝒂𝟓 = 𝟎. 𝟏𝟓 et 𝒃𝟒 = 𝟎. 𝟏𝟐𝟓
𝟐 𝟐
𝒂𝟓 +𝒃𝟓 𝟎.𝟏𝟓+𝟎.𝟏𝟐𝟓
𝒙𝟓 = = = 𝟎. 𝟏𝟑𝟕𝟓; 𝒇(𝟎. 𝟏𝟐𝟓) = −𝟎. 𝟎𝟎𝟑 La solution est 𝒙𝟓 = 𝟎. 𝟏𝟑𝟕𝟓
𝟐 𝟐

7
Méthodes numériques

3 Méthode des approximations successives où du point fixe.

Soit g une fonction définie sur un intervalle [a,b], le point 𝒙

̅ qui vérifie 𝒙
̅ = 𝒈(𝒙
̅)
avec 𝒙
̅ ∈ [𝒂, 𝒃] est dit point fixe de la fonction g.

Cette méthode est basée sur le principe du point fixe d’une fonction, on écrit
l’équation f(x)=0 sous la forme x=g(x), ensuite on cherche le point fixe 𝒙
̅ de la fonction
g. Pour cela on crée la suite 𝒙𝒏+𝟏 = 𝒈(𝒙𝒏 ) ( n=0,1,2….) avec 𝑥0 donnée par dichotomie
par exemple.

On démarre de 𝒙𝟎 pour n=0, on calcule 𝒙𝟏 = 𝒈(𝒙𝟎 ) ensuite n=1, on calcule 𝒙𝟐 =

𝒈(𝒙𝟏 ),..., 𝒙𝒏+𝟏 = 𝒈(𝒙𝒏 ). Sous certaines conditions, la suite {𝒙𝒏 }𝒏=𝟎,∞ converge vers la
solution 𝒙
̅ point fixe de g et solution de l’équation f(x)=0.

Exemple : Ecrire l’équation f(x)=0 sous la forme x = g(x) si 𝒇(𝒙) = 𝒙𝟐 + 𝟑𝒆𝒙 − 𝟏𝟐.

On peut écrire 𝒙 = 𝒈𝟏 (𝒙) = 𝒙𝟐 + 𝟑𝒆𝒙 − 𝟏𝟐 + 𝒙

𝒙 = 𝒈𝟐 (𝒙) = √𝟏𝟐 − 𝟑𝒆𝒙

𝟏𝟐−𝒙𝟐
𝒙 = 𝒈𝟑 (𝒙) = 𝒍𝒏 ( )
𝟑

Pour pouvoir choisir la forme de g adéquate pour le calcul, un critère de

convergence de cette méthode doit être vérifié.

3.1 Critère de convergence et d’arrêt de calculs pour la méthode des

approximations successives.
Soit g une fonction dérivable définie de [a,b] → [a,b] tel que (condition suffisante) :

|𝒈′(𝒙)| ≤ 𝒌 < 𝟏 ∀ 𝒙 ∈ [𝒂, 𝒃]

Alors la suite {𝒙𝒏 }𝒏=𝟎,∞ définie par 𝒙𝒏+𝟏 = 𝒈(𝒙𝒏 ) ( n=0,1,2….) converge
indépendamment de la valeur de 𝒙𝟎 vers l’unique point fixe 𝒙
̅ de g.

Si plusieurs formes de g vérifient cette condition, on aura plusieurs valeurs de

k. On choisit celle avec la valeur minimale de k. En pratique, on calcule 𝒌 =
𝒎𝒂𝒙𝒙∈[𝒂,𝒃] |𝒈′(𝒙)| qui doit être inférieure à l’unité pour que la méthode converge.

On arrête les calculs pour cette méthode lorsque la différence absolue entre
deux itérations successives est inférieure à une certaine précision 𝜀 donnée.
|𝒙𝒏+𝟏 − 𝒙𝒏 | < 𝜺

8
Méthodes numériques

Exemple : Trouver la première racine de l’équation 𝐥𝐧(𝒙) − 𝒙𝟐 + 𝟐 = 𝟎 qui

appartient à [0.1, 0.5] avec une précision ε=0.001. On écrit cette équation sous la
forme x = g(x) et on vérifie les conditions de convergence. On peut écrire :
𝟐 −𝟐
𝒙 = 𝒆𝒙 = 𝒈𝟏 (𝒙) et
𝒙 = √𝒍𝒏(𝒙) + 𝟐 = 𝒈𝟐 (𝒙)
Vérifions la condition de convergence pour cette méthode 𝒌 = 𝒎𝒂𝒙𝒙∈[𝒂,𝒃] |𝒈′(𝒙)|
𝟐 −𝟐
𝒌𝟏 = 𝒎𝒂𝒙𝒙∈[𝟎.𝟏,𝟎.𝟓] |𝒈𝟏 ′(𝒙)| = 𝒎𝒂𝒙𝒙∈[𝟎.𝟏,𝟎.𝟓] |𝟐𝒙𝒆𝒙 | on a 𝒈𝟏 ′(𝒙) strictement croissante
𝟐 −𝟐
donc 𝒌𝟏 = 𝒎𝒂𝒙𝒙=𝟎.𝟓 |𝟐 ∗ 𝟎. 𝟓𝒆𝟎.𝟓 | = 𝟎. 𝟏𝟕𝟒 < 𝟏 cette forme converge.
𝟐
On écrit : 𝒙𝒏+𝟏 = 𝒈𝟏 (𝒙𝒏 ) = 𝒆𝒙𝒏 −𝟐 (n=0,1,2,…..)
Commençons x0=0.3 le milieu de l’intervalle initial donné :
𝟐
𝒏 = 𝟎, 𝒙𝟏 = 𝒈𝟏 (𝒙𝟎 ) = 𝒆𝒙𝟎−𝟐 = 0.148
On calcule |𝑥1 − 𝑥0 | = 0.152 > 𝜀 ;
𝟐
On continue 𝒏 = 𝟏, 𝒙𝟐 = 𝒈𝟏 (𝒙𝟏 ) = 𝒆𝒙𝟏−𝟐 = 0.138.
On calcule |𝑥2 − 𝑥1 | = 0.01 > 𝜀
𝟐
On continue 𝒏 = 𝟐, 𝒙𝟑 = 𝒈𝟏 (𝒙𝟐 ) = 𝒆𝒙𝟐−𝟐 = 0.138.
On calcule |𝑥3 − 𝑥2 | = 0.00 < 𝜀, La solution est 𝒙𝟐 = 𝟎. 𝟏𝟑𝟖

9
Méthodes numériques

4 Méthode de Newton-Raphson
C’est la méthode la plus efficace et la plus utilisée, elle repose sur le
développement de Taylor. Si f(x) est continue et continument dérivable dans le
voisinage de 𝒙
̅ solution de f(x)=0, alors le développement en série de Taylor autour
d’un estimé 𝒙𝒏 proche de 𝒙
̅ s’écrit :

(𝒙
̅−𝒙𝒏 ) ′ ̅−𝒙𝒏 )𝟐 ′′
(𝒙
̅) = 𝒇(𝒙𝒏 ) +
𝒇(𝒙 𝒇 (𝒙𝒏 ) + 𝒇 (𝒙𝒏 ) + ⋯.
𝟏! 𝟐!

Si 𝒙𝒏 est un estimé proche de 𝒙

̅, alors le carré de l’erreur 𝜺𝒏 = 𝒙
̅ − 𝒙𝒏 et les termes de
degrés supérieurs sont négligeable. Sachant que 𝒇(𝒙
̅) = 𝟎 on obtient la relation
approximative :

̅ − 𝒙𝒏 )𝒇′ (𝒙𝒏 ) ≈ 𝟎
𝒇(𝒙𝒏 ) + (𝒙

𝒇(𝒙𝒏 )
Donc ̅ = 𝒙𝒏 −
𝒙
𝒇′ (𝒙𝒏 )

On peut écrire la (𝒏 + 𝟏)𝒆𝒎𝒆 itération approximant 𝑥̅ est :

𝒇(𝒙 )
𝒙𝒏+𝟏 = 𝒙𝒏 − 𝒇′ (𝒙𝒏 ) (n=0,1,2,…..)
𝒏

Cette suite, si elle converge, doit converger vers la solution 𝒙

̅ de f(x)=0. On remarque
que f’(x) doit être non nulle.

4.1 Critère de convergence de la méthode de Newton-Raphson

Soit une fonction f définie sur [a,b] telle que :

i. 𝒇(𝒂)𝒇(𝒃) < 𝟎
ii. 𝒇′(𝒙) 𝒆𝒕 𝒇′′(𝒙) sont non nulles et gardent un signe constant sur
l’intervalle donné.

4.2 Critère d’arrêt de calcul pour la méthode de Newton-Raphson

Si la condition de convergence est vérifiée, le procédé itératif doit converger.
Cela veut dire que chaque nouvelle itération est meilleure que la précédente, de ce
fait on peut dire que si on a une précision 𝜺, on arrête le calcul lorsque la différence
absolue entre deux approximations successives est inférieure à la précision donnée.
C’est-à-dire :

|𝒙𝒏+𝟏 − 𝒙𝒏 | ≤ 𝜺
Si cette condition est vérifiée on prend 𝒙𝒏+𝟏 comme solution de f(x)=0.

10
Méthodes numériques

Exemple : Trouver la première racine de l’équation 𝒇(𝒙) = 𝐥𝐧(𝒙) − 𝒙𝟐 + 𝟐 = 𝟎

qui appartient à [0.1, 0.5] avec une précision ε=0.0001. On calcule la dérivée
première et seconde de f et on vérifie les conditions de convergence.
𝟏
On a : 𝒇′(𝒙) = 𝒙 − 𝟐𝒙 qui est strictement décroissante et positive sur l’intervalle donné.
𝟏
𝒇′(𝒙) > 𝟎 et 𝒇′′ (𝒙) = − 𝒙𝟐 − 𝟐, 𝒇′′(𝒙) < 𝟎 sur l’intervalle donné. La condition de
𝒇(𝒙 ) 𝐥𝐧(𝒙𝒏 )−𝒙𝒏 𝟐 +𝟐
convergence est vérifiée. On écrit donc : 𝒙𝒏+𝟏 = 𝒙𝒏 − 𝒇′ (𝒙𝒏 ) = 𝒙𝒏 − 𝟏 , (n=0,1,2,…..).
𝒏 −𝟐𝒙𝒏
𝒙𝒏

Commençons x0=0.3 le milieu de l’intervalle initial donné :

𝐥𝐧(𝒙𝟎 ) − 𝒙𝟎 𝟐 + 𝟐
𝒏 = 𝟎, 𝒙𝟏 = 𝒙𝟎 − = 𝟎. 𝟎𝟒𝟏𝟕
𝟏
𝒙𝟎 − 𝟐𝒙 𝟎

On calcule |𝑥1 − 𝑥0 | > 𝜀 ;

On continue 𝒏 = 𝟏, 𝒙𝟐 = 0.0910.
On calcule |𝑥2 − 𝑥1 | => 𝜀
On continue 𝒏 = 𝟐, 𝒙𝟑 = 0.1285.
On calcule |𝑥3 − 𝑥2 | > 𝜀 ;
On continue 𝒏 = 𝟑, 𝒙𝟒 = 0.1376.
On calcule |𝑥4 − 𝑥3 | => 𝜀
On continue 𝒏 = 𝟒, 𝒙𝟓 = 0.1379.
On calcule |𝑥5 − 𝑥4 | => 𝜀
On continue 𝒏 = 𝟓, 𝒙𝟔 = 𝟎. 𝟏𝟑𝟕𝟗. La solution est 𝒙𝟒 = 𝟎. 𝟏𝟑𝟕𝟗
Remarques :

• La méthode de bissection est inconditionnellement convergente, son

inconvénient est sa lenteur pour obtenir la solution avec une grande
précision. Elle peut servir pour démarrer d’autres méthodes plus
performantes.
• La méthode des approximations successives est plus rapide que celle de
bissection à condition qu’elle converge.
• La méthode de Newton-Raphson est la plus rapide, elle permet d’obtenir
des solutions très précises en un nombre réduit d’itérations.

11
Méthodes numériques

Chapitre II : Interpolation polynômiale

Soit par exemple une expérience où on enregistre la distance parcourue par

un objet en fonction du temps, les résultats sont donnés dans le tableau suivant :
t(sec) 0 1 2 3 4
X(m) 0 5 15 0 3

On veut par exemple calculer la position de l’objet au temps t=2.5 sec ou la

vitesse de l’objet à un temps donné. Pour cela, il faut avoir une forme analytique de
x en fonction de t, X(t). Cette forme doit au moins coïncider avec les points donnés
𝟒 𝒅𝑿(𝒕)
dans le tableau. Ensuite on peut calculer X(2.5), ∫𝟎 𝑿(𝒕)𝒅𝒕 où bien 𝒗(𝒕) = 𝒅𝒕
.

Dans ce chapitre, on va considérer l’approximation de X(t) par une forme

polynômiale c’est-à-dire :
𝑿(𝒕) = 𝒂𝟎 + 𝒂𝟏 𝒕 + 𝒂𝟐 𝒕𝟐 + ⋯ + 𝒂𝒏 𝒕𝒏
Avec 𝒂𝒊 (𝒊 = 𝟎, 𝒏) sont des coefficients à déterminer.
Les polynômes que nous allons étudier différent seulement par la façon de
déterminer les coefficients 𝒂𝒊 (𝒊 = 𝟎, 𝒏), car pour un tableau de valeurs données le
polynôme d’interpolation est unique.

2.1 Polynôme d’interpolation de Lagrange.

Soient (n+1) points distincts 𝒙𝟎 , 𝒙𝟏 , … … , 𝒙𝒏 et f une fonction dont les valeurs

sont 𝒇(𝒙𝟎 ), 𝒇(𝒙𝟏 ), … … , 𝒇(𝒙𝒏 ). Alors, il existe un seul polynôme de degré inférieur ou égal
à n et qui coïncide avec les points d’interpolation, i e :

𝒇(𝒙𝒌 ) = 𝑷𝒏 (𝒙𝒌 ), 𝒌 = 𝟎, 𝟏, 𝟐, … … , 𝒏
Ce polynôme est donné par :

𝒏
𝑷𝒏 (𝒙) = ∑ 𝒇(𝒙𝒊 )𝑳𝒊 (𝒙) = 𝒇(𝒙𝟎 )𝑳𝟎 (𝒙) + 𝒇(𝒙𝟏 )𝑳𝟏 (𝒙) + ⋯ … + 𝒇(𝒙𝒏 )𝑳𝒏 (𝒙)
𝒊=𝟎

(𝒙−𝒙𝒊 ) (𝒙−𝒙𝟎 ) (𝒙−𝒙𝟏 ) (𝒙−𝒙 ) (𝒙−𝒙 ) (𝒙−𝒙𝒏 )

Avec 𝑳𝒌 (𝒙) = ∑𝒏𝒊=𝟎,𝒊≠𝒌 = … . . (𝒙 𝒌−𝟏 ) (𝒙 𝒌+𝟏 ) … (k=0,..,n)
(𝒙𝒌 −𝒙𝒊 ) (𝒙𝒌 −𝒙𝟎 ) (𝒙𝒌 −𝒙𝟏 ) 𝒌 −𝒙𝒌−𝟏 𝒌 −𝒙𝒌+𝟏 (𝒙 𝒌 −𝒙𝒏 )

𝑳𝒌 (𝒙) sont dits coefficients polynômes de Lagrange, ils sont orthogonaux c’est-à-dire
𝑳𝒌 (𝒙𝒋 ) = 𝟎 et 𝑳𝒌 (𝒙𝒌 ) = 𝟏.

12
Méthodes numériques

Exemple : Reprenons la table donnée au début du chapitre et essayons de calculer

le polynôme de Lagrange pour cette table. Notons que pour n+1 points le degré du
polynôme est inférieure ou égal à n. Pour notre cas on a 5 points, cela nous donne
un polynôme de de gré inférieur ou égal à 4.

𝟒
𝑿(𝒕) ≈ 𝑷𝟒 (𝒕) = ∑ 𝒇(𝒕𝒊 )𝑳𝒊 (𝒕) = 𝒇(𝒕𝟎 )𝑳𝟎 (𝒕) + 𝒇(𝒕𝟏 )𝑳𝟏 (𝒕) + 𝒇(𝒕𝟐 )𝑳𝟐 (𝒕) + 𝒇(𝒕𝟑 )𝑳𝟑 (𝒕) + 𝒇(𝒕𝟒 )𝑳𝟒 (𝒕)
𝒊=𝟎

Avec les coefficients 𝑓(𝑡𝑖 ) sont les valeurs de X(ti) aux points donnés ti, on remplace
et on écrit donc :

𝑿(𝒕) ≈ 𝑷𝟒 (𝒕) = 𝟎 ∗ 𝑳𝟎 (𝒕) + 𝟓 ∗ 𝑳𝟏 (𝒕) + 𝟏𝟓 ∗ 𝑳𝟐 (𝒕) + 𝟎 ∗ 𝑳𝟑 (𝒕) + 𝟑 ∗ 𝑳𝟒 (𝒕)

Ensuite, on calcule les coefficients polynômes de Lagrange :

(𝒕−𝒕𝒊 ) (𝒕−𝒕𝟏 ) (𝒕−𝒕𝟐 ) (𝒕−𝒕𝟑 ) (𝒕−𝒕𝟒 )

𝑳𝟎 (𝒕) = ∑𝟒𝒊=𝟎,𝒊≠𝟎 (𝒕 = (𝒕 Noter bien qu’il est inutile de calculer les
𝒌 −𝒕𝒊 ) 𝟎 −𝒕𝟏 ) (𝒕𝟎 −𝒕𝟐 ) (𝒕𝟎 −𝒕𝟑 ) (𝒕𝟎 −𝒕𝟒 )

coefficients polynômes 𝑳𝟎 (𝒕) et 𝑳𝟑 (𝒕) car ils seront multipliés par zéro dans le
remplacement.

(𝒕−𝒕𝒊 ) (𝒕−𝒕𝟎 ) (𝒕−𝒕𝟐 ) (𝒕−𝒕𝟑 ) (𝒕−𝒕𝟒 ) (𝒕−𝟎) (𝒕−𝟐) (𝒕−𝟑) (𝒕−𝟒) 𝟏

𝑳𝟏 (𝒕) = ∑𝟒𝒊=𝟎,𝒊≠𝟏 = = = − (𝒕𝟒 − 𝟗𝒕𝟑 + 𝟐𝟔𝒕𝟐 − 𝟐𝟒𝒕)
(𝒕𝒌 −𝒕𝒊 ) (𝒕𝟏 −𝒕𝟎 ) (𝒕𝟏 −𝒕𝟐 ) (𝒕𝟏 −𝒕𝟑 ) (𝒕𝟏 −𝒕𝟒 ) (𝟏−𝟎) (𝟏−𝟐) (𝟏−𝟑) (𝟏−𝟒) 𝟔

(𝒕−𝒕𝒊 ) (𝒕−𝒕𝟎 ) (𝒕−𝒕𝟏 ) (𝒕−𝒕𝟑 ) (𝒕−𝒕𝟒 ) (𝒕−𝟎) (𝒕−𝟏) (𝒕−𝟑) (𝒕−𝟒) 𝟏

𝑳𝟐 (𝒕) = ∑𝟒𝒊=𝟎,𝒊≠𝟐 = = = (𝒕𝟒 − 𝟖𝒕𝟑 + 𝟏𝟗𝒕𝟐 − 𝟏𝟐𝒕)
(𝒕𝒌 −𝒕𝒊 ) (𝒕𝟐 −𝒕𝟎 ) (𝒕𝟐 −𝒕𝟏 ) (𝒕𝟐 −𝒕𝟑 ) (𝒕𝟐 −𝒕𝟒 ) (𝟐−𝟎) (𝟐−𝟏) (𝟐−𝟑) (𝟐−𝟒) 𝟒

(𝒕−𝒕𝒊 ) (𝒕−𝒕𝟎 ) (𝒕−𝒕𝟏 ) (𝒕−𝒕𝟐 ) (𝒕−𝒕𝟑 ) (𝒕−𝟎) (𝒕−𝟏) (𝒕−𝟐) (𝒕−𝟑) 𝟏

𝑳𝟒 (𝒕) = ∑𝟒𝒊=𝟎,𝒊≠𝟒 = = = ( 𝒕𝟒 − 𝟔𝒕𝟑 + 𝟏𝟏𝒕𝟐 − 𝟔𝒕)
(𝒕𝒌 −𝒕𝒊 ) (𝒕𝟒 −𝒕𝟎 ) (𝒕𝟒 −𝒕𝟏 ) (𝒕𝟒 −𝒕𝟐 ) (𝒕𝟒 −𝒕𝟑 ) (𝟒−𝟎) (𝟒−𝟏) (𝟒−𝟐) (𝟒−𝟑) 𝟐𝟒

Finalement on remplace les coefficients polynômes et on obtient :

𝑿(𝒕) ≈ 𝑷𝟒 (𝒕) = −𝟐𝟓. 𝟕𝟓𝒕 + 𝟓𝟎. 𝟗𝟓𝟖𝟑𝟑𝒕𝟐 − 𝟐𝟑. 𝟐𝟓𝒕𝟑 + 𝟑. 𝟎𝟒𝟏𝟔𝟕𝒕𝟒

13
Méthodes numériques

2.2 Polynôme d’interpolation de Newton

On a vu que le polynôme de Lagrange utilise (n+1) coefficients polynômes qui

sont eux aussi des polynômes de degré inférieur ou égal à n. Le calcul de ces
coefficients polynômes est aussi une tâche délicate, c’est pourquoi il est intéressant
d’utiliser une autre formulation plus souple ; c’est le polynôme de Newton.

Le calcul du polynôme de Newton commence par la construction d’un

polynôme de degré 1, 𝑷𝟏 (𝒙) qui passe par les deux premiers points. Ensuite, ce
dernier sera utilisé pour calculer un autre de degré 2, 𝑷𝟐 (𝒙) qui passe par les trois
premiers points et ainsi de suite jusqu’au polynôme final de degré inférieur ou égal à
n, 𝑷𝒏 (𝒙). On a la relation de récurrence suivante entre deux polynômes successifs
𝑷𝒊−𝟏 (𝒙) et 𝑷𝒊 (𝒙) (i=2,3,..,n+1):

𝑷𝟏 (𝒙) = 𝒂𝟎 + 𝒂𝟏 (𝒙 − 𝒙𝟎 )
𝑷𝟐 (𝒙) = 𝑷𝟏 (𝒙) + 𝒂𝟐 (𝒙 − 𝒙𝟎 )(𝒙 − 𝒙𝟏 )
𝑷𝟑 (𝒙) = 𝑷𝟐 (𝒙) + 𝒂𝟑 (𝒙 − 𝒙𝟎 )(𝒙 − 𝒙𝟏 )(𝒙 − 𝒙𝟐 )
………………………………………..
{𝑷𝒏 (𝒙) = 𝑷𝒏−𝟏 (𝒙) + 𝒂𝒏 (𝒙 − 𝒙𝟎 )(𝒙 − 𝒙𝟏 ) … (𝒙 − 𝒙𝒏−𝟏 )

On remarque que les coefficients 𝒂𝒌 (k=0,…,n) sont les éléments essentiels

dans le calcul des polynômes de Newton. Ces coefficients sont les différences divisées
d’ordre k de la fonction f.

2.3 Calcul des différences divisées d’une fonction f.

Les différences divisées d’une fonction f basée sur les points 𝒙𝟎 , 𝒙𝟏 , … … , 𝒙𝒏
𝒚𝒊
sont données par : 𝒂𝒌 = 𝒇[𝒙𝟎 , 𝒙𝟏 , … … , 𝒙𝒌 ] = ∑𝒌𝒊=𝟎 ∏𝒌
𝒋=𝟎,𝒋≠𝒊(𝒙𝒊 −𝒙𝒋 )

En pratique et pour le nombre limité de points, les différences divisées sont calculées
en utilisant un tableau qui a la forme suivante :

𝒙𝒌 𝒇(𝒙𝒌 ) = 𝒇[𝒙𝒌 ] DD1 DD2 DD3 DD4

𝒙𝟎 𝒇[𝒙𝟎 ]
𝒙𝟏 𝒇[𝒙𝟏 ] 𝒇[𝒙𝟎 , 𝒙𝟏 ] 𝒇[𝒙 , 𝒙 , 𝒙 ]
𝟎 𝟏 𝟐

𝒙𝟐 𝒇[𝒙𝟐 ] 𝒇[𝒙𝟏 , 𝒙𝟐 ] 𝒇[𝒙 , 𝒙 , 𝒙 ] 𝒇[𝒙𝟎 , 𝒙𝟏 , 𝒙𝟐 , 𝒙𝟑 ] 𝒇[𝒙 , 𝒙 , 𝒙 , 𝒙 , 𝒙 ]

𝟏 𝟐 𝟑 𝟎 𝟏 𝟐 𝟑 𝟒

𝒙𝟑 𝒇[𝒙𝟑 ] 𝒇[𝒙𝟐 , 𝒙𝟑 ] 𝒇[𝒙 , 𝒙 , 𝒙 ] 𝒇[𝒙𝟏 , 𝒙𝟐 , 𝒙𝟑 , 𝒙𝟒 ]

𝟐 𝟑 𝟒

𝒙𝟒 𝒇[𝒙𝟒 ] 𝒇[𝒙𝟑 , 𝒙𝟒 ]

14
Méthodes numériques ST L 2 S4 Dr. MAMERI A. Dép. GM, FSSA, ULBM

𝑓[𝑥1 ]−𝑓[𝑥0 ] 𝑓[𝑥2 ]−𝑓[𝑥1 ]

Avec : 𝑓[𝑥0 , 𝑥1 ] = , 𝑓[𝑥1 , 𝑥2 ] = , … … ..
𝑥1 −𝑥0 𝑥2 −𝑥1

𝑓[𝑥1 , 𝑥2 ] − 𝑓[𝑥0 , 𝑥1 ] 𝑓[𝑥2 , 𝑥3 ] − 𝑓[𝑥1 , 𝑥2 ]

𝑓[𝑥0 , 𝑥1, 𝑥2 ] = , 𝑓[𝑥1 , 𝑥2, 𝑥3 ] = ,……..
𝑥2 − 𝑥0 𝑥3 − 𝑥1
𝑓[𝑥1 , 𝑥2 , 𝑥3 ] − 𝑓[𝑥0 , 𝑥1 , 𝑥2 ] 𝑓[𝑥2 , 𝑥3 , 𝑥4 ] − 𝑓[𝑥1 , 𝑥2 , 𝑥3 ]
𝑓[𝑥0 , 𝑥1 , 𝑥2 , 𝑥3 ] = , 𝑓[𝑥1 , 𝑥2 , 𝑥3 , 𝑥4 ] =
𝑥3 − 𝑥0 𝑥4 − 𝑥1

𝑓[𝑥1 ,𝑥2 ,𝑥3 ,𝑥4 ]−𝑓[𝑥0 ,𝑥1 ,𝑥2 ,𝑥3 ]

𝑓[𝑥0 , 𝑥1 , 𝑥2 , 𝑥3 , 𝑥4 ] =
𝑥4 −𝑥0

2.4 Erreur d’interpolation

C’est l’erreur commise lorsqu’on remplace la fonction f par le polynôme

d’interpolation équivalent. Elle est notée par 𝜺(𝒙) car elle varie d’un point à un autre
dans l’intervalle d’interpolation. Cette erreur doit être nulle aux points
d’interpolation, 𝜺(𝒙𝒊 ) = 𝟎, (i=0,…n).

Si la fonction f est continue et (n+1) fois dérivable sur l’intervalle d’interpolation

[𝒂 = 𝒙𝟎 , 𝒃 = 𝒙𝒏 ]. Alors pour tout 𝒙 ∈ [𝒂, 𝒃] il existe 𝒛 ∈ [𝒂, 𝒃] tel que :

𝒏 (𝒙 − 𝒙𝒊 ) (𝒏+𝟏)
𝜺(𝒙) = |𝒇(𝒙) − 𝑷𝒏 (𝒙)| = ∏ 𝒇 (𝒛)
𝒊=𝟎 (𝒏 + 𝟏)!

Si |𝒇(𝒏+𝟏) (𝒛)| ≤ 𝑴 ∀ 𝒛 ∈ [𝒂, 𝒃] on peut écrire :

𝒏 (𝒙 − 𝒙𝒊 )
𝜺(𝒙) = |𝒇(𝒙) − 𝑷𝒏 (𝒙)| ≤ ∏ 𝑴
𝒊=𝟎 (𝒏 + 𝟏)!

Dans ce cas M est majorant de la fonction 𝒇(𝒏+𝟏) (𝒙) sur l’intervalle [a,b].

Exemple : Reprenons la table donnée au début du chapitre et essayons de calculer

le polynôme de Newton pour cette table. Notons que pour n+1 points le degré du
polynôme est inférieure ou égal à n. Pour notre cas on a 5 points, cela nous donne
un polynôme de de gré inférieur ou égal à 4. Ecrivons les polynômes de Newton :
𝑷𝟏 (𝒕) = 𝒂𝟎 + 𝒂𝟏 (𝒕 − 𝒕𝟎 )
𝑷𝟐 (𝒕) = 𝑷𝟏 (𝒕) + 𝒂𝟐 (𝒕 − 𝒕𝟎 )(𝒕 − 𝒕𝟏 )
𝑷𝟑 (𝒕) = 𝑷𝟐 (𝒕) + 𝒂𝟑 (𝒕 − 𝒕𝟎 )(𝒕 − 𝒕𝟏 )(𝒕 − 𝒕𝟐 )
𝑷𝟒 (𝒕) = 𝑷𝟑 (𝒕) + 𝒂𝟒 (𝒕 − 𝒕𝟎 )(𝒕 − 𝒕𝟏 )(𝒕 − 𝒕𝟐 )(𝒕 − 𝒕𝟑 )
{ Les 𝒂𝒊 sont les différences diviséees d′ ordre 𝒊

Calculons la table des différences divisées.

𝒕𝒌 𝑓(𝑡𝑘 ) = 𝑓[𝑡𝑘 ] DD1 DD2 DD3 DD4
0 𝟎 = 𝒂𝟎
𝟓 = 𝒂𝟏
1 5 𝟐. 𝟓 = 𝒂𝟐
10 −𝟓 = 𝒂𝟑
2 15 −12.5 𝟑. 𝟎𝟒𝟏𝟔𝟕 = 𝒂𝟒
−15 9 7.1667
3 0
3
4 3

15
Méthodes numériques

Remplaçant les 𝑎𝑖 et les 𝑡𝑖 par leurs valeurs dans les polynômes de Newton, on trouve :

𝑷𝟏 (𝒕) = 𝒂𝟎 + 𝒂𝟏 (𝒕 − 𝒕𝟎 ) = 𝟎 + 𝟓(𝒕 − 𝟎) = 𝟓𝒕

𝑷𝟐 (𝒕) = 𝑷𝟏 (𝒕) + 𝒂𝟐 (𝒕 − 𝒕𝟎 )(𝒕 − 𝒕𝟏 ) = 𝟓𝒕 + 𝟐. 𝟓(𝒕 − 𝟎)(𝒕 − 𝟏) = 𝟐. 𝟓𝒕𝟐 + 𝟐. 𝟓𝒕

𝑷𝟑 (𝒕) = 𝑷𝟐 (𝒕) + 𝒂𝟑 (𝒕 − 𝒕𝟎 )(𝒕 − 𝒕𝟏 )(𝒕 − 𝒕𝟐 ) = 𝟐. 𝟓𝒕𝟐 + 𝟐. 𝟓𝒕 − 𝟓(𝒕 − 𝟎)(𝒕 − 𝟏)(𝒕 − 𝟐) = −𝟓𝒕𝟑 + 𝟏𝟕. 𝟓𝒕𝟐 − 𝟕. 𝟓𝒕

𝑷𝟒 (𝒕) = 𝑷𝟑 (𝒕) + 𝒂𝟒 (𝒕 − 𝒕𝟎 )(𝒕 − 𝒕𝟏 )(𝒕 − 𝒕𝟐 )(𝒕 − 𝒕𝟑 ) = −𝟓𝒕𝟑 + 𝟏𝟕. 𝟓𝒕𝟐 − 𝟔𝒕 + 𝟑. 𝟎𝟒𝟏𝟕(𝒕 − 𝟎)(𝒕 − 𝟏)(𝒕 − 𝟐)(𝒕 − 𝟑)

𝑷𝟒 (𝒕) = 3.04167𝑡4 − 23.25𝑡3 + 50.95833𝑡2 − 25.75𝑡 C’est le même polynôme que celui de
Lagrange.

20
P2 P1
(2,15)
15

10

(1,5)
5
Pi(t)

(4,3)
(0,0)
(3,0)
0
P4
Points d'interpolation
-5

P3
-10
0 1 2 3 4
t

Fig. 2.1. Tracés des polynômes d’interpolation de Newton

16
Méthodes numériques

Chapitre III : Intégration Numérique

Dans ce chapitre, on va étudier quelques méthodes approximatives pour le

calcul des intégrales limitées. Aussi ces méthodes permettent le calcul des intégrales
qui n’ont pas de solutions directes ou analytiques. On peut aussi calculer l’intégrale
d’une fonction donnée sous forme tabulaire ou discrète.

3.1 Formule du trapèze

Cette formule est très simple, elle permet de remplacer la courbe f(x) de la
fonction à intégrer par une ligne droite qui relie les points (a,f(a)) et (b,f(b)) ce qui
donne un trapèze (Fig. 3.1).

2.0

1.5

(b,f(b))
1.0 f(x)
f

0.5 (a,f(a))

0.0
a b
h
-0.5
0.0 0.5 1.0 1.5 2.0
x

Fig. 3.1 : Méthode du trapèze

L’intégrale est donc remplacée par la surface du trapèze :

𝒃 𝒉
𝒔 = ∫𝒂 𝒇(𝒙) = 𝟐 (𝒇(𝒂) + 𝒇(𝒃))

Avec h=b-a est dit pas d’intégration. On peut remarquer qu’il y’a une différence
importante entre la courbe de la fonction et la ligne droite, cela veut qu’on commît
une erreur de calcul. Pour minimiser cette erreur, on utilise une autre forme plus
adaptée de cette formule.

17
Méthodes numériques

3.1.1 Formule du trapèze généralisée.

On divise l’intervalle [a,b] en plusieurs sous intervalles égaux et on applique la
formule du trapèze à chaque sous intervalle (Fig. 3.2). On a donc les sous intervalles
[𝒂 = 𝒙𝟎 , 𝒙𝟏 ] ∪ [𝒙𝟏 , 𝒙𝟐 ] ∪ … . .∪ [𝒙𝒏−𝟏 , 𝒙𝒏 = 𝒃], l’application de la formule du trapèze
donne :
𝑏
ℎ ℎ ℎ ℎ
∫ 𝑓(𝑥) = (𝑓(𝑥0 ) + 𝑓(𝑥1 )) + (𝑓(𝑥1 ) + 𝑓(𝑥2 )) + (𝑓(𝑥2 ) + 𝑓(𝑥3 )) + ⋯ + (𝑓(𝑥𝑛−1 ) + 𝑓(𝑥𝑛 ))
𝑎 2 2 2 2
𝑏 𝑛−1 𝑛−1
ℎ ℎ
∫ 𝑓(𝑥) = (𝑓(𝑥0 ) + 2 ∑ 𝑓(𝑥𝑖 ) + 𝑓(𝑥𝑛 )) = (𝑓0 + 2 ∑ 𝑓𝑖 + 𝑓𝑛 )
𝑎 2 𝑖=1 2 𝑖=1

Fig. 3.2. Méthode du trapèze généralisée

3.1.2 Erreur d’intégration

C’est la différence entre l’intégrale exacte de la fonction et celle calculée par la

méthode du trapèze, elle est notée par R(f).
𝒙 𝒉 𝒃−𝒂 𝟐 ′′
𝑹(𝒇) = ∫𝒙 𝒏 𝒇(𝒙)𝒅𝒙 − 𝟐 (𝒇𝟎 + 𝟐 ∑𝒏−𝟏
𝒊=𝟏 𝒇𝒊 + 𝒇𝒏 ) = − 𝟏𝟐
𝒉 𝒇 (𝒛) avec 𝒛 ∈ [𝒂, 𝒃]
𝟎

1 2
Exemple : Soit à calculer l’intégrale ∫0 𝑒 −𝑥 𝑑𝑥 avec une précision de 0.001 par la
méthode du trapèze. On doit premièrement trouver le nombre de division à faire pour
𝒃−𝒂 𝟐 ′′
obtenir cette précision. L’erreur d’intégration s’écrit 𝑹(𝒇) = − 𝟏𝟐
𝒉 𝒇 (𝒛) sa valeur
absolue doit être inférieure ou égale à la précision donnée (0.001), c.à.d. :

𝒃 − 𝒂 𝟐 ′′
|𝑹(𝒇)| = |− 𝒉 𝒇 (𝒛)| ≤ 𝟎. 𝟎𝟎𝟏
𝟏𝟐

2 2
La fonction 𝑓(𝑥) = 𝑒 −𝑥 donc sa dérivée seconde est 𝑓"(𝑥) = 2(2𝑥 2 − 1)𝑒 −𝑥 , cette
fonction est strictement croissante dans l’intervalle donné (Fig. 3.3).

18
Méthodes numériques

f"(x)
y

-1

-2
0,0 0,5 1,0
x

Fig. 3.3. Courbe de la dérivée

2
seconde de 𝑒 −𝑥

On calcule

𝑀 = 𝑚𝑎𝑥|𝒇′′ (𝒛)| = 𝟐 à 𝒙 = 𝟎

𝒃−𝒂 𝟐 ′′
donc |𝑹(𝒇)| = |−
𝟏𝟐
𝒉 𝒇 (𝒛)| ≤ 𝟎. 𝟎𝟎𝟏

𝟏𝟐∗𝟎.𝟎𝟎𝟏 𝟏
d’où 𝒉 ≤ √ (𝟏−𝟎)∗𝟐 = 𝟎. 𝟎𝟕𝟕𝟒 donc n = = 𝟏𝟐. 𝟗𝟏 on prend 13 divisions.
𝟎.𝟎𝟕𝟕𝟒

1
Le pas d’intégration ℎ = 13.

1 12
−𝑥 2
1 2 𝑖 2 2
∫ 𝑒 𝑑𝑥 = (𝑒 −0 + 2 ∑ 𝑒 −(13) + 𝑒 −1 ) = 0.74646
0 2 ∗ 13
𝑖=1

19
Méthodes numériques

3.2 Formules de Simpson

Dans cette formule on ne remplace pas la fonction par une droite mais par une
parabole de degré n inférieure ou égale à deux. Cette dernière doit passer par trois
points (𝒙𝟎 , 𝒚𝟎 ), (𝒙𝟏 , 𝒚𝟏 )𝐞𝐭 (𝒙𝟐 , 𝒚𝟐 ) ce qui fait que cette méthode n’est applicable que pour
un nombre pair de tranches (une tranche c’est l’intervalle entre deux points) (Fig.
3.4.). La formule de Simpson s’écrit :
𝒃
𝒉
∫ 𝒇(𝒙) ≅ (𝒇(𝒙𝟎 ) + 𝟒𝒇(𝒙𝟏 ) + 𝒇(𝒙𝟐 ))
𝒂 𝟑

Fig. 3.4. Méthode de Simpson

Si on généralise la formule de Simpson pour 2n sous intervalles avec un pas

𝒃−𝒂
d’intégration 𝒉= 𝟐𝒏
, 𝒂 = 𝒙𝟎 < 𝒙𝟏 < ⋯ … . < 𝒙𝟐𝒏 = 𝒃 et 𝒙𝒌 = 𝒂 + 𝒉𝒌 pour
k=0,1,2,……,2n.

La formule de Simpson généralisée s’écrit :

𝒃
𝒉
∫ 𝒇(𝒙) ≅ (𝒇(𝒙𝟎 ) + 𝟐 ∑ 𝒇(𝒙𝒊 ) + 𝟒 ∑ 𝒇(𝒙𝒊 ) + 𝒇(𝒙𝟐𝒏 ))
𝒂 𝟑 𝒊 𝒑𝒂𝒊𝒓 𝒊 𝒊𝒎𝒑𝒂𝒊𝒓

L’erreur d’interpolation de la formule de Simpson s’écrit :

𝒃−𝒂
𝑹(𝒇) = − 𝟏𝟖𝟎 𝒉𝟒 𝒇(𝟒) (𝒛) avec 𝒛 ∈ [𝒂, 𝒃]

20
Méthodes numériques

1 2
Exemple : Soit à calculer l’intégrale ∫0 𝑒 −𝑥 𝑑𝑥 avec une précision de 0.001 par la
méthode de Simpson. On doit premièrement trouver le nombre de division à faire
pour obtenir cette précision.
L’erreur d’intégration s’écrit :

(𝒃−𝒂)
𝑹(𝒇) = − 𝟏𝟖𝟎
𝒉𝟒 𝒇(𝟒) (𝒛)
sa valeur absolue doit être inférieure ou égale à la précision donnée (0.001), c.à.d. :

(𝒃 − 𝒂) 𝟒 (𝟒)
|𝑹(𝒇)| = |− 𝒉 𝒇 (𝒛)| ≤ 𝟎. 𝟎𝟎𝟏
𝟏𝟖𝟎

2 2
La fonction 𝑓(𝑥) = 𝑒 −𝑥 sa dérivée quatrième est 𝑓 (4) (𝑥) = (16𝑥 4 − 48𝑥 2 + 12)𝑒 −𝑥 ,
cette fonction est non monotone dans l’intervalle donné. On calcule son maximum
par le traceur Origin (Fig. 3.3.).

𝑀 = 𝑚𝑎𝑥|𝒇(𝟒) (𝒛)| = 𝟏𝟐 à 𝒙 = 𝟎.

𝟒 𝟏𝟖𝟎∗𝟎.𝟎𝟎𝟏
D’ou 𝒉 ≤ √ (𝟏−𝟎)∗𝟏𝟐 = 𝟎. 𝟑𝟓

𝟏 1
donc 𝟐𝒏 = = 𝟐. 𝟖𝟓 on prend 𝟒 divisions, le pas d’intégration ℎ = = 0.25.
𝟎.𝟑𝟓 4

1 2 0.25 −02 2 2 2 2
On trouve : ∫0 𝑒 −𝑥 𝑑𝑥 = 3
(𝑒 + 4(𝑒 −0.25 + 𝑒 −0.75 ) + 2𝑒 −0.5 + 𝑒 −1 ) = 0.7469

21
Méthodes numériques

3.3 Méthode de quadrature

Cette méthode permet l’élaboration de formules d’intégration numériques en

se basant sur le polynôme de Lagrange. On peut par exemple retrouver les formules
du trapèze et de Simpson par cette méthode ou bien de construire d’autres formules
d’intégration plus performantes. Dans cette méthode on remplace la fonction par un
polynôme de Lagrange puis on intègre le polynôme trouvé, on écrit donc :

𝒇(𝒙) ≅ 𝑷𝒏 (𝒙) = ∑𝒏𝒊=𝟎 𝒇𝒊 𝑳𝒊 (𝒙)

𝒃 𝒃 𝒃
On intègre : ∫𝒂 𝒇(𝒙)𝒅𝒙 ≅ ∫𝒂 𝑷𝒏 (𝒙)𝒅𝒙 = ∑𝒏𝒊=𝟎 𝒇𝒊 ∫𝒂 𝑳𝒊 (𝒙)𝒅𝒙

𝒃
Si on pose : 𝑨𝒊 = ∫𝒂 𝑳𝒊 (𝒙)𝒅𝒙 pour i=0,1,…,n

𝒃
On obtient la formule de quadrature : ∫𝒂 𝒇(𝒙)𝒅𝒙 ≅ ∑𝒏𝒊=𝟎 𝒇𝒊 𝑨𝒊

On doit maintenant choisir la forme de la fonction f(x), dans notre cas on prend :

𝒇(𝒙) = 𝒙𝒌 avec k=0,1,2,….,n

𝒃 𝒃𝒌+𝟏 −𝒂𝒌+𝟏
On remplace dans l’intégrale : ∫𝒂 𝒙𝒌 𝒅𝒙 ≅ ∑𝒏𝒊=𝟎 𝒇𝒊 𝑨𝒊 = 𝒌+𝟏
k=0,1,2,….,n

En variant i de 0 à n pour chaque valeur de k, on aura:

𝒃𝟏 −𝒂𝟏
pour k=0 : 𝒙𝟎𝟎 𝑨𝟎 + 𝒙𝟎𝟏 𝑨𝟏 + ⋯ … + 𝒙𝟎𝒏 𝑨𝟎 =
𝟏
𝒃𝟐 −𝒂𝟐
pour k=1 : 𝒙𝟏𝟎 𝑨𝟎 + 𝒙𝟏𝟏 𝑨𝟏 + ⋯ … + 𝒙𝟏𝒏 𝑨𝟎 = 𝟐
𝒃𝟑 −𝒂𝟑
pour k=2 : 𝒙𝟐𝟎 𝑨𝟎 + 𝒙𝟐𝟏 𝑨𝟏 + ⋯ … + 𝒙𝟐𝒏 𝑨𝟎 =
𝟑
……………………………………………………………………………………

𝒃𝒏+𝟏 −𝒂𝒏+𝟏
pour k=n : 𝒙𝒏𝟎 𝑨𝟎 + 𝒙𝒏𝟏 𝑨𝟏 + ⋯ … + 𝒙𝒏𝒏 𝑨𝟎 = 𝒏+𝟏

D’où on obtient pour toutes les valeurs de i et k le système suivant :

𝟏 𝟏 𝟏 … … 𝟏 𝑨𝟎 𝑰𝟎
𝒙𝟎 𝒙𝟏 𝒙𝟐 … . . 𝒙𝒏 𝑨𝟏 𝑰𝟏
𝒃𝒌+𝟏 −𝒂𝒌+𝟏
𝟐 𝟐 𝟐
𝒙𝟎 𝒙𝟏 𝒙𝟐 … . 𝒙𝒏 𝟐 𝑨 𝟐 = 𝑰𝟐 avec 𝑰𝒌 =
𝒌+𝟏
…………………. . .
[ 𝒙𝒏𝟎 𝒙𝒏𝟏 𝒙𝒏𝟐 … . 𝒙𝒏𝒏 ] [𝑨𝒏 ] [𝑰𝒏 ]

Le déterminant de la matrice du système trouvé est dit de ‘’Von-Dermonde’’ il est non

nul d’où la solution de ce système existe et unique (𝑨𝟎 , 𝑨𝟏 , … . . , 𝑨𝒏 ).

22
Méthodes numériques

1 2
Exemple : Soit à calculer l’intégrale ∫0 𝑒 −𝑥 𝑑𝑥 avec une formule qui a la forme
1 2
suivante : ∫0 𝑒 −𝑥 𝑑𝑥 = 𝑨𝟎 𝒇(𝟎) + 𝑨𝟏 𝒇(𝟎. 𝟐𝟓) + 𝑨𝟐 𝒇(𝟎. 𝟓) + 𝑨𝟑 𝒇(𝟎. 𝟕𝟓) + 𝑨𝟐 𝒇(𝟏. 𝟎𝟎). Dans ce
cas on cherche premièrement les constantes 𝐴𝑖 puis on calcule l’intégrale. Ces
constantes sont données par le système d’équations suivant :

1 1 1 1 1 𝑨𝟎 𝟏
0 0.25 0.5 0.75 1 𝑨𝟏 𝟎. 𝟓
0 0.0625 0.25 0.5626 1 𝑨𝟐 = 𝟎. 𝟑𝟑𝟑
0 0.015625 0.125 0.421875 1 𝑨𝟑 𝟎. 𝟐𝟓𝟎
[0 0.00390625 0.0625 0.31640625 1] [𝑨𝟒 ] [𝟎. 𝟐𝟎𝟎]

La solution du système est (0.0691, 0.3812, 0.1052, 0.3694, 0.0751) donc :

1
2 𝟐 𝟐 𝟐 𝟐
∫ 𝑒 −𝑥 𝑑𝑥 = 0.0691𝒆𝟎 + 0.3812𝒆−𝟎.𝟐𝟓 + 0.1052𝒆−𝟎.𝟓 + 0.3694𝒆−𝟎.𝟕𝟓 + 0.0751𝒆−𝟏 = 𝟎. 𝟕𝟒𝟕𝟐
0

23
Méthodes numériques

Chapitre VI : Résolution numérique des équations

différentielles ordinaires : Problème de Cauchy

Les équations différentielles sont utilisées dans la modélisation mathématique

de la quasi- totalité des phénomènes physiques. Une équation différentielle est une
relation entre une variable et sa ou ses dérivées de divers ordres. Dans ce chapitre,
on va considérer les équations différentielles ordinaires du premier ordre avec une
condition initiale imposée (problème de Cauchy).
Le problème de Cauchy est défini par la solution de l’équation différentielle
suivante :
𝒅𝒚(𝒕)
= 𝒚′ = 𝒇(𝒕, 𝒚(𝒕))
{ 𝒅𝒕 𝒚(𝒕𝟎 ) = 𝒚𝟎 est une condition initiale imposée.
𝒚(𝒕𝟎 ) = 𝒚𝟎
On remarque que la solution doit obligatoirement passer par le point (𝒕𝟎 , 𝒚𝟎 ).
La solution du problème de Cauchy existe et unique si la fonction 𝒇(𝒕, 𝒚(𝒕))
satisfait la condition de Lipschitz en y sur un rectangle R définit par 𝒂 ≤ 𝒕 ≤ 𝒃 et
𝒄 ≤ 𝒚 ≤ 𝒅. Cette condition nécessite que |𝒇(𝒕, 𝒚𝟏 ) − 𝒇(𝒕, 𝒚𝟐 )| ≤ 𝑳|𝒚𝟏 − 𝒚𝟐 |.
𝝏𝒇(𝒕,𝒚)
En pratique pour vérifier cette condition, on calcule 𝑴𝒂𝒙 | 𝝏𝒚
| ≤ 𝑳 sur R

4.1 Méthode d’Euler

Soit un intervalle [a,b] sur lequel on cherche la solution d’un problème de
Cauchy
𝒅𝒚(𝒕)
= 𝒚′ = 𝒇(𝒕, 𝒚(𝒕))
{ 𝒅𝒕
𝒚(𝒕𝟎 = 𝒂) = 𝒚𝟎
La fonction 𝒇(𝒕, 𝒚(𝒕)) vérifie la condition de Lipschitz en y sur le rectangle R.
La première étape consiste à diviser l’intervalle donné en n points équidistants
𝒃−𝒂
ce qui donne un pas d’intégration 𝒉 = 𝒏
, le point d’abscisse 𝒕𝒊 est donné par 𝒕𝒊 = 𝒂 +

𝒊𝒉 (pour i=1,2,…,n).

a=t0 t1 t2 ………….. tn-1 b=tn t

24
Méthodes numériques

On va donc résoudre le problème sur l’intervalle [𝒂 = 𝒕𝟎 , 𝒃 = 𝒕𝒏 ] avec 𝒚(𝒕𝟎 = 𝒂) =

𝒚𝟎 . Si les fonctions 𝑦(𝒕), 𝒚′ (𝒕) et 𝒚′′(𝒕) sont continues, on peut écrire le développement
en série de Taylor pour 𝒚(𝒕) au voisinage de 𝒕𝟎 .

(𝒕 − 𝒕𝟎 ) (𝒕 − 𝒕𝟎 )𝟐
𝒚(𝒕) = 𝒚(𝒕𝟎 ) + 𝒚′ (𝒕𝟎 ) + 𝒚′′ (𝒕𝟎 ) +⋯
𝟏! 𝟐!
Sachant que 𝒚′ (𝒕𝟎 ) = 𝒇(𝒕𝟎 , 𝒚(𝒕𝟎 )) et 𝒉 = 𝒕𝟏 − 𝒕𝟎 ,
𝒉𝟐
Ecrivons 𝒚(𝒕𝟏 ) = 𝒚(𝒕𝟎 ) + 𝒇(𝒕𝟎 , 𝒚(𝒕𝟎 ))𝒉 + 𝒚′′ (𝒕𝟎 ) 𝟐! + ⋯

Si h est suffisamment petit alors on peut négliger 𝒉𝟐 et donc les termes d’ordre deux
et plus, on obtient donc :

𝒚(𝒕𝟏 ) = 𝒚(𝒕𝟎 ) + 𝒉𝒇(𝒕𝟎 , 𝒚(𝒕𝟎 )) → 𝒚𝟏 = 𝒚𝟎 + 𝒉𝒇(𝒕𝟎 , 𝒚𝟎 )

C’est l’approximation d’Euler d’ordre 1, en répétant le procédé on génère une

séquence de points 𝒚𝟏 , 𝒚𝟐 , 𝒚𝟑 , … … , 𝒚𝒏−𝟏 , 𝒚𝒏 approximant 𝒚 = 𝒚(𝒕), en général :

𝒕 = 𝒂 + 𝒊𝒉, 𝒊 = 𝟎, 𝟏, 𝟐, … . . , 𝒏 − 𝟏
{ 𝒊
𝒚𝒊+𝟏 = 𝒚𝒊 + 𝒉𝒇(𝒕𝒊 , 𝒚𝒊 ), 𝒚𝟎 = 𝒚(𝒕𝟎 )

Exemple : Résoudre le problème de Cauchy suivant par la méthode d’Euler en

prenant un pas d’intégration h=0.25.

𝑦 ′ = 2 − 𝑡𝑦 2 𝑡 ∈ [0,1], 𝑦 ∈ [0,1]
{
𝑦(0) = 1

Vérifions la condition de Lipchitz sur le rectangle R défini par 𝑡 ∈ [0,1], 𝑦 ∈ [0,1].

𝝏𝒇(𝒕,𝒚) 𝝏(2−𝑡𝑦 2 )
Max | 𝝏𝒚
| = Max | 𝝏𝒚
| = |−𝟐𝒕𝒚| = 𝟐 < 𝑳 Condition vérifiée

On divise l’intervalle de t, [0,1] avec une pas h=0.25, soit 0, 0.25, 0.50, 0.75 et 1.

t0 =0 t1=0.25 t2 =0.50 t3 =0.75 t4 =1

y0 =y(0)=1 y1 y2 y3 y4

𝒕𝒊 = 𝒊𝒉, 𝒊 = 𝟎, 𝟏, 𝟐 et 𝟑
On écrit : {
𝒚𝒊+𝟏 = 𝒚𝒊 + 𝒉𝒇(𝒕𝒊 , 𝒚𝒊 ) = 𝒚𝒊 + 𝟎. 𝟐𝟓(𝟐 − 𝒕𝒊 𝒚𝟐𝒊 ), 𝒚𝟎 = 𝟏

25
Méthodes numériques

𝒊 = 𝟎, 𝒚𝟏 = 𝒚𝟎 + 𝒉𝒇(𝒕𝟎 , 𝒚𝟎 ) = 𝟏 + 𝟎. 𝟐𝟓(𝟐 − 𝟎 ∗ 𝟏𝟐 ) = 𝟏. 𝟓
𝒊 = 𝟏, 𝒚𝟐 = 𝒚𝟏 + 𝒉𝒇(𝒕𝟏 , 𝒚𝟏 ) = 𝟏. 𝟓 + 𝟎. 𝟐𝟓(𝟐 − 𝟎. 𝟐𝟓 ∗ 𝟏. 𝟓𝟐 ) = 𝟏. 𝟖𝟓𝟗𝟒
𝒊 = 𝟐, 𝒚𝟑 = 𝒚𝟐 + 𝒉𝒇(𝒕𝟐 , 𝒚𝟐 ) = 𝟏. 𝟖𝟓𝟗 + 𝟎. 𝟐𝟓(𝟐 − 𝟎. 𝟓 ∗ 𝟏. 𝟖𝟓𝟗𝟐 ) = 𝟏. 𝟗𝟐𝟕𝟐
𝒊 = 𝟑, 𝒚𝟒 = 𝒚𝟑 + 𝒉𝒇(𝒕𝟑 , 𝒚𝟑 ) = 𝟏. 𝟗𝟐𝟕 + 𝟎. 𝟐𝟓(𝟐 − 𝟎. 𝟕𝟓 ∗ 𝟏. 𝟗𝟐𝟕𝟐 ) = 𝟏. 𝟕𝟑𝟎𝟖

La solution exacte est :

𝟑 𝟑
𝟐 √𝟐𝒄𝟏 𝒕𝑨𝒊(𝟐𝟐⁄𝟑 𝒕)+𝟐 √𝟐𝒕𝑩𝒊(𝟐𝟐⁄𝟑 𝒕)
𝒚(𝒕) = Cas de l’exemple
𝟐𝒕(𝒄𝟏 𝑨′𝒊 (𝟐𝟐⁄𝟑 𝒕)+𝑩′𝒊 (𝟐𝟐⁄𝟑 𝒕))

Avec Ai et Bi sont les fonctions

d’Airy et d’Airy de seconde espèce.

Fig. 4.1. Solutions y(t) pour différentes

conditions initiales y(0).

26
Méthodes numériques

4.2 Méthode d’Euler améliorée(Huns)

D’après le développement en série de Taylor, on a :

(𝒕 − 𝒕𝟎 ) (𝒕 − 𝒕𝟎 )𝟐 (𝒕 − 𝒕𝟎 )𝟑
𝒚(𝒕) = 𝒚(𝒕𝟎 ) + 𝒚′ (𝒕𝟎 ) ′′ (𝒕 )
+𝒚 𝟎 +𝒚′′′ (𝒕 )
𝟎 +⋯
𝟏! 𝟐! 𝟑!
Si on prend les trois premiers termes et on néglige les autres d’ordre supérieur, on
obtient :

(𝒕 − 𝒕𝟎 ) (𝒕 − 𝒕𝟎 )𝟐
𝒚(𝒕) ≅ 𝒚(𝒕𝟎 ) + 𝒚′ (𝒕𝟎 ) + 𝒚′′ (𝒕𝟎 ) +⋯
𝟏! 𝟐!
En remplaçant 𝒕 par 𝒕𝟏 cela donne :

(𝒕𝟏 − 𝒕𝟎 ) (𝒕𝟏 − 𝒕𝟎 )𝟐 𝒉𝟐
𝒚(𝒕𝟏 ) = 𝒚(𝒕𝟎 ) + 𝒚′ (𝒕𝟎 ) + 𝒚′′ (𝒕𝟎 ) = 𝒚(𝒕𝟎 ) + 𝒚′ (𝒕𝟎 )𝒉 + 𝒚′′ (𝒕𝟎 )
𝟏! 𝟐! 𝟐!
𝒚′(𝒕𝟏 )−𝒚′(𝒕𝟎 ) 𝒉
Or 𝒚′′ (𝒕𝟎 ) = on aura 𝒚(𝒕𝟏 ) = 𝒚(𝒕𝟎 ) + [𝒇(𝒕𝟏 , 𝒚𝟏 ) + 𝒇(𝒕𝟎 , 𝒚𝟎 )]
𝒉 𝟐

En général la formule d’Euler améliorée s’écrit :

𝒉
𝒚𝒊+𝟏 = 𝒚𝒊 + [𝒇(𝒕𝒊+𝟏 , 𝒚𝑬𝒊+𝟏 ) + 𝒇(𝒕𝒊 , 𝒚𝒊 )]
𝟐
On remarque que cette formule donne 𝒚𝒊+𝟏 en fonction de 𝒚𝑬𝒊+𝟏 qui doit être calculée
par la méthode d’Euler.

Exemple : Résoudre le problème de Cauchy précédant par la méthode d’Euler

améliorée en prenant un pas d’intégration h=0.25.
On a :
𝒕𝒊 = 𝒊𝒉, 𝒊 = 𝟎, 𝟏, 𝟐 et 𝟑
{ 𝒉 𝒉
𝒚𝒊+𝟏 = 𝒚𝒊 + [𝒇(𝒕𝒊+𝟏 , 𝒚𝑬𝒊+𝟏 ) + 𝒇(𝒕𝒊 , 𝒚𝒊 )] = 𝒚𝒊 + [(𝟐 − 𝒕𝒊 𝒚𝟐𝒊 ) + (𝟐 − 𝒕𝒊+𝟏 𝒚𝟐𝒊+𝟏𝑬 )], 𝒚𝟎 = 𝟏
𝟐 𝟐
Avec : 𝒚𝑬𝒊+𝟏 = 𝒚𝒊 + 𝒉(𝟐 − 𝒕𝒊 𝒚𝟐𝒊 ) la solution obtenue par la méthode d’Euler. On
𝒉 𝟐
remplace 𝒚𝑬𝒊+𝟏 , on aura : 𝒚𝒊+𝟏 = 𝒚𝒊 + [(𝟐 − 𝒕𝒊 𝒚𝟐𝒊 ) + (𝟐 − 𝒕𝒊+𝟏 (𝒚𝒊 + 𝒉(𝟐 − 𝒕𝒊 𝒚𝟐𝒊 )) )],
𝟐

𝒚𝟏 = 𝟏. 𝟒𝟐𝟗𝟕, 𝒚𝟐 = 𝟏. 𝟔𝟔𝟐𝟗, 𝒚𝟑 = 𝟏. 𝟔𝟖𝟎𝟓, 𝒚𝟒 = 𝟏. 𝟓𝟕𝟓𝟎

27
 DEUXIEME PARTIE

Introduction à la programmation en langage Python

1. Présentation

Le Langage Python

Python est un langage de programmation (au même titre que le C, C++, fortran, java . . .),
développé en 1989. Ses principales caractéristiques sont les suivantes :
— «open-source» : son utilisation est gratuite et les fichiers sources sont disponibles et modi-
fiables ;
— simple et très lisible ;
— doté d’une bibliothèque de base très fournie ;
— importante quantité de bibliothèques disponibles : pour le calcul scientifique, les statis-
tiques, les bases de données, la visualisation . . . ;
— grande portabilité : indépendant vis à vis du système d’exploitation (linux, windows, Ma-
cOS) ;
— orienté objet ;
— typage dynamique : le typage (association à une variable de son type et allocation zone
mémoire en conséquence) est fait automatiquement lors de l’exécution du programme, ce
qui permet une grande flexibilité et rapidité de programmation, mais qui se paye par une
surconsommation de mémoire et une perte de performance ;
— présente un support pour l’intégration d’autres langages.

Comment faire fonctionner le code source ?

Il existe deux techniques principales pour traduire un code source en langage machine :
— la compilation : une application tierce, appelée compilateur, transforme les lignes de code en
un fichier exécutable en langage machine. A chaque fois que l’on apporte une modification au
programme, il faut recompiler avant de voir le résultat.
— l’interprétation : un interpréteur s’occupe de traduire ligne par ligne le programme en langage
machine. Ce type de langage offre une plus grande commodité pour le développement, mais
les exécutions sont souvent plus lentes.
Dans le cas de Python, on peut admettre pour commencer qu’il s’agit d’un langage interprété,
qui fait appel des modules compilés. Pour les opérations algorithmiques coûteuses, le langage
Python peut s’interfacer à des bibliothèques écrites en langage de bas niveau comme le langage C.

Les différentes versions

Il existe deux versions de Python : 2.7 et 3.3. La version 3.3 n’est pas une simple amélioration
de la version 2.2. Attention, toutes les librairies Python n’ont pas effectué la migration de 2.7 à 3.3.

L’interpréteur

Dans un terminal, taper python (interpréteur classique) ou ipython (interpréteur plus évolué)
pour accéder à un interpréteur Python. Vous pouvez maintenant taper dans ce terminal des ins-
tructions Python qui seront exécutées.

28
 Utilisations de Python et librairies

— web : Django, Zope, Plone,. . .

— bases de données : MySQL, Oracle,. . .
— réseaux : TwistedMatrix, PyRO, VTK,. . .
— représentation graphique : matplotlib, VTK,. . .
— calcul scientifique : numpy, scipy, . . .

Le mode programmation

Il s’agit d’écrire dans un fichier une succession d’instructions qui ne seront éffectuées que lorsque
vous lancerez l’exécution du programme. Cela permet tout d’abord de sauvegarder les commandes
qui pourront être utilisées ultérieurement, et d’autre part d’organiser un programme, sous-forme
de fichier principal, modules, fonctions . . ..
Le fichier à exécuter devra avoir l’extension .py, et devra contenir en première ligne le chemin
pour accéder au compilateur Python, ainsi que l’encodage :
#!/usr/bin/env python
# -*- coding: utf-8 -*-

Il pourra être exécuté en lançant dans un terminal la commande python [Link].

L’environnement de programmation Spyder (sous Anaconda)

L’environnement Spyder permet de réaliser des programmes informatiques écrits avec le lan-
gage Python. Il est disponible avec la distribution Anaconda, qui présente de nombreux avan-
tages, notamment celui d’être simple à installer. Son téléchargement se fait à l’adresse suivante :
https ://[Link]/cshop/anaconda/
Une fois la distribution Anaconda téléchargée et installée, on peut commencer à lancer Spyder en
tapant Spyder dans un terminal.

Au lancement de Spyder, apparaît une fenêtre partagée en deux zones. La zone en-bas à droite,
appelée console (shell en anglais), est celle où l’on peut travailler de façon interactive avec l’in-
terprteur Python. La zone à gauche est un éditeur de texte, spécialement conçu pour écrire des
programmes dans le langage Python.

2. Types et opérations de base

Le langage Python est orienté objet, c’est-à-dire qu’il permet de créer des objets, en définissant
des attributs et des fonctions qui leur sont propres. Cependant, certains objets sont pré-définis
dans le langage. Nous allons voir à présent les plus importants.

29
(a) les nombres et les booléens
— entiers (32 bits)
type : int
— réels (64 bits)
type : float
exemples de valeurs : 4. 5.1 1.23e-6
— complexes
type : complex
exemples de valeurs : 3+4j 3+4J
— booléens type : bool
exemples de valeurs : True False

(b) opérations de base

— affectation

>>> i = 3 # i vaut 3
>>> a, k=True, 3.14159
>>> k=r=2.15
>>> x=complex(3,4)
— affichage

>>> i
3
>>> print(i)
3
— Opérateurs addition, soustraction, multiplication, division

+, -, *, \, \%,
— Opérateurs puissance, valeur absolue

**, pow, abs, \dots

— Opérateurs de comparaison

==, is, !=, is not, >, >=, <, <=

— Opérateurs logiques

or, and, not

— Conversion

>>> int(3.1415)
3
>>> float(3)
3.

(c) les chaînes de caractères

Une chaîne de caractères (string en anglais) est un objet de la classe (ou de type) str. Une
chaîne de caractères peut être définie de plusieurs façons :
>>> "je suis une chaine"
’je suis une chaine’
>>> ’je suis une chaine’
’je suis une chaine’

30
 >>> ’pour prendre l\’apostrophe’
’pour prendre l’ apostrophe’
>>> "pour prendre l’apostrophe"
"pour prendre l’apostrophe"
>>> " " "ecrire
sur
plusieurs
lignes" " "
’ecrire\nsur\nplusieurs\nlignes’

— concaténation
On peut mettre plusieurs chaînes de caractères bout à bout avec l’opérateur binaire de
concaténation, noté +.
>>> s = ’i vaut’
>>> i = 1
>>> print( s+i )
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: cannot concatenate ’str’ and ’int’ objects
>>> print( s + " %d %s "%(i, "m."))
i vaut 1 m.
>>> print( s + ’ ’ + str(i))
i vaut 1
>>> print(’*-’ * 5)
*-*-*-*-*-

— accès aux caractères

Les caractères qui composent une chaîne sont numérotés à partir de zéro. On peut y accéder
individuellement en faisant suivre le nom de la chaîne d’un entier encadré par une paire
de crochets :
>>> "bonjour"[3]; "bonjour"[-1]
’j’
’r’
>>> "bonjour"[2:]; "bonjour"[:3]; "bonjour"[3:5]
’njour’
’bon’
’jo’
>>> "bonjour"[-1::-1];
’ruojnob’

— méthodes propres
— len(s) : renvoie la taille d’une chaîne,
— [Link] : recherche une sous-chaîne dans la chaîne,
— [Link] : enlève les espaces de fin,
— [Link] : remplace une chaîne par une autre,
— ...

(d) les listes

Une liste consiste en une succession ordonnée d’objets, qui ne doivent pas nécessairement être
du même type. Les termes d’une liste sont numérotés à partir de 0 (comme pour les chaînes
de caractères).
— initialisation

31
 >>> []; list();
[]
[]
>>> [1,2,3,4,5]; [’point’,’triangle’,’quad’];
[1, 2, 3, 4, 5]
[’point’, ’triangle’, ’quad’]
>>> [1,4,’mesh’,4,’triangle’,[’point’,6]];
[1, 4, ’mesh’, 4, ’triangle’, [’point’, 6]]
>>> range(10)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> range(2,10,2)
[2, 4, 6, 8]

— modification
Contrairement aux chaînes de caractères, on peut modifier les éléments d’une liste :
>>> l=[1,2,3,4,5]
>>> l[2:]=[2,2,2]
>>> l
[1, 2, 2, 2, 2]

— concaténation
>>> [0]*7
[0, 0, 0, 0, 0, 0, 0]
>>> L1, L2 = [1,2,3], [4,5]
>>> L1
[1, 2, 3]
>>> L2
[4, 5]
>>> L1+L2
[1, 2, 3, 4, 5]

— méthodes propres
— len(L) : renvoie la taille de la liste L,
— [Link] : trie la liste L,
— [Link] : ajoute un élément à la fin de la liste L,
— [Link] : inverse la liste L,
— [Link] : recherche un élément dans la liste L,
— [Link] : retire un élément de la liste L,
— [Link] : retire le dernier élément de la liste L,
— ...

(e) copie d’un objet

>>> L = [’Dans’,’python’,’tout’,’est’,’objet’]
>>> T = L
>>> T[4] = ’bon’
>>> T
[’Dans’, ’python’, ’tout’, ’est’, ’bon’]
>>> L
[’Dans’, ’python’, ’tout’, ’est’, ’bon’]
>>> L=T[:]
>>> L[4]=’objet’
>>> T;L
[’Dans’, ’python’, ’tout’, ’est’, ’bon’]

32
 [’Dans’, ’python’, ’tout’, ’est’, ’objet’]

(f) quelques remarques importantes

— en Python, tout est objet,

— une chaîne de caractères est immuable, tandis qu’une liste est muable,
— type permet de connaître le type d’un objet,
— id permet de connaître l’adresse d’un objet,
— eval permet d’évaluer une chaîne de caractères.

3. Commentaires
# ceci est un commentaire

4. Noms de variables

Python fait la distinction entre minuscules et majuscules.

Conseil : donner des noms significatifs aux variables et aux fonctions.

5. Les structures de contrôle

(a) L’indentation

Les fonctions Python n’ont pas de begin ou end explicites, ni d’accolades qui pourraient
marquer là où commence et où se termine le code de la fonction. Le seul délimiteur est les
deux points (« : ») et l’indentation du code lui-même. Les blocs de code (fonctions, instructions
if, boucles for ou while etc) sont définis par leur indentation. L’indentation démarre le bloc et
la désindentation le termine. Il n’y a pas d’accolades, de crochets ou de mots clés spécifiques.
Cela signifie que les espaces blancs sont significatifs et qu’ils doivent être cohérents. Voici un
exemple :
a = -150
if a < 0:
print(’a est negatif’)

ligne d’en-tête
première instruction du bloc
. . .
dernière instruction du bloc

Fonctionnement par blocs :

Bloc 1
...
Ligne d’en-tête :
Bloc 2
...
Ligne d’en-tête :
Bloc 3
...
Ligne d’en-tête :
Bloc 2 (suite)
...
Bloc 1 (suite)
...

33
(b) Le test de conditions

Le test de condition se fait sous la forme générale suivante :

if < test1 > :
< blocs d’instructions 1>
elif < test2 > :
< blocs d’instructions 2 >
else :
< blocs d’instructions 3 >

Un exemple :

a = 10.
if a > 0:
print(’a est strictement positif’)
if a >= 10:
print (’ a est un nombre ’)
else:
print (’ a est un chiffre ’)
a += 1
elif a is not 0:
print(’ a est strictement negatif ’)
else:
print(’ a est nul ’)

Un autre exemple :

L = [1, 3, 6, 8]
if 9 in L:
print ’9 est dans la liste L’
else:
[Link](9)

(c) La boucle conditionnelle

La forme générale d’une boucle conditionnelle est

while < test1 > :
< blocs d’instructions 1 >
if < test2 > : break
if < test3 > : continue
else :
< blocs d’instructions 2 >
où l’on a utilisé les méthodes suivantes :

break : sort de la boucle sans passer par else,

continue : remonte au début de la boucle,
pass : ne fait rien.
La structure finale else est optionnelle. Elle signifie que l’instruction est lancée si et seulement
si la boucle se termine normalement.

Quelques exemples :

— Boucle infinie :

34
 while 1:
pass
— y est-il premier ?
x = y/2
while x > 1 :
if y%x ==0
print (str(y)+’ est facteur de ’+str(x))
break
x = x - 1
else :
print( str(y)+ ’ est premier’)

(d) La boucle inconditionnelle

La forme générale d’une boucle inconditionnelle est

for < cible > in < objet > :
< blocs d’instructions 1 >
if < test1 > : break
if < test2 > : continue
else :
< blocs d’instructions 2 >

Quelques exemples :

sum = 0
for i in [1, 2, 3, 4] :
sum += 1

prod = 1
for p in range(1, 10) :
prod *= p

s = ’bonjour’
for c in s :
print c,

L = [ x + 10 for x in range(10) ]

6. Les fonctions

La structure générale de définition d’une fonction est la suivante

def < nom fonction > (arg1, arg2,... argN):
...
bloc d’instructions
...
return < valeur( s ) >

Voici quelques exemples :

def table7():
n=1
while n < 11:
print n*7
n+=1

35
 Un fonction qui n’a pas de return renvoie par défaut None.

def table(base):
n=1
while n < 11:
print n*base
n+=1

def table(base, debut=0, fin=11):

print (’Fragment de la table de multiplication par ’+str(base)+’ : ’)
n=debut
l=[]
while n < fin:
print n*base
[Link](n*base)
n+=1
return l

On peut également déclarer une fonction sans connaître ses paramètres :

>>> def f (*args, **kwargs):
... print(args)
... print(kwargs)

>>> f(1,3,’b’,j=1)
(1, 3, ’b’)
{’j’: 1}

Il existe une autre façon de déclarer une fonction de plusieurs paramètres :

lambda argument 1, ... , argument N : expression utilisant les arguments

Exemple :
>>> f = lambda x, i : x**i
>>> f(2,4)
16

7. Les modules

Un module est un fichier comprenant un ensemble de définitions et d’instructions compréhensibles

par Python. Il permet d’étendre les fonctionnalités du langage. Voici tout d’abord un exemple de
module permettant l’utilisation des nombres de Fibonacci.

Fichier [Link]
# Module nombres de Fibonacci
def print_fib(n) :
" " "
ecrit la serie de Fibonacci jusqu’a n
" " "
a, b = 0, 1
while b < n:
print(b),

36
 a, b = b, a + b
print

def print_fib(n) :
" " "
retourne la serie de Fibonacci jusqu’a n
" " "
result, a, b = [], 0, 1
while b < n:
[Link](b),
a, b = b, a + b
return result

Utilisation du module [Link]

>>> import fibo
>>> fibo.print_fib(1000)
1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987
>>> fibo.list_fib(100)
[1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89]

L’importation de modules

Il existe plusieurs façons d’importer un module :

— import fibo
— import fibo as f
— from fibo import print_ fib, list_ fib
— from fibo import * (importe tous les noms sauf les variables et les fonctions privées)
Le module math

Ce module fournit un ensemble de fonctions mathématiques pour les réels :

pi
sqrt
cos, sin, tan, acos, ...
...
Tapez import math pour importer le module, puis help(math) pour avoir toutes les informations
sur le module : chemin, fonctions, constantes.

8. Le module NumPy

NumPy est un outil performant pour la manipulation de tableaux à plusieurs dimensions. Il

ajoute en effet le type array, qui est similaire à une liste, mais dont tous les éléments sont du
même type : des entiers, des flottants ou des booléens.

Le module NumPy possède des fonctions basiques en algèbre linéaire, ainsi que pour les transfor-
mées de Fourier.

(a) Création d’un tableau dont on connaît la taille

>>> import numpy as np

>>> a = [Link](4)
>>> a

37
 array([ 0., 0., 0., 0.])
>>> nx, ny = 2, 2
>>> a=[Link]((nx,ny))
>>> a
array( [ [ 0., 0. ] ,
[ 0., 0.] ] )

Un tableau peut être multidimensionnel, comme ici, de dimension 3 :

>>> a=[Link]((nx,ny,3))
>>> a
array( [ [ [ 0., 0., 0. ],
[ 0., 0., 0. ] ],

[ [ 0., 0., 0. ],
[ 0., 0., 0. ] ] ] )
Mais nous nous limiterons dans ce cours aux tableaux uni et bi-dimensionnels.

Il existe également les fonctions [Link], [Link], [Link], [Link], ...

Par exemple [Link](3) est la matrice identité d’ordre 3, [Link] est un tableau vide,
[Link](5) est le vecteur [1 1 1 1 1] et [Link](3,2) est la matrice à 3 lignes et 2 colonnes
contenant des 1 sur la diagonale et des zéro partout ailleurs.

Par défaut les éléments d’un tableau sont des float (un réel en double précision) ; mais on
peut donner un deuxième argument qui précise le type (int, complex, bool, ...). Exemple :
>>> [Link](2, dtype=int)

(b) Création d’un tableau avec une séquence de nombre

La fonction linspace(premier, dernier, n) renvoie un tableau unidimensionnel commen-

çant par premier, se terminant par dernier avec n éléments régulièrement espacés. Une
variante est la fonction arange :
>>> a = [Link](-4, 4, 9)
>>> a
array([-4., -3., -2., -1., 0., 1., 2., 3., 4.])
>>> a = [Link](-4, 4, 1)
>>> a
array([-4, -3, -2, -1, 0, 1, 2, 3])

On peut convertir une ou plusieurs listes en un tableau via la commande array :

>>> a = [Link]([1, 2, 3])
>>> a
array([1, 2, 3])
>>> b = [Link](range(10))
>>> b
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> L1, L2 = [1, 2, 3], [4, 5, 6]
>>> a = [Link]([L1, L2])
>>> a
array([[1, 2, 3],
[4, 5, 6]])

(c) Création d’un tableau à partir d’une fonction

38
 >>> def f(x, y):
... return x**2 + [Link](y)
...
>>> a = [Link](f, (2, 3))
>>> a
array([[ 0. , 0.84147098, 0.90929743],
[ 1. , 1.84147098, 1.90929743]])
La fonction f a ici été appliquée aux valeurs x = [0, 1] et y = [0, 1, 2], le résultat est

[ [f (0, 0), f (0, 1), f (0, 2)], [f (1, 0), f (1, 1), f (0, 2)] ].

(d) Manipulations d’un tableau

i. Caractéristiques d’un tableau

— [Link] : retourne les dimensions du tableau
— [Link] : retourne le type des éléments du tableau
— [Link] : retourne le nombre total d’éléments du tableau
— [Link] : retourne la dimension du tableau (1 pour un vecteur, 2 pour une matrice)
ii. Indexation d’un tableau

Comme pour les listes et les chaines de caractères, l’indexation d’un vecteur (ou tableau
de dimension 1) commence à 0. Pour les matrices (ou tableaux de dimension 2), le premier
index se réfère à la ligne, le deuxième à la colonne. Quelques exemples :

>>> L1, L2 = [1, -2, 3], [-4, 5, 6]

>>> a = [Link]([L1, L2])
>>> a[1, 2] # extrait l’élément en 2ème ligne et 3ème colonne
6

>>> a[:, 1] # extrait toute la deuxième colonne

array([-2, 5])

>>> a[1,0:2:2] # extrait les éléments d’indice [dbut=0, pas=2, fin=2]

# de la première ligne, le dernier élément n’est pas inclus
array([-4])

>>> a[:, -1:0:-1] # extrait les éléments d’indice [dbut=-1, pas=-1, fin=0]
# de la première colonne, le dernier élément n’est pas inclus
array([[3, -2],
[6, 5]])

>>> a[a < 0] # extrait les éléments négatifs

array([-2, -4])

iii. Copie d’un tableau

Un tableau est un objet. Si l’on affecte un tableau A à un autre tableau B, A et B font ré-
férence au même objet. En modifiant l’un on modifie donc automatiquement l’autre. Pour
éviter cela on peut faire une copie du tableau A dans le tableau B moyennant la fonction
copy(). Ainsi, A et B seront deux tableaux identiques, mais ils ne feront pas référence au
même objet.

39
 >>> a = [Link](1, 5, 5)
>>> b = a
>>> c = [Link]()
>>> b[1] = 9
>>> a; b; c;
array([ 1., 9., 3., 4., 5.])
array([ 1., 9., 3., 4., 5.])
array([ 1., 2., 3., 4., 5.])

Une façon de faire une copie de tableau sans utiliser la méthode copy() est la suivante :
>>> d = [Link]([Link], [Link])
>>> d[:] = b
>>> d
array([ 1., 2., 3., 4., 5.])

iv. Redimensionnement d’un tableau

Voici comment modifier les dimensions d’un tableau :

>>> a = [Link](1, 10, 10)

>>> a
array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])
>>> [Link] = (2, 5)
>>> a
array([[ 1., 2., 3., 4., 5.],
[ 6., 7., 8., 9., 10.]])
>>> [Link] = ([Link],)
>>> a
array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])
>>> [Link](2, 5)
array([[ 1., 2., 3., 4., 5.],
[ 6., 7., 8., 9., 10.]])
>>> a
array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])
Attention : la fonction reshape ne modifie pas l’objet, elle crée une nouvelle vue.
La fonction flatten renvoie une vue d’un tableau bi-dimensionnel en tableau uni-dimensionnel.
>>> a = [Link](1, 10, 10)
>>> [Link] = (2, 5)
>>> a
array([[ 1., 2., 3., 4., 5.],
[ 6., 7., 8., 9., 10.]])
>>> [Link]
array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])

v. Boucles sur les tableaux

>>> a = [Link]((2, 3))
>>> for i in range([Link][0]):
... for j in range([Link][1]):
... a[i, j] = (i + 1)*(j + 1)
>>> print a
[[ 1. 2. 3.]
[ 2. 4. 6.]]

39
 >>> for e in a:
... print e
[ 1. 2. 3.]
[ 2. 4. 6.]

(e) Calculs avec des tableaux

Pour bien des calculs il est possible d’éviter d’utiliser des boucles (qui sont très consom-
matrices de temps de calcul). On peut faire directement les calculs sur des tableaux. Ceux-ci
sont faits via des fonctions C, un langage de programmation bas niveau, et sont donc plus
rapides. C’est ce qu’on appelle «vectoriser» un programme.

Voici un exemple, dans lequel on veut calculer b = 3a − 1 :

>>> a = [Link](0, 1, 1E+06)

>>> %timeit b = 3*a -1
100 loops, best of 3: 10.7 ms per loop
>>> b = [Link](1E+06)
>>> %timeit for i in xrange([Link]): b[i] = 3*a[i] - 1
10 loops, best of 3: 1.31 s per loop

Voici un autre exemple de calcul fait directement sur un tableau :

>>> def f1(x):
... return [Link](-x*x)*[Link](1+x*[Link](x))
...
>>> x = [Link](0, 1, 1e6)
>>> a = f1(x)

(f) Produit matriciel

Le calcul matriciel se fait avec la fonction dot().

Attention, si le deuxième argument est un vecteur ligne, il sera transformé si besoin en vecteur
colonne, par transposition. Ceci est dû au fait qu’il est plus simple de définir un vecteur ligne,
par exemple v=[Link]([0,1,2]), qu’un vecteur colonne, ici v=[Link]([[0],[1],[2]]).
Exemple :
>>> A=[Link]([[1,1,1],[0,1,1],[0,0,1]]
>>> v=[Link]([0,1,2])
>>> [Link](v,A)
array([0, 1, 3])
>>> [Link](A,v)
array([3, 3, 2])
On peut aussi transposer préalablement un vecteur ligne :
>>> v=[Link]([[0,1,2]])
>>> v=[Link](v)
>>> v
Attention, l’opération A**2 correspond à une élévation au carré terme à terme. Pour lever une
matrice au carré il faut taper [Link](A,A).

Pour lever une matrice carrée à une puissance n il faut faire une boucle :
>>> B=[Link]()
>>> for i in range(1,n):
... B = [Link](A,B)

40
 A la sortie de cette boucle B contiendra An .

(g) Produit scalaire

Il faut utiliser la fonction vdot().

(h) Le sous-module linalg

Ce module propose des méthodes numériques pour inverser une matrice, calculer son déter-
minant, résoudre un système linéaire . . .

9. Importer et exporter des données avec Python

(a) Méthode manuelle

Pour ouvrir et fermer un fichier nommé [Link] on crée une variable fichier de
type file à travers laquelle on pourra accéder au fichier.
fichier = open (’[Link]’,’r’)
Le paramètre «r» indique qu’on accède au fichier en mode lecture. Pour accéder au fichier en
mode écriture il faut utiliser le paramètre «w». A la fin des opérations sur le fichier, on le
ferme en appelant la fonction close :
[Link]()
Voici trois façons différentes de parcourir les lignes d’un fichier ouvert en mode écriture :
for ligne in fichier :
La variable ligne est une chaîne de caractère qui prendra les valeurs successives des lignes du
fichier de lecture (sous la forme d’une chaîne de caractères).
ligne = [Link]()
A chaque fois que la fonction readline de fichier est appelée, la ligne suivante du fichier est
lue et constitue la valeur de retour de la fonction.
lignes = [Link] ()
Cette fonction lit toutes les lignes d’un coup et la valeur de retour est une liste de chaînes de
caractères. Attention, si le fichier est gros cette méthode est à éviter car elle bloque instan-
tannément une grande place mémoire. Il est dans ce cas préférable d’utiliser une méthode de
lecture caractère par caractère (grâce à la fonction read()) ou ligne par ligne.

Enfin, la fonction [Link](chaine) permet d’écrire la chaîne de caractères en argu-

ment dans un fichier ouvert en mode écriture.

(b) Méthode automatisée

La fonction [Link] lit les lignes du fichier en argument et renvoie une matrice dont
les lignes correspondent aux lignes du fichier et les colonnes aux différentes valeurs délimitées
par un espace.

Parmi les options les plus utiles, skiprows est le nombre de lignes à ne pas lire dans l’en–tête
(par défaut 0) et delimiter remplace le délimiteur (espace) par la chaîne que l’on souhaite.

Exemple : le fichier [Link] contient les lignes :

41
 chat, chien
3,4
2,5
1,8
On veut supprimer la première ligne qui contient le nom des colonnes mais pas de valeurs
numériques, et on déclare que le délimiteur des valeurs est la virgule :
t = [Link] ("[Link]", skiprows=1, delimiter=’,’)
t
renvoie
array([ [ 3., 4.],
[ 2., 5.],
[ 1., 8.]])

10. Représentation discrète d’une fonction avec Matplotlib

Pour tracer la courbe représentative d’une fonction f avec Matplotlib, il faut tout d’abord la
“discrétiser” ; c’est-à-dire définir une liste de points (xi , f (xi )) qui va permettre d’approcher le
graphe de la fonction par une ligne brisée. Bien sûr, plus on augmente le nombre de points, plus
la ligne brisée est “proche” du graphe (en un certain sens).

Plus précisément, pour représenter le graphe d’une fonction réelle f définie sur un intervalle
[a, b], on commence par construire un vecteur X, discrétisant l’intervalle [a, b], en considérant un
ensemble de points équidistants dans [a, b]. Pour ce faire, on se donne un naturel N grand et
on considère X=[Link](a,b,N) (ou, ce qui revient au même, on pose h = Nb−a −1 et on
considère
X=[Link](a,b+h,h)).

On construit ensuite le vecteur Y dont les composantes sont les images des Xi par f et on
trace la ligne brise reliant tous les points (Xi , Yi ).

(a) La fonction plot du module [Link]

import [Link] as plt

import numpy as np

x = [Link](0., 2*[Link], 100)

[Link](x, [Link](x/2)*[Link](5*x), ’-ro’)

[Link](’Fonction $f(x)=e^{x/2} cos(5x)$’)

[Link](’$x$’)
[Link](1, 15, ’Courbe’, fontsize=22)
[Link](’1D_exemple.pdf’)
[Link]()

Cette fonction permet de tracer des lignes brisées. Si X et Y sont deux vecteurs-lignes (ou
vecteurs-colonnes) de même taille n, alors [Link](X,Y) permet de tracer la ligne brisée qui
relie les points de coordonnées (X(i),Y(i)) pour i = 1, . . . , n.

Pour connaître toutes les options, le mieux est de se référer à la documentation de Matplotlib.

42
 25
Fonction f(x) = ex/2 cos(5x)
20
15 Courbe
10
5
0
5
10
15
20
0 1 2 3 4 5 6 7
x

Exemple : tracé d’un carré

Coordonnées des sommets : X = (−1; 1; 1; −1; −1) et Y = (−1; −1; 1; 1; −1).

X=[Link]([-1, 1, 1,-1,-1])
Y=[Link]([-1, -1, 1,1,-1])
[Link](X,Y,’r’,lw=3)
[Link](’equal’)
[Link]([-2,2,-2,2])

[Link](’Carre de cote 2’)

[Link](’carre_exemple.pdf’)
[Link]()

Carre de cote 2
1.0

0.5

0.0

0.5

1.0

2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0

(b) Tracé de plusieurs lignes brisées

43
 On peut superposer plusieurs courbes dans le même graphique.
import [Link] as plt
import numpy as np
X = [Link](0, 2*[Link], 256)
Ycos = [Link](X)
Ysin = [Link](X)
[Link](X,Ycos,’b’)
[Link](X,Ysin,’r’)
[Link]()

1.0
cos(x)
sin(x)
0.5

0.0

0.5

1.0
0 1 2 3 4 5 6 7
x

(c) Création de plusieurs graphiques dans une même fenêtre

La commande subplot permet de découper la fenêtre graphique en plusieurs sous-fenêtres.

Voici un exemple d’utilisation de cette commande.

x = [Link](-4., 4., 25)

[Link](2, 2, 1)
[Link](x, -x**4 + x**3 + x**2 + 1, ’o-’)
[Link]("$f$")
[Link](2, 2, 2)
[Link](x, -4*x**3 + 3*x**2 + 2*x, ’-’)
[Link]("$f^{(1)}$")
[Link](2, 2, 3)
[Link](x, -12*x**2 + 6*x + 2, ’--’)
[Link]("$f^{(2)}$")
[Link](2, 2, 4)
[Link](x, -24*x + 6, ’:’)
[Link]("$f^{(3)}$")
plt.tight_layout()
[Link]("1D_subplot2.pdf")
[Link]()

44
 50
f 300 f(1)
0 200
50
100 100
150
200 0
250 100
300
350 200
4 3 2 1 0 1 2 3 4 4 3 2 1 0 1 2 3 4

50 f(2) 150 f(3)

0 100
50 50
100
150 0
200 50
250 100
4 3 2 1 0 1 2 3 4 4 3 2 1 0 1 2 3 4

Matplotlib dispose d’autres commandes de tracé comme loglog, polar...

45
Méthodes numériques

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

C’est la méthode la plus précise et la plus utilisée, elle est d’ordre quatre.
L’intervalle [a,b] est divisé en n sous intervalles de largeur h, cette formule s’écrit :
𝒉
𝒚𝒊+𝟏 = 𝒚𝒊 + 𝟔 (𝑲𝟏 + 𝟐(𝑲𝟐 + 𝑲𝟑 ) + 𝑲𝟒 )

𝑲𝟏 = 𝒇(𝒕𝒊 , 𝒚𝒊 )
𝒉 𝒉
𝑲𝟐 = 𝒇 (𝒕𝒊 + , 𝒚𝒊 + 𝑲𝟏 )
𝟐 𝟐
𝒉 𝒉
𝑲𝟑 = 𝒇 (𝒕𝒊 + , 𝒚𝒊 + 𝑲𝟐 )
𝟐 𝟐
𝑲𝟒 = 𝒇(𝒕𝒊 + 𝒉, 𝒚𝒊 + 𝒉𝑲𝟑 )

Exemple : Résoudre le problème de Cauchy précédant par la méthode de Runge-

Kutta d’ordre 4 en prenant un pas d’intégration h=0.25.

𝒉
𝒚𝒊+𝟏 = 𝒚𝒊 + 𝟔 (𝑲𝟏 + 𝟐(𝑲𝟐 + 𝑲𝟑 ) + 𝑲𝟒 )

𝑲𝟏 = 𝒇(𝒕𝒊 , 𝒚𝒊 ) = 𝟐 − 𝒕𝒊 𝒚𝟐𝒊
𝟐
𝒉 𝒉 𝒉 𝒉
𝑲𝟐 = 𝒇 (𝒕𝒊 + , 𝒚𝒊 + 𝑲𝟏 ) = 𝟐 − (𝒕𝒊 + ) (𝒚𝒊 + 𝑲𝟏 )
𝟐 𝟐 𝟐 𝟐
𝟐
𝒉 𝒉 𝒉 𝒉
𝑲𝟑 = 𝒇 (𝒕𝒊 + , 𝒚𝒊 + 𝑲𝟐 ) = 𝟐 − (𝒕𝒊 + ) (𝒚𝒊 + 𝑲𝟐 )
𝟐 𝟐 𝟐 𝟐
𝑲𝟒 = 𝒇(𝒕𝒊 + 𝒉, 𝒚𝒊 + 𝒉𝑲𝟑 ) = 𝟐 − (𝒕𝒊 + 𝒉)(𝒚𝒊 + 𝒉𝑲𝟑 )𝟐
i=0, 𝑲𝟏 = 𝟐, 𝑲𝟐 = 𝟏. 𝟖𝟎𝟒𝟕, 𝑲𝟑 = 𝟏. 𝟖𝟏𝟐𝟐, 𝑲𝟒 = 𝟏. 𝟒𝟕𝟐𝟐, 𝒚𝟏 = 𝟏. 𝟒𝟒𝟔𝟏
y2=1.7028, y3=1.7317, y4=1.6148

2.0

1.8

1.6
Euler Ordre1
y

Euler A Ordre2
1.4 RK4 Ordre4
RK6 Ordre6

1.2

1.0
0.0 0.2 0.4 0.6 0.8 1.0
t

Fig. 4.2. Comparaison entre différentes méthodes

46