博雅课微分方程模型

动态模型

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

分析对象特征的变化规律

预报对象特征的未来性态

研究控制对象特征的手段

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

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

按照内在规律或用类比法建立微分方程

传染病模型问题

描述传染病的传播过程

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

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

预防传染病蔓延的手段

按照传播过程的一般规律，用机理分析方法建立模型

已感染人数(病人)i(t)

每个病人每天有效接触(足以使人致病)人数为

模型1假设若有效接触的是病人，则不能使病人数增加必须区分已感染者(病人)和未感染者(健康人)建模?模型2区分已感染者(病人)和未感染者(健康人)假设1）总人数N不变，病人和健康人的比例分别为2）每个病人每天有效接触人数为

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

~日接触率SI模型模型21/2tmii010ttm~传染病高潮到来时刻

(日接触率)tm

Logistic模型病人可以治愈！?t=tm,di/dt最大dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始条件’,‘自变量’)

记号:在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分.任何D后所跟的字母为因变量，自变量可以指定或由系统规则选定为确省.

Matlab求微分方程的解析解Matlab命令dsolve('Dy=lam*y*(1-y)','y(0)=i0‘)结果：1/(1-exp(-lam*t)*(-1+i0)/i0)例一例二输入命令:y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')结果为:y=3*exp(-2*x)*sin(5*x)例三输入命令：

[x,y]=dsolve('Dx=2*x-3*y','Dy=4*x-5*y')结果：x=3/4*C1*exp(-2*t)+C2*exp(-t)y=C1*exp(-2*t)+C2*exp(-t)

模型3传染病无免疫性——病人治愈成为健康人，健康人可再次被感染增加假设SIS模型3）病人每天治愈的比例为

~日治愈率建模

~日接触率1/

~感染期

~一个感染期内每个病人的有效接触人数，称为接触数。模型3i0i0接触数

=1~阈值感染期内有效接触感染的健康者人数不超过病人数1-1/

i0模型2(SI模型)如何看作模型3(SIS模型)的特例idi/dt01>10ti>11-1/

i0t

1di/dt<0模型4传染病有免疫性——病人治愈后即移出感染系统，称移出者SIR模型假设1）总人数N不变，病人、健康人和移出者的比例分别为2）病人的日接触率

,日治愈率

,

接触数

=/建模需建立的两个方程模型4SIR模型无法求出的解析解在相平面上研究解的性质模型4消去dtSIR模型相轨线的定义域相轨线11si0D在D内作相轨线的图形，进行分析si101D模型4SIR模型相轨线及其分析传染病蔓延传染病不蔓延s(t)单调减

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

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

i(t)单调降至01/~阈值P3P4P2S0模型4SIR模型预防传染病蔓延的手段

(日接触率)卫生水平

(日治愈率)

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

的估计

降低s0提高r0

提高阈值1/

降低

(=

/

)

,

群体免疫模型4SIR模型被传染人数的估计记被传染人数比例x<<s0i0P1

i0

0,s0

1

小,s0

1提高阈值1/

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

=

模型验证1905年到1906年印度孟买瘟疫数据这些数据相当于给出了dr/dt由得将作Taylor展开前三项在初始值r0=0的前提下，图形结果Matlab程序1.、建立M文件ill1.mfunctiony=ill1(b,t)y=b(1)./(cosh(b(2)*t-b(3)).*cosh(b(2)*t-b(3)))2、建立M文件ill1solve.mdata=[70,65,…];t=[4:32];plot(t,data);holdonb0=[800,0.5,3];b=nlinfit(t,data,‘ill1',b0)data2=ill1(b,t);plot(t,data2,'r')方法一t=[4:32];data=[70,65,…]plot(t,data);holdonfun=@(b,t)b(1)./(cosh(b(2)*t-b(3)).*cosh(b(2)*t-b(3)));b0=[200,0.5,3];b=nlinfit(t,data,fun,b0);data2=ill1(b,t);plot(t,data2,'r')

方法二微分方程的数值解Matlab求常微分方程的数值解

[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：分别为设定的相对误差和绝对误差.1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意:2、ts=[0:50];x0=[0.02,0.98];[t,x]=ode45(‘ill2’,ts,x0);plot(t,x(:,1),‘r’,t,x(:,2),‘b’)pause;plot(x(:,2),x(:,1))解

1、建立m-文件ill2.m如下：

functiondx=ill2(t,x)a=1;b=0.3;dx=[a*x