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

[Ryan Gutenkunst]

Hans Fangohr wrote:
> 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.
 >  <snip>
> 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?
>  <snip>
> 
> Many thanks,
> 
> Hans
> 

Hi Hans,

In our work with biochemical networks we've encountered a very similar 
requirement. (We, however, want to stop the integration, change some 
parameters, than start again from where we left off.) Our solution isn't 
particularly pretty, but it does give the time and variable values when 
the condition is fulfilled exactly.

Consider the integration of dy/dt where y can be a vector. We want to 
know exactly when the condition c(y) = 0 is satisfied. We integrate 
dy/dt from 0 to some time T, requesting a reasonably fine grid of times. 
Then we go back and calculate c(y, t) for each time point in our 
trajectory. If two adjacent time points (t1 and t2) have different signs 
for c(y, t), we know our condition was fulfilled sometime in between 
those times.

Now we can integrate backward from t2 to find the exact time the 
condition was fulfilled. To do so, we make a change of variables from t 
to c. Namely, we integrate the equations defined by:
   [dy/dc = dy/dt / dc/dt, dt/dc = 1 / dc/dt]
from c = c(t2) to c = 0. (Note that this is one more equation that we 
were integrating before.) The initial conditions are [y(t2), 0]. This 
integration terminates with the values [y(tc), t2 - tc] where tc is the 
time when c crossed zero. We can then insert these values into our 
trajectory, and terminate the integration process if we want.

So there are some obvious drawbacks to this:
  1) The condition needs to be expressed in a continuous way, and we 
need to know dc/dt analytically. This usually isn't so bad since it's 
probably simply related to dy/dt. (For example, we could write your r <= 
0 as c = r - 0 then dc/dt = dr/dt which we already know how to calculate.)
  2) We have to integrate all the way out to T first, before we can 
check whether a condition fired. So if N conditions fire, we're doing 
about N times as much integration work.
  3) If the sign of c(y) changes but then changes back before the next 
time point we sampled in the initial integration, we'll miss that firing 
of the condition.

Nevertheless, it's pretty easy to write in pure Python, and it works 
without access to the integrator's internal approximation of the 
function. Doing this exactly and efficiently would require low-level 
access to the guts of the integrator. There do exist integrators that 
can handle this sort of condition-checking internally, but as far as I 
can tell the one odeint wraps isn't among them.

For now, it's not worth it to me to wrap up another package, especially 
since that would be yet another thing for our users to have to install. 
If there's interest in moving scipy over to a more sophisticated 
integrator, I'd be happy to help, but I'm still somewhat of a newbie and 
have no clue wrt FORTRAN, so I couldn't do it alone. I also expect that 
it might break a lot of existing code if we aren't careful with outputs.

Anyways, I hope my explanation of our solution is reasonably clear. If 
folks are interested, I'll clean-up, comment, and post our code. And, of 
course, if anyone has a better solution, I'm dying to hear it. :-)

Cheers,
Ryan

-- 
Ryan Gutenkunst               |
Cornell Dept. of Physics      |       "It is not the mountain
                               |        we conquer but ourselves."
Clark 535 / (607)255-6068     |        -- Sir Edmund Hillary
AIM: JepettoRNG               |
          http://www.physics.cornell.edu/~rgutenkunst/