solve_ivp（ODE 求解器）的实现

Question

这是我关于堆栈溢出的第一篇文章。所以我在 SciPy 的文档中遇到了这个 solve_ivp 求解器的例子。正在解决的问题如下：

大炮在撞击时向上发射并伴随终端事件。事件的终端和方向字段是通过猴子修补函数来应用的。这里 y[0] 是位置，y[1] 是速度。弹丸从位置 0 开始，速度 +10。请注意，积分永远不会达到 t=100，因为事件是终端事件。

文档中的代码如下：

>>> def upward_cannon(t, y): return [y[1], -0.5]
>>> def hit_ground(t, y): return y[0]
>>> hit_ground.terminal = True
>>> hit_ground.direction = -1
>>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)
>>> print(sol.t_events)
[array([40.])]
>>> print(sol.t)
[0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]

我已经使用这个求解器来求解其他微分方程。此外，我了解 terminal 和 direction 字段的用法。但是仅对于这个示例，我无法理解函数 upward_cannon() 是如何使用return [y[1],-0.5] 工作的。对应于print(sol.y) 的输出如下：

[[ 0.00000000e+00  9.99897510e-04  1.09985977e-02  1.10958105e-01
   1.10780373e+00  1.08013149e+01  8.02419261e+01 -1.42108547e-14]
 [ 1.00000000e+01  9.99995000e+00  9.99945005e+00  9.99445055e+00
   9.94445555e+00  9.44450555e+00  4.44500550e+00 -1.00000000e+01]]

由于代码中没有提到底层微分方程，求解器如何为y 生成上述值？我知道return [y[1],-0.5] 正在做某事，但我无法解释它在做什么。

Answer 1

你写的

由于代码中没有提到底层微分方程......

没有明确提到，但其实有一个微分方程。设 h(t) 为炮弹的高度。微分方程为

h''(t) = -0.5

也就是说，加速度是恒定的并且是向下的。这是Newton's second law，假设重力恒定。本例假设万有引力常数与质量之比为0.5。

初始条件为

h(0) = 0,  h'(0) = 10

方程 h''(t) = -0.5 是一个（微不足道的！）二阶微分方程。我们可以用普通的定积分来解决它，但为了示例，我们使用 ODE 求解器。要使用solve_ivp（或 SciPy 中的任何其他 ODE 求解器），我们必须将二阶 DE 转换为一阶方程组。设 y0(t) = h(t) 和 y1(t) = h'(t)。那么

y0'(t) = h'(t) = y1(t)
y1'(t) = h''(t) = -0.5

那是在

中实现的系统

def upward_cannon(t, y):
    return [y[1], -0.5]