seir方程求解matlab,欧拉数值法求解SEIR传染病动力学模型（MATLAB）

clc %%具体的MATLAB程序

clear

close all

%% 参数

Q=1;    % 潜伏者相对于感染者传播能力的比值

A=1/14;  %隔离时间14天

G=1/7 ; % 潜伏期向感染者的转化速率潜伏期为7d

c=2;  %接触率

p=1;   %有效接触系数

Si=0.13;  %感染者隔离速率

Sq=0.13;  % 隔离

YI=0.007;

YH=0.014;

B=2.05*10^(-9);

q=10^(-6);

a=2.7*10^(-4);

h=0.01;  %求解的步长

T=50 ;  %预测50天的

t=0:h:T;

S=59170000;% 当地总人口数1

E=4007;

I=524*1.5;

SQ=2776;

EQ=400;

H=1+EQ;

R=31;

DS=zeros(1,length(t));     DS1=zeros(1,length(t));     % 参数保存地方

DE=zeros(1,length(t));   DE1=zeros(1,length(t));

DI=zeros(1,length(t));   DI1=zeros(1,length(t));

DSQ=zeros(1,length(t));  DSQ1=zeros(1,length(t));

DEQ=zeros(1,length(t));   DEQ1=zeros(1,length(t));

DH=zeros(1,length(t));   DH1=zeros(1,length(t));

DR=zeros(1,length(t));   DR1=zeros(1,length(t));

DS(1)=S;%参数初始化

DE(1)=E;

DI(1)=I;

DSQ(1)=SQ;

DEQ(1)=EQ;

DH(1)=H;

DR(1)=R;

for i=2:length(t)

DS1(i-1)=-(p*c*B+p*c*q*(1-B))*( DS(i-1) )*( DI(i-1)+1*DE(i-1)  )+A*DSQ(i-1);

DS(i)=DS(i-1)+DS1(i-1)*h;

DE1(i-1)=p*c*B*(1-q)*DS(i-1)*( DI(i-1)+1*DE(i-1) )-G*DE(i-1);

DE(i)=DE(i-1)+DE1(i-1)*h;

DI1(i-1)=G*DE(i-1)-(Si+a+YI)*DI(i-1);

DI(i)=DI(i-1)+DI1(i-1)*h;

DSQ1(i-1)=p*c*q*(1-B)*DS(i-1)*( DI(i-1)+1*DE(i-1) )-A*DSQ(i-1);

DSQ(i)=DSQ(i-1)+DSQ1(i-1)*h;

DEQ1(i-1)=p*c*B*q*DS(i-1)*( DI(i-1)+1*DE(i-1) )-G*DE(i-1);

DEQ(i)= DEQ(i-1)+ DEQ1(i-1)*h;

DH1(i-1)=Si*DI(i-1)+Sq*DEQ(i-1)-(a+YH)*DH(i-1);

DH(i)=DH(i-1)+DH1(i-1)*h;

DR1(i-1)=(YI)*DI(i-1)+YH*DH(i-1);

DR(i)=DR(i-1)+DR1(i-1)*h;

end

subplot(3,3,1);

plot(t,DS); hold on;title('易感染总人数');

subplot(3,3,2);

plot(t,DE); hold on;title('潜伏者');

subplot(3,3,3);

plot(t,DI); hold on; title('感染人数');

subplot(3,3,4);

plot(t,DSQ); hold on; title('隔离易感者');

subplot(3,3,5);

plot(t,DEQ); hold on ; title('隔离潜伏者');

subplot(3,3,6);

plot(t,DH); hold on; title('住院患者');

subplot(3,3,7);

plot(t,DR); hold on ; title('康复人群');