5.3 Integração Numérica — Método de Simpson Adaptativo
1. Introdução
A regra de Simpson adaptativa é um algoritmo recursivo de integração numérica que usa uma estimativa do erro para calcular uma integral definida através de múltiplas aplicações do método de Simpson. Como em muitas abordagens numéricas, se o erro estimado ultrapassa um erro tolerado, o algoritmo divide o intervalo de integração em dois sub-intervalos e reaplica o método de Simpson nesses intervalos de forma independente, somando os resultados. Essa técnica acaba por ser mais eficiente que outros métodos de integração por convergir quadraticamente ao resultado.
2. O Método
Seja o resultado da aproximação da integral utilizando a regra de Simpson — ou — no intervalo . Assim,
O critério de parada da sub-divisão dos intervalos é calculado pela seguinte equação
com
</center>
onde é o ponto médio do intervalo e é o erro tolerado para o intervalo.
3. Implementação
A implementação utilizando as fórmulas indutivas é
def adaptative_simpson(f, xi, xe, errto=10e-10):
"""
Numerical integration using the Adaptative Simpson's Rule.
integral = adaptative_simpson(f, 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 simpson_rule(g, a, b):
middle_point = (a + b) / 2.0
h = abs(a - b) / 6.0
n = g(a) + 4.0 * g(middle_point) + g(b)
integral = h * n
return integral
def apply_simpson_rule(g, a, b):
middle_point = (a + b) / 2.0
left = simpson_rule(g, a, middle_point)
right = simpson_rule(g, middle_point, b)
return (left, middle_point, right)
def recursion(g, a, b, errno, integral):
(left, middle, right) = apply_simpson_rule(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, simpson_rule(f, xi, xe))
return integral
Um exemplo de utilização dessa função pode ser obtido em http://www.sawp.com.br/code/integral/integral_adaptative_simpson.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/.