[SciPy-user] Problems with odeint

Víctor Martínez-Moll victor.martinez at uib.es
Fri Jun 9 04:04:59 CDT 2006


Hi all,

I've been a SciLab user for some time and I'm evaluating SciPy as a 
development tool.

The first thing I tried is to solve a simple second order diferential 
equation using odeint(). The problem is that depending on the function I 
want to integrate I get nice results, but for most of them I get simply 
nothing or nonsense answers. Is not a problem of the function having a 
strange behaviour or having singularity points. For example if I try to 
solve:
     d2y/dt2 = 1-sin(y)
either I get nothing or wrong solutions (the best thing I got was 
setting:hmin=0.01,atol=.001), while If I do about the same procedure in 
SciLab I get a nice and smooth set of curves. The strangest thing is 
that if I use exactly the same procedure to solve:
     d2y/dt2 = 1-y
then I get the right solution, which seems to indicate that I'm doing 
the right thing (although of course I know I'm not because I do not 
belive that odeint is not able to solve such a silly thing).

I've only checked it with the last enthon distribution I found: 
enthon-python2.4-1.0.0.beta2.exe

The simple procedure I wrote in Python and its equivalent in SciLab that 
does the right thing in are:

######################################
## The Python one                   ##
######################################

from scipy import *
from matplotlib.pylab import *

def dwdt(w,t):
     return [w[1],1.0-sin(w[0])]

t=arange(0.0,2.0*pi,.01)

ww = integrate.odeint(dwdt,[0.0,0.0],t,hmin=0.01,atol=.001)

y = ww[:,0]
dy =ww[:,1]
ddy = 1.0-sin(ww[:,0])

plot(t,y,label='y')
plot(t,dy,label='dy')
plot(t,ddy,label='ddy')
legend()
show()

######################################
## The SciLab one                   ##
######################################

function result = dwdt(t,w)
   result(1) = w(2)
   result(2) = 1.0-sin(w(1))
endfunction

t = 0.0:0.01:2.0*%pi;

ww = ode([0.0;0.0],0.0,t,dwdt);

y = ww(1,:);
dy = ww(2,:);
ddy=1.0-sin(ww(1,:));

xset("window", 0)
xbasc()
plot2d(t,[y;dy;ddy]',leg="y at dy@ddy")

########################################

Any ideas or suggestions will be wellcome.


-- 
Víctor Martínez Moll           | Universitat de les Illes Balears
Departament de Física          | Edifici Mateu Orfila
Àrea d'Enginyeria Mecànica     | E-07122, Palma de Mallorca, SPAIN
e-mail: victor.martinez at uib.es | Tel:34-971171374 Fax:34-971173426



More information about the SciPy-user mailing list