【数学建模】Matlab实现SEIR模型
2024-05-03 07:04:51
1974年Hoppensteadt首先在文[1]中建立和研究了具有年龄结构的传染病模型．至今，具有年龄结构的 传染病模型的研究已有许多成果（见[2]-[5]等），但这些模型大多不考虑染病年龄、潜伏期等对疾病传播 的影响．这就不能准确的描述某些具有较长潜伏期和病程的传染病（麻疹、肺结核等）的传播,因此建立 和研究具有生理年龄、潜伏期和染病年龄的传染病模型具有重要意义．目前，同时具有生理年龄、染病 年龄、潜伏期的传染病模型研究结果较少，只有一些特殊问题的结果，如文[6]-[7]研究了具有生理年龄和 染病年龄的传染病模型,文[8]-[9]建立和研究了同时具有潜伏期、染病年龄和生理年龄的无免疫型SEIS传 染病模型及其解的适定性．本文针对麻疹、肺结核等疾病,将人群分为S（易感类）、E（潜伏类）、I（染病 类）、R（治愈类）四类人群，利用K-M仓室模型原理（见文[10]）,易建立同时具有生理年龄、染病年龄、潜 伏期的SEIR模型 。
这里总结了五个模型，分别是SI模型，SIS模型，SIR模型，SIRS模型，SEIR模型。
这几种模型的特点先介绍一下。
首先定义SEIR：
S为易感者 (Susceptible)，指未得病者，但缺乏免疫能力，与感染者接触后容易受到感染；
E为暴露者 (Exposed)，指接触过感染者，但暂无能力传染给其他人的人，对潜伏期长的传染病适用；
I为感病者 (Infective)，指染上传染病的人，可以传播给 S 类成员，将其变为 E 类或 I 类成员；
R为康复者 (Recovered)，指被隔离或因病愈而具有免疫力的人。如免疫期有限，R 类成员可以重新变为 S 类。
一、SI模型
该模型只考虑易感者和感病者，感病者不断去感染易感者。
随着时间推移，该模型感染者越来越多直到所有人都感染。
其微分方程为：
其中beta为感染率。
二、SIS模型
该模型依然只考虑易感者和感病者，感病者不断去感染易感者，这里感病者会得到治疗恢复成易感者，不过恢复后依然可能得病。
随着时间推移，该模型感染者和易感者会达到动态平衡。
其微分方程为：
其中beta为感染率，gamma为治愈率。
三、SIR模型
该模型考虑易感者、感病者与康复者，其中感病者不断感染易感者，而感病者又不断接受治疗成为康复者，康复者因为得到抗体不会再成为易感者。
随着时间推移，该模型康复者越来越多，最终所有人都成为康复者。
其微分方程为：
其中beta为感染率，gamma为治愈率。
四、SIRS模型
该模型考虑易感者、感病者与康复者，其中感病者不断感染易感者，而感病者不断接受治疗成为康复者，康复者获得抗体能够抵抗一段时间，不过最终还是会成为易感者。
该模型和SIS模型很像，区别就是康复者能够抵抗一段时间，也就是有一定的复感率，而SIS模型的复感率为1。
随着时间推移，该模型同样会达到一个动态平衡。
其微分方程为：
其中beta为感染率，gamma为治愈率，alpha为复感率。
五、SEIR模型
该模型四种参与者全部考虑，其中感病者不断感染易感者，易感者得到病毒会成为潜伏者，潜伏者依然能够使易感者感染，潜伏者有一定概率自己痊愈，感染者接受治疗也能够痊愈，痊愈后不再能够感染该病。
随着时间推移，该模型最终也是所有人都会成为康复者。
其微分方程为：
其中beta为感染率，gamma1为潜伏期康复率，gamma2为患者康复率，alpha为潜伏期转阳率。
function Message_Spread_Mode
tic
load 'Data\Link.txt'; %读入连接矩阵
% load '\Data\Point_X.txt'; %读入横坐标
% load '\Data\Point_Y.txt'; %读入纵坐标
%-------------------------------------------------------------------------%
%状态分布及状态转移概率SEIR
%0:易感状态S（Susceptible） P_0_1; （P_0_3:预免疫系数）
%1:潜伏状态E（Exposed） P_1_0；P_1_2；P_1_3
%2:染病状态I（Infected） P_2_0；P_2_3
%3:免疫状态R（Recovered） P_3_0
%-------------------------------------------------------------------------%
%计算各用户节点的度
De=sum(Link); %用户节点的度
%------------——————----参数设置与说明--------------------------------%
[M N]=size(Link); %连接矩阵的规模
I_E=0.6; %潜伏期E用户的传染强度
I_I=0.9; %发病期I用户的传染强度
lamda=sum(De)/M; %用户单位时间内平均发送信息的数量
%P_m1：用户预免疫系数
%State：用户所处状态State=zeros(1,M);0:表示易感状态（Susceptible）
%---------------------------------1---------------------------------------%
%先讨论用户预免疫系数P_m1对病毒传播的影响
TimeStep=50;%input('短信网络内病毒传播模拟时间：');
P_m1=[0.1,0.5,0.9]; %用户预免疫系数
% State=zeros(TimeStep,M); %用户的状态
G_t=5; %G_t：用户的免疫持续时间，反映了病毒的变异频率
F_t=5; %F_t：用户从发现病毒到杀毒并升级病毒库的时间
for i=1:length(P_m1)TimeLong_F=zeros(1,M); %用户处于染病期的时间长短TimeLong_E=zeros(1,M); %用户处于潜伏期的时间长短Sta=zeros(1,M); %用户的状态 %进行预免疫设定for j=1:Mif rand(1)<=P_m1(i)Sta(j)=3; %进入免疫状态TimeLong_E(j)=1; %出入潜伏期的时间为1elsecontinue;endend%状态转换%初始随机选择一个节点为病源点（此时不能选处于免疫状态的点）%问题：节点度大小存在差别，可能模拟出来的结果有出于% 为避免这个问题，我们取度最大的节点为病源节点，如果已免疫，则选次大的，一次下去[Number,Sta]=Select_Infected_Point(M,Sta,De);%Number：病源节点%State ：确定病源节点以后的节点状态矩阵State=zeros(TimeStep,M);Number_State=zeros(4,TimeStep); %用户处于个状态的统计数量for t=1:TimeStepif t==1State(t,:)=Sta;else%模拟每个用户节点的状态for j=1:M%判断用户节点处于什么状态，然后根据其状态确定其转变情况if State(t-1,j)==0 %此时处于易感状态0，可能向潜伏期转移Num=Select_Number_Near(j,Link); %找出节点j的邻居节点P=zeros(1,length(Num)); %邻居节点感染该节点的概率for k=1:length(Num)if State(t-1,Num(k))==1 %节点处于潜伏期E（1）P(k)=I_E/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./...(factorial([1:De(Num(k))]-1)));elseif State(t-1,Num(k))==2 %节点处于染病期I（2）P(k)=I_I/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./(factorial([1:De(Num(k))]-1)));elsecontinue;endendendP_0_1=max(P); %节点感染病毒的概率if rand<=P_0_1 %此时节点进入潜伏期State(t,j)=1;elseState(t,j)=State(t-1,j); endelseif State(t-1,j)==1 %此时处于潜伏状态E，可能向易感S，染病I和免疫R转移if rand<=1/(1+exp(-De(j))) %向染病状态I转移 State(t,j)=2;TimeLong_F(j)=TimeLong_F(j)+1; %用户j处于染病状态的时间长短 elseif rand<=1/(1+exp(-De(j))) %向易感状态S转移 State(t,j)=0;elseif rand<=1/(1+exp(-De(j))) %向免疫状态R转移State(t,j)=3;TimeLong_E(j)=TimeLong_E(j)+1; %免疫时间增加1elseState(t,j)=State(t-1,j); %状态不变，依然为潜伏期E（1）endendendelseif State(t-1,j)==2 %此时处于欺染病状态I，可能向易感S，免疫R转移if TimeLong_F(j)<=F_t %表示此时用户不对病毒进行任何处理State(t,j)=State(t-1,j); %此时用户维持在原状态ITimeLong_F(j)=TimeLong_F(j)+2;else%此时用户对进行杀毒并升级病毒库，进入免疫状态RState(t,j)=3;TimeLong_F(j)=0; %处于感染期（中毒状态）的时间长度TimeLong_E(j)=1; %进入免疫期的时间长度endelse%此时用户处于免疫期if TimeLong_E<=G_t %病毒此时并未突变，维持原状态R（免疫状态）State(t,j)=State(t-1,j);TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加elseif rand<=1/G_t %病毒突变，状态转移为易感状态SState(t,j)=0;TimeLong_E(j)=0;else%此时用户状态依然不变State(t,j)=State(t-1,j);TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加endendendendendendend%统计各状态的节点数量Number_State(1,t)=sum(State(t,:)==0);%处于易感状态S的总节点数量Number_State(2,t)=sum(State(t,:)==1);%处于易感状态E的总节点数量Number_State(3,t)=sum(State(t,:)==2);%处于易感状态I的总节点数量Number_State(4,t)=sum(State(t,:)==3);%处于易感状态R的总节点数量figure(i)if rem(t,3)==0plot([t-1 t],[Number_State(1,t-1) Number_State(1,t)],'md-'),hold onplot([t-1 t],[Number_State(2,t-1) Number_State(2,t)],'gh:'),hold onplot([t-1 t],[Number_State(3,t-1) Number_State(3,t)],'bs-.'),hold onplot([t-1 t],[Number_State(4,t-1) Number_State(4,t)],'k.-'),hold onelsecontinue;endlegend('易感状态Susceptible','潜伏状态Exposed','染病状态Infected','免疫状态Recovered')xlabel('模拟时间')ylabel('各状态的用户数量')end
end
P_m1=0.3; %用户预免疫系数
% State=zeros(TimeStep,M); %用户的状态
% G_t=5; %G_t：用户的免疫持续时间，反映了病毒的变异频率
G_t=[1,5,9];
F_t=5; %F_t：用户从发现病毒到杀毒并升级病毒库的时间
for i=1:length(G_t)TimeLong_F=zeros(1,M); %用户处于染病期的时间长短TimeLong_E=zeros(1,M); %用户处于潜伏期的时间长短Sta=zeros(1,M); %用户的状态 %进行预免疫设定for j=1:Mif rand(1)<=P_m1Sta(j)=3; %进入免疫状态TimeLong_E(j)=1; %出入潜伏期的时间为1elsecontinue;endend%状态转换%初始随机选择一个节点为病源点（此时不能选处于免疫状态的点）%问题：节点度大小存在差别，可能模拟出来的结果有出于% 为避免这个问题，我们取度最大的节点为病源节点，如果已免疫，则选次大的，一次下去[Number,Sta]=Select_Infected_Point(M,Sta,De);%Number：病源节点%State ：确定病源节点以后的节点状态矩阵State=zeros(TimeStep,M);Number_State=zeros(4,TimeStep); %用户处于个状态的统计数量for t=1:TimeStepif t==1State(t,:)=Sta;else%模拟每个用户节点的状态for j=1:M%判断用户节点处于什么状态，然后根据其状态确定其转变情况if State(t-1,j)==0 %此时处于易感状态0，可能向潜伏期转移Num=Select_Number_Near(j,Link); %找出节点j的邻居节点P=zeros(1,length(Num)); %邻居节点感染该节点的概率for k=1:length(Num)if State(t-1,Num(k))==1 %节点处于潜伏期E（1）P(k)=I_E/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./...(factorial([1:De(Num(k))]-1)));elseif State(t-1,Num(k))==2 %节点处于染病期I（2）P(k)=I_I/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./...(factorial([1:De(Num(k))]-1)));elsecontinue;endendendP_0_1=max(P); %节点感染病毒的概率if rand<=P_0_1 %此时节点进入潜伏期State(t,j)=1;elseState(t,j)=State(t-1,j); endelseif State(t-1,j)==1 %此时处于潜伏状态E，可能向易感S，染病I和免疫R转移if rand<=1/(1+exp(-De(j))) %向染病状态I转移 State(t,j)=2;TimeLong_F(j)=TimeLong_F(j)+1; %用户j处于染病状态的时间长短 elseif rand<=1/(1+exp(-De(j))) %向易感状态S转移 State(t,j)=0;elseif rand<=1/(1+exp(-De(j))) %向免疫状态R转移State(t,j)=3;TimeLong_E(j)=TimeLong_E(j)+1; %免疫时间增加1elseState(t,j)=State(t-1,j); %状态不变，依然为潜伏期E（1）endendendelseif State(t-1,j)==2 %此时处于欺染病状态I，可能向易感S，免疫R转移if TimeLong_F(j)<=F_t %表示此时用户不对病毒进行任何处理State(t,j)=State(t-1,j); %此时用户维持在原状态ITimeLong_F(j)=TimeLong_F(j)+2;else%此时用户对进行杀毒并升级病毒库，进入免疫状态RState(t,j)=3;TimeLong_F(j)=0; %处于感染期（中毒状态）的时间长度TimeLong_E(j)=1; %进入免疫期的时间长度endelse%此时用户处于免疫期if TimeLong_E<=G_t(i) %病毒此时并未突变，维持原状态R（免疫状态）State(t,j)=State(t-1,j);TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加elseif rand<=1/G_t(i) %病毒突变，状态转移为易感状态SState(t,j)=0;TimeLong_E(j)=0;else%此时用户状态依然不变State(t,j)=State(t-1,j);TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加endendendendendendend%统计各状态的节点数量Number_State(1,t)=sum(State(t,:)==0);%处于易感状态S的总节点数量Number_State(2,t)=sum(State(t,:)==1);%处于易感状态E的总节点数量Number_State(3,t)=sum(State(t,:)==2);%处于易感状态I的总节点数量Number_State(4,t)=sum(State(t,:)==3);%处于易感状态R的总节点数量figure(i+5)if rem(t,3)==0plot([t-1 t],[Number_State(1,t-1) Number_State(1,t)],'md-'),hold onplot([t-1 t],[Number_State(2,t-1) Number_State(2,t)],'gh:'),hold onplot([t-1 t],[Number_State(3,t-1) Number_State(3,t)],'bs-.'),hold onplot([t-1 t],[Number_State(4,t-1) Number_State(4,t)],'k.-'),hold onelsecontinue;endlegend('易感状态Susceptible','潜伏状态Exposed','染病状态Infected','免疫状态Recovered')xlabel('模拟时间')ylabel('各状态的用户数量')end
end
P_m1=0.3; %用户预免疫系数
% State=zeros(TimeStep,M); %用户的状态
% G_t=5; %G_t：用户的免疫持续时间，反映了病毒的变异频率
G_t=5;
F_t=[1,5,9]; %F_t：用户从发现病毒到杀毒并升级病毒库的时间
for i=1:length(F_t)TimeLong_F=zeros(1,M); %用户处于染病期的时间长短TimeLong_E=zeros(1,M); %用户处于潜伏期的时间长短Sta=zeros(1,M); %用户的状态 %进行预免疫设定for j=1:Mif rand(1)<=P_m1Sta(j)=3; %进入免疫状态TimeLong_E(j)=1; %出入潜伏期的时间为1elsecontinue;endend%状态转换%初始随机选择一个节点为病源点（此时不能选处于免疫状态的点）%问题：节点度大小存在差别，可能模拟出来的结果有出于% 为避免这个问题，我们取度最大的节点为病源节点，如果已免疫，则选次大的，一次下去[Number,Sta]=Select_Infected_Point(M,Sta,De);%Number：病源节点%State ：确定病源节点以后的节点状态矩阵State=zeros(TimeStep,M);Number_State=zeros(4,TimeStep); %用户处于个状态的统计数量for t=1:TimeStepif t==1State(t,:)=Sta;else%模拟每个用户节点的状态for j=1:M%判断用户节点处于什么状态，然后根据其状态确定其转变情况if State(t-1,j)==0 %此时处于易感状态0，可能向潜伏期转移Num=Select_Number_Near(j,Link); %找出节点j的邻居节点P=zeros(1,length(Num)); %邻居节点感染该节点的概率for k=1:length(Num)if State(t-1,Num(k))==1 %节点处于潜伏期E（1）P(k)=I_E/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./...(factorial([1:De(Num(k))]-1)));elseif State(t-1,Num(k))==2 %节点处于染病期I（2）P(k)=I_I/De(Num(k))*sum((lamda.^([1:De(Num(k))]).*exp(-lamda))./...(factorial([1:De(Num(k))]-1)));elsecontinue;endendendP_0_1=max(P); %节点感染病毒的概率if rand<=P_0_1 %此时节点进入潜伏期State(t,j)=1;elseState(t,j)=State(t-1,j); endelseif State(t-1,j)==1 %此时处于潜伏状态E，可能向易感S，染病I和免疫R转移if rand<=1/(1+exp(-De(j))) %向染病状态I转移 State(t,j)=2;TimeLong_F(j)=TimeLong_F(j)+1; %用户j处于染病状态的时间长短 elseif rand<=1/(1+exp(-De(j))) %向易感状态S转移 State(t,j)=0;elseif rand<=1/(1+exp(-De(j))) %向免疫状态R转移State(t,j)=3;TimeLong_E(j)=TimeLong_E(j)+1; %免疫时间增加1elseState(t,j)=State(t-1,j); %状态不变，依然为潜伏期E（1）endendendelseif State(t-1,j)==2 %此时处于欺染病状态I，可能向易感S，免疫R转移if TimeLong_F(j)<=F_t(i) %表示此时用户不对病毒进行任何处理State(t,j)=State(t-1,j); %此时用户维持在原状态ITimeLong_F(j)=TimeLong_F(j)+2;else%此时用户对进行杀毒并升级病毒库，进入免疫状态RState(t,j)=3;TimeLong_F(j)=0; %处于感染期（中毒状态）的时间长度TimeLong_E(j)=1; %进入免疫期的时间长度endelse%此时用户处于免疫期if TimeLong_E<=G_t %病毒此时并未突变，维持原状态R（免疫状态）State(t,j)=State(t-1,j);TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加elseif rand<=1/G_t %病毒突变，状态转移为易感状态SState(t,j)=0;TimeLong_E(j)=0;else%此时用户状态依然不变State(t,j)=State(t-1,j);TimeLong_E(j)=TimeLong_E(j)+1; %处于免疫期的时间增加endendendendendendend%统计各状态的节点数量Number_State(1,t)=sum(State(t,:)==0);%处于易感状态S的总节点数量Number_State(2,t)=sum(State(t,:)==1);%处于易感状态E的总节点数量Number_State(3,t)=sum(State(t,:)==2);%处于易感状态I的总节点数量Number_State(4,t)=sum(State(t,:)==3);%处于易感状态R的总节点数量figure(i+10)if rem(t,3)==0plot([t-1 t],[Number_State(1,t-1) Number_State(1,t)],'md-'),hold onplot([t-1 t],[Number_State(2,t-1) Number_State(2,t)],'gh:'),hold onplot([t-1 t],[Number_State(3,t-1) Number_State(3,t)],'bs-.'),hold onplot([t-1 t],[Number_State(4,t-1) Number_State(4,t)],'k.-'),hold onelsecontinue;endlegend('易感状态Susceptible','潜伏状态Exposed','染病状态Infected','免疫状态Recovered')xlabel('模拟时间')ylabel('各状态的人口数量')end
end
toc
【数学建模】Matlab实现SEIR模型相关推荐
- 数学建模matlab 优化模型,数学建模实验中三种优化模型的分析
1. 引言 在20世纪中期,数学建模 [1] 就在欧美国度首次被发现,而在中国的呈现稍晚些,但是大约在80年代初始咱们国家也就有了.它的核心即是创立数学模型 [2] ,使得问题获得最优化的解决.而数学 ...
- dna序列分类数学建模matlab,数学建模DNA序列分类模型(终稿).doc
您所在位置:网站首页 > 海量文档  > 高等教育 > 生物学 数学建模DNA序列分类模型(终稿).doc32页 本文 ...
- Python小白的数学建模课-09.微分方程模型
小白往往听到微分方程就觉得害怕,其实数学建模中的微分方程模型不仅没那么复杂,而且很容易写出高水平的数模论文. 本文介绍微分方程模型的建模与求解,通过常微分方程.常微分方程组.高阶常微分方程 3个案例手 ...
- 数学建模学习记录——数学规划模型
数学建模学习记录--数学规划模型 一.线性规划问题 MatLab中线性规划的标准型 MatLab中求解线性规划的命令 二.整数线性规划问题 三.非线性规划问题 MatLab中非线性规划的标准型 Mat ...
- 数学建模之减肥计划 matlab编程,数学建模matlab算法大全第02章整数规划.pdf
数学建模matlab算法大全第02章整数规划 第二章 整数规划 §1 概论 1.1 定义 规划中的变量(部分或全部)限制为整数时,称为整数规划.若在线性规划模型中, 变量限制为整数,则称为整数线性规划 ...
- 视频教程-MATLAB数学建模-Matlab
MATLAB数学建模 图像和算法等领域有多年研究和项目经验:指导发表科技核心期刊经验丰富:多次指导数学建模爱好者参赛. 宋星星 ¥12.00 立即订阅 扫码下载「CSDN程序员学院APP」,1000+ ...
- matlab已知飞机下轮廓线,数学建模matlab插值及拟合.ppt
数学建模matlab插值及拟合 2013数学建模培训 插值与拟合 插值与拟合属数值分析中函数逼近内容.在数学建模竞赛中,插值与拟合是一种常用的数据分析手段,被公认为建模中的十大算法之一. 本节首先通过 ...
- 【数学建模】常用微分方程模型 + 详细手写公式推导 + Matlab代码实现
文章目录 一.学习内容 二.学习时间 三.学习产出 3.1 微分方程基本概念 3.2 微分方程在数学建模中的应用 3.3 微分方程常用模型 3.3.1 人口增长模型 3.3.1.1 指数增长模型(马尔 ...
- 模型 趋势预测 数学建模 matlab,经济发展趋势预测数学建模.doc
经济发展趋势预测数学建模.doc 经济发展趋势预测 摘要 本题给出了从1978年到2009年该国的GDP与工业值.建筑业及农林渔业的变化的数据,对于问题1,需建立国内生产总值, Y=732.2776+ ...
最新文章
- (5.0) Tomcat 8 源码, 初始化 bootstrap
- （转）jQuery第五课：Ajax
- 不错html5画布效果，可惜网站不需要。亲们，可以看一下~
- 批量删除数据库中有特定开始字符的表、视图和存储过程
- mongodb 入门笔记
- 48-如何实现unix2dos功能
- Win11控制面板里面怎么找到系统安全？
- 卷成这样，非逼我用RTX 3090?（深度学习GPU平台种草
- 项目方案-标书-文档等等编写规范
- SQL 优化极简法则，你掌握几个？
- CGVAE -＞ Delinker -＞ DeepCoy
- 远程打开其他电脑的computer management
- python Pandas文件读写
- multisim页面不够大_微信订阅号页面改版，iOS 14.2存在严重耗电问题，贾跃亭再成被执行人，特斯拉回应被判欺诈，这就是今天的其他大新闻！...
- SSH婚庆策划系统案例
- 跟着小皮老师了解Go语言LiteIDE详细使用教程❤
- 0基础小白，如何入门数学建模？
- 【中国银联】数据挖掘笔试+三面面经
- 计算机毕业设计(50)java小程序毕设作品之校园浴室预约小程序系
- 错误java sql_java.sql.SQLSyntaxErrorException：语法错误：遇到“(”行