2
$\begingroup$

I'm using Python scipy.integrate.ode, and I want to stop my integration at a certain condition. So I use the integrator "dopri5" and use the method "set_solout" to specify a function for my stop condition.

Unfortunately, in some cases, the program says my problem is stiff, and exits. (The message is "UserWarning: dopri5: problem is probably stiff (interrupted)"). So I try another integrator. The documentation mentions 3 other integrators ("vode", "zvode" and "lsoda") which might be suitable for stiff problems. (There is also "dop853", but that exits in the same way as "dopri5").

Here's my problem: none of these integrators support the "set_selout" method. So: is there another way to stop any of these integrators on a condition?

EDIT (a few hours later): I got around this with repeated calls of the "integrate" method (with integrator set to "lsoda") with a short time step; and checking my end condition after each call. Ugly, but gets the job done. I'd still like to know if there's a better way to do this.

$\endgroup$
  • 2
    $\begingroup$ The functionality you're after goes by the name "event detection" or "rootfinding". As far as I can tell, it is not supported in Scipy but is supported by some solvers in pysundials. $\endgroup$ – David Ketcheson Jan 19 '17 at 6:47
2
$\begingroup$

I don't know of a better way to handle this with SciPy since I don't think it has event handling. But if you're willing to venture beyond SciPy, the following software have the capability, either documented as event handling or rootfinding:

  • Sundials' CVODE
  • MATLAB's ode23 and ode15s
  • Julia's DifferentialEquations.jl Sundials wrapper and Rosenbrock methods
  • Mathematica's NDSolve
$\endgroup$
  • $\begingroup$ Thanks. My hack works for now, but if I return to the problem, I'll use Julia or Matlab. (Probably Julia, which I kind of like). $\endgroup$ – Peter B Jan 24 '17 at 5:06
1
$\begingroup$

If you want to stick to the Python scientific family, you could opt for Assimulo. Assimulo is a wrapper around a lot of ODE integrators, providing a common interface. If you happen to be running on Windows, you can find wheels at Christoph Gohlke's page.

Assimulo allows to handle state events to allow for discontinuities in the RHS function of your ODE but you can use them also to stop the iteration when a certain condition is fulfilled.

The procedure is to define a function 'state_events' that analyzes whether an event took place and a second function 'handle_event' to act accordingly (changing the state equations, changing the variables, stopping the iteration, ...).

The example in the Assimulo package is a pendulum that hits a wall. When it hits the wall, the state changes in the sense that the direction should be reversed and some energy should be lost.

I modified the example that comes with Assimulo, using a simple exponential decay problem where I asked to stop the integration when the concentration of the species drops below a threshold value. See the code below.

#!/usr/bin/env python 
# -*- coding: utf-8 -*-
import numpy as N
import pylab as P
from assimulo.problem import Explicit_Problem
from assimulo.solvers.sundials import CVode
from assimulo.exception import TerminateSimulation


def run_example():

    def decay(t,y):

        yd_0 = -0.5*y[0]

        return N.array([yd_0])

    def state_events(t,y,sw):
        """
        This is our function that keep track of our events, when the sign
        of any of the events has changed, we have an event.
        """

        return N.array([y[0]-50.0])


    def handle_event(solver, event_info):
        """
        Event handling. This functions is called when Assimulo finds an event as
        specified by the event functions.
        """
        state_info = event_info[0] #We are only interested in state events info

        if state_info[0] != 0: #Check if the first event function have been triggered
            raise TerminateSimulation


    #Initial values
    y0 = [100.0] #Initial states
    t0 = 0.0             #Initial time

    #Create an Assimulo Problem
    mod = Explicit_Problem(decay, y0, t0)

    mod.state_events = state_events #Sets the state events to the problem
    mod.handle_event = handle_event #Sets the event handling to the problem
    mod.name = 'Simple decay'   #Sets the name of the problem

    #Create an Assimulo solver (CVode)
    sim = CVode(mod)

    #Specifies options 
    sim.discr = 'BDF'     #Sets the discretization method
    sim.iter = 'Newton' #Sets the iteration method
    sim.rtol = 1.e-8        #Sets the relative tolerance
    sim.atol = 1.e-6        #Sets the absolute tolerance

    #Simulation
    ncp = 200     #Number of communication points
    tfinal = 10.0 #Final time

    t,y = sim.simulate(tfinal, ncp) #Simulate

    #Print event information
    sim.print_event_data()

    #Plot
    P.plot(t,y)
    P.show()

if __name__=='__main__':
    run_example()
$\endgroup$
0
$\begingroup$

Event detection in python:

The solution could be https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

From the documentation: ‘RK45’ or ‘RK23’ method for non-stiff problems and ‘Radau’ or ‘BDF’ for stiff problems.

$\endgroup$

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.