L'intégration multidimensionnelle intégrale dans scipy

Motivation: j'ai un multidimensionnelle intégrale, qui pour être complet, j'ai reproduit ci-dessous. Il s'agit du calcul du second coefficient de viriel quand il y a une importante anisotropie:

L'intégration multidimensionnelle intégrale dans scipy

Ici W est une fonction de toutes les variables. C'est une fonction connue, qui, je peux définir une fonction python.

Programmation Question: Comment puis-je obtenir scipy à intégrer cette expression? Je pensais d'enchaînement de deux triples quadruples ( scipy.s'intégrer.tplquad ), mais je suis inquiet au sujet de la performance et de la précision. Est-il la plus grande dimension de l'intégration sur scipy, que l'on peut gérer un nombre arbitraire de imbriqués les intégrales? Si non, quelle est la meilleure façon de le faire?

Vous avez peut-être mieux essayer de Sympy.

OriginalL'auteur Hooked | 2012-12-28

  1. 21

    Avec une plus-dimensionnelle partie intégrante de ce genre, les méthodes de monte carlo sont souvent une technique utile - ils convergent sur la réponse à l'inverse de la racine carrée du nombre de fonction des évaluations, ce qui est mieux pour dimension supérieure, alors vous aurez généralement obtenir de même assez sophistiqué méthodes adaptatives (sauf si vous savez quelque chose de très spécifique au sujet de votre intégrande - symétries qui peuvent être exploités, etc.)

    La mcint package effectue de monte-carlo, de l'intégration: en cours d'exécution avec un non-trivial W qui est néanmoins intégrable si nous connaissons la réponse, nous obtenons (notez que j'ai tronqué r à partir de [0,1); vous aurez à faire une sorte de journal de transformation ou de quelque chose pour obtenir la demi-illimitée de domaine en quelque chose de souple, pour la plupart des intégrateurs numériques):

    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 " "

    De course donne

    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 %

    Vous pourriez grandement de la vitesse par la vectorisation de la génération de nombres aléatoires, etc.

    Bien sûr, vous pouvez la chaîne de la triple intégrales comme vous le suggérez:

    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, "%"

    qui est lent, mais donne de très bons résultats pour ce cas simple. Ce qui est mieux est d'aller pour en venir à la complexité de votre W est et ce que vos exigences de précision sont. Simple (rapide à évaluer) W avec une grande précision va vous pousser à ce genre de méthode; compliqué (de lent à évaluer) W modérées avec des exigences de précision va vous pousser vers MC techniques.

    Result =  1632.10498552  estimated error =  3.59054059995e-11
Known result =  1632.10498552  error =  4.54747350886e-13  =  2.7862628625e-14 %
    Merci! Je vais jeter un oeil à mcint et voir si cela fonctionne mieux que mon ad-hoc MC méthode que j'ai maintenant.
    +1 Bon d'apprendre de nouvelles choses 🙂

    OriginalL'auteur Jonathan Dursi

  2. 6

    Je vais juste faire quelques commentaires d'ordre général sur comment bien faire ce genre d'intégrale, mais ce conseil n'est pas spécifique à scipy (trop long pour un commentaire, même s'il n'est pas une réponse).

    Je ne sais pas votre cas d'utilisation, c'est à dire si vous êtes satisfait avec une "bonne" réponse avec quelques chiffres de précision qui pourrait être obtenu de manière simple à l'aide de Monte-Carlo comme indiqué dans Jonathan Dursi de réponse, ou si vous voulez vraiment pousser la précision numérique autant que possible.

    Que j'ai effectué d'analyse de Monte-Carlo et en quadrature calculs de coefficients du viriel moi-même. Si vous voulez faire les intégrales avec précision, alors il ya quelques choses que vous devriez faire:

    1. Tenter de faire le plus grand nombre des intégrales exactement que possible; il se peut bien que l'intégration de certains de vos coordonnées est assez simple.

    2. Envisager de transformer vos variables de l'intégration, de sorte que la fonction à intégrer est aussi lisse que possible. (Cela aide pour les deux Monte-Carlo et en quadrature).

    3. De Monte-Carlo, l'utilisation de l'importance de l'échantillonnage pour une meilleure convergence.

    4. Pour la quadrature, avec 7 intégrales c'est peut-être possible d'obtenir vraiment une convergence rapide à l'aide de tanh-sinh quadrature. Si vous pouvez l'obtenir jusqu'à 5 intégrales, alors vous devriez être en mesure d'obtenir jusqu'à 10 chiffres de précision pour votre intégrale. Je recommande fortement mathtool /ARPREC à cet effet, disponible a partir de David Bailey page d'accueil: http://www.davidhbailey.com/

    Merci pour l'entrée. Avez-vous l'esprit sur l'élaboration de #2? A priori comment puis-je savoir ce qu'est un bon transformer? Puisque vous avez fait ce type de calculs avant, toute entrée supplémentaire serait appréciée.

    OriginalL'auteur Nathan

  3. 1

    D'abord à vous dire que je ne suis pas bon en math, donc s'il vous plaît être gentil. De toute façon, voici mon essai:

    Notez que dans votre question, il y a 6 variables, mais 7 intégrales!?

    Dans Python à l'aide de 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)))

    et voici le résultat: [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*}

    Vous pouvez jouer un peu plus pour votre question 😉

    Qui dirait qu'elle est en train de faire calcul symbolique, c'est à dire votre W est une fonction linéaire des variables d'entrée, d'où le résultat exact. Pour moi W est non-linéaire et non exprimable comme une fonction mathématique, mais comme le résultat d'un autre calcul (donc défini comme une fonction python). Il est exact que je ne devrait avoir que 6 intégrales, je dois avoir emporté TeXing.

    OriginalL'auteur Developer

  4. 1

    Jonathan Dursi a fait une très bonne réponse. Je vais juste ajouter à sa réponse.

    Maintenant scipy.integrate a une fonction nommée nquad que l'on peut effectuer un multi-dimensionnelle partie intégrante sans tracas. Voir ce lien pour plus d'informations. Ci-dessous, nous calculons l'intégrale à l'aide de nquad avec Jonathan exemple:

    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)

    Le résultat est plus précis que le enchaînés tplquad:

    >>> print(diff)
0.0

    OriginalL'auteur Dongchen Zou