La méthode de Simpson est une méthode d’intégration qui permet d’obtenir une approcimation d’une intégrale à l’aide de paraboles.
Méthode d’intégration de Simpson: introduction
J’ai déjà rédigé quelques articles sur d’autres méthodes d’approximation d’une intégrale (méthode des rectangles, des trapèzes, du point médian), qui font appel à des polygones.
La méthode de Simpson va plus loin: elle approxime localement une courbe par un morceau de parabole.
Prenons une fonction $f$ définie sur un intervalle $[a;b]$. On subdivise cet intervalle en intervalles $[x_0=a;x_2]$, $[x_2;x_4]$, …, $[x_{2n-2};x_{2n}=b]$. Sur chaque intervalle $[x_{2k};x_{2(k+1)}]$, on choisit alors trois points qui ont pour coordonnées:
- $(x_{2k};f(x_{2k}))$
- $(x_{2k+1};f(x_{2k+1}))$, avec $m=\frac{x_{2k}+x_{2(k+1)}}{2}$
- $(x_{2(k+1)};f(x_{2(k+1)}))$
L’idée consiste alors à construire la parabole passant par ces trois points. L’aire sous cette parabole servira d’approximation de $\displaystyle\int_{x_{2k}}^{x_{2(k+1)}} f(x)\text{d}x$.
Comme on peut le voir sur le schéma, l’intervalle $[a;b]$ est subdivisé en un nombre pair d’intervalles.
Méthode de Simpson: la formule
Polynôme d’interpolation de Lagrange
Avant de regarder ce que donne la formule engendrée par la méthode de Simpson, il faut parler de ce polynôme.
Considérons $n+1$ points de coordonnées $(x_k;f(x_k))$, $0 \leq k \leq n$. Si on souhaite construire un polynôme dont la courbe représentative passe par tous ces points, on utilise la formule de Lagrange:$$L(x)=\sum_{i=0}^n f(x_i)\prod_{\substack{k=0\\k\neq i}}^n \frac{x-x_k}{x_k-x_i}.$$
Appliqué à 3 points $x_k$, $x_{k+1}$ et $x_{k+2}$, cela donne:$$\begin{align*}P_k(x)&=\frac{x-x_{k+1}}{x_k-x_{k+1}}\times \frac{x-x_{k+2}}{x_k-x_{k+2}}f(x_{k})\\& + \frac{x-x_{k}}{x_{k+1}-x_{k}}\times \frac{x-x_{k+2}}{x_{k+1}-x_{k+2}}f(x_{k+1}) \\ & + \frac{x-x_{k}}{x_{k+2}-x_{k}}\times \frac{x-x_{k+1}}{x_{k+2}-x_{k+1}}f(x_{k+2}).\end{align*}$$
N’oublions pas que $x_{k+1}-x_k=x_{k+2}-x_{k+1}=h$, ce qui simplifie l’écriture du polynôme de Lagrange:$$\begin{align*}P_k(x)&=\frac{(x-x_{k}-h)(x-x_k-2h)}{2h^2}f(x_k) + \frac{(x-x_{k})(x-x_k-2h)}{h^2}f(x_{k+1})\\& + \frac{(x-x_{k})(x-x_k-h)}{2h^2}f(x_{k+2}).\end{align*}$$
La formule (composite) de Simpson
On va poser $t=x-x_k$. Alors, le polynôme devient:$$P_k(t)=\frac{(t-h)(t-2h)}{2h^2}f(x_k)-\frac{t(t-2h)}{h^2}f(x_{k+1})+\frac{t(t-h)}{2h^2}f(x_{k+2}).$$ De plus,
- quand $x=x_k$, $t=0$;
- quand $x=x_{k+2}$, $t=x_{k+2}-x_k=2h$;
- $\text{d}t = \text{d}x$.
Ainsi,$$\begin{align*}\int_{x_{k}}^{x_{k+2}}P_k(x)\text{d}x & = \frac{f(x_k)}{2h^2}\int_0^{2h} (t-h)(t-2h) ~\text{d}t – \frac{f(x_{k+1})}{h^2}\int_0^{2h} t(t-2h)~\text{d}t \\& + \frac{f(x_{k+2})}{2h^2}\int_0^{2h} t(t-h)~\text{d}t\\\int_{x_{k}}^{x_{k+2}}P_k(x)\text{d}x&= \frac{f(x_k)}{2h^2}\int_0^{2h} (t^2-3ht+2h^2) ~\text{d}t\\ & – \frac{f(x_{k+1})}{h^2}\int_0^{2h} (t^2-2ht)~\text{d}t \\& + \frac{f(x_{k+2})}{2h^2}\int_0^{2h} (t^2-ht)~\text{d}t\\ \int_{x_{k}}^{x_{k+2}}P_k(x)\text{d}x&=\frac{f(x_k)}{2h^2}\left[ \frac{t^3}{3}-\frac{3ht^2}{2}+2h^2t\right]_0^{2h} \\ & – \frac{f(x_{k+1})}{h^2}\left[ \frac{t^3}{3}-ht^2\right]_0^{2h} \\ & + \frac{f(x_{k+2})}{2h^2}\left[ \frac{t^3}{3} – \frac{ht^2}{2} \right]_0^{2h}\\ \int_{x_{k}}^{x_{k+2}}P_k(x)\text{d}x&= \frac{h}{3}f(x_k) + \frac{4h}{3}f(x_{k+1}) + \frac{h}{3}f(x_{k+2})\\ & = \frac{h}{3}\big( f(x_k) + 4f(x_{k+1}) + f(x_{k+2})\big)\end{align*}$$
À ce stade, on n’oublie pas que $h=\frac{b-a}{2n}$.
On arrive alors à l’approximation:$$\boxed{\int_a^b f(x)~\text{d}x \approx \frac{b-a}{6n}\sum_{k=0}^{n-1}\big( f(x_k) + 4f(x_{k+1}) + f(x_{k+2})\big)}$$
Implémentation en Python
from math import log, exp
def f(x):
return 5 * log(x) * exp(-x/2)
def simpson(a,b,n):
S = 0
h = (b-a)/(2*n)
for k in range(n):
x0 = a + 2*k*h
x1 = x0 + h
x2 = x0 + 2*h
S += f(x0) + 4*f(x1) + f(x2)
return h * S / 3
>>> print( simpson(1,10,10) )
5.4286825969998596
