The (brief) documentation for scipy.integrate.ode
says that two methods (dopri5
and dop853
) have stepsize control and dense output. Looking at the examples and the code itself, I can only see a very simple way to get output from an integrator. Namely, it looks like you just step the integrator forward by some fixed dt, get the function value(s) at that time, and repeat.
My problem has pretty variable timescales, so I'd like to just get the values at whatever time steps it needs to evaluate to achieve the required tolerances. That is, early on, things are changing slowly, so the output time steps can be big. But as things get interesting, the output time steps have to be smaller. I don't actually want dense output at equal intervals, I just want the time steps the adaptive function uses.
A related notion (almost the opposite) is "dense output", whereby the steps taken are as large as the stepper cares to take, but the values of the function are interpolated (usually with accuracy comparable to the accuracy of the stepper) to whatever you want. The fortran underlying scipy.integrate.ode
is apparently capable of this, but ode
does not have the interface. odeint
, on the other hand, is based on a different code, and does evidently do dense output. (You can output every time your right-hand-side is called to see when that happens, and see that it has nothing to do with the output times.)
So I could still take advantage of adaptivity, as long as I could decide on the output time steps I want ahead of time. Unfortunately, for my favorite system, I don't even know what the approximate timescales are as functions of time, until I run the integration. So I'll have to combine the idea of taking one integrator step with this notion of dense output.
Apparently, scipy 1.0.0 introduced support for dense output through a new interface. In particular, they recommend moving away from scipy.integrate.odeint
and towards scipy.integrate.solve_ivp
, which as a keyword dense_output
. If set to True
, the returned object has an attribute sol
that you can call with an array of times, which then returns the integrated functions values at those times. That still doesn't solve the problem for this question, but it is useful in many cases.
odeint came first and is uses lsoda from the FORTRAN package odepack to solve ODEs. solve_ivp is a more general solution that lets use decide which integrator to use to solve ODEs. If you define the method param as method='LSODA' then this will use the same integrator as odeint .
y = odeint(model, y0, t) model: Function name that returns derivative values at requested y and t values as dydt = model(y,t) y0: Initial conditions of the differential states. t: Time points at which the solution should be reported.
The SciPy function odeint returns a matrix which has k columns, one for each of the k variables, and a row for each time point specified in t.
Since SciPy 0.13.0,
The intermediate results from the
dopri
family of ODE solvers can now be accessed by asolout
callback function.
import numpy as np from scipy.integrate import ode import matplotlib.pyplot as plt def logistic(t, y, r): return r * y * (1.0 - y) r = .01 t0 = 0 y0 = 1e-5 t1 = 5000.0 backend = 'dopri5' # backend = 'dop853' solver = ode(logistic).set_integrator(backend) sol = [] def solout(t, y): sol.append([t, *y]) solver.set_solout(solout) solver.set_initial_value(y0, t0).set_f_params(r) solver.integrate(t1) sol = np.array(sol) plt.plot(sol[:,0], sol[:,1], 'b.-') plt.show()
Result:
The result seems to be slightly different from Tim D's, although they both use the same backend. I suspect this having to do with FSAL property of dopri5
. In Tim's approach, I think the result k7 from the seventh stage is discarded, so k1 is calculated afresh.
Note: There's a known bug with set_solout not working if you set it after setting initial values. It was fixed as of SciPy 0.17.0.
I've been looking at this to try to get the same result. It turns out you can use a hack to get the step-by-step results by setting nsteps=1 in the ode instantiation. It will generate a UserWarning at every step (this can be caught and suppressed).
import numpy as np from scipy.integrate import ode import matplotlib.pyplot as plt import warnings def logistic(t, y, r): return r * y * (1.0 - y) r = .01 t0 = 0 y0 = 1e-5 t1 = 5000.0 #backend = 'vode' backend = 'dopri5' #backend = 'dop853' solver = ode(logistic).set_integrator(backend, nsteps=1) solver.set_initial_value(y0, t0).set_f_params(r) # suppress Fortran-printed warning solver._integrator.iwork[2] = -1 sol = [] warnings.filterwarnings("ignore", category=UserWarning) while solver.t < t1: solver.integrate(t1, step=True) sol.append([solver.t, solver.y]) warnings.resetwarnings() sol = np.array(sol) plt.plot(sol[:,0], sol[:,1], 'b.-') plt.show()
result:
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With