其他分享
首页 > 其他分享> > 数模-微分方程(SI模型及其四种拓展)

数模-微分方程(SI模型及其四种拓展)

作者:互联网

SI模型

image
image

代码

fun1.m

function dx=fun1(t,x)   % 大家可以修改里面的参数,来看结果的变化
    global TOTAL_N   % 定义总人数为全局变量
    beta = 0.1;  % 易感染者与已感染者接触且被传染的强度
    dx = zeros(2,1);  % x(1)表示S  x(2)表示I
    dx(1) = - beta*x(1)*x(2)/TOTAL_N;  
    dx(2) = beta*x(1)*x(2)/TOTAL_N;
end

code.m

%% 最简单的SI模型
%% 直接求解析解求不出来
dsolve('Dx1=-0.1*x1*x2/1000','Dx2=0.1*x1*x2/1000','x1(0)=999,x2(0)=1','t')

%% 根据S+I=N的关系,消去一个变量后就可以求出来了
x1 = dsolve('Dx1=-0.1*x1*(1000-x1)/1000','x1(0)=999','t')  % S的数量
x2 = 1000-x1  % I的数量
figure(1)
fplot(x1,[1,200],'r')
hold on
fplot(x2,[1,200],'b')
legend('易感染者S','患者I')

%% 为了和后面更加复杂的模型的求解统一,我们这里就求数值解
clc;clear
global TOTAL_N   % 定义总人数为全局变量(可以在子函数中使用)
TOTAL_N = 1000;  % 总人数
i0 = 1; % 初始时刻患者(已感染者)的人数
[t,x]=ode45('fun1',[1:200],[TOTAL_N-i0 i0]); 
x = round(x);  % 可以对x进行四舍五入(人数为整数)
% % 注意:四舍五入后人口数加起来可能不等于总人数了,但这个误差可以忽略。
figure(1)
 % x的第一列是易感染者S的数量,x的第二列是患者I的数量
plot(t,x(:,1),'r-',t,x(:,2),'b-','Linewidth',1.5)  % 线的宽度设置为1.5
legend('易感染者S','患者I')

结果

image

SI模型扩展1

image

代码

fun2.m

function dx=fun2(t,x)   % 大家可以修改里面的参数,来看结果的变化
    global TOTAL_N   % 定义总人数为全局变量
    beta = 0.1;  % 易感染者与已感染者接触且被传染的强度
    if t > 60
        beta = beta/10; % 第60期后禁止大规模聚会,使得传染强度beta缩小为原来的10倍
    end
%     beta = 0.1 - 0.001*t
%     if beta < 0
%         beta = 0;
%     end
    dx = zeros(2,1);  % x(1)表示S  x(2)表示I
    dx(1) = - beta*x(1)*x(2)/TOTAL_N;  
    dx(2) = beta*x(1)*x(2)/TOTAL_N;
end

code.m

%% 考虑某种使得参数beta降低的因素(例如禁止大规模聚会、采取隔离措施等)
% 第60期后禁止大规模聚会,使得传染强度beta缩小为原来的10倍
clc;clear
global TOTAL_N   % 定义总人数为全局变量(可以在子函数中使用)
TOTAL_N = 1000;  % 总人数
i0 = 1; % 初始时刻患者(已感染者)的人数
[t,x]=ode45('fun2',[1:200],[TOTAL_N-i0 i0]); 
x = round(x);  % 对x进行四舍五入(人数为整数)
figure(2)
plot(t,x(:,1),'r-',t,x(:,2),'b-','Linewidth',1.5)  % x的第一列是易感染者S的数量,x的第二列是患者I的数量
legend('易感染者S','患者I')

% 考虑1000期
[t,x]=ode45('fun2',[1:1000],[TOTAL_N-i0 i0]); 
x = round(x);  % 对x进行四舍五入(人数为整数)
figure(3)
plot(t,x(:,1),'r-',t,x(:,2),'b-','Linewidth',1.5)  % x的第一列是易感染者S的数量,x的第二列是患者I的数量
legend('易感染者S','患者I')

% 两张图画到一起
[t,x1]=ode45('fun1',[1:500],[TOTAL_N-i0 i0]);   % 原来调用fun1,返回x1
[t,x2]=ode45('fun2',[1:500],[TOTAL_N-i0 i0]);   % 禁止大规模聚会后调用fun2,返回x2
figure(4)
plot(t,x1(:,1),'r-',t,x1(:,2),'b-','Linewidth',1.5)  
hold on  % 接着在上面那个图形上面画图
plot(t,x2(:,1),'r--',t,x2(:,2),'b--','Linewidth',1.5)   % 两个小横线--表示绘制的为虚线
axis([0 500 0 1000])  % 设置坐标轴范围,x轴为0-500 y轴为0-1000
legend('原来S','原来I','现在S','现在I')  % 在图中可以手动拖动图例的位置

结果

200期

image

考虑1000期

image

与fun1的基本模型放在一起

image

SI模型扩展2

image

代码

fun3.m

function dx=fun3(t,x)   % 大家可以修改里面的参数,来看结果的变化
    beta = 0.1;  % 易感染者与已感染者接触且被传染的强度
    u = 0.002;  % 人口出生率
    v = 0.001;  % 自然死亡率
    dx = zeros(3,1);  % x(1)表示易感染者S  x(2)表示已感染者I  x(3)表示自然死亡人数ND
    dx(1) = - beta*x(1)*x(2)/(x(1)+x(2)) +u*(x(1)+x(2)) - v*x(1);  
    dx(2) = beta*x(1)*x(2)/(x(1)+x(2)) - v*x(2);
    dx(3) = v*x(1) + v*x(2);
end

code.m

%% 增加人口自然出生率和死亡率,但不考虑疾病的死亡率
clc;clear
TOTAL_N = 1000;  % 总人数
i0 = 1; % 初始时刻患者(已感染者)的人数
[t,x]=ode45('fun3',[1:200],[TOTAL_N-i0 i0 0]);  % 别忘了给ND初始值0
x = round(x);  % 对x进行四舍五入(人数为整数)
figure(5)
% x的第一列是易感染者S的数量,x的第二列是患者I的数量, x的第三列是自然死亡人数ND
plot(t,x(:,1),'r-',t,x(:,2),'b-',t,x(:,3),'k-','Linewidth',1.5)   
legend('易感染者S','患者I','自然死亡人数ND')

结果

image

SI模型扩展3

image

代码

fun4.m

function dx=fun4(t,x)   % 大家可以修改里面的参数,来看结果的变化
    beta = 0.1;  % 易感染者与已感染者接触且被传染的强度
    d = 0.01;  % 得病的死亡率
    dx = zeros(3,1);  % x(1)表示易感染者S  x(2)表示已感染者I  x(3)表示患病死亡人数ID
    dx(1) = - beta*x(1)*x(2)/(x(1)+x(2)) ;  
    dx(2) = beta*x(1)*x(2)/(x(1)+x(2)) -d*x(2);
    dx(3) = d*x(2);
end

code.m

%% 不考虑人口自然出生率和死亡率,只考虑疾病的死亡率
clc;clear
TOTAL_N = 1000;  % 总人数
i0 = 1; % 初始时刻患者(已感染者)的人数
[t,x]=ode45('fun4',[1:500],[TOTAL_N-i0 i0 0]);  % 别忘了给ID初始值0
x = round(x);  % 对x进行四舍五入(人数为整数)
figure(6)
% x的第一列是易感染者S的数量,x的第二列是患者I的数量, x的第三列是患病死亡人数ID
plot(t,x(:,1),'r-',t,x(:,2),'b-',t,x(:,3),'k-','Linewidth',1.5)   
legend('易感染者S','患者I','患病死亡人数ID')

结果

image

SI模型扩展4

image

代码

fun5.m

function dx=fun5(t,x)   % 大家可以修改里面的参数,来看结果的变化
    beta = 0.1;  % 易感染者与已感染者接触且被传染的强度
    d = 0.01;  % 得病的死亡率
    u = 0.002;  % 人口出生率
    v = 0.001;  % 自然死亡率
    % x(1)表示易感染者S  x(2)表示已感染者I  x(3)表示患病死亡人数ID  x(4)表示自然死亡人数ND
    dx = zeros(4,1);  
    dx(1) = - beta*x(1)*x(2)/(x(1)+x(2))+u*(x(1)+x(2)) - v*x(1) ;  
    dx(2) = beta*x(1)*x(2)/(x(1)+x(2)) -d*x(2)- v*x(2);
    dx(3) = d*x(2);
    dx(4) = v*x(1) + v*x(2);
end

code.m

%% 同时考虑人口自然出生率和死亡率和疾病的死亡率
clc;clear
TOTAL_N = 1000;  % 总人数
i0 = 1; % 初始时刻患者(已感染者)的人数
[t,x]=ode45('fun5',[1:500],[TOTAL_N-i0 i0 0 0]);  % 别忘了给ID和ND的初始值都为0
x = round(x);  % 对x进行四舍五入(人数为整数)
figure(7)
% x的第一列是易感染者S的数量,x的第二列是患者I的数量
% x的第三列是患病死亡人数ID  ,x的第四列是自然死亡人数ND
plot(t,x(:,1),'r-',t,x(:,2),'b-',t,x(:,3),'k-',t,x(:,4),'g-','Linewidth',1.5)  
legend('易感染者S','患者I','患病死亡人数ID','自然死亡人数ND')

结果

image

标签:i0,感染者,beta,SI,数模,dx,微分方程,TOTAL,人数
来源: https://www.cnblogs.com/jgg54335/p/15185111.html