在枕中积分多维积分

Question

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

​

​

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

编程问题：如何让scipy集成这个表达式？我想把两个三重四角(scipy.integrate.tplquad)链接在一起，但我担心性能和准确性。在scipy中是否有一个高维积分器，它可以处理任意数量的嵌套积分？如果没有，那么最好的方法是什么呢？

Answer 1

使用像这样的高维积分，monte方法通常是一种有用的技术--它们以函数求值数的倒数平方根的形式收敛到答案上，这对高维的计算更好，这样你通常就会摆脱即使是相当复杂的自适应方法(除非你知道一些关于积分对称性的非常具体的东西，可以利用，等等)。

mcint包执行monte集成:运行一个可积的非平凡的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

乔纳森·杜西给出了一个很好的答案。我只想补充一下他的答案。

现在，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

我只想就如何准确地进行这种积分提出几点一般性意见，但这个建议并不是专门针对for的(对评论来说太长了，尽管它不是一个答案)。

我不知道你的用例，即你是否满意一个“好的”答案的几位数的准确性，可以直接获得蒙特卡洛在乔纳森杜西的答案概述，或你是否真的想推动数字的准确性尽可能。

我自己也做过分析，蒙特卡罗和正交计算。如果您想要精确地进行积分，那么您应该做一些事情：

1. 尝试尽可能准确地执行多个积分；很可能在您的一些坐标中的集成非常简单。
2. 考虑转换您的积分变量，使integrand尽可能平滑。(这对蒙特卡罗和正交都有帮助)。
3. 对于蒙特卡罗，使用重要抽样进行最佳收敛。
4. 对于求积，对于7个积分，使用tanh-sinh求积可以得到非常快的收敛。如果你能把它降到5积分，那么你的积分就能得到10位数的精度。为此，我强烈推荐mathtool / ARPREC，可在David的主页上查阅：http://www.davidhbailey.com/。