python - 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

你写的

Since there is no underlying differential equation mentioned in the code...

没有明确提及，但实际上，存在一个微分方程。设 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]