三重积分，Python，约束积分集

Question

我想计算一个球体和一个无限长圆柱体在某个距离b处的交集体积，于是决定用一个简单的Python脚本来实现。我的要求是计算时间少于1秒，并且结果要有超过3位有效数字。

我的想法是这样的：我们把半径为R的球体放在原点，然后把半径为R'的圆柱体放在z轴上，从(b,0,0)开始延伸。接着，我们在球体内部进行积分，使用一个阶跃函数：如果我们在圆柱体内部，返回1；如果不在，则返回0。这样，我们就可以在同时在球体和圆柱体内部的区域进行积分，也就是计算它们的交集。

我尝试使用scipy.integrate.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

在处理一个在两个区域内是常数的间断函数时，进行三重自适应数值积分是个非常糟糕的主意，尤其是如果你想要速度或准确性的话。

我建议一个更好的方法是从分析上简化这个问题。

可以通过变换将圆柱体与某个轴对齐。这样做会把球体移动到一个不在原点的位置。

接下来，找出球体与圆柱体在这个轴上相交的范围。

然后在这个轴的变量上进行积分。在轴上任何固定值的相交区域，实际上就是两个圆的相交面积，这个面积可以通过三角函数和一些简单的计算得出。

最后，你将得到一个精确的结果，几乎不需要计算时间。

Answer 2

假设你能把你的数据处理得很好，让球体的中心点在 [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 3

我按照eat的建议，用简单的蒙特卡罗积分法解决了这个问题，但我的实现速度太慢了。我的需求增加了，所以我按照woodchips的建议，重新从数学上整理了这个问题。

简单来说，我把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]