[SciPy-user] stopping criterion for integrate.odeint?
Fernando.Perez at colorado.edu
Wed Mar 23 18:34:53 CST 2005
Ryan Gutenkunst wrote:
> 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.
In case anyone wants a full reference:
- M. Henon, Physica 5D, 412 (1982).
This is a standard trick in the construction of Poincare maps in dynamical
systems. The above paper by Henon is the place where I first read about it,
though it is something simple enough that it has probably been 'reinvented'
numerous times and in many different contexts.
More information about the SciPy-user