Accessing earlier values in odeint

Question

I'm having some trouble using odeint function from scipy. I'm translating a discrete system into a continuous one, but some equation in the discrete model requires that I access to the previous value of a variable that I'm currently integrating. How could I translate this behaviour?

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]

I've posted a snippet of the code, the problem is with the denominator of the first two equations. Thanks in advance.

Answer 1

From the interaction terms you can directly read off the ODE system \begin{align} \dot S &= \nu\frac{S·I}{N-H-Q-D} \\ \dot I &=-\nu·\frac{S·I}{N-H-Q-D} - \gamma·I - \alpha·I - \psi·I \\ \dot Q &= \alpha·I - \eta_q·Q - \kappa_h·Q + \kappa_q·H \\ \dot H &= \psi·I - \eta_h·H + \kappa_h·Q - \kappa_q·H - \zeta·H \\ \dot R &= \eta_q·Q + \eta_h·H \end{align} If you now solve that with the Euler method with step size $1$, you get the discrete method back. However, if you simualate that with any other method, that is, get closer to the exact ODE solution, the result might be distinctly different. Compare Question about Euler's Method and the SIR epidemic model using a spreadsheet for a similar question (only backwards, from ODE to large-step Euler) on the basic Zombie Apocalypse model.