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

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.

possible duplicate (identified by blaz) in http://stackoverflow.com/questions/12926393/using-adaptive-step-sizes-with-scipy-integrate-ode — Lutz Lehmann, Oct 03 '15 at 18:39

score 0 · Answer 1 · edited May 23 '17 at 11:52

0

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, 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.

edited May 23 '17 at 11:52

Community

1
1

answered Oct 03 '15 at 11:29

Lutz Lehmann

25,219
2
22
51

Hey LutzL thank you for your anwser. I set the evaluation of `f=force(t)`inside the `MassSpring_with_force`. However if you try running the code it still does not work. If you chose i.e. the `vode` integrator, the integration process runs as expected. – blaz Oct 03 '15 at 14:47
Set `rtol=1e-12`, it is possible that the indicated stepsize is sufficient for `rtol=1e-6` in the order 5 dopri5 method. – Lutz Lehmann Oct 03 '15 at 15:27
I still get three iterations. – blaz Oct 03 '15 at 17:03
Try increasing `nsteps`, as that is the complaint in the warning. However, 500 inner iteration steps would relate to a too small step size for any tame problem. – Lutz Lehmann Oct 03 '15 at 17:20
Here is a similar problem: http://stackoverflow.com/questions/12926393/using-adaptive-step-sizes-with-scipy-integrate-ode – blaz Oct 03 '15 at 17:48
setting the `dt` as you recommended works, but this makes a constant step-size – blaz Oct 03 '15 at 17:51
1

And does the method in the similar question of setting `nsteps=1` and suppressing the warnings work for you? – Lutz Lehmann Oct 03 '15 at 18:38
Yes, this works. However, the `state_ode_f.successful()` is always `False`. – blaz Oct 04 '15 at 08:24
1

Yes of course it is false. You are preventing the algorithm to reach the requested end time `t_final`, it is a hack, an inefficient one as it is with all the back'n'forth from python to fortran. But if it works for you, it is fine. – Lutz Lehmann Oct 04 '15 at 10:27

Stepsize control of dopri5 integrator

1 Answers1