模拟僵尸末日¶
| 日期 | 2011-04-24（最后修改），2011-04-24（创建） |
|---|
此示例演示了如何使用 SciPy 求解一阶常微分方程组。请注意，可以使用 SciPy 求解 N 阶方程，方法是将其转换为一阶方程组。在这个轻松的示例中，可以使用常微分方程组来模拟“僵尸入侵”，使用Munz 等人 2009 年中指定的方程。
该系统表示为
dS/dt = P - BSZ - dS dZ/dt = BSZ + GR - ASZ dR/dt = dS + ASZ - GR
使用以下符号
- S：易感受害者人数
- Z：僵尸人数
- R：被“杀死”的人数
- P：人口出生率
- d：自然死亡的可能性
- B：“僵尸病”传播的可能性（活人变成僵尸）
- G：死者复活成僵尸的可能性
- A：僵尸被完全消灭的可能性
这涉及求解由以下公式给出的常微分方程组：dy/dt = f(y, t)
其中y = [S, Z, R]。
用于解决此系统的代码如下
在 [1]
# zombie apocalypse modeling
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
plt.ion()
plt.rcParams['figure.figsize'] = 10, 8
P = 0 # birth rate
d = 0.0001 # natural death percent (per day)
B = 0.0095 # transmission percent (per day)
G = 0.0001 # resurect percent (per day)
A = 0.0001 # destroy percent (per day)
# solve the system dy/dt = f(y, t)
def f(y, t):
Si = y[0]
Zi = y[1]
Ri = y[2]
# the model equations (see Munz et al. 2009)
f0 = P - B*Si*Zi - d*Si
f1 = B*Si*Zi + G*Ri - A*Si*Zi
f2 = d*Si + A*Si*Zi - G*Ri
return [f0, f1, f2]
# initial conditions
S0 = 500. # initial population
Z0 = 0 # initial zombie population
R0 = 0 # initial death population
y0 = [S0, Z0, R0] # initial condition vector
t = np.linspace(0, 5., 1000) # time grid
# solve the DEs
soln = odeint(f, y0, t)
S = soln[:, 0]
Z = soln[:, 1]
R = soln[:, 2]
# plot results
plt.figure()
plt.plot(t, S, label='Living')
plt.plot(t, Z, label='Zombies')
plt.xlabel('Days from outbreak')
plt.ylabel('Population')
plt.title('Zombie Apocalypse - No Init. Dead Pop.; No New Births.')
plt.legend(loc=0)
# change the initial conditions
R0 = 0.01*S0 # 1% of initial pop is dead
y0 = [S0, Z0, R0]
# solve the DEs
soln = odeint(f, y0, t)
S = soln[:, 0]
Z = soln[:, 1]
R = soln[:, 2]
plt.figure()
plt.plot(t, S, label='Living')
plt.plot(t, Z, label='Zombies')
plt.xlabel('Days from outbreak')
plt.ylabel('Population')
plt.title('Zombie Apocalypse - 1% Init. Pop. is Dead; No New Births.')
plt.legend(loc=0)
# change the initial conditions
R0 = 0.01*S0 # 1% of initial pop is dead
P = 10 # 10 new births daily
y0 = [S0, Z0, R0]
# solve the DEs
soln = odeint(f, y0, t)
S = soln[:, 0]
Z = soln[:, 1]
R = soln[:, 2]
plt.figure()
plt.plot(t, S, label='Living')
plt.plot(t, Z, label='Zombies')
plt.xlabel('Days from outbreak')
plt.ylabel('Population')
plt.title('Zombie Apocalypse - 1% Init. Pop. is Dead; 10 Daily Births')
plt.legend(loc=0)
输出[1]
章节作者：ChristopherCampo