访问 odeint 中的早期值

Question

我在使用 scipy 中的 odeint 函数时遇到了一些问题。我正在将离散系统转换为连续系统，但离散模型中的某些方程要求我访问我当前正在积分的变量的先前值。我该如何翻译这种行为？

import numpy as np
days_of_prediction = 15
N = 100
discrete_S0 = np.zeros((days_of_prediction, 1))
discrete_I0 = np.zeros((days_of_prediction, 1))
discrete_Q0 = np.zeros((days_of_prediction, 1))
discrete_H0 = np.zeros((days_of_prediction, 1))
discrete_D0 = np.zeros((days_of_prediction, 1))
discrete_S0[0] = 99
discrete_I0[0] = 1
discrete_Q0[0] = 0
discrete_H0[0] = 0
discrete_D0[0] = 0
v=0.1
alpha = 0.3
gamma = 1/21
psi = 0.2
k_h=0.1
k_q=0.1
eta_h=0.3
eta_q=0.3
for t in range(days_of_prediction - 1):
    discrete_S0[t + 1] = discrete_S0[t] - v * discrete_S0[t] * 
discrete_I0[t] / (N - discrete_Q0[t] - discrete_H0[t] - discrete_D0[t])
    discrete_I0[t + 1] = discrete_I0[t] + v * discrete_S0[t] * discrete_I0[t] / (N - discrete_Q0[t] - discrete_H0[t] - discrete_D0[t]) - gamma * discrete_I0[t] - alpha * discrete_I0[t] - psi * discrete_I0[t]
    discrete_Q0[t + 1] = discrete_Q0[t] + alpha * discrete_I0[t] - eta_q * discrete_Q0[t] - k_h *discrete_Q0[t] + k_q * discrete_H0[t]
    discrete_H0[t + 1] = discrete_H0[t] + psi * discrete_I0[t] - eta_h * discrete_H0[t] + k_h * discrete_Q0[t] - k_q * discrete_H0[t] - zeta * discrete_H0[t]
    discrete_R0[t + 1] = discrete_R0[t] + eta_q * discrete_Q0[t] + eta_h * discrete_H0[t]

我已经发布了代码的 sn-p，问题在于前两个方程的分母。 提前致谢。

Answer 1

在这样的方程系统中，一些变量的先前值需要在其他变量的演化方程中，你可以定义你的函数如下：

def fun(RHS, t):

    # get initial boundary condition values
    discrete_S0 = RHS[0]
    discrete_I0 = RHS[1]
    discrete_Q0 = RHS[2]
    discrete_H0 = RHS[3]
    discrete_D0 = RHS[4]

    # calculte rate of respective variables
    discrete_S0dt = - v * discrete_S0 * discrete_I0 / (N - discrete_Q0 - discrete_H0 - discrete_D0)
    discrete_I0dt = v * discrete_S0 * discrete_I0 / (N - discrete_Q0 - discrete_H0 - discrete_D0) - gamma * discrete_I0 - alpha * discrete_I0 - psi * discrete_I0
    discrete_Q0dt = alpha * discrete_I0 - eta_q * discrete_Q0 - k_h *discrete_Q0 + k_q * discrete_H0
    discrete_H0dt = psi * discrete_I0 - eta_h * discrete_H0 + k_h * discrete_Q0 - k_q * discrete_H0 - zeta * discrete_H0
    discrete_D0dt = eta_q * discrete_Q0 + eta_h * discrete_H0

    # Left-hand side of ODE
    LHS = np.zeros([5,])

    LHS[0] = discrete_S0dt
    LHS[1] = discrete_I0dt
    LHS[2] = discrete_Q0dt
    LHS[3] = discrete_H0dt
    LHS[4] = discrete_D0dt

    return LHS

注意：检查所有变量的速率方程（discrete_S0dt、discrete_I0dt 等），我可能写错了。请验证并更正自己。

之后，你可以解决它（根据你的边界条件）如下：

import numpy as np
from scipy.integrate import odeint

y0 = [99, 1, 0, 0, 0]
t = np.linspace(0,13,14)

res = odeint(fun, y0, t)

这里 y0 是函数 fun 在 t=0 时定义的所有变量的初始边界条件。这就是变量 t 从 0 开始的原因。

另外，你可以得到所有变量的结果如下：

res[:,0]
res[:,1]
res[:,2]
res[:,3]
res[:,4]