[SciPy-user] Problems with odeint

[Víctor Martínez-Moll]

Hi again,

Sorry for bothering again but as my message arrived just before the list 
problem and I've received no answer I'm wondering if maybe some message 
was lost.

So, again, any suggestions on how to solve this simple 2nd order 
differential equation with Scipy? Or, what the hell I did wrong to get 
such strangre results?

Thanks in advance and best regards.

Victor

PS: If you think this is not the right forum or you think I could find 
the answer to my question reading any particular documentation please 
let me know too. I've read some of the documentation I found online but 
maybe I've missed something.

En/na Víctor Martínez-Moll ha escrit:
> 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