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 产生消极的结果