下载文档原格式
河北大学《数学模型》实验实验报告一、实验目的二、实验要求1.实验7-1 传染病模型2（ SI模型）——画di/dt~ i曲线图（参考教材 p137-138）传染病模型 2（ SI 模型）：;di/dt=ki(1-i),i(0)=i其中， i(t)是第 t 天病人在总人数中所占的比例。
λ是每个病人每天有效接触的平均人数（日接触率）。
i0是初始时刻（ t=0）病人的比例。
取 k=0.1，画出 di/dt~ i 曲线图，求 i 为何值时di/dt达到最大值，并在曲线图上标注。
试编写一个 m 文件来实现。
参考程序运行结果（在图形窗口菜单选择 Edit/Copy Figure，复制图形）：[提示]1）画曲线图用 fplot 函数，调用格式如下：fplot(fun,lims)fun 必须为一个 M 文件的函数名或对变量 x 的可执行字符串。
若 lims 取[xmin xmax]，则 x 轴被限制在此区间上。
若 lims 取[xmin xmax ymin ymax]，则 y 轴也被限制。
本题可用fplot('0.1*x*(1-x)',[0 1.1 0 0.03]);2）求最大值用求解边界约束条件下的非线性最小化函数 fminbnd，调用格式如下：x=fminbnd(‘fun’,x1,x2)fun 必须为一个 M 文件的函数名或对变量 x 的可执行字符串。
返回自变量 x 在区间 x1<x<x2 上函数取最小值时的 x 值。
本题可用x=fminbnd('-0.1*x*(1-x)',0,1)y=0.1*x*(1-x)4）指示最大值坐标用线性绘图函数plot，调用格式如下：plot(x1,y1,’颜色线型数据点图标’, x2,y2,’颜色线型数据点图标’,…) 说明参见《数学实验》 p225本题可用hold on; %在上面的同一张图上画线（同坐标系）plot([0,x],[y,y],':',[x,x],[0,y],':');3）图形的标注使用文本标注函数 text，调用格式如下：格式 1text(x,y,文本标识内容,’HorizontalAlignment’,’字符串 1’)x,y 给定标注文本在图中添加的位置。
3.12传染病模型摘要：本文是一个对传染病的研究问题。
通过把一般把传染病流行范围内的人群分成三类：S类，易感者(Susceptible)，指未得病者，但缺乏免疫能力，与感染者接触后容易受到感染；I类，感病者(Infective)，指染上传染病的人，它可以传播给S类成员；R类，移出者(Removal)，指被隔离或因病愈而具有免疫力的人。
建立数学模型用极限和微积分等数学方法对传染病传播规律进行研究。
关键词：传染病极限和微积分正文1 传染病〔Infectious Diseases〕是由各种病原体引起的能在人与人、动物与动物或人与动物之间相互传播的一类疾病。
病原体中大部分是微生物，小部分为寄生虫，寄生虫引起者又称寄生虫病。
有些传染病，防疫部门必须及时掌握其发病情况，及时采取对策，因此发现后应按规定时间及时向当地防疫部门报告，称为法定传染病。
中国目前的法定传染病有甲、乙、丙3类，共37种医学科学的发展已经能够有效地预防和控制许多传染病，天花在世界范围内被消灭，鼠疫、霍乱等传染病得到控制。
但是仍然有一些传染病暴发或流行，危害人们的健康和生命。
在发展中国家，传染病的流行仍十分严重；即使在发达国家，一些常见的传染病也未绝迹，而新的传染病还会出现，如爱滋病（AIDS）等。
有些传染病传染很快，导致很高的致残率，危害极大，因而对传染病在人群中传染过程的定量研究具有重要的现实意义。
传染病流行过程的研究与其他学科有所不同，不能通过在人群中实验的方式获得科学数据。
事实上，在人群中作传染病实验是极不人道的。
所以有关传染病的数据、资料只能从已有的传染病流行的报告中获取。
这些数据往往不够全面，难以根据这些数据来准确地确定某些参数，只能大概估计其范围。
基于上述原因，利用数学建模与计算机仿真便成为研究传染病流行过程的有效途径之一。
2问题提出上世纪初，瘟疫还经常在世界的某些地区流行，被传染的人数与哪些因素有关？如何预报传染病高潮的到来？为什么同一地区一种传染病每次流行时，被传染的人数大致不变？3 模型分析社会、经济、文化、风俗习惯等因素都会影响传染病的传播，而最直接的因素是：传染者的数量及其在人群中的分布、被传染者的数量、传播形式、传播能力、免疫能力等，在建立模型时不可能考虑所有因素，只能抓住关键的因素，采用合理的假设，进行简化。
传染病模型摘要当今社会，人们开始意识到通过定量地研究传染病的传播规律,建立传染病的传播模型，可以为预测和控制传染病提供可靠、足够的信息。
本文利用微分方程稳定性理论对传统传染病动力学建模方式进行综述,且针对甲流,SARS等新生传染病模型进行建模和分析。
不同类型的传染病的传播过程有其各自不同的特点，我们不是从医学的角度一一分析各种传染病的传播，而是从一般的传播机理分析建立各种模型，如简单模型，SI模型，SIS模型，SIR模型等。
本文中，我们应用传染病动力学模型来描述疾病发展变化的过程和传播规律，运用联立微分方程组体现疫情发展过程中各类人的内在因果联系,并在此基础上建立方程求解算法。
然后，通过借助Matlab程序拟合出与实际较为符合的曲线并进行了疫情预测，评估各种控制措施的效果,从而不断完善文中的模型。
本文由简到难、全面地评价了该模型的合理性与实用性，而后对模型和数据也做了较为扼要的分析，进一步改进了模型的不妥之处。
同时，在对问题进行较为全面评价的基础上又引入更为全面合理的假设，运用双线性函数模型对卫生部的措施进行了评价并给出建议，做好模型的完善与优化工作.关键词：传染病模型,简单模型，SI,SIS，SIR,微分方程，Matlab。
一、问题重述有一种传染病（如SARS、甲型H1N1）正在流行，现在希望建立适当的数学模型,利用已经掌握的一些数据资料对该传染病进行有效地研究,以期对其传播蔓延进行必要的控制,减少人民生命财产的损失。
考虑如下的几个问题,建立适当的数学模型，并进行一定的比较分析和评价展望.1、不考虑环境的限制，设单位时间内感染人数的增长率是常数，建立模型求t时刻的感染人数。
2、假设单位时间内感染人数的增长率是感染人数的线性函数，最大感染时的增长率为零。
建立模型求t时刻的感染人数。
3、假设总人口可分为传染病患者和易感染者，易感染者因与患病者接触而得病,而患病者会因治愈而减少且对该传染病具有很强的免疫功能，建立模型分析t时刻患病者与易感染者的关系,并对传染情况（如流行趋势，是否最终消灭)进行预测.二、问题分析1、这是一个涉及传染病传播情况的实际问题，其中涉及传染病感染人数随时间的变化情况及一些初始资料，可通过建立相应的微分方程模型加以解决.2、问题表述中已给出了各子问题的一些相应的假设。
实验：传染病模型 Si 模型 问题建立基于以下两个假设的模型，求出平衡点，给出参数，图示模型曲线。
一、模型假设（1）不考虑人口的出生、死亡流动等群动力因素。
人口始终保持一个常数N 。
人群分为易感染者和已感染者，t 时刻这两类人在总人口所占比例分别记作i(t)s(t),（2）一个病人一旦与易感者接触就必然具有一定的传染力。
设每个病人每天有效接触的平均的人数为β，当易感染者与病人接触时就会变成病人。
二、建立模型根据假设，每个病人可以使()s t β个健康人染病，因为病人数为()Ni t ，所以每天共有()()s t Ni t β 于是得Nsi dtdiNβ= 又因为1)()(=+t i t s再记初始时刻（0=t ）的病人比例为0i ，则0)0(),1(**i i i i dt di=-=β解得*01()11(1)*ti t e i β-=+-三、求解平衡点 0)1(**=-=i i dtdiβ，得平衡点1,021==i i 设)1(**)(i i i F -=β 易得 0)'1(,0)'0(<-==βF F 故01=i 不稳定，12=i 稳定四、模型求解Xt表示t时刻病人数，x0表示初始病人数，a表示日接触率>> syms Xt x t x0 a>> [Xt]=dsolve('Dx-a*x*(1-x)','x(0)=x0')Xt =1/(1-exp(-a*t)*(-1+x0)/x0) %求出Xt的符号表达设a=5，x0=0.01则>> a=5;x0=0.01;Xt=1/(1-exp(-a*t)*(-1+x0)/x0)Xt =1/(1+99*exp(-5*t)) %求出Xt的解析解>> fplot('1/(1+99*exp(-5*t))',[0,2]) %做出Xt的图像Xt的图像Sis模型问题建立基于以下三个假设的模型，求出平衡点，给出参数，图示模型曲线。
传染病模型模型假设1. 总人数N 不变.人群分为健康者、病人和病愈免疫的移出者( Re m oved)三类,称SIR 模型.三类人在总人数N 中占的比例分别记作s( t), i( t)和r( t).2. 病人的日接触率为λ,日治愈率为μ(与SI 模型相同),传染期接触数为σ= λ/μ.模型构成由假设1 显然有s( t) + i( t) + r( t) = 1 (12)根据条件2 方程(8)仍成立.对于病愈免疫的移出者而言应有Nd rd t= μNi (13)再记初始时刻的健康者和病人的比例分别是s0 ( s0 > 0)和i0 ( i0 > 0)(不妨设移出者的初始值r0 = 0),则由(8),(12),(13)式, SIR 模型的方程可以写作(14)d i/ d t =λsi - μi, i(0) = i0ds /d t = - λsi, s(0) = s0数值计算在方程(14)中设λ= 1,μ= 0. 3, i(0) = 0. 02, s(0) = 0.98.MA TLAB软件编程function y=ill(t,x)a=1;b=0.3y=[a*x(1)*x(2)-b*x(1),-a*x(1)*x(2)]'输入ts=0:50x0=[0.02,0.98][t,x]=ode45('ill',ts,x0);[t,x]结果为ans =0 0.0200 0.98001.0000 0.0390 0.95252.0000 0.0732 0.90193.0000 0.1285 0.81694.0000 0.2033 0.69275.0000 0.2795 0.54386.0000 0.3312 0.39957.0000 0.3444 0.28398.0000 0.3247 0.20279.0000 0.2863 0.149310.0000 0.2418 0.114511.0000 0.1986 0.091712.0000 0.1599 0.076713.0000 0.1272 0.066514.0000 0.1004 0.059315.0000 0.0787 0.054316.0000 0.0614 0.050717.0000 0.0478 0.048018.0000 0.0371 0.046019.0000 0.0287 0.044520.0000 0.0223 0.043421.0000 0.0172 0.042622.0000 0.0133 0.041923.0000 0.0103 0.041524.0000 0.0079 0.041125.0000 0.0061 0.040826.0000 0.0047 0.040627.0000 0.0036 0.040428.0000 0.0028 0.040329.0000 0.0022 0.040230.0000 0.0017 0.040131.0000 0.0013 0.040032.0000 0.0010 0.040033.0000 0.0008 0.040034.0000 0.0006 0.039935.0000 0.0005 0.039936.0000 0.0004 0.039937.0000 0.0003 0.039938.0000 0.0002 0.039939.0000 0.0002 0.039940.0000 0.0001 0.039941.0000 0.0001 0.039942.0000 0.0001 0.039943.0000 0.0001 0.039944.0000 0.0000 0.039845.0000 0.0000 0.039846.0000 0.0000 0.039847.0000 0.0000 0.039848.0000 0.0000 0.039849.0000 0.0000 0.039850.0000 0.0000 0.0398plot(t,x(:,1),t,x(:,2)),grid,pause图1plot(x(:,2),x(:,1)),grid图2结果：输出的简明计算结果列入表1, i( t), s( t)的图形见图1,图2 是i～ s 的图形,称为相轨线,初值i(0) = 0.02, s(0) = 0.98 相当于图2中的P0 点,随着t的增加,( s, i)沿轨线自右向左运动.由表1、图1、图2 可以看出, i( t)由初值增长至约t = 7 时达到最大值,然后减少, t→∞, i→0; s( t)则单调减少, t →∞, s→0.0398.。
河北大学《数学模型》实验实验报告一、实验目的二、实验要求1.实验7-1 传染病模型2（ SI模型）——画di/dt~ i曲线图（参考教材 p137-138）传染病模型 2（ SI 模型）：;di/dt=ki(1-i),i(0)=i其中， i(t)是第 t 天病人在总人数中所占的比例。
λ是每个病人每天有效接触的平均人数（日接触率）。
i0是初始时刻（ t=0）病人的比例。
取 k=0.1，画出 di/dt~ i 曲线图，求 i 为何值时di/dt达到最大值，并在曲线图上标注。
试编写一个 m 文件来实现。
参考程序运行结果（在图形窗口菜单选择 Edit/Copy Figure，复制图形）：[提示]1）画曲线图用 fplot 函数，调用格式如下：fplot(fun,lims)fun 必须为一个 M 文件的函数名或对变量 x 的可执行字符串。
若 lims 取[xmin xmax]，则 x 轴被限制在此区间上。
若 lims 取[xmin xmax ymin ymax]，则 y 轴也被限制。
本题可用fplot('0.1*x*(1-x)',[0 1.1 0 0.03]);2）求最大值用求解边界约束条件下的非线性最小化函数 fminbnd，调用格式如下：x=fminbnd(‘fun’,x1,x2)fun 必须为一个 M 文件的函数名或对变量 x 的可执行字符串。
返回自变量 x 在区间 x1<x<x2 上函数取最小值时的 x 值。
本题可用x=fminbnd('-0.1*x*(1-x)',0,1)y=0.1*x*(1-x)4）指示最大值坐标用线性绘图函数plot，调用格式如下：plot(x1,y1,’颜色线型数据点图标’, x2,y2,’颜色线型数据点图标’,…) 说明参见《数学实验》 p225本题可用hold on; %在上面的同一张图上画线（同坐标系）plot([0,x],[y,y],':',[x,x],[0,y],':');3）图形的标注使用文本标注函数 text，调用格式如下：格式 1text(x,y,文本标识内容,’HorizontalAlignment’,’字符串 1’)x,y 给定标注文本在图中添加的位置。
’HorizontalAlignment’为水平控制属性，控制文本标识起点位于点(x,y)同一水平线上。
’字符串 1’为水平控制属性值，取三个值之一：‘left’，点(x,y)位于文本标识的左边。
‘center’，点(x,y)位于文本标识的中心点。
‘right’，点(x,y)位于文本标识的右边。
格式 2text(x,y, 文本标识内容,’VerticalAlignment’,’字符串 2’)x,y 给定标注文本在图中添加的位置。
’VerticalAlignment’为垂直控制属性，控制文本标识起点位于点(x,y)同一垂直线上。
’字符串 1’为垂直控制属性值，取四个值之一：‘middle’，’top’，’cap’，’baseline’，’bottom’。
（对应位置可在命令窗口应用确定）本题可用text(0,y,'(di/dt)m','VerticalAlignment','bottom');text(x,-0.001,num2str(x),'HorizontalAlignment','center');4）坐标轴标注调用函数 xlabel， ylabel 和 title本题可用title('SI模型di/dt~i曲线');xlabel('i');ylabel('di/dt');2.实验7-2 传染病模型2（ SI模型）——画i~t曲线图（参考教材 p137-138）传染病模型 2（ SI 模型）：;di/dt=ki(1-i),i(0)=i其中，i(t)是第t 天病人在总人数中所占的比例。
k 是每个病人每天有效接触的平均人数（日接触率）。
i0是初始时刻（t=0）病人的比例求出微分方程的解析解i(t)，画出如下所示的i~t 曲线（i(0)=0.15, k=0.2,t=0~30）。
试编写一个 m 文件来实现。
（在图形窗口菜单选择Edit/Copy Figure，复制图形）[提示]1）求解微分方程常微分方程符号解用函数 dsolve,调用格式如下：dsolve(‘equ1’,’equ2’,…,’变量名’)以代表微分方程及初始条件的符号方程为输入参数，多个方程或初始条件可在一个输入变量内联立输入，且以逗号分隔。
默认的独立变量为 t，也可把 t 变为其他的符号变量。
字符 D 代表对独立变量的微分，通常指 d/dt。
本题可用x=dsolve(‘Dx=k*x*(1-x)’,’x(0)=x0’)2) 画出 i~t 曲线（ i(0)=0.15, λ=0.2, t=0~30）用 for 循环，函数 length, eval, plot, axis, title, xlabel, ylabel3.实验7-3 传染病模型3（ SIS模型）——画di/dt~ i曲线图（参考教材 p138-139）已知传染病模型 3（ SIS 模型）：di/dt=-λ i[i-(1-1/σ )],i(0)=i其中，i(t)是第t 天病人在总人数中所占的比例。
λ是每个病人每天有效接触的平均人数（日接触率）。
i0是初始时刻（t=0）病人的比例。
σ是整个传染期内每个病人有效接触的平均人数（接触数）。
取λ=0.1，σ =1.5，画出如下所示的di/dt~ i曲线图。
试编写一个 m 文件来实现。
（在图形窗口菜单选择 Edit/Copy Figure，复制图形）[提示]用fplot函数画出di/dt~ i曲线图；在上图上用plot函数画一条过原点的水平用title, xlabel, ylabel标注。
4.实验7-4 传染病模型3（ SIS模型）——画i~t曲线图（参考教材 p138-139）已知传染病模型 3（ SIS 模型）：di/dt=-λ i[i-(1-1/σ )],i(0)=i其中，i(t)是第t 天病人在总人数中所占的比例。
λ 是每个病人每天有效接触的平均人数（日接触率）。
i0是初始时刻（t=0）病人的比例。
σ 是整个传染期内每个病人有效接触的平均人数（接触数）。
实验要求 :求出微分方程的解析解i(t)。
取λ =0.2, σ =3, t=0~40，画出如下所示的图形。
试编写一个m 文件来实现。
其中蓝色实线为i(0)=0.2 时的i~t 曲线（第1 条）；黑色虚点线为过点（0, 1-1/σ ）的水平线（第 2 条）；红色虚线为i(0)=0.9 时的i~t 曲线（第3 条）。
[提示]图例标注可用legend('i(0)=0.2','1-1/¦σ ','i(0)=0.9');5.实验7-5 传染病模型4（ SIR模型）（参考教材p140-141）SIR 模型的方程 :di/dt=λ si-μ i i(0)=i0ds/dt=-λ si s(0)=s0实验要求：1．设λ =1，μ =0.3，i(0)=0.02，s(0)=0.98。
输入 p139 的程序，并修改程序中的[t,x]，使得输出的数据格式如下（提示：取 4 位小数，使用四舍五入取整函数 round，矩阵剪裁和拼接）：ans =Columns 1 through 60 1 2 3 4 50.02 0.039 0.0732 0.1285 0.2033 0.27950.98 0.9525 0.9019 0.8169 0.6927 0.5438Columns 7 through 126 7 8 9 10 150.3312 0.3444 0.3247 0.2863 0.2418 0.07870.3995 0.2839 0.2027 0.1493 0.1145 0.0543Columns 13 through 1820 25 30 35 40 45三、实验内容1.实验7-1 传染病模型2（ SI模型）——画di/dt~ i曲线图在matlab中建立M文件fun1.m代码如下：function y=fun(x)k=0.1;y=k*x*[1-x];Fun2.m代码如下：function y=fun(x)k=0.1;y=-k*x*[1-x];在命令行输入以下代码：fplot('fun1',[0 1.1 0 0.03]);x=fminbnd('fun2',0,1);y=0.1*x*(1-x);hold on;plot([0,x],[y,y],'-',[x,x],[0,y],'-');text(0,y,'(di/dt)m','VerticalAlignment','bottom');text(x,-0.001,num2str(x),'HorizontalAlignment','center'); title('SI模型di/dt~i曲线');xlabel('i');ylabel('di/dt');hold off2.实验7-2 传染病模型2（ SI模型）——画i~t曲线图在matlab中建立M文件fun22.m代码如下：k=0.2;x0=0.15;x=dsolve('Dx=k*x*(1-x)','x(0)=x0');tt=linspace(0,31,1001);for i=1:1001t=tt(i);xx(i)=eval(x);endplot(tt,xx)axis([0,31,0,1.1]);title('图1 SI模型i~t曲线');xlabel('t(天)');ylabel('i（病人所占比例）');在命令行输入以下代码：fun22;3.实验7-3 传染病模型3（ SIS模型）——画di/dt~ i曲线图在matlab中建立M文件fun3.m代码如下：function y=fun(x)a=0.1;b=1.5;y=-a*x*[x-(1-1/b)];在命令行输入以下代码：fplot('fun3',[0 0.4 -0.0005 0.003]);x=fminbnd('fun3',0,1);title('SIS模型di/dt~i曲线');xlabel('i');ylabel('di/dt');>> hold on>> plot([0,0.4],[0,0])4.实验7-4 传染病模型3（ SIS模型）——画i~t曲线图在matlab中建立M文件fun4.m代码如下：function y=fun(x)x=dsolve('Dx=-0.2*x*(x-(1-1/3))','x(0)=0.2');tt=linspace(0,41,1001);for i=1:1001t=tt(i);xx(i)=eval(x);endplot(tt,xx);hold on;plot([0,40],[1-1/3,1-1/3],'-k');x=dsolve('Dx=-0.2*x*(x-(1-1/3))','x(0)=0.9');tt=linspace(0,41,1001);for i=1:1001t=tt(i);xx(i)=eval(x);endplot(tt,xx,'-r');axis([0,40,0,1]);title('图1 SI模型i~t曲线(λ =0.2, σ =3)');xlabel('t(天)');ylabel('i（病人所占比例）');legend('i(0)=0.2','1-1/σ ','i(0)=0.9');在命令行输入以下代码：fun4;5.实验7-5 传染病模型4（ SIR模型）在matlab中建立M文件fun5.m代码如下：function y=fun(t,x)a=1;b=0.3;y=[a*x(1)*x(2)-b*x(1),-a*x(1)*x(2)]';在命令行输入以下代码：>> ts=0:50;>> x0=[0.02,0.98];>> [t,x]=ode45('fun5',ts,x0);>> plot(t,x(:,1),t,x(:,2)),grid,pause>> plot(x(:,2),x(:,1)),grid,四、实验结果及其分析1.实验7-1 传染病模型2（ SI模型）——画di/dt~ i曲线图分析：,这时病人增加得在最快，可以认为是医院的门诊当i=1/2时di/dt达到最大值(di/dt)m量最大的一天，预示着传染病高潮的到来，是医疗卫生部门关注的时刻。
