Stepsize control of dopri5 integrator
Question: I am trying to solve a simple example with the `dopri5` integrator in
`scipy.integrate.ode`. As the documentation states

> This is an explicit runge-kutta method of order (4)5 due to Dormand & Prince
> (with stepsize control and dense output).

this should work. So here is my example:

    
    
    import numpy as np
    from scipy.integrate import ode
    import matplotlib.pyplot as plt
    
    
    def MassSpring_with_force(t, state):
        """ Simple 1DOF dynamics model: m ddx(t) + k x(t) = f(t)"""
        # unpack the state vector
        x = state[0]
        xd = state[1]
    
        # these are our constants
        k = 2.5 # Newtons per metre
        m = 1.5 # Kilograms
    
        # force
        f = force(t)
    
        # compute acceleration xdd
        xdd = ( ( -k*x + f) / m )
    
        # return the two state derivatives
        return [xd, xdd]
    
    
    def force(t):
        """ Excitation force """
        f0 = 1  # force amplitude [N]
        freq = 20  # frequency[Hz]
        omega = 2 * np.pi *freq  # angular frequency [rad/s]
        return f0 * np.sin(omega*t)
    
    
    # Time range
    t_start = 0
    t_final = 1
    
    # Main program
    state_ode_f = ode(MassSpring_with_force)
    
    state_ode_f.set_integrator('dopri5', rtol=1e-6, nsteps=500,
                           first_step=1e-6, max_step=1e-3)
    
    state2 = [0.0, 0.0]  # initial conditions
    state_ode_f.set_initial_value(state2, 0)
    
    sol = np.array([[t_start, state2[0], state2[1]]], dtype=float)
    
    print("Time\t\t Timestep\t dx\t\t ddx\t\t state_ode_f.successful()")
    
    while state_ode_f.t < (t_final):
        state_ode_f.integrate(t_final, step=True)
        sol = np.append(sol, [[state_ode_f.t, state_ode_f.y[0], state_ode_f.y[1]]], axis=0)
        print("{0:0.8f}\t {1:0.4e} \t{2:10.3e}\t {3:0.3e}\t {4}".format(
                state_ode_f.t, sol[-1, 0]- sol[-2, 0], state_ode_f.y[0], state_ode_f.y[1], state_ode_f.successful()))
    

The result I get is:

    
    
    Time         Timestep    dx      ddx         state_ode_f.successful()
    0.49763822   4.9764e-01      2.475e-03   -8.258e-04  False
    0.99863822   5.0100e-01      3.955e-03   -3.754e-03  False
    1.00000000   1.3618e-03      3.950e-03   -3.840e-03  False
    

with a warning:

> c:\python34\lib\site-packages\scipy\integrate_ode.py:1018: UserWarning:
> dopri5: larger nmax is needed self.messages.get(idid, 'Unexpected idid=%s' %
> idid))

The result is incorect. If I run the same code with `vode` integrator, I get
the expected result.

**Edit**

A similar issue is described here: [Using adaptive step sizes with
scipy.integrate.ode](http://stackoverflow.com/questions/12926393/using-
adaptive-step-sizes-with-scipy-integrate-ode)

The suggested solution recommends setting `nsteps=1`, which solves the ODE
correctly and with step-size control. However the integrator returns
`state_ode_f.successful()` as `False`.


Answer: No, there is nothing wrong. You are telling the integrator to perform an
integration step to `t_final` and it performs that step. Internal steps of the
integrator are not reported.

* * *

The sensible thing to do is to give the desired sampling points as input of
the algorithm, set for example `dt=0.1` and use

    
    
    state_ode_f.integrate( min(state_ode_f.t+dt, t_final) )
    

There is no single-`step` method in `dopri5`, only `vode` has it defined, see
the source code
<https://github.com/scipy/scipy/blob/v0.14.0/scipy/integrate/_ode.py#L376>,
this could account for the observed differences.

As you found in [Using adaptive step sizes with
scipy.integrate.ode](http://stackoverflow.com/questions/12926393/using-
adaptive-step-sizes-with-scipy-integrate-ode), one can force single-step
behavior by setting the iteration bound `nsteps=1`. This will produce a
warning every time, so one has to suppress these specific warnings to see a
sensible result.

* * *

You should not use a parameter (which is a constant for the integration
interval) for a time-dependent force. Use inside `MassSpring_with_force` the
evaluation `f=force(t)`. Possibly you could pass the function handle of
`force` as parameter.

