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

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

SIR模型

程式碼

fun1.m

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

code.m

%% 最簡單的SIR模型
clc;clear
N = 1000;  % 總人數
i0 = 1; % 初始時刻患者（已感染者）的人數
[t,x]=ode45('fun1',[1:500],[N-i0 i0 0]);   % 記得給康復者R一個初始值
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')

結果

分別取康復率為0.01,0.02一直到0.05,把患病人數I的結果放到一個圖形上(相對於上面的模型)

程式碼

fun2.m

function dx=fun2(t,x)   % 大家可以修改裡面的引數，來看結果的變化
    global GAMMA % 定義康復率為全域性變數
    beta = 0.1;  % 易感染者與已感染者接觸且被傳染的強度
    dx = zeros(3,1);  % x(1)表示S  x(2)表示I  x(3)表示R
    C = x(1)+x(2);  % 傳染病系統中的有效人群,也就是課件中的N' = S+I
    dx(1) = - beta*x(1)*x(2)/C;  
    dx(2) = beta*x(1)*x(2)/C - GAMMA*x(2);
    dx(3) = GAMMA*x(2);
end

code.m

%% 分別取康復率為0.01,0.02一直到0.05,把患病人數I的結果放到一個圖形上
clc;clear
global GAMMA % 定義康復率為全域性變數(我習慣用大寫字母來定義)
N = 1000;  % 總人數
i0 = 1; % 初始時刻患者（已感染者）的人數
% https://ww2.mathworks.cn/help/matlab/ref/colormap.html
c = jet(5);  % jet是matlab自帶的顏色圖陣列，等會畫圖我們要用不同顏色區分
% 每行一個0-1之間三元素 RGB 向量。第一列紅色強度。第二列綠色強度。第三列藍色強度。
for i = 1:5
    GAMMA = 0.01*i;   % 在迴圈中來改變康復率GAMMA的值
    [t,x]=ode45('fun2',[1:500],[N-i0 i0 0]);   % 記得給康復者R一個初始值
    x = round(x);  % 對x進行四捨五入(人數為整數)
    % x的第一列是易感染者S的數量，x的第二列是患者I的數量，x的第三列是康復者R的數量
    figure(2)
    % 這裡我們就只繪製患病人數I的圖形了
    plot(t, x(:,2), '-', 'color', c(i,:),'DisplayName',num2str(GAMMA),'Linewidth',1.5)    
    % 'color' 後面的c(i,:)是前面生成的一個RGB向量，可以用來指定這個曲線的顏色
    % 'DisplayName' 後面需要加上當前圖例對應的文字
    hold on
end
legend show    % 要加上這句話，否則圖例不會顯示
% 等價於直接使用命令： legend('0.01','0.02','0.03','0.04','0.05')

結果

不定義全域性變數，我們可以用另外一種方法（這種方法的通用性更強）

由於我們引數多加了一個GAMMA，所以我們要改為函式控制代碼的方式

程式碼

newfun2.m

% % 不定義全域性變數的寫法
function dx=newfun2(t,x,GAMMA)   % 把GAMMA當成引數傳入進來
    beta = 0.1;  % 易感染者與已感染者接觸且被傳染的強度
    dx = zeros(3,1);  % x(1)表示S  x(2)表示I  x(3)表示R
    C = x(1)+x(2);  % 傳染病系統中的有效人群,也就是課件中的N' = S+I
    dx(1) = - beta*x(1)*x(2)/C;  
    dx(2) = beta*x(1)*x(2)/C - GAMMA*x(2);
    dx(3) = GAMMA*x(2);
end

code.m

% % 不定義全域性變數，我們可以用另外一種方法（這種方法的通用性更強）
clc;clear
N = 1000;  % 總人數
i0 = 1; % 初始時刻患者（已感染者）的人數
c = jet(5);  
for i = 1:5
    GAMMA = 0.01*i;   % 在迴圈中來改變康復率GAMMA的值
    [t,x]=ode45(@(t,x) newfun2(t,x,GAMMA),[1:500],[N-i0 i0 0]);  
    x = round(x);  % 對x進行四捨五入(人數為整數)
    % x的第一列是易感染者S的數量，x的第二列是患者I的數量，x的第三列是康復者R的數量
    figure(2)
    % 這裡我們就只繪製患病人數I的圖形了
    plot(t, x(:,2), '-', 'color', c(i,:),'DisplayName',num2str(GAMMA),'Linewidth',1.5) 
    hold on
end
legend show    % 要加上這句話，否則圖例不會顯示

結果

SIR模型的擴充套件1

程式碼

fun3.m

function dx=fun3(t,x)   % 大家可以修改裡面的引數，來看結果的變化
    beta = 0.1;  % 易感染者與已感染者接觸且被傳染的強度
    gamma = 0.01;  % 康復率
    if t > 100
        gamma = gamma * 10;  % 第100期後研發了疫苗，使得康復率增加為原來的10倍
    end
    dx = zeros(3,1);  % x(1)表示S  x(2)表示I  x(3)表示R
    C = x(1)+x(2);  % 傳染病系統中的有效人群,也就是課件中的N' = S+I
    dx(1) = - beta*x(1)*x(2)/C;  
    dx(2) = beta*x(1)*x(2)/C - gamma*x(2);
    dx(3) = gamma*x(2); 
end

code.m

%% 考慮某種使得康復率gamma增加的因素（例如研發了疫苗、醫療裝置升級等）
% 第100期後研發了疫苗，使得康復率增加為原來的10倍
clc;clear
N = 1000;  % 總人數
i0 = 1; % 初始時刻患者（已感染者）的人數
[t,x]=ode45('fun3',[1:200],[N-i0 i0 0]);   % 記得給康復者R一個初始值
x = round(x);  % 對x進行四捨五入(人數為整數)
figure(3)
% 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')

結果

SIR模型的擴充套件2

程式碼

fun4.m

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

code.m

%% 考慮疾病的死亡率的SIR模型
clc;clear
N = 1000;  % 總人數
i0 = 1; % 初始時刻患者（已感染者）的人數
[t,x]=ode45('fun4',[1:400],[N-i0 i0 0 0]);   % 記得給患病死亡人數ID一個初始值
x = round(x);  % 對x進行四捨五入(人數為整數)
figure(4)
% x的第一列是易感染者S的數量，x的第二列是患者I的數量
% x的第三列是康復者R的數量, x的第四列是患病死亡人數ID的數量
plot(t,x(:,1),'r-',t,x(:,2),'b-',t,x(:,3),'g-',t,x(:,4),'k-','Linewidth',1.5)   
legend('易感染者S','患者I','康復者R','患病死亡人數ID')