Integral pela fórmula de Simpson

Iremos calcular o valor de uma integral definida

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

A solução descrita aqui foi publicada em uma das dissertações de Thomas Simpson em 1743.

Fórmula de Simpson

Seja $n$ um número natural. Dividimos o segmento de integração $[a, b]$ em $2n$ partes iguais:

$$x_i = a + i h, ~~ i = 0 \ldots 2n,$$ $$h = \frac {b-a} {2n}.$$

Agora calculamos a integral separadamente em cada um dos segmentos $[x_ {2i-2}, x_ {2i}]$, $i = 1 \ldots n$, e, em seguida, adicionamos todos os valores.

Portanto, suponha que consideremos o próximo segmento $[x_ {2i-2}, x_ {2i}], i = 1 \ldots n$. Substitua a função $f(x)$ por uma parábola $P(x)$ passando por 3 pontos $(x_ {2i-2}, x_ {2i-1}, x_ {2i})$. Essa parábola sempre existe e é única; pode ser encontrado analiticamente. Por exemplo, podemos construí-lo usando a interpolação polinomial de Lagrange. A única coisa que resta a fazer é integrar esse polinômio. Se você fizer isso para uma função geral $f$, você receberá essa expressão:

$$\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} $$

Adicionando esses valores em todos os segmentos, obtemos a fórmula final de Simpson:

$$\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} $$

Erro

O erro na aproximação de uma integral pela fórmula de Simpson é

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

onde $\xi$ é um número entre $a$ e $b$.

O erro é assintoticamente proporcional a $(b-a)^5$. No entanto, as derivações acima sugerem um erro proporcional a $(b-a)^4$. A regra de Simpson ganha uma ordem extra porque os pontos nos quais o integrando é avaliado são distribuídos simetricamente no intervalo $[a, b]$.

Implementação

Aqui, $f(x)$ é uma função definida pelo usuário.

const int N = 1000 * 1000; // número de etapas 

double simpson_integration(double a, double b){
    double h = (b - a) / N;
    double s = f(a) + f(b); // a = x_0     b = x_2n
    for (int i = 1; i <= N - 1; ++i) { // refere-se a fórmula final de Simpson
        double x = a + h * i;
        s += f(x) * ((i & 1) ? 4 : 2);
    }
    s *= h / 3;
    return s;
}

Problemas