Как узнать количество восприимчивых, зараженных и выздоровевших за время = 50, где S(50), I(50), R(50)? (СЭР МОДЕЛЬ)

Эта ссылка может помочь вам построить модель модели ODE, производные от SIR также здесь у меня есть код для вас:

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()

Вывод будет таким:

Самый простой способ получить значение в момент времени 50 — это вычислить значение в момент времени 50. Когда вы вычисляете данные за 100 дней с примерно 10 точками данных в день, отразите это в построении временного массива.

time = np.linspace(0,days,10*days+1)

Note that linspace(a,b,N) produces N nodes that have between them a step of size (b-a)/(N-1).

Затем вы получаете данные за 50-й день с временным индексом 500 (и последующие 9).

Для этой медленно движущейся системы и относительно небольшого временного шага вы получите приемлемую точность с помощью реализованного метода порядка 1, но вы получите лучшую точность с помощью метода более высокого порядка, такого как RK4.

Вам необходимо везде применять связанные обновления ко всем компонентам. Это требует чередования шагов RK4, которые у вас есть, например, исправленного шага.

   S_k2 = derivada_S(time[i] + (h/2), I[i] + (h/2)*I_k1, S[i] + (h/2)*S_k1)

Требует, чтобы значение I_k1 было предварительно вычислено. Также обратите внимание, что h должен быть фактором к S_k1, его не следует добавлять.

Всего должно получиться

for i in range(data-1):
    S_k1 = derivada_S(time[i], I[i], S[i])
    I_k1 = derivada_I(time[i], I[i], S[i])
    R_k1 = derivada_R(time[i], I[i])

    S_k2 = derivada_S(time[i] + (1/2)*h, I[i] + (h/2)*I_k1, S[i] + (h/2)*S_k1)
    I_k2 = derivada_I(time[i] + (1/2)*h, I[i] + (h/2)*I_k1, S[i] + (h/2)*S_k1)
    R_k2 = derivada_R(time[i] + (1/2)*h, I[i] + (h/2)*I_k1)

    S_k3 = derivada_S(time[i] + (h/2), I[i] + (h/2)*I_k2, S[i] + (h/2)*S_k2)
    I_k3 = derivada_I(time[i] + (h/2), I[i] + (h/2)*I_k2, S[i] + (h/2)*S_k2)
    R_k3 = derivada_R(time[i] + (h/2), I[i] + (h/2)*I_k2)

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

Обратите внимание, что h является фактором для I_k1, S_k1 и т. д. У вас есть сумма.

Замена только этого фрагмента кода дает сюжет

Но перед этим есть еще одна проблема. Вы определили временной шаг как 0.05, так что t=50 достигается на последнем месте. Поскольку система автономна, содержимое массива времени не имеет значения, но маркировка оси x должна делиться на 2. Значения, которые вам нужны, на самом деле являются последними значениями, вычисленными с помощью data = 10*time_k+1.

S[-1]=0.10483, I[-1]=8.11098e-05, R[-1]=0.89509

Чтобы предыдущее обсуждение оставалось в силе, вы также можете установить h=t[1]-t[0], чтобы t=50 достигалось посередине в i=500.

Вы можете использовать интегратор, доступный на scipy.integrate.solve_ivp, и с ним использовать метод Рунге-Кутты четвертого порядка (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
'''

График вывода