首页 文章

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

提问于
浏览 次
3

Basically ...我需要一种在微分方程中包含相移的方法 . 也就是说，我在我的系统函数的定义中返回了类似Y（t-3）的dY / dt . 像这个微分方程：

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) ?

基本上就是这样 . 我真的希望有一个简单的答案，但我似乎无法跟踪一个 .

Specifically ...我试图用数字计算由以下函数定义的系统的解决方案：

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，但这不是理想的，因为它消除了我可以“即插即用”任何通用公式 .

python scipy ode odeint

3 回答

  • 3
    回复于
  • 1

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

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

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

    另一种解决方案是以符号方式进行（如果可以，您可以尝试SymPy） .

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

    第四种解决方案可能是某些模拟器docs推荐的：使用 c2d() 并在离散时间内模拟整个模型并将过去的变量存储在列表/数组中（无插值，但您可能需要使用小步骤以获得更高的精度） .

    第五个解决方案是使用Padé approximation来表示你的模型's delay (might working depending on your case). There' s a pade() function在python-control中恰好接近这个 .

    回复于
  • 2

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

    如果你这样做的话，你最终只得到 Y(t-3) 的近似值，但是如果时间间隔足够接近，这种方法对你来说可能是令人满意的 . 毕竟，由数值ODE求解器计算的 Y(t) 值也是近似值 .

    回复于

相关问题