Mètode de Simpson

De Viquipèdia
Dreceres ràpides: navegació, cerca
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 és un mètode d'integració numèrica, una aproximació numèrica de la integral definida. Específicament, és la següent aproximació:

 \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 a Thomas Simpson.[1]

Obtenció[modifica | modifica el codi]

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

Polinomi d'interpolació de segon grau[modifica | modifica el codi]

Una forma d'obtenir-la es basa a substituir la funció a integrar f(x) per la funció polinòmica de segon grau P(x) que té els mateixos valors que 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,

 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

 \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[modifica | modifica el codi]

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

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

I el mètode trapezial

 T = \tfrac12 (b-a) (f(a)+f(b)).

Els errors d'aquestes aproximacions són

 -\tfrac1{24} (b-a)^3 f''(a) + O((b-a)^4) \quad\text{i}\quad \tfrac1{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 mitja aritmètica ponderada

 \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[modifica | modifica el codi]

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

 \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[modifica | modifica el codi]

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

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

on \xi és algún nombre entre a i b.[3]

L'error és (asimptòticament) proporcional a (b-a)^5. En canvi, els procediments d'obtenció emprats anteriorment, suggerien una error proporcional a (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[modifica | modifica el codi]

Si l'interval d'integració [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 [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 [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 complert. D'aquesta mena d'aproximació sen diu el mètode de Simpson compost.

Suposeu que l'interval [a, b] es parteix en n subintervals, amb n un nombre senar. Llavors, el mètode de Simpson compost ve donat per

\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 x_i=a+ih per i=0, 1, ..., n-1, n amb h=(b-a)/n; en particular, x_0=a i x_n=b. La fórmula de més amunt també es pot escriure com

\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 comes pel mètode de Simpson compost està fitat (en valor absolut) per

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

on h és el "pas", donat per h=(b-a)/n.[4]

Aquesta formulació parteix l'interval [a,b] en subintervals de igual 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[modifica | modifica el codi]

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 dóna una aproximació notablement millor.

Notes[modifica | modifica el codi]

  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[modifica | modifica el codi]

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. 2nd edition. John Wiley & Sons, 1989. ISBN 0-471-50023-2. 
  • Burden, Richard L. and Faires, J. Douglas. Numerical Analysis. 7th edition. 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[modifica | modifica el codi]