使用带有 scipy odeint 的 numpy 数组答案

【问题标题】：Using numpy arrays with scipy odeint使用带有 scipy odeint 的 numpy 数组
【发布时间】：2013-10-01 12:14:17
【问题描述】：

我正在使用 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”的值。因此，例如，我不想只有一个 = -2.0，而是：

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

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

【问题讨论】：

标签： python arrays scipy odeint

【解决方案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 代码，翻译完全自动化。

【讨论】：