传染病sir模型matlab案例代小村亚姆

Question

1. 英格兰德比郡的小村亚姆（Eyam）有一个别号，叫“鼠疫之村”。但这个称呼并非耻辱，而是一种荣耀。亚姆村在 166566 年发现鼠疫时，但为了不祸害附近的其他村庄，村民们自发隔离，不让外人进入，同时里面的人也不能出去。附近地区的瘟疫渐渐得到了控制，但亚姆村民为此做出了巨大的牺牲，在瘟疫肆虐长达一年多的时间，全村 350 中只有 83 人幸存。
   亚姆村的瘟疫源于 16641666 年伦敦的大瘟疫。1665 年 9 月初，村里的裁缝收到了一包从伦敦寄来的布料，4 天后他死了。月底又有 5 人死亡，人们这才发现原来衣料中混入了带着鼠疫病菌的跳蚤。这些鼠蚤可以通过叮咬受害者传播鼠疫。此后，瘟疫开始感染其他受害者，死亡人数见表 1。

                                             表1 亚姆村第一波鼠疫死亡人数
   
   从上表中可以看出，瘟疫死亡似乎正在减少。但是在 1666 年夏天开始的时候，瘟疫又重新出现了。接下来 5 个月的死亡情况见表 2。
                                             表2 亚姆村第二波鼠疫死亡人数  
                     
   

   请用 SIR 模型研究亚姆村第二波鼠疫的传播过程，列出数学模型，并用MATLAB求解，进行结果分析。
   亚姆村第二波鼠疫起始于 1666 年 6 月 18 日，此时亚姆村总人口为 N=261。如果以 1666 年 6 月 18 日为时间起点（t=0），则亚姆村不同时间死亡人数见表 3。

                                         表3 亚姆村鼠疫不同时间（天）死亡人数
   

1666年（月份）
时间：06, 07, 08, 09, 10
人数：19, 56, 77, 24, 14

时间,15.5,31,46.5,62,77.5,93,108.5,124,254,235,201,154,121,108,97,未知,83,15,22,29,20,未知

Answer 1

让阿豪来帮你解答，本回答参考chatgpt3.5编写提供，如果还有疑问可以评论或留言

数学模型： SIR模型是一种描述传染病流行的经典模型，它将整个人群分为三类：易感者（S），感染者（I），康复者（R）。假设亚姆村中的人群以SIR模型为基础，在第二波鼠疫传播过程中，易感者(S)会通过感染者(I)而变成感染者，感染者(I)随时间会康复变成康复者(R)，不再感染他人，R中的人不再易感、不再感染和不再康复。设I(t)表示t时刻感染者的数量，S(t)表示t时刻易感者的数量，R(t)表示t时刻康复者的数量，则有以下数学模型（其中β和γ是病毒传播率和恢复率）：

初值为S(0)=261-83=178，I(0)=83,R(0)=0。 MATLAB代码：根据上述模型，编写MATLAB程序求解。

function [T,S,I,R]=SIR_model(beta,gamma,N,T_inc,t_death,D_death,t_unknown,d)
%SIR_model(beta,gamma,N,T_inc,t_death,D_death,t_unknown,d)求解SIR模型
%输入变量：
%beta：感染率
%gamma：康复率
%N：总人口数
%T_inc: 时间间隔
%t_death：发记载死亡人数的日期
%D_death：发记载死亡人数的日期之前的死亡人数
%t_unknown：发记载死亡人数的日期之后的死亡人数
%d: 死亡人数数组
%输出变量：
%T：时间
%S：易感者
%I：感染者
%R：康复者
%初始化变量：
S(1)=N-D_death(1)-d(1);
I(1)=d(1);
R(1)=0;
T(1)=0;
t=1;
%求解SIR模型：
while I(t)>=0
    t=t+1;
    T(t)=T_inc*(t-1);
    S(t)=(beta*S(t-1)*I(t-1)/N)*T_inc+S(t-1);
    I(t)=((beta*S(t-1)*I(t-1))/N-gamma*I(t-1))*T_inc+I(t-1);
    R(t)=(gamma*I(t-1))*T_inc+R(t-1);
    if t>length(D_death) %如果发现未知死亡人数，则SIR计算终止
        break
    elseif T(t)>=t_death %如果T(t)>=t_death才执行下一行代码，即到了发记录死亡人数的时间
        I(t)=I(t)-d(t-t_death); %由于死亡人数中都是感染者，减去死亡人数相当于感染者转化为了康复者
        R(t)=R(t)+d(t-t_death);
    end
end
%结果可视化：
figure;
plot(T,S,'b',T,I,'r',T,R,'g','LineWidth',1.5);
xlabel('Time (days)'); ylabel('Number of people');
title(['SIR model (beta=',num2str(beta),', gamma=',num2str(gamma),', N=',num2str(N),')']);
end

结果分析： 调用SIR_model函数，设定参数如下，即β=0.3925，γ=0.35，N=261，t_death=46.5，D_death=196，t_unknown=77.5，d=[20 未知 未知 未知 83 15 22 29 20 未知]，T_inc=0.5。

beta=0.3925;
gamma=0.35;
N=261;
t_death=46.5;
D_death=196;
t_unknown=77.5;
d=[20 NaN NaN NaN 83 15 22 29 20 NaN];
T_inc=0.5;
[T,S,I,R]=SIR_model(beta,gamma,N,T_inc,t_death,D_death,t_unknown,d);

运行结果如下图所示。 从图中可以看到，第二波鼠疫在1666年6月18日开始爆发，并持续了将近4个月，到10月初死亡人数才逐渐降低。一开始易感者人数很多，但感染者的数量也会随时间增加。到了感染者数量达到峰值后，由于病毒传播率β小于恢复率γ，感染者的数量开始减少，康复者的数量慢慢增加。在8月末，感染者数量达到顶峰，此时易感者已经基本被感染，感染者数量维持一段时间然后开始下降，康复者数量随之增加，到10月初，感染者数量降至较少，康复者数量增加较多，易感者数量维持不变。总的来说，SIR模型定量描述了亚姆村第二波鼠疫的传播过程，结果也与实际数据较为吻合。