利用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小时内删除。

作 者 | 郭志龙

编 辑 | 郭志龙

校 对 | 郭志龙