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

Posted 2023-02-18

【中文标题】在 Python 中求解递归微分方程组

【英文标题】：Solve system of recursive differential equation in Python

【发布时间】：2021-08-02 19:30:51

【问题描述】：

所以我试图在 Python 中求解以下微分方程组。

System of differential equations

如您所见，对于 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, z0, 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, h2h3, 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, 0, 0,    0], 
              [0, 0, 0, 0, 0,    0],
              [0, 0, 0, 0, h3a0, 0],
              [0, 0, 0, 0, 0,    0],
              [0, 0, 0, 0, 0,    0],
              [0, 0, 0, 0, 0,    0],
              ])


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,1,1,1,1,1,0,0,0,0,0,0]
t = np.linspace(0,90,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,1,1,1,1]
t = np.linspace(0,90)
z = odeint(modeln0manually, u0, t)

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

【参考方案1】：

您的系统是耦合的，即使是三角形的。所以最紧凑的方法就是把它作为一个耦合系统来解决

A = np.array([[-3.273, 3.2], [3.041, -3.041]])
B = np.array([[0, 0.071], [0, 0]])

def model0n(u,t):
    Ra = u.reshape([-1,2])
    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,1,0,0]
t = np.linspace(0,90,90+1)
u = odeint(model0n,u0,t)

【讨论】：

感谢您的回答，但是我的编码技能在这里可能有点不足。在使用 odeint 语句时尝试您的解决方案会出错，因为 Python 无法将 u 重塑为具有 n 行和 2 列的向量，n=1 除外，这导致它仅返回问题的初始条件。不过，我可能从您的解决方案中误解了某些内容。

u 需要包含R(0|t) 到R(n|t) 的2*(n+1) 条目，它同时评估所有组件DE 的右侧。我使它更灵活，没有外部尺寸n。请注意，这还没有测试，确实循环索引是错误的。

现在它已经过测试，它产生了结果。现在增加n 意味着向u0 添加更多条目，始终是偶数。

您提供的解决方案很有魅力，谢谢！但是，当将模型扩展为具有六个方程的系统时，使 A 和 B 成为 6x6 矩阵，我又遇到了一些麻烦。我用 [-1,6] 重塑 Ra 并提供了新矩阵，但它不会产生 n > 0 的值和 n = 0 的错误概率。如果我在我的问题中首先返回并使用 a_ij*z[j] 放置六个方程，它适用于 n=0。扩展模型时我缺少什么吗？这对我来说似乎很直接。我已将我的新代码放入 A 和 B 中具有任意值的问题中

会不会是B 中的尾随逗号？我在代码中没有发现其他异常情况。也许我错过了.flatten() 的一些内容，也可以使用reshape([-1]) 来解构状态数据。 // 确切的错误信息是什么？也可以再次检查np.matmul（@ 使用的操作），但这应该只是实现通常的矩阵乘法。