# How to get ODE solution at specified time points?

The code below basically illustrates my problem. It is a test code for a pendulum. I solve it using a method suggested on https://stackoverflow.com/questions/12926393/using-adaptive-step-sizes-with-scipy-integrate-ode). Now I want to plot points at a range of evenly spaced times to show the fact the dynamic nature of the pendulum swings. For example I want to get the solution and plot the corresponding solution state space points at the times

[0,0.2,0.4,0.6,0.8,1.0]


Is there any automatic way Python can do this (using interpolation of the internal solution points which are obtained during the ODE integration)?

Here is the main code:

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
import warnings

def f(t,y):
l = 1
m = 1
d = 1
g = 9.8
return [y, -np.sin(y)*g/l-y*d/m]

y0, t0 = [np.pi/2, 0], 0
t1 = 500

backend = 'vode'

solver = ode(f).set_integrator(backend, nsteps=1)
solver.set_initial_value(y0, t0)
# suppress Fortran-printed warning
solver._integrator.iwork = -1

y, t = [], []
warnings.filterwarnings("ignore", category=UserWarning)
while solver.t < t1:
solver.integrate(t1, step=True)
y.append(solver.y)
t.append(solver.t)
warnings.resetwarnings()
y = np.array(y)
t = np.array(t)

plt.plot(y[:,0], y[:,1], 'b-')
plt.plot(y[0::5,0], y[0::5,1], 'b.')

plt.show()


GertVdE's answer covers the user-facing implementation well.

Essentially, you want two different things to happen:

• the integrator should take as large a time step as possible
• you want output at regularly-spaced points in time anyway

Here's how that works in practice:

• The integrator has its own internal time step information. If you were to use VODE in "one-step" mode (calling the Fortran library DVODE, or its successor, CVODE, in SUNDIALS), it would take as large a step as possible (subject to constraints on stability and error). It would then output that step. If you were to do this repeatedly, your output would be the time steps of the integrator. However, that output is not what you want -- you want regularly spaced output.
• The integrator can also interpolate between internal time steps using the polynomial approximation constructed by the integrator. In Hairer and Wanner, this functionality is called "dense output". So when you query output from the integrator (as it generates it in time) between internal time steps, what it is really doing is interpolating between these time steps to get the output.
• The interface GertVdE describes seems to use "normal" mode, which does dense output automatically and takes time steps internally, because most users care about the output they get, and aren't concerned about the numerical details in the implementation, so long as they work correctly.

You can see all of these details in the SUNDIALS Usage Notes and in the various SUNDIALS Users Guides, which cover the theory behind the integrators as well as their use. Although these sources of documentation cover the successors to VODE (CVODE and CVODES), the theory should be the same, and the interfaces should be similar, except that CVODE and CVODES have additional functionality.

• So "dense output" could be used, for example, to re-interpolate later on a new grid of timepoints, or to find "events" such as closest approach. This is not available (yet) in SciPy? So far I've found these, which suggest this might be available in the future? (1), (2), (3), (4)
– uhoh
Oct 24, 2015 at 6:24

The example given in the docs of scipy.integrate.ode hints at using successive integration between your points of interest:

>>> r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True)
>>> r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
>>> t1 = 10
>>> dt = 1
>>> while r.successful() and r.t < t1:
>>>     r.integrate(r.t+dt)
>>>     print("%g %g" % (r.t, r.y))


The scipy.integrate.odeint function takes an input argument 't': A sequence of time points for which to solve for y. The initial value point should be the first element of this sequence.

• I'm not sure if this applies to my case: I don't give a time step to the integrator, please see my code. I just want a way to interpolate my solution for specified times. May 13, 2014 at 18:20
• Jack, you need to integrate from point to point that you are interested in. So the call to integrate takes as argument your times in your vector May 13, 2014 at 18:31
• Thanks for your answer. But I need the integrator to determine the suitable step size, because I don't want a uniform spacing of solutions because my system will change very rapidly in some regions, requiring many small time steps, and in regions where the system changes slowly, larger time steps can be used for efficiency. Please let me know if I have misunderstood something. May 13, 2014 at 18:37
• The integrator will choose a suitable step size based on error control in between the integration boundaries. You are in control of these boundaries and can (ab-)use them to get information at certain moments in time. May 13, 2014 at 18:45
• Ok, so in your example above, choosing a step size of 1 will give me 10 evenly distributed points in time. But lets say that most of the movement of my system happens between time 9 and 10. Then I will have a bunch of useless points from time 0-9 and then very little information (if I were to plot my output say) between time 9-10. May 13, 2014 at 18:51

### New since SciPy 1.0.0 (late 2017):

Another solution is to solve with solve_ivp and use the dense_output option which allows interpolating between solution steps:

import numpy as np
from scipy.integrate import ode, solve_ivp
import matplotlib.pyplot as plt
import warnings

def f(t,y):
l = 1
m = 1
d = 1
g = 9.8
return [y, -np.sin(y)*g/l-y*d/m]

y0, t0 = [np.pi/2, 0], 0
t1 = 100

t = np.linspace(t0,t1,2000)
sol = solve_ivp(f, (t0, t1), y0, dense_output=True)
y_real = sol.y
y_interp = sol.sol(t)

plt.clf()
plt.plot(y_real[0,:],y_real[1,:],'or')
plt.plot(y_interp[0,:],y_interp[1,:],'.-b')
plt.show()