如何在时间= 50 ( S(50)，I(50)，R(50) )中找出易感、感染和康复个体的数量？(SIR模型)

Question

如何在时间= 50 ( S(50)，I(50)，R(50) )中找出易感、感染和康复个体的数量？(SIR模型)

# Equações diferenciais e suas condições iniciais
h = 0.05
beta = 0.8
nu = 0.3125

def derivada_S(time,I,S):
    return -beta*I*S

def derivada_I(time,I,S):
    return beta*I*S - nu*I

def derivada_R(time,I):
    return nu*I

S0 = 0.99
I0 = 0.01
R0 = 0.0

time_0 = 0.0
time_k = 100
data = 1000

# vetor representativo do tempo
time = np.linspace(time_0,time_k,data)

S = np.zeros(data)
I = np.zeros(data)
R = np.zeros(data)

S[0] = S0
I[0] = I0
R[0] = R0

for i in range(data-1):
    S_k1 = derivada_S(time[i], I[i], S[i])
    S_k2 = derivada_S(time[i] + (1/2)*h, I[i], S[i] + h + (1/2)*S_k1)
    S_k3 = derivada_S(time[i] + (1/2)*h, I[i], S[i] + h + (1/2)*S_k2)
    S_k4 = derivada_S(time[i] + h,  I[i], S[i] + h + S_k3)
    
    S[i+1] = S[i] + (h/6)*(S_k1 + 2*S_k2 + 2*S_k3 + S_k4)

    I_k1 = derivada_I(time[i], I[i], S[i])
    I_k2 = derivada_I(time[i] + (1/2)*h, I[i], S[i] + h + (1/2)*I_k1)
    I_k3 = derivada_I(time[i] + (1/2)*h, I[i], S[i] + h + (1/2)*I_k2)
    I_k4 = derivada_I(time[i] + h,  I[i], S[i] + h + I_k3)
    
    I[i+1] = I[i] + (h/6)*(I_k1 + 2*I_k2 + 2*I_k3 + I_k4)
    
    R_k1 = derivada_R(time[i], I[i])
    R_k2 = derivada_R(time[i] + (1/2)*h, I[i])
    R_k3 = derivada_R(time[i] + (1/2)*h, I[i])
    R_k4 = derivada_R(time[i] + h, I[i])
    
    R[i+1] = R[i] + (h/6)*(R_k1 + 2*R_k2 + 2*R_k3 + R_k4)

plt.figure(figsize=(8,6))
plt.plot(time,S, label = 'S')
plt.plot(time,I, label = 'I')
plt.plot(time,R, label = 'R')
plt.xlabel('tempo (t)')
plt.ylabel('Susceptível, Infectado e Recuperado')
plt.grid()
plt.legend()
plt.show()

我正在解决一个大学问题，python应用Runge的第四阶，但我不知道如何收集时间= 50的数据。

Answer 1

这个链接可能会帮助您构建模型SIR衍生ODE模型，这里我也为您提供了代码：

import numpy as np
import matplotlib.pyplot as plt

Beta = 1.00205
Gamma = 0.23000
N = 1000

def func_S(t,I,S):
    return - Beta*I*S/N

def func_I(t,I,S):
    return Beta*I*S/N - Gamma*I

def func_R(t,I):
    return Gamma*I


# physical parameters
I0 = 1
R0 = 0
S0 = N - I0 - R0
t0 = 0
tn = 50



# Numerical Parameters
ndata = 1000



t = np.linspace(t0,tn,ndata)
h = t[2] - t[1]

S = np.zeros(ndata)
I = np.zeros(ndata)
R = np.zeros(ndata)

S[0] = S0
I[0] = I0
R[0] = R0


for i in range(ndata-1):
    k1 = func_S(t[i], I[i], S[i])
    k2 = func_S(t[i]+0.5*h, I[i], S[i]+h+0.5*k1)
    k3 = func_S(t[i]+0.5*h, I[i], S[i]+h+0.5*k2)
    k4 = func_S(t[i]+h, I[i], S[i]+h+k3)
    
    S[i+1] = S[i] + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
    
    kk1 = func_I(t[i], I[i], S[i])
    kk2 = func_I(t[i]+0.5*h, I[i], S[i]+h+0.5*kk1)
    kk3 = func_I(t[i]+0.5*h, I[i], S[i]+h+0.5*kk2)
    kk4 = func_I(t[i]+h, I[i], S[i]+h+kk3)
    
    I[i+1] = I[i] + (h/6)*(kk1 + 2*kk2 + 2*kk3 + kk4)
    
    l1 = func_R(t[i], I[i])
    l2 = func_R(t[i]+0.5*h, I[i])
    l3 = func_R(t[i]+0.5*h, I[i])
    l4 = func_R(t[i]+h, I[i])
    
    R[i+1] = R[i] + (h/6)*(l1 + 2*l2 + 2*l3 + l4)
    
    
plt.figure(1)
plt.plot(t,S)
plt.plot(t,I)
plt.plot(t,R)
plt.show()

输出如下：

​

Answer 2

您可以使用scipy.integrate.solve_ivp上可用的积分器，并使用四阶Runge方法(DOP853、RK23、RK45和Radau)。

##########################################
# AUTHOR  : CARLOS DUARDO DA SILVA LIMA  #
# DATE    : 12/01/2022                   #
# LANGUAGE: python                       #
# IDE     : GOOGLE COLAB                 #
# PROBLEM : MODEL SIR                    #
##########################################

import numpy as np
from scipy.integrate import odeint, solve_ivp, RK45
import matplotlib.pyplot as plt

t_i = 0.0 # START TIME
t_f = 50.0 # FINAL TIME
N   = 1000

#t = np.linspace(t_i,t_f,N)
t_span = np.array([t_i,t_f])

# INITIAL CONDITIONS OF THE SOR MODEL
S0 = 0.99
I0 = 0.01
R0 = 0.0
r0 = np.array([S0,I0,R0])

# ORDINARY DIFFERENTIAL EQUATIONS OF THE SIR MODEL
def SIR(t,y,b,k):
  s,i,r = y
  ode1 = -b*s*i
  ode2 = b*s*i-k*i
  ode3 = k*i
  return np.array([ode1,ode2,ode3])

# INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS (FOURTH ORDER RUNGE-KUTTA, RADAU)
#sol_solve_ivp = solve_ivp(SIR,t_span,y0 = r0,method='Radau', rtol=1E-09, atol=1e-09, args = (0.8,0.3125))
sol_solve_ivp = solve_ivp(SIR,t_span,y0 = r0,method='RK45', rtol=1E-09, atol=1e-09, args = (0.8,0.3125))

# T, S, I, R FUNCTIONS
t_= sol_solve_ivp.t
s = sol_solve_ivp.y[0, :]
i = sol_solve_ivp.y[1, :]
r = sol_solve_ivp.y[2, :]

# GRAPHIC
plt.figure(1)
plt.style.use('dark_background')
plt.figure(figsize = (8,8))
plt.plot(t_,s,'c-',t_,i,'g-',t_,r,'y-',lw=1.5)
#plt.title(r'$\frac{dS(t)}{dt} = -bs(t)i(t)$, $\frac{dI(t)}{dt} = bs(t)i(t)-ki(t)$ and $\frac{dR(t)}{dt} = ki(t)$')
plt.title(r'SIR Model', color = 'm')
plt.xlabel(r'$t(t)$', color = 'm')
plt.ylabel(r'$S(t)$, $I(t)$ and $R(t)$', color = 'm')
plt.legend(['S', 'I', 'R'], shadow=True)
plt.grid(lw = 0.95,color = 'white',linestyle = '--')
plt.show()

''' SEARCH WEBSITES
https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology
https://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-the-differential-equation-model
'''

输出图