数学建模sir模型python_Python实现简单的SI传播模型

Python实现简单的SI传播模型

9 W1 N' t% O% P% ]9 M#SI疾病传播模型的原理9 k0 t. @& x2 |7 h+ M9 ~* b9 C* X

在经典的传染病模型中，种群(Population)内N个个体的状态可分为如下几类$ K8 m: {+ j- `0 e

4 M) O A4 `: V1 O& E U) H" }易感状态(Susceptible)。一个个体在感染前是处于易感状态的，即该个体有可能被邻居个体感染。

% v* X0 d2 ?8 y, Y' ] p H易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。，即该个体还会以一定概率感染其邻居个体。, y$ b) A b, V. w& q

移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态，当一个个体经历过一个完整的感染周期后，该个体就不再被感染，因此就可以不再考虑改革提。, i& W- F! [! g' H1 U6 c( M

SI传播模型是最简单的疾病传播模型，模型中的所有个体都只可能处于两个状态中的一个

9 C6 U' U" x6 t/ B即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。$ F% e2 @6 C% L) }7 r9 N' _5 I2 q3 R

在给定时刻t，令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目，显然有

- z+ e; B6 f: v* q6 N+ PS(t)+I(t)恒等于N，这里，N是个体总数。随着时间t的增长，易感个体与感染个体的接触6 E; w& R. s9 L' J* f3 T

会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β，疾病仅在; ?+ z5 r) `9 l$ D) P

感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t，易感个体的比例为S(t)/N，感染个体的数量为I(t)，一次，易感个体的数量将以如下变化率减少4 Z* N! G# z3 o3 \4 a" u

ds/dt = -β*S(t)I(t)/N

! d6 e4 y2 L) D2 j同时，感染个体的数量会以与易感个体相反的变化率增加，

" S# S( X" Z; Ids/dt = βS(t)*I(t)/N

+ O! N, ?! ]% c2 ^/ j分别将时刻t处于易感状态和感染状态的个体所占比例记为，

: U7 ^! p. B* Bs(t)=S(t)/N5 M3 X$ M) b) T+ Q; }( {

i(t)=I(t)/N6 X. B/ Z; n N+ c9 E2 D7 w

显然有，- O) u @3 M% u5 f

s(t)+i(t)恒等于1，此时之前的公式可以记做9 O( Y" B) P! N- k! s) N2 C

ds/dt=-βsi

2 y$ d( k, J6 n8 R& v* Wdi/dt=βsi

; E P+ N3 F3 S; B" M即

% P0 x7 P7 N. B+ l7 l7 }di/dt=βi(1-i)( t& @1 G* g* I

上式也成为Logistic增长方程式(Logistic growth equation),

4 `3 I( K/ {1 f1 _& G/ O. B方程的解和图像如图9 x- j# @" g( i6 J* k. _

2020-4-18 16:14 上传

下载附件 (213.53 KB)

) _& _+ h% s3 K" C5 u& p( e! b! J代码和相关文件以及环境链接：链接：https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ

% e+ k$ p) k6 `" e: L" S提取码：z448

% {* z: G7 ?4 q; [2 L

0 E' l' j9 l* P' m9 A) M9 ]$ z2 ~0 Q/ a0 {. p+ @

'''

4 t9 E" Q4 }) ?实验环境Python2.7.13,igraph包，cairo包，numpy包

4 t$ q# o U9 X% A V'''

1 Z: ?2 ]& |7 ~. ?4 T# -*- coding:utf8 -*

9 y- ~/ K% b2 |, R% L2 Hfrom igraph import *

6 j, Q1 U% X3 y# y% b: F) Z6 yimport numpy as numpy0 ^2 F' z G9 q

from numpy import *% J7 {( ?- t, C: Q

import random

# m; n: M/ ~9 D* k0 T

$ d! ~7 Q* J' W: ^: ?, Tdef len_arr(infected_array,nodes_num):#获取感染数组长度

8 `3 i3 }4 ]; P* H: s5 {3 j1 M len_value=0#初始化长度

7 i' z1 I9 X1 y4 P' w len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)

/ U8 f% x: M7 F( j' r) h return len_value# r0 y7 D* h2 R6 m2 A

. f0 R; M/ g- B- m* t* P4 Ng=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)( O0 M" J7 l A7 q6 z" `; f9 m

summary(g)# Y2 B& L4 O' g. L( K! V1 d% o8 K' w

nodes_num=g.vcount()#统计图中的结点个数

$ h) k' c/ |9 y' t3 O X/ P, knet_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat- F. J& f% `7 P* U. A" Q

g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色+ @7 _% ^/ u" C K

a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a

* ]8 @% M O3 l- D1 mnodes_state=matrix(a).T#nodes_state通过转置a矩阵创建，用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重！！！)6 Q6 X4 V; M" ?

#第一列是节点编号，第二列是节点状态，感染状态用-2表示，第三列是节点感染的时间# K# {# _$ }0 q$ p

print(nodes_state)

% g( W$ H) ~5 s/ F( cinfected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染 34代表网络节点数

! R! {1 `# x$ s2 `print(infected_array)9 Z8 Q, F( y4 v+ q/ f% ?

+ G& x! n E7 w! s/ ?$ ^( Y* {0 J- l8 ?infe_rate=1#传播率(感染率) 1代表邻接点100%被感染

/ j% y' j, i( y: b, Lset_time=2#传播次数(感染次数) 2次% K$ z( F( _8 l

source_seed=1#感染源位置

' O, N- _0 a1 N' `1 m* A3 W$ c% wnodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1

9 p, A/ e; ?8 Pnodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态

! Q; s7 u. |7 M E& @1 E( onodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1

$ D% P5 A+ v0 H2 p. _5 Hg.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红% y8 [& ], b% j& e2 o O

infected_array[0]=source_seed#将感染源的位置存入被感染节点列表3 G8 B$ m n8 ]4 r# }

plot(g)#绘制

, B( {& B0 c, l+ {: i* K j" Y7 F, E' A

stop=False#感染过程结束的标记+ l) r8 x9 t" y$ s) P

temp_time=0#第几次感染

, Z+ ~. x+ m+ `( u# H4 k" rtemp_len=0#本轮的感染源数量初始化- |" ], w# z4 `( y* r, } N) ~/ j

, _4 w( P2 ~2 j. M( v8 Xwhile not stop:# P& P1 e4 q* I: i

i=0#记录让每个感染源都传播一次

; m* m: y3 E( L1 ]* N4 m if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行7 U1 i ]& ?5 Y8 T

temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量

$ F/ _3 K- C4 T4 v$ ]9 a- R while i

7 e0 y- j, V2 v' c# F- Y temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间' y* b% g' S0 X1 k

nei_count=0#下一轮可以被感染到的节点数量

# m y# n' q* j- Y# ^ #生成下一轮可能被感染的节点的集合nei_arr

- ~0 T f! ?+ G/ U2 F3 Q for j in range(nodes_num):#遍历节点

8 f+ |5 c7 i3 k- T: R if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染4 x. p0 y: A3 d$ e v' b+ V$ D

nei_count=nei_count+1#下一轮可以被感染到的节点数量++. A+ m8 d+ D) }/ k" a

nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染. Q7 W, K9 L6 e" `' V8 B! r

t=0

+ {3 [: n4 t- F3 w$ n: V. u for j in range(nodes_num):

' f- L5 b8 h' E5 O. c* {) F if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:3 H1 |& L J; m( w

nei_arr[t]=j+15 V) f/ a' b2 \ L7 Q! r2 M

t=t+1

$ k: q* r6 v/ v- d/ l; i. y ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组4 ^# E6 x& B4 u9 c0 W

#random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象

. z/ R* l; S3 ^& ?0 l if len(ran_infe_arr)>0:#存在需要被感染的节点; Z6 Z4 P: N2 g2 I; X1 U

t=0#让ran_infe_arr内每个感染源都被感染9 I8 t3 O& r. z0 i

while t

* T q6 V& c+ s* { nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态, u: `% L* Q7 }4 ?8 k

nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间

# F* M2 g* ^4 K1 U& r* g infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中$ L8 T6 h: O" A0 ~+ ~, X" e

g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色7 v+ q+ C8 w$ Y2 \, d. @

plot(g)#绘制

" v; C# H# J4 @ t=t+1. u# }& F+ g4 h) ^

i=i+1

: j2 \" ^/ ?$ w if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染

3 b0 h B) w" ?; `1 C. Z: Z stop=True - ~& r2 G+ Y. r1 L% Z, ^

- x% y, C3 f1 |( y/ S0 P# C0 R3 ?& ^. ?7 p

视频演示bilibili传送门- g E( t+ L# c9 l+ e$ Z; K7 p

效果图

# S) s2 @3 S! x8 E' C

1 M! Z+ M0 e$ l

) G6 z1 P: i8 n! J

2020-4-18 16:15 上传

下载附件 (174.22 KB)

0 y+ G# M/ P; X

; L/ |( N2 h; G+ _

2020-4-18 16:15 上传

下载附件 (22.44 KB)

6 Z7 t, V4 P* R) D% V5 ]1 J# C# a! r1 q+ e* \- y( j' i9 f

2020-4-18 16:16 上传

下载附件 (42.85 KB)

: \" L" @1 l" _$ K' E* N. ]

\1 ?. b( N; E0 L0 H0 @, T `* y% z

2020-4-18 16:16 上传

下载附件 (17.44 KB)

9 B$ n8 d% E+ x. c3 W

( ^+ Y& c- b6 {5 }

2020-4-18 16:16 上传

下载附件 (19.78 KB)

/ g- X8 I: a! v n% k+ C

* m# e; D, V& X6 i3 G% ?$ U

2020-4-18 16:16 上传

下载附件 (21.16 KB)

————————————————7 z' i) m% X. e7 H- t5 \: D

原文链接：https://blog.csdn.net/wdays83892469/article/details/80878862

! N% ^7 P5 m- E

8 _) @$ w5 B9 `; a1 T9 U- y1 @; X+ a( S( \6 O0 J6 N! f$ @

原文链接：https://blog.csdn.net/weixin_42297843/article/details/113719185