[SciPy-user] stopping criterion for integrate.odeint?

[Hans Fangohr]

Dear all,

we have been using odeint in a number of situations and it works
great. Now we'd like to interrupt the integration of the set of ODEs
once a conditions is fulfilled that depends on the solution that is
being computed.

Here is an example. Suppose we want to solve the ODEs associated with an 
object being thrown straight into the sky (i.e. against the gravitational 
force). Assume we can describe this by

r''(t) = -g

where r'' is the second derivative of the (scalar) position with
respect to time.

We convert this into two 1st order ODEs

v' = -g  and
r' = v

and use odeint to solve these. (So far, so good, please find attached
below the complete program demonstrating this).

Here is the complication: Assume we want to interrupt the integration
if r<=0, e.g. when the object hits the ground. Is there an elegant way
of doing this?

The source code below should demonstrate what we try to do.

A few more comments:

(i) the actual example I am trying to solve is more complicated (yes, I 
can compute the time when the object hits the ground analytically in this 
case ;-) but that's not the point here)

(ii) options I considered include:

(a) integrating the set of ODEs in small steps dt by calling odeint 
repeatedly. Feasible but not beautiful as the final solution will depend 
on the size of this dt.

(b) raising an exception in the rhs function (doesn't work)

(c) setting all derivatives to zero in the rhs function (seems to do what 
I want since r and v don't change anymore, but is wasting CPU cycles 
because the integration doesn't stop).

Okay, I hope I outlined by problem. Doe anyone have any advice?

Many thanks,

Hans


------------------------------------------------------------------------------------------------------------
import Numeric as N
from scipy.integrate import odeint

def rhs(y,t):
     """expect y = N.array([r,v])
     Right-Hand-Side function to solve 2nd order ODE r''= -g
     (with g being a constant).

     returns dr/dt = v
             dv/dt = a = F/m = -g and

     """

     r,v = y

     drdt = v
     dvdt = -9.81 #m/s^2

     # this is what I don't know how to do elegantly
     if r<=0:
         print "We'd like to stop integrating now"

     return N.array([drdt,dvdt])

y = odeint(rhs, N.array([0,10]), N.arrayrange(0,3,0.1))


try:
     import pylab
     pylab.plot(y[:,0])
     pylab.show()
except ImportError:
     pass                                 #can't import matplotlib