python - indefinite - scipy help
Integrando una integral multidimensional en scipy (3)
Motivación: Tengo una integral multidimensional, que para completar he reproducido a continuación. Viene del cálculo del segundo coeficiente virial cuando hay anisotropía significativa:
Aquí W es una función de todas las variables. Es una función conocida, una para la que puedo definir una función python.
Pregunta de programación: ¿Cómo puedo obtener scipy para integrar esta expresión? Estaba pensando en encadenar dos cuádruples cuádruples ( scipy.integrate.tplquad ) juntos, pero me preocupa el rendimiento y la precisión. ¿Hay un integrador dimensional más alto 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 cual es mejor para una dimensión más alta, por lo general saldrán de métodos de adaptación bastante sofisticados (a menos que sepa algo muy específico sobre 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, no obstante, es integrable, así que sabemos la respuesta que obtenemos (nótese que he truncado r para ser desde [0,1); Tendrás que hacer algún tipo de transformación de logs o algo para que ese dominio semi-ilimitado se convierta 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ías acelerarlo enormemente al vectorizar la generación de números aleatorios, etc.
Por supuesto, puede encadenar las integrales triples como sugiere:
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. W simple (rápido de evaluar) con alta precisión lo empujará a este tipo de método; complicado (lento de evaluar) W con requisitos de precisión moderada lo empujará hacia las técnicas MC.
Result = 1632.10498552 estimated error = 3.59054059995e-11
Known result = 1632.10498552 error = 4.54747350886e-13 = 2.7862628625e-14 %
Primero, decir que no soy tan bueno en matemáticas, así que por favor sé amable. De todos modos, aquí está mi intento:
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: [código LateX]
/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*}
Puede jugar un poco más por su pregunta;)
Haré un par de comentarios generales acerca de cómo hacer exactamente este tipo de integral, pero este consejo no es específico para scipy (demasiado tiempo para un comentario, aunque no sea una respuesta).
No conozco su caso de uso, es decir, si está satisfecho con una "buena" respuesta con unos pocos dígitos de precisión que podría obtenerse directamente usando Monte Carlo como se describe en la respuesta de Jonathan Dursi, o si realmente desea presionar el numérico precisión en la medida de lo posible.
Realicé cálculos analíticos, de Monte Carlo y en cuadratura de coeficientes viriales. Si quiere hacer las integrales con precisión, entonces hay algunas cosas que debe hacer:
Intente realizar tantas integrales exactamente como sea posible; bien puede ser que la integración en algunas de sus coordenadas sea bastante simple.
Considere transformar sus variables de integración para que el integrando sea lo más sencillo posible. (Esto ayuda tanto a Monte Carlo como a la cuadratura).
Para Monte Carlo, use el muestreo de importancia para la mejor convergencia.
Para la cuadratura, con 7 integrales, es posible lograr una convergencia realmente rápida usando la cuadratura tanh-sinh. Si puede bajarlo a 5 integrales, entonces debería poder obtener 10 dígitos de precisión para su integral. Recomiendo mathtool / ARPREC para este propósito, disponible en la página principal de David Bailey: http://www.davidhbailey.com/