How to get ODE solution at specified time points?

Question

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[1], -np.sin(y[0])*g/l-y[1]*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[2] = -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()

Answer 1

See also: Python: Grid with step control ODE solver

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.

Answer 2

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.

Answer 3

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[1], -np.sin(y[0])*g/l-y[1]*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()