[SciPy-user] integrate.odeint and event handling

Rob Clewley rob.clewley@gmail....
Wed Jan 21 14:33:26 CST 2009


Let's be clear about the expected functionality of this posted code...

> So your example would read something like... (if not tested this btw...)
>
> # ------------------------------
> import numpy as np
> from scipy import integrate
>
> def df(f,t):
>    return f-2.
>
> def df_stop(f,t):
>    return f < 0.0
>
> f0 = 1.
> t0 = 0.
> t_max = 5.
> nout = 100
> ts = np.linspace(t0,t_max,nout)
>
>
> fs = [f0,]
> df_continue = True
> i = 0
> while df_continue:
>     f = integrate.odeint(df,fs[i],[ts[i],ts[i+1]])
>     i+=1
>     if i==nout-1:
>         df_continue = False
>     elif df_stop(f[1][0],ts[i+1]):
>         df_continue = False
>     else:
>         fs.append( f[1][0] )
>
> fs = np.array( fs )
>
>

This won't stop integration at the actual time that the event occurred
(the OP said he wants to stop when f=0 and I am assuming he means to
some significant accuracy) - it only stops at some time after the
event occurred, up to an error of the fixed step size. The whole point
of the lsodar and pydstool routines is to be able to have an
integration that stops precisely when an event occurs, up to a
predetermined error tolerance. In this code, you would have to
re-integrate between the last two time points (the one before and the
one after the event) at much smaller time steps to discover where the
event is more accurately. This is efficiently done in the other codes.

-Rob


More information about the SciPy-user mailing list