[SciPy-user] integrate.odeint and event handling

[James]

Chris,

The way I would go about this is more akin to how you would have to use
ODEPACK in Fortran. The odeint function takes a list of output timesteps,
but the Fortran routine is called once for each desired output and uses the
previous call as the initial conditions for the next.

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 )


# ------------------------------
>
>
You could probably integrate the output time conditions into dt_stop by
using a fixed timestep to make it a bit cleaner.

Cheers,
James

On Wed, Jan 21, 2009 at 5:44 PM, Christopher W. MacMinn <cmac at mit.edu>wrote:

> Hello -
>
> I was wondering if integrate.odeint offers any event handling
> capabilities.
>
> For example, say that I want to solve the simple ODE df/dt=f-2 with
> f(t=0)=1.  Also, say that I want to stop the integration when f=0,
> maybe because I don't care about negative values of f, or maybe
> because what I really want to know is the value of t when f=0.  (The
> analytical solution is f(t) = 2-exp(t), and f=0 at t=ln(2).)
>
> The MATLAB code below will produce the desired solution.  It will stop
> integration when f=0, and I believe it will also integrate with some
> care near f=0.  Additionally, MATLAB returns the vector of times at
> which the solution is evaluated, so I can easily grab the value of t
> when f=0.
>
> % ---------------------------------------------------------
> function [ts,fs] = ode_with_events()
>
>    function dfdt = df(t,f)
>        dfdt = f-2.;
>    end
>
>    function [value,isterminal,direction] = df_events(t,f)
>        value = f;
>        isterminal = 1;
>        direction = 0;
>    end
>
>    f0 = 1.;
>    t0 = 0.;
>    t_max = 5.;
>    options = odeset('events', at df_events);
>    [ts,fs] = ode45(@df,[t0,t_max],f0,options);
>
> end
> % ---------------------------------------------------------
>
>
> The Python code below integrates the ODE just fine, but is there a way
> to get the "event" functionality described above?
>
> # ---------------------------------------------------------
> import numpy as np
> from scipy import integrate
>
> def ode_with_events():
>
>    def df(f,t):
>        return f-2.
>
>    f0 = 1.
>    t0 = 0.
>    t_max = 5.
>    ts = np.linspace(t0,t_max,100)
>    fs = integrate.odeint(df,f0,ts)
>
>    return ts,fs
> # ---------------------------------------------------------
>
>
> Thanks!
>
> Best, Chris
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20090121/0c43bcbd/attachment.html>