5.4 Integração Numérica — Integração de Romberg

1. Introdução

A integração de Romberg é um método proposto em 1955, e é utilizado para calcular integrais numéricas com alta velocidade de convergência. Esta técnica é semelhante à regra dos trapézios, uma vez que utiliza esse método sucessivas vezes para calcular o valor da integral. Contudo, difere-se da regra dos trapézios por aplicar a extrapolação de Richardson para acelerar a convergência da aproximação. O método de Romberg é uma forma geral de calcular as fórmulas de Newton-Cotes, uma vez que avalia os integrandos em pontos igualmente espaçados, onde os integrandos precisam ter derivadas contínuas.

 

2. Extrapolação de Richardson

Assim como apresentado anteriormente no post sobre a extrapolação de Richardson, onde mostramos como acelerar a convergência de derivadas numéricas, da mesma forma podemos utilizar esse método para aumentar a acurácia das funções de integração numérica. Desta maneira, a estimativa da integral \(I \) e o erro associado \(E \) com a aplicação múltipla da regra do trapézio podem ser representados de um modo geral por

\( I = I(h) + E(h) \)

onde \(I \) é o valor exato da integral, \(I(h) \) é a aproximação a partir da aplicação da regra do trapézio (em função do passo \(h \) ) e \(E(h) \) é o erro de truncamento. Se fizermos duas estimativas independentes usando passos distintos \(h_1 \) e \(h_2 \) , temos que

\( I(h_1) + E(h_1) = I(h_2) + E(h_2) \)

Como o erro da aplicação da regra do trapézio é
\( E \approx – \dfrac{b-a}{12} h^2 \bar{f}” \)

onde \(\bar{f}” \) é a derivada média da função nos pontos avaliados. Supondo que \(\bar{f}” \) é constante, independente do tamanho do passo, a equação acima pode ser usada para determinar a razão entre os erros:

\( \dfrac{E(h_1)}{E(h_2)} \approx \dfrac{h_1^2}{h_2^2} \)

ou melhor
\( E(h_1) \approx E(h_2) \left( \dfrac{h_1}{h_2} \right)^2 \)

que pode ser substituida na equação das integrais:
\( I(h_1) + E(h_2) \left( \dfrac{h_1}{h_2} \right)^2 \approx I(h_2) + E(h_2) \)

isolando \(E(h_2) \) , temos
\( E(h_2) \approx \dfrac{I(h_1) – I(h_2)}{1 – (h_1 / h_2)^2} \)

Logo, obteremos uma estimativa do erro de truncamento em termos das estimativas da integral e seus tamanhos de passo. Essa estimativa pode então ser substituída em
\( I = I(h_2) + E(h_2) \)

para fornecer uma estimativa mais acurada
\( I \approx I(h_2) + \dfrac{1}{(h_1/h_2)^2 – 1} \left[ I(h_2) – I(h_1) \right] \)

Para o caso em que \(h_2 = h_1/2 \) , temos que
\( I \approx \dfrac{4}{3} I(h_2) – \dfrac{1}{3} I(h_1) \)

 

3. O Método de Romberg

Observe que na última equação acima, os coeficientes somam 1. Isso equivale à fatores de peso de uma média ponderada. Assim, conforme a precisão aumenta, podemos colocar peso relativamente maior na estimativa superior da integral. Essas ponderações são expressas pela seguinte equação:

\( I_{j, k} = \dfrac{4^{k-1} I_{j+1, k-1} – I_{j, k-1}}{4^{k-1}-1} \)

onde \(I_{j+1,k-1} \) e \(I_{j, k-1} \) são as integrais mais e menos acuradas, respectivamente, e \(I_{j, k} \) é a integral melhorada. O índice \(k\) representa a respectiva fórmula de Newton-Cotes, isto é, para \(k=1 \) temos a regra do trapézio, para \(k=2 \) temos a regra de Simpson, e assim por diante. O índice \(j \) serve para indicar a estimativa mais acurada da menos acurada (quanto maior \(j \) , mais acurada). Assim, o método pode ser definido indutivamente da seguinte forma

  • \(R(0,0) = \dfrac{1}{2} \left( b-a \right) \left( f(a) + f(b) \right)\)
  • \(R(n,0) = \dfrac{1}{2} R(n-1, 0) + h_n \sum_{k=1}^{2^{n-1}} f(a + h_n (2k-1))\)
  • \(R(n,m) = \dfrac{1}{4^m – 1} \left( 4^m R(n,m-1) – R(n-1, m-1) \right)\)

onde \(a \) e \(b \) são os limites do intervalo de integração, \(n \geq 1 \) , $m \geq 1 \( e \) h_n = \dfrac{b-1}{2^n}$.

 

4. Implementação

A implementação utilizando as fórmulas indutivas é

def romberg_slow(f, xi, xe, ni=20, mi=2):
    """
    Numerical integration using Romberg's Integration.
 
    integral = romberg_slow(f, xi, xe, ni=20, mi=2)
 
    INPUT:
      * f: function to be integrated
      * xi: beginning of integration interval
      * xe: end of integration interval
      * ni: initial partitions n
      * mi: initial order m
 
    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
    """
    def R00():
        return 0.5 * (xe - xi) * (f(xi) + f(xe))
 
    def Rn0(n):
        hn = (xe - xi) / (2.0 ** n)
        g = lambda k: f(xi + hn * (2.0 * k - 1.0))
        summ = hn * sum([g(k) for k in xrange(1, 2 ** (n - 1))])
        return 0.5 * R(n - 1, 0) + summ
 
    def Rnm(n, m):
        frac = 1.0 / ((4.0 ** m) - 1.0)
        rec1 = R(n, m - 1)
        rec2 = R(n - 1, m - 1)
        rec = (4.0 ** m) * rec1 - rec2
        return frac * rec
 
    def R(n, m):
        if n == 0 and m == 0:
            return R00()
        elif n > 0 and m == 0:
            return Rn0(n)
        else:
            return Rnm(n, m)
    return R(ni, mi)

Contudo, essa implementação ingênua pode ser muito custosa para muitos \(n \) ou ordens de integração ( \(m \) ) muito grandes. Assim, uma versão mais eficiente do método de Romberg é listada abaixo:

def romberg_fast(f, xi, xe, n=20):
    """
    Numerical integration using Romberg's Integration.
 
    integral = romberg_fast(f, xi, xe, n=20)
 
    INPUT:
      * f: function to be integrated
      * xi: beginning of integration interval
      * xe: end of integration interval
      * n: number of iterations
 
    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
    """
    def calc_trapezoids(s, k, i):
        m = 2 ** (k - 1)
        si = trapezoid(f, xi, xe, m)
        return (s[i], si)
 
    def improve_convergence(s, var, i, k):
        sk = (4 ** (i - 1) * s[i - 1] - var) / (4 ** (i - 1) - 1)
        return (sk, s[i], s[k])
 
    def compute():
        var = 0.0
        s = [1.0 for i in xrange(0, n)]
        for k in xrange(1, n):
            for i in xrange(1, k + 1):
                if i == 1:
                    (var, s[i]) = calc_trapezoids(s, k, i)
                else:
                    (s[k], var, s[i]) = improve_convergence(s, var, i, k)
        return s[n - 1]
 
    return compute()

Um exemplo de utilização dessas funções pode ser obtido em http://www.sawp.com.br/code/integral/integral_romberg.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).