最近在研究熱力圖，發現了一篇可能有用的Python模擬疫情擴散的文章，可以部分模擬熱力圖，整篇文章原文內容如下：
瘟疫蔓延，連芬蘭都難以幸免
受傑森的《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 import numpy as np 2 import math 3 import matplotlib.pyplot as plt 4 5 from matplotlib import rcParams 6 import matplotlib.image as mpimg 7 rcParams['font.family'] = 'serif' 8 rcParams['font.size'] = 16 9 rcParams['figure.figsize'] = 12, 8 10 from PIL import Image
適當的beta和gamma值就能夠摧毀大半江山
1 beta = 0.010 2 gamma = 1
還記得導數的定義么？當導數已知，假設Δt很小的情況下，經過重新整理，它可以用來近似預測函數的下一個取值，我們已經聲明過u′(t)。
回想前面：
我們把函數(u(t +△t))在下一個時間步記為
，
表示當前時間步。
這種方法叫做歐拉法，代碼如下：
1 def euler_step(u, f, dt): 2 return u + dt * f(u)
我們需要函數f(u)。友好的numpy提供了簡潔的數組操作。我可能會在另一篇文章中回顧它，因為它們太強大了，需要更多的解釋，但現在這樣就能達到效果：
1 def f(u): 2 S = u[0] 3 I = u[1] 4 R = u[2] 5 6 new = np.array([-beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] + 7 S[0:-2, 1:-1]*I[0:-2, 1:-1] + 8 S[2:, 1:-1]*I[2:, 1:-1] + 9 S[1:-1, 0:-2]*I[1:-1, 0:-2] + 10 S[1:-1, 2:]*I[1:-1, 2:]), 11 beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] + 12 S[0:-2, 1:-1]*I[0:-2, 1:-1] + 13 S[2:, 1:-1]*I[2:, 1:-1] + 14 S[1:-1, 0:-2]*I[1:-1, 0:-2] + 15 S[1:-1, 2:]*I[1:-1, 2:]) - gamma*I[1:-1, 1:-1], 16 gamma*I[1:-1, 1:-1] 17 ]) 18 19 padding = np.zeros_like(u) 20 padding[:,1:-1,1:-1] = new 21 padding[0][padding[0] < 0] = 0 22 padding[0][padding[0] > 255] = 255 23 padding[1][padding[1] < 0] = 0 24 padding[1][padding[1] > 255] = 255 25 padding[2][padding[2] < 0] = 0 26 padding[2][padding[2] > 255] = 255 27 28 return padding
導入北歐國家的人口密度圖並進行下采樣，以便較快地得到結果
1 from PIL import Image 2 img = Image.open('popdens2.png') 3 img = img.resize((img.size[0]/2,img.size[1]/2)) 4 img = 255 - np.asarray(img) 5 imgplot = plt.imshow(img) 6 imgplot.set_interpolation('nearest')
北歐國家的人口密度圖（未包含丹麥）
S矩陣，也就是易感個體，應該近似於人口密度。感染者初始值是0，我們把斯德哥爾摩作為第一感染源。
1 S_0 = img[:,:,1] 2 I_0 = np.zeros_like(S_0) 3 I_0[309,170] = 1 # patient zero
因為還沒人死亡，所以把矩陣也置為0.
1 R_0 = np.zeros_like(S_0)
接着初始化模擬時長等。
1 T = 900 # final time 2 dt = 1 # time increment 3 N = int(T/dt) + 1 # number of time-steps 4 t = np.linspace(0.0, T, N) # time discretization 5 6 # initialize the array containing the solution for each time-step 7 u = np.empty((N, 3, S_0.shape[0], S_0.shape[1])) 8 u[0][0] = S_0 9 u[0][1] = I_0 10 u[0][2] = R_0
我們需要自定義一個顏色表，這樣才能將感染矩陣顯示在地圖上。
1 import matplotlib.cm as cm 2 theCM = cm.get_cmap("Reds") 3 theCM._init() 4 alphas = np.abs(np.linspace(0, 1, theCM.N)) 5 theCM._lut[:-3,-1] = alphas
下面坐下來欣賞吧…
1 for n in range(N-1): 2 u[n+1] = euler_step(u[n], f, dt)
讓我們再做一下圖像渲染，把它做成gif，每個人都喜歡gifs！
1 from images2gif import writeGif 2 3 keyFrames = [] 4 frames = 60.0 5 6 for i in range(0, N-1, int(N/frames)): 7 imgplot = plt.imshow(img, vmin=0, vmax=255) 8 imgplot.set_interpolation("nearest") 9 imgplot = plt.imshow(u[i][1], vmin=0, cmap=theCM) 10 imgplot.set_interpolation("nearest") 11 filename = "outbreak" + str(i) + ".png" 12 plt.savefig(filename) 13 keyFrames.append(filename) 14 15 images = [Image.open(fn) for fn in keyFrames] 16 gifFilename = "outbreak.gif" 17 writeGif(gifFilename, images, duration=0.3) 18 plt.clf()
瘟疫蔓延的gif圖，連芬蘭都難以幸免
看，唯一安全的地方是人口密度不太高的北部地區，動畫中連芬蘭最后都被感染了，現在你明白了吧。
如果你想了解更多微分方程的求解，溫馨向您推薦LorenaABarba 的實用數值方法的Python實現課程，你可以從中學習所有實數的數值求解方法，而不是本文中簡單的這種。
更新：你可以在這里找到Ipython版本的筆記。
實際運行上面的例子的代碼是存在問題的，往往會編譯不通過，下面是經過調整后的代碼例子，可以正常運行：
上圖是代碼需要的圖片popdens2.png
具體的代碼：
import numpy as np import math import matplotlib.pyplot as plt 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 = 0.010 gamma = 1 def euler_step(u, f, dt): return u + dt * f(u) 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 Image img = 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 = img[:,:,1] I_0 = np.zeros_like(S_0) I_0[309,170] = 1 # patient zero R_0 = np.zeros_like(S_0) T = 900 # final time dt = 1 # time increment N = int(T/dt) + 1 # number of time-steps t = np.linspace(0.0, T, N) # time discretization # initialize the array containing the solution for each time-step u = np.empty((N, 3, S_0.shape[0], S_0.shape[1])) u[0][0] = S_0 u[0][1] = I_0 u[0][2] = R_0 import matplotlib.cm as cm theCM = 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) 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) import imageio def create_gif(image_list, gif_name,gif_duration=0.1): frames = [] for image_name in image_list: frames.append(imageio.imread(image_name)) # Save them as frames into a gif imageio.mimsave(gif_name, frames, 'GIF', duration = gif_duration) image_list=keyFrames gif_name = 'zombie.gif' create_gif(image_list, gif_name,gif_duration=0.3) plt.clf()