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

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

введите описание изображения здесь

Здесь W – функция всех переменных. Это известная функция, для которой я могу определить функцию python.

Вопрос программирования. Как мне получить scipy для интеграции этого выражения? Я думал о соединении двух тройных квадов ( scipy.integrate.tplquad ) вместе, но я беспокоюсь о производительности и точности. Есть ли более высокий размерный интегратор в scipy , который может обрабатывать произвольное число вложенных интегралов? Если нет, то каков наилучший способ сделать это?

3 Solutions collect form web for “Интегрирование многомерного интеграла в scipy”

С таким высокомерным интегралом, как это, методы 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 % 

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

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

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

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

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

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

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

Сначала сказать, что я не так хорош в математике, поэтому, пожалуйста, будьте добрыми. Во всяком случае, вот моя попытка:
Обратите внимание, что в вашем вопросе есть 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]

 \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*} 

Вы можете немного поиграть в свой вопрос;)

Python - лучший язык программирования в мире.