四十米的大刀,你有几秒逃生

最近看了李永乐老师的课程,其中有一期是“抽出四十米的大刀,你有几秒钟的逃生时间”。起源于一个网上小视频,爸爸辅导女儿的作业,在填写单位的时候,女儿给铅笔填单位时,填了米。爸爸说,铅笔的单位怎么可能说米呢,信不信我抽出四十米的大刀,女儿气急的说,你的四十米的大刀刚好削我15米的铅笔,哈哈哈。详细可以看看那期视频。

问题

我们一起来看看这个脑洞,假设爸爸抽取L米长度的大刀,大刀只在重力g的作用下下落,一端握在手里,另一侧则是转动下落。这是一个典型的刚体转动问题,详细过程可以参考视频。这里给出最终的积分公式:

t = \sqrt{\frac{L}{3g}}\int_0 ^\theta \frac{\rm d \theta}{\sqrt{\sin \theta_0 – \sin \theta}}

抽出40米的大刀,从30度开始下落,逃生的时间是1.78秒。我们就来算一下这个定积分,当然是编程求解了。

直接调用函数求解

使用数值积分的函数scipyintegrate函数。完整程序如下:

from scipy import integrate, pi, sin, sqrt, linspace

# 初始参数
theta = pi / 6  # 初始角度
L = 40  # 大刀长度
g = 9.8  # 重力加速度

f = lambda x: sqrt(L / (3 * g)) / (sqrt(sin(theta) - sin(x)))
print(integrate.quad(f, 0, theta))

# 输出结果为:
# (1.777013730182824, 2.5051516416851882e-11)
# 第一个是结果,第二个是误差

编写积分函数

我们知道一维定积分就是曲线下的面积,于是我们可以根据定义来计算这个定积分。

画出函数图像,用矩形来代替曲线下的面积,当划分越细的时候,矩形的面积和就越接近曲线下的面积。我们只需要计算这些矩形的和就可以了。

from scipy import pi, sin, sqrt, linspace

# 初始参数
theta = pi / 6  # 初始角度
L = 40  # 大刀长度
g = 9.8  # 重力加速度
N = 10000  # 划分的


def compute(theta, N):
    x = linspace(0, theta, N)[:-1]  # 分割x,舍弃最后一个x坐标
    f = lambda x: sqrt(L / (3 * g)) / (sqrt(sin(theta) - sin(x)))
    y = f(x)
    delta_x = theta / N  # 等距离划分,每个划分的长度
    return delta_x * sum(y)  # 矩形的求和


print(compute(theta, N=N))
# 1.7652852253068916

当划分次数N=10000时,结果是1.765,误差还是比较大的,可以进一步增大N的值,使得划分更加细,误差会不断减小。分别尝试了一下几种分割。

x 1e3 1e4 1e5 1e6 1e7
Y 1.735 1.765 1.774 1.777 1.778

最终我们得到非常接近1.78s的结果。如果大砍刀的初始角度是60度,那么逃生时间是多少呢?我们只需要改一下theta=pi/3,运行程序就能得到3.036秒(N=1e5)。

总结

可能看过课程的小朋友会好奇最后的1.78秒到底是怎么得到,今天我们就用编程求解了这个积分表达式,并且分别使用了两种方法,第一种是别人写好的函数,计算的速度和精度都很好;第二种方法是根据积分的定义,比较粗糙的求解。从这里可以看出编程是一个不断优化计算方法(算法)的过程,有兴趣的小朋友可以自己尝试。

关注机器学习和算法的码农,喜欢编程和读书
文章已创建 74

发表评论

电子邮件地址不会被公开。

相关文章

开始在上面输入您的搜索词,然后按回车进行搜索。按ESC取消。

返回顶部