Calcul Numérique

Ce cours est consacré à la présentation de quelques méthodes usuelles d’analyse numérique. Il aborde
des techniques de dérivation (intégration) numériques, la résolution des équations non linéaires et la
résolution des équations différentielles ordinaires.

1. Préliminaires
Les nombres décimaux sont représentés en Python par une mantisse (partie fractionnaire) et un exposant.
Considérons la représentation de 2100 en entier et en notation scientifique :

[1]: x = 2**100
y = 2.**100
print("x={} et y={}".format(x,y))

x=1267650600228229401496703205376 et y=1.2676506002282294e+30
La mantisse de 2100 est alors 1.2676506002282294 et son exposant est +30. Pour éviter de manipuler
des représentations différentes d’un même nombre, on convient de placer la virgule décimale après le
premier chiffre non nul, d’où le nom notation scientifique. Les nombres réels sont représentés en Python
avec le type double précision de la norme IEEE754. Les valeurs possibles pour l’exposant sont limités ;
approximativement, en deçà de 10−323 le nombre est arrondi à 0 ; au delà de 10308 , il y a dépassement de
capacité : le nombre est considéré comme infini (inf = ∞).

[2]: x = 1e-324
y = 1e-323
z = 1e308
w = 1e309
print("x={} ; y={} ; z={} et w={}".format(x, y, z, w))

x=0.0 ; y=1e-323 ; z=1e+308 et w=inf

II existe également des limitations pour la mantisse ; les flottants sont représentés avec 15 significatifs
(Exercice : taper un algorithme pour le vérifier). Plus précisément, le module sys permet d’accéder à la
constante epsilon (ε) qui donne la différence entre 1.0 et le flottant suivant le plus proche :

[3]: from sys import float_info as fi

x = fi.epsilon
print("x={}".format(x))

x=2.220446049250313e-16
 Calcul Numérique 2. Résolution d’équation non linéaire 2

Notons que cette constante peut-être utilisée pour écrire une fonction qui teste l’égalité entre deux flottants
dans la limite de la précision de l’interpréteur. En plus, un nombre réel x sera donc stocké en mémoire
sous forme arrondie x̃ vérifiant :
x̃ = x(1 + εx ) avec |εx | 6 ε
Le nombre εx est alors appelé erreur relative et le nombre xεx erreur absolue.
En conclusion, les problèmes d’arrondis sont à prendre en compte dès que l’on travaille avec des nombres
à virgule flottante. En général, pour éviter les pertes de précision, on essayera autant que peut se faire de
respecter les règles suivantes : éviter d’additionner ou de soustraire deux nombres d’ordres de grandeur
très différente ; éviter de soustraire deux nombres presque égaux.
De plus, pour calculer une somme de termes ayant des ordres de grandeur très différents, on applique le
principe dit photo de classe : les petits devant, les grands derrière. Par exemple, si on souhaite donner
+∞
X 1 π4
une approximation décimale de ζ(4) = dont la valeur exacte est , on constate qu’on obtient
n=1
n 4 90
une meilleure précision en commençant par sommer les termes les plus petits :

[4]: from math import *

def zeta(n,p):
return sum(1/i**p for i in range(1, n+1))

def zeta_inverse(n,p):
return sum(1/i**p for i in range(n, 0, -1))

n, p = 2**18, 4
print('L'erreur du plus petit au plus grand :',
"{: .16f}".format(abs(zeta_inverse(n, p) - pi**4/90)))

print('L'erreur du plus grand au plus petit :',

"{: .16f}".format(abs(zeta(n, p) - pi**4/90)))

L'erreur du plus petit au plus grand : 0.0000000000000002

L'erreur du plus grand au plus petit : 0.0000000000002769

2. Résolution d’équation non linéaire

Pour résoudre numériquement une équation non linéaire de la forme f (x) = 0, il faut au préalable étudier
les variations de la fonction pour isoler les zéros.

Théorème 1.
Soit f : [a, b] → R une fonction continue sur un segment.
Si f (a) · f (b) 6 0, alors il existe ` ∈ [a, b] tel que f (`) = 0.

On suppose dans la suite que cette étude est déjà réalisée et que l’on dispose d’un intervalle [a, b] sur
lequel la fonction est strictement monotone, et tel que l’on ait f (a) ∗ f (b) < 0. D’après le théorème de la
bijection monotone, on en déduit que f admet un unique zéro dans ]a, b[. Pour programmer les différentes
méthodes numériques de résolution, nous allons utiliser la classe parente suivante :

[5]: # *******************************************
# Classe Parente
# *******************************************
Calcul Numérique 2. Résolution d’équation non linéaire 3

class Nonlineaire:
def __init__(self, approx, eps=1e-14, maxIter=10**7):
self.approx = approx
self.eps = eps
self.maxIter = maxIter # borne max pour limiter les itérations

def resolution(self, f):

raise NotImplementedError()

2.1. Méthode de dichotomie

On part d’une fonction f : [a, b] → R continue, avec a < b, et f (a) · f (b) 6 0. Voici la première étape de
la construction : on regarde le signe de la valeur de la fonction f appliquée au point milieu a+b 2 .

2 ) 6 0, alors il existe c ∈ [a, 2 ] tel que f (c) = 0.

• Si f (a) · f ( a+b a+b

• Si f (a) · f ( a+b
2 ) > 0, cela implique que f ( 2 ) · f (b) 6 0, et alors il existe c ∈ [ 2 , b] tel que
a+b a+b

f (c) = 0.
y
y

a+b
a 2
b x
f ( a+b )>0
2
f ( a+b
2
)<0
a
a+b b x
2

Nous avons obtenu un intervalle de longueur moitié dans lequel l’équation (f (x) = 0) admet une solution.
On itère alors le procédé pour diviser de nouveau l’intervalle en deux.
Voici le processus complet :
• Au rang 0 : On pose a0 = a, b0 = b. Il existe une solution x0 de l’équation (f (x) = 0) dans l’intervalle
[a0 , b0 ].
• Au rang 1 :
— Si f (a0 ) · f ( a0 +b
2 ) 6 0, alors on pose a1 = a0 et b1 =
0
2 ,
a0 +b0

— sinon on pose a1 = a0 +b0

2 et b1 = b.
— Dans les deux cas, il existe une solution x1 de l’équation (f (x) = 0) dans l’intervalle [a1 , b1 ].
• ...
• Au rang n : supposons construit un intervalle [an , bn ], de longueur 2n ,
b−a
et contenant une solution
xn de l’équation (f (x) = 0). Alors :
— Si f (an ) · f ( an +b
2
n
) 6 0, alors on pose an+1 = an et bn+1 = an +bn
2 ,
— sinon on pose an+1 = an +bn
2 et bn+1 = bn .
— Dans les deux cas, il existe une solution xn+1 de l’équation (f (x) = 0) dans l’intervalle
[an+1 , bn+1 ].
L’erreur absolue de l’approximation de l par an ou bn est alors majorée par bn − an = b−a
2n ; et l’erreur
relative est de l’ordre de bn|c+a
n|
n
.
 Calcul Numérique 2. Résolution d’équation non linéaire 4

[6]: # *******************************************
# Méthode de dichotomie
# *******************************************
class Dichotomie(Nonlineaire):
def subdivise(self, a, b, f):
return 0.5 * (a+b)

def resolution(self, f):

[a, b] = self.approx
a, b = min(a, b), max(a, b)
eps, maxIter = self.eps, self.maxIter
for i in range(maxIter):
c = self.subdivise(a, b, f)
if b-a < eps:
return c
if f(a) * f(c) < 0:
b = c
else:
a = c
return c

2.2. Méthode de régula-falsi

Cette méthode est une variante de la précédente : au lieu de choisir, pour découper l’intervalle [an , bn ], le
milieu cn , on considère l’intersection de la corde joignant les points An (an , f (an )) et Bn (bn , f (bn )) avec
l’axe des abscisses. On pose alors
bn − an
cn = an − f (an )
f (bn ) − f (an )
La convergence est souvent plus rapide qu’avec la méthode de la dichotomie. Pour écrire une classe
RegulaFalsi, il suffit de la faire hériter des méthodes de la classe parente Dichotomie et en surchargeant
le méthode subdivise (c’est un exemple de polymorphisme) :

[7]: # *******************************************
# Méthode regula falsi
# *******************************************
class RegulaFalsi(Dichotomie):
def subdivise(self, a, b, f):
return a - f(a) * (b - a) / (f(b) - f(a))

2.3. Méthode du point fixe

On se propose de résoudre une équation de la forme
f (x) = x.
Les solutions de cette équation sont appelés des points fixes de f ; géométriquement, ce sont les abscisses
des points d’intersection entre la courbe représentative de f et la première bissectrice des axes. Une
méthode fructueuse pour résoudre ce type d’équation est de procéder par itération, à condition toutefois
que soient vérifiées par exemple les hypothèses du théorème du point fixe, dont voici l’énoncé :

Proposition 1.
Soit I un intervalle fermé stable pour une application contractante f de rapport de Lipschitz k ∈ [0, 1[ ;
alors
 Calcul Numérique 2. Résolution d’équation non linéaire 5

• f possède un unique point fixe l ∈ l ;

x0 = α

• pour tout α ∈ I, la suite et on a la majoration
∀n ∈ N xn+1 = f (xn ) converge vers l
∀n ∈ N, |xn − l| 6 k n |x0 − l|

Les points fixes peuvent être classés en différentes catégories suivant la valeur de f 0 (l) en points fixes
attractifs (|f 0 (l)| < 1), répulsifs (|f 0 (l)| > 1), super-attractifs (|f 0 (l)| = 0), indifférents (|f 0 (l)| = 1).

[8]: # *******************************************
# Méthode du point fixe
# *******************************************
class PointFixe(Nonlineaire):
def resolution(self, f):
eps, maxIter = self.eps, self.maxIter
x0, x1 = self.approx, f(self.approx)
for i in range(maxIter):
if abs(x0 - x1) < eps:
return x1
x0, x1 = x1, f(x1)
return x1

2.4. Méthode de Newton

La méthode de Newton est une des méthodes les plus utilisées pour la résolution des équations nonlinéaires ;
en effet, d’une part, elle converge généralement très vite, d’autre part, elle se généralise aux systèmes non
linéaires.
Le principe en est le suivant : partant d’une valeur u0 , on trace la tangente à la courbe représentative de
f (u0 )
f au point M0 (u0 , f (u0 )) ; cette droite coupe l’axe des abscisses en un point d’abscisse u1 = u0 − 0 .
f (u0 )

f (un )

un
un+1

On réitère alors l’opération tant que deux valeurs consécutives sont distantes d’au moins la valeur de
seuil epsilon (ε). Ceci revient à appliquer la méthode du point fixe avec la fonction
f (x)
ϕ(x) = x − .
f 0 (x)
 Calcul Numérique 2. Résolution d’équation non linéaire 6

[9]: # *******************************************
# Méthode de Newton
# *******************************************
class Newton(PointFixe):
def resolution(self, f, fp):
return PointFixe.resolution(self, lambda x:x-f(x)/fp(x))

2.5. Méthode de la sécante

Dans le cas où l’évaluation numérique de la dérivée de f s’avère complexe, on peut modifier la méthode de
Newton en remplaçant le calcul de la dérivée par un taux d’accroissements ; on obtient alors la méthode
de la sécante, qui est une méthode à deux points :
xn − xn−1
xn+1 = xn − f (xn ) ∗
f (xn ) − f (xn−1 )
y

a a0 a00 b
x
A00
A0
A

Proposition 2.
Soit f : [a, b] → R une fonction continue, strictement croissante et convexe telle que f (a) 6 0, f (b) > 0.
Alors la suite définie par
b − an
a0 = a et an+1 = an − f (an )
f (b) − f (an )
est croissante et converge vers la solution ` de (f (x) = 0).

L’hypothèse f convexe signifie exactement que pour tout x, x0 dans [a, b] la sécante (ou corde) entre
(x, f (x)) et (x0 , f (x0 )) est au-dessus du graphe de f .
y

(x0 , f (x0 ))
x
x0 x

(x, f (x))

Dans ce cas, la convergence est moins rapide que dans le cas de la méthode de Newton, mais le résultat
peut être obtenu plus rapidement si le calcul de f 0 est très long.
 Calcul Numérique 2. Résolution d’équation non linéaire 7

[10]: # *******************************************
# Méthode de la sécante
# *******************************************
class Secante(Nonlineaire):
def resolution(self, f):
eps, maxIter = self.eps, self.maxIter
[x0, x1] = self.approx
for i in range(maxIter):
x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))
if abs(x0 - x2) < eps:
return x2
x0, x1 = x1, x2
return x2

2.6. Comparatif des différentes méthodes

Comparons les différentes méthodes de résolution sur la recherche d’une approximation décimale de
√
2. D’une part, on peut appliquer les méthodes précédentes (hormis celle du point fixe) à l’équation
x2 − 2 = 0. D’autre part, on peut appliquer la méthode du point fixe en réécrivant cette équation sous la
1 a √
x = f (x) avec f (x) = (x + ) (Exercice : prouver que f est stable sur ] a, +∞[ ; y est 12 -lipschitzienne
2 x √
et possède sur cet intervalle pour unique point fixe le réel a, qui est un point super-attractif).

[11]: # *******************************************
# Comparatifs des différentes méthodes
# *******************************************
from math import *
print('-'*62)
print("{:>7}|{:^10}|{:^10}|{:^10}|{:^10}|{:^10}"
.format("MaxIter", "Dichotomie", "RegulaFalsi",
"Secante", "PointFixe", "Newton"))
for i in range(1, 9):
di = Dichotomie([1, 2], maxIter=i)
rg = RegulaFalsi([1, 2], maxIter=i)
se = Secante([1, 2], maxIter=i)
pf = PointFixe(2, maxIter=i)
nw = Newton(2, maxIter=i)
comp = [eval(m).resolution(lambda x:x**2-2) for m in ['di', 'rg', 'se']]
comp.append(pf.resolution(lambda x:0.5*(x+2/x)))
comp.append(nw.resolution(lambda x:x**2-2, lambda x:2*x))
comp = [abs(sqrt(2) - c) for c in comp]
print(' '*5, '-'*56)
print("{:>7d}|{: .3e}|{: .4e}|{: .3e}|{: .3e}|{: .3e}"
.format(i, comp[0], comp[1], comp[2], comp[3], comp[4]))
print(' '*5, '-'*56)

--------------------------------------------------------------
MaxIter|Dichotomie|RegulaFalsi| Secante |PointFixe | Newton
--------------------------------------------------------
1| 8.579e-02| 8.0880e-02| 8.088e-02| 2.453e-03| 2.453e-03
--------------------------------------------------------
2| 1.642e-01| 1.4358e-02| 1.421e-02| 2.124e-06| 2.124e-06
--------------------------------------------------------
Calcul Numérique 3. Dérivation numérique 8

3| 2.892e-01| 4.2046e-04| 4.206e-04| 1.595e-12| 1.595e-12

--------------------------------------------------------
4| 2.267e-01| 4.2046e-04| 2.124e-06| 2.220e-16| 0.000e+00
--------------------------------------------------------
5| 1.955e-01| 4.2046e-04| 3.158e-10| 2.220e-16| 2.220e-16
--------------------------------------------------------
6| 1.798e-01| 4.2046e-04| 2.220e-16| 2.220e-16| 2.220e-16
--------------------------------------------------------
7| 1.720e-01| 4.2046e-04| 0.000e+00| 2.220e-16| 2.220e-16
--------------------------------------------------------
8| 1.681e-01| 4.2046e-04| 2.220e-16| 2.220e-16| 2.220e-16
--------------------------------------------------------
La méthode du point fixe converge ici à une vitesse vertigineuse du fait que le point fixe est super-attractif.
À noter que même si la méthode de Newton permet en général d’obtenir une convergence quadratique,
un mauvais choix de la valeur initiale peut provoquer la divergence de cette méthode. D’où l’importance
d’une étude préalable soignée de la fonction f .
Travaux pratiques 1.
1. Faire trois itérations de la méthode de la bissection pour les fonctions suivantes et à partir des intervalles indiqués.
Déterminer le nombre d’itérations nécessaires pour obtenir une solution dont le chiffre des millièmes est significatif.
(a) f (x) = x2 | sin x| − 4.1 dans l’intervalle [0, 4]
(b) f (x) = x6 − x − 1 dans l’intervalle [1, 2]
2. Calculer les points fixes des fonctions suivantes et vérifier s’ils sont attractifs ou répulsifs.
(a) f (x) = arcsin(x)
(b) f (x) = 5 + x − x2
3. Faire 5 itérations de la méthode de la sécante pour les fonctions :
(a) f (x) = x3 − 2x2 − 5, avec x0 = 3,
(b) f (x) = x6 − x − 1, avec x0 = 1.5.
Dans chaque cas, prendre le x0 donné et poser x1 = x0 + 1.
4. On considère l’équation :
th(x) = x2 − 1
(a) Montrer que cette équation possède deux racines assez proches de 0.
(b) En partant de x0 = 2, écrire un programme afin de comparer expérimentalement le comportement des deux
méthodes itératives suivantes.
p  
xn+1 = 1 + th(xn ), xn+1 = arcsinh cosh(xn ) ∗ (x2n − 1) .
(c) Que se passe-t-il si on part de x0 = 0 ?
(d) Déterminer expérimentalement un intervalle I = [0, a] pour lequel la première méthode itérative converge
quel que soit x0 ∈ I.
5. Évaluer la quantité :
q
3
p3 √
3
s= 3+ 3 + 3 + ...
Suggestion : mettre cette équation au cube et obtenir une équation de la forme f (s) = 0. Résoudre cette dernière
à l’aide de la méthode de Newton à partir de s0 = 1.

3. Dérivation numérique
Lorsqu’une fonction est trop compliquée pour que le calcul analytique de sa dérivée soit pratique, on peut
obtenir une valeur approchée de la dérivée en calculant un taux d’accroissement τx (h) :
f (x + h) − f (x)
f 0 (x) ≈ := τx (h)
h
En théorie, plus h est petit, plus τx (h) est proche de f 0 (x). Mais, eu égard aux problèmes d’arrondi, il ne
faut pas prendre h trop petit. Pour une fonction suffisamment régulière, il est donc judicieux de choisir
√
une valeur de h qui soit de l’ordre de ε, c’est-à-dire de l’ordre de 10−8 .
 Calcul Numérique 3. Dérivation numérique 9

L’écriture d’une classe pour la dérivation numérique peut être :

[12]: # ***********************************
# Classe Parente
# ***********************************
class Diff:
def __init__(self, f, h=1e-8):
self.f = f
self.h = h
# ***********************************
# Approx. en avant d'ordre 1
# ***********************************
class DiffAvant1(Diff):
def __call__(self, x):
f, h = self.f, self.h
return (f(x+h)-f(x))/(h)

Comparons pour différentes valeurs de h les valeurs données pour la dérivée de la fonction sinus en 0 :

[13]: # ***********************************
# Exemple : dérivée de sinus en zéro
# ***********************************
from math import *
print('- '*15); print('{: ^8} | {: ^20}'.format(' h ', 'cos(0)')); print('- '*15)
for n in range(1, 13):
x, h = 0, 10**(-n); cos = DiffAvant1(sin, h)
print('{: 7.1e} | {: .15f}'.format(h, cos(x)))

- - - - - - - - - - - - - - -
h | cos(0)
- - - - - - - - - - - - - - -
1.0e-01 | 0.998334166468282
1.0e-02 | 0.999983333416666
1.0e-03 | 0.999999833333342
1.0e-04 | 0.999999998333333
1.0e-05 | 0.999999999983333
1.0e-06 | 0.999999999999833
1.0e-07 | 0.999999999999998
1.0e-08 | 1.000000000000000
1.0e-09 | 1.000000000000000
1.0e-10 | 1.000000000000000
1.0e-11 | 1.000000000000000
1.0e-12 | 1.000000000000000
On constate qu’à partir de h = 10−8 la valeur donnée est exacte. Notons qu’il existe d’autres formules
d’approximation de la dérivée, basées sur des calcul de différences finies et qui permettent d’obtenir un
ordre d’approximation plus élevé :
f (x) − f (x − h)
f 0 (x) ≈ approximation en arrière d’ordre 1
h
f (x + h) − f (x − h)
f 0 (x) ≈ approximation centrée d’ordre 2
2h
4f (x + h) − f (x + 2h) − 3f (x)
f 0 (x) ≈ approximation en avant d’ordre 2
2h
−4f (x − h) + f (x − 2h) + 3f (x)
f 0 (x) ≈ approximation en arrière d’ordre 2
2h
 Calcul Numérique 3. Dérivation numérique 10

etc.
Pour les programmer, il suffit d’ajouter une classe pour chacune dans le fichier parent.py :

[14]: # ***********************************
# Approx. centré d'ordre 2
# ***********************************
class DiffCentre2(Diff):
def __call__(self, x):
f, h = self.f, self.h
return (f(x+h)-f(x-h))/(2*h)

Comparons, par exemple, l’efficacité de la formule de la différence centrée d’ordre 2 pour le calcul de la
dérivée de la fonction racine carrée :

[15]: # ***********************************
# Exemple : comparaison
# ***********************************
def f(x): return sqrt(1+x)

print('- '*27)
print('{: ^8} | {: ^20} | {: ^20}'
.format(' h ', 'DiffAvant1', 'DiffCentre2'))
print('- '*27)
for n in range(1, 13):
x, h = 0, 10**(-n)
g1 = DiffAvant1(f, h)
g2 = DiffCentre2(f, h)
print('{: 7.1e} | {: .17f} | {: .17f}'.format(h, g1(x), g2(x)))

- - - - - - - - - - - - - - - - - - - - - - - - - - -
h | DiffAvant1 | DiffCentre2
- - - - - - - - - - - - - - - - - - - - - - - - - - -
1.0e-01 | 0.48808848170151631 | 0.50062775059818931
1.0e-02 | 0.49875621120889502 | 0.50000625027344925
1.0e-03 | 0.49987506246096380 | 0.50000006250000562
1.0e-04 | 0.49998750062396624 | 0.50000000062444538
1.0e-05 | 0.49999875000317212 | 0.50000000000327560
1.0e-06 | 0.49999987505877641 | 0.50000000001437783
1.0e-07 | 0.49999998807948032 | 0.50000000029193359
1.0e-08 | 0.49999999696126451 | 0.50000000251237964
1.0e-09 | 0.50000004137018550 | 0.50000004137018550
1.0e-10 | 0.50000004137018550 | 0.50000004137018550
1.0e-11 | 0.50000004137018550 | 0.50000004137018550
1.0e-12 | 0.50004445029117051 | 0.50004445029117051
Cette fois-ci, on voit apparaître très nettement la perte de précision lorsque h est trop petit.
Travaux pratiques 2.
1. Évaluer la dérivée de f (x) = ex en x = 0 à l’aide des différences avant et arrière d’ordre 2. Prendre h = 0.05 et
h = 0.025 et calculer le rapport des erreurs commises.
2. À partir des données suivantes :
Fonction tabulée
x 1.00 1.01 1.02
f (x) 1.27 1.32 1.38
calculer les approximations de f 0 (1.005), f 0 (1.015) et f 00 (1.01) à l’aide des formules de différences centrées. Si les
 Calcul Numérique 4. Intégration numérique 11

valeurs de f (x) sont précises à ±0.005, quelle est l’erreur maximale sur chacune des approximations ?
√

3. Soit f (x) = 12 x 1 − x2 + arcsin(x) . Pour x0 = 0.5, calculer des approximations de f 0 (x0 ) par la formule de
différence arrière d’ordre 2 avec h = 0.1, h = 0.05, h = 0.025, h = 0.0125 et h = 0, 00625. Dans chacun des cas,
calculer l’erreur exacte en valeur absolue. Cette erreur se comporte-telle en O(h2 ) ?
4. Obtenir la formule de différence centrée d’ordre 4 permettant de faire l’approximation de f 00 (x).

4. Intégration numérique
Dans cette section, nous nous intéresserons aux algorithmes permettant le calcul numérique d’intégrale.
Le principe général commun aux méthodes que nous allons présenter est simple : découper l’intervalle
d’intégration en petits morceaux, et approcher sur chacun de ces petits intervalles la courbe représentative
de la fonction f par une courbe très simple pour laquelle le calcul d’aire est facile. Avant de décrire
quelques méthodes d’intégration numérique, on peut commencer par écrire une super-classe (la classe
Integration).

[16]: # ************************************
# Classe Parente
# ************************************
class Integration:
def __init__(self, a, b, n):
self.a, self.b, self.n = a, b, n
self.poids_pivots = self.quadrature()

def quadrature(self):
raise NotlmplementedError

def integrale(self, f):

return sum(w * f(x) for w, x in self.poids_pivots)

4.1. Méthode des rectangles.

Soit f une fonction continue sur un segment [a, b]. Pour calculer son intégrale approchée par la méthode
b−a
des rectangles, on pose h = et on pose xk = a + k ∗ h, pour tout entier k allant de 0 à n (avec
n
n
X
a0 = a et an = b). On pose ensuite Sn (f ) = h f (xk ).
k=0

y y = ex

x
0 n = 10 1

En théorie plus n est grand plus l’approximation est précise (mais plus le pas est petit !).
 Calcul Numérique 4. Intégration numérique 12

[17]: # ************************************
# Méthode des rectangles
# ************************************
class Rectangles(Integration):
def quadrature(self):
a , b, n = self.a, self.b, self.n
h = (b - a) / n
return [(h, a + k * h) for k in range(n)]

4.2. Méthodes des trapèzes

Le principe général est exactement le même que celui de la méthode des rectangles, mais on approche
cette fois-ci la courbe sur le segment [xk , xk+1 ] par le segment de droite reliant les deux points de la
courbe d’abscisses xk et xk+1 , ce qui revient à calculer une somme d’aires de trapèzes :
n−1
X f (xk ) + f (xk+1 )
Sn (f ) = h .
2
k=0

Exercice : Montrer qu’on peut utiliser des rectangles pour calculer l’aire du trapèze.

[18]: # ************************************
# Méthode des trapèzes
# ************************************
class Trapezes(Integration):
def quadrature(self):
a , b, n = self.a, self.b, self.n
h = (b - a) / n
return [(h, a + (k + 0.5)* h) for k in range(n)]

4.3. Méthode de Gauss

Les méthodes de quadrature de Gauss consistent à choisir non seulement les poids, mais aussi les pivots,
pour que la méthode soit aussi précise que possible. Pour trouver la répartition optimale des pivots, il faut
faire le détour par les polynômes orthogonaux (dont l’exposé théorique ne sera pas détaillé ici !). Pour
appliquer cette méthode au calcul d’une intégrale sur un segment quelconque, on effectue un changement
de variable affine :
n
b + a b + a
Z b
b − a 1 b − a b−aX b − a
Z
f (x)dx = f x+ dx ≈ wj f x+
a 2 −1 2 2 2 j=1 2 2
 Calcul Numérique 4. Intégration numérique 13

[19]: # ************************************
# Méthode de Gauss-Legendre
# ************************************
class GaussLegendre(Integration):
def gaussPivots(self, m, tol=1e-14):
def legendre(t, m): # calcule le m-ième polynôme de Legendre et son dérivé
p0, p1 = 1, t
for k in range(1, m):
p = ((2*k+1)*t*p1 - k*p0) / (1+k)
p0, p1 = p1, p
dp = m*(p0 - t*p1) / (1 - t**2)
return p, dp

w = [0 for i in range(m)] # poids

x = [0 for i in range(m)] # pivots
nRoots = (m + 1) // 2 # nombres de racines positives
for i in range(nRoots):
t = cos(pi*(i + 0.75) / (m + 0.5)) # racine approchée
for j in range(30):
p, dp = legendre(t, m)
dt = -p / dp; t = t + dt # méthode de Newton
if abs(dt) < tol:
x[i], x[m-i-1] = t, -t
w[i] = 2 / (1 - t**2) / dp**2; w[m-i-1] = w [i]
break
return zip(w, x)

def quadrature(self):
a , b, n = self.a, self.b, self.n
return [[(b - a)*w*0.5, (a + b)*0.5 + (a - b)*0.5*x]
for w, x in self.gaussPivots(n)]

4.4. Comparatif des différentes méthodes

Il ne nous reste plus qu’à comparer la méthode de Gauss-Legendre aux méthodes précédentes, ce que
nous faisons en cherchant à calculer une valeur approchée du nombre π à l’aide de la formule
 Z 21 p √ 
3
π = 12 1 − t dt −
2
0 8

[20]: # ************************************
# Comparatif des trois méthodes
# ************************************
from math import *
print('- '*27)
print("{:>3} | {:^15} | {:^14} | {:^10}".format("n",

"rectangles", "trapèzes", "Gauss-Legendre"))

print('- '*27)
for i in range(1, 8):
r, t = Rectangles(0, 0.5, 2**i), Trapezes(0, 0.5, 2**i)
g = GaussLegendre(0, 0.5, 2**i)
Calcul Numérique 4. Intégration numérique 14

comp = [eval(m).integrale(lambda x:sqrt(1-x**2))

for m in ['r', 't', 'g']]
comp = [abs(pi - 12 * (c - sqrt(3) / 8)) for c in comp]
print("{0:>3d} | {1[0]: .8e} | {1[1]: .8e}"
"| {1[2]: .5e}".format(2**i, comp))

- - - - - - - - - - - - - - - - - - - - - - - - - - -
n | rectangles | trapèzes | Gauss-Legendre
- - - - - - - - - - - - - - - - - - - - - - - - - - -
2 | 1.65068645e-01 | 1.78757927e-02| 4.34283e-04
4 | 9.14722187e-02 | 4.49974053e-03| 1.08247e-07
8 | 4.79859796e-02 | 1.12695453e-03| 7.99361e-15
16 | 2.45564671e-02 | 2.81866525e-04| 3.19744e-14
32 | 1.24191668e-02 | 7.04746518e-05| 1.55431e-14
64 | 6.24482072e-03 | 1.76191647e-05| 2.66454e-14
128 | 3.13121994e-03 | 4.40482253e-06| 3.64153e-14
On remarque également que la précision de la méthode de Gauss-Legendre est excellente avec seulement
8 points de quadrature, mais c’est au prix d’une programmation plus ardue. Au-delà de 8 points, les
erreurs d’arrondi dans la détermination des pivots et des points en gâtent la précision.

Travaux pratiques 3.
Z 2
1. Utiliser une méthode numérique pour évaluer ln(x)dx.
0
Z π
2
2. À l’aide d’une certaine méthode d’intégration numérique, on a évalué sin xdx en utilisant 3 valeurs différentes
0
de h. On a obtenu les résultats suivants.
Valeurs de I
h I
0.1 1.001 235
0.2 1.009 872
0.3 1.078 979
Compte tenu de la valeur exacte de I, déduire l’ordre de convergence de la quadrature employée.
3. Soit une fonction aléatoire X suivant une loi normale. La probabilité que X soit inférieure ou égale à x (notée
P(X 6 x)) est donnée par la fonction :
Z x
1
t2
P X6x = √ e− 2 dt
2π −∞
t2
Comme la fonction e− 2 n’a pas de primitive, on calcule cette probabilité pour différentes valeurs de x et l’on
garde les résultats dans des tables comme celle-ci :
Valeurs de P(X 6 x)
x P(X 6 x)
1.0 0.841 3447
1.1 0.864 3339
1.2 0.884 9303
1.3 0.903 1995
1.4 0.919 2433
(a) Calculer, toujours en vous servant de la table fournie, la dérivée P0 (X 6 1.2) avec une précision d’ordre 4.
Comparer avec la valeur exacte de cette dérivée.
(b) Calculer, toujours en vous servant de la table fournie, la dérivée P00 (X 6 1.2) avec une précision d’ordre 2
 Calcul Numérique 5. Résolution numérique d’équations différentielle 15

5. Résolution numérique d’équations différentielle

5.1. Position du problème

Dans cette section, on s’intéresse aux équations différentielles d’ordre 1 dont la fonction inconnue est soit
à valeurs réelles, soit à valeurs vectorielles.
Soit U un ouvert de R × Rn et f : U −→ R une application continue sur l’ouvert U . Alors pour (t0 , X0 )
de U , on appelle solution du problème de Cauchy :
x (t) = f (t, x(t))
 0
(1)
x(t0 ) = x0
tout couple (I, x) où I est un intervalle contenant t0 et x : I −→ Rn une solution sur I de l’équation
différentielle x0 (t) = f (t, x(t)) telle que x(t0 ) = x0 .
L’équation différentielle qui apparaît dans (1) sera appelée équation différentielle scalaire si n = 1 , et
système différentiel sinon.
Si l’application f : U −→ Rn est de classe C n sur l’ouvert U , alors, d’après le théorème de Cauchy-
Lipschitz , pour (t0 , x0 ) ∈ U , il existe une unique solution maximale au problème de Cauchy (1). Dans la
plupart des cas, on ne sait pas résoudre explicitement le problème de Cauchy ; d’où la nécessité de mettre
au point des méthodes numériques de résolution approchée d’un tel problème.
Nous allons créer une classe ODESolver qui servira d’interface pour la programmation des différents
schémas numériques. Cette classe possède comme attribut la fonction f , le pas h, les conditions initiales
(t0 , x0 ) ; de plus, elle programme la résolution de l’équation de départ sur un intervalle I = [a, b] ; pour ce
faire, il faut parcourir l’intervalle I à partir de t0 vers la droite jusqu’à atteindre b, puis à partir de t0
vers la gauche jusqu’à atteindre a.

[21]: # ************************************
# Fonction auxiliaire
# ************************************
def terminate(u, t, step_no):
eps = 1.0E-6 # nombre assez petit
return abs(u[step_no]) < eps # assez proche de zero?
# ************************************
# Classe parente
# ************************************
class ODESolver:
def __init__(self, f):
self.f = lambda u, t: np.asarray(f(u, t), float)

def advance(self):
raise NotImplementedError

def set_initial_condition(self, U0):

if isinstance(U0, (float,int)): # EDO scalaire
self.neq = 1
U0 = float(U0)
else: # Système EDO
U0 = np.asarray(U0)
self.neq = U0.size
self.U0 = U0
 Calcul Numérique 5. Résolution numérique d’équations différentielle 16

def solve(self, time_points, terminate=None):

if terminate is None:
terminate = lambda u, t, step_no: False

self.t = np.asarray(time_points)
n = self.t.size
if self.neq == 1: # EDO scalaire
self.u = np.zeros(n)
else: # Système EDO
self.u = np.zeros((n,self.neq))

# Supposons que self.t[0] est celui de self.U0

self.u[0] = self.U0

# La boucle temporelle
for k in range(n-1):
self.k = k
self.u[k+1] = self.advance()
if terminate(self.u, self.t, self.k+1):
break # terminer la boucle sur k
return self.u[:k+2], self.t[:k+2]

5.2. Méthode d’Euler

La méthode d’Euler est la méthode la plus simple pour résoudre numériquement une équation différentielle.
Elle présente un réel intérêt théorique, puisqu’elle peut être utilisée pour démontrer le théorème de
Cauchy-Lipschitz. Toutefois, son intérêt pratique est limité par sa faible précision.
x(t + h) − x(t)
L’idée d’Euler consiste à utiliser l’approximation : x0 (t) ≈ pour h assez petit. Autrement
h
dit : x(t + h) ≈ x(t) + x (t) ∗ h = x(t) + f (t, x(t)) ∗ h. On construit ainsi une suite de points de la manière
0

suivante :
∀k ∈ N, tk+1 = tk + h et xk+1 = xk + f (t, xk ) ∗ h
Voici une manière de programmer ce schéma numérique :

[22]: # ***********************************
# Euler Forward
# ***********************************
class ForwardEuler(ODESolver):
def advance(self):
u, f, k, t = self.u, self.f, self.k, self.t
h = t[k+1] - t[k]
u_new = u[k] + h*f(u[k], t[k])
return u_new

Testons l’efficacité de cette méthode pour le problème : x0 = x avec x(0) = 1 dont la solution est ex .

[23]: # ***********************************
# Test via y'=y avec y(0)=1 sur [0,3]
# ***********************************
import numpy as np
import matplotlib.pyplot as plt
from math import *
 Calcul Numérique 5. Résolution numérique d’équations différentielle 17

def f(u, t):

return u

T = 3
for h in 0.8, 0.5, 0.2, 0.1: # Pour différents dt (pas)
n = int(round(T/h))
solver = ForwardEuler(f)
solver.set_initial_condition(1)
u, t = solver.solve(np.linspace(0, T, n+1))
plt.plot(t, u, label = 'pas=%g' % h)
plt.legend(loc="upper left")

x = np.linspace(0, 3, 31)
plt.plot(x, np.exp(x), color='black',label='exacte')
plt.legend(loc="upper left")
plt.grid()

Il faut bien comprendre qu’à chaque étape, on repart dans la direction d’une solution exacte de l’équation
différentielle, mais qui n’est pas celle qui est du problème Cauchy initial. Sur le graphe ci-dessous, on trace
la solution exacte pour la condition initiale x(0) = 1 et les courbes intégrales qui passent par (tk , xk ).

[24]: # ***********************************
# Courbes intégrales
# ***********************************
T = 3
h = 0.1 # pas fixe
n = int(round(T/h))
solver = ForwardEuler(f)
for j in [0.2, 0.4, 0.6, 0.8, 1.2, 1.4, 1.6, 1.8]: # points variables
solver.set_initial_condition(j)
u, t = solver.solve(np.linspace(0, T, n+1))
plt.plot(t, u, color = 'red')
x = np.linspace(0, 3, 31)
 Calcul Numérique 5. Résolution numérique d’équations différentielle 18

plt.plot(x, np.exp(x), color='black')

plt.show()

Pour juger de la qualité d’une méthode (ou schéma) numérique de résolution d’équations différentielles, il
faut prendre en compte plusieurs critères.
L’erreur de consistance donne l’ordre de grandeur de l’erreur effectué à chaque pas. Elle mesure l’erreur
qu’entraîne le fait d’approcher le nombre dérivé par un taux d’accroissement. À l’aide de la formule de
Taylor-Lagrange, on peut montrer que dans le cas de la méthode d’Euler, l’erreur de consistance est
dominée par h2 ; on dit alors que cette méthode est d’ordre 1. Ainsi, d’un point de vue théorique, plus le
pas est petit, meilleure sera l’approximation. Un schéma numérique est dit consistant si la somme des
erreurs de consistance tend 0 quand h → 0 (ce qui est le cas ici).
La stabilité contrôle la différence entre deux solutions approchées correspondant à deux conditions initiales
voisines : un schéma est dit stable si un petit écart entre les conditions initiales x(t0 ) = x0 et x̃(t0 ) = x0 +ε
et de petites erreurs d’arrondi dans le calcul récurrent des provoquent une erreur finale |x(t0 ) − x̃(t0 )|
contrôlable . La méthode d’Euler est stable.
Un schéma est dit convergent lorsque l’erreur globale (qui est le maximum des écarts entre la solution
exacte et la solution approchée) tend vers 0 lorsque le pas tend vers 0. En fait, pour que le schéma soit
convergent, il suffit que la méthode soit consistante et stable.
Enfin pour la mise en application du schéma, il faut aussi prendre en compte l’influence des erreurs
d’arrondi. En effet, afin de minimiser l’erreur globale théorique, on pourrait être tenté d’appliquer la
méthode d’Euler avec un pas très petit, mais ce faisant, outre que le temps de calcul deviendrait irréaliste,
les erreurs d’arrondi feraient diverger la solution approchée très rapidement. La question de trouver de
manière théorique le pas optimal peut s’avérer un problème épineux.
Dans le script suivant, on applique la méthode d’Euler sur l’intervalle [0, 1] avec la condition initiale
x(0) = 1. Pour minimiser le temps de calcul, l’algorithme est réécrit sous la forme suivante.

[25]: # ***********************************
# Exemple
# ***********************************
from math import exp
 Calcul Numérique 5. Résolution numérique d’équations différentielle 19

print('-'*54)
print("{: ^8} | {: ^12} | {: ^12} | {: ^13}".format('h', 'x(t)', 'exp(t)',␣
,→'erreur'))

print('-'*54)
for k in range(1, 10): # Pour k>=6, l'algo est trop lent
h, t, x = 10**(-k), 0, 1
while t < 1:
t, x = t + h, x * (h + 1)
print("{0:>8g} | {1:>.10f} | {2:>.10f} | {3:>.10f}".format(h, x, exp(t), exp(t)␣
,→- x))

------------------------------------------------------
h | x(t) | exp(t) | erreur
------------------------------------------------------
0.1 | 2.8531167061 | 3.0041660239 | 0.1510493178
0.01 | 2.7048138294 | 2.7182818285 | 0.0134679990
0.001 | 2.7169239322 | 2.7182818285 | 0.0013578962
0.0001 | 2.7184177414 | 2.7185536702 | 0.0001359288
1e-05 | 2.7182954199 | 2.7183090114 | 0.0000135915
1e-06 | 2.7182804691 | 2.7182818285 | 0.0000013594
1e-07 | 2.7182819660 | 2.7182820996 | 0.0000001336
1e-08 | 2.7182817983 | 2.7182818347 | 0.0000000363
1e-09 | 2.7182820738 | 2.7182818299 | -0.0000002438
Pour un pas h 6 10−9 le temps de calcul devient trop important, de plus, on voit poindre l’effet des
erreurs d’arrondis, puisque la solution approchée est devenu supérieure à l’exponentielle, ce qui est
mathématiquement faux, en vertu de l’inégalité :
 1 n
∀n ∈ N, 1+ 6e
n
En général, il ne suffit pas qu’un schéma numérique soit convergent pour qu’il donne de bons résultats
sur n’importe quelle équation différentielle. Encore faut-il que le problème soit bien posé (hypothèses de
Cauchy-Lipschitz doivent être vérifiées !) et bien conditionné (une petite erreur sur l’un des des termes
aura des répercussions allant en s’amplifiant).

5.3. Méthode de Runge-Kutta d’ordre 4

Si, la méthode d’Euler présente un intérêt théorique, on préfère en pratique des méthodes d’ordre plus élevé.
Parmi la multitude des schémas numériques celle qui présente le meilleur rapport précision/complexité
est certainement celle de Runge-Kutta d’ordre 4. En voici le schéma pour tout entier k allant de 0 à N ,

= h ∗ f (tn , xn )

 k1
tk+1 = tk + h = h ∗ f tn + h2 , xn + k21
 

k2

où
xk+1 = xk + 16 (k1 + 2k2 + 2k3 + k4 ) = h ∗ f tn + h2 , xn + k22


 k3

= h ∗ f (tn + h, xn + k3 )

k4
Sans plus attendre, écrivons ce schéma.

[26]: # ***********************************
# Méthode de Runge-Kutta d'ordre 4
# ***********************************
class RungeKutta4(ODESolver):
def advance(self):
 Calcul Numérique 5. Résolution numérique d’équations différentielle 20

u, f, k, t = self.u, self.f, self.k, self.t

dt = t[k+1] - t[k]
dt2 = dt/2.0
K1 = dt*f(u[k], t[k])
K2 = dt*f(u[k] + 0.5*K1, t[k] + dt2)
K3 = dt*f(u[k] + 0.5*K2, t[k] + dt2)
K4 = dt*f(u[k] + K3, t[k] + dt)
u_new = u[k] + (1/6.0)*(K1 + 2*K2 + 2*K3 + K4)
return u_new

Effectuons une comparaison des méthodes d’Euler et de Runge-Kutta en les appliquant au problème de
Cauchy : x0 = x et x(0) = 1 sur l’intervalle [0, 1].

[27]: # *******************************************************
# Testons graphiquement ForwardEuler vs RungeKutta4
# *******************************************************
# figure()
# legends = []
T = 3
dt = 0.5
n = int(round(T/dt))
t_points = np.linspace(0, T, n+1)
for solver_class in RungeKutta4, ForwardEuler:
solver = solver_class(f)
solver.set_initial_condition(1)
u, t = solver.solve(t_points)
plt.plot(t, u, label = '%s' % solver_class.__name__)
plt.legend(loc="upper left")
# hold('on')
x = np.linspace(0, 3, 31)
plt.plot(x, np.exp(x), color='black',label='exacte')
plt.legend(loc="upper left")
plt.show()
 Calcul Numérique 5. Résolution numérique d’équations différentielle 21

[28]: # **************************************
# Comparaison pour différents pas
# *************************************
from math import exp
print('-'*51)
print("{: ^10} | {: ^18} | {: ^17}".format('h', 'erreur Euler', 'erreur RK4'))
print('-'*51)
T = 1
for i in range(1, 6):
h = 10**(-i)
n = int(round(T/h))
t_points = np.linspace(0, T, n+1)
solverEuler = ForwardEuler(f)
solverRK4 = RungeKutta4(f)
solverEuler.set_initial_condition(1)
solverRK4.set_initial_condition(1)
u_e, t_r = solverEuler.solve(t_points)
u_r, t_r = solverRK4.solve(t_points)
print("{0:>10g} | {1:>.12e} | {2:>.10e}"
.format(h, exp(1) - u_e[-1], exp(1) - u_r[-1]))

---------------------------------------------------
h | erreur Euler | erreur RK4
---------------------------------------------------
0.1 | 1.245393683590e-01 | 2.0843238793e-06
0.01 | 1.346799903752e-02 | 2.2464297089e-10
0.001 | 1.357896223148e-03 | 2.1760371283e-14
0.0001 | 1.359016338212e-04 | -8.4376949872e-15
1e-05 | 1.359128448319e-05 | -1.2878587086e-14
On constate que pour h = 0.001, la méthode d’Euler ne donne que trois décimales significatives du nombre
e, alors que la méthode de Runge-Kutta en donne quatorze ! On remarque également qu’il est inutile
d’appliquer la méthode de Runge-Kutta avec un pas h > 0.001, puisque dans ce cas les erreurs d’arrondi
prennent le dessus par rapport au gain théorique de précision.