5.2 Integração Numérica — Regras de Simpson

1. Introdução

A regra do trapézio, conforme apresentada no post anterior, possui três características que a tornam inapropriada em alguns casos. A primeira propriedade que devemos salientar se dá na sua utilização em relação ao tamanho da aplicação. Isto é, para aplicações individuais em funções bem comportadas, a utilização da regra do trapézio é computacionalmente boa com acurácia suficiente. Contudo, se for necessário uma acurácia muito alta, a regra do trapézio exige um grande esforço computacional, havendo o risco dos cálculos envolvidos estourar a precisão do processador. Ainda que o esforço computacional possa ser desprezível nos computadores atuais, ele pode impactar no desempenho da aplicação quando há um intenso cálculo de integrais ou a própria função a ser integrada é muito complexa.

Os erros de arrendondamento podem limitar nossa habilidade de determinar integrais. Isso decorre tanto da palavra do computador quanto da grande quantidade de cálculos em ponto flutuante, que podem ir acumulando resíduos ao longo da computação. Assim, além de aplicarmos a regra do trapézio com segmentos menores, outra forma de obter uma estimativa mais acurada de uma integral é utilizar polinômios de grau mais alto para ligar os pontos. As fórmulas polinomiais para aproximação dessas integrais são chamadas de regras de Simpson.

 

2. Regra 1/3 de Simpson

A regra 1/3 de Simpson é obtida quando um polinômio interpolador de segundo grau é utilizado. Isto é, \(f_n(x) = f_2(x)\):

\( I = \int_{a}^{b} f(x) dx \approx \int_{a}^{b} f_2(x) dx \)

Se \(a \) e \(b \) forem designados por \(x_0 \) e \(x_2 \) e se \(f_2(x) \) for representado por um polinômio de Lagrange, temos que a integral acima é

\( I = \int_{x_0}^{x_2} \left[ \dfrac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} f(x_0) + \dfrac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} f(x_1) + \dfrac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} f(x_2) \right] ~ dx \)

Essa integral simplificada resulta em
\( I = \dfrac{h}{3} \left[ f(x_0) + 4 f(x_1) + f(x_2) \right] \)

onde \(h=\dfrac{b-a}{2} = \dfrac{x_2-x_0}{2} \) . Essa equação é conhecida como a regra 1/3 de Simpson, sendo a segunda fórmula de Newton-Cotes.

 

2.1. Regra 1/3 em múltiplos intervalos

Assim como fizemos com a regra do trapézio, dividindo o intervalo de integração em diversos subintervalos, fazemos com a regra de Simpson para melhorarmos a estimativa da integral. Assim, ao dividirmos o intervalo \([a.b] \) em \(n \) segmentos, temos

\( h = \dfrac{b-a}{n} \)

Tomando \(b=x_n \) e \(a=x_0 \) , a integral é representada como
\( I = \int_{x_0}^{x_2} f(x) dx + \int_{x_2}^{x_4} f(x) dx + \cdots + \int_{x_{n-2}}^{x_n} f(x) dx \)

onde \(f(x) \) é substituído em cada parte dessa integral por
\( f(x) = \dfrac{(x-x_{i+1})(x-x_{i+2})}{(x_i-x_{i+1})(x_{i}-x_{i+2})} f(x_i) + \dfrac{(x-x_i)(x-x_{i+2})}{(x_{i+1}-x_i)(x_{i+1}-x_{i+2})} f(x_{i+1}) + \dfrac{(x-x_i)(x-x_{i+1})}{(x_{i+2}-x_i)(x_{i+2}-x_{i+1})} f(x_{i+2}) \)

Assim, ao substituirmos a equação acima em cada integral individual, temos que
\( I \approx 2h \dfrac{f(x_0)+ 4 f(x_1)+f(x_2)}{6} + 2h \dfrac{f(x_2) + 4 f(x_3) + f(x_4)}{6} + \cdots 2h \dfrac{f(x_{n-2}) + 4 f(x_{n-1}) + f(x_n)}{6} \)

que pode ser expressa de uma forma mais simplificada:
\( I \approx \dfrac{(b-a)}{3n} \left[ f(x_0) + 4 \sum_{i=1,3,5}^{n-1} f(x_i) + 2 \sum_{j=2,4,6}^{n-2} f(x_j) + f(x_n) \right] \)

 

3. Regra 3/8 de Simpson

De forma análoga à regra 1/3 e à regra do trapézio, a regra 3/8 consiste em tomar um polinômio aproximador. Nesse caso, tomamos a terceira fórmula de Newton-Cotes, \(f_n(x) = f_3(x)\):

\( I \int_{a}^{b} f(x) dx \approx \int_{a}^{b} f_3(x) dx \)

que fornece a seguinte fórmula da regra 3/8 de Simpson
\( I \approx \dfrac{3h}{8} \left[ f(x_0) + 3 f(x_1) + 3 f(x_2) + f(x_3) \right] \)

A figura abaixo demonstra a diferença da utilização dos métodos do trapézio (esquerda), regra de Simpson 1/3 (meio) e Simpson 3/8 (direita) para a função \(f(x) = x^x – x!\) utilizando 2 partições. Como podemos ver, a escolha do método de interpolação tem grande implicação na acurácia da aproximação.



numerical

 

4. Implementação

Regra de Simpson 1/3:

def simpson13(f, xi, xe, n=10000):
    """
    Numerical integration using the Simpson 1/3 Rule.
 
    integral = simpson13(f, xi, xe, n=10000)
 
    INPUT:
      * f: function to be integrated
      * xi: beginning of integration interval
      * xe: end of integration interval
      * n: number of interval divisions
 
    return: \int_{xi}^{xe} f(x)
 
    Author: Pedro Garcia [sawp@sawp.com.br]
    see: http://www.sawp.com.br
 
    License: Creative Commons
             http://creativecommons.org/licenses/by-nc-nd/2.5/br/
 
    Jul 2012
    """
    h = float(xe - xi) / n
    s = lambda i: i % 2 == 0
    fxi = lambda i: f(xi + i * h)
    g = lambda i: s(i) * (2.0 * fxi(i)) + (not s(i)) * (4.0 * fxi(i))
    sumatory = sum(map(g, xrange(1, n)))
    integral = (fxi(0) + sumatory + fxi(n)) * (h / 3.0)
    return integral

Regra de Simpson 3/8:

def simpson38(f, xi, xe, n=10001):
    """
    Numerical integration using the Simpson 3/8 Rule.
 
    integral = simpson38(f, xi, xe, n=10000)
 
    INPUT:
      * f: function to be integrated
      * xi: beginning of integration interval
      * xe: end of integration interval
      * n: number of interval divisions
 
    return: \int_{xi}^{xe} f(x)
 
    Author: Pedro Garcia [sawp@sawp.com.br]
    see: http://www.sawp.com.br
 
    License: Creative Commons
             http://creativecommons.org/licenses/by-nc-nd/2.5/br/
 
    Jul 2012
    """
    h = float(xe - xi) / n
    s = lambda i: i % 3 == 0
    fxi = lambda i: f(xi + i * h)
    g = lambda i: s(i) * (2.0 * fxi(i)) + (not s(i)) * (3.0 * fxi(i))
    sumatory = sum(map(g, xrange(1, n)))
    integral = (fxi(0) + sumatory + fxi(n)) * (3.0 * h / 8.0)
    return integral

Um exemplo de utilização dessas funções pode ser obtido em http://www.sawp.com.br/code/integral/integrals_simpson.py.

 

5. Copyright

Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.

References

[1] Anthony Ralston and Philip Rabinowitz, A First Course in Numerical Analysis (2nd ed.), McGraw-Hill and Dover, (2001).

[2] N.B Franco, Cálculo Numérico, Pearson Prentice Hall (2006).