圆柱体直径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