數模-微分方程(SEIR模型)
阿新 • • 發佈:2022-05-09
模型
程式碼1
fun1.m
function dx=fun1(t,x) % 大家可以修改裡面的引數,來看結果的變化 beta = 0.1; % 易感染者與已感染者接觸且被傳染的強度 sigma = 0.2; % 潛伏者轉換為感染者的速率 gamma = 0.02; % 康復率 dx = zeros(4,1); % x(1)表示S x(2)表示E x(3)表示I x(4)表示R C = x(1)+x(2)+x(3); % 傳染病系統中的有效人群,也就是課件中的N' = S+E+I dx(1) = - beta*x(1)*x(3)/C; dx(2) = beta*x(1)*x(3)/C - sigma*x(2); dx(3) = sigma*x(2) - gamma*x(3); dx(4) = gamma*x(3); end
code.m
%% 最簡單的SEIR模型: E不具有傳染性 clc;clear N = 1000; % 總人數 i0 = 1; % 初始時刻患者(已感染者)的人數 [t,x]=ode45('fun1',[1:500],[N-i0 0 i0 0]); % x = round(x); % 對x進行四捨五入(人數為整數) figure(1) % x的第一列是易感染者S的數量,x的第二列是潛伏者E的數量 % x的第三列是患者I的數量, x的第四列是康復者R的數量 plot(t,x(:,1),'r-',t,x(:,2),'m-',t,x(:,3),'b-',t,x(:,4),'g-','Linewidth',1.5) legend('易感染者S','潛伏者E','患者I','康復者R')
結果1
程式碼2
fun2.m
function dx=fun2(t,x) % 大家可以修改裡面的引數,來看結果的變化 beta1 = 0.1; % 易感染者與已感染者接觸且被傳染的強度 beta2 = 0.05; % 易感染者與潛伏者接觸且被傳染的強度 sigma = 0.2; % 潛伏者轉換為感染者的速率 gamma = 0.02; % 康復率 dx = zeros(4,1); % x(1)表示S x(2)表示E x(3)表示I x(4)表示R C = x(1)+x(2)+x(3); % 傳染病系統中的有效人群,也就是課件中的N' = S+E+I dx(1) = - beta1*x(1)*x(3)/C - beta2*x(1)*x(2)/C; dx(2) = beta1*x(1)*x(3)/C + beta2*x(1)*x(2)/C - sigma*x(2); dx(3) = sigma*x(2) - gamma*x(3); dx(4) = gamma*x(3); end
code.m
%% E具有傳染性
clc;clear
N = 1000; % 總人數
i0 = 1; % 初始時刻患者(已感染者)的人數
[t,x]=ode45('fun2',[1:500],[N-i0 0 i0 0]);
% x = round(x); % 對x進行四捨五入(人數為整數)
figure(2)
% x的第一列是易感染者S的數量,x的第二列是潛伏者E的數量
% x的第三列是患者I的數量, x的第四列是康復者R的數量
plot(t,x(:,1),'r-',t,x(:,2),'m-',t,x(:,3),'b-',t,x(:,4),'g-','Linewidth',1.5)
legend('易感染者S','潛伏者E','患者I','康復者R')
% 把兩個圖放到一起看看區別
figure(3)
[t,x1]=ode45('fun1',[1:500],[N-i0 0 i0 0]);
plot(t,x1(:,1),'r-',t,x1(:,2),'m-',t,x1(:,3),'b-',t,x1(:,4),'g-','Linewidth',1.5)
hold on
[t,x2]=ode45('fun2',[1:500],[N-i0 0 i0 0]);
plot(t,x2(:,1),'r--',t,x2(:,2),'m--',t,x2(:,3),'b--',t,x2(:,4),'g--','Linewidth',1.5)
legend('易感染者S','潛伏者E','患者I','康復者R')