传染病模型的元胞自动机实现(附程序)

传染病模型的元胞自动机实现clc;clear all;close all;TT0=[];TT1=[];TT2=[];TT3=[];TT4=[];TT5=[];TT6=[]; %========构建元胞网络============n=100; %网格规模L=73; %仿真天数，365/5=73%=============参数设定===========p1=0.001; %0-1（易感者-潜伏者）yita=0.001; %5-0（治愈着-易感者）Step=10; %迁移步长m=0.1; %迁移比例popu=round(m*n); %每次迁移的人数round(X)朝最近的方向取整，对X中的每个元素朝最近的方向取整T1=1.4; %1-3（潜伏者-隔离者时间步长）?T2=2; %2-3(染病着-隔离者时间步长) ?T3=3; %隔离时间步长T4=73; %免疫时间步长T5=15; %接种疫苗时间步长(10-14天表现出症状) ?q1=0.3; %潜伏者强制隔离阀值?q2=0.5; %染病着强制隔离阀值?q3=0.6; %隔离者强制治愈阀值q4=0.9; %接种疫苗力度?q5=0.9; %6-5（接疫苗者-治愈者）q6=0.5; %1-2q7=0.8; %5-0lambda1=0.01075; %患者的死亡概率lambda2=0.00525; %隔离者的死亡概率%% 初始化元胞people=zeros(n); %初始化人群people(rand(n)<p1)=1; %初始化人群rand(n)生成一个n*n阶随机矩阵,元素均匀分布在区间（0,1）上tt0=[]; tt1=[];tt2=[];tt3=[];tt4=[];tt5=[];tt6=[];step1=zeros(n); %状态变化后，时间清0 step2=zeros(n);step3=zeros(n);step4=zeros(n);step5=zeros(n);for k=1:1S=people; %记录当前人群状态S易感者fi=rand(n); %患者的感染力Ri=randn(n); %易感者抵抗力，randn(n)生成n阶随机矩阵，元素服从均值为0，方差为1的标准正态分布Ri(S~=0)=0; %仅易感者有抵抗力fi(S~=1)=0; %仅患者有传染力%===各步长记录===step1(S==1)=step1(S==1)+1; % 1 表示刷新一次step2(S==2)=step2(S==2)+1;step3(S==3)=step3(S==3)+1;step4(S==5)=step4(S==5)+1;step5(S==0)=step5(S==0)+1;%===各步长清零===step1(S~=1)=0;step2(S~=2)=0;step3(S~=3)=0;step4(S~=5)=0;step5(S~=0)=0;%===各传染概率===%每个元胞共有8种（上下左右，东南、西南、东北、西北）邻居，每种邻居构成一个矩阵ne1=sqrt(1.*fi([n 1:n-1],1:n).*(1-Ri));ne2=sqrt(1.*fi([2:n1],1:n).*(1-Ri));ne3=sqrt(1.*fi(1:n,[n1:n-1]).*(1-Ri));ne4=sqrt(1.*fi(1:n,[2:n 1]).*(1-Ri));ne5=sqrt(0.7.*fi([n 1:n-1],[2:n 1]).*(1-Ri));ne6=sqrt(0.7.*fi([2:n 1],[2:n 1]).*(1-Ri));ne7=sqrt(0.7.*fi([2:n 1],[n 1:n-1]).*(1-Ri));ne8=sqrt(0.7.*fi([n 1:n-1],[n 1:n-1]).*(1-Ri));%===选出上述矩阵中各个位置的最大值===for i=1:nfor j=1:n%Pij(i,j)=ne1(i,j)+ne2(i,j)+ne3(i,j)+ne4(i,j)+...% ne5(i,j)+ne6(i,j)+ne7(i,j)+ne8(i,j);Pij(i,j)=max([ne1(i,j),ne2(i,j),ne3(i,j),ne4(i,j),...ne5(i,j),ne6(i,j),ne7(i,j),ne8(i,j)] );endend%===演化过程===people(S==0&rand(n)<Pij)=1; %%0-1people(S==0&step5>=T5&rand(n)<q4)=6; %%0-6 if S==6&rand(n)<q5people(S==6&rand(n)<q5)=5; %%6-5elsepeople(S==6)=0;endpeople(S==1&rand(n)<q6)=2; %%1-2people((S==1&step1>=T1)|(S==1&step1<=T1&rand(n)< q1))=3; %1-3 潜伏者变成隔离者people((S==2&step2>=T2)|(S==2&step2<=T2&rand(n) <q2))=3; %%2-3people(S==2&rand(n)<lambda1)=4;%2-4患者变成死亡者people(S==3&step3>=T3|(S==3&step3<=T3&rand(n) <q3))=5; %3-5people(S==3&rand(n)<lambda2)=4; %3-4people(S==5&step4>=T4&rand(n)<q7)=0; %5-0%部分治愈者变成易感者% people(S==6&rand(n)>=q5)=0; %%6-0%==%=迁移===mo1=randperm(n); %随机顺序randperm(n)返回一个随机排列的整数1:n mo2=randperm(n);move1=[mo1(1:popu);mo2(1:popu)]; % 每次移动元胞的位置（i,j）move2=ceil((rand(2,popu)-0.5)*Step); %迁移长度%2个随机数，移动用的ceil朝正无穷大方向取整，rand(n,m)生成一个n*m阶随机矩阵move2=move1+move2;move2(move2>n)=move2(move2>n)-n; % 减去move2(move2<1)=move2(move2<1)+n;%for i=1:popuif people(move1(1,i),move1(2,i))~=4 && people(move2(1,i),move2(2,i))~=4pc=people(move1(1,i),move1(2,i));people(move1(1,i),move1(2,i))=people(move2(1,i),move2( 2,i));people(move2(1,i),move2(2,i))=pc; %变的是状态，位置不变endend%===求各类人数和===t0=length(find(people==0)); % length(X)计算向量或矩阵的长度,参量X是向量，直接返回其长度；X为非空矩阵，返回矩阵X的行数和列数值中的较大值t1=length(find(people==1)); % find(k) 返回按行检索稀疏矩阵X中非零元素的"位置索引"向量k 146页t2=length(find(people==2));t3=length(find(people==3));t4=length(find(people==4));t5=length(find(people==5));t6=length(find(people==6));tt0=[tt0 t0];tt1=[tt1 t1];tt2=[tt2 t2];tt3=[tt3 t3];tt4=[tt4 t4];tt5=[tt5 t5];tt6=[tt6 t6];% % ==================画图过程================image(people*20);title(sprintf('Recovered = %.0f,Death = %.0f',t5,t4),'fontsize',15);xlabel(sprintf('t in %.0f day',k),'fontsize',13);drawnowendfor k=2:LS=people; %记录当前人群状态S易感者fi=rand(n); %患者的感染力Ri=randn(n); %易感者抵抗力，randn(n)生成n阶随机矩阵，元素服从均值为0，方差为1的标准正态分布Ri(S~=0)=0; %仅易感者有抵抗力fi(S~=1)=0; %仅患者有传染力%===各步长记录===step1(S==1)=step1(S==1)+1; % 1 表示刷新一次step2(S==2)=step2(S==2)+1;step3(S==3)=step3(S==3)+1;step4(S==5)=step4(S==5)+1;step5(S==0)=step5(S==0)+1;%===各步长清零===step1(S~=1)=0;step2(S~=2)=0;step3(S~=3)=0;step4(S~=5)=0;step5(S~=0)=0;%===各传染概率===%每个元胞共有8种（上下左右，东南、西南、东北、西北）邻居，每种邻居构成一个矩阵ne1=sqrt(1.*fi([n 1:n-1],1:n).*(1-Ri));ne2=sqrt(1.*fi([2:n 1],1:n).*(1-Ri));ne3=sqrt(1.*fi(1:n,[n1:n-1]).*(1-Ri));ne4=sqrt(1.*fi(1:n,[2:n 1]).*(1-Ri));ne5=sqrt(0.7.*fi([n 1:n-1],[2:n 1]).*(1-Ri));ne6=sqrt(0.7.*fi([2:n 1],[2:n 1]).*(1-Ri));ne7=sqrt(0.7.*fi([2:n 1],[n1:n-1]).*(1-Ri));ne8=sqrt(0.7.*fi([n 1:n-1],[n 1:n-1]).*(1-Ri));%===选出上述矩阵中各个位置的最大值===for i=1:nfor j=1:n%Pij(i,j)=ne1(i,j)+ne2(i,j)+ne3(i,j)+ne4(i,j)+...% ne5(i,j)+ne6(i,j)+ne7(i,j)+ne8(i,j);Pij(i,j)=max([ne1(i,j),ne2(i,j),ne3(i,j),ne4(i,j),...ne5(i,j),ne6(i,j),ne7(i,j),ne8(i,j)] );endend%===演化过程===people(S==0&rand(n)<Pij)=1; %%0-1people(S==1&rand(n)<q6)=2; %%1-2people((S==1&step1>=T1)|(S==1&step1<=T1&rand(n)< q1))=3; %1-3 潜伏者变成隔离者people((S==2&step2>=T2)|(S==2&step2<=T2&rand(n) <q2))=3; %%2-3people(S==2&rand(n)<lambda1)=4;%2-4患者变成死亡者people(S==3&step3>=T3|(S==3&step3<=T3&rand(n) <q3))=5; %3-5people(S==3&rand(n)<lambda2)=4; %3-4people(S==5&step4>=T4&rand(n)<q7)=0; %5-0%部分治愈者变成易感者% people(S==6&rand(n)>=q5)=0; %%6-0 %==%=迁移===mo1=randperm(n); %随机顺序randperm(n)返回一个随机排列的整数1:n mo2=randperm(n);move1=[mo1(1:popu);mo2(1:popu)]; % 每次移动元胞的位置（i,j）move2=ceil((rand(2,popu)-0.5)*Step); %迁移长度%2个随机数，移动用的ceil朝正无穷大方向取整，rand(n,m)生成一个n*m阶随机矩阵move2=move1+move2;move2(move2>n)=move2(move2>n)-n; % 减去move2(move2<1)=move2(move2<1)+n;%for i=1:popuif people(move1(1,i),move1(2,i))~=4 && people(move2(1,i),move2(2,i))~=4pc=people(move1(1,i),move1(2,i));people(move1(1,i),move1(2,i))=people(move2(1,i),move2( 2,i));people(move2(1,i),move2(2,i))=pc; %变的是状态，位置不变endend%===求各类人数和===t0=length(find(people==0)); % length(X)计算向量或矩阵的长度,参量X是向量，直接返回其长度；X为非空矩阵，返回矩阵X的行数和列数值中的较大值t1=length(find(people==1)); % find(k) 返回按行检索稀疏矩阵X中非零元素的"位置索引"向量k 146页t2=length(find(people==2));t3=length(find(people==3));t4=length(find(people==4));t5=length(find(people==5));t6=length(find(people==6));tt0=[tt0 t0];tt1=[tt1 t1];tt2=[tt2 t2];tt3=[tt3 t3];tt4=[tt4 t4];tt5=[tt5 t5];tt6=[tt6 t6];% % ==================画图过程================image(people*20);title(sprintf('Recovered = %.0f,Death = %.0f',t5,t4),'fontsize',15);xlabel(sprintf('t in %.0f day',k),'fontsize',13);drawnowendx=1:L;plot(x,tt1,x,tt2,x,tt3,x,tt4)legend('潜伏者','染病者','隔离者','死亡者')。