• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

SIR模型实现(matlab)

原作者: [db:作者] 来自: [db:来源] 收藏 邀请
 
 
 

 

clc
clear
close all;
A = 0.4;
B = 0.1;
I = 0.4;
S = 0.5;

%ode
tspan = [0 50];
y0 = [I S];
[t, y] = ode45(@(t,y)odefun(t,y,A,B), tspan, y0);
r = 1-y(:,1)-y(:,2);

%euler
n = size(r,1);
h = 50 / (n-1);
t_0 = [0:h:50]';
y_i = zeros(n,1);
y_s = zeros(n,1);
y_i(1) = I;
y_s(1) = S;
for i = 1:n-1
y_i(i+1) = h*[A*y_i(i)*y_s(i) - B*y_i(i)]+y_i(i);
y_s(i+1) = h*[-A*y_i(i)*y_s(i)]+y_s(i);
end
r_0 = 1 - y_i(:,1) - y_s(:,1);

%画图
subplot(2,2,1);
plot(t,y(:,1),'-o',t,y(:,2),'-.',t,r,'g');
hold on;
legend('生病人数:i(t)','健康人数:s(t)','移除人数:r(t)','Location','Best');
ylabel('占人口比例%');
xlabel('时间t');
str = ['接触数λ/μ:',num2str(A/B),' 初始生病人数:',num2str(I),',初始健康人数:',num2str(S)];
text(15,0.4,str,'FontSize',10);
title('SIR模型(ode)');


subplot(2,2,2);
plot(t_0,y_i,'-o',t_0,y_s,'-.',t_0,r_0,'g');
hold on;
legend('生病人数:i(t)','健康人数:s(t)','移除人数:r(t)','Location','Best');
ylabel('占人口比例%');
xlabel('时间t');
str = ['接触数λ/μ:',num2str(A/B),' 初始生病人数:',num2str(I),',初始健康人数:',num2str(S)];
text(15,0.4,str,'FontSize',10);
title('SIR模型(euler)');

subplot(2,2,3);
plot(t_0,y_i,'r-',t,y(:,1),'-.');
diff = sum(abs(y_i - y(:,1)));
str1 = ['生病人数对比图i(t), 误差:',num2str(diff)];
title(str1);
legend('euler','ode','Location','Best');
ylabel('占人口比例%');
xlabel('时间t');

subplot(2,2,4);
plot(t_0,y_s,'r-',t,y(:,2),'-.');
diff = sum(abs(y_s - y(:,2)));
str1 = ['健康人数对比图s(t), 误差:',num2str(diff)];
title(str1);
legend('euler','ode','Location','Best');
ylabel('占人口比例%');
xlabel('时间t');

 


function dydt = odefun(t,y,A,B)
dydt = zeros(2,1);
dydt(1) = A*y(1)*y(2) - B*y(1);
dydt(2) = -A*y(1)*y(2);
end

结果

 

 

 


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
Delphi 7~XE系列升级安装Indy10.6发布时间:2022-07-18
下一篇:
delphi 利用 InterlockedCompareExchange 实现主线程维一锁等待发布时间:2022-07-18
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap