ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

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

2022-05-08 10:33:13  阅读:286  来源: 互联网

标签:i0 感染者 beta SI 数模 dx 微分方程 TOTAL 人数


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

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有