数模-微分方程（SI模型及其四种拓展）
作者：互联网
SI模型
代码
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')
结果
SI模型扩展1
代码
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期
考虑1000期
与fun1的基本模型放在一起
SI模型扩展2
代码
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')
结果
SI模型扩展3
代码
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')
结果
SI模型扩展4
代码
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')
结果
标签：i0,感染者,beta,SI,数模,dx,微分方程,TOTAL,人数 来源： https://www.cnblogs.com/jgg54335/p/15185111.html