Интегрирование многомерного интеграла в scipy

Ответ 1

С таким высокомерным интегралом, как это, методы monte carlo часто являются полезной техникой - они сходятся на ответе как обратный квадратный корень от числа оценок функций, что лучше для более высокого измерения, тогда вы обычно получаете из-за достаточно сложных адаптивных методов (если вы не знаете что-то очень специфическое в отношении вашего подынтегрального выражения - симметрии, которые можно использовать и т.д.)

Пакет mcint выполняет интегральную интеграцию carlo: работает с нетривиальным W, который тем не менее интегрируется, поэтому мы знаем, ответим нам (обратите внимание, что я усекал r от [0,1]; вам придется сделать какой-то лог-преобразование или что-то, чтобы получить этот полубесконечный домен в нечто приемлемое для большинства числовых интеграторов):

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

Запуск дает

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 %

Вы могли бы значительно ускорить это путем векторизации генерации случайных чисел и т.д.

Конечно, вы можете связать тройные интегралы, как вы предлагаете:

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

который является медленным, но дает очень хорошие результаты для этого простого случая. Который лучше будет свести к тому, насколько сложна ваша W и каковы ваши требования к точности. Простая (быстрая оценка) W с высокой точностью подтолкнет вас к такому методу; сложный (медленно оцениваемый) W с умеренными требованиями к точности подтолкнет вас к методам MC.

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

Ответ 2

Я просто сделаю пару общих замечаний о том, как точно сделать такой интеграл, но этот совет не является специфичным для scipy (слишком длинный для комментария, хотя это не ответ).

Я не знаю вашего случая использования, то есть, удовлетворены ли вы "хорошим" ответом с помощью нескольких цифр точности, которые можно было бы получить прямо с помощью Монте-Карло, как указано в ответе Джонатана Дурси, или вы действительно хотите насколько возможно, нажмите цифровую точность.

Я сам выполнил аналитические, Монте-Карло и квадратурные вычисления вириальных коэффициентов. Если вы хотите точно выполнить интегралы, то вам нужно сделать несколько вещей:

• Попытка выполнить как можно больше интегралов так, насколько это возможно; вполне возможно, что интеграция в некоторые из ваших координат довольно проста.

• Попробуйте преобразовать переменные интеграции так, чтобы подынтегральное выражение было максимально гладким. (Это помогает как для Монте-Карло, так и для квадратуры).

• Для Монте-Карло используйте выборку важности для лучшей конвергенции.

  <квадратичный, с 7 интегралами, возможно, можно получить действительно быструю сходимость, используя квадратичную формулу tanh-sinh. Если вы можете получить его до 5 интегралов, тогда вы сможете получить 10 цифр точности для своего интеграла. Я очень рекомендую mathtool/ARPREC для этой цели, доступный на главной странице David Bailey:

Ответ 3

Прежде всего сказать, что я не так хорош в математике, поэтому, пожалуйста, будьте добрыми. Во всяком случае, вот моя попытка:
Обратите внимание, что в вашем вопросе есть 6 переменные, но 7 интегралы!?
В Python с помощью 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)))  

и вот результат: [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*}

Вы можете сыграть немного больше для своего вопроса;)