need help with ODE solver from scipy

T.Crane trevis.crane at gmail.com
Tue Apr 24 11:53:15 EDT 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




More information about the Python-list mailing list