Integración por la fórmula de Simpson

Vamos a calcular el valor de una integral definida

$$\int_a ^ b f (x) dx$$

La solución descrita aquí se publicó en una de las disertaciones de Thomas Simpson en 1743.

Fórmula de Simpson

Sea $n$ algún número natural. Dividimos el segmento de integración $[a, b]$ en $2n$ partes iguales:

$$x_i = a + i h, ~~ i = 0 \ldots 2n,$$

$$h = \frac {b-a} {2n}.$$

Ahora calculamos la integral por separado en cada uno de los segmentos $[x_ {2i-2}, x_ {2i}]$, $i = 1 \ldots n$, y luego sumamos todos los valores.

Así, supongamos que consideramos el siguiente segmento $[x_ {2i-2}, x_ {2i}], i = 1 \ldots n$. Reemplazamos la función $f(x)$ en él por una parábola $P(x)$ que pasa por 3 puntos $(x_ {2i-2}, x_ {2i-1}, x_ {2i})$. Tal parábola siempre existe y es única; se puede hallar analíticamente. Por ejemplo podríamos construirla usando la interpolación polinómica de Lagrange. Lo único que queda por hacer es integrar este polinomio. Si se hace esto para una función general $f$, se recibe una expresión notablemente simple:

$$\int_{x_ {2i-2}} ^ {x_ {2i}} f (x) ~dx \approx \int_{x_ {2i-2}} ^ {x_ {2i}} P (x) ~dx = \left(f(x_{2i-2}) + 4f(x_{2i-1})+(f(x_{2i})\right)\frac {h} {3}$$

Sumando estos valores sobre todos los segmentos, obtenemos la fórmula de Simpson final:

$$\int_a ^ b f (x) dx \approx \left(f (x_0) + 4 f (x_1) + 2 f (x_2) + 4f(x_3) + 2 f(x_4) + \ldots + 4 f(x_{2N-1}) + f(x_{2N}) \right)\frac {h} {3}$$

Error

El error al aproximar una integral por la fórmula de Simpson es

$$-\tfrac{1}{90} \left(\tfrac{b-a}{2}\right)^5 f^{(4)}(\xi)$$

donde $\xi$ es algún número entre $a$ y $b$.

El error es asintóticamente proporcional a $(b-a)^5$. Sin embargo, las derivaciones de arriba sugieren un error proporcional a $(b-a)^4$. La regla de Simpson gana un orden extra porque los puntos en los que se evalúa el integrando están distribuidos simétricamente en el intervalo $[a, b]$.

Implementación

Aquí, $f(x)$ es alguna función definida por el usuario.

const int N = 1000 * 1000; // number of steps (already multiplied by 2)
 
double simpson_integration(double a, double b){
    double h = (b - a) / N;
    double s = f(a) + f(b); // a = x_0 and b = x_2n
    for (int i = 1; i <= N - 1; ++i) { // Refer to final Simpson's formula
        double x = a + h * i;
        s += f(x) * ((i & 1) ? 4 : 2);
    }
    s *= h / 3;
    return s;
}

Problemas de práctica

• Latin American Regionals 2012 - Environment Protection