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 integración numérica que aparece impreso, ^[2] aunque los métodos adaptativos más modernos Actualmente se prefieren generalmente los 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 interpolatoria que es exacta cuando el integrando es un polinomio de grado tres o inferior. 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 . Entonces, la estimación obtenida es exacta para polinomios de grado cinco o menos. ${\displaystyle S(a,m)+S(m,b)}$ ${\displaystyle S(a,b)}$ ${\displaystyle [S(a,m)+S(m,b)-S(a,b)]/15}$

Procedimiento matemático

Definición de términos

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 dado por la regla de Simpson son las estimaciones de , y respectivamente, y es la tolerancia máxima de error 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}$

Pasos procesales

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}}$

Consideración numérica

Algunas entradas no lograrán converger rápidamente en el método de Simpson adaptativo, lo que provocará un desbordamiento insuficiente de la tolerancia y producirá 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 ambas (como Kuncir). El tamaño del intervalo también puede aproximarse al épsilon de la máquina local , dando 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 ) . Los criterios de terminación originales pasarían entonces a ser 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 un cambio en la tolerancia ε ₀ a ε′ _{0 .}
  • Las rutinas recursivas ahora necesitan 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 mayor recursividad en un intervalo:
    1. Si parece que se ha alcanzado el redondeo, cambie E' a D ( a , m ) . ^[a]
    2. De lo contrario, ajuste E' a max( 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 ε′ ₀ efectivo en todo el intervalo:
  • Registre cada x _i en el que E' se cambie en una matriz de ( x _i , ε _i ' ) pares. La primera entrada debe ser ( a , ε′ ₀ ) .
  • El ε _eff real es la media aritmética de todos los ε′ ₀ , 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á deshabilitado.
  • No se debe utilizar el término de corrección.

La maniobra de elevación é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 código 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 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 permanece lineal con respecto a la capa de recursividad.

Pitón

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

de  __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 y también devuelve m y f(m) para reutilizar """ m = ( a + b ) / 2 fm = f ( m ) retorno ( m , fm , abs ( b - a ) / 6 * ( fa + 4 * fm + fb ))                           def  _quad_asr ( f ,  a ,  fa ,  b ,  fb ,  eps ,  entero ,  m ,  fm ): """  Implementación recursiva eficiente de la regla de Simpson adaptativa.  Se conservan los valores de función al inicio, medio y final de los intervalos.  "" " lm , flm , izquierda = _quad_simpsons_mem ( f , a , fa , m , fm ) rm , frm , derecha = _quad_simpsons_mem ( f , m , fm , b , fb ) delta = izquierda + derecha - entero si abs ( delta ) < = 15 * eps : regresar izquierda + derecha + delta / 15 regresar _quad_asr ( f , a , fa , m , fm , eps / 2 , izquierda , lm , flm ) + \ _quad_asr ( f , m , fm , b , fb , eps / 2 , derecha , rm , frm )                                                            def  quad_asr ( f ,  a ,  b ,  eps ): """Integrar f de a a b usando la regla de Simpson adaptativa con error máximo de eps.""" fa , fb = f ( a ), f ( b ) m , fm , entero = _quad_simpsons_mem ( f , a , fa , b , fb ) devuelve _quad_asr ( f , a , fa , b , fb , eps , entero , m , fm )                         de  matemáticas  importar  sin imprimir ( quad_asr ( sin ,  0 ,  1 ,  1e-09 ))

C

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

#include <math.h>  // incluye archivo para fabs y sin #include <stdio.h> // incluye archivo para printf y perror #include <errno.h>   /** Regla de Simpson adaptable, núcleo recursivo */ float adaptiveSimpsonsAux ( float ( * f )( float ), float a , float b , float eps , float entera , float fa , float fb , float fm , int rec ) { float m = ( a + b ) / 2 , h = ( b - a ) / 2 ; flotador lm = ( a + m ) / 2 , rm = ( m + b ) / 2 ; // serio problema numérico: no convergerá if (( eps / 2 == eps ) || ( a == lm )) { errno = EDOM ; regresar entero ; } flotador flm = f ( lm ), frm = f ( rm ); flotar izquierda = ( h / 6 ) * ( fa + 4 * flm + fm ); flotar a la derecha = ( h / 6 ) * ( fm + 4 * frm + fb ); flotador 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 izquierda + derecha + ( delta ) / 15 ; devolver adaptiveSimpsonsAux ( f , a , m , eps / 2 , izquierda , fa , fm , flm , rec -1 ) + adaptiveSimpsonsAux ( f , m , b , eps / 2 , derecha , fm , fb , frm , rec -1 ) ; }                                              /** Envoltorio de reglas de Simpson adaptable * (completa las evaluaciones de funciones en caché) */ float adaptiveSimpsons ( float ( * f )( float ), // función ptr para integrar float a , float b , // intervalo [a,b] float epsilon , // tolerancia a errores int maxRecDepth ) { // límite de recursividad errno = 0 ; flotar h = b - a ; si ( h == 0 ) devuelve 0 ; flotador fa = f ( a ), fb = f ( b ), fm = f (( a + b ) / 2 ); flotador S = ( h / 6 ) * ( fa + 4 * fm + fb ); devolver adaptiveSimpsonsAux ( f , a , b , épsilon , S , fa , fb , fm , maxRecDepth ); }                                                            /** ejemplo de uso */ #include <stdlib.h> // para el ejemplo hostil (función rand) static int callcnt = 0 ; flotante estático sinfc ( flotante x ) { callcnt ++ ; devolver pecado ( x ); } flotador estático frand48c ( flotante x ) { callcnt ++ ; devolver drand48 (); } int main () { // Sea I 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" ); // ¿Fue exitoso? (la profundidad = 1 es demasiado superficial) printf ( "(%d evaluaciones) \n " , callcnt );                                          llamada = 0 ; sarand48 ( 0 ); I = Simpsons adaptativos ( frand48c , 0 , 0,25 , 1e-5 , 25 ); // una función aleatoria printf ( "integrar(frand48, 0, 0.25) = %lf \n " , I ); perror ( "adaptativoSimpsons" ); // no convergerá printf ( "(%d evaluaciones) \n " , callcnt ); devolver 0 ; }                   

Esta implementación se ha incorporado a un trazador de rayos C++ destinado a la simulación de rayos X con láser en el Laboratorio Nacional de 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 integración 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 conductual. La función exportada calcula la integral indeterminada para alguna función dada f . ^[5]

;; -------------------------------------------------- --------------------- ;; interfaz, con contrato ( proporcionar/contratar [adaptive-simpson ( ->i (( f ( -> real? real? )) ( L real? ) ( R ( L ) ( y/c real? ( >/c L ) ))) ( #:épsilon ( ε real? )) ( r real? )) ] )                    ;; -------------------------------------------------- --------------------- ;; implementación ( definir ( adaptive-simpson f L R #:epsilon [ε . 000000001] ) ( definir f@L ( f L )) ( definir f@R ( f R )) ( definir-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-values ​​( 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 ( define é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 medio intervalo (1 func 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" integrándose de izquierda a derecha y ajustando el ancho del intervalo según sea necesario. ^[2]
• El algoritmo 103 de Kuncir (1962) es el integrador recursivo, bisector y adaptativo 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 sentencia goto . Protege contra el desbordamiento de anchos de intervalo (bucle BB) y cancela 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 ⁽²⁾ 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 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, se utiliza para la terminación, por lo que en su lugar se utiliza una mejora de 1963 . ^{[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 con la bisección para reducir los costos computacionales (Modificaciones 1+2, coincidiendo con el integrador Kuncir) y mejora las estimaciones de error de McKeeman 1962/63 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 ε al mínimo permitido por la precisión actual y devolver el nuevo error. ^[3]

Notas

1. El 4i original solo menciona elevar E'. Sin embargo, un texto posterior mencionó que se puede bajar en grandes pasos.
2. Esto probablemente también se aplica a los desbordamientos del ancho de intervalo/tolerancia en el caso simplista.

Bibliografía

1. GF Kuncir (1962), "Algoritmo 103: integrador de reglas de Simpson", Comunicaciones de la ACM , 5 (6): 347, doi : 10.1145/367766.368179
2. Para un integrador adaptativo no recursivo anterior que recuerda más a los solucionadores de EDO , consulte S. Henriksson (1961), "Contribución n.° 2: integración numérica de Simpson con longitud de paso variable", BIT Numerical Mathematics , 1 : 290
3. 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
4. Berrill, Mark A. "RayTrace-miniapp: src/AtomicModel/interp.hpp · de5e8229bccf60ae5c1c5bab14f861dc0326d5f9". ORNL GitLab .
5. Felleisen, Matías. "[raqueta] integración adaptativa de Simpson". "Usuarios" de la lista de correo de Racket . Consultado el 26 de septiembre de 2018 .
6. McKeeman, William Marshall (1 de diciembre de 1962). "Algoritmo 145: integración numérica adaptativa según la regla de Simpson". Comunicaciones de la ACM . 5 (12): 604. doi : 10.1145/355580.369102 .

enlaces externos

• Módulo para la regla de Simpson adaptativa