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

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

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

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

Warning: 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.

代码:

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() 

请您参考如下方法:

假设您能够以球体原点位于 [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