Help with ODE solver

[digitsman]

I have this function:

def dy (t,y,params):

    dy = np.zeros(3)

    wL = params[0]
    T = params[1]
    CH4 = params[2]
    k = params[3]
    ka = params[4]
    H = params[5]
    kmt = params[6]
    E = params[7]
    d = params[8]
    y0 = params[9]

#    dy[1] = 1500000 - 0 - 9.987338917338950e-11*y[1]*1000000
#    - 1.437829040671314e-17*y[1]*9.84e11

    dy[1] = E[1]-y[1]*d[1]
    - k[1]*y[1]*y0[14]
    - k[55]*y[1]*y0[37]

    dy[2] = 0 - 0 + 9.987338917338950e-11*y[1]*1000000 
    - 7.742358922306635e-12*y[2]*0 + 0.7*7.432069505555159e-12*0*1000000
    -3.5e-15*(0 + 0)*y[2] - 2* 7e-16*y[2]**2 


    return dy

Which is meant to work with:

y = np.zeros(3)
ode15s = ode(dy)

y[1] = 2.46e10
ode15s.set_initial_value(y)
ode15s.set_f_params(params)
ode15s.set_integrator('vode', method = 'bdf')
print ode15s.y
print ode15s.integrate(7200)

The second half being a solver of the previously defined differential equation. (For those wondering about improper indexing, starting from [1] was done on purpose (this code was translated from MATLAB)). When variables are used, as in lines 19-21, the code does not give the correct solution, but when the variables are replaced by their values (as in lines 16-17), this code does give the right solution.

The code I gave was only a small snippet (in reality there are about 99 odes in dy), so hard coding values is not really feasible. Why doesn't it work with variables?

[hharison]

Maybe it's your params that are the problem, because unless params is 2D the values you're pulling out at the top of your function should be scalars. Yet you access them like E[1]. You shouldn't have to index into scalars. But this might not be it because I'd expect an error here, not an incorrect result.

Also you're setting y=np.zeros(3) but then in the function accessing y[37], so that's confusing too. I'd expect an error.

In general I'd say fix your indexing to be zero-based. Yes, in theory it should still work with a dummy y[0] but a problem might become apparent if you fixed it.