CODEWITHSUNDEEP
pythonnumpyscipyode
blaz
blaz

Reputation: 4258

Stepsize control of dopri5 integrator

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.

Upvotes: 2

Views: 2046

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

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.

Upvotes: 0

Related Questions