我需要一个算法来快速计算一个立体的体积,这个形状由一个函数指定,给定一个点P(x,y,z),如果P是立体中的一个点则返回1,否则返回0。
我尝试使用numpy进行以下测试:
import numpy
from scipy.integrate import *
def integrand(x,y,z):
if x**2. + y**2. + z**2. <=1.:
return 1.
else:
return 0.
g=lambda x: -2.
f=lambda x: 2.
q=lambda x,y: -2.
r=lambda x,y: 2.
I=tplquad(integrand,-2.,2.,g,f,q,r)
print I
但是,它失败了,并显示以下错误:
警告(来自警告模块): 文件“C:\ Python27 \ lib \ site-packages \ scipy \ integrate \ quadpack.py”,第321行 warnings.warn(msg,IntegrationWarning) 积分警告:已达到最大子区间数(50)。 如果增加限制没有改善建议分析被积函数以确定困难所在位置。 如果可以确定本地困难的位置(奇异性,不连续),则可能从将区间分割并在子范围上调用积分器中获得。 可能应使用专用积分器。
警告(来自警告模块): 文件“C:\ Python27 \ lib \ site-packages \ scipy \ integrate \ quadpack.py”,第321行 warnings.warn(msg,IntegrationWarning) 积分警告:算法不收敛。 在外推表中检测到舍入误差 假定无法实现所请求的公差,并且返回的结果(如果full_output = 1)是 可以获得的最佳结果。
警告(来自警告模块): 文件“C:\ Python27 \ lib \ site-packages \ scipy \ integrate \ quadpack.py”,第321行 warnings.warn(msg,IntegrationWarning) 积分警告:检测到舍入误差,阻止 达到所需的公差。 错误可能是 被低估。
警告(来自警告模块): 文件“C:\ Python27 \ lib \ site-packages \ scipy \ integrate \ quadpack.py”,第321行 warnings.warn(msg,IntegrationWarning) 积分警告:积分可能发散或收敛缓慢。
因此,自然地,我寻找了“专用积分器”,但是找不到任何可以满足我的需求的积分器。
然后,我尝试使用蒙特卡罗方法编写自己的积分,并使用相同的形状进行了测试。
import random
# Monte Carlo Method
def get_volume(f,(x0,x1),(y0,y1),(z0,z1),prec=0.001,init_sample=5000):
xr=(x0,x1)
yr=(y0,y1)
zr=(z0,z1)
vdomain=(x1-x0)*(y1-y0)*(z1-z0)
def rand((p0,p1)):
return p0+random.random()*(p1-p0)
vol=0.
points=0.
s=0. # sum part of variance of f
err=0.
percent=0
while err>prec or points<init_sample:
p=(rand(xr),rand(yr),rand(zr))
rpoint=f(p)
vol+=rpoint
points+=1
s+=(rpoint-vol/points)**2
if points>1:
err=vdomain*(((1./(points-1.))*s)**0.5)/(points**0.5)
if err>0:
if int(100.*prec/err)>=percent+1:
percent=int(100.*prec/err)
print percent,'% complete\n error:',err
print int(points),'points used.'
return vdomain*vol/points
f=lambda (x,y,z): ((x**2)+(y**2)<=4.) and ((z**2)<=9.) and ((x**2)+(y**2)>=0.25)
print get_volume(f,(-2.,2.),(-2.,2.),(-2.,2.))
但是这个程序运行得太慢了。对于这个程序,我将使用大约100次左右的数值积分,并且我还要在更大的形状上进行计算,如果按照现在的速度进行计算,需要几分钟甚至一两个小时,更不用说我想要比2位小数更好的精度了。
我尝试过实现MISER Monte Carlo方法,但是遇到了一些困难,而且我仍然不确定它会快多少。
因此,我想问是否有任何库可以完成我的要求,或者是否有任何更好的算法可以比同等精度下运行速度快几倍。欢迎提出任何建议,因为我已经花了相当长的时间在这个问题上了。
编辑:
如果我无法在Python中解决这个问题,我愿意转换到任何其他既可编译又具有相对易用的GUI功能的语言。欢迎提出任何建议。