瘟疫蔓延,连芬兰都难以幸免
受杰森的《Almost Looks Like Work》启发,我来展示一些病毒传播模型。需要注意的是这个模型并不反映现实情况,因此不要误以为是西非可怕的传染病。相反,它更应该被看做是某种虚构的僵尸爆发现象。那么,让我们进入主题。
这就是SIR模型,其中字母S、I和R反映的是在僵尸疫情中,个体可能处于的不同状态。
- S 代表易感群体,即健康个体中潜在的可能转变的数量。
- I 代表染病群体,即僵尸数量。
- R 代表移除量,即因死亡而退出游戏的僵尸数量,或者感染后又转回人类的数量。但对与僵尸不存在治愈者,所以我们就不要自我愚弄了(如果要把SIR模型应用到流感传染中,还是有治愈者的)。
至于β(beta)和γ(gamma):
- β(beta)表示疾病的传染性程度,只要被咬就会感染。
- γ(gamma)表示从僵尸走向死亡的速率,取决于僵尸猎人的平均工作速率,当然,这不是一个完美的模型,请对我保持耐心。
S′=−βIS告诉我们健康者变成僵尸的速率,S′是对时间的导数。
I′=βIS−γI告诉我们感染者是如何增加的,以及行尸进入移除态速率(双关语)。
R′=γI只是加上(gamma I),这一项在前面的等式中是负的。
上面的模型没有考虑S/I/R的空间分布,下面来修正一下!
一种方法是把瑞典和北欧国家分割成网格,每个单元可以感染邻近单元,描述如下:
其中对于单元,
和
是它周围的四个单元。(不要因为对角单元而脑疲劳,我们需要我们的大脑不被吃掉)。
初始化一些东东。
1 2 3 4 5 6 7 8 9 10 11 |
import numpy as np import math import matplotlib.pyplot as plt %matplotlib inline from matplotlib import rcParams import matplotlib.image as mpimg rcParams['font.family'] = 'serif' rcParams['font.size'] = 16 rcParams['figure.figsize'] = 12, 8 from PIL import Image |
适当的beta和gamma值就能够摧毁大半江山
1 2 |
beta = 0.010 gamma = 1 |
还记得导数的定义么?当导数已知,假设Δt很小的情况下,经过重新整理,它可以用来近似预测函数的下一个取值,我们已经声明过u′(t)。
回想前面:
我们把函数(u(t +△t))在下一个时间步记为,
表示当前时间步。
这种方法叫做欧拉法,代码如下:
1 2 |
def euler_step(u, f, dt): return u + dt * f(u) |
我们需要函数f(u)。友好的numpy提供了简洁的数组操作。我可能会在另一篇文章中回顾它,因为它们太强大了,需要更多的解释,但现在这样就能达到效果:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
def f(u): S = u[0] I = u[1] R = u[2] new = np.array([-beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] + S[0:-2, 1:-1]*I[0:-2, 1:-1] + S[2:, 1:-1]*I[2:, 1:-1] + S[1:-1, 0:-2]*I[1:-1, 0:-2] + S[1:-1, 2:]*I[1:-1, 2:]), beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] + S[0:-2, 1:-1]*I[0:-2, 1:-1] + S[2:, 1:-1]*I[2:, 1:-1] + S[1:-1, 0:-2]*I[1:-1, 0on-o">:-2]*I[1:-1, 0 ![]() 瘟疫蔓延,连芬兰都难以幸免 受杰森的《Almost Looks Like Work》启发,我来展示一些病毒传播模型。需要注意的是这个模型并不反映现实情况,因此不要误以为是西非可怕的传染病。相反,它更应该被看做是某种虚构的僵尸爆发现象。那么,让我们进入主题。 这就是SIR模型,其中字母S、I和R反映的是在僵尸疫情中,个体可能处于的不同状态。
至于β(beta)和γ(gamma):
S′=−βIS告诉我们健康者变成僵尸的速率,S′是对时间的导数。 I′=βIS−γI告诉我们感染者是如何增加的,以及行尸进入移除态速率(双关语)。 R′=γI只是加上(gamma I),这一项在前面的等式中是负的。 上面的模型没有考虑S/I/R的空间分布,下面来修正一下! 一种方法是把瑞典和北欧国家分割成网格,每个单元可以感染邻近单元,描述如下: 其中对于单元 初始化一些东东。
适当的beta和gamma值就能够摧毁大半江山
还记得导数的定义么?当导数已知,假设Δt很小的情况下,经过重新整理,它可以用来近似预测函数的下一个取值,我们已经声明过u′(t)。 回想前面: 我们把函数(u(t +△t))在下一个时间步记为 这种方法叫做欧拉法,代码如下:
我们需要函数f(u)。友好的numpy提供了简洁的数组操作。我可能会在另一篇文章中回顾它,因为它们太强大了,需要更多的解释,但现在这样就能达到效果:
|