标签：SIR i0 感染者 beta 数模 gamma dx 微分方程 GAMMA
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')
结果
标签：SIR,i0,感染者,beta,数模,gamma,dx,微分方程,GAMMA 来源： https://www.cnblogs.com/jgg54335/p/15185565.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。