3
$\begingroup$

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()
$\endgroup$
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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) $\endgroup$ – uhoh Oct 24 '15 at 6:24
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Dipole May 13 '14 at 18:20
  • $\begingroup$ 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 $\endgroup$ – GertVdE May 13 '14 at 18:31
  • $\begingroup$ 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. $\endgroup$ – Dipole May 13 '14 at 18:37
  • $\begingroup$ 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. $\endgroup$ – GertVdE May 13 '14 at 18:45
  • $\begingroup$ 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. $\endgroup$ – Dipole May 13 '14 at 18:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.