基本上......我需要一种在我的微分方程中包含相移的方法。也就是说,我的系统函数定义中有返回 dY/dt 之类的 Y(t-3)。就像这个微分方程:
dY/dt = a*Y(t) + b*Y(t-tau)
现在,如果我尝试将其编写为传递给 scipy.odeint 的系统定义函数,我会迷失方向:
def eqtnSystem(A,t):
Y = A
a = 1
b = 5
tau = 3
return a*Y + b*??? # how do I Y(t-tau) ?
基本上就是这样。我真的希望有一个简单的答案,但我似乎无法找到一个。
具体来说......我正在尝试以数值方式计算由以下函数定义的系统的解决方案:
def etaFunc(A,t):
#...definition of all those constants is here...
return array([(gamma[0,0]*xi(t-theta[0])[0] - eta[0] + zeta[0])/tau[0],\
(gamma[1,1]*xi(t-theta[1])[1] - eta[1] + zeta[1])/tau[1],\
(gamma[2,2]*xi(t-theta[2])[2] - eta[2] + zeta[2])/tau[2],\
( beta[3,0]*pastEta(t-theta[3])[0] \
+ beta[3,1]*pastEta(t-theta[4])[1] \
+ beta[3,2]*pastEta(t-theta[5])[2] -eta[3]+ zeta[3])/tau[3],\
( beta[4,3]*pastEta(t-theta[6])[3] \
+ beta[4,2]*pastEta(t-theta[7])[2] - eta[4] + zeta[4])/tau[4]])
然后这个函数随后被赋予 odeint,如下所示:
ETA = integrate.odeint(etaFunc,initCond,time)
然后我可以像这样取出 ETA 的每个单独组件(如 eta_0)ETA[:,0]
:
我在这里遇到的问题是pastEta(t-theta[?])
. 目前,这是一个试图找到已经计算的 eta 值的函数(对于 whenstart_time < t-theta[?] < t
和theta[?] > 0
。这不是很好。
我看到在这种情况下,我可以单独找到 eta 的每个组件,然后获取先前计算的 eta 组件(eta_0、eta_1、eta_2)的计算值来计算 eta_3,并且类似地计算 eta_4,但这并不理想,因为它剥夺了我可以“即插即用”任何通用公式。