在具有相移变量的方程上使用 scipy odeint

Question

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

我看到在这种情况下，我可以单独找到 eta 的每个组件，然后获取先前计算的 eta 组件（eta_0、eta_1、eta_2）的计算值来计算 eta_3 和类似的 eta_4，但这并不理想，因为它带走了我能够“即插即用”任何通用公式。

Answer 1

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

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

• http://zulko.wordpress.com/2013/03/01/delay-differential-equations-easy-with-python/（使用scipy.integrate.ode）
• https://pypi.python.org/pypi/PyDDE
• http://pydelay.sourceforge.net/

Answer 2

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

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

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

另一种解决方案是象征性地这样做（如果可以，您可以尝试SymPy）。

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

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

第五个解决方案是使用Padé approximation 来表示您的模型的延迟（根据您的情况可能会起作用）。在 python-control 中有一个 pade() function 来近似这个。

Answer 3

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

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