Integrando una integral multidimensional en scipy.

Motivación: tengo una integral multidimensional, que para completarla, he reproducido a continuación. Proviene del cálculo del segundo coeficiente virial cuando existe una anisotropía significativa:

introduzca la descripción de la imagen aquí

Aquí W es una función de todas las variables. Es una función conocida, una para la que puedo definir una función de python.

Pregunta de progtwigción: ¿Cómo obtengo scipy para integrar esta expresión? Estaba pensando en encadenar dos quads triples ( scipy.integrate.tplquad ), pero me preocupa el rendimiento y la precisión. ¿Hay un integrador dimensional superior en scipy , uno que pueda manejar un número arbitrario de integrales anidadas? Si no, ¿cuál es la mejor manera de hacer esto?

Con una integral de mayor dimensión como esta, los métodos de monte carlo son a menudo una técnica útil: convergen en la respuesta como la raíz cuadrada inversa del número de evaluaciones de funciones, lo que es mejor para una dimensión más alta que, por lo general, saldrá de la par. Métodos adaptativos bastante sofisticados (a menos que sepa algo muy específico acerca de su integrando: simetrías que pueden ser explotadas, etc.)

El paquete mcint realiza una integración de monte carlo: se ejecuta con un W no trivial que, sin embargo, es integrable, por lo que sabemos la respuesta que obtenemos (tenga en cuenta que he truncado r para ser de [0,1); Tendrá que hacer algún tipo de transformación de registro o algo para convertir ese dominio semi-ilimitado en algo manejable para la mayoría de los integradores numéricos):

 import mcint import random import math def w(r, theta, phi, alpha, beta, gamma): return(-math.log(theta * beta)) def integrand(x): r = x[0] theta = x[1] alpha = x[2] beta = x[3] gamma = x[4] phi = x[5] k = 1. T = 1. ww = w(r, theta, phi, alpha, beta, gamma) return (math.exp(-ww/(k*T)) - 1.)*r*r*math.sin(beta)*math.sin(theta) def sampler(): while True: r = random.uniform(0.,1.) theta = random.uniform(0.,2.*math.pi) alpha = random.uniform(0.,2.*math.pi) beta = random.uniform(0.,2.*math.pi) gamma = random.uniform(0.,2.*math.pi) phi = random.uniform(0.,math.pi) yield (r, theta, alpha, beta, gamma, phi) domainsize = math.pow(2*math.pi,4)*math.pi*1 expected = 16*math.pow(math.pi,5)/3. for nmc in [1000, 10000, 100000, 1000000, 10000000, 100000000]: random.seed(1) result, error = mcint.integrate(integrand, sampler(), measure=domainsize, n=nmc) diff = abs(result - expected) print "Using n = ", nmc print "Result = ", result, "estimated error = ", error print "Known result = ", expected, " error = ", diff, " = ", 100.*diff/expected, "%" print " " 

Correr da

 Using n = 1000 Result = 1654.19633236 estimated error = 399.360391622 Known result = 1632.10498552 error = 22.0913468345 = 1.35354937522 % Using n = 10000 Result = 1634.88583778 estimated error = 128.824988953 Known result = 1632.10498552 error = 2.78085225405 = 0.170384397984 % Using n = 100000 Result = 1646.72936 estimated error = 41.3384733174 Known result = 1632.10498552 error = 14.6243744747 = 0.8960437352 % Using n = 1000000 Result = 1640.67189792 estimated error = 13.0282663003 Known result = 1632.10498552 error = 8.56691239895 = 0.524899591322 % Using n = 10000000 Result = 1635.52135088 estimated error = 4.12131562436 Known result = 1632.10498552 error = 3.41636536248 = 0.209322647304 % Using n = 100000000 Result = 1631.5982799 estimated error = 1.30214644297 Known result = 1632.10498552 error = 0.506705620147 = 0.0310461413109 % 

Podría acelerar enormemente esto vectorizando la generación de números aleatorios, etc.

Por supuesto, puede encadenar las integrales triples como sugiera:

 import numpy import scipy.integrate import math def w(r, theta, phi, alpha, beta, gamma): return(-math.log(theta * beta)) def integrand(phi, alpha, gamma, r, theta, beta): ww = w(r, theta, phi, alpha, beta, gamma) k = 1. T = 1. return (math.exp(-ww/(k*T)) - 1.)*r*r*math.sin(beta)*math.sin(theta) # limits of integration def zero(x, y=0): return 0. def one(x, y=0): return 1. def pi(x, y=0): return math.pi def twopi(x, y=0): return 2.*math.pi # integrate over phi [0, Pi), alpha [0, 2 Pi), gamma [0, 2 Pi) def secondIntegrals(r, theta, beta): res, err = scipy.integrate.tplquad(integrand, 0., 2.*math.pi, zero, twopi, zero, pi, args=(r, theta, beta)) return res # integrate over r [0, 1), beta [0, 2 Pi), theta [0, 2 Pi) def integral(): return scipy.integrate.tplquad(secondIntegrals, 0., 2.*math.pi, zero, twopi, zero, one) expected = 16*math.pow(math.pi,5)/3. result, err = integral() diff = abs(result - expected) print "Result = ", result, " estimated error = ", err print "Known result = ", expected, " error = ", diff, " = ", 100.*diff/expected, "%" 

que es lento pero da muy buenos resultados para este caso simple. Lo que es mejor se reducirá a lo complicado que es su W y cuáles son sus requisitos de precisión. Simple (rápido de evaluar) W con alta precisión lo empujará a este tipo de método; La W complicada (lenta de evaluar) con requisitos de precisión moderada lo empujará hacia técnicas de MC.

 Result = 1632.10498552 estimated error = 3.59054059995e-11 Known result = 1632.10498552 error = 4.54747350886e-13 = 2.7862628625e-14 % 

Solo haré un par de comentarios generales sobre cómo hacer con precisión este tipo de integral, pero este consejo no es específico para scipy (demasiado largo para un comentario, aunque no sea una respuesta).

No conozco su caso de uso, es decir, si está satisfecho con una respuesta “buena” con unos pocos dígitos de precisión que podría obtenerse utilizando Monte Carlo como se describe en la respuesta de Jonathan Dursi, o si realmente quiere presionar el botón numérico. Precisión en la medida de lo posible.

He realizado cálculos analíticos, de Monte Carlo y en cuadratura de coeficientes viriales. Si quieres hacer las integrales con precisión, hay algunas cosas que debes hacer:

  1. Intente realizar tantas integrales como sea posible; bien puede ser que la integración en algunas de sus coordenadas sea bastante simple.

  2. Considere la posibilidad de transformar sus variables de integración para que el integrando sea lo más suave posible. (Esto ayuda tanto para Monte Carlo como para la cuadratura).

  3. Para Monte Carlo, use el muestreo de importancia para la mejor convergencia.

  4. Para la cuadratura, con 7 integrales puede ser posible obtener una convergencia realmente rápida utilizando la cuadratura de tanh-sinh. Si puede reducirlo a 5 integrales, entonces debería poder obtener 10 dígitos de precisión para su integral. Recomiendo altamente mathtool / ARPREC para este propósito, disponible en la página de inicio de David Bailey: http://www.davidhbailey.com/

Jonathan Dursi ha hecho una muy buena respuesta. Acabo de añadir a su respuesta.

Ahora scipy.integrate tiene una función llamada nquad que permite realizar una integral multidimensional sin problemas. Vea este enlace para más información. A continuación, calculamos la integral utilizando nquad con el ejemplo de Jonathan:

 from scipy import integrate import math import numpy as np def w(r, theta, phi, alpha, beta, gamma): return(-math.log(theta * beta)) def integrand(r, theta, phi, alpha, beta, gamma): ww = w(r, theta, phi, alpha, beta, gamma) k = 1. T = 1. return (math.exp(-ww/(k*T)) - 1.)*r*r*math.sin(beta)*math.sin(theta) result, error = integrate.nquad(integrand, [[0, 1], # r [0, 2 * math.pi], # theta [0, math.pi], # phi [0, 2 * math.pi], # alpha [0, 2 * math.pi], # beta [0, 2 * math.pi]]) # gamma expected = 16*math.pow(math.pi,5)/3 diff = abs(result - expected) 

El resultado es más preciso que el tplquad encadenado:

 >>> print(diff) 0.0 

Primero decir que no soy tan bueno en matemáticas, así que por favor sea amable. De todos modos, aquí está mi bash:
Tenga en cuenta que en su pregunta hay 6 variables pero 7 integrales.
En Python usando Sympy :

 >>> r,theta,phi,alpha,beta,gamma,W,k,T = symbols('r theta phi alpha beta gamma W k T') >>> W = r+theta+phi+alpha+beta+gamma >>> Integral((exp(-W/(k*T))-1)*r**2*sin(beta)*sin(theta),(r,(0,2*pi)),(theta,(0,pi)),(phi,(0,2*pi)),(alpha,(0,2*pi)),(beta,(0,pi)),(gamma,(0,pi))) >>> integrate((exp(-W)-1)*r**2*sin(beta)*sin(theta),(r,(0,2*pi)),(theta,(0,pi)),(phi,(0,2*pi)),(alpha,(0,2*pi)),(beta,(0,pi)),(gamma,(0,pi))) 

Y aquí está el resultado: [LateX code]

 \begin{equation*}- \frac{128}{3} \pi^{6} - \frac{\pi^{2}}{e^{2 \pi}} - \frac{\pi}{e^{2 \pi}} - \frac{2}{e^{2 \pi}} - \frac{\pi^{2}}{e^{3 \pi}} - \frac{\pi}{e^{3 \pi}} - \frac{2}{e^{3 \pi}} - 3 \frac{\pi^{2}}{e^{6 \pi}} - 3 \frac{\pi}{e^{6 \pi}} - \frac{2}{e^{6 \pi}} - 3 \frac{\pi^{2}}{e^{7 \pi}} - 3 \frac{\pi}{e^{7 \pi}} - \frac{2}{e^{7 \pi}} + \frac{1}{2 e^{9 \pi}} + \frac{\pi}{e^{9 \pi}} + \frac{\pi^{2}}{e^{9 \pi}} + \frac{1}{2 e^{8 \pi}} + \frac{\pi}{e^{8 \pi}} + \frac{\pi^{2}}{e^{8 \pi}} + \frac{3}{e^{5 \pi}} + 3 \frac{\pi}{e^{5 \pi}} + 3 \frac{\pi^{2}}{e^{5 \pi}} + \frac{3}{e^{4 \pi}} + 3 \frac{\pi}{e^{4 \pi}} + 3 \frac{\pi^{2}}{e^{4 \pi}} + \frac{1}{2 e^{\pi}} + \frac{1}{2}\end{equation*} 

Puedes jugar un poco más para tu pregunta;)