[SciPy-user] stopping criterion for integrate.odeint?
Hans Fangohr
H.FANGOHR at soton.ac.uk
Wed Mar 23 17:28:53 CST 2005
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
More information about the SciPy-user
mailing list