Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Imitate ode45 function from MATLAB in Python

I am wondering how to export MATLAB function ode45 to python. According to the documentation is should be as follows:

 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")

The results are completely different, Matlab returns different dimensions than Python.

like image 590
Migui Mag Avatar asked Jan 24 '18 17:01

Migui Mag


People also ask

How do you write an ODE45 function in Matlab?

[ t , y ] = ode45( odefun , tspan , y0 ) , where tspan = [t0 tf] , integrates the system of differential equations y ' = f ( t , y ) from t0 to tf with initial conditions y0 . Each row in the solution array y corresponds to a value returned in column vector t .

Is ode23 faster than ODE45?

ode23 is a three-stage, third-order, Runge-Kutta method. ode45 is a six-stage, fifth-order, Runge-Kutta method. ode45 does more work per step than ode23, but can take much larger steps.

What integration method does ODE45 use?

ode23 and ode45 are automatic step-size Runge-Kutta-Fehlberg integration methods. ode23 uses a simple second and third order pair of formulas for medium accuracy and ode45 uses a fourth and fifth order pair for higher accuracy.

Does Simulink use ODE45?

The ODE45 is a variable step solver. Simulink solves the problem as a fixed step solver. The ODE45 can replicate the results of Simulink if a limit is imposed on the maximum step that it can take using the odeset function.


1 Answers

As @LutzL mentioned, you can use the newer API, solve_ivp.

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

If t_eval is not specified, then you won't have one record per one timestamp, which is mostly the cases I assume.

Another side note is that for odeint and often other integrators, the output array is a ndarray of a shape of [len(time), len(states)], however for solve_ivp, the output is a list(length of state vector) of 1-dimension ndarray(which length is equal to t_eval).

So you have to merge it if you want the same order. You can do so by:

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

First you need to reshape to make it a [n,1] array, and merge it horizontally. Hope this helps!

like image 155
extraymond Avatar answered Sep 18 '22 10:09

extraymond