Simulation Python du modèle épidémique (modèle Kermack-McKendrick)
Épidémie
Citation: https://kotobank.jp/word/%E3%82%A8%E3%83%94%E3%83%87%E3%83%9F%E3%83%83%E3%82%AF-683127
> Épidémie En médecine et en santé publique, les personnes atteintes d'une certaine maladie font des prédictions normales dans une zone ou un groupe donné.
> Dépassement d'un montant important. Une maladie infectieuse telle que la grippe est répandue dans une zone spécifique.
> L'état dans lequel cela se produit simultanément dans diverses parties du monde est appelé pandémie.
Modèle Kermack-McKendrick
http://mathworld.wolfram.com/Kermack-McKendrickModel.html
Un modèle qui représente l'évolution du nombre de personnes infectées par des maladies infectieuses dans une population fermée. Ci-après, il est représenté par une équation différentielle simultanée.
\begin{eqnarray}
\frac{dx(x)}{dt}&=&-kx(t)y(t)\\
\frac{dy(t)}{dt}&=&kx(t)y(t)-ly(t)\\
\frac{dz(t)}{dt}&=&ly(t)
\end{eqnarray}
$ x (t) ≥ 0 $: nombre de personnes en bonne santé, $ y (t) ≥ 0 $: nombre de personnes malades, $ z (t) ≥ 0 $: nombre de personnes rétablies, $ k ≥ 0 $ : Taux d'infection, $ l ≥ 0 $: Taux de guérison
En tant qu'interprétation ・ Ensemble 1: Le taux de variation du nombre de personnes en bonne santé $ x $ est proportionnel au nombre de personnes en bonne santé x au nombre de personnes malades. Le coefficient de proportionnalité est de $ -k $. Puisque $ k, x, y ≥ 0 $, le nombre de personnes en bonne santé n'augmentera pas. Si le nombre de personnes malades devient 0, le nombre de personnes en bonne santé sera maintenu constant. · Type 2: Le taux de changement du nombre de malades $ y $ est déterminé par le taux de changement du nombre de personnes en bonne santé et la contribution du nombre de malades. Plus le taux de variation du nombre de personnes en bonne santé est élevé ou plus le nombre de personnes malades est élevé, plus le taux de variation du nombre de personnes malades est faible. Si le signe $ kx (t) -l $ est positif, le nombre de malades augmentera, et s'il est négatif, il diminuera. ・ Type 3: Le taux de variation du nombre de personnes qui ont récupéré $ z $ est proportionnel au nombre de personnes malades. La constante de proportionnalité est $ l $. Plus il y a de malades, plus le nombre de personnes qui guérissent augmentera rapidement.
Un point fixe
\begin{eqnarray}
\frac{dx(x)}{dt}&=&-kx(t)y(t)=0…①\\
\frac{dy(t)}{dt}&=&kx(t)y(t)-ly(t)=0…②\\
\frac{dz(t)}{dt}&=&ly(t)=0…③
\end{eqnarray}
Le point qui devient.
À partir de ② $ x $ = $ \ frac {l} {k} $ De ① et ③ $ y = 0 $
Quantité de stockage
\begin{eqnarray}
\frac{dx(x)}{dy(t)}&=&\frac{-kx(t)y(t)}{kx(t)y(t)-ly(t)}\\
&=&\frac{-kx(t)}{kx(t)-l}\\
(\frac{l}{k}\frac{1}{x}-1)dx(t)&=&dy(t)\\
\int{(\frac{l}{k}\frac{1}{x}-1)dx(t)}&=&\int{dy(t)}\\
\frac{l}{k}logx(t)-x(t)&=&y(t)+E
\end{eqnarray}
$ \ Frac {l} {k} logx (t) -x (t) -y (t) $ est la quantité de stockage (elle prend une valeur constante même si l'heure change).
Simulation par calcul numérique
Méthode de calcul
http://qiita.com/popondeli/items/391810fd727cefb190e7 La méthode de Rungecutta décrite dans est étendue et utilisée pour les équations différentielles simultanées.
·Code source
import numpy as np
import matplotlib.pyplot as plt
import math
k = 1.0
l = 1.0
DELTA_T = 0.001
MAX_T = 100.0
def dx(x, y, z, t):
return -k * x * y
def dy(x, y, z, t):
return k * x * y - l * y
def dz(x, y, z, t):
return l * y
def calc_conceved_quantity(x, y, z, t):
return (l/k)*math.log(x) - x - y
def runge_kutta(init_x, init_y, init_z, init_t):
x, y, z, t = init_x, init_y, init_z, init_t
history = [[x, y, z, t, calc_conceved_quantity(x, y, z, t)]]
cnt = 0
while t < MAX_T:
cnt += 1
k1_x = DELTA_T*dx(x, y, z, t)
k1_y = DELTA_T*dy(x, y, z, t)
k1_z = DELTA_T*dz(x, y, z, t)
k2_x = DELTA_T*dx(x + k1_x/2, y + k1_y/2, z + k1_z/2, t + DELTA_T/2)
k2_y = DELTA_T*dy(x + k1_x/2, y + k1_y/2, z + k1_z/2, t + DELTA_T/2)
k2_z = DELTA_T*dz(x + k1_x/2, y + k1_y/2, z + k1_z/2, t + DELTA_T/2)
k3_x = DELTA_T*dx(x + k2_x/2, y + k2_y/2, z + k2_z/2, t + DELTA_T/2)
k3_y = DELTA_T*dy(x + k2_x/2, y + k2_y/2, z + k2_z/2, t + DELTA_T/2)
k3_z = DELTA_T*dz(x + k2_x/2, y + k2_y/2, z + k2_z/2, t + DELTA_T/2)
k4_x = DELTA_T*dx(x + k3_x, y + k3_y, z + k3_z, t + DELTA_T)
k4_y = DELTA_T*dy(x + k3_x, y + k3_y, z + k3_z, t + DELTA_T)
k4_z = DELTA_T*dz(x + k3_x, y + k3_y, z + k3_z, t + DELTA_T)
x += (k1_x + 2*k2_x + 2*k3_x + k4_x)/6
y += (k1_y + 2*k2_y + 2*k3_y + k4_y)/6
z += (k1_z + 2*k2_z + 2*k3_z + k4_z)/6
t += DELTA_T
x = max(0, x)
y = max(0, y)
z = max(0, z)
if cnt % 100 == 0:
history.append([x, y, z, t, calc_conceved_quantity(x, y, z, t)])
return history
if __name__ == '__main__':
#Calcul numérique par la méthode Rungekutta
history = np.array(runge_kutta(init_x = 1, init_y = 0.1, init_z = 0, init_t = 0))
x_min, x_max = min(history[:,0]), max(history[:,0])
y_min, y_max = min(history[:,1]), max(history[:,1])
z_min, z_max = min(history[:,2]), max(history[:,2])
t_min, t_max = 0, MAX_T
#Tracer le diagramme de phase x vs y
plt.subplot(4, 1, 1)
plt.title("x vs y")
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.scatter(history[:,0], history[:,1])
# x(Nombre de personnes en bonne santé)Tracez le changement d'heure de
plt.subplot(4, 1, 2)
plt.title("t vs x")
plt.xlim(t_min, t_max)
plt.ylim(0, x_max)
plt.scatter(history[:,3], history[:,0])
# y(Nombre de malades)Tracez le changement d'heure de
plt.subplot(4, 1, 3)
plt.title("t vs y")
plt.xlim(t_min, t_max)
plt.ylim(0, y_max)
plt.scatter(history[:,3], history[:,1])
#Tracer la quantité de stockage au fil du temps
plt.subplot(4, 1, 4)
plt.title(u"t vs conserved quantity")
plt.xlim(t_min, t_max)
plt.scatter(history[:,3], history[:,4])
plt.show()
$ k $ et $ l $ valent 1,0.
Résultats dus aux différences de valeurs initiales x et y
Diagramme de phase X vs y lorsque la simulation est répétée 100 fois avec les valeurs initiales de $ x et y $ étant donné une valeur aléatoire de 0 à 1.
Code source
import numpy as np
import matplotlib.pyplot as plt
import math
import random
k = 1.0
l = 1.0
DELTA_T = 0.001
MAX_T = 100.0
def dx(x, y, z, t):
return -k * x * y
def dy(x, y, z, t):
return k * x * y - l * y
def dz(x, y, z, t):
return l * y
def calc_conceved_quantity(x, y, z, t):
return (l/k)*math.log(x) - x - y
def runge_kutta(init_x, init_y, init_z, init_t):
x, y, z, t = init_x, init_y, init_z, init_t
history = [[x, y, z, t, calc_conceved_quantity(x, y, z, t)]]
cnt = 0
while t < MAX_T:
cnt += 1
k1_x = DELTA_T*dx(x, y, z, t)
k1_y = DELTA_T*dy(x, y, z, t)
k1_z = DELTA_T*dz(x, y, z, t)
k2_x = DELTA_T*dx(x + k1_x/2, y + k1_y/2, z + k1_z/2, t + DELTA_T/2)
k2_y = DELTA_T*dy(x + k1_x/2, y + k1_y/2, z + k1_z/2, t + DELTA_T/2)
k2_z = DELTA_T*dz(x + k1_x/2, y + k1_y/2, z + k1_z/2, t + DELTA_T/2)
k3_x = DELTA_T*dx(x + k2_x/2, y + k2_y/2, z + k2_z/2, t + DELTA_T/2)
k3_y = DELTA_T*dy(x + k2_x/2, y + k2_y/2, z + k2_z/2, t + DELTA_T/2)
k3_z = DELTA_T*dz(x + k2_x/2, y + k2_y/2, z + k2_z/2, t + DELTA_T/2)
k4_x = DELTA_T*dx(x + k3_x, y + k3_y, z + k3_z, t + DELTA_T)
k4_y = DELTA_T*dy(x + k3_x, y + k3_y, z + k3_z, t + DELTA_T)
k4_z = DELTA_T*dz(x + k3_x, y + k3_y, z + k3_z, t + DELTA_T)
x += (k1_x + 2*k2_x + 2*k3_x + k4_x)/6
y += (k1_y + 2*k2_y + 2*k3_y + k4_y)/6
z += (k1_z + 2*k2_z + 2*k3_z + k4_z)/6
t += DELTA_T
x = max(0, x)
y = max(0, y)
z = max(0, z)
if cnt % 100 == 0:
history.append([x, y, z, t, calc_conceved_quantity(x, y, z, t)])
return history
if __name__ == '__main__':
for i in range(0, 100):
init_x = random.random()
init_y = random.random()
#Calcul numérique par la méthode Rungekutta
history = np.array(runge_kutta(init_x, init_y, init_z = 0, init_t = 0))
x_min, x_max = min(history[:,0]), max(history[:,0])
y_min, y_max = min(history[:,1]), max(history[:,1])
z_min, z_max = min(history[:,2]), max(history[:,2])
t_min, t_max = 0, MAX_T
#Tracer le diagramme de phase x vs y
plt.title("x vs y")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.scatter(history[:,0], history[:,1])
plt.show()
En dessinant une courbe, elle tombe à $ y = 0 $ (= nombre de personnes infectées 0).
Résultats dus aux différences de paramètres k et l (taux d'infection, taux de récupération)
Taux de récupération plus élevé
La chute à $ y = 0 $ devient raide
Augmentation du taux d'infection
Une fois infecté, le nombre de personnes infectées augmente et diminue à 0
Recommended Posts