python - 球坐标中的蒙特卡罗积分

Question

我尝试在 6 维中实现简单的蒙特卡罗积分。积分仅涵盖两个 3D 球体，在下文中球体的坐标标记为 r1 和 r2。 当使用笛卡尔坐标并忽略球体之外的任何东西时，积分效果很好。当被积函数取决于 r1 和 r2 之间的角度时，使用球坐标会失败。

看起来 r1 和 r2 很可能指向相同的方向，而我希望对齐是完全随机的。

用于生成球面坐标的变换可以在这里找到:http://de.mathworks.com/help/matlab/math/numbers-placed-randomly-within-volume-of-sphere.html

以下代码使用两个不同的被积函数以及基于笛卡尔和球坐标的蒙特卡罗积分来说明这一点:

import numpy as np
from math import *

N=1e5
R=2
INVERT = False # if this is set to True, the results are correct
INTEGRAND = None


def spherical2cartesian(r, cos_theta, cos_phi):
    sin_theta = sqrt(1-cos_theta**2)
    sin_phi = sqrt(1-cos_phi**2)
    return np.array([r*sin_theta*cos_phi, r*sin_theta*sin_phi, r*cos_theta])

# Integrands (cartesian)

def integrand_sum_sqr(r1,r2):
    r1=np.linalg.norm(r1)
    r2=np.linalg.norm(r2)
    if r1>R or r2>R:
        return 0.0;
    return r1**2 + r2**2

def integrand_diff_sqr(r1,r2):
    delta = r1-r2
    r1=np.linalg.norm(r1)
    r2=np.linalg.norm(r2)
    if r1>R or r2>R:
        return 0.0;
    return np.linalg.norm(delta)**2

# integrand for spherical integration (calls cartesian version)

def integrand_spherical(r1,r2):
    # see the following link for the transformation used: http://de.mathworks.com/help/matlab/math/numbers-placed-randomly-within-volume-of-sphere.html
    r1 = spherical2cartesian((3*r1[0])**(1.0/3.0), r1[1], cos(r1[2]))
    r2 = spherical2cartesian((3*r2[0])**(1.0/3.0), r2[1], cos(r2[2]))
    if INVERT:
        s = (-1)**np.random.choice(2,6)
        r1=s[0:3]*r1
        r2=s[3:6]*r2
    return INTEGRAND(r1,r2)

# monte carlo integration routines

def monte_carlo(name,func,samples,V):
    results=np.array([func(x[0:3],x[3:6]) for x in samples])
    avg = results.mean()*V
    std = results.std(ddof=1)*V/sqrt(len(samples))
    print name,": ",avg," +- ",std

def monte_carlo_cartesian():
    V=(2*R)**6
    samples = np.random.rand(N,6)
    samples = R*(2*samples-1)
    monte_carlo('cartesian',INTEGRAND,samples,V)

def monte_carlo_spherical():
    V=(4.0/3.0*R**3*pi)**2
    samples = np.random.rand(6,N)
    samples = np.array([R**3*samples[0]/3.0, 2*samples[1]-1, 2*pi*samples[2], R**3*samples[3]/3.0, 2*samples[4]-1, 2*pi*samples[5]])
    samples = samples.T
    monte_carlo('spherical',integrand_spherical,samples,V)

# run all functions with all different monte carlo versions
print "Integrating sum_sqr, expected:",32.0/15.0*R**8*pi**2
INTEGRAND=integrand_sum_sqr
monte_carlo_cartesian()
monte_carlo_spherical()

print "Integrating diff_sqr, expected:",32.0/15.0*R**8*pi**2
INTEGRAND=integrand_diff_sqr
monte_carlo_cartesian()
monte_carlo_spherical()

典型输出:

Integrating sum_sqr, expected: 5390.11995025
cartesian :  5406.6226422  +-  29.5405030567
spherical :  5392.72811794  +-  5.23871574928
Integrating diff_sqr, expected: 5390.11995025
cartesian :  5360.20055643  +-  34.8851044924
spherical :  4141.88319573  +-  9.69351527901

最后的积分显然是错误的。 为什么 r1 和 r2 相关？我该如何解决这个问题？

Answer 1

上面代码的问题在下面一行

sin_phi = sqrt(1-cos_phi**2)

这只会产生正结果，而 sin(phi) 会产生 phi > pi 的负结果