Chapitre 9 : Méthodes Numériques avec Python

IPEST La Marsa

A. AMMAR

31 mai 2021

A. Ammar Représentation de l’information 1 / 19

Différentiation numérique

Soit f une application réelle dérivable sur un intervalle ouvert ]a, b[ non vide. La définition du
nombre de la dérivée f 0 (x ) en x ∈]a, b[ fournit implicitement une méthode de calcul approché de
f 0 (x ).

différence décentrée à droite

Le quotient différentiel est appelé différence décentrée à
droite pour h (pas de l’intervalle) suffisamment petit

f (x + h) − f (x )
f 0 (x ) ≈
h

différence décentrée à gauche

On peut tout aussi bien choisir une différence décentrée à
gauche
f (x ) − f (x − h)
f 0 (x ) ≈
h

Remarque
Dans le cas où f n’est dérivable qu’à droite (resp. gauche) en x , on devra s’obliger à faire tendre h
vers zéro par valeurs supérieures (resp. inférieures).

A. Ammar Représentation de l’information 2 / 19

Différentiation numérique

Différence finie
Si une fonction admet une dérivée à droite et une dérivée à gauche , alors elle admet une dérivée
symétrique appelée aussi différence finie centrée
f (a + h) − f (a)
fd0 (a) = lim
h→0+ h
f (a + h) − f (a) f (a − h) − f (a)
fg0 (a) = lim = lim
h→0 − h h→0+ −h
f (a + h) − f (a) + f (a) − f (a − h)
fs0 (a) = lim
h→0 2h
Ainsi, dans le cas de bonnes fonctions, régulières, pour des raisons qu’il est aisé de vérifier
géométriquement, une meilleure approximation de f 0 (x ) est donnée par la formule de dérivée
symétrique :
f (x + h) − f (x − h)
f 0 (x ) ≈
2h

Application en Python
Voir TD 13 : Exercice 2 . . .

A. Ammar Représentation de l’information 3 / 19

Recherche des zéros d’une fonction : Recherche dichotomique

Principe
a+b
Le principe de la recherche dichotomique consiste à calculer f ( ). Suivant le signe de cette
2
quantité l’une des deux relations est nécessairement vérifiée :
a+b a+b
f (a) × f ( )≤0 ou f( ) × f (b) ≤ 0
2 2
a+b
h i
Le premier cas assure l’existence d’une racine dans l’intervalle a, f ( ) , le second dans
2
a+b
h i
l’intervalle f ( ), b , ramenant dans chacun des deux cas la recherche à un intervalle
2
d’amplitude moitié moins grande.

Algorithme
L’algorithme de recherche consiste donc à itérer
deux suites (un )n∈N et (vn )n∈N définies par les
valeurs initiales u0 = a et v0 = b et la relation de
récurrence :

(un , m) si f (un ) × f (m) ≤ 0
(un , vv ) =
(m, vn ) sinon.

un + vn
avec m = .
2
A. Ammar Représentation de l’information 4 / 19
 Recherche des zéros d’une fonction : Recherche dichotomique
La recherche dichotomique nécessite quatre paramètres : la fonction f , les valeurs de a et b et la
valeur de  (qui par défaut sera égal à 10−12 ).

1 def dicho(f, a, b, epsilon=1e−12):
2 if f(a) ∗ f(b) > 0:
3 return None
4 u, v = a, b
5 while abs(v − u) > 2 ∗ epsilon:
6 w = (u + v) / 2
7 if f(u) ∗ f(w) <= 0:
8 v = w
9 else:
10 u = w
11 return (u + v) / 2
 

Utilisons cet algorithme pour chercher l’unique racine de la fonction sinus entre 3 et 4 :
>>> from numpy import sin
>>> dicho(sin, 3, 4)
3.141592653589214
√
Obtenons maintenant une approximation de 2 :
>>> dicho(lambda x: x * x - 2, 1, 2)
1.4142135623724243
A. Ammar Représentation de l’information 5 / 19
Recherche des zéros d’une fonction : Recherche dichotomique

Utilisation du module scipy

scipy est un ensemble de modules spécialement dédiés au calcul numérique ; il n’est donc pas
étonnant d’y trouver une fonction permettant de réaliser la recherche dichotomique que nous
venons d’étudier. Cette fonction se nomme bisect et se trouve dans le module [Link] .
Comme la majorité des autres fonctions de ce module ne nous intéressent pas pour l’instant, on se
contente de charger en mémoire cette unique fonction :
>>> from [Link] import bisect
Avant de l’utiliser, lisons la documentation qui lui est associée :
>>> help(bisect)
Utilisons cette fonction pour chercher l’unique racine de la fonction sinus entre 3 et 4 :
>>> bisect(sin, 3, 4)
3.141592653589214

A. Ammar Représentation de l’information 6 / 19

Recherche des zéros d’une fonction : Méthode de Newton-Raphson

Principe de la méthode de Newton

• La méthode déjà étudiée : dichotomie, déterminent un zéro de f en lequel se produit un
changement de signe en construisant une suite décroissante d’intervalles contenant ce
dernier.
• La méthode de Newton-Raphson (également appelée méthode de Newton) s’affranchit de
cette contrainte, en procédant par approximations successives pour construire une suite
(xn )n∈N qui converge en général vers un zéro de f : l’équation f (x ) = 0 est remplacée par
l’équation affine φn (x ) = 0, où φn est l’équation de la tangente à la fonction f au point
d’abscisse xn .
Cette tangente a pour équation
y = f (xn ) + f 0 (xn ) × (xn+1 − xn ) donc
f (xn )
φn (xn+1 ) = 0 ⇐⇒ xn+1 = xn − 0 .
f (xn )
La méthode de Newton-Raphson consiste donc
en l’itération d’une suite (xn )n∈N définie par la
donnée d’une valeur x0 et de la relation de
récurrence :
f (xn )
xn+1 = xn −
f 0 (xn )

A. Ammar Représentation de l’information 7 / 19

 Recherche des zéros d’une fonction : Méthode de Newton-Raphson

Mise en œuvre pratique

La recherche dichotomique nécessite quatre paramètres : la fonction f et sa dérivée df , la valeur
initiale x et la valeur de  (qui par défaut sera égal à 10−12 ).

1 def NewtonRaphson(f,df,x,epsilon=1e−12):
2 X = f(x)
3 n = 0
4 while abs(X) > epsilon:
5 Xpoint = df(x)
6 x = x − X/Xpoint
7 n +=1
8 X = f(x)
9 return x, n
 

Utilisons cet algorithme pour chercher l’unique racine de la fonction sinus avec la valeur initiale
x =3:
>>> from numpy import sin, cos
>>> NewtonRaphson(sin,cos, 3)
(3.141592653589793, 3)

A. Ammar Représentation de l’information 8 / 19

Recherche des zéros d’une fonction : Méthode de Newton-Raphson

Utilisation du module [Link]

Dans le paragraphe précédent, nous avons appris que dans le module [Link] se trouve
une fonction nommée bissect qui applique la méthode dichotomique pour trouver un zéro d’une
fonction. Dans ce même module se trouve une fonction nommée newton qui, comme son nom
l’indique, applique la méthode de Newton-Raphson.
Avant de l’utiliser, lisons la documentation qui lui est associée :
>>> help(newton)
Utilisons cette fonction pour chercher l’unique racine de la fonction sinus entre 3 et 4 :
>>> newton(sin, 3, cos)
3.141592653589214

Exercice
Voir TD 13 : Exercice 3 . . .

A. Ammar Représentation de l’information 9 / 19

Valeur approchée d’une intégrale : Méthode des rectangles

Principe
• on fait coïncider le sommet haut gauche du rectangle avec la courbe : c’est la méthode des
rectangles à gauche,
• on fait coïncider le sommet haut droit du rectangle avec la courbe : c’est la méthode des
rectangles à droite,
• on fait coïncider le milieu du coté haut du rectangle avec la courbe : c’est la méthode du
point milieu.

A. Ammar Représentation de l’information 10 / 19

Valeur approchée d’une intégrale : Méthode des rectangles

Algorithmes
1 méthode des rectangles à gauche, on obtient :
Z b n−1
X
f (x )dx ≈ h f (xi ) ; avec xi = a + ih
a i=0

2 méthode des rectangles à droite, on obtient :

Z b n−1
X
f (x )dx ≈ h f (xi ) ; avec xi = a + (i + 1)h
a i=0

3 méthode du point milieu, on obtient :

Z b n−1
X 1
f (x )dx ≈ h f (xi) ; avec xi = a + (i + )h
a
2
i=0

A. Ammar Représentation de l’information 11 / 19

Valeur approchée d’une intégrale : Méthode des trapèzes

Principe et algorithme
L’idée de la méthode des trapèzes est d’approcher l’aire sous la courbe par des trapèzes dont l’aire
est facilement calculables.
1 Choisir le nombre n de trapèzes qu’on veut
sous la courbe.
2 on répartit sur [a, b] n + 1 points de la
façon équitable suivante : x0 = a et
xn+1 = xn + h, h = (b − a)/n.
3 l’aire de chaque petit trapèze est
(base1 + base2) × hauteur
Ai = =
2
h
× [f (xi ) + f (xi−1 )].
2
4 Nous obtenons l’aire recherchée en
sommant l’aire de tous les trapèzes entre a
et b, ce qui nous donne :
Z b n−1
h X
f (x )dx ≈ (f (a) + f (b)) + h f (xi )
a
2
i=1

avec xi = a + ih

A. Ammar Représentation de l’information 12 / 19

Valeur approchée d’une intégrale : Méthode des rectangles

Algorithmes
1 méthode des rectangles à gauche, on obtient :
Z b n−1
X
f (x )dx ≈ h f (xi ) ; avec xi = a + ih
a i=0

2 méthode des rectangles à droite, on obtient :

Z b n−1
X
f (x )dx ≈ h f (xi ) ; avec xi = a + (i + 1)h
a i=0

3 méthode du point milieu, on obtient :

Z b n−1
X 1
f (x )dx ≈ h f (xi) ; avec xi = a + (i + )h
a
2
i=0

A. Ammar Représentation de l’information 13 / 19

 Valeur approchée d’une intégrale : Méthode de Simpson

Principe et algorithme
Dans la méthode des trapèzes, nous avons en fait interpolé f (x ) par une droite entre les points i
et i + h de l’intervalle. Dans la méthode de Simpson, nous n’allons plus interpoler par une droite
mais par un polynôme de degré 2, ce qui va améliorer notre précision.
Plaçons nous autour d’un point x0 appartenant à
l’intervalle [a, b], dans la maille de calcul x0 − h et x0 + h.
Pour un accroissement (x − x0 ), le développement de
Taylor limité au second ordre nous donne :

f (x ) = f (x0 )+(x −x0 )f 0 (x0 )+(1/2)(x −x0 )2 f ”(x0 )+O((x −x0 )3 ).

Nous savons que f 0 (x0 ) = [f (x0 + h) − f (x0 − h)] /2h et

que f ”(x0 ) = [f (x0 + h) − 2f (x0 ) + f (x0 − h)] /h2 . Si
nous remplaçons ces valeurs dans le développement limité
et que l’on intègre entre x0 − h et x0 + h, on obtient l’aire
h
élémentaire (f (x0 + h) + 4f (x0 ) + f (x0 − h)).
3
L’intégrale recherchée s’obtient en sommant toutes les
aires élémentaires :
Z b n−1 n−2
h X X
f (x )dx ≈ [(f (a) + f (b)) + 4 f (xi ) + 2 f (xi )]
a
3
i=1 i=2

avec xi = a + ih.
A. Ammar Représentation de l’information 14 / 19
Valeur approchée d’une intégrale : module [Link]

Le module [Link] contient une fonction nommée quad qui permet le calcul approché
d’une intégrale. Cette fonction, très efficace, utilise différentes routines (la plupart sont basées sur
des méthodes de Gauss, ces méthodes sont hors programme) de manière à obtenir la meilleure
approximation possible. Dans son utilisation la plus simple, cette fonction prend trois paramètres
Rb :
f , a et b et renvoie un tuple (I, e) où I est une valeur approchée de l’intégrale I(f ) = a f (x )dx
et e une estimation de l’erreur |I(f ) − I|.
On notera que les valeurs a et b peuvent éventuellement être égales à ±∞.
>>> import numpy as np
>>> from [Link] import quad
>>> quad(lambda x: [Link](x), 0, [Link])
(2.0, 2.220446049250313e-14)
>>> quad(lambda x: [Link](-x*x), -[Link], [Link])
(1.7724538509055159, 1.4202636780944923e-08)
Rπ
Le premier exemple calcule une valeur approchée de l’intégrale 0
sin(x )dx = 2 ; le second une
R +∞ 2 √
valeur approchée de l’intégrale de Gauss −∞ e −x dx = π

>>> [Link]([Link])
1.7724538509055159

A. Ammar Représentation de l’information 15 / 19

Méthodes d’Euler pour la résolution numérique des équations différentielles
ordinaires

Présentation du problème
Nous allons dans un premier temps nous intéresser aux équations différentielles que l’on peut
mettre sous la forme :
ẋ = f (x (t), t)
où f est une fonction définie sur une partie U de R2 , à valeurs dans R.
Une solution de cette équation différentielle est une fonction x de classe C 1 définie sur un certain
intervalle I de R et à valeurs dans R vérifiant :
• (i) ∀t ∈ I ; (x (t), t) ∈ U ;
• (ii) ∀t ∈ I ; ẋ = f (x (t), t).
Nous allons adjoindre à cette équation différentielle une condition initiale sous la forme d’un
couple (x0 , t0 ) ∈ U et chercher à résoudre le problème de Cauchy suivant :

ẋ = f (x (t), t)
x (t0 ) = x0

Sous certaines conditions sur f que nous ne détaillerons pas, ce problème admet une unique
solution, que nous allons chercher à déterminer numériquement.

A. Ammar Représentation de l’information 16 / 19

Méthodes d’Euler pour la résolution numérique des équations différentielles
ordinaires
Les méthodes que nous allons étudier consistent à subdiviser l’intervalle de temps [t0 , t0 + T ] en
n + 1 points t0 < t1 < · · · < tn = t0 + T puis à approcher la relation :
Z tk+1 Z tk+1
x (tk+1 ) − x (tk ) = ẋ dt = f (x (t), t)dt
tk tk

Méthode d’Euler explicite (progressive)

En posant hk = tk+1 − tk , ceci conduit à définir

une suite de valeurs x0 , x1 , . . . , xn à partir de la
condition initiale x0 et de la relation de
récurrence :

∀k ∈ [0, n − 1], xk+1 = xk + hk f (xk , tk )

On observera qu’en général, seul le premier
point x0 de cette méthode est une valeur
exacte ; les autres points sont calculés à
partir de l’approximation précédente, ce qui progressive
peut conduire la valeur calculée xk à s’écarter
de plus en plus de la valeur exacte x (tk ).
A. Ammar Représentation de l’information 17 / 19
Méthodes d’Euler pour la résolution numérique des équations différentielles
ordinaires

Méthode d’Euler implicite (rétrograde)

La méthode d’Euler implicite
Rt consiste à
approcher l’intégrale t k+1 f (x (t), t)dt par la
k
méthode du rectangle droit, ce qui conduit à
définir la suite (x0 , x1 , . . . , xn ) par les relations : rétrograde

∀k ∈ [0, n − 1], xk+1 = xk + hk f (xk+1 , tk+1 )

On observe que cette relation ne procure pas une

relation explicite de xk+1 puisque ce terme est
aussi présent dans le second membre.
Dans la pratique, la méthode d’Euler implicite se
révèle souvent plus stable que la méthode
explicite : elle est moins précise à court terme,
mais diverge moins rapidement de la solution
exacte que la méthode explicite.

A. Ammar Représentation de l’information 18 / 19

 Méthodes d’Euler pour la résolution numérique des équations différentielles
ordinaires

Module [Link]
La fonction odeint du module [Link] est dédiée à la résolution des systèmes
numériques, mais s’applique aussi aux équations scalaires.

ẏ = 1 + cos(t) ∗ y
y (t0 ) = 0

 25
1 import [Link] as plt
2 import [Link] as sc 20
3 import numpy as np
4 def f(y, t): 15
5 return 1 + [Link](t) ∗ y

y(t)
6 t = [Link](0, 10, 100) 10
7 y = [Link](f, 0, t)
8 [Link](t, y, "−ro") 5
9 [Link]("t"); [Link]("y(t)")
0 [Link]()
  0
0 2 4 6 8 10
t

A. Ammar Représentation de l’information 19 / 19