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

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

# 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, применяя четвертый порядок Рунге-Кутты, но я не знаю, как собрать данные для времени = 50.

«Вам нужно решить связанную систему как связанную систему», см. /runge-kutta-problems-in‌​-js/… и многие другие.

Lutz Lehmann 21.11.2022 06:57
Мутабельность и переработка объектов в Python
Мутабельность и переработка объектов в Python
Объекты являются основной конструкцией любого языка ООП, и каждый язык определяет свой собственный синтаксис для их создания, обновления и...
Другой маршрут в Flask Python
Другой маршрут в Flask Python
Flask - это фреймворк, который поддерживает веб-приложения. В этой статье я покажу, как мы можем использовать @app .route в flask, чтобы иметь другую...
14 Задание: Типы данных и структуры данных Python для DevOps
14 Задание: Типы данных и структуры данных Python для DevOps
Проверить тип данных используемой переменной, мы можем просто написать: your_variable=100
Python PyPDF2 - запись метаданных PDF
Python PyPDF2 - запись метаданных PDF
Python скрипт, который будет записывать метаданные в PDF файл, для этого мы будем использовать PDF ридер из библиотеки PyPDF2 . PyPDF2 - это...
Переменные, типы данных и операторы в Python
Переменные, типы данных и операторы в Python
В Python переменные используются как место для хранения значений. Пример переменной формы:
Почему Python - идеальный выбор для проекта AI и ML
Почему Python - идеальный выбор для проекта AI и ML
Блог, которым поделился Harikrishna Kundariya в нашем сообществе Developer Nation Community.
1
1
56
3
Перейти к ответу Данный вопрос помечен как решенный

Ответы 3

Эта ссылка может помочь вам построить модель модели 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()

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

Вы повторяете ошибку вопроса, реализуете сложный метод первого порядка, а не РК4.

Lutz Lehmann 21.11.2022 06:58

Даже не метод первого порядка, S[i]+h+0.5*S_k1 должно быть S[i]+h*0.5*S_k1, а I_k2 = func_I(t[i]+0.5*h, I[i], S[i]+h+0.5*I_k1) должно быть I_k2 = func_I(t[i]+0.5*h, I[i]+h*0.5*I_k1, S[i]+h*0.5*S_k1) и так далее.

Lutz Lehmann 21.11.2022 14:15
Ответ принят как подходящий

Самый простой способ получить значение в момент времени 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.

Спасибо за ответ! Я изменил переменную времени на time = np.linspace(0, time_k, data), где time_k = 100 и data = int(time_k*10 + 1). Для этого я взял S[500] = 0,053, I[500] = 0,016, R[500] = 1,0044. Это нормально? Я думал, что он не может иметь значения больше 1.

ksa 21.11.2022 13:38

Это происходит, если вы реализуете не RK4, а слегка связанный метод порядка 1. Так что это не приятно, но и не удивительно.

Lutz Lehmann 21.11.2022 13:45

Вы можете использовать интегратор, доступный на 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
'''

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

Другие вопросы по теме