对于程序,我需要一种算法来非常快速地计算实体的体积。该形状由一个函数指定,给定一个点P(x,y,z),如果P是实体的一个点,则返回1,如果P不是实体的一个点,则返回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
但它失败,给我以下错误:
因此,自然地,我寻找“特殊用途的集成商”,但找不到能满足我需求的任何产品。
然后,我尝试使用蒙特卡洛方法编写自己的积分,并以相同的形状对其进行测试:
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蒙特卡洛方法,但遇到了一些困难,但我仍不确定它会快多少。
因此,我要问是否有任何库可以满足我的要求,或者是否有更好的算法可以工作几倍(对于相同的精度)。任何建议都是值得欢迎的,因为我已经为此工作了一段时间了。
编辑:
如果我无法使用Python进行此工作,那么我愿意切换到任何其他既可编译又具有相对简单的GUI功能的语言。欢迎任何建议。
最佳答案
正如其他人已经指出的那样,很难找到 bool 函数给定的域的数量。您可以使用pygalmesh(属于我的一个小项目),它位于CGAL顶部,并为您提供四面体网格。这个
import numpy
import pygalmesh
import meshplex
class Custom(pygalmesh.DomainBase):
def __init__(self):
super(Custom, self).__init__()
return
def eval(self, x):
return (x[0]**2 + x[1]**2 + x[2]**2) - 1.0
def get_bounding_sphere_squared_radius(self):
return 2.0
mesh = pygalmesh.generate_mesh(Custom(), cell_size=1.0e-1)
给你
从那里开始,您可以使用各种程序包来提取体积。一种可能性:meshplex(我动物园里的另一种):
import meshplex
mp = meshplex.MeshTetra(mesh.points, mesh.cells["tetra"])
print(numpy.sum(mp.cell_volumes))
给
4.161777
足够接近真实值
4/3 pi = 4.18879020478...
。如果要提高精度,请减少上述网格生成中的cell_size
。