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 o resultado da aproximação da integral utilizando o método de Romberg no intervalo . Assim, o resultado da integral será
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/.