1

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

我的想法是这样的：我们放置半径为 R 的球体，使其中心位于原点，我们放置半径为 R' 的圆柱体，使其轴在 z 中从 (b,0,0 ）。我们在球体上积分，使用一个阶跃函数，如果我们在圆柱体内部，则返回 1，如果不是，则返回 0，因此在受球体和圆柱体内部约束的集合上积分 1，即相交。

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

警告：已达到最大细分数 (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()
4

4 回答 4

3

假设您能够转换和缩放数据，使得球体的原点位于其中[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
于 2011-04-26T15:36:55.933 回答
2

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

我建议一个更好的主意是从分析上减少问题。

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

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

在该轴变量上积分。沿轴的任何固定值的相交面积只是两个圆的相交面积，而这反过来又可以使用三角函数和一点点努力简单地计算出来。

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

于 2011-04-26T14:39:48.357 回答
1

正如 eat 所建议的，我使用简单的 MC 集成解决了它，但我的实现速度很慢。我的要求增加了。因此，正如木片所建议的那样，我在数学上重新表述了这个问题。

基本上，我将 x 的极限表述为 z 和 y 的函数，并将 y 表述为 z 的函数。然后，i 本质上是在交集上积分 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]
于 2011-04-27T16:16:10.217 回答
0

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

我会为域生成一个四面体网格，然后将单元体积相加。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为卷。减小单元或边缘尺寸以获得更高的精度。

于 2017-03-24T13:02:32.063 回答