Using numpy arrays with scipy odeint

Question

I'm using scipy to solve a system of ordinary differential equations. For simplicity, take my code to be:

import scipy as sp
import numpy as np
from scipy.integrate import odeint
from numpy import array

def deriv(y,t): # return derivatives of the array y
 a = -2.0
 b = -0.1
 return array([ y[1], a*y[0]+b*y[1] ])
time = np.linspace(0.0,10.0,100)
yinit = array([0.0005,0.2]) # initial values
y = odeint(deriv,yinit,time)

But now I want to solve this system for several values of the constant "a". So instead of having just a = -2.0 for example, I'd like to have:

a = array([[-2.0,-1.5,-1.,-0.5]])

and solve the system for each value of a. Is there a way to do this without having to loop over each element of the array? Can I do it all at once?

Answer 1

You can reformulate your system of equations by writing two equations for each value of a. One way to do it is

from scipy.integrate import odeint
from numpy import array,linspace,tile,empty_like
a = array([-2.0,-1.5,-1.,-0.5])
b = array([-.1,-.1,-.1,-.1])

yinit = tile(array([0.0005,0.2]),a.size)

def deriv(y,_):
    dy = empty_like(y)
    dy[0::2]=y[1::2]
    dy[1::2]=y[::2]*a+b*y[1::2]
    return dy
time = linspace(0.0,10.0,100)
y = odeint(deriv,yinit,time)

you will find in y[:,0] the solution for y for a=-2, in y[:,2] the solution for a=-1.5, and so on and so forth with y[:,-1] being the derivative of y for a=-.5.

Anyway, if you want to vectorize this to gain speed, you may be interested in a library that converts your python code to c, and then compiles it. I personally use pydelay because I need to simulate time delays as well, and would recommend it. You don't even have to deal with the c code, the translation is completely automated.