So I've written a function to use a particular solver type to integrate a couple system of ODEs. For some reason, which I can't quite get a handle on, the function I wrote to define the RHS of the ODE is returning a tuple when run inside the integrating function (or so the error indicates), but whenever I test the array it returns separately, the array seems to be...well... and array, which is supposed to work fine.
def solve_decay_system(t, u_0, solver_type="dopri5", K_1=0.0, K_2=0.0, K_3=0.0):
def decay(u, K_1 = K_1, K_2=K_2, K_3=K_3):
dydt = numpy.array([-K_1 * u[0], K_1 * u[0] - K_2 * u[1], K_2 * u[1] - K_3*u[2], K_3*u[2]])
dydt = numpy.array(dydt)
print dydt, type(dydt)
return dydt
#return [1, 2, 3, 4]
u = numpy.empty((4, t.size))
u_0 = numpy.array([1.0, 0.0, 0.0, 0.0])
u[0, 0] = u_0[0]
u[1, 0] = u_0[1]
u[2, 0] = u_0[2]
u[3, 0] = u_0[3]
integrator = integrate.ode(decay)
integrator.set_integrator(solver_type)
integrator.set_initial_value(u[:, 0])
integrator.set_f_params(u, K_1, K_2, K_3)
for (n, t_n) in enumerate(t[1:]):
print n
integrator.integrate(t_n)
if not integrator.successful():
break
u[:, n + 1] = integrator.y
return u
Any suggestions as to how to resolve this would be very helpful. Here is the resultant error, for reference:
ValueError Traceback (most recent call last)
<ipython-input-230-8f9be46c3d95> in <module>()
----> 1 solve_decay_system(t, u_0, solver_type="dopri5", K_1=0.0, K_2=0.0, K_3=0.0)
<ipython-input-229-a05ddca0f334> in solve_decay_system(t, u_0, solver_type, K_1, K_2, K_3)
24 for (n, t_n) in enumerate(t[1:]):
25 print n
---> 26 integrator.integrate(t_n)
27 if not integrator.successful():
28 break
C:\Users\Ben\Anaconda2\lib\site-packages\scipy\integrate\_ode.pyc in integrate(self, t, step, relax)
409 except SystemError:
410 # f2py issue with tuple returns, see ticket 1187.
--> 411 raise ValueError('Function to integrate must not return a tuple.')
412
413 return self._y
ValueError: Function to integrate must not return a tuple.
The derivative function expects (time, state,...) as first arguments,
def decay(t, u, K_1 = K_1, K_2=K_2, K_3=K_3):
even if the ODE is autonomous.
The state is not a parameter, esp. if you give the list of states as parameter.
integrator.set_f_params(K_1, K_2, K_3)
With these changes your code should run. (Up to n=8.)
Even if you just return a list
def decay(t, u, K_1 = K_1, K_2=K_2, K_3=K_3):
return [-K_1 * u[0], K_1 * u[0] - K_2 * u[1], K_2 * u[1] - K_3*u[2], K_3*u[2]]