2021-07-22MATLAB基于元胞自动机模型的传染病扩散模型

MATLAB基于元胞自动机模型的传染病扩散模型

蓝色小点代表未感染过的健康人，蓝色圆圈代表康复者，红色星号代表感染者

基本思路：

1. 地图矩阵可以分为三类，第一类是健康人矩阵，第二类是感染者矩阵，第三类是康复者矩阵。健康人感染后则落入感染者矩阵，感染者康复后则升上康复者矩阵，且康复者对病毒具有免疫能力。允许多个人在同一点。
2. 通过建立三维矩阵Track（在x轴和y轴之上加入时间轴形成三维）来追踪已感染的人的感染时间，从而判定康复。
3. 通过设置NewMap和NewPatientMap以及NewSurvivalMap三个过渡矩阵来更新每次的随机移动，防止先更新的人在遍历过程中再次被选中导致一个回合内多次移动。
4. 第一次遍历让三类人随机走动并判定死亡，第二次遍历只判定感染，最后一次遍历只判定康复。这样安排能尽量减少遍历次数并避免bug产生。

%将地图分为三种，第一种是健康人地图，第二种是病人地图，第三种是康复者地图，多个人可以在同一个点。 %初始化 pplnum = 1000; %地图内总人数 V = 1; %所有人的移动速度 R = 5; %病人传染健康人的范围半径 P = 0.6; %健康人处在风险区域被感染的概率 T = 100; %康复所需时间 D = 0.01; %病人每步的死亡概率 Size = 201; %地图大小 Map = zeros(Size); %正常人分布地图 NewMap = zeros(Size);%正常人分布地图过渡矩阵，消除循环中部分人移动多次的bug PatientMap = zeros(Size); %病人分布地图 NewPatientMap = zeros(Size); %病人分布地图过渡矩阵 SurvivalMap = zeros(Size); %康复者地图 NewSurvivalMap = zeros(Size); %康复者过渡矩阵 direction=[0 -V; 0 V;V 0; -V 0]; %上下左右四个方向 DangerousArea = []; %存储会被感染的邻域 NewPatientNum = 1; %存储本轮新产生的病人数量 Track = zeros(Size,Size,T); %记录病人行踪，方便判定康复 Time = 0; %总循环次数，用来记录动画 %使人群随机分布 for i = 1:pplnum-1 xi = randperm(Size,1); yi = randperm(Size,1); Map(xi,yi) = Map(xi,yi)+1; end %随机选择零号病人位置 for i = 1:1 PatientMap(randperm(Size,1),randperm(Size,1)) = 1; end Track(:,:,1) = PatientMap; %寻找会被感染的区域 [X,Y] = meshgrid(-R:R); for i = 1:(2*R+1)^2 if X(i)^2 + Y(i)^2 <= R^2 DangerousArea = [DangerousArea;X(i),Y(i)]; end end %主程序 while sum(PatientMap,"all") > 0 %在病人清零前一直循环 Time = Time+1; %循环次数增加，这一步是方便动画帧的捕捉 %在第一次遍历中，完成健康人、康复者、病人随机走动以及病人死亡判定 for i = 1:Size %建立for循环嵌套，遍历图中所有点 for k = 1:Size %健康人随机走动 for num = 1:Map(i,k) r1 = randi(4); %四个方向随机选一个走 NewY = min(max(i+direction(r1,1),1),Size); %新的坐标 NewX = min(max(k+direction(r1,2),1),Size); NewMap(NewY,NewX) = NewMap(NewY,NewX)+1; %人移到新位置 end %康复者随机走动 for num = 1:SurvivalMap(i,k) r3 = randi(4); %四个方向随机选一个走 NewY3 = min(max(i+direction(r3,1),1),Size); %新的坐标 NewX3 = min(max(k+direction(r3,2),1),Size); NewSurvivalMap(NewY3,NewX3) = NewSurvivalMap(NewY3,NewX3)+1;%人移到新位置 end %让本轮感染前的病人随机走动 for num = 1:PatientMap(i,k) if rand<D %死了就直接删除Track数据，不进行随机移动就不会进入NewPatientMap for n3 = T:-1:1 %倒序在Track中查找 if Track(i,k,n3) ~= 0 Track(i,k,n3) = Track(i,k,n3)-1; break end end else %没死再进行随机移动 r2 = randi(4); %四个方向随机选一个走 NewY2 = min(max(i+direction(r2,1),1),Size); %新的坐标 NewX2 = min(max(k+direction(r2,2),1),Size); NewPatientMap(NewY2,NewX2) = NewPatientMap(NewY2,NewX2)+1;%人移到新位置 for n2 = T:-1:1 %倒序，从最早感染的开始查找 if Track(i,k,n2) ~= 0 Track(i,k,n2) = Track(i,k,n2)-1; %改变轨迹位置 Track(NewY2,NewX2,n2) = Track(NewY2,NewX2,n2)+1; break end end end end end end %位置数据已经全部转移至过渡矩阵，因此以下操作全部在New地图上进行 %第二次遍历只完成判定感染 for i = 1:Size %建立for循环嵌套，遍历图中所有点 for k = 1:Size if NewPatientMap(i,k) ~= 0 %如果此处有病人 for n = 1:size(DangerousArea,1) %选取该病人周围空间 npy = min(max(i+DangerousArea(n,1),1),Size); npx = min(max(k+DangerousArea(n,2),1),Size); %选取周围空间的一点 for n1 = 1:NewMap(npy,npx) %如果该点有健康人，则遍历健康人 if rand<P %有P的概率此健康人感染 NewPatientMap(npy,npx) = NewPatientMap(npy,npx)+1; %加入新病人 NewMap(npy,npx) = NewMap(npy,npx)-1; %从健康人地图剔除新病人 Track(npy,npx,1) = Track(npy,npx,1)+1; %记录新病人位置 end end end end end end %第三次遍历判定康复 for i = 1:Size %建立for循环嵌套，遍历图中所有点 for k = 1:Size for n8 = 1:Track(i,k,end) NewSurvivalMap(i,k) = NewSurvivalMap(i,k)+1; %增加康复者 NewPatientMap(i,k) = NewPatientMap(i,k)-1; %减少病人 end end end %其他操作 Map = NewMap; %将新地图赋值到原地图 PatientMap = NewPatientMap; SurvivalMap = NewSurvivalMap; NewMap = zeros(Size); %清空过渡矩阵 NewPatientMap = zeros(Size); NewSurvivalMap = zeros(Size); Track(:,:,2:end) = Track(:,:,1:end-1); %将轨迹整体后移一个单位时间 Track(:,:,1) = zeros(Size); [pply,pplx] = find(Map); %找到当前人的坐标 [ptty,pttx] = find(PatientMap); [svvy,svvx] = find(SurvivalMap); plot(pplx-1,pply-1,'b.') %绘图 axis([-10 Size+9 -10 Size+9]); hold on plot(pttx-1,ptty-1,'r*') plot(svvx-1,svvy-1,'bo') hold off M(:,Time) = getframe; clf end movie(M,1,50) disp(sum(SurvivalMap,'all')+sum(Map,'all')) %统计幸存者 

笔者第一次编程，希望各位也能多多批评，给出建议。