3

基本上......我需要一种在我的微分方程中包含相移的方法。也就是说,我的系统函数定义中有返回 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[?] < ttheta[?] > 0。这不是很好。

我看到在这种情况下,我可以单独找到 eta 的每个组件,然后获取先前计算的 eta 组件(eta_0、eta_1、eta_2)的计算值来计算 eta_3,并且类似地计算 eta_4,但这并不理想,因为它剥夺了我可以“即插即用”任何通用公式。

4

3 回答 3

3

有许多现有的库和示例可以做到这一点。

http://www.google.fi/search?q=python+delay+differential+equation给了我:

于 2013-04-04T15:53:34.267 回答
2

延迟不完全是线性函数。通常的阶跃延迟在拉普拉斯域中表示为e**(a*s)/s,其中a是延迟。

这意味着除非您有一些解决方法,否则“正常”ODE 求解器将无法工作。通常这种解决方法不是很容易做到,因为对于僵硬的问题,您通常无法使用足够好的近似值进行插值。

无论如何,其中一个解决方案是使用其他答案中发布的库。

另一种解决方案是象征性地进行(如果可以,您可以尝试SymPy)。

第三种解决方案是存储过去的结果并进行插值以找到您需要的确切过去(可能不够好)。

第四种解决方案可能是某些模拟器文档所推荐的:c2d()在离散时间使用和模拟整个模型,并将过去的变量存储在列表/数组中(没有插值,但您可能需要使用小步骤以获得更好的准确性)。

第五个解决方案是使用Padé 近似来表示模型的延迟(可能会根据您的情况起作用)。python-control 中有一个pade()函数可以精确地近似这个。

于 2015-07-18T02:07:22.477 回答
1

一种方法是在开始时间和结束时间之间integrate.odeint()运行integrate.odeint()许多短时间间隔,在列表中的每个短时间间隔之后存储时间值和输出 Y 值。例如,每次需要时,您都可以使用scipy.interpolate.interp1d()在列表中插入 Y 值Y(t-3)

当然,如果你这样做,你只会得到一个近似值Y(t-3),但如果时间间隔足够接近,这种方法可能会让你满意。毕竟,Y(t)数值 ODE 求解器计算的值也是近似值。

于 2013-04-04T05:59:14.180 回答