Mètode de Simpson

El mètode de Simpson es pot obtenir aproximant l'integrand f(x) (en blau) per un polinomi d'interpolació de segon grau P(x) (en vermell).

En càlcul numèric, el mètode de Simpson o regla de Simpson és un mètode d'integració numèrica que dona una aproximació numèrica de la integral definida. Específicament, és la següent aproximació:

${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right].}$

Rep el seu nom en honor de Thomas Simpson.^[1]

Obtenció

El mètode de Simpson es pot obtenir de diverses maneres.

Polinomi d'interpolació de segon grau

Una forma d'obtenir-la es basa a substituir la funció a integrar ${\displaystyle f(x)}$ per la funció polinòmica de segon grau ${\displaystyle P(x)}$ que té els mateixos valors que ${\displaystyle f(x)}$ als extrems a i b i al mig m = (a+b) / 2. Es pot fer servir la interpolació polinòmica de Lagrange per a trobar una expressió d'aquest polinomi,

${\displaystyle P(x)=f(a){\frac {(x-m)(x-b)}{(a-m)(a-b)}}+f(m){\frac {(x-a)(x-b)}{(m-a)(m-b)}}+f(b){\frac {(x-a)(x-m)}{(b-a)(b-m)}}.}$

Un càlcul fàcil (però pesat) mostra que

${\displaystyle \int _{a}^{b}P(x)\,dx={\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right].}$^[2]

Mitjana aritmètica ponderada dels mètodes trapezial i rectangular

Una altra forma d'obtenir el mètode de Simpson és calculant la mitjana aritmètica ponderada de dues aproximacions més simples: el mètode rectangular

${\displaystyle M=(b-a)f\left({\frac {a+b}{2}}\right)}$

I el mètode trapezial

${\displaystyle T={\tfrac {1}{2}}(b-a)(f(a)+f(b)).}$

Els errors d'aquestes aproximacions són

${\displaystyle -{\tfrac {1}{24}}(b-a)^{3}f''(a)+O((b-a)^{4})\quad {\text{i}}\quad {\tfrac {1}{12}}(b-a)^{3}f''(a)+O((b-a)^{4}),}$

respectivament. D'aquí en resulta que el terme principal de l'error s'anul·la si es pren la mitjana aritmètica ponderada

${\displaystyle {\frac {2M+T}{3}}.}$

Aquesta mitjana aritmètica ponderada és exactament el mètode de Simpson.

Emprant una altra aproximació (per exemple, el mètode trapezial amb doble quantitat de punts), és possible d'obtenir una mitjana aritmètica ponderada adequada per eliminar un altre terme de l'error. Aquest és el mètode de Romberg.

Coeficients indeterminats

El tercer mètode comença a partir de la hipòtesi que

${\displaystyle \int _{a}^{b}f(x)\,dx\approx \alpha f(a)+\beta f\left({\frac {a+b}{2}}\right)+\gamma f(b).}$

Els coeficients α, β i γ es poden determinar imposant que aquesta aproximació sigui exacta pel cas de polinomis de segon grau. Això porta al mètode de Simpson.

Error

L'error en aproximar una integral amb el mètode de Simpson és

${\displaystyle -{\frac {(b-a)^{5}}{2880}}f^{(4)}(\xi ),}$

on ${\displaystyle \xi }$ és algun nombre entre ${\displaystyle a}$ i ${\displaystyle b}$.^[3]

L'error és (asimptòticament) proporcional a ${\displaystyle (b-a)^{5}}$. En canvi, els procediments d'obtenció emprats anteriorment, suggerien una error proporcional a ${\displaystyle (b-a)^{4}}$. El mètode de Simpson guanya un ordre extra perquè els punts en què s'avalua l'integrand, estan distribuïts simètricament en l'interval [a, b].

Mètode de Simpson compost

Si l'interval d'integració ${\displaystyle [a,b]}$ és "petit" en algun sentit, llavors el mètode de Simpson donarà una aproximació adequada de la integral exacta. Petit, el que en realitat significa és que la funció a integrar sigui relativament suau sobre l'interval ${\displaystyle [a,b]}$. Per aquest tipus de funcions, una interpolació amb un polinomi de segon grau suau, com el que es fa servir en el mètode de Simpson, donarà bons resultats.

En canvi, sovint el cas és que la funció que s'està intentant integrar no és suau sobre l'interval. Típicament, això significa que la funció és altament oscil·lant o que no té derivades a certs punts. En aquests casos, el mètode de Simpson pot donar resultats molt pobres. Una forma habitual de manejar aquest problema és a base de trencar l'interval ${\displaystyle [a,b]}$ en un determinat nombre de subintervals petits. Llavors s'aplica el mèode de Simpson a cada subinterval, i se sumen els resultats per a obtenir una aproximació de la integral sobre l'inteval complet. D'aquesta mena d'aproximació sen diu el mètode de Simpson compost.

Suposeu que l'interval ${\displaystyle [a,b]}$ es parteix en ${\displaystyle n}$ subintervals, amb ${\displaystyle n}$ un nombre parell. Llavors, el mètode de Simpson compost dona l'aproximació

${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {h}{3}}{\bigg [}f(x_{0})+2\sum _{j=1}^{n/2-1}f(x_{2j})+4\sum _{j=1}^{n/2}f(x_{2j-1})+f(x_{n}){\bigg ]},}$

on ${\displaystyle x_{i}=a+ih}$ per ${\displaystyle i=0,1,...,n-1,n}$ amb ${\displaystyle h=(b-a)/n}$; en particular, ${\displaystyle x_{0}=a}$ i ${\displaystyle x_{n}=b}$. La fórmula anterior també es pot escriure com

${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {h}{3}}{\bigg [}f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+2f(x_{4})+\cdots +4f(x_{n-1})+f(x_{n}){\bigg ]}.}$

L'error comès pel mètode de Simpson compost està fitat (en valor absolut) per

${\displaystyle {\frac {h^{4}}{180}}(b-a)\max _{\xi \in [a,b]}|f^{(4)}(\xi )|,}$

on ${\displaystyle h}$ és el pas, donat per ${\displaystyle h=(b-a)/n.}$^[4]

Aquesta formulació parteix l'interval ${\displaystyle [a,b]}$ en subintervals de la mateixa longitud. A la pràctica, sovint és avantatjós emprar subintervals de diferents longituds, i concentrar els esforços en els llocs on l'integrand té menys bon comportament. Això porta al mètode de Simpson adaptatiu.

Implementació en Python del mètode de Simpson

Aquí hi ha una implementació del mètode de Simpson en Python.

def metode_de_simpson(f, a, b):
 "Aproxima la integral definida de f des de a fins a b pel mètode de Simpson."
 c = (a + b) / 2.0
 h3 = abs(b - a) / 6.0
 return h3 * (f(a) + 4.0*f(c) + f(b))

# Càlcul de la integral de sin(x) des de 0 fins a 1
from math import sin
print metode_de_simpson(sin, 0, 1)

Aquesta és la versió del mètode de Simpson compost, també en Python.

def metode_de_simpson_compost(f,a,b,n):
 "Aproxima la integral definida de f des de a fins a b amb el mètode de Simpson
 compost, dividint l'interval en n parts."
 dx = (float(b)-a)/n
 i = 0
 fks = 0.0
 while i<=n:
 xk = a+i*dx
 if i==0 or i==n:
 fk = f(xk)
 elif i%2 == 1:
 fk = 4*f(xk)
 else:
 fk = 2*f(xk)
 fks += fk
 i += 1
 return (dx/3)*fks

# Càlcul de la integral de sin(x) des de 0 fins a 1
from math import sin
print metode_de_simpson_compost(sin, 0, 1, 10000)

Integrant sin(x) des de 0 fins a 1 amb el primer codi dona 0.459862189871, el segon dona 0.459697694132 mentre que el verdader valor és 1 − cos(1) = 0.45969769413... La versió composta del mètode dona una aproximació notablement millor.

Notes

1. Süli and Mayers, §7.2
2. Atkinson, p. 256; Süli and Mayers, §7.2
3. Atkinson, equation (5.1.15); Süli and Mayers, Theorem 7.2
4. Atkinson, pp. 257+258; Süli and Mayers, §7.5

Referències

El mètode de Simpson s'esmenta en molts llibre de text de càlcul numèric:

• Atkinson, Kendall A.. An Introduction to Numerical Analysis. 2a ed.. John Wiley & Sons, 1989. ISBN 0-471-50023-2. 
• Burden, Richard L. and Faires, J. Douglas. Numerical Analysis. 7a ed.. Brooks/Cole, 2000. ISBN 0-534-38216-9. 
• Süli, Endre and Mayers, David. An Introduction to Numerical Analysis. Cambridge University Press, 2003. ISBN 0-521-81026-4 (hardback), ISBN 0-521-00794-1 (paperback). 

Enllaços externs

• Simpson's Rule for Numerical Integration
• Application of Simpson's Rule - Earthwork Excavation Arxivat 2007-10-18 a Wayback Machine.
• Simpson's 1/3rd rule of integration - Notes, PPT, Mathcad, Matlab, Mathematica, Maple