Python数值积分求区域体积

Question

我需要一个算法来快速计算一个立体的体积，这个形状由一个函数指定，给定一个点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功能的语言。欢迎提出任何建议。

- Lambda

---

值得注意的是，形状的“大小”（净体积）对于数值积分或蒙特卡罗方法并不重要。然而，形状边界的复杂性将改变某些算法的收敛性。 - Hooked

@Hooked，形状的大小决定了积分区域的大小，该区域越大，需要更多的点来获得相同的精度。 - Lambda

所有形状都可以缩放以适应某个边界框。从那里开始，一个简单的MC采样算法将以相同的时间常数运行，无论净体积如何。我认为我们在定义精度方面有所不同 - 通常我认为误差是相对于尺寸而言的，而不是绝对误差。例如，如果我要积分一个苹果和地球的体积，对于其中一个，我不会在意我偏离了几米... - Hooked

@Hooked 那很有道理；我指的是绝对误差，但正如你所提到的，相对误差可能更合理。 - Lambda

Answer 1

正如其他人已经指出的那样，找到由布尔函数给定的域的体积很困难。您可以使用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。

- Nico Schlömer

Answer 2

你的函数不是连续函数，我认为这样做积分会很困难。

那么怎么办呢：

import numpy as np

def sphere(x,y,z):
    return x**2 + y**2 + z**2 <= 1

x, y, z = np.random.uniform(-2, 2, (3, 2000000))

sphere(x, y, z).mean() * (4**3), 4/3.0*np.pi

输出：

(4.1930560000000003, 4.1887902047863905)

或者VTK：

from tvtk.api import tvtk

n = 151
r = 2.0
x0, x1 = -r, r
y0, y1 = -r, r
z0, z1 = -r, r
X,Y,Z = np.mgrid[x0:x1:n*1j, y0:y1:n*1j, z0:z1:n*1j]
s = sphere(X, Y, Z)

img = tvtk.ImageData(spacing=((x1-x0)/(n-1), (y1-y0)/(n-1), (z1-z0)/(n-1)), 
                     origin=(x0, y0, z0), dimensions=(n, n, n))

img.point_data.scalars = s.astype(float).ravel()

blur = tvtk.ImageGaussianSmooth(input=img) 
blur.set_standard_deviation(1)
contours = tvtk.ContourFilter(input = blur.output) 

contours.set_value(0, 0.5)
mp = tvtk.MassProperties(input = contours.output)
mp.volume, mp.surface_area

输出：

4.186006622559839, 12.621690438955586

- HYRY

---

如果函数是矢量化的，那么您的numpy示例是高效的方法。对于某些情况，您可以在积分中使用numexpr使其更快，这将并行评估表达式。 - Davidmh

Answer 3

如果没有给出一点边界提示，这似乎非常困难：

import numpy as np
from scipy.integrate import *

def integrand(z,y,x):
    return 1. if x**2 + y**2 + z**2 <= 1. else 0.

g=lambda x: -2
h=lambda x: 2
q=lambda x,y: -np.sqrt(max(0, 1-x**2-y**2))
r=lambda x,y: np.sqrt(max(0, 1-x**2-y**2))
I=tplquad(integrand,-2.,2.,g,h,q,r)

print I

- mtrbean

---

我刚试着运行了你的代码，它可以工作！但是我想问一下：为什么你的代码能够正常工作，而我的代码却输出了几个“警告”？你只是基本上改变了一个维度中的积分限制，那么它不还是在圆柱体内进行积分，并在边界点处出现问题吗？ - Lambda