常微分方程模型演示文稿

常微分方程模型演示文稿當前第1頁\共有73頁\編于星期二\22點（優(yōu)選）第九講常微分方程模型當前第2頁\共有73頁\編于星期二\22點Malthus模型預測美國人口當前第3頁\共有73頁\編于星期二\22點Malthus模型預測美國人口當前第4頁\共有73頁\編于星期二\22點Malthus模型預測的優(yōu)缺點優(yōu)點短期預報比較準確缺點不適合中長期預報原因預報時假設人口增長率r

為常數(shù)。沒有考慮環(huán)境對人口增長的制約作用。當前第5頁\共有73頁\編于星期二\22點2.阻滯增長模型假設人口增長率r(t)是t時刻人口x(t)的減函數(shù)：其中，xm

為考慮到受自然資源和環(huán)境條件限制所能容納的最大人口數(shù)量（稱最大人口容量）

模型假設模型建立當前第6頁\共有73頁\編于星期二\22點模型分析（定性分析）人口將遞減并趨向于xm！人口將始終保持xm不變！人口將遞增并趨向于xm！

無論在哪種情況下，人口最終將趨向于最大人口容量！模型求解當前第7頁\共有73頁\編于星期二\22點人口增長率達到最大值當前第8頁\共有73頁\編于星期二\22點阻滯增長模型預測美國人口當前第9頁\共有73頁\編于星期二\22點阻滯增長模型預測美國人口當前第10頁\共有73頁\編于星期二\22點阻滯增長模型預測的優(yōu)缺點優(yōu)點中期預報比較準確缺點理論上很好，實用性不強原因預報時假設固有人口增長率r

以及最大人口容量xm為定值。實際上這兩個參數(shù)（特別是xm

）很難確定，而且會隨著社會發(fā)展情況變化而變化。前面圖中曲線末端分叉就是由于這個原因。當前第11頁\共有73頁\編于星期二\22點6.2

藥物在體內(nèi)的分布與排除藥物進入機體形成血藥濃度(單位體積血液的藥物量)血藥濃度需保持在一定范圍內(nèi)——給藥方案設計

藥物在體內(nèi)吸收、分布和排除過程——藥物動力學

建立房室模型——藥物動力學的基本步驟房室——機體的一部分，藥物在一個房室內(nèi)均勻分布(血藥濃度為常數(shù))，在房室間按一定規(guī)律轉移本節(jié)討論二室模型——中心室(心、肺、腎等)和周邊室(四肢、肌肉等)當前第12頁\共有73頁\編于星期二\22點中心室周邊室給藥排出模型假設中心室(1)和周邊室(2),容積不變藥物在房室間轉移速率及向體外排除

速率，與該室血藥濃度成正比藥物從體外進入中心室，在二室間

相互轉移,從中心室排出體外模型建立當前第13頁\共有73頁\編于星期二\22點線性常系數(shù)非齊次方程對應齊次方程通解模型建立當前第14頁\共有73頁\編于星期二\22點幾種常見的給藥方式1.快速靜脈注射t=0瞬時注射劑量d的藥物進入中心室,血藥濃度立即為d/V1給藥速率f

(t)和初始條件當前第15頁\共有73頁\編于星期二\22點2.恒速靜脈滴注t>T時,c1(t)和c2(t)按指數(shù)規(guī)律衰減趨于零藥物以恒定速率k進入中心室0Tt￡￡當前第16頁\共有73頁\編于星期二\22點吸收室中心室3.口服或肌肉注射相當于藥物(劑量d)先進入吸收室，吸收后再進入中心室吸收室藥量x0(t)當前第17頁\共有73頁\編于星期二\22點參數(shù)估計各種給藥方式下的c1(t),c2(t)

取決于參數(shù)k12,k21,k13,V1,V2以快速靜脈注射為例,在ti(i=1,2,n)測得c1(ti)由較大的用最小二乘法定A,由較小的用最小二乘法定B,當前第18頁\共有73頁\編于星期二\22點參數(shù)估計法一進入中心室的藥物全部排除當前第19頁\共有73頁\編于星期二\22點參數(shù)估計法二%構造非線性擬合函數(shù)[TWOEXPS.M]functionE=twoexps(a,x,y)x=x(:);y=y(:);Y=a(1)*exp(-a(3)*x)+a(2)*exp(-a(4)*x);E=sum((y-Y).^2)當前第20頁\共有73頁\編于星期二\22點a0=[10011];options=optimset('fminsearch');options.TolX=0.01;options.Display='off';a=fminsearch(@ps,a0,options,x,y)

a=[112.23780.18232.1773]當前第21頁\共有73頁\編于星期二\22點6.3傳染病模型問題描述傳染病的傳播過程分析受感染人數(shù)的變化規(guī)律預報傳染病高潮到來的時刻預防傳染病蔓延的手段按照傳播過程的一般規(guī)律，用機理分析方法建立模型當前第22頁\共有73頁\編于星期二\22點已感染人數(shù)(病人)i(t)每個病人每天有效接觸(足以使人致病)人數(shù)為模型1假設若有效接觸的是病人，則不能使病人數(shù)增加必須區(qū)分已感染者(病人)和未感染者(健康人)建模?當前第23頁\共有73頁\編于星期二\22點模型2區(qū)分已感染者(病人)和未感染者(健康人)假設1）總人數(shù)N不變，病人和健康人的比例分別為

2）每個病人每天有效接觸人數(shù)為,且使接觸的健康人致病建模~日接觸率SI模型當前第24頁\共有73頁\編于星期二\22點模型21/2tmii010ttm~傳染病高潮到來時刻

(日接觸率)tmLogistic模型病人可以治愈！?t=tm,di/dt最大當前第25頁\共有73頁\編于星期二\22點模型3傳染病無免疫性——病人治愈成為健康人，健康人可再次被感染增加假設SIS模型3）病人每天治愈的比例為~日治愈率建模~日接觸率1/~感染期

~一個感染期內(nèi)每個病人的有效接觸人數(shù)，稱為接觸數(shù)。當前第26頁\共有73頁\編于星期二\22點模型3i0i0接觸數(shù)=1~閾值感染期內(nèi)有效接觸感染的健康者人數(shù)不超過病人數(shù)1-1/i0模型2(SI模型)如何看作模型3(SIS模型)的特例idi/dt01>10ti>11-1/i0t1di/dt<0當前第27頁\共有73頁\編于星期二\22點模型4傳染病有免疫性——病人治愈后即移出感染系統(tǒng)，稱移出者SIR模型假設1）總人數(shù)N不變，病人、健康人和移出者的比例分別為2）病人的日接觸率

,日治愈率,

接觸數(shù)=/建模需建立的兩個方程當前第28頁\共有73頁\編于星期二\22點模型4SIR模型無法求出的解析解在相平面上研究解的性質當前第29頁\共有73頁\編于星期二\22點模型4消去dtSIR模型相軌線的定義域相軌線11si0D在D內(nèi)作相軌線的圖形，進行分析當前第30頁\共有73頁\編于星期二\22點si101D模型4SIR模型相軌線及其分析傳染病蔓延傳染病不蔓延s(t)單調減相軌線的方向P1s0imP1:s0>1/

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

i(t)單調降至01/

~閾值P3P4P2S0當前第31頁\共有73頁\編于星期二\22點模型4SIR模型預防傳染病蔓延的手段

(日接觸率)衛(wèi)生水平(日治愈率)醫(yī)療水平傳染病不蔓延的條件——s0<1/的估計降低s0提高r0提高閾值1/降低(=/),群體免疫當前第32頁\共有73頁\編于星期二\22點SIR模型被傳染人數(shù)的估計法一記被傳染人數(shù)比例x<<s0i0P1i00,s01小,s01提高閾值

1/

降低被傳染人數(shù)比例xs0-1/=當前第33頁\共有73頁\編于星期二\22點被傳染人數(shù)的估計法二X=fzero(‘x-1.2*log(x/0.96)-0.99’,0.5)

X=0.8651當前第34頁\共有73頁\編于星期二\22點6.4多種群生態(tài)數(shù)學模型意大利生物學家D’Ancona曾致力于魚類種群相互制約關系的研究，他從第一次世界大戰(zhàn)期間,地中海各港口捕獲的幾種魚類捕獲量百分比的資料中，發(fā)現(xiàn)鯊魚等的比例有明顯增加（見下表），而供其捕食的食用魚的百分比卻明顯下降.顯然戰(zhàn)爭使捕魚量下降，食用魚增加，鯊魚等也隨之增加，但為何鯊魚的比例大幅增加呢？他無法解釋這個現(xiàn)象，于是求助于其岳父，著名的意大利數(shù)學家V.Volterra，希望建立一個食餌—捕食系統(tǒng)的數(shù)學模型，定性或定量地回答這個問題.當前第35頁\共有73頁\編于星期二\22點

該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關系，沒有考慮食餌和捕食者自身的阻滯作用，是Volterra提出的最簡單的模型.模型（一）不考慮捕獲當前第36頁\共有73頁\編于星期二\22點定理

Volterra微分方程組對應初值問題的解是周期函數(shù)，且解的周期平均值為

當前第37頁\共有73頁\編于星期二\22點首先，建立m-文件shier.m如下：

functiondx=shier(t,x)dx=zeros(2,1);dx(1)=x(1)*(1-0.1*x(2));dx(2)=x(2)*(-0.5+0.02*x(1));其次，建立主程序shark.m如下：

[t,x]=ode45('shier',[015],[252]);plot(t,x(:,1),'-',t,x(:,2),'*')plot(x(:,1),x(:,2))ToMatlab(shark)當前第38頁\共有73頁\編于星期二\22點求解結果：左圖反映了x1（t）與x2（t）的關系?？梢圆聹y：x1（t）與x2（t）都是周期函數(shù)。當前第39頁\共有73頁\編于星期二\22點模型（二）考慮人工捕獲設表示捕獲能力的系數(shù)為e，相當于食餌的自然增長率由a降為a-e，捕食者的死亡率由c增為c+e設戰(zhàn)前捕獲能力系數(shù)e=0.3,戰(zhàn)爭中降為e=0.1,則戰(zhàn)前與戰(zhàn)爭中的模型分別為:Volterra原理當前第40頁\共有73頁\編于星期二\22點模型求解:1、分別用m-文件shier1.m和shier2.m定義上述兩個方程2、建立主程序shark1.m,求解兩個方程，并畫出兩種情況下鯊魚數(shù)在魚類總數(shù)中所占比例x2(t)/[x1(t)+x2(t)]ToMatlab(shark1)實線為戰(zhàn)前的鯊魚比例，“*”線為戰(zhàn)爭中的鯊魚比例結論：戰(zhàn)爭中鯊魚的比例比戰(zhàn)前高！當前第41頁\共有73頁\編于星期二\22點functiony=shier(t,x)r=1;d=0.5;a=0.1;b=0.02;y=diag([r-a*x(2),-d+b*x(1)])*x;shier.mts=0::0.1:35;x0=[25,2];[t,x]=ode45('shier',ts,x0);[t,x],plot(t,x),grid,gtext('x1(t)'),gtext('x2(t)'),pause,plot(x(:,1),x(:,2)),grid,xlabel('x1'),ylabel('x2')shiyan42注：ts中終值(=15)和步長=(0.1)的確定當前第42頁\共有73頁\編于星期二\22點當前第43頁\共有73頁\編于星期二\22點計算結果（數(shù)值，圖形）x(t),y(t)是周期函數(shù)，相圖(x,y)是封閉曲線；觀察，猜測x(t),y(t)的周期約為10.7;xmax=99.3,xmin=2.0,ymax=28.4,ymin=2.0.用數(shù)值積分可算出x(t)一周期的平均值為25,y(t)一周期的平均值為10.當前第44頁\共有73頁\編于星期二\22點6.5其它生態(tài)數(shù)學模型存在一大類生態(tài)模型源于對Volterra模型的改造模型1考慮食餌種群與外界有遷入或遷出（1978）外界有食餌遷入外界有食餌遷出也可以表示人工干預，如投放或捕獲模型討論食餌－捕食者模型當前第45頁\共有73頁\編于星期二\22點模型2考慮食餌種群內(nèi)部存在生存競爭G.Bojadziev表示當沒有捕食者存在時食餌種群的環(huán)境容納量模型3考慮食餌和捕食者種群內(nèi)部都存在生存競爭張錦炎（1979）表示當沒有捕食者存在時食餌種群的環(huán)境容納量當前第46頁\共有73頁\編于星期二\22點模型4考慮雙方內(nèi)部都存在生存競爭，且捕食者另有食物來源E.C.Pielou平衡點：當前第47頁\共有73頁\編于星期二\22點模型兩種群競爭模型競爭模型競爭排斥原理（CompetitionExclutionLaw）

多個種群依靠同一個生存資源而生活，如果生活在同一個地理空間，

獵取相同食物或營養(yǎng)物。在有限的相同生存資源條件下，如果存在競爭關系，

它們必然相互排斥，展開激烈的生存競爭。結局是競爭力較弱的種群滅絕，

競爭力最強的種群達到其環(huán)境容納量。

當前第48頁\共有73頁\編于星期二\22點模型討論種群y

最終將被滅絕,種群x最終趨于最大容量種群x

最終將被滅絕,種群y最終趨于最大容量存在過正平衡點的一條分界線，將第一象限分成

種群x和種群y

的兩個吸引域。種群x和種群y

最終達到穩(wěn)定的正平衡態(tài)競爭排斥原理

是針對前三種情形得出的結論，

第四種情況極為罕見。當前第49頁\共有73頁\編于星期二\22點%M函數(shù)functiondy=cwf1(t,y)dy=zeros(2,1);dy(1)=0.1*y(1)*(1-0.001*y(1)-0.0008*y(2));dy(2)=0.2*y(2)*(1-0.0012*y(1)-0.001*y(2));當前第50頁\共有73頁\編于星期二\22點

主程序[T,X]=ode45(‘cwf1’,[0200],[200200]);[T,Y]=ode45(‘cwf1’,[0200],[500200]);[T,Z]=ode45(‘cwf1’,[0200],[1200500]);Plot(X(:,1),X(:,2),(Y(:,1),Y(:,2),(Z(:,1),Z(:,2))當前第51頁\共有73頁\編于星期二\22點

當前第52頁\共有73頁\編于星期二\22點模型兩種群互惠模型互惠模型研究多個種群之間相互依賴、共生現(xiàn)象。模型討論模型有三個平衡點，分別為當前第53頁\共有73頁\編于星期二\22點兩種群最終達到穩(wěn)定平衡態(tài)兩種群共生P3為正平衡點P3穩(wěn)定兩種群最終達不到穩(wěn)定平衡態(tài)P3不穩(wěn)定當前第54頁\共有73頁\編于星期二\22點%m程序functiondy=cwf2(t,y)

dy=zeros(2,1);dy(1)=0.1*y(1)*(1-0.001*y(1)+0.0005*y(2));dy(2)=0.2*y(2)*(-1+0.0015*y(1)-0.001*y(2));當前第55頁\共有73頁\編于星期二\22點[T,X]=ode45(‘cwf2’,[0,100],[1600,2800]);[T,Y]=ode45(‘cwf2’,[0,100],[1200,2500]);[T,Z]=ode45(‘cwf2’,[0,100],[2500,2200]);Plot(X(:,1),X(:,2),Y(:,1),Y(:,2),Z(:,1),Z(:,2))Text(2000,2600,2000,2600,'{\sigma_1<1,\sigma_2>1,\sigma_1\sigma_2<1}')當前第56頁\共有73頁\編于星期二\22點當前第57頁\共有73頁\編于星期二\22點6.6常微分方程的數(shù)值解及實驗（一）常微分方程數(shù)值解的定義在生產(chǎn)和科研中所處理的微分方程往往很復雜且大多得不出一般解。而在實際上對初值問題，一般是要求得到解在若干個點上滿足規(guī)定精確度的近似值，或者得到一個滿足精確度要求的便于計算的表達式。因此，研究常微分方程的數(shù)值解法是十分必要的。當前第58頁\共有73頁\編于星期二\22點歐拉公式1、用差商代替導數(shù)若步長h較小，則有故有公式：此即歐拉法（向前歐拉法）。當前第59頁\共有73頁\編于星期二\22點2、使用數(shù)值積分對方程y’=f(x,y),兩邊由xi到xi+1積分，并利用梯形公式，有：實際應用時，與歐拉公式結合使用：此即改進的歐拉法。故有公式：當前第60頁\共有73頁\編于星期二\22點例1求解初值問題當前第61頁\共有73頁\編于星期二\22點當前第62頁\共有73頁\編于星期二\22點龍格—庫特方法

考慮微分中值定理