在scipy中集成多维积分

Question

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

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

编程问题:如何scipy集成此表达式？我想把两个三重四边形链接scipy.integrate.tplquad在一起,但我担心性能和准确性.是否有更高维的积分器scipy,可以处理任意数量的嵌套积分？如果没有,最好的方法是什么？

Answer 1

对于像这样的高维积分,monte carlo方法通常是一种有用的技术 - 它们作为函数评估数的倒数平方根收敛于答案,这对于更高维度更好,那么你通常会得到更多相当复杂的自适应方法(除非你知道关于你的被积函数的一些非常具体的东西 - 可以被利用的对称性等)

该mcint包执行蒙特卡洛积分:用不平凡的运行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 %

Answer 2

Jonathan Dursi 给出了很好的答案。我只会补充他的回答。

现在scipy.integrate有一个名为的函数nquad，可以轻松执行多维积分。有关更多信息，请参阅此链接。下面我们使用nquad乔纳森的例子计算积分：

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)

结果比链式更准确tplquad：

>>> print(diff)
0.0

Answer 3

我将就如何准确地完成这种积分做一些一般性的评论,但这个建议并不是针对scipy的(对于评论来说太长了,即使它不是答案).

我不知道你的用例,也就是说你是否满意于一个准确的几位数的"好"答案,可以使用Jonathan Dursi的答案中概述的蒙特卡罗直接获得,或者你是否真的想推动数字准确性尽可能高.

我自己进行了维数系数的解析,蒙特卡罗和正交计算.如果你想准确地进行积分,那么你应该做一些事情:

1. 尝试尽可能多地执行积分; 很可能你的一些坐标中的集成非常简单.

2. 考虑转换集成变量,以使被积函数尽可能平滑.(这有助于蒙特卡罗和正交).

3. 对于蒙特卡罗,使用重要性采样以获得最佳收敛.

4. 对于正交,使用7个积分,可以使用tanh-sinh积分实现真正快速的收敛.如果你可以将它归结为5个积分,那么你应该能够得到10个精度的数字作为积分.我强烈建议使用mathtool/ARPREC,David David Bailey的主页:http://www.davidhbailey.com/