博雅課微分方程模型

動(dòng)態(tài)模型

描述對象特征隨時(shí)間(空間)的演變過程

分析對象特征的變化規(guī)律

預(yù)報(bào)對象特征的未來性態(tài)

研究控制對象特征的手段

根據(jù)函數(shù)及其變化率之間的關(guān)系確定函數(shù)微分方程建模

根據(jù)建模目的和問題分析作出簡化假設(shè)

按照內(nèi)在規(guī)律或用類比法建立微分方程

傳染病模型問題

描述傳染病的傳播過程

分析受感染人數(shù)的變化規(guī)律

預(yù)報(bào)傳染病高潮到來的時(shí)刻

預(yù)防傳染病蔓延的手段

按照傳播過程的一般規(guī)律，用機(jī)理分析方法建立模型

已感染人數(shù)(病人)i(t)

每個(gè)病人每天有效接觸(足以使人致病)人數(shù)為

模型1假設(shè)若有效接觸的是病人，則不能使病人數(shù)增加必須區(qū)分已感染者(病人)和未感染者(健康人)建模?模型2區(qū)分已感染者(病人)和未感染者(健康人)假設(shè)1）總?cè)藬?shù)N不變，病人和健康人的比例分別為2）每個(gè)病人每天有效接觸人數(shù)為

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

~日接觸率SI模型模型21/2tmii010ttm~傳染病高潮到來時(shí)刻

(日接觸率)tm

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

記號(hào):在表達(dá)微分方程時(shí)，用字母D表示求微分，D2、D3等表示求高階微分.任何D后所跟的字母為因變量，自變量可以指定或由系統(tǒng)規(guī)則選定為確省.

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

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

模型3傳染病無免疫性——病人治愈成為健康人，健康人可再次被感染增加假設(shè)SIS模型3）病人每天治愈的比例為

~日治愈率建模

~日接觸率1/

~感染期

~一個(gè)感染期內(nèi)每個(gè)病人的有效接觸人數(shù)，稱為接觸數(shù)。模型3i0i0接觸數(shù)

=1~閾值感染期內(nèi)有效接觸感染的健康者人數(shù)不超過病人數(shù)1-1/

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

i0t

1di/dt<0模型4傳染病有免疫性——病人治愈后即移出感染系統(tǒng)，稱移出者SIR模型假設(shè)1）總?cè)藬?shù)N不變，病人、健康人和移出者的比例分別為2）病人的日接觸率

,日治愈率

,

接觸數(shù)

=/建模需建立的兩個(gè)方程模型4SIR模型無法求出的解析解在相平面上研究解的性質(zhì)模型4消去dtSIR模型相軌線的定義域相軌線11si0D在D內(nèi)作相軌線的圖形，進(jìn)行分析si101D模型4SIR模型相軌線及其分析傳染病蔓延傳染病不蔓延s(t)單調(diào)減

相軌線的方向P1s0imP1:s0>1/

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

i(t)單調(diào)降至01/~閾值P3P4P2S0模型4SIR模型預(yù)防傳染病蔓延的手段

(日接觸率)衛(wèi)生水平

(日治愈率)

醫(yī)療水平傳染病不蔓延的條件——s0<1/

的估計(jì)

降低s0提高r0

提高閾值1/

降低

(=

/

)

,

群體免疫模型4SIR模型被傳染人數(shù)的估計(jì)記被傳染人數(shù)比例x<<s0i0P1

i0

0,s0

1

小,s0

1提高閾值1/

降低被傳染人數(shù)比例xs0-1/

=

模型驗(yàn)證1905年到1906年印度孟買瘟疫數(shù)據(jù)這些數(shù)據(jù)相當(dāng)于給出了dr/dt由得將作Taylor展開前三項(xiàng)在初始值r0=0的前提下，圖形結(jié)果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')

方法二微分方程的數(shù)值解Matlab求常微分方程的數(shù)值解

[t，x]=solver（’f’,ts,x0,options）ode45ode23ode113ode15sode23s由待解方程寫成的m-文件名ts=[t0,tf]，t0,tf為自變量的初值和終值函數(shù)的初值自變量值函數(shù)值用于設(shè)定誤差限(缺省時(shí)設(shè)定相對誤差10-3,絕對誤差10-6),命令為：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分別為設(shè)定的相對誤差和絕對誤差.1、在解n個(gè)未知函數(shù)的方程組時(shí)，x0和x均為n維向量，m-文件中的待解方程組應(yīng)以x的分量形式寫成.2、使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組.注意: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