受杰森的《Almost Looks Like Work》啟發，我來展示一些病毒傳播模型。需要注意的是這個模型并不反映現實情況，因此不要誤以為是西非可怕的傳染病。相反，它更應該被看做是某種虛構的僵尸爆發現象。那么，讓我們進入主題。
這就是SIR模型，其中字母S、I和R反映的是在僵尸疫情中，個體可能處于的不同狀態。
上面的模型沒有考慮S/I/R的空間分布，下面來修正一下！
一種方法是把瑞典和北歐國家分割成網格，每個單元可以感染鄰近單元，描述如下：
其中對于單元，和是它周圍的四個單元。（不要因為對角單元而腦疲勞，我們需要我們的大腦不被吃掉）。
初始化一些東東。
import numpy as npimport mathimport matplotlib.pyplot as plt %matplotlib inlinefrom matplotlib import rcParamsimport matplotlib.image as mpimgrcParams['font.family'] = 'serif'rcParams['font.size'] = 16rcParams['figure.figsize'] = 12, 8from PIL import Image
適當的beta和gamma值就能夠摧毀大半江山
beta = 0.010gamma = 1
還記得導數的定義么？當導數已知，假設Δt很小的情況下，經過重新整理，它可以用來近似預測函數的下一個取值，我們已經聲明過u′(t)。
初始化一些東東。
import numpy as npimport mathimport matplotlib.pyplot as plt %matplotlib inlinefrom matplotlib import rcParamsimport matplotlib.image as mpimgrcParams['font.family'] = 'serif'rcParams['font.size'] = 16rcParams['figure.figsize'] = 12, 8from PIL import Image
適當的beta和gamma值就能夠摧毀大半江山
beta = 0.010gamma = 1
還記得導數的定義么？當導數已知，假設Δt很小的情況下，經過重新整理，它可以用來近似預測函數的下一個取值，我們已經聲明過u′(t)。
這種方法叫做歐拉法，代碼如下：
def euler_step(u, f, dt): return u + dt * f(u)
我們需要函數f(u)。友好的numpy提供了簡潔的數組操作。我可能會在另一篇文章中回顧它，因為它們太強大了，需要更多的解釋，但現在這樣就能達到效果：
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, 0:-2] + S[1:-1, 2:]*I[1:-1, 2:]) - gamma*I[1:-1, 1:-1], gamma*I[1:-1, 1:-1] ]) padding = np.zeros_like(u) padding[:,1:-1,1:-1] = new padding[0][padding[0] < 0] = 0 padding[0][padding[0] > 255] = 255 padding[1][padding[1] < 0] = 0 padding[1][padding[1] > 255] = 255 padding[2][padding[2] < 0] = 0 padding[2][padding[2] > 255] = 255 return padding
導入北歐國家的人口密度圖并進行下采樣，以便較快地得到結果
from PIL import Imageimg = Image.open('popdens2.png')img = img.resize((img.size[0]/2,img.size[1]/2))img = 255 - np.asarray(img)imgplot = plt.imshow(img)imgplot.set_interpolation('nearest')
北歐國家的人口密度圖（未包含丹麥）
S矩陣，也就是易感個體，應該近似于人口密度。感染者初始值是0，我們把斯德哥爾摩作為第一感染源。
S_0 = img[:,:,1]I_0 = np.zeros_like(S_0)I_0[309,170] = 1 # patient zero
因為還沒人死亡，所以把矩陣也置為0.
R_0 = np.zeros_like(S_0)
接著初始化模擬時長等。
T = 900 # final timedt = 1 # time incrementN = int(T/dt) + 1 # number of time-stepst = np.linspace(0.0, T, N) # time discretization # initialize the array containing the solution for each time-stepu = np.empty((N, 3, S_0.shape[0], S_0.shape[1]))u[0][0] = S_0u[0][1] = I_0u[0][2] = R_0
我們需要自定義一個顏色表，這樣才能將感染矩陣顯示在地圖上。
import matplotlib.cm as cmtheCM = cm.get_cmap("Reds")theCM._init()alphas = np.abs(np.linspace(0, 1, theCM.N))theCM._lut[:-3,-1] = alphas下面坐下來欣賞吧…
for n in range(N-1): u[n+1] = euler_step(u[n], f, dt)
讓我們再做一下圖像渲染，把它做成gif，每個人都喜歡gifs！
from images2gif import writeGif keyFrames = []frames = 60.0 for i in range(0, N-1, int(N/frames)): imgplot = plt.imshow(img, vmin=0, vmax=255) imgplot.set_interpolation("nearest") imgplot = plt.imshow(u[i][1], vmin=0, cmap=theCM) imgplot.set_interpolation("nearest") filename = "outbreak" + str(i) + ".png" plt.savefig(filename) keyFrames.append(filename) images = [Image.open(fn) for fn in keyFrames]gifFilename = "outbreak.gif"writeGif(gifFilename, images, duration=0.3)plt.clf()新聞熱點
疑難解答
圖片精選
