在 Python 中求解递归微分方程组

问题描述

所以我试图用 Python 解决以下微分方程组。

如您所见，对于 {0,1,2,3,...} 中的每个 n，系统依赖于前一个系统的解决方案。

我尝试解决 n=0 的系统，并找到了一个解决方案 R(0|t)，我可以将其插入 R(1|t) 中，Python 可以毫无问题地解决系统。我已将解 R(0|t) 定义为 r0(t) 并按如下方式实现了 n=1 的解：

def model(z,t):
    dxdt = -3.273*z[0] + 3.2*z[1] + r0(t)
    dydt = 3.041*z[0] - 3.041*z[1]
    dzdt = [dxdt,dydt]
    return dzdt

z0 = [0,0]

t = np.linspace(0,90,90)

z1 = odeint(model,z0,t)

但是，我想通过在求解 n 时调用 n-1 系统的解决方案来概括这个解决方案。由于微分方程只有矩阵右上角不为零的项，我们只需要从先前的解中得到 z1 的解。我试过了

def model0(z,t):
    dxdt = -3.273*z[0] + 3.2*z[1] 
    dydt = 3.041*z[0] - 3.041*z[1]
    dzdt = [dxdt,dydt]
    return dzdt

z0 = [1,1]

t = np.linspace(0,90)

def model1(z,t):
    dxdt = -3.273*z[0] + 3.2*z[1] + 0.071*odeint(model0,t)[t,1]
    dydt = 3.041*z[0] - 3.041*z[1]
    dzdt = [dxdt,dydt]
    return dzdt


z1 = [0,0]


z = odeint(model1,z1,t)

没有任何运气。有没有人在 Python 中解决这些递归的 odes 系统的经验？

提前致谢。

更新了 6x6 矩阵和 6 个函数的代码：


A = np.array([[h1h1,h1h2,h1h3,h1a1,h1a2,h1a3],[h2h1,h2h2,h2h3,h2a1,h2a2,h2a3],[h3h1,h3h3,h3a1,h3a2,h3a3],[a1h1,a1h2,a1h3,a1a1,a1a2,a1a3],[a2h1,a2h2,a2h3,a2a1,a2a2,a2a3],[a3h1,a3h2,a3h3,a3a1,a3a2,a3a3]
              ])


B = np.array([[0,0],[0,h3a0,])


def model0n(u,t):
    Ra = u.reshape([-1,6])
    n = len(Ra) - 1
    dRa = np.zeros(Ra.shape)
    dRa[0] = A @ Ra[0]
    for i in range(1,n+1): 
        dRa[i] = A @ Ra[i] + B @ Ra[i-1]
    return dRa.flatten()

u0 = [1,0]
t = np.linspace(0,90+1)

u = odeint(model0n,u0,t)

上面的结果是 u[:,0] 的下图： Plot for u[:,0] which is supposed to be probabilities

对于 n=0，它提供执行矩阵乘积“手动”的结果：


def modeln0manually(z,t):
    d1dt = h1h1*z[0] + h1h2 * z[1] + h1h3*z[2] + h1a1*z[3] + h1a2*z[4] + h1a3*z[5]
    d2dt = h2h1*z[0] + h2h2 * z[1] + h2h3*z[2] + h2a1*z[3] + h2a2*z[4] + h2a3*z[5]
    d3dt = h3h1*z[0] + h3h2 * z[1] + h3h3*z[2] + h3a1*z[3] + h3a2*z[4] + h3a3*z[5]
    d4dt = a1h1*z[0] + a1h2 * z[1] + a1h3*z[2] + a1a1*z[3] + a1a2*z[4] + a1a3*z[5]
    d5dt = a2h1*z[0] + a2h2 * z[1] + a2h3*z[2] + a2a1*z[3] + a2a2*z[4] + a2a3*z[5]
    d6dt = a3h1*z[0] + a3h2 * z[1] + a3h3*z[2] + a3a1*z[3] + a3a2*z[4] + a3a3*z[5]
    drdt = [d1dt,d2dt,d3dt,d4dt,d5dt,d6dt]    
    return drdt    


u0 = [1,1]
t = np.linspace(0,90)
z = odeint(modeln0manually,t)

导致 u[:,0] 的情节： Plot of u[:,0] as it is supposed to be

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

differential-equations odeint python recursion