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