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()
