SIS病毒模型解析及模拟

解析

在SI模型基础上增加感染者 I 恢复健康的概率\(\gamma\)。
恢复后仍有可能被再次感染。得到

常微分方程

\( \frac{\mathrm{d}}{\mathrm{d}t}I=r\beta I\frac{S}{N}-\gamma I \\ \frac{\mathrm{d}}{\mathrm{d}t}S=-r\beta I\frac{S}{N}+\gamma I\)

现假设总人数N=100，感染者接触易感者的人数r=2，传染概率\(\beta\)=0.03，恢复概率\(\gamma\)=0.01，初始感染者为1即\( I\left(0\right)=1\)。得到

解

\(I\left(t\right)=\frac{250}{3\,{\left({\mathrm{e}}^{\mathrm{log}\left(\frac{247}{3}\right)-\frac{t}{20}} +1\right)}} \\ S\left(t\right)=\frac{50\,{\left({\mathrm{e}}^{\frac{t}{20}+\mathrm{log}\left(\frac{3}{247}\right)} +6\right)}}{3\,{\left({\mathrm{e}}^{\frac{t}{20}+\mathrm{log}\left(\frac{3}{247}\right)} +1\right)}}\)

为观察其变化率对其求导，得

导数

\(\frac{\mathrm{d}}{\textrm{d}t}I=\frac{25\,{\mathrm{e}}^{\mathrm{log}\left(\frac{247}{3}\right)-\frac{t}{20}} }{6\,{{\left({\mathrm{e}}^{\mathrm{log}\left(\frac{247}{3}\right)-\frac{t}{20}} +1\right)}}^2 } \\ \frac{\mathrm{d}}{\textrm{d}t}S=\frac{5\,{\mathrm{e}}^{\frac{t}{20}+\mathrm{log}\left(\frac{3}{247}\right)} }{6\,{\left({\mathrm{e}}^{\frac{t}{20}+\mathrm{log}\left(\frac{3}{247}\right)} +1\right)}}-\frac{5\,{\mathrm{e}}^{\frac{t}{20}+\mathrm{log}\left(\frac{3}{247}\right)} \,{\left({\mathrm{e}}^{\frac{t}{20}+\mathrm{log}\left(\frac{3}{247}\right)} +6\right)}}{6\,{{\left({\mathrm{e}}^{\frac{t}{20}+\mathrm{log}\left(\frac{3}{247}\right)} +1\right)}}^2 }\)

为找出增长速度最快的时间，求其导数最大值，即

最快增长时间点

\(t={I^{\prime } \left(t\right)}_{\mathrm{max}} =20\,\mathrm{log}\left(247\right)-20\,\mathrm{log}\left(3\right)\approx 88\)

即此时\(N\approx 42\)

为分析何时进入稳态，故解导函数等于零（实际上该数学模型永不为0，这里以\( {10}^{-3}\)代替），得

平衡点

\(t=20\,\mathrm{log}\left(\frac{247}{3}\right)-20\,\mathrm{log}\left(\frac{6247}{3}-\frac{50\,\sqrt{5}\,\sqrt{3122}}{3}\right)\approx 255\)

图解

通解

解析解

\( I\left(t\right)=\frac{N\,{\left(\gamma -\beta \,r\right)}}{{\mathrm{e}}^{\mathrm{log}\left(\frac{N\,\gamma +I_0 \,\beta \,r-N\,\beta \,r}{I_0 }\right)+\gamma \,t-\beta \,r\,t} -\beta \,r} \\ S\left(t\right)=-\frac{N\,{\left(\gamma -{\mathrm{e}}^{\mathrm{log}\left(\frac{N\,\gamma -S_0 \,\beta \,r}{N-S_0 }\right)+\gamma \,t-\beta \,r\,t} \right)}}{{\mathrm{e}}^{\mathrm{log}\left(\frac{N\,\gamma -S_0 \,\beta \,r}{N-S_0 }\right)+\gamma \,t-\beta \,r\,t} -\beta \,r}\)

导数

\(\frac{d}{\textrm{d}t}I=-\frac{N\,{\mathrm{e}}^{\mathrm{log}\left(\frac{N\,\gamma +I_0 \,\beta \,r-N\,\beta \,r}{I_0 }\right)+\gamma \,t-\beta \,r\,t} \,{{\left(\gamma -\beta \,r\right)}}^2 }{{{\left({\mathrm{e}}^{\mathrm{log}\left(\frac{N\,\gamma +I_0 \,\beta \,r-N\,\beta \,r}{I_0 }\right)+\gamma \,t-\beta \,r\,t} -\beta \,r\right)}}^2 } \\ \\ 
\begin{array}{*{35}{l}}
\frac{d}{\text{d}t}S=\frac{N\,{{\sigma }_{1}}\,\left( \gamma -\beta \,r \right)}{{{\sigma }_{1}}-\beta \,r}+\frac{N\,{{\sigma }_{1}}\,\left( \gamma -{{\sigma }_{1}} \right)\,\left( \gamma -\beta \,r \right)}{{{\left( {{\sigma }_{1}}-\beta \,r \right)}^{2}}} \\
其中 \\
{{\sigma }_{1}}={{\text{e}}^{\log \left( \frac{N\,\gamma -{{S}_{0}}\,\beta \,r}{N-{{S}_{0}}} \right)+\gamma \,t-\beta \,r\,t}} \\
\end{array}\)

增长最快时间点

\( t=\frac{\mathrm{log}\left(-\beta \,r\right)-\mathrm{log}\left(\frac{N\,\gamma +I_0 \,\beta \,r-N\,\beta \,r}{I_0 }\right)+2\,\pi \,k\,\mathrm{i}}{\gamma -\beta \,r} \\
\left(\begin{array}{c}
N\not= 0\wedge {\mathrm{e}}^{\mathrm{log}\left(\frac{N\,\gamma +I_0 \,\beta \,r-N\,\beta \,r}{I_0 }\right)+\frac{\gamma \,\sigma_1 }{\gamma -\beta \,r}-\frac{\beta \,r\,\sigma_1 }{\gamma -\beta \,r}} \not= \beta \,r\wedge k\in \mathbb{Z}\wedge \gamma \not= \beta \,r\\
\sigma_1 =\mathrm{log}\left(-\beta \,r\right)-\mathrm{log}\left(\frac{N\,\gamma +I_0 \,\beta \,r-N\,\beta \,r}{I_0 }\right)+2\,\pi \,k\,\mathrm{i}
\end{array}\right)\)

平衡点

取\(0={10}^{-3}\)

\(\begin{align}
& t=\frac{\log \left( \beta \,r-{{\sigma }_{3}}-{{\sigma }_{2}}-500\,N\,\gamma \,{{\sigma }_{6}}+{{\sigma }_{5}}+{{\sigma }_{4}} \right)-{{\sigma }_{1}}+2\,\pi \,k\,\text{i}}{\gamma -\beta \,r} \\
& 其中 \\
& {{\sigma }_{1}}=\log \left( \frac{N\,\gamma +{{I}_{0}}\,\beta \,r-N\,\beta \,r}{{{I}_{0}}} \right) \\
& {{\sigma }_{2}}=500\,N\,{{\beta }^{2}}\,{{r}^{2}} \\
& {{\sigma }_{3}}=500\,N\,{{\gamma }^{2}} \\
& {{\sigma }_{4}}=1000\,N\,\beta \,\gamma \,r \\
& {{\sigma }_{5}}=500\,N\,\beta \,r\,{{\sigma }_{6}} \\
& {{\sigma }_{6}}=\sqrt{-\frac{-250\,N\,{{\beta }^{2}}\,{{r}^{2}}+500\,N\,\beta \,\gamma \,r+\beta \,r-250\,N\,{{\gamma }^{2}}}{250\,N}} \\
\end{align}\)

查看所有模型及说明：点击病毒传染的动态扩散模拟及数学模型

网络模拟

更新中.....

图示

Matlab代码

笛卡尔坐标系模拟

图示(动态过程未给出)

Matlab代码

clear
figure
N = 100;        %人口总数
I = 1;          %传染者
S = N - I;      %易感者

R = 2;          %感染者接触易感者的人数
B = 0.03;       %传染概率
Y = 0.01;       %康复概率
days = 300;
Dynamic = false;%设置为true显示动态模拟过程

Sdata=[];

state=zeros(1,N);
index=randperm(N,I);
state(index)=1;
if Dynamic
    axis([0 N 0 N]);
    hold on
end
for i=1:days
    xpos=randperm(N);
    ypos=randperm(N);
    if sum(state)==100
        break
    end
    for j=1:N
        if state(j)~=1
            continue
        end
        if rand()<=Y
            state(j)=0;
            continue
        end
        dis=sqrt((xpos-xpos(j)).^2+(ypos-ypos(j)).^2);
        l=0;
        for k=1:N
            if state(k)==1
                continue
            end
            l=l+1;
            peo(l).dis=dis(k);
            peo(l).num=k;
        end
        T = struct2table(peo);
        sortedT = sortrows(T,'dis');
        sortedS = table2struct(sortedT);
        tp=sortedT{1:R,2};
        for k=1:length(tp)
            if rand()<=B
                state(tp(k))=1;
            end
        end
    end
    Sdata(i)=sum(state);
    %a=-8;
    %b=8;
    %xpos=xpos+ a + (b-a).*rand(1,N);
    %ypos=ypos+ a + (b-a).*rand(1,N);
    if Dynamic && rem(i,10)==0
        phd=scatter(xpos,ypos,[],state,"filled");
        title(['第' num2str(i) '天'])
        drawnow
        delete(phd)
    end
end

figure
set(gcf,'visible',true)
MaxFigure
days=length(Sdata);
plot(1:days,Sdata,1:days,N-Sdata)
hold on
T = 1:days+100;
for idx = 1:length(T)-1
    if S(idx)<0
        S(idx)=0;
    elseif S(idx)>N
        S(idx)=N;
    end
    if I(idx)<0
        I(idx)=0;
    elseif I(idx)>N
        I(idx)=N;
    end
    S(idx+1) = S(idx) - R*B*I(idx)*S(idx)/N + Y*I(idx);
    I(idx+1) = I(idx) + R*B*I(idx)*S(idx)/N - Y*I(idx);
end
plot(T,S,T,I);
legend('模拟传染者','模拟易感者','易感者','传染者');
xlabel('天');ylabel('人数')
title('SIS模型')Code language: Matlab (matlab)