蒙特卡罗在球坐标系中的积分

我试图在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是相关的？我该如何解决这个问题？

来源

2016-06-24 val

与上述代码的问题是在下面的行

sin_phi = sqrt(1-cos_phi**2)

这仅产生正结果，而sin(phi)产生阴性结果phi > pi

来源

2016-06-27 13:20:37 val

蒙特卡罗在球坐标系中的积分

回答

相关问题