suggestions for improving scipy.integrate.odeint

Dear all, dear Travis (author of odeint?), let me first say that scipy is a very useful extension to Python and that the combination of Python, scipy (and Numerix) and -- in our case -- pylab, rivals MATLAB for our purposes. Over the last few months, I have been using scipy and in particular scipy.integrate.odeint in a introductory lecture course on numeric computation for some engineering degrees. Thanks to the creativity of the students (they have the most fantastic ideas you could not dream of) I have come across two minor issues with odeint, which I'd like to point out here. If these could be resolved (in the long term), then this would make odeint much easier to use for beginners. I'll explain why as I go along. 1. The initial value (array) can change while odeint carries out the integration. This source code #---------------------------- import Numeric from scipy.integrate import odeint def rhs(y,t): return -y[0] y0 = Numeric.array([1.0]) print "Initial value:",y0 ts = Numeric.arange(0,10,1) ys = odeint( rhs, y0, ts ) print "Initial value:",y0 #---------------------------- produces the following output: #=========================== Initial value: [ 1.] Initial value: [ 0.00012341] #=========================== I can kind-of-see why this may happen (without digging in the source). However, this is highly confusing to the students because y0 is an input argument to the odeint function and should not change when odeint is called. (There is, of course, the issue of passing mutable objects such as lists to functions but I'd rather not start explaining that to non-computer science students in their very first lectures on Python. Also, if y0 is a list (rather than an array), this problems does not exist: if y0 is a list, then it is not modified by running odeint!) 2. The sequence of values t for which odeint returns the numerical solution y(t) cannot be a list. Here is an example: #---------------------------- import Numeric from scipy.integrate import odeint def rhs(y,t): return -y[0] y0 = Numeric.array([1.0]) print "this works" ts = Numeric.arange(0,10,1) ys = odeint( rhs, y0, ts ) print "this doesn't" ts = range(10) ys = odeint( rhs, y0, ts ) #---------------------------- The output of this code is: #=========================== this works this doesn't Traceback (most recent call last): File "<stdin>", line 15, in ? File "/usr/lib/python2.3/site-packages/scipy/integrate/odepack.py", line 114, in odeint t = t.copy() AttributeError: 'list' object has no attribute 'copy' #=========================== Since y0 can be a list, it would be consistent to accept ts to be a list, too. (This is in fact very useful: by posing problems carefully, students can use odeint 'from lecture one' even if they haven't learned much about data types yet.) I hope this is useful information. Thanks again for providing such a nice wrapper around ODEPACK. Cheers, Hans

Dear Hans, On Tue, 19 Apr 2005, Hans Fangohr wrote: [...]
1. The initial value (array) can change while odeint carries out the integration. This source code
#---------------------------- import Numeric from scipy.integrate import odeint
def rhs(y,t): return -y[0]
y0 = Numeric.array([1.0]) print "Initial value:",y0 ts = Numeric.arange(0,10,1) ys = odeint( rhs, y0, ts ) print "Initial value:",y0 #----------------------------
produces the following output: #=========================== Initial value: [ 1.] Initial value: [ 0.00012341] #===========================
I can kind-of-see why this may happen (without digging in the source). However, this is highly confusing to the students because y0 is an input argument to the odeint function and should not change when odeint is called.
I think this is actually done on purpose, as y0 contains the last entry of ys, print ys[-1] [ 0.00012341] This allows to restart odeint with this initial value to continue the integration. We stumbled across the same some time ago (I thought this was documented somewhere, but it seems it is not, at least not in the doc-string to odeint). Best, Arnd

Hi Arnd,
On Tue, 19 Apr 2005, Hans Fangohr wrote:
[...]
1. The initial value (array) can change while odeint carries out the integration. This source code
#---------------------------- import Numeric from scipy.integrate import odeint
def rhs(y,t): return -y[0]
y0 = Numeric.array([1.0]) print "Initial value:",y0 ts = Numeric.arange(0,10,1) ys = odeint( rhs, y0, ts ) print "Initial value:",y0 #----------------------------
produces the following output: #=========================== Initial value: [ 1.] Initial value: [ 0.00012341] #===========================
I can kind-of-see why this may happen (without digging in the source). However, this is highly confusing to the students because y0 is an input argument to the odeint function and should not change when odeint is called.
I think this is actually done on purpose, as y0 contains the last entry of ys, print ys[-1] [ 0.00012341] This allows to restart odeint with this initial value to continue the integration.
I can see what you are saying. However, I would argue that this 'obscure' feature of returning the last value of ys (indirectly) in y0 is not 'pythonic' because it is not clear that this would happen. Note also, that y0 does not change if the type of y0 is a list (!) rather than an Numeric.array (as in my example). If the intention is that y0 is always the same as ys[-1] then this should be consistent independent of the type of y0 (although --- as I said before --- I can't see being a good idea yet).
We stumbled across the same some time ago (I thought this was documented somewhere, but it seems it is not, at least not in the doc-string to odeint). Correct :)
Thanks for the feedback, Hans
participants (2)
-
Arnd Baecker
-
Hans Fangohr