美文网首页
基于MATLAB的新冠病毒传播元胞自动机模拟

基于MATLAB的新冠病毒传播元胞自动机模拟

作者: 奈何缘浅wyj | 来源:发表于2020-09-26 11:26 被阅读0次

    鉴于2020年2月的新型冠状病毒的爆发,病毒在人群中的传播过程引起了人们的关注,这里用元胞自动机简单地模拟病毒在人群中传播的过程。

    数据的初始化

    以下代码设置GUI按键以及初始化人群分布,这里我把人群的分布作为二维上的正态分布。

    计算二维散点的正态分布可以用matlab中的mvnrnd函数

    具体用法为(以双变量为例):先构造该两个变量的均值向量,再构造协方差矩阵,输入为mvnrnd(均值向量,协方差矩阵,想得到的点的个数),该函数的输出为n行2列的矩阵,每一行中的两个元素分别代表所得的正态分布的点的行下标,列下标(假设在一个以均值点为中心无限大的矩阵中的二维正态分布投点)。

    % 设置GUI按键
    plotbutton=uicontrol('style','pushbutton','string','运行', 'fontsize',12, 'position',[150,400,50,20], 'callback', 'run=1;');
    erasebutton=uicontrol('style','pushbutton','string','停止','fontsize',12,'position',[250,400,50,20],'callback','freeze=1;');
    quitbutton=uicontrol('style','pushbutton','string','退出','fontsize',12,'position',[350,400,50,20],'callback','stop=1;close;');
    number = uicontrol('style','text','string','1','fontsize',12, 'position',[20,400,50,20]);
    total_num=100;
    n=100;
    mu=[50,50];% 均值向量
    Sigma=[200 0;0 200];% 协方差矩阵
    r=mvnrnd(mu,Sigma,total_num);
    r=fix(r);
    z = zeros(n,n);
    k=sub2ind(size(z),r(:,1),r(:,2));%得到人的坐标的序列索引
    healthy = z;inf=z;carrier=z;cardate=z;%初始化三种类型的人群的矩阵以及携带日期的矩阵
    first=randi([10,90]);
    healthy(k)=200;healthy(k(first))=0;%一开始没有感染者
    carrier(k(first))=200;cardate(k(first))=1;%初始化携带者天数矩阵
    % 建立图像句柄
    
    

    主体循环部分

    % 主事件循环
    stop= 0; run = 1;freeze = 0; 
    while stop==0
        if run==1
            %先进行感染判断
          [x11,y11]=find(inf~=0);
          m=[];n=[];
          for i=1:length(x11)
              m=[m;x11(i)-2:x11(i)+2];
              n=[n;y11(i)-2:y11(i)+2];
              m(m<1)=1;m(m>100)=100;
              n(n<1)=1;n(n>100)=100;
          end
          %m的每一行为距离感染者为2以内的行坐标
          %n的每一行为距离感染者为2以内的列坐标
          %接下来判断这个以感染者为中心5*5的矩阵中是否有健康人
          x22=[];y22=[];fm=0;
          if ~isempty(m)
          fm=length(m(:,1));
        end
          for i=1:fm
            mat=healthy(m(i,:),n(i,:));%某个感染者附近的坐标
            [x21,y21]=find(healthy(m(i,:),n(i,:))~=0);%healthy(m(i,:),n(i,:)是5*5的矩阵,其第几行就代表m那一行的第几个元素,n同理对应列下标
            x22=[x22,m(i,x21)];y22=[y22,n(i,y21)];
          end
    
          for i=1:length(x22)
          pro=rand;a=pro>0.6;b=pro<0.6;
          healthy(x22(i),y22(i))=200*a;
          carrier(x22(i),y22(i))=200*b;
          cardate(x22(i),y22(i))=b;
          end
          %得到感染者附近的健康人坐标,并以一定的概率使其成为携带者
          %接下来进行携带者变为感染者,以及携带天数增加
          date=find(cardate==4);
          carrier(date)=0;cardate(date)=0;inf(date)=200;
          cardate(cardate~=0)=cardate(cardate~=0)+1;
    
          % 现在来进行人群的移动
    
          [x1_,y1_]=find(inf~=0);infnum=length(x1_);
          inf=z;
          for i=1: infnum
          x1_(i)=randi([x1_(i)-7,x1_(i)+7]);y1_(i)=randi([y1_(i)-7,y1_(i)+7]);
          while x1_(i)<3||x1_(i)>97||y1_(i)<3||y1_(i)>97   %若碰壁则重新取点
              x1_(i)=randi([x1_(i)-7,x1_(i)+7]);y1_(i)=randi([y1_(i)-7,y1_(i)+7]);
              if x1_(i)<0%因为重新取点可能导致不停的往一个方向加或减从而永远也无法跳出while所以对于过分的点直接让其等于一个正常值
                  x1_(i)=3;
              elseif y1_(i)<0
                  y1_(i)=3;
              elseif x1_(i)>100
                  x1_(i)=97;
              elseif y1_(i)>100
                  y1_(i)=97;         
              end
          end
          end%感染者的移动
          idx1=sub2ind(size(z),x1_,y1_);%得出感染者移动后的坐标索引
          inf(idx1)=200;%将感染者的位置标注在矩阵上
    
          [x2,y2]=find(healthy~=0);healnum=length(x2);%得到健康人的矩阵行列索引,以及健康人的数量
          healthy = z;
          for i=1: healnum   
          x2(i)=randi([x2(i)-7,x2(i)+7]);y2(i)=randi([y2(i)-7,y2(i)+7]);%健康人随机运动
          while x2(i)<3||x2(i)>97||y2(i)<3||y2(i)>97||~isempty(find(idx1==sub2ind(size(z),x2(i),y2(i)), 1))%检测健康人是否与感染者的位置重叠若重叠则换位置
              x2(i)=randi([x2(i)-7,x2(i)+7]);y2(i)=randi([y2(i)-7,y2(i)+7]);
              if x2(i)<0
                  x2(i)=3;
              elseif y2(i)<0
                  y2(i)=3;
              elseif x2(i)>100
                  x2(i)=97;
              elseif y2(i)>100
                  y2(i)=97;         
              end
          end
          end
          idx2=sub2ind(size(z),x2,y2);%得到健康人移动后的下标
          healthy(idx2)=200;
    
          [x3,y3]=find(carrier~=0);carnum=length(x3);
          dateidx=sub2ind(size(z),x3,y3);
          cardate_=cardate;
          cardate=z;
          carrier=z;
          for i=1: carnum   
          x3(i)=randi([x3(i)-7,x3(i)+7]);y3(i)=randi([y3(i)-7,y3(i)+7]);%携带者随机运动
          while x3(i)<3||x3(i)>97||y3(i)<3||y3(i)>97||~isempty(find(idx1==sub2ind(size(z),x3(i),y3(i)), 1))||~isempty(find(idx2==sub2ind(size(z),x3(i),y3(i)), 1))%检测携带者是否与感染者和健康人的位置重叠若重叠则重新取点
              x3(i)=randi([x3(i)-7,x3(i)+7]);y3(i)=randi([y3(i)-7,y3(i)+7]);
              if x3(i)<0
                  x3(i)=3;
              elseif y3(i)<0
                  y3(i)=3;
              elseif x3(i)>100
                  x3(i)=97;
              elseif y1(i)>100
                  y3(i)=97;         
              end
          end
          end
          idx3=sub2ind(size(z),x3,y3);%得到携带者移动后的下标索引
          cardate(idx3)=cardate_(dateidx);
          carrier(idx3)=200;
    
            image(cat(3,inf,healthy,carrier));
            stepnumber = 1 + str2double(get(number,'string'));%更新迭代步数,每迭代一次就加一
            set(number,'string',num2str(stepnumber))%将加一后的迭代代数更新
           pause(0.7)
    
        end
        if freeze==1
            run = 0;
            freeze = 0;
        end
    end
    
    

    可以看出如果没有采取任何措施,新冠肺炎的传播在前期是越来越快的,但是当感染者数量接近总人数的时候,增长速率明显减缓。

    资源传送门

    • 关注【做一个柔情的程序猿】公众号
    • 在【做一个柔情的程序猿】公众号后台回复 【python资料】【2020秋招】 即可获取相应的惊喜哦!

    「❤️ 感谢大家」

    • 点赞支持下吧,让更多的人也能看到这篇内容(收藏不点赞,都是耍流氓 -_-)
    • 欢迎在留言区与我分享你的想法,也欢迎你在留言区记录你的思考过程。

    相关文章

      网友评论

          本文标题:基于MATLAB的新冠病毒传播元胞自动机模拟

          本文链接:https://www.haomeiwen.com/subject/lqscuktx.html