python - 将 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)

但现在我想针对常量“a”的多个值求解此系统。因此，例如，我不希望只有 a = -2.0，而是:

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

并针对 a 的每个值求解系统。有没有一种方法可以做到这一点而不必遍历数组的每个元素？我可以一次完成吗？

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 代码，翻译是完全自动化的。