動機:多次元積分がありますが、完全を期すために以下に再掲します。これは、有意な異方性がある場合の 2 番目のビリアル係数の計算から得られます。
ここで、W はすべての変数の関数です。これは既知の関数であり、python 関数を定義できる関数です。
プログラミングの質問:scipyこの式を統合するにはどうすればよいですか? トリプルクワッド（scipy.integrate.tplquad）を2つつなげようと思っていたのですが、性能や精度が心配です。scipyに、任意の数のネストされた積分を処理できる高次元の積分器はありますか? そうでない場合、これを行う最善の方法は何ですか?
このような高次元の積分では、モンテカルロ法はしばしば有用な手法です - 関数評価の数の逆平方根として答えに収束します。かなり洗練された適応法 (被積分関数について非常に具体的なことを知っていない限り - 悪用できる対称性など)
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 %
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)
結果は chained よりも正確ですtplquad:
>>> print(diff)
0.0
この種の積分を正確に行う方法について、いくつかの一般的なコメントを作成しますが、このアドバイスは scipy に固有のものではありません (回答ではありませんが、コメントには長すぎます)。
あなたのユースケース、つまり、ジョナサン・ドゥルシの回答で概説されているように、モンテカルロを使用して直接取得できる数桁の精度の「良い」回答に満足しているかどうか、または本当に数値をプッシュしたいかどうかはわかりません可能な限り正確に。
私は自分でビリアル係数の分析、モンテカルロ、および直交計算を実行しました。積分を正確に行いたい場合は、次のことを行う必要があります。
できるだけ多くの積分を正確に実行しようとします。いくつかの座標での統合は非常に簡単かもしれません。
被積分関数ができるだけ滑らかになるように、積分変数を変換することを検討してください。(これは、モンテカルロと求積法の両方に役立ちます)。
モンテカルロの場合、最適な収束のために重要度サンプリングを使用します。
求積の場合、7 つの積分を使用すると、tanh-sinh 求積を使用して非常に高速な収束が得られる可能性があります。積分を 5 まで下げることができれば、積分の精度を数十桁にすることができるはずです。この目的には、David Bailey のホームページから入手できる mathtool / ARPREC を強くお勧めします: http://www.davidhbailey.com/
最初に言っておきますが、私は数学があまり得意ではないので、親切にしてください。とにかく、ここに私の試みがあります:
あなたの質問には6つの変数がありますが、 7つの積分があることに注意してください!? 使用
中： PythonSympy
>>> 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*}
あなたはあなたの質問のためにもう少し遊ぶかもしれません;）