[SciPy-user] suggestions for improving scipy.integrate.odeint

[Hans Fangohr]

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