Integrazione di un integrale multidimensionale in scipy

Question

Motivazione: Ho un integrale multidimensionale, che per completezza ho riprodotto qui sotto. Viene dal calcolo del secondo coefficiente del viriale quando v'è una significativa anisotropia:

Qui W è una funzione di tutte le variabili. È una funzione nota, una per la quale posso definire una funzione python.

Domanda di programmazione: Come si ottiene scipy per integrare questa espressione? Stavo pensando di concatenare due triple quad (scipy.integrate.tplquad) insieme, ma sono preoccupato per le prestazioni e la precisione. Esiste un integratore di dimensioni superiori in scipy, uno in grado di gestire un numero arbitrario di integrali nidificati? In caso contrario, qual è il modo migliore per farlo?

Answer 1

Con un superiore che integrali dimensionali come questo, i metodi di monte carlo sono spesso una tecnica utile - convergono sulla risposta come la radice quadrata inversa del numero di valutazioni di funzione, che è meglio per una dimensione più alta quindi generalmente uscire anche abbastanza sofisticati metodi adattivi (a meno che non si sa qualcosa di molto specifiche sul integrando - simmetrie che possono essere sfruttate, etc.)

Il pacchetto mcint esegue un'integrazione Monte Carlo: in esecuzione con un non banale W che è comunque integrabile, quindi conosciamo la risposta che otteniamo (notare che ho troncato r per essere da [0,1); dovrete fare una sorta di registro di trasformare o qualcosa per ottenere quel dominio semi-illimitata in qualcosa trattabili per la maggior parte integratori numerici):

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

Correre dà

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 % 

Si potrebbe accelerare notevolmente questo da vettorizzare la generazione di numeri casuali, ecc

Naturalmente, è possibile concatenare le integrali tripli come lei suggerisce:

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

che è lento ma fornisce ottimi risultati per questo semplice caso. Qual è il modo migliore per capire quanto è complicato il tuo W e quali sono i tuoi requisiti di precisione. W semplice (veloce da valutare) con alta precisione ti spingerà a questo tipo di metodo; W complicato (lento a valutare) con requisiti di precisione moderati vi spingerà verso le tecniche MC.

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

Answer 2

In primo luogo per dire che non sono così bravo in matematica, quindi per favore sii gentile. Comunque, ecco il mio tentativo:
Nota che nella tua domanda ci sono variabili ma integrali !?
In Python utilizzando 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))) 

e qui è il risultato: [codice 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*} 

Puoi giocare un po 'di più per la tua domanda;)

Answer 3

mi limiterò a fare un paio di osservazioni di carattere generale su come farlo con precisione questo tipo di integrale, ma questo consiglio non è specifico per SciPy (troppo lungo per un commento, anche se non è una risposta) .

Non conosco il tuo caso d'uso, vale a dire se sei soddisfatto di una risposta "buona" con poche cifre di precisione che possono essere ottenute direttamente usando Monte Carlo come delineato nella risposta di Jonathan Dursi, o se vuoi davvero per spingere la precisione numerica il più lontano possibile.

Ho eseguito io stesso calcoli analitici, Monte Carlo e in quadratura di coefficienti virali. Se si vuole fare gli integrali con precisione, poi ci sono alcune cose che si dovrebbe fare:

1. tenta di eseguire il maggior numero di integrali esattamente come possibili; potrebbe essere che l'integrazione in alcune delle tue coordinate sia abbastanza semplice.

2. Considerare la trasformazione delle variabili di integrazione in modo che l'integrando sia il più fluido possibile. (Questo aiuta sia per Monte Carlo che per la quadratura).

3. Per Monte Carlo, utilizzare il campionamento di importanza per la migliore convergenza.

4. Per la quadratura, con 7 integrali, è possibile ottenere una convergenza davvero rapida utilizzando la quadratura tanh-sinh. Se riesci a ottenere 5 integrali, dovresti essere in grado di ottenere 10 cifre di precisione per il tuo integrale. Mi raccomando mathtool/ARPREC per questo scopo, disponibile presso l'homepage di David Bailey: http://www.davidhbailey.com/