Python-Simulation des Epidemiemodells (Kermack-McKendrick-Modell)
Epidemie
Zitat: https://kotobank.jp/word/%E3%82%A8%E3%83%94%E3%83%87%E3%83%9F%E3%83%83%E3%82%AF-683127
> Epidemie In der medizinischen und öffentlichen Gesundheit machen Menschen mit einer bestimmten Krankheit in einem bestimmten Bereich oder einer bestimmten Gruppe normale Vorhersagen.
> Eine große Menge überschreiten. Eine Infektionskrankheit wie Influenza ist in einem bestimmten Gebiet weit verbreitet.
> Der Zustand, in dem dies gleichzeitig in verschiedenen Teilen der Welt geschieht, wird als Pandemie bezeichnet.
Kermack-McKendrick-Modell
http://mathworld.wolfram.com/Kermack-McKendrickModel.html
Ein Modell, das Veränderungen in der Anzahl der mit Infektionskrankheiten infizierten Menschen in einer geschlossenen Population darstellt. Im Folgenden wird es durch eine simultane Differentialgleichung dargestellt.
\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 $: Anzahl gesunder Menschen, $ y (t) ≥ 0 $: Anzahl kranker Menschen, $ z (t) ≥ 0 $: Anzahl wiedergefundener Menschen, $ k ≥ 0 $ : Infektionsrate, $ l ≥ 0 $: Wiederherstellungsrate
Als Interpretation · Set 1: Die Änderungsrate der Anzahl gesunder Menschen $ x $ ist proportional zur Anzahl gesunder Menschen x der Anzahl kranker Menschen. Der Proportionalitätskoeffizient beträgt $ -k $. Da $ k, x, y ≥ 0 $ ist, wird die Anzahl gesunder Menschen nicht zunehmen. Wenn die Anzahl der Kranken 0 wird, wird die Anzahl der Gesunden konstant gehalten. · Typ 2: Die Änderungsrate der Anzahl der Kranken $ y $ wird durch die Änderungsrate der Anzahl der gesunden Menschen und den Beitrag der Anzahl der Kranken bestimmt. Je größer die Änderungsrate der Anzahl gesunder Menschen oder je größer die Anzahl der Kranken ist, desto geringer ist die Änderungsrate der Anzahl der Kranken. Wenn das Vorzeichen von $ kx (t) -l $ positiv ist, nimmt die Anzahl der Kranken zu, und wenn es negativ ist, nimmt es ab. ・ Typ 3: Die Änderungsrate der Anzahl der Personen, die $ z $ wiedererlangt haben, ist proportional zur Anzahl der Kranken. Die Proportionalitätskonstante ist $ l $. Je mehr Kranke es gibt, desto schneller steigt die Zahl der Menschen, die sich erholen.
Fixpunkt
\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}
Der Punkt, der wird.
Von ② $ x $ = $ \ frac {l} {k} $ Von ① und ③ $ y = 0 $
Lagermenge
\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) $ ist die Speichermenge (sie nimmt einen konstanten Wert an, auch wenn sich die Zeit ändert).
Simulation durch numerische Berechnung
Berechnungsmethode
http://qiita.com/popondeli/items/391810fd727cefb190e7 Die in beschriebene Rungecutta-Methode wird erweitert und für simultane Differentialgleichungen verwendet.
·Quellcode
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__':
#Numerische Berechnung nach der Rungekutta-Methode
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
#Zeichnen Sie das x gegen y-Phasendiagramm
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(Anzahl gesunder Menschen)Zeichnen Sie die Zeitänderung von
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(Anzahl der Kranken)Zeichnen Sie die Zeitänderung von
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])
#Zeichnen Sie die zeitliche Änderung der Speichermenge
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 $ und $ l $ sind 1.0.
Ergebnisse aufgrund unterschiedlicher x- und y-Anfangswerte
X-gegen-y-Phasendiagramm, wenn die Simulation 100 Mal wiederholt wird, wobei die Anfangswerte von $ x und y $ einen zufälligen Wert von 0 bis 1 erhalten.
Quellcode
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()
#Numerische Berechnung nach der Rungekutta-Methode
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
#Zeichnen Sie das x gegen y-Phasendiagramm
plt.title("x vs y")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.scatter(history[:,0], history[:,1])
plt.show()
Beim Zeichnen einer Kurve fällt sie auf $ y = 0 $ (= Anzahl der infizierten Personen 0).
Ergebnisse aufgrund unterschiedlicher Parameter k und l (Infektionsrate, Wiederherstellungsrate)
Höhere Wiederherstellungsrate
Der Abfall auf $ y = 0 $ wird steil
Erhöhte Infektionsrate
Einmal infiziert, steigt und sinkt die Anzahl der infizierten Personen auf 0