MATLAB求解微分方程组—以一种传染病的动力学模型求解为例

利用ode45（）函数求解微分方程组

一种可自愈的传染病，未患病的易感人群 S 以一定概率感染之后，成为潜伏期感染者 E，之后部分潜伏期感染者发展成为显性感染者 I，另一部分成为隐性感染者 A，显性和隐性感染者在病愈后成为恢复者 R。随疾病的自由传播，不同类型的人群数量可用如下动力学微分方程来描述：

假设当疾病开始传播 D 天后，防疫部门发现疫情，并开始采取限制措施。现 在有两种限制措施：（1）对显性感染者 I 采取严格隔离措施直至痊愈，阻断显性感染者与其他 人群的接触和传播。采取这一措施后人群数量的变化为：

（2）采用持续 Q 日的所有人群居家隔离措施，但隔离期满后直接恢复原有 的自由接触状态，而未对人员进行全面检测以发现感染者。隔离期间的人群数量变化为：

编程求解

主程序

clc;
clear all;
close all;
x0 = [1665;0;2; 0; 0];%初值
tspan = [0 7];%设置变量范围
tspan1 = [1 14];
tspan2 = [2 28];
[t,x] = ode45(@myfun,tspan,x0);
[t1,x1] = ode45(@myfun,tspan1,x0);
[t2,x2] = ode45(@myfun,tspan2,x0);
%画图
figure;
subplot(3,1,1);
plot(t,x(:,1),'black');
title('持续7日所有人隔离政策') 
legend('为患病的感染人群S')
xlabel('t/天数')
subplot(3,1,2);
plot(t1,x1(:,1),'r');
xlabel('t/天数')
title('持续14日所有人隔离政策') 
legend('为患病的感染人群S')
subplot(3,1,3);
plot(t2,x2(:,1),'b');
xlabel('t/天数')
title('持续28日所有人隔离政策') 
legend('为患病的感染人群S')
figure;
subplot(3,1,1);
plot(t,x(:,2),'g');
xlabel('t/天数')
title('持续7日所有人隔离政策') 
legend('潜伏期感染者E')
subplot(3,1,2);
plot(t1,x1(:,2),'r');
title('持续14日所有人隔离政策') 
legend('潜伏期感染者E')
xlabel('t/天数')
subplot(3,1,3);
plot(t2,x2(:,2),'b');
xlabel('t/天数')
title('持续28日所有人隔离政策') 
legend('潜伏期感染者E')
figure;
subplot(3,1,1);
plot(t,x(:,3),'b');
xlabel('t/天数')
title('持续7日所有人隔离政策') 
legend('显性感染者I')
subplot(3,1,2);
plot(t1,x1(:,3),'r');
xlabel('t/天数')
title('持续14日所有人隔离政策') 
legend('显性感染者I')
subplot(3,1,3);
plot(t2,x2(:,3),'g');
xlabel('t/天数')
title('持续28日所有人隔离政策') 
legend('显性感染者I')
figure;
subplot(3,1,1);
plot(t,x(:,4),'y');
xlabel('t/天数')
title('持续7日所有人隔离政策') 
legend('隐性感染者A')
subplot(3,1,2);
plot(t1,x1(:,4),'r');
xlabel('t/天数')
title('持续14日所有人隔离政策') 
legend('隐性感染者A')
subplot(3,1,3);
plot(t2,x2(:,4),'b');
xlabel('t/天数')
title('持续28日所有人隔离政策') 
legend('隐性感染者A')
figure;
subplot(3,1,1);
plot(t,x(:,5),'r');
title('持续7日所有人隔离政策') 
legend('恢复者R')
xlabel('t/天数')
subplot(3,1,2);
plot(t1,x1(:,5),'g');
xlabel('t/天数')
title('持续14日所有人隔离政策') 
legend('恢复者R')
subplot(3,1,3);
plot(t2,x2(:,5),'b');
xlabel('t/天数')
title('持续14日所有人隔离政策') 
legend('恢复者R')




tspan = [0 7];%设置变量范围
tspan1 = [1 14];
tspan2 = [2 28];
[t3,x3] = ode45(@myfun1,tspan,x0);
[t4,x4] = ode45(@myfun1,tspan1,x0);
[t5,x5] = ode45(@myfun1,tspan2,x0);
%画图
figure;
subplot(3,1,1);
plot(t3,x3(:,1),'black');
title('7日显性感染者I痊愈') 
legend('为患病的感染人群S')
xlabel('t/天数')
subplot(3,1,2);
plot(t4,x4(:,1),'r');
xlabel('t/天数')
title('14日显性感染者I痊愈') 
legend('为患病的感染人群S')
subplot(3,1,3);
plot(t5,x5(:,1),'b');
xlabel('t/天数')
title('28日显性感染者I痊愈') 
legend('为患病的感染人群S')
figure;
subplot(3,1,1);
plot(t3,x3(:,2),'g');
xlabel('t/天数')
title('7日显性感染者I痊愈') 
legend('潜伏期感染者E')
subplot(3,1,2);
plot(t4,x4(:,2),'r');
title('14日显性感染者I痊愈') 
legend('潜伏期感染者E')
xlabel('t/天数')
subplot(3,1,3);
plot(t5,x5(:,2),'b');
xlabel('t/天数')
title('28日显性感染者I痊愈') 
legend('潜伏期感染者E')
figure;
subplot(3,1,1);
plot(t3,x3(:,3),'b');
xlabel('t/天数')
title('7日显性感染者I痊愈') 
legend('显性感染者I')
subplot(3,1,2);
plot(t4,x4(:,3),'r');
xlabel('t/天数')
title('14日显性感染者I痊愈') 
legend('显性感染者I')
subplot(3,1,3);
plot(t5,x5(:,3),'g');
xlabel('t/天数')
title('28日显性感染者I痊愈')  
legend('显性感染者I')
figure;
subplot(3,1,1);
plot(t3,x3(:,4),'y');
xlabel('t/天数')
title('7日显性感染者I痊愈') 
legend('隐性感染者A')
subplot(3,1,2);
plot(t4,x4(:,4),'r');
xlabel('t/天数')
title('14日显性感染者I痊愈') 
legend('隐性感染者A')
subplot(3,1,3);
plot(t5,x5(:,4),'b');
xlabel('t/天数')
title('28日显性感染者I痊愈') 
legend('隐性感染者A')
figure;
subplot(3,1,1);
plot(t3,x3(:,5),'r');
title('7日显性感染者I痊愈') 
legend('恢复者R')
xlabel('t/天数')
subplot(3,1,2);
plot(t4,x4(:,5),'g');
xlabel('t/天数')
title('14日显性感染者I痊愈') 
legend('恢复者R')
subplot(3,1,3);
plot(t5,x5(:,5),'b');
xlabel('t/天数')
title('28日显性感染者I痊愈') 
legend('恢复者R')
xlswrite('x.xlsx',x);
xlswrite('x1.xlsx',x1);
xlswrite('x2.xlsx',x2);
xlswrite('x3.xlsx',x3);
xlswrite('x4.xlsx',x4);
xlswrite('x5.xlsx',x5);

微分方程组程序1

function dydt = myfun1(t,y)
%y(1) = S y(2) = E  y(3) = I ,y(4) = A y(5) = R
%对显性感染者I采取隔离措施直至痊愈
%初始化参数
beta = 0.00105;
k = 0.5;
p = 0.14;
omiga = 0.53;
omiga1 = 0.83;
gama = 0.23;
gama1 = 0.24;
%定义函数
dydt = zeros(5,1);%初始化
dydt(1) = -beta*y(1)*k*y(1);
dydt(2) = beta*y(1)*k*y(1)-(1-p)*omiga*y(2)-p*omiga1*y(2);
dydt(3) = (1-p)*omiga*y(2)-gama*y(3);
dydt(4) = p*omiga1*y(2)-gama1*y(4);
dydt(5) = gama *y(3)+gama1*y(4);
end

微分方程组程序2

function dydt = myfun(t,y)
%y(1) = S y(2) = E  y(3) = I ,y(4) = A y(5) = R
%采用持续Q日政策
%初始化参数
beta = 0.00105;
k = 0.5;
p = 0.14;
omiga = 0.53;
omiga1 = 0.83;
gama = 0.23;
gama1 = 0.24;
%定义函数
dydt = zeros(5,1);%初始化
dydt(1) = 0;
dydt(2) = -(1-p)*omiga*y(2)-p*omiga1*y(2);
dydt(3) = (1-p)*omiga*y(2)-gama*y(3);
dydt(4) = p*omiga1*y(2)-gama1*y(4);
dydt(5) = gama *y(3)+gama1*y(4);
end

结果

本文内容来源于网络，仅供参考学习，如内容、图片有任何版权问题，请联系处理，24小时内删除。

作 者 | 郭志龙
编 辑 | 郭志龙
校 对 | 郭志龙