El método adaptativo de Simpson , también llamado regla adaptativa de Simpson , 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 métodos adaptativos más modernos basados en la cuadratura de Gauss-Kronrod y la cuadratura de Clenshaw-Curtis son ahora generalmente preferidos. El método adaptativo de Simpson usa una estimación del error que obtenemos al calcular una integral definida usando 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 adaptativo de Simpson 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 .
Un criterio para determinar cuándo dejar de subdividir un intervalo, sugerido por JN Lyness, [3] es
dónde es un intervalo con punto medio , , , y son las estimaciones dadas por la regla de Simpson sobre los intervalos correspondientes y es la tolerancia deseada para el intervalo.
La regla de Simpson es una regla de cuadratura interpolatoria que es exacta cuando el integrando es un polinomio de grado tres o menor. Utilizando la extrapolación de Richardson , la estimación de Simpson más precisa para seis valores de función se combina con la estimación menos precisa para tres valores de función aplicando la corrección . La estimación así obtenida es exacta para polinomios de grado cinco o menos.
Consideración numérica
Algunas entradas no convergerán rápidamente en el método adaptativo de Simpson, lo que provocará que la tolerancia se desborde y produzca un bucle infinito. Los métodos simples para protegerse contra este problema incluyen agregar una limitación de profundidad (como en la muestra 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 , dando a = b .
El artículo de 1969 de Lychee 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 ε 0 .
- Para cada subintervalo [ a , b ] , defina D ( a , b ) , la estimación del error, como. Defina E = 180 ε 0 / ( B - A ) . Los criterios de terminación originales se convertirían entonces en D ≤ E .
- Si D ( a , m ) ≥ D ( a , b ) , se ha alcanzado el nivel de redondeo o se encuentra un cero para f (4) en el intervalo. Es necesario un cambio en la tolerancia ε 0 a ε ' 0 .
- Las rutinas recursivas ahora necesitan devolver un nivel D para el intervalo actual. Una variable rutina-estático E ' = 180 ε' 0 / ( B - A ) se define y se inicializa a E .
- (Modificación 4 i, ii) Si se usa más recursividad en un intervalo:
- Si parece que se ha alcanzado el redondeo, cambie la E ' por D ( a , m ) . [a]
- De lo contrario, ajuste E ' al máximo ( E , D ( a , m )) .
- Es necesario cierto control de los ajustes. Deben inhibirse los aumentos significativos y las disminuciones menores de las tolerancias.
- Para calcular el ε ' 0 efectivo en todo el intervalo:
- Registre cada x i en el que E ' se cambia en una matriz de ( x i , ε i ' ) pares. La primera entrada debe ser ( a , ε ' 0 ) .
- El ε eff real es la media aritmética de todo ε ' 0 , ponderado por el ancho de los intervalos.
- Si la E ' actual 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á desactivado.
- No se debe utilizar el término de corrección.
La maniobra de elevación de épsilon permite que la rutina se utilice en un modo de "mejor esfuerzo": dada una tolerancia inicial cero, la rutina intentará obtener la respuesta más precisa y devolverá un nivel de error real. [3] : 492
Implementaciones de muestra
Una técnica de implementación común que se muestra a continuación es transmitir 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 de nuevo para los subintervalos. Hacerlo 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 permanece lineal con la capa de recursiones.
Pitón
Aquí hay una implementación del método adaptativo de Simpson en Python .
from __future__ import division # python 2 compat # versión adaptativa "estructurada", traducida de Racket def _quad_simpsons_mem ( f , a , fa , b , fb ): "" "Evalúa la regla de Simpson, devolviendo también 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 adaptativa de Simpson. Se conservan los valores de función al principio, 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 = izquierda + derecha - entero si abs ( delta ) < = 15 * eps : return izquierda + derecha + 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 aab usando la regla de Simpson adaptativa con un error máximo de eps." "" fa , fb = f ( a ), f ( b ) m , fm , entero = _quad_simpsons_mem ( f , a , fa , b , fb ) return _quad_asr ( f , a , fa , b , fb , eps , entero , m , fm )desde la importación matemática sin print ( quad_asr ( sin , 0 , 1 , 1e-09 ))
C
Aquí hay una implementación del método adaptativo de Simpson en C99 que evita evaluaciones redundantes de f y cálculos en cuadratura. Incluye las tres defensas "simples" contra problemas numéricos.
#include // incluir archivo para fabs y sin#include // incluir archivo para printf y perror#include / ** Regla adaptativa de Simpson, núcleo recursivo * / float adaptiveSimpsonsAux ( float ( * f ) ( float ), float a , float b , float eps , float entero , float fa , float fb , float fm , int rec ) { float m = ( a + b ) / 2 , h = ( b - a ) / 2 ; flotar lm = ( a + m ) / 2 , rm = ( m + b ) / 2 ; // problema numérico serio: no convergerá si (( eps / 2 == eps ) || ( a == lm )) { errno = EDOM ; volver entero ; } flm flotante = f ( lm ), frm = f ( rm ); flotar a la izquierda = ( h / 6 ) * ( fa + 4 * flm + fm ); flotar a la derecha = ( h / 6 ) * ( fm + 4 * frm + fb ); flotar delta = izquierda + derecha - entero ; if ( rec <= 0 && errno ! = EDOM ) errno = ERANGE ; // límite de profundidad demasiado superficial // 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 , izquierda , fa , fm , flm , rec -1 ) + adaptiveSimpsonsAux ( f , m , b , eps / 2 , derecha , fm , fb , frm , rec -1 ) ; }/ ** Adaptive Simpson's Rule Wrapper * (llena las evaluaciones de funciones en caché) * / float adaptiveSimpsons ( float ( * f ) ( float ), // función ptr para integrar float a , float b , // intervalo [a, b] float épsilon , // tolerancia a errores int maxRecDepth ) { // límite de recursividad errno = 0 ; flotar h = b - a ; si ( h == 0 ) devuelve 0 ; flotar fa = f ( a ), fb = f ( b ), fm = f (( a + b ) / 2 ); flotar S = ( h / 6 ) * ( fa + 4 * fm + fb ); return adaptiveSimpsonsAux ( f , a , b , epsilon , S , fa , fb , fm , maxRecDepth ); }/ ** ejemplo de uso * / #include // para el ejemplo hostil (función rand)static int callcnt = 0 ; static float sinfc ( float x ) { callcnt ++ ; return sinf ( x ); } flotante estático frand48c ( flotante x ) { callcnt ++ ; return drand48 (); } int main () { // Sea yo la integral de sin (x) de 0 a 2 float I = adaptiveSimpsons ( sinfc , 0 , 2 , 1e-5 , 3 ); printf ( "integrar (sinf, 0, 2) =% lf \ n " , I ); // imprime el resultado perror ( "adaptiveSimpsons" ); // ¿Tuvo éxito? (profundidad = 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 ( "integrar (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 envoltorios para integrar en múltiples dimensiones. [4]
Raqueta
A continuación se muestra una implementación del método Simpson adaptativo en Racket con un contrato de software de comportamiento. La función exportada calcula la integral indeterminada para alguna 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 ) ))) ( #: épsilon ( ε 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-1call-to-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 ε entero M f @ M ) ( define-valores ( leftM f @ leftM left * ) ( simpson-1call-to-f f L f @ L M f @ M )) ( definir-valores ( rightM f @ rightM right * ) ( simpson-1call-to-f f M f @ M R f @ R )) ( definir delta * ( - ( + izquierda * derecha * ) entero )) ( cond [ ( <= ( abs delta * ) ( * 15 ε )) ( + izquierda * derecha * ( / delta * 15 )) ] [else ( definir épsilon1 ( / ε 2 )) ( + ( asr f L f @ L M f @ M épsilon1 izquierda * izquierdaM f @ izquierdaM ) ( asr f M f @ M R f @ R épsilon1 derecha * derechaM f @ derechaM )) ] ));; evaluar la mitad de un intervalo (1 función eval) ( definir ( simpson-1call-to-f f L f @ L R f @ R ) ( definir M ( mid L R )) ( definir f @ M ( f M )) ( valores M f @ M ( * ( / ( abs ( - R L )) 6 ) ( + f @ L ( * 4 f @ M ) f @ R ))))( definir ( medio L R ) ( / ( + L R ) 2. ))
Algoritmos relacionados
- 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 y bisector original. El algoritmo 103 consiste en una rutina más grande con una subrutina anidada (bucle AA), que se vuelve recursiva mediante el uso de la instrucció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 recursividad y S (2) es la estimación más precisa. [1]
- El algoritmo 145 de McKeeman (1962) es un integrador igualmente recursivo que divide el intervalo en tres en lugar de dos partes. La recursividad está escrita de una manera más familiar. [6] El algoritmo de 1962, que se considera demasiado cauteloso, utiliza para la terminación, por lo que una mejora de 1963 utiliza en lugar de. [3] : 485,487
- Lyness (1969) es casi el integrador actual. Creado como un conjunto de cuatro modificaciones de McKeeman 1962, reemplaza la trisección con la bisección para reducir los costos computacionales (Modificaciones 1 + 2, coincidiendo con el integrador de Kuncir) y mejora las estimaciones de error de 1962/63 de McKeeman al quinto orden (Modificación 3), en una forma 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]
Notas
Bibliografía
- ^ a b G.F. Kuncir (1962), "Algoritmo 103: Integrador de reglas de Simpson", Comunicaciones del ACM , 5 (6): 347, doi : 10.1145 / 367766.368179
- ^ a b Para un integrador adaptativo no recursivo anterior que recuerda más a los solucionadores de ODE , consulte S. Henriksson (1961), "Contribución n. ° 2: Integración numérica de Simpson con longitud variable de paso", BIT Numerical Mathematics , 1 : 290
- ^ a b c d e f JN Lyness (1969), "Notas sobre la rutina adaptativa de cuadratura de Simpson", Journal of the ACM , 16 (3): 483–495, doi : 10.1145 / 321526.321537
- ^ a b Berrill, Mark A. "RayTrace-miniapp: src / AtomicModel / interp.hpp · de5e8229bccf60ae5c1c5bab14f861dc0326d5f9" . ORNL GitLab .
- ^ Felleisen, Matthias. "[raqueta] integración adaptativa de simpson" . "Usuarios" de la lista de correo de Racket . Consultado el 26 de septiembre de 2018 .
- ^ McKeeman, William Marshall (1 de diciembre de 1962). "Algoritmo 145: integración numérica adaptativa por regla de Simpson". Comunicaciones de la ACM . 5 (12): 604. doi : 10.1145 / 355580.369102 .
enlaces externos
- Módulo para la regla adaptativa de Simpson