Использование numpy массивов с scipy odeint

Question

Я использую scipy для решения системы обыкновенных дифференциальных уравнений. Для простоты, возьмите мой код, чтобы быть:

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) 

Но теперь я хочу, чтобы решить эту систему для нескольких значений константы «а». Так что вместо того, чтобы просто = -2,0, например, я хотел бы иметь:

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

и решить систему для каждого значения а. Есть ли способ сделать это без необходимости перебирать каждый элемент массива? Могу ли я сделать все сразу?

Answer 1

Вы можете переформулировать свою систему уравнений, написав два уравнения для каждого значения a. Один из способов сделать это является

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) 

вы найдете y[:,0] решение для y для a=-2, в y[:,2] решение для a=-1.5, и так далее, и так далее с y[:,-1] является производной y для a=-.5.

В любом случае, если вы хотите векторизовать это, чтобы получить скорость, вам может быть интересна библиотека, которая преобразует ваш код на python в c и затем скомпилирует его. Я лично использую pydelay, потому что мне тоже нужно моделировать временные задержки, и я бы порекомендовал его. Вам даже не нужно иметь дело с кодом c, перевод полностью автоматизирован.