战胜疫情——简单的冠状病毒模型（Matlab代码实现）

目录

1 知识回顾

2 岂曰无衣，与子同裳

3 简单的冠状病毒模型（Matlab代码实现）

---

1 知识回顾

新冠病毒传播模拟（Matlab实现）

 

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