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);
　　微分方程组程序1function 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
　　微分方程组程序2function 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小时内删除。
　　作 者 | 郭志龙
　　编 辑 | 郭志龙
　　校 对 | 郭志龙