2

我正在使用 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 的每个值求解系统。有没有办法做到这一点而不必遍历数组的每个元素？我可以一次性完成吗？

4

1 回答 1

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]for 的解决方案中y，a=-2在fory[:,2]的解决方案中a=-1.5，等等，都是fory[:,-1]的导数。ya=-.5

无论如何，如果您想对其进行矢量化以提高速度，您可能会对将您的 python 代码转换为 c 的库感兴趣，然后对其进行编译。我个人使用pydelay是因为我也需要模拟时间延迟，并且会推荐它。您甚至不必处理 c 代码，翻译是完全自动化的。

于 2013-09-25T21:55:32.973 回答