Método de Simpson adaptativo

El método de Simpson adaptativo , también llamado regla de Simpson adaptativa , es un método de integración numérica propuesto por GF Kuncir en 1962. ^[1] Es probablemente el primer algoritmo adaptativo recursivo para la integración numérica que aparece impreso, ^[2] aunque ahora se prefieren generalmente métodos adaptativos más modernos basados ​​en la cuadratura de Gauss-Kronrod y la cuadratura de Clenshaw-Curtis . El método de Simpson adaptativo utiliza una estimación del error que obtenemos al calcular una integral definida utilizando la regla de Simpson . Si el error excede una tolerancia especificada por el usuario, el algoritmo requiere subdividir el intervalo de integración en dos y aplicar el método de Simpson adaptativo a cada subintervalo de manera recursiva. La técnica suele ser mucho más eficiente que la regla de Simpson compuesta , ya que utiliza menos evaluaciones de funciones en lugares donde la función está bien aproximada por una función cúbica .

La regla de Simpson es una regla de cuadratura interpolativa que es exacta cuando el integrando es un polinomio de grado tres o menor. Mediante la extrapolación de Richardson , la estimación más precisa de Simpson para seis valores de función se combina con la estimación menos precisa para tres valores de función aplicando la corrección . Por lo tanto, la estimación obtenida es exacta para polinomios de grado cinco o menor. ${\displaystyle S(a,m)+S(m,b)}$${\displaystyle S(a,b)}$${\displaystyle [S(a,m)+S(m,b)-S(a,b)]/15}$

Un criterio para determinar cuándo dejar de subdividir un intervalo, sugerido por JN Lyness, ^[3] es $${\displaystyle |S(a,m)+S(m,b)-S(a,b)|<15\varepsilon }$$

donde es un intervalo con punto medio , mientras que , , y dadas por la regla de Simpson son las estimaciones de , , y respectivamente, y es la tolerancia de error máxima deseada para el intervalo. ${\displaystyle [a,b]}$${\displaystyle m={\frac {a+b}{2}}}$${\displaystyle S(a,b)\,\!}$${\displaystyle S(a,m)\,\!}$${\displaystyle S(m,b)}$${\textstyle \int _{a}^{b}f(x)\,dx}$${\textstyle \int _{a}^{m}f(x)\,dx}$${\textstyle \int _{m}^{b}f(x)\,dx}$${\displaystyle \varepsilon }$

Nota, . ${\displaystyle \varepsilon _{i+1}=\varepsilon _{i}/2}$

Para realizar el método de Simpson adaptativo, haga lo siguiente: si , agregue y a la suma de las reglas de Simpson que se utilizan para aproximar la integral; de lo contrario, realice la misma operación con y en lugar de . ${\displaystyle |S(a,m)+S(m,b)-S(a,b)|<15\varepsilon _{i}}$${\displaystyle S(a,m)}$${\displaystyle S(m,b)}$${\textstyle \left|S{\left(a,{\frac {m-a}{2}}\right)}+S{\left({\frac {m-a}{2}},m\right)}-S(a,m)\right|<15\varepsilon _{i+1}}$${\textstyle \left|S{\left(m,{\frac {b-m}{2}}\right)}+S{\left({\frac {b-m}{2}},b\right)}-S(m,b)\right|<15\varepsilon _{i+1}}$${\displaystyle |S(a,m)+S(m,b)-S(a,b)|<15\varepsilon _{i}}$

Algunas entradas no convergerán rápidamente en el método adaptativo de Simpson, lo que provocará que la tolerancia se desborde y se produzca un bucle infinito. Algunos métodos simples para protegerse contra este problema incluyen agregar una limitación de profundidad (como en el ejemplo C y en McKeeman), verificar que ε /2 ≠ ε en aritmética de punto flotante, o ambos (como Kuncir). El tamaño del intervalo también puede acercarse al épsilon de la máquina local , lo que da a = b .

El artículo de Lyness de 1969 incluye una "Modificación 4" que aborda este problema de una manera más concreta: ^{[3] : 490–2}

• Sea el intervalo inicial [ A , B ] . Sea la tolerancia original ε ₀ .
• Para cada subintervalo [ a , b ] , defina D ( a , b ) , la estimación del error, como . Defina E = 180 ε ₀ / ( B - A ) . El criterio de terminación original sería entonces D ≤ E .${\textstyle {\frac {12}{a-b}}[S^{(2)}(a,b)-S(a,b)]={\frac {1}{4}}[f(a)-4f(lm)+6f(m)-4f(rm)+f(b)]}$
• Si D ( a , m ) ≥ D ( a , b ) , se ha alcanzado el nivel de redondeo o se encuentra un cero para f ⁽⁴⁾ en el intervalo. Es necesario cambiar la tolerancia ε ₀ a ε′ _{0 .}
  • Las rutinas recursivas ahora deben devolver un nivel D para el intervalo actual. Se define una variable estática de rutina E' = 180 ε' ₀ / ( B - A ) y se inicializa en E .
  • (Modificación 4 i, ii) Si se utiliza más recursión en un intervalo:
    1. Si parece que se ha alcanzado el redondeo, cambie E' por D ( a , m ) . ^[a]
    2. De lo contrario, ajuste E' a max( E , D ( a , m )) .
  • Es necesario un cierto control de los ajustes. Se deben evitar aumentos significativos y pequeñas disminuciones de las tolerancias.
• Para calcular el ε′ ₀ efectivo en todo el intervalo:
  • Registra cada x _i en el que E' se transforma en una matriz de pares ( x _i , ε _i ' ) . La primera entrada debe ser ( a , ε′ ₀ ) .
  • La ε _eff real es la media aritmética de todos los ε′ ₀ , ponderada por el ancho de los intervalos.
• Si la corriente E' para un intervalo es mayor que E , entonces la aceleración/corrección de quinto orden no se aplicaría: ^[b]
  • El factor "15" en los criterios de terminación está deshabilitado.
  • No se debe utilizar el término de corrección.

La maniobra de aumento de épsilon permite utilizar la rutina en un modo de "máximo esfuerzo": dada una tolerancia inicial de cero, la rutina intentará obtener la respuesta más precisa y devolverá un nivel de error real. ^{[3] : 492}

Una técnica de implementación común que se muestra a continuación es pasar f( a ), f( b ), f( m ) junto con el intervalo [ a , b ] . Estos valores, utilizados para evaluar S ( a , b ) en el nivel principal, se utilizarán nuevamente para los subintervalos. Al hacerlo, se reduce el costo de cada llamada recursiva de 6 a 2 evaluaciones de la función de entrada. El tamaño del espacio de pila utilizado se mantiene lineal a la capa de recursiones.

Aquí hay una implementación del método de Simpson adaptativo en Python .

de  __future__  import  division  # compatibilidad con python 2 
# versión adaptativa "estructurada", traducida de Racket 
def  _quad_simpsons_mem ( f ,  a ,  fa ,  b ,  fb ): 
"""Evalúa la regla de Simpson, y también devuelve m y f(m) para reutilizar""" m = ( a + b ) / 2 fm = f ( m ) return ( m , fm , abs ( b - a ) / 6 * ( fa + 4 * fm + fb ))    
          
      
                   

def  _quad_asr ( f ,  a ,  fa ,  b ,  fb ,  eps ,  whole ,  m ,  fm ): 
"""     Implementación recursiva eficiente de la regla de Simpson adaptativa.     Se conservan los valores de la función al inicio, medio y final de los intervalos.     """ lm , flm , left = _quad_simpsons_mem ( f , a , fa , m , fm ) rm , frm , right = _quad_simpsons_mem ( f , m , fm , b , fb ) delta = left + right - whole if abs ( delta ) <= 15 * eps : return left + right + delta / 15 return _quad_asr ( f , a , fa , m , fm , eps / 2 , left , lm , flm ) + \
            _quad_asr ( f , m , fm , b , fb , eps / 2 , derecha , rm , frm )    



             
            
          
         
               
                       

def  quad_asr ( f ,  a ,  b ,  eps ): 
"""Integra f de a a b usando la regla de Simpson adaptativa con un error máximo de eps.""" fa , fb = f ( a ), f ( b ) m , fm , whole = _quad_simpsons_mem ( f , a , fa , b , fb ) return _quad_asr ( f , a , fa , b , fb , eps , whole , m , fm )    
        
            
             

desde  matemáticas  importar  sin 
imprimir ( quad_asr ( sin ,  0 ,  1 ,  1e-09 ))

Aquí se presenta una implementación del método adaptativo de Simpson en C99 que evita evaluaciones redundantes de f y cálculos de cuadratura. Incluye las tres defensas "simples" contra problemas numéricos.

#include <math.h>   // archivo de inclusión para fabs y sin #include <stdio.h> // archivo de inclusión para printf y perror #include <errno.h> 
 
 

/** Regla de Simpson adaptativa, núcleo recursivo */ 
float adaptiveSimpsonsAux ( float ( * f )( float ), float a , float b , float eps , float whole , float fa , float fb , float fm , int rec ) { float m = ( a + b ) / 2 , h = ( b - a ) / 2 ; float lm = ( a + m ) / 2 , rm = ( m + b ) / 2 ; // problema numérico serio: no convergerá si (( eps / 2 == eps ) || ( a == lm )) { errno = EDOM ; return whole ; } float flm = f ( lm ), frm = f ( rm ); float left = ( h / 6 ) * ( fa + 4 * flm + fm ); flotante derecha = ( h / 6 ) * ( fm + 4 * frm + fb ); flotante delta = izquierda + derecha - entero ;        
                                    
                   
                 
    
                  
               
              
             
           

    if ( rec <= 0 && errno != EDOM ) errno = ERANGE ; // límite de profundidad demasiado bajo // Lyness 1969 + extrapolación de Richardson; ver artículo if ( rec <= 0 || fabs ( delta ) <= 15 * eps ) return left + right + ( delta ) / 15 ; return adaptiveSimpsonsAux ( f , a , m , eps / 2 , left , fa , fm , flm , rec -1 ) + adaptiveSimpsonsAux ( f , m , b , eps / 2 , right , fm , fb , frm , rec -1 ); }            
    
           
             
               
                   


/** Envoltorio de reglas de Simpson adaptable 
* (rellena las evaluaciones de funciones almacenadas en caché) */ 
float adaptiveSimpsons ( float ( * f )( float ), // función ptr para integrar float a , float b , // intervalo [a,b] float epsilon , // tolerancia de error int maxRecDepth ) { // límite de recursión errno = 0 ; float h = b - a ; if ( h == 0 ) return 0 ; float fa = f ( a ), fb = f ( b ), fm = f (( a + b ) / 2 ); float S = ( h / 6 ) * ( fa + 4 * fm + fb ); return adaptiveSimpsonsAux ( f , a , b , epsilon , S , fa , fb , fm , maxRecDepth ); }       
                                
                                 
                              
      
         
         
               
           
             


/** ejemplo de uso */ 
#include <stdlib.h> // para el ejemplo hostil (función rand) static int callcnt = 0 ; static float sinfc ( float x ) { callcnt ++ ; return sinf ( x ); } static float frand48c ( float x ) { callcnt ++ ; return drand48 (); } int main () { // Sea I la integral de sin(x) de 0 a 2 float I = adaptiveSimpsons ( sinfc , 0 , 2 , 1e-5 , 3 ); printf ( "integrate(sinf, 0, 2) = %lf \n " , I ); // imprima el resultado perror ( "adaptiveSimpsons" ); // ¿Fue exitoso? (depth=1 es demasiado superficial) printf ( "(%d evaluaciones) \n " , callcnt ); 
    
         
         
  
    
           
        
                       
     

    callcnt = 0 ; srand48 ( 0 ); I = adaptiveSimpsons ( frand48c , 0 , 0.25 , 1e-5 , 25 ); // una función aleatoria printf ( "integrate(frand48, 0, 0.25) = %lf \n " , I ); perror ( "adaptiveSimpsons" ); // no convergerá printf ( "(%d evaluaciones) \n " , callcnt ); return 0 ; }   
           
     
                            
     
     

Esta implementación se ha incorporado a un trazador de rayos C++ destinado a la simulación de láser de rayos X en el Laboratorio Nacional Oak Ridge ^[4] , entre otros proyectos. La versión ORNL se ha mejorado con un contador de llamadas, plantillas para diferentes tipos de datos y contenedores para la integración en múltiples dimensiones. ^[4]

A continuación se muestra una implementación del método Simpson adaptativo en Racket con un contrato de software conductual. La función exportada calcula la integral indeterminada para una función dada f . ^[5]

;; ----------------------------------------------------------------------------- 
;; interfaz, con contrato 
( proporcionar/contrato 
[adaptive-simpson ( ->i (( f ( -> ¿ real? ¿real? )) ( L ¿ real? ) ( R ( L ) ( y/c ¿real? ( >/c L )))) ( #:epsilon ( ε ¿ real? )) ( r ¿real? )) ] )  
               
         
        

;; ----------------------------------------------------------------------------- 
;; implementación 
( define ( adaptive-simpson f L R #:epsilon [ε . 000000001] ) ( define f@L ( f L )) ( define f@R ( f R )) ( define-valores ( M f@M entero ) ( simpson-1llamada-a-f f L f@L R f@R )) ( asr f L f@L R f@R ε entero M f@M ))       
     
     
           
           

;; la implementación "eficiente" 
( 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 )) ] ))          
              
           
        
  
               
         
                       
                      

;; evaluar la mitad de un intervalo (1 func eval) 
( 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 ))))      
      
     
                 

( definir ( medio L R ) ( / ( + L R ) 2. ))        

• Henriksson (1961) es una variante no recursiva de la regla de Simpson. Se "adapta" integrando de izquierda a derecha y ajustando el ancho del intervalo según sea necesario. ^[2]
• El algoritmo 103 de Kuncir (1962) es el integrador adaptativo recursivo bisectriz original. El algoritmo 103 consta de una rutina más grande con una subrutina anidada (bucle AA), que se vuelve recursiva mediante el uso de la declaración goto . Protege contra el desbordamiento de los anchos de intervalo (bucle BB) y se interrumpe tan pronto como se excede el eps especificado por el usuario. El criterio de terminación es , donde n es el nivel actual de recursión y S ⁽²⁾ es la estimación más precisa. ^[1]${\displaystyle |S^{(2)}(a,b)-S(a,b)|<2^{-n}\epsilon \,}$
• El algoritmo 145 de McKeeman (1962) es un integrador recursivo similar que divide el intervalo en tres partes en lugar de dos. La recursión está escrita de una manera más familiar. ^[6] El algoritmo de 1962, que se consideró demasiado cauteloso, utiliza para la terminación, por lo que una mejora de 1963 utiliza en su lugar. ^{[3] : 485, 487}${\displaystyle |S^{(3)}(a,b)-S(a,b)|<3^{-n}\epsilon \,}$${\displaystyle {\sqrt {3}}^{\,-n}\epsilon }$
• Lyness (1969) es casi el integrador actual. Creado como un conjunto de cuatro modificaciones de McKeeman 1962, reemplaza la trisección por la bisección para reducir los costos computacionales (Modificaciones 1+2, coincidentes con el integrador de Kuncir) y mejora las estimaciones de error de McKeeman 1962/63 al quinto orden (Modificación 3), de una manera relacionada con la regla de Boole y el método de Romberg . ^{[3] : 489} La Modificación 4, no implementada aquí, contiene disposiciones para el error de redondeo que permite elevar el ε al mínimo permitido por la precisión actual y devolver el nuevo error. ^[3]

• Módulo para la regla de Simpson adaptativa