scipy计算三重积分

圆柱体直径4cm，长20cm，两端各有一个球，每个球直径20cm，密度为5600kg/m3计算rotational inertia matrix .坐标原点为圆柱的中心，沿长度的方向为z轴。

先在草稿纸上写出直角坐标的积分公式，将其写入程序即可。计算三重积分用到的函数是tplquad，此外还有二重积分，定积分，n重积分。f为被积函数，cyl为求圆柱的积分，sphereup求上球的积分，spheredown求下球的积分，返回一个元组，第一个元素是积分值，第二个元素是精度。

from scipy import integrate
f = lambda x, y, z: x**2+y**2

cyl=integrate.tplquad(f, -0.1, 0.1, 
             lambda y: -0.02, lambda y: 0.02,
     lambda x, y: -np.sqrt(0.0004-y**2), lambda x, y: np.sqrt(0.0004-y**2))
sphereup=integrate.tplquad(f, 0.1, 0.3, 
         lambda z: -np.sqrt(0.40*z - z**2 - 0.03), 
         lambda z: np.sqrt(0.40*z - z**2 - 0.03),
     lambda z,y: -np.sqrt(0.40*z - z**2 - 0.03 - y**2), 
     lambda z,y: np.sqrt(0.40*z - z**2 - 0.03 - y**2))
spheredown=integrate.tplquad(f, -0.3, -0.1, 
         lambda z: -np.sqrt(-0.40*z - z**2 - 0.03), 
         lambda z: np.sqrt(-0.40*z - z**2 - 0.03),
     lambda z,y: -np.sqrt(-0.40*z - z**2 -0.03 - y**2), 
     lambda z,y: np.sqrt(-0.40*z - z**2 - 0.03 - y**2))

print((cyl[0]+sphereup[0]+spheredown[0])*5600)

 

当把被积函数改为1时，即为求圆柱体和球体的体积：

f=lambda x,y,z:1

 此时计算的结果为

print(cyl[0],sphereup[0],spheredown[0])

0.0002513274461536556

0.004188790206670827

0.004188790206670827

用球体、圆柱体体积公式计算结果为：

In[5]：4/3*np.pi*0.001
Out[5]: 0.0041887902047863905

In[6]: np.pi*0.02**2*0.2
Out[6]: 0.0002513274122871835