need help with ODE solver from scipy

Discussion in 'Python' started by T.Crane, Apr 24, 2007.

  1. T.Crane

    T.Crane Guest

    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
    T.Crane, Apr 24, 2007
    #1
    1. Advertising

  2. T.Crane

    Robert Kern Guest

    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
    Robert Kern, Apr 24, 2007
    #2
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Paxcal
    Replies:
    0
    Views:
    1,828
    Paxcal
    Feb 11, 2004
  2. Robert Oschler

    An ode to re.finditer()

    Robert Oschler, Aug 1, 2004, in forum: Python
    Replies:
    0
    Views:
    335
    Robert Oschler
    Aug 1, 2004
  3. Ron Alford

    Ode to Urllib2

    Ron Alford, Oct 14, 2004, in forum: Python
    Replies:
    1
    Views:
    299
    Avi Berkovich
    Oct 15, 2004
  4. Jamey Cribbs

    An Ode To My Two Loves

    Jamey Cribbs, Feb 5, 2005, in forum: Python
    Replies:
    7
    Views:
    324
    Mike Meyer
    Feb 9, 2005
  5. vml
    Replies:
    3
    Views:
    686
Loading...

Share This Page