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

传染病模型的元胞自动机实现

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-2

q7=0.8; %5-0

lambda1=0.01075; %患者的死亡概率

lambda2=0.00525; %隔离者的死亡概率

%% 初始化元胞

people=zeros(n); %初始化人群

people(rand(n)

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:1

S=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,[n

1: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:n

for 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)] );

end

end

%===演化过程===

people(S==0&rand(n)

people(S==0&step5>=T5&rand(n)

people(S==6&rand(n)

else

people(S==6)=0;

end

people(S==1&rand(n)

people((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)

people(S==2&rand(n)

%2-4患者变成死亡者

people(S==3&step3>=T3|(S==3&step3<=T3&rand(n)

people(S==3&rand(n)

people(S==5&step4>=T4&rand(n)

%部分治愈者变成易感者

% 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:popu

if people(move1(1,i),move1(2,i))~=4 && people(move2(1,i),move2(2,i))~=4

pc=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; %变的是状态，位置不变

end

end

%===求各类人数和===

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);

drawnow

end

for k=2:L

S=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,[n

1: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:n

for 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)] );

end

end

%===演化过程===

people(S==0&rand(n)

people(S==1&rand(n)

people((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)

people(S==2&rand(n)

%2-4患者变成死亡者

people(S==3&step3>=T3|(S==3&step3<=T3&rand(n)

people(S==3&rand(n)

people(S==5&step4>=T4&rand(n)

%部分治愈者变成易感者

% 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:popu

if people(move1(1,i),move1(2,i))~=4 && people(move2(1,i),move2(2,i))~=4

pc=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; %变的

是状态，位置不变

end

end

%===求各类人数和===

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);

drawnow

end

x=1:L;

plot(x,tt1,x,tt2,x,tt3,x,tt4)

legend('潜伏者','染病者','隔离者','死亡者')