动机：我有一个多维积分，为了完整起见，我在下面复制了它。它来自于存在显着各向异性时的第二维里系数的计算：
这里 W 是所有变量的函数。这是一个已知函数，我可以为它定义一个 python 函数。
编程问题：如何scipy整合这个表达式？我正在考虑将两个三重四边形 ( scipy.integrate.tplquad) 链接在一起，但我担心性能和准确性。中是否有更高维的积分器scipy，可以处理任意数量的嵌套积分？如果没有，最好的方法是什么？
问问题
18183 次
4 回答
32
对于像这样的高维积分，蒙特卡洛方法通常是一种有用的技术 - 它们将答案收敛为函数评估次数的平方根倒数，这对于高维更好，然后你通常会得到偶数相当复杂的自适应方法（除非你知道关于你的被积函数的一些非常具体的东西——可以利用的对称性等）
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 %
于 2012-12-28T20:46:07.960 回答
9
Jonathan Dursi 给出了一个非常好的答案。我将补充他的答案。
现在scipy.integrate有一个名为的函数nquad，可以毫不费力地执行多维积分。有关更多信息，请参阅此链接。nquad下面我们使用Jonathan 的示例计算积分：
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
于 2018-10-23T21:42:10.193 回答
6
我将就如何准确地进行这种积分做一些一般性的评论，但这个建议并不特定于 scipy（评论太长，即使它不是答案）。
我不知道您的用例，即您是否对具有几位数准确度的“好”答案感到满意，这可以使用 Jonathan Dursi 的回答中概述的蒙特卡洛直接获得，或者您是否真的想推动数字尽可能准确。
我自己对维里系数进行了分析、蒙特卡罗和正交计算。如果你想准确地做积分，那么你应该做一些事情：
尝试尽可能多地进行积分；很可能在您的某些坐标中集成非常简单。
考虑转换积分变量，使被积函数尽可能平滑。（这对蒙特卡洛和正交都有帮助）。
对于 Monte Carlo，使用重要性采样以获得最佳收敛。
对于求积，使用 7 个积分，使用 tanh-sinh 求积可能获得真正快速的收敛。如果您可以将其降低到 5 个积分，那么您应该能够为您的积分获得 10 位数的精度。为此，我强烈推荐 mathtool / ARPREC，可从 David Bailey 的主页获得：http ://www.davidhbailey.com/
于 2013-01-07T20:56:58.050 回答
1
首先要说我数学不是很好，所以请善待。无论如何，这是我的尝试：
请注意，在您的问题中有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*}
你可以为你的问题多玩一点；）
于 2012-12-28T17:35:03.510 回答