在 scipy 中积分多维积分 [英] Integrating a multidimensional integral in scipy

问题描述

动机:我有一个多维积分，为了完整起见，我在下面复制了它.它来自于存在显着各向异性时的第二维里系数的计算:

Motivation: I have a multidimensional integral, which for completeness I have reproduced below. It comes from the computation of the second virial coefficient when there is significant anisotropy:

这里 W 是所有变量的函数.这是一个已知的函数，我可以为它定义一个 python 函数.

Here W is a function of all the variables. It is a known function, one which I can define a python function for.

编程问题: 如何让 scipy 整合这个表达式?我正在考虑链接两个三重四边形(scipy.integrate.tplquad) 在一起，但我担心性能和准确性.scipy 中是否有更高维的积分器，可以处理任意数量的嵌套积分?如果没有，最好的方法是什么?

Programming Question: How do I get scipy to integrate this expression? I was thinking of chaining two triple quads (scipy.integrate.tplquad) together, but I'm worried about performance and accuracy. Is there a higher dimensional integrator in scipy, one that can handle an arbitrary number of nested integrals? If not, what is the best way to do this?

推荐答案

对于像这样的高维积分，蒙特卡罗方法通常是一种有用的技术——它们收敛于函数数量的平方根的倒数评估，这对于更高维度更好，那么您通常会摆脱甚至相当复杂的自适应方法(除非您对被积函数非常了解 - 可以利用的对称性等)

With a higher-dimensional integral like this, monte carlo methods are often a useful technique - they converge on the answer as the inverse square root of the number of function evaluations, which is better for higher dimension then you'll generally get out of even fairly sophisticated adaptive methods (unless you know something very specific about your integrand - symmetries that can be exploited, etc.)

mcint 包执行蒙特卡罗集成:使用非平凡的<代码运行>W 仍然是可积的，所以我们知道我们得到的答案(注意我已经截断了 r 来自 [0,1)；您必须进行某种对数变换或其他操作，才能将该半无界域转换为大多数数值积分器易于处理的域):

The mcint package performs a monte carlo integration: running with a non-trivial W that is nonetheless integrable so we know the answer we get (note that I've truncated r to be from [0,1); you'll have to do some sort of log transform or something to get that semi-unbounded domain into something tractable for most numerical integrators):

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 %

您可以通过向量化随机数生成等来大大加快速度.

You could greatly speed this up by vectorizing the random number generation, etc.

当然，您可以按照您的建议链接三重积分:

Of course, you can chain the triple integrals as you suggest:

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 技术.

which is slow but gives very good results for this simple case. Which is better is going to come down to how complicated your W is and what your accuracy requirements are. Simple (fast to evaluate) W with high accuracy will push you to this sort of method; complicated (slow to evaluate) W with moderate accuracy requirements will push you towards MC techniques.

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

这篇关于在 scipy 中积分多维积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！