Integración numérica de integral multidimensional con límites conocidos

12

Tengo una integral impropia (bidimensional)

I=AW(x,y)F(x,y)dxdy

donde el dominio de integración es menor que , pero restringido aún más por . Como y son suaves yx = [ - 1 , 1 ] y = [ - 1 , 1 ] F ( x , y ) > 0 F W W 0Ax=[1,1]y=[1,1]F(x,y)>0FWW0en los límites, la relación posterior implica que el integrando puede ser singular en los límites. Sin embargo, el integrando es finito. Yo, hasta ahora, calculo esta integral con integración numérica anidada. Esto es exitoso pero lento. Busco un método más apropiado (más rápido) para abordar la integral, tal vez un método de Montecarlo. Pero necesito uno que no ponga puntos en el límite del dominio no cúbico A y tome el límite de la integral impropia correctamente. ¿Puede una transformación integral ayudar a esta expresión general? Tenga en cuenta que puedo resolver para como una función de e incluso calcular para algunas funciones especiales de peso .y x I W ( x , y )F(x,y)yxIW(x,y)

Highsciguy
fuente
¿Puedes ser un poco más específico sobre los métodos que has usado hasta ahora? ¿Qué rutinas específicas has usado de forma anidada? Además, ¿está dentro de , es decir, ¿las raíces de solo están en el límite? A F ( x , y )F(x,y)0AF(x,y)
Pedro
El algoritmo GSL QAGS: gnu.org/software/gsl/manual/html_node/… . ¡Gracias por las ediciones (no vi la opción de escribir ecuaciones)!
highsciguy

Respuestas:

7

Descargo de responsabilidad: escribí mi tesis doctoral sobre la cuadratura adaptativa, por lo que esta respuesta estará sesgada hacia mi propio trabajo.

QAGS de GSL es el antiguo integrador QUADPACK , y no es completamente robusto, especialmente en presencia de singularidades. Esto generalmente lleva a los usuarios a solicitar muchos más dígitos de precisión de los que realmente necesitan, lo que hace que la integración sea bastante costosa.

Si está utilizando GSL, es posible que desee probar mi propio código, CQUAD , que se describe en este documento . Está diseñado para hacer frente a las singularidades, tanto en los bordes del intervalo como dentro del dominio. Tenga en cuenta que la estimación del error es bastante sólida, por lo que solo solicite tantos dígitos como realmente necesite.

Con respecto a la integración de Monte-Carlo, depende del tipo de precisión que esté buscando. Tampoco estoy muy seguro de qué tan bien funcionará cerca de las singularidades.

Pedro
fuente
Definitivamente voy a echar un vistazo a esto porque será más sencillo implementarlo. De hecho, experimenté que la rutina QAGS no era súper estable para este problema.
highsciguy
¿Hay alguna manera de influir en la aparición de 'GSL_EDIVERGE'? Parece aparecer para algunos parámetros.
highsciguy
@highsciguy: el algoritmo devuelve GSL_EDIVERGE cuando cree que la integral no es finita. Si pudiera darme un ejemplo para el que falla, podría echarle un vistazo más de cerca.
Pedro
Es un poco difícil aislar una rutina simple, ya que está incrustada en un código genérico para integrales n-dimensinales. Veré ... Pero para y fijo, 1 / sqrt (F (x, y)) debería comportarse como 1 / sqrt (x) cuando x se acerca a los ceros de F (x, y) ya que F (x, y) entonces se puede escribir como un polinomio en x. Pero podría ser que el comportamiento 1 / sqrt (x) comience tarde. También podría ser que la precisión numérica del integrando no sea demasiado buena.
highsciguy
1
@highsciguy: Sí, esta es una mala idea. La mayoría de las reglas de cuadratura suponen que el integrando tiene cierto grado de suavidad, y si lo establece en cero a partir de algún punto arbitrario, está introduciendo una discontinuidad. ¡Obtendrá resultados mucho mejores si usa el intervalo real!
Pedro
5

Los métodos de Monte Carlo en general no pueden competir con la cuadratura adaptativa a menos que tenga una integral de alta dimensión donde no pueda permitirse la explosión combinatoria de los puntos de la cuadratura con la dimensión.

La razón es relativamente fácil de entender. Tomemos, por ejemplo, solo donde es la dimensión del problema. Digamos, por simplicidad, que subdivides cada dimensión en intervalos intermedios, es decir, obtienes celdas de hipercubos en total. Supongamos además que usa una fórmula de Gauss con puntos de Gauss, solo como un ejemplo. Entonces tiene puntos de cuadratura en total, y debido a que los puntos Gauss le proporcionan precisión de orden , , su precisión general en función de los puntos de evaluación será [0,1]nf(x)dnxnMMnkN=(kM)nk(2k1)e=O(h5)=O(M(2k1))

e=O(N(2k1)/n).
Por otro lado, los métodos de Monte Carlo generalmente solo le proporcionan convergencia de errores como que es peor que para cualquier fórmula de Gauss con al menos puntos por intervalo. La razón es relativamente simple de entender: la cuadratura de Gauss elige los puntos de interpolación de una manera inteligente, Monte Carlo no lo es. No se puede esperar nada útil de esto último. (Por supuesto, hay situaciones en las que la cuadratura gaussiana es difícil: por ejemplo, en su caso donde el dominio de integración tiene una forma irregular; pero en ese caso, es probable que aún sea mejor hacer una integración adaptativa o similar).k > n / 4 + 1 / 2
e=O(N1/2)
k>n/4+1/2

Ahora, hay problemas prácticos (de estabilidad) con la integración con más de, digamos, 8 o 10 puntos por intervalo. Entonces, si quieres , entonces no puedes ir más allá de . Por otro lado, en ese caso, incluso elegir un solo intervalo por dirección ( ) produce puntos de integración, mucho más de lo que podría en una evaluación de por vida. En otras palabras, siempre que pueda evaluar suficientes puntos de integración, la cuadratura en las subdivisiones de su dominio de integración es siempre el enfoque más eficiente. En los casos en los que tiene una integral de alta dimensión para la que no puede evaluar los puntos de integración en una sola subdivisión, las personas usan los métodos de Monte Carlo a pesar de su peor orden de convergencia.n = 30 M = 1 N = 8 30k8n=30M=1N=830

Wolfgang Bangerth
fuente
1

Pruebe una cuadratura doble exponencial anidada (consulte las implementaciones de Ooura ). Esta técnica utiliza una transformación variable que hace que el integrando transformado se comporte muy suavemente en los límites y es muy eficiente para manejar singularidades en el límite. También hay una muy buena lista de referencias en la cuadratura DE en su sitio web.

GertVdE
fuente