В пакете scipy.integrate
есть целая пачка функций для интегрирования обыкновенных дифференциальных уравнений: https://docs.scipy.org/doc/scipy/reference/integrate.html#solving-initial-value-problems-for-ode-systems
Универсальная обёртка - функция solve_ivp
Пример: предмет подбрасывают вверх, начальная высота 10 метров, скорость 1 метр в секунду. Решение вычисляется для интервала t \in [0,3]
, используется метод по-умолчанию явный Рунге-Кутты пятого порядка
import scipy.integrate as spi
import numpy as np
def a(x, t):
# Притяжение не зависит ни от высоты, ни от времени
return -9.81
def diff_equation(t, xv):
# Первый параметр функции - время, второй - вектор [высота, скорость]
x, v = xv
# dx/dt = v, dv/dt = a(x,t)
return [v, a(x, t)]
x0 = 10
v0 = 1
# Вызов решателя, интервал 0-3 сек, начальные условия, набор 10 точек t от 0 до 3
solution = spi.solve_ivp(diff_equation, [0, 3], [x0, v0], t_eval=np.linspace(0,3, 10))
if solution.success:
print("t: ", solution.t)
print("x: ", solution.y[0])
print("v: ", solution.y[1])
else:
print("НЕШМОГЛА... ", solution)
Результат
t: [0. 0.33333333 0.66666667 1. 1.33333333 1.66666667
2. 2.33333333 2.66666667 3. ]
x: [ 10. 9.78833333 8.48666667 6.095 2.61333333
-1.95833333 -7.62 -14.37166667 -22.21333333 -31.145 ]
v: [ 1. -2.27 -5.54 -8.81 -12.08 -15.35 -18.62 -21.89 -25.16 -28.43]
Если 10 отсечек вам мало, то поменяйте значение параметра t_eval
Например, t_eval=np.arange(0,3.001, 0.01)
задаёт набор моментов времени от 0
до 3
секунд с шагом 0.01
секунда.
Набор точек в t_eval
не влияет на точность вычислений, это лишь команда решателю для каких отсечек вывести значения.