Imitate ode45 function from MATLAB in Python

Question

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.

Take a look at [this example](https://stackoverflow.com/questions/16239678/how-to-use-dorpi5-or-dop853-in-python/16240484#16240484) in an old answer of mine. — Warren Weckesser, Jan 24 '18 at 17:47
`vdp1([0,20],[2,0])` is an array, the result of passing two lists to your function. `ode` expects a function, such as `vdp1` itself. In the MATLAB you pass `@vdp1`, not `vdpt1([0 20],[2 0])` to `ode45. — hpaulj, Jan 24 '18 at 18:24

score 9 · Answer 1 · edited Sep 24 '19 at 05:20

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!

Does this do variable step sizes like ode45? – crobar Sep 15 '18 at 20:20 — crobar, Sep 15 '18 at 20:20

score 4 · Accepted Answer · 2018-01-25T18:06:46.703

The interface of integrate.ode is not as intuitive as of a simpler method odeint which, however, does not support choosing an ODE integrator. The main difference is that ode does not run a loop for you; if you need a solution at a bunch of points, you have to say at what points, and compute it one point at a time.

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

Or use the new API with the `odeint`-like interface in [`scipy.integrate.RK45`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html) or [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html). — Lutz Lehmann, Jan 26 '18 at 07:44
@LutzL Please post an answer, which would give this approach more visibility that a comment. I, for one, didn't know about these API. — , Jan 26 '18 at 21:18

score 3 · Answer 3 · edited Apr 17 '22 at 02:00

3

The function scipy.integrate.solve_ivp uses the method RK45 by default, similar the method used by Matlab's function ODE45 as both use the Dormand-Pierce formulas with fourth-order method accuracy.

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

Ordinary Differential Equations (solve_ivp)

edited Apr 17 '22 at 02:00

Lost Fool

19
5

answered Jan 27 '21 at 23:56

David Contreras

169
1
3

1

Note that Matlab's ode45 has a "Refine" option that defaults to 4. That is, per computed point it adds 3 additional interpolated points from the segment to the return vectors. There is no such option in python solve_ivp, so that with the same tolerances leading to approximately the same work ode45 will return more points that give smoother looking plots. – Lutz Lehmann Jan 28 '21 at 12:18

Imitate ode45 function from MATLAB in Python

3 Answers3

Linked