[SciPy-user] unexpected behavior in odeint
Chris Myers
myers at tc.cornell.edu
Fri Sep 7 14:25:17 CDT 2001
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.
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).
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).
Any insight here would be appreciated.
Chris
==========================================================================
Chris Myers
Cornell Theory Center
--------------------------------------------------------------------------
636 Rhodes Hall email: myers at tc.cornell.edu
Cornell University phone: (607) 255-5894 / fax: (607) 254-8888
Ithaca, NY 14853 http://www.tc.cornell.edu/~myers
--------------------------------------------------------------------------
"To thine own self be blue." - Polonious Funk
==========================================================================
More information about the SciPy-user
mailing list