Skip to main content
 首页 » 编程设计

python之将 numpy 数组与 scipy odeint 一起使用

2025年05月04日14itcoder

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

请您参考如下方法:

您可以通过为 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 代码,翻译是完全自动化的。