數模-微分方程（SIRS模型）

阿新 • • 發佈：2022-05-09

SIRS模型

程式碼1

fun1.m

function dx=fun1(t,x)   % 大家可以修改裡面的引數，來看結果的變化
    global TOTAL_N   % 定義總人數為全域性變數
    beta = 0.1;  % 易感染者與已感染者接觸且被傳染的強度
    gamma = 0.03;  % 康復率
    alpha = 0.01; % 康復者R再次轉變為易感者S的轉移率
    dx = zeros(3,1);  % x(1)表示S  x(2)表示I x(3)表示R
    dx(1) = - beta*x(1)*x(2)/TOTAL_N + alpha*x(3) ;  
    dx(2) = beta*x(1)*x(2)/TOTAL_N - gamma*x(2);
    dx(3) = gamma*x(2) - alpha*x(3) ;  
end

code.m

%% SIRS模型
clc;clear
global TOTAL_N   % 定義總人數為全域性變數（可以在子函式中使用）
TOTAL_N = 1000;  % 總人數
i0 = 1; % 初始時刻患者（已感染者）的人數
[t,x]=ode45('fun1',[1:500],[TOTAL_N-i0 i0 0]); 
x = round(x);  % 對x進行四捨五入(人數為整數)
figure(1)
% x的第一列是易感染者S的數量，x的第二列是患者I的數量 x的第三列是康復者R的數量
plot(t,x(:,1),'r-',t,x(:,2),'b-',t,x(:,3),'g-','Linewidth',1.5)     
legend('易感染者S','患者I','康復者R')

結果1

程式碼2

fun2.m

function dx=fun2(t,x)   % 大家可以修改裡面的引數，來看結果的變化
    beta = 0.1;  % 易感染者與已感染者接觸且被傳染的強度
    gamma = 0.03;  % 康復率
    alpha = 0.01; % 康復者R再次轉變為易感者S的轉移率
    d = 0.005; % 因病的死亡率
    C = x(1)+x(2)+x(3);  % 傳染病系統中的有效人群,也就是課件中的N' = S+I+R
    dx = zeros(4,1);  % x(1)表示S  x(2)表示I x(3)表示R x(4)表示ID
    dx(1) = - beta*x(1)*x(2)/C + alpha*x(3) ;  
    dx(2) = beta*x(1)*x(2)/C - gamma*x(2)-d*x(2);
    dx(3) = gamma*x(2) - alpha*x(3) ;  
    dx(4) = d*x(2);
end

code.m

%%  引入疾病的死亡率d後的SIRS模型
clc;clear
N = 1000;  % 總人數
i0 = 1; % 初始時刻患者（已感染者）的人數
figure(2)
% 第一列是易感染者S的數量，第二列是患者I的數量
% 第三列是康復者R的數量, 第四列是患病死亡人數ID的數量
[t,x1]=ode45('fun2',[1:500],[N-i0 i0 0 0]);    % 別忘了初始化
x1 = round(x1); 
subplot(2,2,1)  % 用subplot函式畫子圖
plot(t,x1(:,1),'r-',t,x1(:,2),'b-',t,x1(:,3),'g-',t,x1(:,4),'k-','Linewidth',1.5)   
title('500期')

[t,x2]=ode45('fun2',[1:2000],[N-i0 i0 0 0]);  
x2 = round(x2);  
subplot(2,2,2)
plot(t,x2(:,1),'r-',t,x2(:,2),'b-',t,x2(:,3),'g-',t,x2(:,4),'k-','Linewidth',1.5)    
title('2000期')

[t,x3]=ode45('fun2',[1:8000],[N-i0 i0 0 0]);  
x3 = round(x3);  
subplot(2,2,[3,4])
plot(t,x3(:,1),'r-',t,x3(:,2),'b-',t,x3(:,3),'g-',t,x3(:,4),'k-','Linewidth',1.5)    
title('8000期')
legend('易感染者S','患者I','康復者R','患病死亡人數ID')

% 知識點：
% subplot是將多個圖畫到一個平面上的工具。
% 使用方法：subplot（m,n,p）
% 其中，m表示是圖排成m行，n表示圖排成n列
% 也就是整個figure中每一行n個子圖，一共m行
% p表示下面這個圖所在的位置，p=1表示從左到右從上到下的第一個位置。
% p也可以為某一行的連續位置，就像我們上方的[3,4]一樣

結果2

程式碼3

fun3.m

function dx=fun3(t,x)   % 大家可以修改裡面的引數，來看結果的變化
    beta = 0.1;  % 易感染者與已感染者接觸且被傳染的強度
    gamma1 = 0.01;  % 完全康復率
    gamma2 = 0.01;  % 暫時康復率
    alpha = 0.01; % 暫時康復者R2再次轉變為易感者S的轉移率
    d = 0.005; % 因病的死亡率
    C = x(1)+x(2)+x(4);  % 傳染病系統中的有效人群,也就是課件中的N' = S+I+R2
    dx = zeros(5,1);  % x(1)表示S  x(2)表示I x(3)表示R1 x(4)表示R2 x(5)表示ID
    dx(1) = - beta*x(1)*x(2)/C + alpha*x(4) ;  
    dx(2) = beta*x(1)*x(2)/C - gamma1*x(2)-gamma2*x(2)-d*x(2);
    dx(3) = gamma1*x(2) ;  
    dx(4) = gamma2*x(2) - alpha*x(4) ;  
    dx(5) = d*x(2);
end

code.m

%%  引入完全康復者和暫時康復者的模型
clc;clear
N = 1000;  % 總人數
i0 = 1; % 初始時刻患者（已感染者）的人數
figure(3)
% 第一列是易感染者S的數量，第二列是患者I的數量
% 第三列是完全康復者R1的數量, 第四列是暫時康復者R1的數量
% 第五列是患病死亡人數ID的數量
[t,x]=ode45('fun3',[1:1000],[N-i0 i0 0 0 0]);   % 別忘了初始化
x = round(x); 
plot(t,x(:,1),'r-',t,x(:,2),'b-',t,x(:,3),'g-',t,x(:,4),'m-',t,x(:,5),'k-','Linewidth',1.5)     
legend('易感染者S','患者I','完全康復者R1','暫時康復者R2','患病死亡人數ID')