20 抛物体

到目前为止，我们使用的微分方程是一阶的，这意味着它们只涉及一阶导数。在本章中，我们将注意力转向二阶常微分方程，它可能涉及一阶和二阶导数。

我们将重新访问第1章中的“自由落体的硬币”示例，并使用run_ode_solver分别在考虑和不考虑空气阻力的情况下查找硬币下降时的位置和速度。

20.1 牛顿第二运动定律

一阶常微分方程组可以写作

其中 G 是 x 和 y 的某个函数（参见http://modsimpy.com/ode). 二阶常微分方程组可写作

其中 H 是 x , y 和 dy/dx 的函数。

在这一章里，我们将用到最著名并且很实用的牛顿第二运动定律，其公式是一个二阶常微分方程：

其中F是力或一组力的总和，m 是运动物体的质量，a 是它的加速度。

牛顿定律可能看起来不像微分方程，除非我们意识到加速度 a 是位置 y 对时间 t 的二阶导数

所以牛顿第二定律也可以写成

这绝对是二阶常微分方程。一般来说，F 可以是时间、位置和速度的函数。

当然，这个“定律”实际上是一个模型，因为它是对现实世界的简化。虽然这通常是正确的：

• 它只适用于 m 为常数的情况。如果质量取决于时间、位置或速度，我们必须使用牛顿定律的更一般形式（见http://modsimpy.com/varmass).
• 对于非常微观的事物，它不是一个好的模型，而另一个模型，量子力学可以更好地描述它们。
• 对于高速运动的物体，它并不是一个很好的模型，可以用另一个模型，相对论力学来描述。

然而，对于质量恒定、以中等速度运动的中型物体，牛顿模型是非常有用的。如果我们能量化作用在这样一个物体上的力，我们就可以预测它将如何移动。

20.2 抛硬币

作为第一个例子，让我们回想在1.1节中讨论过的帝国大厦掉下来的硬币。我们将实施两种模式的系统：一种是不考虑空气阻力，另一种考虑空气阻力。

考虑到帝国大厦高381米，假设硬币是从静止状态掉落的，初始条件是：

init = State(y=381 * m,
            v=0 * m/s)

其中 y 是地平面以上的高度，v 是速度。m 和 s 是 Pint 里提供的单位。

m = UNITS.meter
s = UNITS.second

唯一的系统参数是重力加速度：

g = 9.8 * m/s**2

此外，我们还将指定模拟的持续时间和步长：

t_end = 10 * s
dt = 0.1 * s

有了这些参数，时间步数是100，这对于许多问题来说是很好的。一旦我们有了解决方案，我们将增加步骤的数量，看看它对结果有什么影响。
我们需要一个系统来存储参数：

system = System(init=init, g=g, t_end=t_end, dt=dt)

现在我们需要一个斜率函数，这就是问题的症结所在。如我们所见，run_ode_solver可以解一阶常微分方程组，但牛顿定律是二阶常微分方程组。但是，如果我们认识到这一点

1. 速度 v 是位置的导数，dy/dt
2. 加速度 a 是速度的导数，dv/dt

我们可以将牛顿定律改写为一阶常微分方程组：

并且我们可以把这些等式转换成一个斜率函数：

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

第一个参数state包含硬币的位置和速度。最后一个参数system包含系统参数g，它是由于重力引起的加速度的大小。

第二个参数t是时间。该斜率函数中不使用它，因为模型中的任何因素都不是时间相关的（见第18.1节）。我之所以包含它，是因为这个函数将由run_ode_solver调用，它总是提供相同的参数，不管它们是否需要。

函数的其余部分是微分方程的直接转换，替换a=-g，这表示重力引起的加速度朝着y减小的方向发展。slope_func返回一个包含两个导数的序列。
在调用run_ode_solver之前，最好使用初始条件测试斜率函数：

dydt, dvdt = slpoe_func(system.init, 0, system)

20_1.png

图20.1: 硬币的高度与时间的关系，没有空气阻力

速度为 0 m/s，加速度为 9.8 m/s^2。现在我们将run_ode_solver称为如下：

results, details = run_ode_solver(system, slope_func)

结果是一个有两列的时间框架(TimeFrame)：y 代表硬币的高度；v 代表它的速度。

我们可以像这样画出结果图：

def plot_position(results):
    plot(results.y)
    decorate(xlabel='Time (s)',
            ylabel='Position (m)')

图20.1显示了结果。由于加速度恒定，速度线性增加，位置二次减少；因此，高度曲线是抛物线。

results.y的最后一个值是-109m，这意味着我们运行模拟的时间太长了。解决这个问题的一种方法是利用结果来估计硬币落在地面上的时间。

ModSim库提供交叉，它接受一个TimeSeries和一个值，并在序列通过该值时返回一个时间序列。我们可以这样找到硬币高度为 0 的时间：

t_crossings = crossings(results.y, 0)

结果是一个只有一个值的数组，8.818s。现在，我们可以用t_end=8.818再次运行模拟，但是还有更好的方法。

20.3 事件

作为一个选项，run_ode_solver可以使用一个event_function，该函数检测到一个“事件”，如硬币击中地面，那就结束模拟。

事件函数采用与斜率函数、state、t和system相同的参数。它们应该返回一个在事件发生时变为0的值。以下是一个事件函数，用于检测硬币撞到地面上：

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

返回值是硬币的高度 y，当硬币碰到地面时，它变为0。

我们把事件函数像这样传递给run_ode_solver：

results, details = run_ode_solver(system, slope_func,
                                  events=event_func)

然后我们可以得到在空中的时间和最终速度，如下所示：

t_sidewalk = get_last_label(results) * s
v_sidewalk = get_last_value(results.v)

在你继续之前，你可能想读一下本章的笔记本，chap20.ipynb，并练习一下。有关下载和运行代码的说明，请参阅Section 0.4。

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

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

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

未经允许，请勿转载。