Quadratura de Simpson adaptativa

Quadratura de Simpson adaptativa, também chamada de Regra adaptativa de Simpson, é um método de Integração numérica, proposto por G.F. Kuncir em 1962.^[1] Este é, provavelmente, o primeiro algoritmo recursivo adaptativo para integração numérica impresso,^[2] porém, métodos adaptativos mais modernos, baseados em Quadraturas de Gauss-Kronrod e Quadratura de Clenshaw–Curtis, são geralmente preferidos. O método da Quadratura de Simpson Adaptativa usa uma estimativa do erro obtido ao calcular integrais definidas pelo Método de Simpson. Se o erro excede a tolerância definida pelo usuário, o algoritmo subdivide os intervalos em dois e aplica a Quadratura de Simpson adaptativa em cara intervalo, de maneira recursiva. Usualmente, a técnica é muito mais eficiente que a Regra composta de Simpson., uma vez que realiza menos avaliações em pontos em que a função é bem aproximada por uma Função cúbica.

Um critério sobre quando parar de subdividir um intervalo, sugerido por J.N. Lyness,^[3] é

|S(a,c)+S(c,b)-S(a,b)|/15<\epsilon \,

onde $[a,b]\,\!$ é um intervalo centrado em $c\,\!$ , $S(a,b)\,\!$ , $S(a,c)\,\!$ , e onde $S(c,b)\,\!$ são estimativas dadas pelo Método de Simpson nos intervalos correspondentes e $\epsilon \,\!$ é a tolerância desejada para o intervalo.

O Método de Simpson é uma quadratura interpoladora que apresenta resultados exatos quando o integrante é um polinômio de grau três ou menor. Usando Extrapolação de Richardson, as estimativas de Simpson mais precisas de $S(a,c)+S(c,b)\,$ para seis valores da função são combinadas com estimativas menos precisas $S(a,b)\,$ para três valores da função, aplicando a correção $[S(a,c)+S(c,b)-S(a,b)]/15\,$ . Assim, a estimativa obtida é exata para polinômios de grau cinco ou menos.

Exemplos de Implementação

Python

Segue uma implementação da Quadratura de Simpson adaptativa em Python. Observe que este é um código explicativo, em detrimento de sua eficiência. Todas as chamadas para recursive_asr implicam seis avaliações da função. Para uso, recomenda-se uma modificação, de tal maneira que um mínimo de duas avaliações de função sejam realizadas.

def simpsons_rule(f,a,b):
    c = (a+b) / 2.0
    h3 = abs(b-a) / 6.0
    return h3*(f(a) + 4.0*f(c) + f(b))

def recursive_asr(f,a,b,eps,whole):
    "Recursive implementation of adaptive Simpson's rule."
    c = (a+b) / 2.0
    left = simpsons_rule(f,a,c)
    right = simpsons_rule(f,c,b)
    if abs(left + right - whole) <= 15*eps:
        return left + right + (left + right - whole)/15.0
    return recursive_asr(f,a,c,eps/2.0,left) + recursive_asr(f,c,b,eps/2.0,right)

def adaptive_simpsons_rule(f,a,b,eps):
    "Calculate integral of f from a to b with max error of eps."
    return recursive_asr(f,a,b,eps,simpsons_rule(f,a,b))

from math import sin
print adaptive_simpsons_rule(sin,0,1,.000000001)

C

Segue também uma implementação da Quadratura de Simpson adaptativa em C99 que evita avaliações redundantes de f e cálculos de quadraturas.

#include <math.h>  // include file for fabs and sin
#include <stdio.h> // include file for printf
 
//
// Função recursiva auxiliar parar a função adaptiveSimpsons() abaixo
//                                                                                                 
double adaptiveSimpsonsAux(double (*f)(double), double a, double b, double epsilon,                 
                         double S, double fa, double fb, double fc, int bottom) {                 
  double c = (a + b)/2, h = b - a;                                                                  
  double d = (a + c)/2, e = (c + b)/2;                                                              
  double fd = f(d), fe = f(e);                                                                      
  double Sleft = (h/12)*(fa + 4*fd + fc);                                                           
  double Sright = (h/12)*(fc + 4*fe + fb);                                                          
  double S2 = Sleft + Sright;                                                                       
  if (bottom <= 0 || fabs(S2 - S) <= 15*epsilon)                                                    
    return S2 + (S2 - S)/15;                                                                        
  return adaptiveSimpsonsAux(f, a, c, epsilon/2, Sleft,  fa, fc, fd, bottom-1) +                    
         adaptiveSimpsonsAux(f, c, b, epsilon/2, Sright, fc, fb, fe, bottom-1);                     
}

//
// Quadratura de Simpson Adaptativa
//
double adaptiveSimpsons(double (*f)(double),   // ptr para a função
                           double a, double b,  // intervalo [a,b]
                           double epsilon,  // tolerancia de erro
                           int maxRecursionDepth) {   // recursão cap        
  double c = (a + b)/2, h = b - a;                                                                  
  double fa = f(a), fb = f(b), fc = f(c);                                                           
  double S = (h/6)*(fa + 4*fc + fb);                                                                
  return adaptiveSimpsonsAux(f, a, b, epsilon, S, fa, fb, fc, maxRecursionDepth);                   
}

int main(){
 double I = adaptiveSimpsons(sin, 0, 1, 0.000000001, 10); // calcular a integral de sin(x)
                                                          // de 0 a 1 e armazenar na
                                                          // nova variável I
 printf("I = %lf\n",I); // print the result
 return 0;
}

Racket

Seguem também uma implementação da Quadratura de Simpson Adaptativa em Racket. A função exportada calcula a integral indeterminada para algumas funções f dadas.

;; -----------------------------------------------------------------------------
;; interface

;; Método de Simpson para aproximara integral
(define (simpson f L R)
  (* (/ (- R L) 6) (+ (f L) (* 4 (f (mid L R))) (f R))))

(provide/contract
 [adaptive-simpson 
  (->i ((f (-> real? real?)) (L real?) (R  (L) (and/c real? (>/c L))))
       (#:epsilon (ε real?))
       (r real?))]
 [step (-> real? real?)])

;; -----------------------------------------------------------------------------
;; implementação

(define (adaptive-simpson f L R #:epsilon [ε .000000001])
  (define f@L (f L))
  (define f@R (f R))
  (define-values (M f@M whole) (simpson-1call-to-f f L f@L R f@R))
  (asr f L f@L R f@R ε whole M f@M))

;; eficiência computacional: 2 funções chamadas por passo 
(define (asr f L f@L R f@R ε whole M f@M)
  (define-values (leftM  f@leftM  left*)  (simpson-1call-to-f f L f@L M f@M))
  (define-values (rightM f@rightM right*) (simpson-1call-to-f f M f@M R f@R))
  (define delta* (- (+ left* right*) whole))
  (cond
    [(<= (abs delta*) (* 15 ε)) (+ left* right* (/ delta* 15))]
    [else (define epsilon1 (/ ε 2))
          (+ (asr f L f@L M f@M epsilon1 left*  leftM  f@leftM) 
             (asr f M f@M R f@R epsilon1 right* rightM f@rightM))]))

(define (simpson-1call-to-f f L f@L R f@R)
  (define M (mid L R))
  (define f@M (f M))
  (values M f@M (* (/ (abs (- R L)) 6) (+ f@L (* 4 f@M) f@R))))

(define (mid L R) (/ (+ L R) 2.))

O código é um resumo de um módulo "#lang racket" que inclui uma linha (require rackunit).

Bibliografia

↑ G.F. Kuncir (1962), «Algorithm 103: Simpson's rule integrator», Communications of the ACM, 5 (6): 347
↑ For an earlier, non-recursive adaptive integrator more reminiscent of ODE solvers, see S. Henriksson (1961), «Contribution no. 2: Simpson numerical integration with variable length of step», BIT Numerical Mathematics, 1: 290
↑ J.N. Lyness (1969), «Notes on the adaptive Simpson quadrature routine», Journal of the ACM, 16 (3): 483–495