7

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.

Improve this question
2
  • Take a look at this example in an old answer of mine. Commented Jan 24, 2018 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. Commented Jan 24, 2018 at 18:24

3 Answers 3

Reset to default
9

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!

Improve this answer
Sign up to request clarification or add additional context in comments.

2 Comments

Thx for cleanup.
Does this do variable step sizes like ode45?
4

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

solution

Improve this answer

2 Comments

Or use the new API with the odeint-like interface in scipy.integrate.RK45 or scipy.integrate.solve_ivp.
@LutzL Please post an answer, which would give this approach more visibility that a comment. I, for one, didn't know about these API.
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)

Improve this answer

2 Comments

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.
Nevertheless, I am still obtaining different results in MATLAB using ode45 with reltol = 1e-6 and in Python solve_ivp using RK45 with reltol = 1e-6. Do you have any clues about why this might be happening?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.