Имитация функции ode45 из MATLAB в Python

Question

Мне интересно, как экспортировать функцию MATLAB ode45 в Python. Согласно документации это должно быть следующим:

 MATLAB:  [t,y]=ode45(@vdp1,[0 20],[2 0]);

 Python:  import numpy as np
          def  vdp1(t,y):
              dydt= np.array([y[1], (1-y[0]**2)*y[1]-y[0]])
              return dydt
          import scipy integrate 
          l=scipy.integrate.ode(vdp1([0,20],[2,0])).set_integrator("dopri5")

Результаты совершенно разные, Matlab возвращает разные измерения, чем Python.

Answer 1

Интерфейс integrate.ode не такой интуитивно понятный, как у простого метода odeint, который, однако, не поддерживает выбор интегратора ODE. Основное отличие состоит в том, что ode не запускает цикл для вас; если вам нужно решение по множеству точек, вы должны сказать, в каких точках, и вычислять его по одной точке за раз.

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

def vdp1(t, y):
    return np.array([y[1], (1 - y[0]**2)*y[1] - y[0]])
t0, t1 = 0, 20                # start and end
t = np.linspace(t0, t1, 100)  # the points of evaluation of solution
y0 = [2, 0]                   # initial value
y = np.zeros((len(t), len(y0)))   # array for solution
y[0, :] = y0
r = integrate.ode(vdp1).set_integrator("dopri5")  # choice of method
r.set_initial_value(y0, t0)   # initial values
for i in range(1, t.size):
   y[i, :] = r.integrate(t[i]) # get one more value, add it to the array
   if not r.successful():
       raise RuntimeError("Could not integrate")
plt.plot(t, y)
plt.show()

Answer 2

Как уже упоминалось @LutzL, вы можете использовать новый API.

results = solve_ivp(obj_func, t_span, y0, t_eval = time_series)

Если t_eval не указывается, тогда у вас не будет одной записи на одну временную метку, что в большинстве случаев я предполагаю.

Другое примечание заключается в том, что для odeint и часто другие интеграторы, выходной массив является ndarray формы [len(time), len(states)]Однако для solve_ivpвыходной list(length of state vector) 1-мерного ndarray(длина которого равна t_eval).

Таким образом, вы должны объединить его, если вы хотите такой же порядок. Вы можете сделать это:

Y =results
merged = np.hstack([i.reshape(-1,1) for i in Y.y])

Сначала вам нужно изменить форму, чтобы сделать его [n,1] массив, и слить его по горизонтали. Надеюсь это поможет!

Answer 3

Функция scipy.integrate.solve_ivp по умолчанию использует метод RK45 , аналогичный методу, используемому функцией ODE45 в Matlab, поскольку оба используют формулы Дорманда-Пирса с точностью метода четвертого порядка.

      vdp1 = @(T,Y) [Y(2); (1 - Y(1)^2) * Y(2) - Y(1)];
[T,Y] = ode45 (vdp1, [0, 20], [2, 0]);

      from scipy.integrate import solve_ivp

vdp1 = lambda T,Y: [Y[1], (1 - Y[0]**2) * Y[1] - Y[0]]
sol = solve_ivp (vdp1, [0, 20], [2, 0])

T = sol.t
Y = sol.y

Обыкновенные дифференциальные уравнения (resolve_ivp)