对于程序,我需要一种算法来非常快速地计算实体的体积。该形状由一个函数指定,给定一个点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)

给你

python - 区域体积的Python数值积分-LMLPHP

从那里开始,您可以使用各种程序包来提取体积。一种可能性: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

10-05 20:50
查看更多