Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > Python > need help with ODE solver from scipy

Reply
Thread Tools

need help with ODE solver from scipy

 
 
T.Crane
Guest
Posts: n/a
 
      04-24-2007
Hi,

OK, I'm trying to figure out how to use the ODE solver
(scipy.integrate.ode.ode). Here's what I'm doing (in iPython)

y0 = [0,1,1]
dt = 0.01
tEnd = 12
t0 = 0
Y = zeros([tEnd/dt, 3])

As an aside, I've used this assignment for Y in the past and it
worked. When I tried it this morning I got a TypeError and a message
saying I needed to use an integer. So, instead...

Y = zeros([int(tEnd/dt), 3])
T = zero([int(tEnd/dt)])
index = 0
def foo(t,y):
dy = zeros([3])
dy[0] = y[1]*y[2]
dy[1] = -y[0]*y[2]
dy[2] = -0.51*y[0]*y[1]
return dy

r = ode(foo).set_integrator('vode').set_initial_value( y0,t0)
while r.successful() and r.t < tEnd:
r.integrate(r.t+dt)
Y[index] = r.y
T[index] = r.t
index += 1

As a further aside, this system of three coupled linear ODE's is from
an example in MATLAB's documentation on their function ode45, so I
know what I 'should' get. The while loop and call to ode is adapted
from scipy's (very limited) documentation on scipy.integrate.ode.ode.

At the end of this while loop I've gotten a few different results
depending on what I have no idea. This morning another TypeError
exception was thrown, saying:

TypeError: Array can not be safely cast to required type.

Although, to be fair this is after the output from one iteration:

array([0.01, 1. , 1. ])

So, clearly this isn't working right. Does anyone have any experience
using this function or anything they can contribute?

thanks,
trevis

 
Reply With Quote
 
 
 
 
Robert Kern
Guest
Posts: n/a
 
      04-24-2007
T.Crane wrote:
> Hi,
>
> OK, I'm trying to figure out how to use the ODE solver
> (scipy.integrate.ode.ode).


You will get more help from the scipy-user list than you will here.

http://www.scipy.org/Mailing_Lists

> Here's what I'm doing (in iPython)
>
> y0 = [0,1,1]
> dt = 0.01
> tEnd = 12
> t0 = 0
> Y = zeros([tEnd/dt, 3])
>
> As an aside, I've used this assignment for Y in the past and it
> worked. When I tried it this morning I got a TypeError and a message
> saying I needed to use an integer.


Well, the current behavior is correct. When you give us a float instead of an
int for a dimension, we shouldn't guess what you want. Particularly, in this
case, truncating is incorrect. I ran your code and eventually wound up with an
IndexError because your iteration ran past the end of Y and T.

> So, instead...
>
> Y = zeros([int(tEnd/dt), 3])
> T = zero([int(tEnd/dt)])


Also, a better way to do this is simply to leave these as lists and simply
append to them. That way, you don't have to manage indices.

> index = 0
> def foo(t,y):
> dy = zeros([3])
> dy[0] = y[1]*y[2]
> dy[1] = -y[0]*y[2]
> dy[2] = -0.51*y[0]*y[1]
> return dy
>
> r = ode(foo).set_integrator('vode').set_initial_value( y0,t0)
> while r.successful() and r.t < tEnd:
> r.integrate(r.t+dt)
> Y[index] = r.y
> T[index] = r.t
> index += 1
>
> As a further aside, this system of three coupled linear ODE's is from
> an example in MATLAB's documentation on their function ode45, so I
> know what I 'should' get. The while loop and call to ode is adapted
> from scipy's (very limited) documentation on scipy.integrate.ode.ode.
>
> At the end of this while loop I've gotten a few different results
> depending on what I have no idea. This morning another TypeError
> exception was thrown, saying:
>
> TypeError: Array can not be safely cast to required type.
>
> Although, to be fair this is after the output from one iteration:
>
> array([0.01, 1. , 1. ])
>
> So, clearly this isn't working right. Does anyone have any experience
> using this function or anything they can contribute?


I tried your example and it worked for me (after fixing the IndexError). Could
you come to scipy-user with the complete code that you ran (the one above isn't
the one you ran; there is at least one typo), and the full traceback of the error?

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco

 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
python,win32com,scipy and vb 6 : no module named scipy vml Python 3 05-03-2007 07:57 AM
OT: Ode To Colorectal Surgeon [updated] Verbal Funderburk MCSE 6 03-06-2006 04:18 PM
An ode to re.finditer() Robert Oschler Python 0 07-31-2004 11:57 PM
Scipy question about scipy.optimize.leastsq Paxcal Python 0 02-11-2004 07:52 PM
Ode to a Slit Steve Young Digital Photography 7 10-22-2003 03:05 PM



Advertisments