Мне интересно, как экспортировать функцию 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.
Интерфейс 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()
odeint-подобным интерфейсом в scipy.integrate.RK45 или scipy.integrate.solve_ivp.
- person Lutz Lehmann; 26.01.2018
Как упоминал @LutzL, вы можете использовать более новый API, solve_ivp.
results = solve_ivp(obj_func, t_span, y0, t_eval = time_series)
Если t_eval не указано, то у вас не будет одной записи на одну временную метку, что я предполагаю в большинстве случаев.
Поэтому вам нужно объединить его, если вы хотите тот же порядок. Вы можете сделать это:
Y =results
merged = np.hstack([i.reshape(-1,1) for i in Y.y])
Сначала вам нужно изменить форму, чтобы сделать его массивом [n,1], и объединить его по горизонтали. Надеюсь это поможет!
Функция scipy.integrate.solve_ivp по умолчанию использует метод RK45, аналогичный методу, используемому функцией Matlab ODE45, так как обе используют метод Дорманда-Пирса. формулы с точностью метода четвертого порядка.
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
Обыкновенные дифференциальные уравнения (solve_ivp)
vdp1([0,20],[2,0])— это массив, результат передачи двух списков в вашу функцию.odeожидает функцию, такую какvdp1. В MATLAB вы передаете@vdp1, а неvdpt1([0 20],[2 0])в `ode45. - person hpaulj schedule 24.01.2018