Реализация solve_ivp (ODE решатель)

Question

Это мой первый пост о переполнении стека. Так что я наткнулся на этот пример решателя solve_ivp в документации SciPy. Решаемая проблема заключается в следующем:

Пушка стреляла вверх с конечным событием при ударе. Поля терминала и направления события применяются путем исправления функции. Здесь 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]

Я использовал этот решатель для решения других дифференциальные уравнения. Кроме того, я понимаю использование полей терминал и направление . Но только для этого примера я не могу понять, как работает функция 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

То есть ускорение постоянное и в направлении вниз. Это второй закон Ньютона с предположением о постоянной гравитационной силе. В примере предполагается, что отношение гравитационной постоянной и массы равно 0,5.

Начальные условия

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

Уравнение h '' (t) = -0,5 представляет собой (тривиально !) дифференциальное уравнение второго порядка. Мы могли бы решить это с помощью простой старой определенной интеграции, но для примера мы используем решатель ODE. Чтобы использовать solve_ivp (или любой из других решателей ODE в SciPy), мы должны преобразовать 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]