Python-Simulation of the Epidemic Model (Kermack-McKendrick Model)
Epidemic
Quote: https://kotobank.jp/word/%E3%82%A8%E3%83%94%E3%83%87%E3%83%9F%E3%83%83%E3%82%AF-683127
> Epidemic In medical and public health, people with a certain disease make normal predictions in a certain area or population.
> Exceeding a large amount. An infectious disease such as influenza is prevalent in a specific area.
> The state in which this occurs simultaneously in various parts of the world is called a pandemic.
Kermack-McKendrick model
http://mathworld.wolfram.com/Kermack-McKendrickModel.html
A model that represents changes in the number of people infected with infectious diseases in a closed population. Hereinafter, it is represented by simultaneous differential equations.
\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 $: number of healthy people, $ y (t) ≥0 $: number of sick people, $ z (t) ≥0 $: number of recovered people, $ k ≥0 $ : Infection rate, $ l ≥ 0 $: Recovery rate
As an interpretation ・ Set 1: The rate of change in the number of healthy people $ x $ is proportional to the number of healthy people x the number of sick people. The proportional coefficient is $ -k $. Since $ k, x, y ≥ 0 $, the number of healthy people will not increase. If the number of sick people becomes 0, the number of healthy people will be kept constant. ・ Type 2: The rate of change in the number of sick people $ y $ is determined by the rate of change in the number of healthy people and the contribution of the number of sick people. The greater the rate of change in the number of healthy people or the greater the number of sick people, the lower the rate of change in the number of sick people. If the sign of $ kx (t) -l $ is positive, the number of sick people will increase, and if it is negative, it will decrease. ・ Formula 3: The rate of change in the number of people who have recovered $ z $ is proportional to the number of sick people. The constant of proportionality is $ l $. The more sick people there are, the faster the number of people who recover will increase.
Fixed point
\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}
The point that becomes.
From ② $ x $ = $ \ frac {l} {k} $ From ① and ③ $ y = 0 $
Conserved quantity
\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) $ is the conserved quantity (it takes a constant value even if the time changes).
Simulation by numerical calculation
Method of calculation
http://qiita.com/popondeli/items/391810fd727cefb190e7 The Runge-Kutta method described in is extended and used for simultaneous differential equations.
·Source code
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__':
#Numerical calculation by Runge-Kutta method
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
#Plot x vs y phase diagram
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(Number of healthy people)Plot the time change of
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(Number of sick people)Plot the time change of
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])
#Plot the time change of the conserved quantity
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 $ and $ l $ are 1.0.
Results due to differences in x and y initial values
X vs y phase diagram when the simulation is repeated 100 times with the initial values of $ x and y $ given a random value from 0 to 1.
Source code
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()
#Numerical calculation by Runge-Kutta method
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
#Plot x vs y phase diagram
plt.title("x vs y")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.scatter(history[:,0], history[:,1])
plt.show()
While drawing a curve, it falls to $ y = 0 $ (= number of infected people 0).
Results due to differences in parameters k and l (infection rate, recovery rate)
Higher recovery rate
The fall to $ y = 0 $ becomes steep
Increased infection rate
Once infected, the number of infected people increases and decreases to 0.