5.5 Integração Numérica — Integração de Romberg Adaptativa

1. Introdução

No post sobre a integração de Romberg, vimos que este método utiliza uma estratégia para aumentar a acurácia da regra do trapézio através de sucessivas aplicações. Também apresentamos o Método de Simpson Adaptativo, que é um algoritmo recursivo de integração numérica que utiliza múltiplas aplicações do método de Simpson, dividindo o intervalo de integração em vários outros sub-intervalos. Ambas estratégias podem ser adaptadas para aprimorar a convergência do método do trapézio. Contudo, podemos substituir o método de Simpson pelo método de Romberg, o que permite um aumento na velocidade de convergência ainda mais significante no algoritmo adaptativo.

 

2. O Método

Seja \(S(a,b) \) o resultado da aproximação da integral utilizando o método de Romberg no intervalo \([a,b] \) . Assim, o resultado da integral será

\( I = S(a,b) = \int_a^b f(x) dx \)

Onde o critério de parada para a sub-divisão dos intervalos é o mesmo utilizando no Simpson adaptativo.

 

3. Implementação

As funções que implementam o método de Romberg adaptativo são:

def trapezoid(f, xi, xe, n=10000):
    """
    Numerical integration using the Trapezoid Rule.
 
    integral = trapezoid(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
    sumatory = sum([f(xi + i * h) for i in xrange(1, n)])
    trapezoids = f(xi) + 2.0 * sumatory + f(xi + n * h)
    integral = trapezoids * (h / 2.0)
    return integral
def romberg(f, xi, xe, n=10):
    """
    Numerical integration using Romberg's Integration.
 
    integral = romberg(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()
def adaptative_romberg(f, xi, xe, errto=10e-10):
    """
    Numerical integration using the Adaptative Romberg's.
 
    integral = adaptative_rombergf, xi, xe, errto=10e-10)
 
    INPUT:
      * f: function to be integrated
      * xi: beginning of integration interval
      * xe: end of integration interval
      * errto: numerical error tolerated in the integration
 
    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 apply_romberg(g, a, b):
        middle_point = (a + b) / 2.0
        left = romberg(g, a, middle_point)
        right = romberg(g, middle_point, b)
        return (left, middle_point, right)
 
    def recursion(g, a, b, errno, integral):
        (left, middle, right) = apply_romberg(g, a, b)
        adap = left + right - integral
        err = abs(adap) / 15.0
        rec = lambda x, y, z: recursion(g, x, y, errno / 2.0, z)
        if err < errno:
            return left + right + adap / 15.0
        else:
            return rec(a, middle, left) + rec(middle, b, right)
    integral = recursion(f, xi, xe, errto, romberg(f, xi, xe))
    return integral

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

 

4. 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).