26 案例学习

26.1 计算工具

在第20章，我们将一个二阶微分方程改写成一个一阶方程，并且像这样使用一个斜率函数解决方程：

def slope_func(state, t, system):
    y, v = state
    g = system.g
    
    dydt = v
    dydt = -g
    
    return dydt, dvdt

我们使用crossings函数去寻找在仿真结果中的过零检测。

然后我们用一个像这样的事件函数：

def event_func(state, t, system):
    y, v = state
    return y

当一个事件发生时，就停止仿真。这时注意到，事件函数和斜率函数使用相同的参数。

在第21章，我们建立了一个空气阻力模型并且使用了一个Params对象，这是一个参数集合：

params = Params(height = 381 * m,
                v_init = 0 * m / s,
                g = 9.8 * m/s**2,
                mass = 2.5e-3 * kg,
                diameter = 19e-3 * m,
                rho = 1.2 * kg/m**3,
                v_term = 18 * m / s)

那么，我们可以用一个新的方式创建一个System对象，就是复制来自Params对象的所有变量并进行增加或改变变量：

    return System(params, area=area, C_d=C_d,
              init=init, t_end=t_end)

我们也可以用gradient函数去估计加速度和给定的速度：

a = gradient(results.v)

第22章引入了Vector对象，它们能够在二维或三维空间中代表一些矢量，例如位移，速度，力和加速度。

A = Vector(3,4) * m

这一章还引入了轨迹图，它可以在二维空间中显示一个对象的路径：

x = results.R.extract('x')
y = results.R.extract('y')

plot(x, y, label='trajectory')

在第23章，我们定义了一个可以计算棒球飞行距离的距离函数，它又可以作为发射角度的函数：

def range_func(angle, params):
    params = Params(params, angle=angle)
    system = make_system(params)
    results, details = run_ode_solver(system, slope_func,
                                       events=event_func)
    x_dist = get_last_value(results.R).x
    return x_dist

然后，我们用maximize寻找使距离最大的发射角度：

bounds = [0, 90] * degree
res = maximize(range_func, bounds, params)

这样，你的工具包就完成了。第24章和第25章引入了旋转物理学，但没有新的计算工具。

26.2 蹦极

假设你想要创造“蹦极扣篮”的最高世界记录，这是一个挑战蹦极者在跳跃的最低点在一杯茶中泡曲奇的特技挑战。在这个视频里有一个例子：http://modsimpy.com/dunk.

既然纪录是70米，那么我们设一个80米的蹦极目标。我们会按照下面的模型假设开始：

• 开始时，将蹦极绳索悬挂在一台起重机上，连接点在一杯茶上方的80米处。

• 在绳子完全伸展之前，它不会对蹦床者施加任何力量。这可能不是一个很好的假设；我们将重新开始。

• 在绳子完全伸展之后，它遵循胡克定律；也就是说，它对线施加一个力，这个力与绳子超出其静止长度的延伸成正比。可见 http://modsimpy.com/hooke.

• 蹦极者的质量为75公斤。

• 跳跃者受到阻力，所以他们的末速度为60米/秒。

我们的目标是设置绳子的长度L，和它的弹簧系数k，这样，跳跃着就会一直落到茶杯的位置并且不超过茶杯！

在这本书的资料里，你会找到一个笔记bungee.ipynb，它包含着本次案例学习的初始代码和练习。

26.3 重写蹦极扣篮

在前面的案例学习中，我们假设在绳被拉伸之前不会对跳跃者施加任何力。这是我们常说的因为绳子与跳跃者一起下落，所以绳子没有作用，但这种直觉是不正确的。当绳子下落时，它会将能量传递给跳跃者。

在http://modsimply.com/bungee，你会找到一篇文章[^1]解释了这一过程，并得出跳跃者的加速度a，作为位移y和速度v的函数：

g为重力加速度，L为绳子长度， 为绳子质量m和跳跃者质量M的比。

如果你不相信模型的正确性，这个视频会为你证实：http://modsimpy.com/drop.

在这本书的资料里，你会找到一个笔记bungee2.ipynb，它包含着本次案例学习的初始代码和练习。当我们改变绳子的质量时，系统的行为会发生怎样的变化？当绳子的质量等于跳跃者的质量时，跳跃最低点的净影响是什么？

[^1]Heck, Uylings, and Kdzierska, “Understanding the physics of bungee jumping”, Physics Education, Volume 45, Number 1, 2010.

26.4 蜘蛛侠

在这个案例研究中，我们将开发一个蜘蛛侠从连接在帝国大厦顶部的有弹性的织带索上摇摆的模型。最开始，蜘蛛侠在邻近的大厦顶部，就像图26.1显示的那样。

26_1.png

图26.1：蜘蛛侠案例研究的初始状态图。

原点O，在帝国大厦的底部。向量H代表织带相对于O处连接大厦的位置。向量P代表蜘蛛侠相对于O处的位置。L表示从连接点到蜘蛛侠的向量。

根据下面的箭头，从O到H，再到L，我们可以知道

那么我们可以像这样计算L:

本次的案例学习目标是：

1.实现预测蜘蛛侠轨迹的模型。

2.选择合适的时间让蜘蛛侠松开蛛丝，这样他在着陆前能走得更远。

3.为蜘蛛侠选择最好的角度跳下大楼，放开蛛网，最大限度地扩大射程。

我们将使用以下参数：

1.根据蜘蛛侠的维基，蜘蛛侠重76千克。

2.假设他的末速度为60米/秒。

3.网的长度为100米。

4.网的初始角度是45度向左下方向。

5.当绳子被拉伸时，蛛网的弹性系数为40N/m，当它被压缩时为0。

在这本书的资料里，你会找到一个笔记spiderman.ipynb，它包含着本次案例学习的初始代码。通读这个笔记并运行代码。代码中使用了minimize，它是一个可以搜索最优参数集的SciPy函数（与minimize_scalar相反，它只能沿着单个轴搜索）。

26.5 小猫

让我们模拟一只小猫在铺开厕纸。参考资料，请看这个视频：http://modsimpy.com/kitten.

小猫和纸辊的相互作用是复杂的。为了简单起见，我们假设小猫用恒定的力向下拉动纸辊的自由端。与此同时，我们将忽略滚轮与滚轴之间的摩擦。

图26.2展示了带有r，F和τ 的纸辊。作为矢量，r的方向是垂直进入纸面的，但我们现在只关注他的大小。

以下使我们将会用到的Params 对象以及参数：

26_2.png

图26.2：一卷厕纸的示意图，标注了力、杠杆臂、以及产生的扭矩。

params = Params(Rmin = 0.02 * m,
                Rmax = 0.055 * m,
                Mcore = 15e-3 * kg,
                Mroll = 215e-3 * kg,
                L = 47 * m,
                tension = 2e-4 * N,
                t_end = 180 * s)

上面出现的Rmin是最小半径，Rmax是最大半径，L是纸的长度，Mcore是纸辊中心卷心筒的质量，Mroll是纸的质量，tension是小猫所施加的力，单位为N，这里选择了一个产生合理结果的值。

在 http://modsimpy.com/moment你可以找到简单几何体的转动惯量，我将纸辊中心卷心筒视为一个“圆柱筒薄壳”，将纸辊视为一个“两端开口的厚壁圆柱管”。

薄壳的转动惯量是mr2 ,，其中m是质量，r是壳的半径。

厚壁管的转动惯量为

式中ρ为物料密度，h为筒高，r2为外径，r1为内径。

因为纸卷的外径会随着小猫的铺开而变化，我们必须计算每个时间点的转动惯量，它是一个以当前半径r为自变量的函数。下面是实现这个功能的函数:

def moment_of_inertia(r, system):
    Mcore, Rmin = system.Mcore, system.Rmin
    rho_h = system.rho_h
    
    Icore = Mcore * Rmin**2
    Iroll = pi * rho_h / 2 * (r**4 - Rmin**4)
    return Icore + Iroll

rho_h是密度和高度的乘积，ρh是单位面积上的质量。rho_h在make_system中计算:

def make_system(params):
    L, Rmax, Rmin = params.L, params.Rmax, params.Rmin
    Mroll = params.Mroll
    
    init = State(theta = 0 * radian,
                 omega = 0 * radian/s,
                 y = L)
    area = pi * (Rmax**2 - Rmin**2)
    rho_h = Mroll / area
    k = (Rmax**2 - Rmin**2) / 2 / L / radian
    
    return System(params, init=init, area=area,
                  rho_h=rho_h, k=k)

make_system同样使用公式24.4计算k.

26_3.png

图26.3：一个悠悠球的示意图，标注了重力与绳子产生的张力、张力的杠杆臂、以及产生的扭矩。

在本书的资料中，你会找到一个名为kitten.ipynb的笔记，其中包含本案例的初始代码。用它来实现这个模型，并检查结果是否可信。

26.6 模拟一个悠悠球

假设你拿着一个溜溜球，有一段绳子缠绕在它的轴上，你放下它，同时保持绳子的一端不动。当重力使溜溜球向下加速时，绳上的张力产生向上的力。由于这个力作用于一个与质心偏移的点上，它产生一个使溜溜球旋转的力矩。

图26.3是溜溜球上的力和由它们产生的扭矩的示意图。外面的阴影区域是溜溜球的球身。里面的阴影区域是卷起来的弦，它的半径随着溜溜球展开而变化

在这个模型中，我们不能独立计算出直线加速度和角加速度;我们需要解一个方程组:

这里的和表示我们在把力和力矩相加起来。

在前面的例子中，线速度和角速度是相关的，因为弦展开的方式:

在这个例子中，线加速度和角加速度有相反的正负号，当悠悠球逆时针旋转时，θ会增加而y，也就是绳子滚动的部分的长度减小了。

对两边求导得到线加速度和角加速度的类似关系:

我们可以更简洁地写为：

这种关系不是普遍适用的自然规律；它只适用于这种特定情况，即一个物体沿着另一个物体滚动而不拉伸或滑动。

由于我们设置问题的方式，y实际上有两种含义:它表示滚动的弦的长度和溜溜球的高度，溜溜球的高度随着溜溜球的下落而减小。同理，a表示滚弦长度和溜溜球高度的加速度。

我们可以通过将线性力相加来计算溜溜球的加速度:

T是正的，因为拉力向上，而mg是负的，因为重力向下。

因为重力作用在质心上，它不产生扭矩，所以唯一的扭矩是由张力产生的:

正(向上)的张力产生正(逆时针)的角加速度。

现在我们有三个方程，三个未知数，T, a，和α，和I, m, g，和r作为已知的量。用手解这些方程是很简单的，但是我们也可以使用SymPy求解:

T, a, alpha, I, m, g, r = symbols('T a alpha I m g r')
eq1 = Eq(a, -r * alpha)
eq2 = Eq(T - m*g, m * a)
eq3 = Eq(T * r, I * alpha)
soln = solve([eq1, eq2, eq3], [T, a, alpha])

结果为：

I∗是增大的转动惯量，I+mr2。要模拟这个系统，我们不需要T;我们可以把a和α直接代入到斜率函数中。

在本书的资料中，你会找到一个名为yoyo.ipynb的笔记，其中包含这些方程的推导过程好本案例的初始代码。用它来实现并测试这个模型。

本书的中文翻译由南开大学医学院智能医学工程专业2018级、2019级的师生完成，方便后续学生学习《Python仿真建模》课程。翻译人员（排名不分前后）：薛淏源、金钰、张雯、张莹睿、赵子雨、李翀、慕振墺、许靖云、李文硕、尹瀛寰、沈纪辰、迪力木拉、樊旭波、商嘉文、赵旭、连煦、杨永新、樊一诺、刘志鑫、彭子豪、马碧婷、吴晓玲、常智星、陈俊帆、高胜寒、韩志恒、刘天翔、张艺潇、刘畅。

整理校订由刘畅完成，如果您发现有翻译不当或者错误，请邮件联系changliu@nankai.edu.cn。

本书中文版不用于商业用途，供大家自由使用。

未经允许，请勿转载。