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 e o erro associado com a aplicação múltipla da regra do trapézio podem ser representados de um modo geral por
onde é o valor exato da integral, é a aproximação a partir da aplicação da regra do trapézio (em função do passo ) e é o erro de truncamento. Se fizermos duas estimativas independentes usando passos distintos e , temos que
Como o erro da aplicação da regra do trapézio é
onde é a derivada média da função nos pontos avaliados. Supondo que é constante, independente do tamanho do passo, a equação acima pode ser usada para determinar a razão entre os erros:
ou melhor
que pode ser substituida na equação das integrais:
isolando , temos
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
para fornecer uma estimativa mais acurada
Para o caso em que , temos que
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:
onde e são as integrais mais e menos acuradas, respectivamente, e é a integral melhorada. O índice representa a respectiva fórmula de Newton-Cotes, isto é, para temos a regra do trapézio, para temos a regra de Simpson, e assim por diante. O índice serve para indicar a estimativa mais acurada da menos acurada (quanto maior , mais acurada). Assim, o método pode ser definido indutivamente da seguinte forma
onde e são os limites do intervalo de integração, , $m \geq 1 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 ou ordens de integração ( ) 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/.