目录
1 知识回顾
?
2 岂曰无衣，与子同裳
位卑未敢忘忧国：
?
3 简单的冠状病毒模型（Matlab代码实现）
?
% Cellular Automata
% Corona
clear all;
n = 200; x = linspace(0,1,n); y = linspace(0,1,n);
p_InitialInfection = 0.05; lambda_0 = 0.25; Tend = 40; lambda_r = 0.15;
lambda_contact = 3;
kmax = 1000;
Dt = Tend / (kmax + 1);
Times = linspace(0,Tend,kmax+1);
TDeath = 5; Tsusc = 10;
ProbabilityInfection = zeros(n,n);
pcontact = 1 - exp(-lambda_contact *Dt);
%[X,YY] = meshgrid(x,y);Y = YY';
%S = round(rand(n,n)); Sinit = S;
S = ones(n,n);
I = zeros(n,n);
R = zeros(n,n);
D = zeros(n,n);
InfectionTime = zeros(n,n);
for i = 1:n
for j = 1:n
if (x(j)-0.5)^2 + (y(i)-0.5)^2 < 0.05
chi = rand;
if chi <= p_InitialInfection
I(i,j) = 1; S(i,j) = 0;
InfectionTime(i,j) = InfectionTime(i,j) + Dt;
end;
end;
end;
end;
Sinit = S;
%S = ones(n,n);
TotalNumberS = zeros(kmax+1,1);
TotalNumberI = TotalNumberS;TotalNumberR = TotalNumberS;
TotalNumberD = TotalNumberS;
TotalNumberS(1) = sum(sum(S))/n^2;
TotalNumberI(1) = sum(sum(I))/n^2;
TotalNumberR(1) = sum(sum(R))/n^2;
TotalNumberD(1) = sum(sum(D))/n^2;
RandomContacts = false;
if RandomContacts
for k = 1:kmax
SumNeigh = zeros(n,n);
Iorig = I;
Sorig = S;
for i = 1:n
for j = 1:n
if i > 1
SumNeigh(i,j) = SumNeigh(i,j) + heaviside(pcontact-rand)*Iorig(i-1,j);
end;
if i < n
SumNeigh(i,j) = SumNeigh(i,j) + heaviside(pcontact-rand)*Iorig(i+1,j);
end;
if j > 1
SumNeigh(i,j) = SumNeigh(i,j) + heaviside(pcontact-rand)*Iorig(i,j-1);
end;
if j < n
SumNeigh(i,j) = SumNeigh(i,j) + heaviside(pcontact-rand)*Iorig(i,j+1);
end;
if (Sorig(i,j) == 1)
ProbabilityInfection(i,j) = 1 - exp(-lambda_0*SumNeigh(i,j)*Dt);
chi = rand;
if chi < ProbabilityInfection(i,j)
S(i,j) = 0;I(i,j) = 1;
Infection(i,j) = Infection(i,j) + Dt;
end;
end;
if (Iorig(i,j) == 1)
ProbabilityRecovery(i,j) = 1 - exp(-lambda_r*Dt);
chi = rand;
if chi < ProbabilityRecovery(i,j)
I(i,j) = 0; R(i,j) = 1;InfectionTime(i,j) = 0;
else
InfectionTime(i,j) = InfectionTime(i,j) + Dt;
if InfectionTime(i,j) > TDeath
I(i,j) = 0; D(i,j) = 1;
end;
end;
end;
end;
end;
TotalNumberS(k+1) = sum(sum(S))/n^2;
TotalNumberI(k+1) = sum(sum(I))/n^2;
TotalNumberR(k+1) = sum(sum(R))/n^2;
TotalNumberD(k+1) = sum(sum(D))/n^2;
end;
else
for k = 1:kmax
SumNeigh = zeros(n,n);
Iorig = I;
Sorig = S;
for i = 1:n
for j = 1:n
if i > 1
SumNeigh(i,j) = SumNeigh(i,j) + Iorig(i-1,j);
end;
if i < n
SumNeigh(i,j) = SumNeigh(i,j) + Iorig(i+1,j);
end;
if j > 1
SumNeigh(i,j) = SumNeigh(i,j) + Iorig(i,j-1);
end;
if j < n
SumNeigh(i,j) = SumNeigh(i,j) + Iorig(i,j+1);
end;
if (Sorig(i,j) == 1)
ProbabilityInfection(i,j) = 1 - exp(-lambda_0*SumNeigh(i,j)*Dt);
chi = rand;
if chi < ProbabilityInfection(i,j)
S(i,j) = 0;I(i,j) = 1;
InfectionTime(i,j) = InfectionTime(i,j) + Dt;
end;
end;
if (Iorig(i,j) == 1)
ProbabilityRecovery(i,j) = 1 - exp(-lambda_r*Dt);
chi = rand;
if chi < ProbabilityRecovery(i,j)
I(i,j) = 0; R(i,j) = 1;InfectionTime(i,j) = 0;
else
InfectionTime(i,j) = InfectionTime(i,j) + Dt;
if InfectionTime(i,j) > TDeath
I(i,j) = 0; D(i,j) = 1;
end;
end;
end;
end;
end;
TotalNumberS(k+1) = sum(sum(S))/n^2;
TotalNumberI(k+1) = sum(sum(I))/n^2;
TotalNumberR(k+1) = sum(sum(R))/n^2;
TotalNumberD(k+1) = sum(sum(D))/n^2;
end;
end;
%figure(1); semilogy(TotalNumber);hold on;
figure(1); plot(Times,TotalNumberS,'k','linewidth',3);hold on;
plot(Times,TotalNumberI,'r','linewidth',3);
plot(Times,TotalNumberR,'g','linewidth',3);
plot(Times,TotalNumberD,'b','linewidth',3);
figure(2)
plot(Times,TotalNumberI,'r','linewidth',3);hold on;
return;
figure
hold on
for i = 1:n
for j = 1:n
if Sinit(i,j) == 1
plot(x(i),y(j),'ko');
else
plot(x(i),y(j),'rx');
end;
end;
end;
hold off
return;
figure
hold on
for i = 1:n
for j = 1:n
if Sinit(i,j) == 1
plot(x(i),y(j),'ko');
else
plot(x(i),y(j),'rx');
end;
end;
end;
hold off
figure
plot(TotalNumber);
return
figure
hold on
for i = 1:n
for j = 1:n
if SumNeigh(i,j) == 1
plot(x(j),y(i),'ro');
elseif SumNeigh(i,j) == 2
plot(x(j),y(i),'ko');
elseif SumNeigh(i,j) == 3
plot(x(j),y(i),'rx');
elseif SumNeigh(i,j) == 4
plot(x(j),y(i),'kx');
elseif SumNeigh(i,j) == 0
plot(x(j),y(i),'bo');
end;
end;
end;
hold off
?
