need help with ODE solver from scipy

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

  1. T.Crane

    T.Crane Guest


    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:
    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?

    T.Crane, Apr 24, 2007
    1. Advertisements

  2. T.Crane

    Robert Kern Guest

    You will get more help from the scipy-user list than you will here.
    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.
    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.
    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
    1. Advertisements

Ask a Question

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

You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.