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

Tengo una integral impropia (bidimensional)

$$I=\int_A \frac{W(x,y)}{F(x,y)}\,\mbox{d}x\mbox{d}y$$

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 ≠ $A$$x=[-1,1]$$y=[-1,1]$$F(x,y)>0$$F$$W$$W \ne 0$en 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)$$y$$x$$I$$W(x,y)$

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.

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á $\int_{[0,1]^n} f(x)\; d^nx$$n$$M$$M^n$$k$$N=(kM)^n$$k$$(2k-1)$$e = {\cal O}(h^5)={\cal O}(M^{-(2k-1)})$

$$e = {\cal O}(N^{-(2k-1)/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 / $$e = {\cal O}(N^{-1/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 $k\le 8$$n=30$$M=1$$N=8^{30}$

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.