[SciPy-user] stopping criterion for integrate.odeint?
rng7 at cornell.edu
Wed Mar 23 18:22:09 CST 2005
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.
> 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?
> Many thanks,
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
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. :-)
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 |
More information about the SciPy-user