9传染病模型与微分方程数值解课件

北京邮电大学数学系1第五章微分方程模型常微分方程的数值解5.1传染病模型北京邮电大学数学系1第五章微分方程模型常微分方程的数值北京邮电大学数学系2常微分方程的数值解及实验在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。高数中微分方程解法在实际中基本不会直接使用(一)常微分方程数值解因此，研究常微分方程的数值解法十分必要。北京邮电大学数学系2常微分方程的数值解及实验在生产和科研中所北京邮电大学数学系31、用差商代替导数

若步长h较小，则有故有公式：此即欧拉法(向前欧拉法).(二)建立数值解法的一些途径对应有隐式欧拉法北京邮电大学数学系31、用差商代替导数若步长h较小，则有北京邮电大学数学系42、使用数值积分对方程y’=f(x,y),两边由xi到xi+1积分,并利用梯形公式,有实际应用时，与欧拉公式结合使用故有公式梯形方法/*trapezoidformula*/此即改进的欧拉法北京邮电大学数学系42、使用数值积分对方程y’=f(x,y)北京邮电大学数学系5

中点欧拉公式/*midpointformula*/中心差商近似导数x0x2x1北京邮电大学数学系5中点欧拉公式/*midpoint63、使用泰勒公式以此方法为基础，有龙格-库塔法、线性多步法等方法库塔三阶方法四阶龙格-库塔公式63、使用泰勒公式以此方法为基础，有龙格-库塔法、线性多步法北京邮电大学数学系74、数值公式的精度

当一个数值公式的截断误差可表示为o(hk)时(k为正整数，h为步长)，称它是一个k阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。北京邮电大学数学系74、数值公式的精度当一个数北京邮电大学数学系8[t，x]=solver(‘f’,ts,x0,options）ode45ode23ode113ode15sode23s由待解方程写成的m-文件名ts=[t0，tf]，t0、tf为自变量的初值和终值函数初值条件自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset(’reltol’,rt,’abstol’,at),rt,at：分别为设定的相对误差和绝对误差.(三)用Matlab软件求常微分方程的数值解[t,x]=ode23(‘f’,ts,x0)3级2阶龙格-库塔公式[t,x]=ode45(@f,ts,x0)5级4阶龙格-库塔公式北京邮电大学数学系8[t，x]=solver(‘f’,北京邮电大学数学系91、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意:北京邮电大学数学系91、在解n个未知函数的方程组时，x0和x北京邮电大学数学系10设取步长，从到用四阶龙格-库塔方法微分方程求解实例求解初值问题h=0.2;ts=0:h:1;y0=1;[t,x]=ode45('dfun1',ts,y0);[t,x],plot(t,x)functiondx=dfun1(x,y)dx=y-2*x/y;01.00000.20001.18320.40001.34160.60001.48320.80001.61251.00001.7321建立m-文件dfun1.m如下输入命令3、结果如图北京邮电大学数学系10设取步长北京邮电大学数学系11解:

令y1=x，y2=y1’,1、建立m-文件dfun2.m如下：functiondx=dfun2(t,y)dx=[y(2);(1-y(1)^2)*y(2)-y(1)];例则微分方程变为一阶微分方程组：[t,y]=ode45(@dfun2,[0,20],[2,0]);[t,y]plot(t,y(:,1),'r-',t,y(:,2),'b.-');holdonplot(y(:,1),y(:,2),'co');holdofflegend('t~x','t~x`','x~y');2、取t0=0,tf=20,输入命令：北京邮电大学数学系11解:令y1=x，y2=y1’,1、北京邮电大学数学系123、结果如图[t,x]=ode45(@dfun2,[0,20],[2,0])plot(t,x(:,1),'r-',t,x(:,2),'b.-');holdonplot(x(:,1),x(:,2),'co');holdofflegend('t~x','t~x`','x~y');北京邮电大学数学系123、结果如图[t,x]=ode45(@北京邮电大学数学系13

描述对象特征随时间(空间)的演变过程

分析对象特征的变化规律

预报对象特征的未来性态

研究控制对象特征的手段

根据函数及其变化率之间的关系确定函数微分方程建模

根据建模目的和问题分析作出简化假设

按照内在规律或用类比法建立微分方程动态模型北京邮电大学数学系13描述对象特征随时间(空间)的演变过程北京邮电大学数学系145.1传染病模型

描述传染病的传播过程

分析受感染人数的变化规律

预报传染病高潮到来的时刻

预防传染病蔓延的手段

本世纪初,瘟疫常在世界上某地流行,随着人类文明的不断进步,很多疾病,诸如天花、霍乱已经得到有效的控制.然而,即使在今天,一些贫穷的发展中国家，仍出现传染病流行的现象,医疗卫生部门的官员与专家所关注的问题是：问题提出感染疾病的人数与哪些因素有关?北京邮电大学数学系145.1传染病模型描述传染病的传北京邮电大学数学系15问题分析

不同类型传染病的传播过程有不同的特点,故不从医学的角度对各种传染病的传播过程一一进行分析，而是按一般的传播机理建立模型．

由于传染病在传播的过程涉及因素较多，在分析问题的过程中，不可能通过一次假设建立完善的数学模型.思路：针对结果中的不合理之处，逐步修改假设，最终得出较好的模型。先做出最简单的假设，对得出的结果进行分析，北京邮电大学数学系15问题分析不同类型传染病的传北京邮电大学数学系16模型一模型假设：（1）一人得病后，久治不愈，人在传染期内不会死亡。(2)假设每个病人每天有效接触(足以使人致病)人数为

设已感染人数(病人)x(t),假设是连续可微函数建模?北京邮电大学数学系16模型一模型假设：(2)假设每个病人每北京邮电大学数学系17举个实例x=0:0.1:10;y=exp(x);plot(x,y,'b-');最初只有1个病人，1个病人一天可传染1个人exp(10)=22026北京邮电大学数学系17举个实例x=0:0.1:10;最初只有北京邮电大学数学系18被传染的机会也减少，于是

将变小。若有效接触的是病人，则不能使病人数增加必须区分已感染者(病人)和未感染者(健康人)模型缺点问题：随着时间的推移，病人的数目将无限增加，这一点与实际情况不符．模型修改的关键：

的变化规律原因：当不考虑传染病期间的出生、死亡和迁移时，一个地区的总人数可视为常数,在传染病流行初期，

较大，因此

应为时间t的函数。随着病人的增多，健康人数减少,北京邮电大学数学系18被传染的机会也减少，于是将变小。若有北京邮电大学数学系19模型2区分未感染者(健康人)和已感染者(病人)假设1）总人数N不变，健康人和病人的比例分别为2）每个病人每天有效接触人数为

,且使接触的健康人致病建模

~日接触率SI模型SusceptibleInfective北京邮电大学数学系19模型2区分未感染者(健康人)和已感染者北京邮电大学数学系20Logistic模型方程的解：传染病患者比例与时间t关系传染病人数的变化率与患者比率i的关系染病人数由开始到高峰并逐渐达到稳定增长速度由低增至最高后降落下来对模型作进一步分析i~t感染病人占一半时传染率最大！北京邮电大学数学系20Logistic模型方程的解：传染病北京邮电大学数学系21模型21/2tmtm~传染病高潮到来时刻

(日接触率)tm

，推迟传染高峰的到来,即改善保健措施,提高卫生水平可推迟传染病高潮到来.t=tm,(i=1/2),di/dt最大病人最多的一天日接触率

表示该地区的卫生水平，越小卫生水平越高。i~t北京邮电大学数学系21模型21/2tmtm~传染病高潮到来时北京邮电大学数学系22模型的缺点缺点：当t→∞时，i(t)→1，这表示所有的人最终都将成为病人，这一点与实际情况也不符原因：这是由假设(1)所导致，没有考虑病人可以治愈及病人病发身亡的情况。思考题：考虑有病人病发身亡的情况，再对模型进行修改。北京邮电大学数学系22模型的缺点缺点：当t→∞时，i(t)北京邮电大学数学系23传染病无免疫性——病人治愈成为健康人，健康人可再次被感染增加假设SIS模型3）病人平均每天治愈总病人数的比例为

~日治愈率模型3每天治愈的病人为μNi

；病人治愈后成为仍可被感染的健康者。健康者和病人在总人数中所占的比例分别为s(t)、i(t)，则：s(t)+i(t)＝1(1/μ称为传染病的平均传染期)北京邮电大学数学系23传染病无免疫性——病人治愈成为健康人，北京邮电大学数学系24建模

~日接触率1/

~感染期解析法可求解该模型方程的解=北京邮电大学数学系24建模~日接触率1/~感染期解北京邮电大学数学系25模型讨论

~一个感染期内每个病人的有效接触人数，称为接触数1-1/

idi/dt01>1i0i00ti>11-1/

阙值北京邮电大学数学系25模型讨论~一个感染期内每个病人北京邮电大学数学系26i0i00ti>11-1/

是因为随着传染期内被传染人数占当时健康人数的比例的增加，当时的病人数所占比例也随之上升当σ增大时，i(∞)也增大，北京邮电大学数学系26i0i00ti>11-1/是因为北京邮电大学数学系27控制有效接触（隔离的效果）将最终消灭传染病。原因：感染期内有效接触使健康人数变成的病人人数不超过把病人治愈的人数。模型2(SI模型)如何看作模型3(SIS模型)的特例i0i0t

1di/dt<0因此接触数

=1~阈值(没有康复的)隔离北京邮电大学数学系27控制有效接触（隔离的效果）将最终消灭传北京邮电大学数学系28传染病有免疫性——病人治愈后即移出感染系统，称移出者SIR模型假设1）总人数N不变，病人、健康人和移出者的比例分别为2）病人的日接触率

,日治愈率

,

接触数

=/建模需建立的两个方程模型4某些传染病如麻疹等，治愈后均有很强的免疫力,所以病愈的人既非健康人，也非病人。SusceptibleInfectiveRemoved北京邮电大学数学系28传染病有免疫性——病人治愈后即移出感染北京邮电大学数学系29SIR模型无法求出的解析解在相平面上研究解的性质is~北京邮电大学数学系29SIR模型无法求出北京邮电大学数学系30该方程组无法得到解析解，只能用数值计算的方法。在讨论方程组解的性质时，通常需要用到两个概念：所谓解曲线就是方程组的解相轨线就是将时间参数t消去后得到的i与s的关系曲线解曲线和相轨线解曲线和相轨线北京邮电大学数学系30该方程组无法得到解析解，只能用数值计算北京邮电大学数学系31s~i模型数值解functiondx=dfill(t,x)a=1;b=0.3;dx=[a*x(2)*x(1)-b*x(1);-a*x(2)*x(1)];ts=0:50;x0=[0.02,0.98];[t,x]=ode45('dfill',ts,x0);[t,x]plot(t,x(:,1),t,x(:,2)),gridpause,plot(x(:,2),x(:,1)),gridis0.02000.98000.03900.95250.07320.90190.12850.81690.27950.54380.33120.3995

0.34440.28390.32470.2027…

0.00010.03990.00000.0398北京邮电大学数学系31s~i模型数值解functiondx北京邮电大学数学系32消去dtSIR模型相轨线的定义域相轨线11si0D在D内作相轨线的图形，进行分析北京邮电大学数学系32消去dtSIR模型相轨线北京邮电大学数学系33模型结果分析1、不论初始条件如何，病人总要消失首先由模型知其次北京邮电大学数学系33模型结果分析1、不论初始条件北京邮电大学数学系34相轨线及其分析si101DSIR模型传染病有蔓延过程传染病不蔓延s(t)单调减

相轨线的方向P1s0imP1:s0>1/

i(t)先升后降至0P2:s0<1/

i(t)单调降至01/

~是传染病蔓延与否的阈值P3P4P2S0北京邮电大学数学系34相轨线及其分析si北京邮电大学数学系35预防传染病蔓延的手段

(日接触率)卫生水平

(日治愈率)

医疗水平传染病不蔓延的条件——s0<1/

降低s0提高r0

提高阈值1/

降低

(=

/

)

,

群体免疫与预防北京邮电大学数学系35预防传染病蔓延的手段(日接触率)北京邮电大学数学系36

的估计相轨线一次传染病结束后，可估计出

~一个感染期内每个病人的有效接触人数，称为接触数北京邮电大学数学系36的估计相轨线一次传染病结束后，可估北京邮电大学数学系37被传染人数的估计记被传染人数比例x<<s0i0P1

i0

0,s0

1

小,s0

1提高阈值1/

降低被传染人数比例xs0-1/

=

当

<