Frédéric Legrand Licence Creative Commons 1

Calculs d’intégrales : méthode de Romberg

1. Introduction
On s’intéresse au calcul approchée d’une intégrale d’une fonction d’une variable réelle,
sur un intervalle fini où la fonction est définie en tout point. On verra dans un premier
temps la méthode des trapèzes itérative, qui consiste à appliquer la méthode des trapèzes
pour des intervalles de plus en petits jusqu’à satisfaire une tolérance demandée. Dans
un second temps, nous aborderons la méthode de Romberg, qui consiste à faire une
extrapolation polynomiale à partir des approximations obtenues par la méthode des
trapèzes pour différentes largeurs d’intervalles. La méthode de Romberg permet, pour des
fonctions assez régulières, d’obtenir une accélération très importante de la convergence.

2. Méthode des trapèzes

2.a. Méthode
On considère une fonction f de la variable réelle x dont on cherche à évaluer numé-
riquement l’intégrale sur un intervalle fini [a; b] :
Z b
I= f (x) dx (1)
a

La fonction f est supposée définie en tout point de l’intervalle [a; b].

Pour obtenir une approximation de cette intégrale, on considère une subdivision de
l’intervalle [a; b] en N intervalles égaux, chacun de largeur h = (b − a)/N . Le nombre N
est une puissance de 2. On pose donc :

N = 2p p = 0, 1, 2, . . .
b−a
hp =
N
x0 = a
xN = b
xk = x0 + kh k = 0, 1, 2, . . . N

La figure suivante montre un exemple pour p = 2 :

Frédéric Legrand Licence Creative Commons 2

f4

f3
f1 f2

f0

x
x0 x1 x2 x3 x4
a b

La formule des trapèzes consiste à approximer l’intégrale sur chaque intervalle par l’aire
du trapèze :
Z xk+1
1 1
f (x) dx ' hfk + h (fk+1 − fk ) = h (fk + fk+1 ) (2)
xk 2 2
On obtient ainsi une approximation de l’intégrale, pour un indice p :
 
1 1
Ip = hp f0 + f1 + f2 + . . . + fN −2 + fN −1 + fN = hp Sp (3)
2 2
Mis à part le premier et le dernier terme, correspondant à x = a et x = b, il s’agit de
calculer la somme des fk .
Cette approximation de l’intégrale étant calculée, on obtient une meilleure approxi-
mation en divisant h par deux, c’est-à-dire en ajoutant les points milieux des intervalles.
La figure suivante montre le passage de p = 2 à p = 3 :
Frédéric Legrand Licence Creative Commons 3

f8

f7
f6
f2 f3 f4 f5
f1

f0

x3 x4 x5 x6 x7 x
x0 x1 x2 x8
a b

On remarque que la somme Sp peut être réutilisée pour calculer Sp+1 . Il suffit d’ajouter
les points milieux, encadrés sur la figure.
Pour calculer une approximation de l’intégrale, on choisit une tolérance  puis on
procède par itérations. La première somme à calculer est :
1
S0 = (f (a) + f (b)) (4)
2
Connaissant la somme Sp , on calcule la somme Sp+1 en ajoutant la somme des fk des
points milieux. L’itération est stoppée lorsque la valeur absolue de la différence des deux
évaluations consécutives de l’intégrale est inférieure à la tolérance :

|Ip+1 − Ip | <  (5)

2.b. Implémentation
Les fonctions sont définies dans une classe :

class Romberg:
def __init__(self,f,a,b):
self.f = f
self.a = a
self.b = b
self.p = 0
self.N = 1
self.S = 0.5*(f(a)+f(b))
self.h = (b-a)
self.I = self.S*self.h
self.I_last = self.I
self.n_eval = 0
def iteration(self):
Frédéric Legrand Licence Creative Commons 4

self.n_eval += self.N
x = self.a+self.h*0.5
somme = self.f(x)
for k in range(self.N-1):
x += self.h
somme += f(x)
self.N *= 2
self.p += 1
self.S += somme
self.h *= 0.5
self.I_last = self.I
self.I = self.h*self.S

def iterations(self,P):
I = [self.I]
h = [self.h]
while self.p<P:
self.iteration()
I.append(self.I)
h.append(self.h)
return (numpy.array(h),numpy.array(I))

def trapezes(self,epsilon):
I = [self.I]
h = [self.h]
self.iteration()
while abs(self.I-self.I_last)>epsilon:
self.iteration()
I.append(self.I)
h.append(self.h)
return (numpy.array(h),numpy.array(I))

def romberg(self,epsilon):
jmax = 20
A=numpy.zeros((jmax+1,jmax+1))
A[0][0] = self.I
self.iteration()
A[1][0] = self.I
correction = (A[1][0]-A[0][0])/3
A[1][1] = A[1][0] + correction
j = 2
while abs(correction) > epsilon:
self.iteration()
A[j][0] = self.I
for i in range(1,j+1):
correction = (A[j][i-1]-A[j-1][i-1])/(4**i-1)
A[j][i] = A[j][i-1] + correction
j += 1
Frédéric Legrand Licence Creative Commons 5

return (A[0:j-1,0:j-1],A[j-1][j-1])

Dans le constructeur, on fournit la fonction à intégrer et les bornes de l’intervalle.

Les variables h,p,N sont initialisées. On calcule la première somme S, la première valeur
de h, et la première estimation de l’intégrale I.
La fonction iteration() calcule la somme Sp+1 et met à jour les différentes variables,
en particulier l’estimation self.I de l’intégrale.
La fonction iterations(P) effectue les itérations jusqu’à un rang P (pour le nombre
p). Elle renvoie la liste des largeurs h et la liste des estimations successives de l’intégrale.
La fonction trapezes(epsilon) effectue les itérations pour une tolérance donnée.
Elle renvoie la liste des largeurs h et la liste des estimations successives de l’intégrale.
Comme premier exemple, on considère l’intégrale suivante :
Z 1
1
x5 dx = (6)
0 6
import numpy
def f(x):
return x**5
romberg = Romberg(f,0,1)
(h,I) = romberg.trapezes(1e-7)

Voici l’erreur :

print(I[-1]-1.0/6)
--> 2.4835268314093994e-08

et le nombre d’évaluations de la fonction f :

print(romberg.n_eval)
--> 4095

On trace l’erreur pour les différentes approximations de l’intégrale, en fonction du

carré de h :

from matplotlib.pyplot import *

figure(figsize=(8,6))
plot(h*h,numpy.absolute(I-1.0/6),"o")
xlabel("$h^2$",fontsize=18)
ylabel("$e$",fontsize=18)
xscale(’log’)
yscale(’log’)
grid()
Frédéric Legrand Licence Creative Commons 6

La pente égale à 1 montre que l’erreur est proportionnelle à h2 .

3. Méthode de Romberg
3.a. Principe
Considérons l’application de la méthode des trapèzes jusqu’au rang p, par exemple
p = 3, et traçons l’estimation de l’intégrale en fonction de h2 .

romberg = Romberg(f,0,1)
(h,I) = romberg.iterations(3)
figure(figsize=(8,6))
plot(h*h,I,"o")
xlabel("$h^2$",fontsize=18)
ylabel("$I$",fontsize=18)
axis([0,1.5,0,1])
grid()
Frédéric Legrand Licence Creative Commons 7

L’estimation la plus précise de l’intégrale est I3 , obtenue pour N = 23 = 8. Son erreur

peut être lue sur le graphique tracé plus haut. On pourrait obtenir une estimation ayant
une erreur plus faible avec p = 4, ce qui nécessiterait le calcul d’une somme de 23 termes,
avec autant d’évaluations de la fonction f . La méthode de Romberg consiste à utiliser
les valeurs déjà calculées I0 , I1 , I2 , I3 pour obtenir une meilleure estimation de l’intégrale,
sans avoir à faire de nouvelles évaluations de la fonction f . Elle repose sur le fait que
l’erreur s’exprime en fonction de h2 . Une estimation de l’intégrale beaucoup plus précise
est donc obtenue en extrapolant les points (h2p , Ip ) déjà calculés jusqu’à h2 = 0.
La méthode de Romberg est une méthode d’accélération de convergence par extra-
polation, appelée aussi extrapolation de Richardson.

3.b. Extrapolation polynomiale

Une extrapolation polynomiale est effectuée en cherchant le polynôme de Lagrange
passant par les points (h2p , Ip ). Ce polynôme est évalué en h2 = 0. On utilise pour cela
l’algorithme de Neville.
L’algorithme de Neville est développé de manière itérative (et non pas récursive).
Les différents termes à calculer sont placés dans une matrice que l’on génère ligne après
ligne :

A0,0
A1,0 A1,1
(7)
A2,0 A2,1 A2,2
A3,0 A3,1 A3,2 A3,3
Le premier terme correspond à la première estimation de l’intégrale A0,0 = I0 . Le premier
élément de la deuxième ligne est la seconde estimation de l’intégrale A1,0 = I1 . À partir
des deux points (h20 , A0,0 ) et (h21 , A1,0 ), on effectue une extrapolation en h2 = 0 par un
Frédéric Legrand Licence Creative Commons 8

polynôme de degré 1. Le résultat de cette extrapolation est donnée par la relation de

récurrence de Neville :

−h21 A0,0 + h20 A1,0 A10 − A0,0

A1,1 = 2 2
= A1,0 +  2 (8)
h0 − h1 h0
h1
−1
La seconde écriture fait apparaître le terme correctif par rapport à la dernière estimation
de l’intégrale. Si ce terme correctif est inférieur à une tolérance  alors on arrête l’itération.
Dans le cas contraire, on calcule A2,0 = I2 puis on effectue une extrapolation polynomiale
de degré 2. Celle-ci se fait en calculant tout d’abord A2,1 , l’extrapolation de degré 1
obtenue avec A1,0 et A2,0 , puis en calculant A2,2 , l’extrapolation de degré 2 calculée à
partir de A1,1 et A2,1 .
L’application de la formule de récurrence de Neville conduit à :
Aj,i−1 − Aj−1,i−1
Aj,i = Aj,i−1 +  2 pour i = 1, 1, . . . , j (9)
hj−i
hj
− 1
Cette formule permet de calculer la ligne j. Le premier élément de cette ligne est
Aj,0 = Ij . Le dernier élément Aj,j est l’extrapolation de degré j. Si le terme correctif par
rapport à Aj,j−1 a une valeur absolue inférieure à la tolérance alors Aj,j est l’estimation
finale de l’intégrale. Dans le cas contraire, on effectue le calcul de la ligne suivante pour
faire une extrapolation de degré j + 1.
En explicitant hi , la formule de récurence s’écrit :
Aj,i−1 − Aj−1,i−1
Aj,i = Aj,i−1 + pour i = 1, 2, . . . , j (10)
4i − 1

3.c. Implémentation et exemples

Une fonction romberg(epsilon) est ajoutée à la classe Romberg. La fonction renvoie
la matrice A et la meilleure estimation de l’intégrale.
On reprend l’exemple précédent :

def f(x):
return x**5
romberg = Romberg(f,0,1)
(h,I) = romberg.trapezes(1e-7)

print(I[-1]-1.0/6)
--> 2.4835268314093994e-08

print(romberg.n_eval)
--> 4095

romberg = Romberg(f,0,1)
(A,I) = romberg.romberg(1e-7)

print(I-1.0/6)
--> 0.0
Frédéric Legrand Licence Creative Commons 9

print(romberg.n_eval)
--> 7

Avec la méthode des Trapèzes itérative, 4095 évaluations de la fonction sont né-
cessaires pour arriver à une erreur de l’ordre de 10−8 . Avec la méthode de Romberg,
seulement 7 évaluations sont nécessaires pour arriver à une erreur inférieure à la préci-
sion des nombres flottants (il s’agit bien sûr d’un cas très favorable). Voici la matrice A,
qui permet de voir les extrapolations effectuées :

print(A)
--> array([[ 0.5 , 0. , 0. ],
[ 0.265625 , 0.1875 , 0. ],
[ 0.19238281, 0.16796875, 0.16666667]])

Le premier élément de la troisième ligne est obtenu avec la formule des trapèzes pour
N = 4. L’extrapolation à 3 points, par un polynôme de degré 2, permet d’obtenir la
valeur exacte (à la précision des flottants près). Sur cet exemple, l’expression de Ip en
fonction de h2 est probablement un polynome de degré 2.
Voyons un autre exemple plus intéressant, où l’on ne connaît pas d’expression analy-
tique de la primitive :
Z z
2
erf(z) = √ exp(−x2 ) dx (11)
π 0
La fonction scipy.special.erf(z) fournit sa valeur. On fait le test pour z = 1.

from scipy.special import erf

def f(x):
return numpy.exp(-x*x)
romberg = Romberg(f,0,1)
(h,I) = romberg.trapezes(1e-7)

print(I[-1]-numpy.sqrt(numpy.pi)/2*erf(1))
--> -1.4618215860018324e-08

print(romberg.n_eval)
--> 2047

romberg = Romberg(f,0,1)
(A,I) = romberg.romberg(1e-7)

print(I-numpy.sqrt(numpy.pi)/2*erf(1))
--> 2.8266744500626828e-10

print(romberg.n_eval)
--> 15

La méthode de Romberg nécessite seulement 15 évaluations de la fonction (et l’erreur

est plus faible), contre 2047 pour la méthode des trapèzes itérative.
Frédéric Legrand Licence Creative Commons 10

print(A)
--> array([[ 0.68393972, 0. , 0. , 0. ],
[ 0.73137025, 0.74718043, 0. , 0. ],
[ 0.7429841 , 0.74685538, 0.74683371, 0. ],
[ 0.74586561, 0.74682612, 0.74682417, 0.74682402]])
Le résultat est obtenu avec une extrapolation par un polynôme de degré 3.
Considérons l’intégrale suivante :
Z 1√
π
1 − x2 dx = (12)
0 4
from scipy.special import erf
def f(x):
return numpy.sqrt(1-x*x)
romberg = Romberg(f,0,1)
(h,I) = romberg.trapezes(1e-7)

print(I[-1]-numpy.pi/4)
--> -4.9563892989823444e-08
print(romberg.n_eval)
--> 32767
romberg = Romberg(f,0,1)
(A,I) = romberg.romberg(1e-7)

print(I-numpy.pi/4)
--> -0.00018949376813126584
print(romberg.n_eval)
--> 63
Pour cette intégrale, l’erreur finale obtenue avec la méthode de Romberg est beaucoup
plus grande que la tolérance demandée, ce qui signifie que l’extrapolation ne fonctionne
pas bien. On peut essayer de réduire la tolérance :
romberg = Romberg(f,0,1)
(A,I) = romberg.romberg(1e-15)

print(I-numpy.pi/4)
--> -4.6233250450278263e-08
print(romberg.n_eval)
--> 16383
On arrive à la même erreur que la méthode des trapèzes, mais le nombre d’évaluations
n’est que deux fois moindre. Pour cette fonction, la méthode de Romberg est donc peu
efficace.
En principe, la méthode de Romberg fonctionne avec une grande efficacité lorsque la
fonction à intégrer est infiniment dérivable sur l’intervalle d’intégration. Ce n’est pas le
cas du dernier exemple, dont la dérivée n’est pas définie en x = 1.