1Rectangles, trapezes, Simpson et leurs erreurs
Calculer \(\displaystyle I = \int_a^b f(t)\,\mathrm dt\) lorsque aucune primitive de \(f\) n'est exprimable a l'aide des fonctions usuelles, ou lorsque \(f\) n'est connue que par des valeurs tabulees.
Exemples classiques sans primitive elementaire : \(\displaystyle\int_0^1 e^{-t^2}\mathrm dt\) (loi normale), \(\displaystyle\int_0^1 \frac{\sin t}{t}\mathrm dt\) (sinus integral), \(\displaystyle\int_0^{\pi/2}\frac{\mathrm dt}{\sqrt{1-k^2\sin^2 t}}\) (integrale elliptique, periode du pendule).
- Subdiviser \([a,b]\) en \(n\) sous-intervalles de longueur \(h = \dfrac{b-a}{n}\), de bornes \(x_k = a + kh\).
- Sur chaque sous-intervalle, remplacer \(f\) par un polynome interpolateur de faible degre, dont on sait integrer exactement.
- Sommer les contributions.
Le degre du polynome determine l'ordre de la methode : constante (rectangles, ordre 1), affine (trapezes, ordre 2), quadratique (Simpson, ordre 4).
Methode des rectangles
On approche \(f\) par une constante sur chaque \([x_k, x_{k+1}]\) :
\[ \int_{x_k}^{x_{k+1}} f \approx h\,f(c_k), \qquad c_k \in \{x_k,\ x_{k+1},\ \tfrac{x_k+x_{k+1}}{2}\}. \]
def rectangles_gauche(f, a, b, n):
"""Approximation par la valeur a gauche. Ordre 1."""
h = (b - a) / n
return h * sum(f(a + k*h) for k in range(n))
def rectangles_droite(f, a, b, n):
h = (b - a) / n
return h * sum(f(a + k*h) for k in range(1, n+1))
def point_milieu(f, a, b, n):
"""Valeur au milieu de chaque sous-intervalle. Ordre 2, gratuit !"""
h = (b - a) / n
return h * sum(f(a + (k + 0.5)*h) for k in range(n))Sur un seul sous-intervalle, la formule de Taylor donne pour le rectangle a gauche \[ \int_{x_k}^{x_k+h} f = h f(x_k) + \frac{h^2}{2}f'(x_k) + O(h^3), \] donc une erreur locale en \(O(h^2)\). En sommant sur \(n = (b-a)/h\) intervalles, l'erreur globale est en \(O(h)\) : la methode est d'ordre 1. Precisement, \[ \left|I - R_n\right| \le \frac{(b-a)h}{2}\max_{[a,b]}|f'|. \]
Le point milieu, lui, est d'ordre 2 : par symetrie, le terme en \(f'\) s'annule, et \[ \left|I - M_n\right| \le \frac{(b-a)h^2}{24}\max|f''|. \] C'est remarquable : au meme cout (une evaluation par intervalle), on gagne un ordre.
Methode des trapezes
On approche \(f\) par le segment joignant \((x_k, f(x_k))\) a \((x_{k+1}, f(x_{k+1}))\) :
\[ \int_{x_k}^{x_{k+1}} f \approx h\,\frac{f(x_k)+f(x_{k+1})}{2}. \]
En sommant et en regroupant, chaque point interieur apparait deux fois :
\[ \boxed{\;T_n = h\left(\frac{f(x_0)}{2} + \sum_{k=1}^{n-1} f(x_k) + \frac{f(x_n)}{2}\right).\;} \]
def trapezes(f, a, b, n):
"""Methode des trapezes. Ordre 2. n+1 evaluations de f."""
h = (b - a) / n
s = (f(a) + f(b)) / 2
for k in range(1, n):
s += f(a + k*h)
return h * sSi \(f\) est de classe \(\mathcal C^2\) sur \([a,b]\), il existe \(\xi \in\;]a,b[\) tel que \[ I - T_n = -\frac{(b-a)h^2}{12}f''(\xi), \qquad\text{donc}\qquad |I - T_n| \le \frac{(b-a)h^2}{12}\max_{[a,b]}|f''|. \] La methode est d'ordre 2 : diviser \(h\) par 2 divise l'erreur par 4.
Sur \([x_k, x_{k+1}]\), on pose \(e(h) = \int_{x_k}^{x_k+h} f - \frac h2(f(x_k)+f(x_k+h))\). On a \(e(0) = 0\), et en derivant deux fois par rapport a \(h\) : \(e''(h) = -\frac h2 f''(x_k+h)\). Une double integration et le theoreme des accroissements finis donnent \(e(h) = -\frac{h^3}{12}f''(\xi_k)\). En sommant sur les \(n\) intervalles et en utilisant le theoreme des valeurs intermediaires pour \(f''\), on obtient la formule globale.
Signe de l'erreur : si \(f\) est convexe (\(f'' \ge 0\)), les cordes sont au-dessus de la courbe, donc \(T_n \ge I\). Les trapezes surestiment une integrale de fonction convexe, et le point milieu la sous-estime.
Methode de Simpson
On approche \(f\) par la parabole passant par \((x_k, f(x_k))\), \((x_{k+1/2}, f(x_{k+1/2}))\) et \((x_{k+1}, f(x_{k+1}))\) :
\[ \int_{x_k}^{x_{k+1}} f \approx \frac{h}{6}\Bigl(f(x_k) + 4f(x_{k+1/2}) + f(x_{k+1})\Bigr). \]
def simpson(f, a, b, n):
"""Methode de Simpson composee. n doit etre PAIR. Ordre 4."""
if n % 2 != 0:
n += 1
h = (b - a) / n
s = f(a) + f(b)
s += 4 * sum(f(a + k*h) for k in range(1, n, 2)) # indices impairs
s += 2 * sum(f(a + k*h) for k in range(2, n, 2)) # indices pairs
return h * s / 3\[ I - S_n = -\frac{(b-a)h^4}{180}f^{(4)}(\xi), \qquad |I - S_n| \le \frac{(b-a)h^4}{180}\max|f^{(4)}|. \] La methode est d'ordre 4 : diviser \(h\) par 2 divise l'erreur par 16.
Fait remarquable : bien que construite sur des paraboles, la methode de Simpson est exacte pour les polynomes de degre 3. Le terme d'erreur en \(f^{(3)}\) s'annule par symetrie ; c'est un degre gagne gratuitement.
assert abs(simpson(lambda x: x**3, 0, 1, 2) - 0.25) < 1e-15 # EXACT
assert abs(simpson(lambda x: x**4, 0, 1, 2) - 0.2) > 1e-6 # plus exactSimpson est la combinaison \(S = \dfrac{T + 2M}{3}\), ou \(T\) est la formule des trapezes et \(M\) celle du point milieu sur le meme decoupage. C'est aussi l'extrapolation de Richardson de la formule des trapezes : \[ S_{2n} = \frac{4T_{2n} - T_n}{3}. \] On elimine ainsi le terme d'erreur en \(h^2\), gagnant deux ordres d'un coup (le developpement d'Euler-Maclaurin des trapezes ne contient que des puissances paires).
2Visualisation interactive et comparaison
Comparaison des methodes
| Methode | Formule locale | Ordre | Erreur | Exacte pour |
|---|---|---|---|---|
| Rectangles (gauche/droite) | \(h f(x_k)\) | 1 | \(\frac{(b-a)h}{2}\max|f'|\) | degre 0 |
| Point milieu | \(h f(x_{k+1/2})\) | 2 | \(\frac{(b-a)h^2}{24}\max|f''|\) | degre 1 |
| Trapezes | \(\frac h2(f_k+f_{k+1})\) | 2 | \(\frac{(b-a)h^2}{12}\max|f''|\) | degre 1 |
| Simpson | \(\frac h6(f_k+4f_{k+1/2}+f_{k+1})\) | 4 | \(\frac{(b-a)h^4}{180}\max|f^{(4)}|\) | degre 3 |
- Le point milieu est deux fois plus precis que les trapezes (constante 1/24 contre 1/12) pour le meme nombre d'evaluations. Il est pourtant moins enseigne.
- Les erreurs du point milieu et des trapezes sont de signes opposes : leur combinaison \(\frac{T + 2M}{3}\) elimine le terme dominant et donne Simpson.
- Une methode a \(m\) points est generalement exacte jusqu'au degre \(m-1\) ; les formules symetriques gagnent un degre supplementaire.
from math import exp, e
I = e - 1
print(f"{'n':>6}{'rectangles':>16}{'milieu':>16}{'trapezes':>16}{'Simpson':>16}")
for n in (2, 4, 8, 16, 32, 64):
r = abs(rectangles_gauche(exp, 0, 1, n) - I)
m = abs(point_milieu(exp, 0, 1, n) - I)
t = abs(trapezes(exp, 0, 1, n) - I)
s = abs(simpson(exp, 0, 1, n) - I)
print(f"{n:>6}{r:>16.3e}{m:>16.3e}{t:>16.3e}{s:>16.3e}") n rectangles milieu trapezes Simpson
2 4.259e-01 2.146e-02 4.312e-02 5.792e-05
4 2.281e-01 5.363e-03 1.076e-02 3.587e-06
8 1.181e-01 1.341e-03 2.688e-03 2.232e-07
16 6.008e-02 3.352e-04 6.719e-04 1.393e-08
32 3.030e-02 8.380e-05 1.680e-04 8.703e-10
64 1.522e-02 2.095e-05 4.199e-05 5.439e-11Les rapports d'une ligne a la suivante valent 2, 4, 4 et 16 : ils confirment les ordres 1, 2, 2 et 4. Avec seulement 64 intervalles, Simpson atteint \(5\cdot10^{-11}\) la ou les rectangles plafonnent a \(1{,}5\cdot10^{-2}\) : neuf ordres de grandeur d'ecart.
3Methodes adaptatives, integrales impropres et applications
Methode adaptative
Un pas uniforme est inefficace si \(f\) varie beaucoup sur une petite portion de l'intervalle et peu ailleurs. Une methode adaptative subdivise recursivement, et ne raffine que la ou l'erreur estimee est trop grande.
L'estimation d'erreur repose sur la comparaison entre une approximation grossiere et une approximation fine : pour Simpson d'ordre 4, \[ |I - S_{\text{fine}}| \approx \frac{|S_{\text{fine}} - S_{\text{grossiere}}|}{15}. \]
def simpson_simple(f, a, b):
"""Simpson sur un seul intervalle : 3 evaluations."""
m = (a + b) / 2
return (b - a) / 6 * (f(a) + 4*f(m) + f(b))
def adaptative(f, a, b, eps=1e-10, profondeur=50):
"""Simpson adaptative recursive. Renvoie (valeur, nombre d'evaluations)."""
compteur = [0]
def g(x):
compteur[0] += 1
return f(x)
def aux(a, b, eps, entier, prof):
m = (a + b) / 2
gauche = simpson_simple(g, a, m)
droite = simpson_simple(g, m, b)
if prof <= 0 or abs(gauche + droite - entier) <= 15 * eps:
# extrapolation de Richardson : on ajoute l'estimation d'erreur
return gauche + droite + (gauche + droite - entier) / 15
return (aux(a, m, eps/2, gauche, prof-1)
+ aux(m, b, eps/2, droite, prof-1))
total = simpson_simple(g, a, b)
return aux(a, b, eps, total, profondeur), compteur[0]
from math import sin, exp, pi
# fonction reguliere : peu d'evaluations
print(adaptative(sin, 0, pi)) # (2.0000000000, ~ 40 evaluations)
# fonction a variation localisee : le raffinement se concentre au bon endroit
f = lambda x: exp(-100 * (x - 0.3)**2) # gaussienne tres piquee
v, ev = adaptative(f, 0, 1, 1e-12)
print(f"valeur {v:.12f}, {ev} evaluations")
# comparaison a pas uniforme
for n in (100, 1000, 10000):
print(f" Simpson uniforme n = {n:>6} : {simpson(f, 0, 1, n):.12f}")Sur une gaussienne piquee, la methode adaptative atteint \(10^{-12}\) avec quelques centaines d'evaluations, la ou une methode uniforme en demanderait des dizaines de milliers : elle place les points la ou ils sont utiles.
Integrales impropres
from math import sqrt, log
# 1. singularite integrable : 1/sqrt(x) sur ]0, 1]
trapezes(lambda x: 1/sqrt(x), 0, 1, 1000) # ZeroDivisionError en x = 0
# 2. borne infinie
trapezes(lambda x: exp(-x*x), 0, float("inf"), 1000) # impossible1. Point milieu : il n'evalue jamais aux bornes, donc tolere une singularite integrable aux extremites.
print(point_milieu(lambda x: 1/sqrt(x), 0, 1, 10**6)) # 1.999... (exact : 2)La convergence est toutefois degradee (ordre \(1/2\) au lieu de 2) car \(f''\) n'est pas bornee.
2. Changement de variable pour absorber la singularite. Pour \(\int_0^1 \frac{\mathrm dx}{\sqrt x}\), on pose \(x = t^2\) :
# int_0^1 dx/sqrt(x) = int_0^1 2 dt : la singularite a DISPARU
print(trapezes(lambda t: 2.0, 0, 1, 10)) # 2.0 exactement3. Troncature ou changement de variable pour une borne infinie. Pour \(\int_0^{+\infty} e^{-x^2}\mathrm dx = \frac{\sqrt\pi}{2}\) :
from math import sqrt, pi, exp
# a) troncature : l'integrande est negligeable au-dela de 6
print(simpson(lambda x: exp(-x*x), 0, 6, 1000), sqrt(pi)/2)
# b) changement de variable x = t/(1-t), qui envoie [0,1[ sur [0, +inf[
def g(t):
if t >= 1: return 0.0
x = t / (1 - t)
return exp(-x*x) / (1 - t)**2 # facteur dx/dt
print(point_milieu(g, 0, 1, 10000), sqrt(pi)/2)La troncature exige de majorer la queue : ici \(\int_6^{+\infty} e^{-x^2}\mathrm dx < \frac{e^{-36}}{12} \approx 10^{-17}\), donc negligeable devant la precision machine.
Application — fonctions speciales
from math import exp, sqrt, pi, sin
def erf(x, n=2000):
"""Fonction d'erreur : erf(x) = 2/sqrt(pi) * int_0^x exp(-t^2) dt."""
return 2/sqrt(pi) * simpson(lambda t: exp(-t*t), 0, x, n)
def Si(x, n=2000):
"""Sinus integral : Si(x) = int_0^x sin(t)/t dt. (prolonge par 1 en 0)"""
f = lambda t: 1.0 if t == 0 else sin(t)/t
return simpson(f, 0, x, n)
def periode_pendule(theta0, L=1.0, g=9.81, n=2000):
"""Periode exacte d'un pendule simple d'amplitude theta0 (integrale elliptique)."""
k = sin(theta0 / 2)
integrale = simpson(lambda t: 1/sqrt(1 - k*k*sin(t)**2), 0, pi/2, n)
return 4 * sqrt(L/g) * integrale
import math
assert abs(erf(1) - math.erf(1)) < 1e-12
print(f"Si(1) = {Si(1):.10f}") # 0.9460830704
print(f"petites oscillations : {2*pi*sqrt(1/9.81):.6f} s")
for a in (10, 30, 60, 90, 150):
T = periode_pendule(math.radians(a))
print(f" amplitude {a:>3} deg : T = {T:.6f} s "
f"(+{100*(T/(2*pi*sqrt(1/9.81)) - 1):.2f} %)")petites oscillations : 2.006066 s
amplitude 10 deg : T = 2.007888 s (+0.09 %)
amplitude 30 deg : T = 2.023420 s (+0.86 %)
amplitude 60 deg : T = 2.077848 s (+3.58 %)
amplitude 90 deg : T = 2.178146 s (+8.58 %)
amplitude 150 deg : T = 2.751594 s (+37.16 %)L'approximation des petites oscillations \(T \approx 2\pi\sqrt{L/g}\), universellement enseignee, n'est valable qu'a 1 % pres jusqu'a 30 degres. A 150 degres, l'erreur atteint 37 % : le calcul numerique donne acces au regime non lineaire, inaccessible analytiquement.
L'essentiel du chapitre
- Rectangles : ordre 1. Point milieu et trapezes : ordre 2. Simpson : ordre 4.
- Doubler \(n\) divise l'erreur par \(2^p\) ou \(p\) est l'ordre : c'est le test de validation.
- Formule des trapezes : \(T_n = h\left(\frac{f_0}{2} + \sum_{k=1}^{n-1}f_k + \frac{f_n}{2}\right)\).
- Simpson est exacte jusqu'au degre 3, et vaut \(\frac{4T_{2n}-T_n}{3}\).
- Trapezes surestiment et point milieu sous-estime une integrale de fonction convexe.
- Le point milieu tolere une singularite aux bornes ; sinon, changer de variable.
- Les methodes adaptatives concentrent les evaluations la ou \(f\) varie.