区域体积的Python数值积分 [英] Python Numerical Integration for Volume of Region

问题描述

对于一个程序，我需要一个算法来非常快速地计算固体的体积.此形状由一个函数指定，给定一个点P(x，y，z)，如果P是实体的一个点，则返回1，如果P不是实体的一个点，则返回0.

For a program, I need an algorithm to very quickly compute the volume of a solid. This shape is specified by a function that, given a point P(x,y,z), returns 1 if P is a point of the solid and 0 if P is not a point of the solid.

我尝试通过以下测试使用numpy:

I have tried using numpy using the following test:

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

但是它没有给我以下错误:

but it fails giving me the following errors:

警告(来自警告模块): 文件"C:\ Python27 \ lib \ site-packages \ scipy \ integrate \ quadpack.py"，第321行 warnings.warn(msg，IntegrationWarning) IntegrationWarning:已达到最大细分数量(50). 如果增加极限没有改善，建议进行分析 为了确定难点.如果一个位置 可以确定局部难度(奇异，不连续) 可能会受益于拆分间隔并调用积分器 在子范围上.也许应该使用专用集成器.

Warning (from warnings module): File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, IntegrationWarning) IntegrationWarning: The maximum number of subdivisions (50) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used.

警告(来自警告模块): 文件"C:\ Python27 \ lib \ site-packages \ scipy \ integrate \ quadpack.py"，第321行 warnings.warn(msg，IntegrationWarning) IntegrationWarning:算法不收敛.检测到舍入错误 在外推表中.假设要求的公差 无法实现，并且返回的结果(如果full_output = 1)为 可以得到的最好的.

Warning (from warnings module): File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, IntegrationWarning) IntegrationWarning: The algorithm does not converge. Roundoff error is detected in the extrapolation table. It is assumed that the requested tolerance cannot be achieved, and that the returned result (if full_output = 1) is the best which can be obtained.

警告(来自警告模块): 文件"C:\ Python27 \ lib \ site-packages \ scipy \ integrate \ quadpack.py"，第321行 warnings.warn(msg，IntegrationWarning) IntegrationWarning:检测到舍入错误的发生，这可以防止 达到要求的公差.错误可能是 被低估了.

Warning (from warnings module): File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, IntegrationWarning) IntegrationWarning: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be underestimated.

警告(来自警告模块): 文件"C:\ Python27 \ lib \ site-packages \ scipy \ integrate \ quadpack.py"，第321行 warnings.warn(msg，IntegrationWarning) IntegrationWarning:积分可能是发散的或缓慢收敛的.

Warning (from warnings module): File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, IntegrationWarning) IntegrationWarning: The integral is probably divergent, or slowly convergent.

因此，自然地，我寻找特殊用途的集成商"，但找不到能满足我需要的任何东西.

So, naturally, I looked for "special-purpose integrators", but could not find any that would do what I needed.

然后，我尝试使用Monte Carlo方法编写自己的积分并以相同的形状对其进行测试:

Then, I tried writing my own integration using the Monte Carlo method and tested it with the same shape:

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位精度更高.

but this works too slowly. For this program I will be using this numerical integration about 100 times or so, and I will also be doing it on larger shapes, which will take minutes if not an hour or two at the rate it goes now, not to mention that I want a better precision than 2 decimal places.

我曾尝试实现MISER蒙特卡洛方法，但遇到了一些困难，但我仍不确定它会快多少.

I have tried implementing a MISER Monte Carlo method, but was having some difficulties and I'm still unsure how much faster it would be.

因此，我要问是否有任何库可以满足我的要求，或者是否有更好的算法可以工作几倍(以相同的精度).任何建议都值得欢迎，因为我已经为此工作了一段时间了.

So, I am asking if there are any libraries that can do what I am asking, or if there are any better algorithms which work several times faster (for the same accuracy). Any suggestions are welcome, as I've been working on this for quite a while now.

如果我无法使用Python进行此工作，则可以切换到任何其他可编译且具有相对简单的GUI功能的语言.欢迎任何建议.

If I cannot get this working in Python, I am open to switching to any other language that is both compilable and has relatively easy GUI functionality. Any suggestions are welcome.

推荐答案

正如其他人已经指出的那样，很难找到布尔函数给定的域的数量.您可以使用 pygalmesh (我的一个小项目)位于 CGAL 的顶部，并为您提供四面体网格.这个

As the others have already noted, finding the volume of domain that's given by a Boolean function is hard. You could used pygalmesh (a small project of mine) which sits on top of CGAL and gives you back a tetrahedral mesh. This

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 (我动物园里的另一个):

From there on out, you can use a variety of packages for extracting the volume. One possibility: meshplex (another one out of my zoo):

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.

which is close enough to the true value 4/3 pi = 4.18879020478.... If want more precision, decrease the cell_size in the mesh generation above.

这篇关于区域体积的Python数值积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！