[SciPy-user] unexpected behavior in odeint

[Chris Myers]

Travis,

Included below is a very simple example demonstrating the problem 
I reported (although this is a different example than the specific 
one I described in my earlier email).  The problem is actually more
prevalent in this example.

As an aside, it would be useful if the 'tcur' array returned in
the infodict had the initial time value at the start of the array;
as it is now, the array starts with the first successful time step,
and therefore the length of the solution array is one longer than the
tcur array (making it more awkward, e.g., to plot y(t) vs. t).

Thanks,

Chris


# -------------------------------------------------------------------------
import Numeric
from scipy.integrate import odeint

y0 = 1.

def f(y,t):
    return Numeric.exp(-0.01*y*t)

t = Numeric.arange(0, 3.0, 0.1)
y, yinfo = odeint(f, y0, t, full_output=1)

# -------------------------------------------------------------------------

For this example, I find:

yinfo['tcur'] = 
array([ 0.11397125,  0.20087814,  0.54850571,  0.54850571,  0.54850571,  0.83819534,
             0.83819534,  0.83819534,  1.12788498,  1.12788498,  1.12788498,
             1.41757461,  1.41757461,  1.41757461,  1.70726425,  1.70726425,
             1.70726425,  1.99695389,  1.99695389,  2.38792355,  2.38792355,
             2.38792355,  2.38792355,  2.77889321,  2.77889321,  2.77889321,
             2.77889321,  3.16986287,  3.16986287])

y = 
array([[ 1.        ],
       [ 1.09994666],
       [ 1.1997735 ],
       [ 1.29946076],
       [ 1.39898883],
       [ 1.49833825],
       [ 1.5974897 ],
       [ 1.69642407],
       [ 1.79512247],
       [ 1.89356621],
       [ 1.99173682],
       [ 2.08961611],
       [ 2.18718616],
       [ 2.28442932],
       [ 2.38132824],
       [ 2.47786589],
       [ 2.57402557],
       [ 2.66979092],
       [ 2.76514591],
       [ 2.8600749 ],
       [ 2.95456261],
       [ 3.04859415],
       [ 3.14215499],
       [ 3.23523105],
       [ 3.32780862],
       [ 3.41987444],
       [ 3.51141562],
       [ 3.60241973],
       [ 3.69287475],
       [ 3.7827691 ]])


On Fri, Sep 07, 2001 at 02:17:42PM -0600, Travis Oliphant wrote:
> 
> >
> > I have a question about the intended behavior of the
> > scipy.integrate.odeint function.
> >
> > The documented process for constructing an ODE solution using the
> > odeint() function is to pass a sequence (NumPy array) of times t at
> > which a solution y(t) is desired.  Looking into the underlying lsoda
> > subroutine, and the __odepack.h wrapper, however, I see that this
> > multiple time-step interface is provided in the Python wrapper; that
> > is, lsoda would appear to be constructed to execute a single time step
> > (over and over again), and the multiple time stepping is a convenience
> > provided in the Python overlayer.
> >
> 
> This is true....
> 
> > The problem I am having is that lsoda does not take a time step
> > exactly as requested, but instead takes a step to ensure that the
> > current time is at least as large as the requested time value.  As a
> > result, in the process of feeding an array of requested time points,
> > lsoda can "outrun" the driver and produce multiple steps with the same
> > value of the time tcur.  For example, I may request a solution at times
> > t=arange(0., 10., 0.1).  odeint guarantees that each step will be at
> > least as large as the requested time, but in some cases it would
> > appear that the current time actually already exceeds the requested
> > time.  In a particular example of mine, I find:
> >
> 
> 
> > t = t_requested[65:68] = array([ 6.5,         6.6,         6.7])
> > tcur = t_actual[65:68] = array([ 6.56972519,  6.70778359,  6.70778359])
> >
> > What we see is that on step 66, odeint has incremented itself to a
> > time which exceeds the requested value for step 67, so that it does
> > not step to a larger time during step 67.  odeint does, however,
> > return a different function value at steps 66 and 67, so that the
> > result function y(t) is multivalued for some t (which is not especially
> > appealing to me).
> 
> The problem that is happening is that I did not copy the array t before
> passing it to the _odepack.odeint C-interface.  This interface hands the
> array t to the lsoda routine which causes your t array to be overwritten
> with possible overshot values for t0.
> 
> My understanding, though is that lsoda uses interpolation to give the
> output value at the requested time.
> 
> So, I think your y output vector corresponds to the correct values but
> your t vector may have been overwritten.  What does t_original compared to
> the ouput look like?
> 
> >
> > Am I misusing this multiple time-step feature of odeint?  Is the
> > behavior that I have describe intended?  Should I be exercising caution
> > in the construction of time-sequences to pass to odeint?  An alternate
> > strategy would be to take the current time returned by lsoda, and
> > set the next requested time based on that (but in such a case, one
> > might get fewer time slices than originally requested).
> 
> I had not intended for the t vector to be overwritten (i don't think it is
> a problem all the time -- but you've uncovered a case).
> 
> Is it possible for you to send an example fprime function, intial
> conditions and time-points so I could try to replicate this behavior.
> 
> Thanks,
> 
> -Travis
> 
> 
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user
>