3d积分，python，积分集约束

Question

我想计算球体和无限圆柱体在一定距离 b 处相交的体积，我想我会使用一个快速而肮脏的 python 脚本来完成它。我的要求是 <1 秒的计算，> 3 位有效数字。

我的想法是这样的： 我们放置半径为 R 的球体，使其中心位于原点，并放置半径为 R' 的圆柱体，使其轴从 (b,0,0) 开始沿 z 方向延伸。我们在球体上进行积分，使用阶跃函数，如果我们在圆柱体内部，则返回 1，如果不在圆柱体内部，则返回 0，从而在球体和圆柱体内部约束的集合（即相交）上积分 1。

我使用 scipy.intigrate.tplquad 尝试过此操作。但没有成功。我认为这是因为阶跃函数的不连续性，因为我收到以下警告。当然，我可能只是做错了。假设我没有犯一些愚蠢的错误，我可以尝试制定相交的范围，从而消除对阶跃函数的需要，但我想我可能会尝试先获得一些反馈。任何人都可以发现任何错误，或者指出一些简单的解决方案。

警告：最大数量 已实现细分 (50)。
如果增加限制不会产生 建议分析改进
被积函数以确定 的困难。如果位置为 局部困难可能是 确定（奇异性， 不连续性）人们可能会 分割间隔的收益 并致电积分商 子范围。也许有特殊用途 应使用积分器。

代码：

from scipy.integrate import tplquad
from math import sqrt


def integrand(z, y, x):
    if Rprim >= (x - b)**2 + y**2:
        return 1.
    else:
        return 0.

def integral():
    return tplquad(integrand, -R, R, 
                   lambda x: -sqrt(R**2 - x**2),          # lower y
                   lambda x: sqrt(R**2 - x**2),           # upper y
                   lambda x,y: -sqrt(R**2 - x**2 - y**2), # lower z
                   lambda x,y: sqrt(R**2 - x**2 - y**2), # upper z
                   epsabs=1.e-01, epsrel=1.e-01
                   )

R=1
Rprim=1
b=0.5
print integral()

Answer 1

假设您能够平移和缩放数据，使得球体的原点位于 [0, 0, 0] 且其半径为 1，那么简单的随机近似可以足够快地给你一个合理的答案。因此，类似的事情可能是一个很好的起点：

import numpy as np

def in_sphere(p, r= 1.):
    return np.sqrt((p** 2).sum(0))<= r

def in_cylinder(p, c, r= 1.):
    m= np.mean(c, 1)[:, None]
    pm= p- m
    d= np.diff(c- m)
    d= d/ np.sqrt(d** 2).sum()
    pp= np.dot(np.dot(d, d.T), pm)
    return np.sqrt(((pp- pm)** 2).sum(0))<= r

def in_sac(p, c, r_c):
    return np.logical_and(in_sphere(p), in_cylinder(p, c, r_c))

if __name__ == '__main__':
    n, c= 1e6, [[0, 1], [0, 1], [0, 1]]
    p= 2* np.random.rand(3, n)- 2
    print (in_sac(p, c, 1).sum()/ n)* 2** 3

Answer 2

对两个域上恒定的不连续函数执行三重自适应数值积分是一个非常糟糕的主意，特别是如果您希望看到速度或准确性。

我建议更好的想法是通过分析来减少问题。

通过变换将圆柱体与轴对齐。这会将球体平移到不在原点的某个点。

现在，找到球体与圆柱体沿该轴相交的极限。

对该轴变量进行积分。沿轴的任何固定值的相交面积就是两个圆的相交面积，而这又可以使用三角学和一点点努力来简单计算。

最后，您将获得准确的结果，几乎不需要计算时间。

Answer 3

我按照 eat 的建议，使用简单的 MC 集成解决了这个问题，但我的实现速度很慢。我的要求增加了。因此，我按照木片的建议，用数学方法重新表述了这个问题。

基本上，我将 x 的极限表示为 z 和 y 的函数，并将 y 表示为 z 的函数。然后，本质上，我使用极限在交点上积分 f(z,y,z)=1。我这样做是因为速度增加，允许我绘制体积与 b 的关系，并且因为它允许我通过相对较小的修改来集成更复杂的函数。

我包含我的代码，以防有人感兴趣。

from scipy.integrate import quad
from math import sqrt
from math import pi

def x_max(y,r):
    return sqrt(r**2-y**2)

def x_min(y,r):
    return max(-sqrt(r**2 - y**2), -sqrt(R**2 - y**2) + b) 

def y_max(r):
    if (R<b and b-R<r) or (R>b and b-R>r):
        return sqrt( R**2 - (R**2-r**2+b**2)**2/(4.*b**2) )
    elif r+R<b:
        return 0.
    else: #r+b<R
        return r

def z_max():
    if R>b:
        return R
    else:
        return sqrt(2.*b*R - b**2) 

def delta_x(y, r):
    return  x_max(y,r) - x_min(y,r)

def int_xy(z):
    r = sqrt(R**2 - z**2)
    return quad(delta_x, 0., y_max(r), args=(r))

def int_xyz():
    return quad(lambda z: int_xy(z)[0], 0., z_max())

R=1.
Rprim=1.
b=0.5
print 4*int_xyz()[0]

Answer 4

首先：您可以手动计算交叉点的体积。如果您不想（或不能）这样做，这里有一个替代方案：

我会为域生成一个四面体网格，然后将单元体积相加。使用 pygalmesh 和 meshplex （均由我自己创作）：

import pygalmesh
import meshplex
import numpy

ball = pygalmesh.Ball([0, 0, 0], 1.0)
cyl = pygalmesh.Cylinder(-1, 1, 0.7, 0.1)
u = pygalmesh.Intersection([ball, cyl])

mesh = pygalmesh.generate_mesh(u, cell_size=0.05, edge_size=0.1)

points = mesh.points
cells = mesh.cells["tetra"]

# kick out unused vertices
uvertices, uidx = numpy.unique(cells, return_inverse=True)
cells = uidx.reshape(cells.shape)
points = points[uvertices]

mp = meshplex.MeshTetra(points, cells)
print(sum(mp.cell_volumes))

这为您提供

并打印 2.6567890958740463 作为体积。减小单元或边缘尺寸以获得更高的精度。