常微分方程模型

1、6.1 人口增長(zhǎng)模型,英國(guó)人口學(xué)家Malthus (1766-1834,模型假設(shè),人口自然增長(zhǎng)率 r 為常數(shù)即單位時(shí)間內(nèi)人口的增長(zhǎng)量與當(dāng)時(shí)的人口呈正比,模型建立,1.指數(shù)增長(zhǎng)模型,人口以幾何級(jí)數(shù)增加,模型分析,人口將按指數(shù)規(guī)律無(wú)限增長(zhǎng),人口將始終保持不變,人口將按指數(shù)規(guī)律減少直至絕滅,人口倍增時(shí)間,模型求解,Malthus模型預(yù)測(cè)美國(guó)人口,Malthus模型預(yù)測(cè)美國(guó)人口,Malthus模型預(yù)測(cè)的優(yōu)缺點(diǎn),優(yōu)點(diǎn),短期預(yù)報(bào)比較準(zhǔn)確,缺點(diǎn),不適合中長(zhǎng)期預(yù)報(bào),原因,預(yù)報(bào)時(shí)假設(shè)人口增長(zhǎng)率 r 為常數(shù)。沒(méi)有考慮環(huán)境對(duì)人口增長(zhǎng)的制約作用,2.阻滯增長(zhǎng)模型,假設(shè)人口增長(zhǎng)率 r(t) 是 t 時(shí)刻人口 x(t)

2、的減函數(shù),其中，xm 為考慮到受自然資源和環(huán)境條件限制所能容納的最大人口數(shù)量,稱最大人口容量,模型假設(shè),模型建立,模型分析（定性分析,人口將遞減并趨向于xm,人口將始終保持xm不變,人口將遞增并趨向于xm,無(wú)論在哪種情況下，人口最終將趨向于最大人口容量,模型求解,人口增長(zhǎng)率達(dá)到最大值,阻滯增長(zhǎng)模型預(yù)測(cè)美國(guó)人口,阻滯增長(zhǎng)模型預(yù)測(cè)美國(guó)人口,阻滯增長(zhǎng)模型預(yù)測(cè)的優(yōu)缺點(diǎn),優(yōu)點(diǎn),中期預(yù)報(bào)比較準(zhǔn)確,缺點(diǎn),理論上很好，實(shí)用性不強(qiáng),原因,預(yù)報(bào)時(shí)假設(shè)固有人口增長(zhǎng)率 r 以及最大人口容量 xm 為定值,實(shí)際上這兩個(gè)參數(shù)（特別是 xm ）很難確定，而且會(huì)隨著社會(huì)發(fā)展情況變化而變化,前面圖中曲線末端分叉就是由于這個(gè)原因

3、,6.2 藥物在體內(nèi)的分布與排除,藥物進(jìn)入機(jī)體形成血藥濃度(單位體積血液的藥物量,血藥濃度需保持在一定范圍內(nèi)給藥方案設(shè)計(jì),藥物在體內(nèi)吸收、分布和排除過(guò)程 藥物動(dòng)力學(xué),建立房室模型藥物動(dòng)力學(xué)的基本步驟,房室機(jī)體的一部分，藥物在一個(gè)房室內(nèi)均勻分布(血藥濃度為常數(shù))，在房室間按一定規(guī)律轉(zhuǎn)移,本節(jié)討論二室模型中心室(心、肺、腎等)和周邊室(四肢、肌肉等,模型假設(shè),中心室(1)和周邊室(2),容積不變,藥物在房室間轉(zhuǎn)移速率及向體外排除 速率，與該室血藥濃度成正比,藥物從體外進(jìn)入中心室，在二室間 相互轉(zhuǎn)移,從中心室排出體外,模型建立,線性常系數(shù)非齊次方程,對(duì)應(yīng)齊次方程通解,模型建立,幾種常見(jiàn)的給藥方式,1

4、.快速靜脈注射,t=0 瞬時(shí)注射劑量d 的藥物進(jìn)入中心室,血藥濃度立即為 d/V1,給藥速率 f (t) 和初始條件,2.恒速靜脈滴注,t T時(shí), c1(t)和 c2(t)按指數(shù)規(guī)律衰減趨于零,3.口服或肌肉注射,相當(dāng)于藥物( 劑量d)先進(jìn)入吸收室，吸收后再進(jìn)入中心室,吸收室藥量x0(t,參數(shù)估計(jì),各種給藥方式下的 c1(t), c2(t) 取決于參數(shù)k12, k21, k13, V1,V2,以快速靜脈注射為例 ,在ti(i=1,2,n)測(cè)得c1(ti,由較大的 用最小二乘法定A,由較小的 用最小二乘法定B,參數(shù)估計(jì)法一,進(jìn)入中心室的藥物全部排除,參數(shù) 估計(jì)法二, 構(gòu)造非線性擬合函數(shù) TWOE

5、XPS.M function E=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,a0=100 1 1 ; options=optimset(fminsearch); options.TolX=0.01; options.Display=off; a=fminsearch(ps,a0,options,x,y,a= 112.2378 0.1823 2.1773,6.3 傳染病模型,問(wèn)題,描述傳染病的傳播過(guò)程,分析受感染人數(shù)的變化規(guī)律,預(yù)報(bào)傳染病高潮到來(lái)的時(shí)刻,預(yù)防傳染病蔓延的手段,

6、按照傳播過(guò)程的一般規(guī)律，用機(jī)理分析方法建立模型,已感染人數(shù) (病人) i(t,每個(gè)病人每天有效接觸(足以使人致病)人數(shù)為,模型1,假設(shè),若有效接觸的是病人，則不能使病人數(shù)增加,建模,模型2,區(qū)分已感染者(病人)和未感染者(健康人,假設(shè),1）總?cè)藬?shù)N不變，病人和健康 人的 比例分別為,2）每個(gè)病人每天有效接觸人數(shù)為, 且使接觸的健康人致病,建模,日 接觸率,SI 模型,模型2,tm傳染病高潮到來(lái)時(shí)刻,(日接觸率) tm,病人可以治愈,t=tm, di/dt 最大,模型3,傳染病無(wú)免疫性病人治愈成為健康人，健康人可再次被感染,增加假設(shè),SIS 模型,3）病人每天治愈的比例為,日治愈率,建模,日接觸

7、率,1/ 感染期,一個(gè)感染期內(nèi)每個(gè)病人的有效接觸人數(shù)，稱為接觸數(shù),模型3,接觸數(shù) =1 閾值,感染期內(nèi)有效接觸感染的健康者人數(shù)不超過(guò)病人數(shù),模型2(SI模型)如何看作模型3(SIS模型)的特例,模型4,傳染病有免疫性病人治愈后即移出感染系統(tǒng)，稱移出者,SIR模型,假設(shè),1）總?cè)藬?shù)N不變，病人、健康人和移出者的比例分別為,2）病人的日接觸率 , 日治愈率, 接觸數(shù) = /,建模,需建立 的兩個(gè)方程,模型4,SIR模型,模型4,SIR模型,相軌線 的定義域,在D內(nèi)作相軌線 的圖形，進(jìn)行分析,模型4,SIR模型,相軌線 及其分析,s(t)單調(diào)減相軌線的方向,P1: s01/ i(t)先升后降至0,P

8、2: s01/ i(t)單調(diào)降至0,1/ 閾值,模型4,SIR模型,預(yù)防傳染病蔓延的手段,(日接觸率) 衛(wèi)生水平,日治愈率) 醫(yī)療水平,傳染病不蔓延的條件s01,的估計(jì),降低 s0,提高 r0,提高閾值 1,SIR模型,被傳染人數(shù)的估計(jì)法一,記被傳染人數(shù)比例,小, s0 1,提高閾值 1/ 降低被傳染人數(shù)比例 x,s0 - 1/ =,被傳染人數(shù)的估計(jì)法二,X=fzero(x-1.2*log(x/0.96)-0.99,0.5) X=0.8651,6.4多種群生態(tài)數(shù)學(xué)模型,意大利生物學(xué)家DAncona曾致力于魚(yú)類種群相互制約關(guān)系的研究，他從第一次世界大戰(zhàn)期間,地中海各港口捕獲的幾種魚(yú)類捕獲量百分比

9、的資料中，發(fā)現(xiàn)鯊魚(yú)等的比例有明顯增加（見(jiàn)下表），而供其捕食的食用魚(yú)的百分比卻明顯下降.顯然戰(zhàn)爭(zhēng)使捕魚(yú)量下降，食用魚(yú)增加，鯊魚(yú)等也隨之增加，但為何鯊魚(yú)的比例大幅增加呢,他無(wú)法解釋這個(gè)現(xiàn)象，于是求助于其岳父，著名的意大利數(shù)學(xué)家V.Volterra，希望建立一個(gè)食餌捕食系統(tǒng)的數(shù)學(xué)模型，定性或定量地回答這個(gè)問(wèn)題,該模型反映了在沒(méi)有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關(guān)系，沒(méi)有考慮食餌和捕食者自身的阻滯作用，是Volterra提出的最簡(jiǎn)單的模型,模型（一）不考慮捕獲,定理 Volterra微分方程組對(duì)應(yīng)初值問(wèn)題,的解,是周期函數(shù)，且解的周期平均值為,首先，建立m-文件shier.m如下： func

10、tion dx=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,0 15,25 2); plot(t,x(:,1),-,t,x(:,2),*) plot(x(:,1),x(:,2,To Matlab(shark,求解結(jié)果,左圖反映了x1（t）與x2（t）的關(guān)系。 可以猜測(cè)： x1（t）與x2（t）都是周期函數(shù),模型（二） 考慮人工捕獲,設(shè)表示捕獲能力的系數(shù)為e，相當(dāng)于食餌的自然增長(zhǎng)率由a降為a-e，捕食者的死亡率由c增為

11、c+e,設(shè)戰(zhàn)前捕獲能力系數(shù)e=0.3, 戰(zhàn)爭(zhēng)中降為e=0.1, 則戰(zhàn)前與戰(zhàn)爭(zhēng)中的模型分別為,Volterra原理,模型求解,1、分別用m-文件shier1.m和shier2.m定義上述兩個(gè)方程,2、建立主程序shark1.m, 求解兩個(gè)方程，并畫(huà)出兩種情況下鯊魚(yú)數(shù)在魚(yú)類總數(shù)中所占比例 x2(t)/x1(t)+x2(t,To Matlab(shark1,實(shí)線為戰(zhàn)前的鯊魚(yú)比例，“*”線為戰(zhàn)爭(zhēng)中的鯊魚(yú)比例,結(jié)論：戰(zhàn)爭(zhēng)中鯊魚(yú)的比例比戰(zhàn)前高,注：ts 中終值(=15)和步長(zhǎng)=(0.1)的確定,計(jì)算結(jié)果（數(shù)值，圖形,x(t),y(t)是周期函數(shù)，相圖(x,y)是封閉曲線,x(t),y(t)的周期約為10.

12、7,xmax=99.3, xmin=2.0, ymax=28.4, ymin=2.0,用數(shù)值積分可算出x(t)一周期的平均值為25, y(t)一周期的平均值為10,6.5其它生態(tài)數(shù)學(xué)模型,存在一大類生態(tài)模型源于對(duì)Volterra模型的改造,模型1,考慮食餌種群與外界有遷入或遷出 G.R.Walsh（1978,外界有食餌遷入,外界有食餌遷出,也可以表示人工干預(yù)，如投放或捕獲,模型討論,食餌捕食者模型,模型2,考慮食餌種群內(nèi)部存在生存競(jìng)爭(zhēng) G.Bojadziev,表示當(dāng)沒(méi)有捕食者存在時(shí)食餌種群的環(huán)境容納量,模型3,考慮食餌和捕食者種群內(nèi)部都存在生存競(jìng)爭(zhēng) 張錦炎（1979,表示當(dāng)沒(méi)有捕食者存在時(shí)食餌

13、種群的環(huán)境容納量,模型4,考慮雙方內(nèi)部都存在生存競(jìng)爭(zhēng)，且捕食者另有食物來(lái)源 E.C.Pielou,平衡點(diǎn),模型,兩種群競(jìng)爭(zhēng)模型,競(jìng)爭(zhēng)模型,競(jìng)爭(zhēng)排斥原理（Competition Exclution Law,多個(gè)種群依靠同一個(gè)生存資源而生活，如果生活在同一個(gè)地理空間，獵取相同食物或營(yíng)養(yǎng)物。在有限的相同生存資源條件下，如果存在競(jìng)爭(zhēng)關(guān)系，它們必然相互排斥，展開(kāi)激烈的生存競(jìng)爭(zhēng)。結(jié)局是競(jìng)爭(zhēng)力較弱的種群滅絕，競(jìng)爭(zhēng)力最強(qiáng)的種群達(dá)到其環(huán)境容納量,模型討論,競(jìng)爭(zhēng)排斥原理 是針對(duì)前三種情形得出的結(jié)論， 第四種情況極為罕見(jiàn), M 函數(shù) function dy=cwf1 (t,y) dy=zeros(2,1); dy(

14、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,主程序 T,X=ode45(cwf1,0 200,200 200); T,Y=ode45(cwf1,0 200,500 200); T,Z=ode45(cwf1,0 200,1200 500); Plot(X(:,1),X(:,2), (Y(:,1),Y(:,2),(Z(:,1),Z(:,2,模型,兩種群互惠模型,互惠模型,研究多個(gè)種群之間相互依賴、共生現(xiàn)象,模型討論,模型有三個(gè)平衡點(diǎn)，分別為,兩種群最終達(dá)到穩(wěn)定平衡態(tài),兩種群共生,P3

15、為正平衡點(diǎn),P3 穩(wěn)定,兩種群最終達(dá)不到穩(wěn)定平衡態(tài),P3 不穩(wěn)定,m程序,function dy=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,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(:

16、,2) Text(2000,2600, 2000,2600,sigma_11,sigma_1sigma_21,6.6常微分方程的數(shù)值解及實(shí)驗(yàn),一）常微分方程數(shù)值解的定義,在生產(chǎn)和科研中所處理的微分方程往往很復(fù)雜且大多得不出一般解。而在實(shí)際上對(duì)初值問(wèn)題，一般是要求得到解在若干個(gè)點(diǎn)上滿足規(guī)定精確度的近似值，或者得到一個(gè)滿足精確度要求的便于計(jì)算的表達(dá)式,因此，研究常微分方程的數(shù)值解法是十分必要的,歐拉公式,1、用差商代替導(dǎo)數(shù),若步長(zhǎng)h較小，則有,故有公式,此即歐拉法（向前歐拉法,2、使用數(shù)值積分,對(duì)方程y=f(x,y), 兩邊由xi到xi+1積分，并利用梯形公式，有,實(shí)際應(yīng)用時(shí)，與歐拉公式結(jié)合使用,

17、此即改進(jìn)的歐拉法,故有公式,例1 求解初值問(wèn)題,龍格庫(kù)特方法,考慮微分中值定理,2階龍格庫(kù)特公式,利用泰勒展式得到,4階龍格庫(kù)特公式,一階方程組與高階方程,Volterra方程解曲線 與Volterra方程相軌線,3、使用泰勒公式,以此方法為基礎(chǔ)，有龍格-庫(kù)塔法、線性多步法等方法,4、數(shù)值公式的精度,當(dāng)一個(gè)數(shù)值公式的截?cái)嗾`差可表示為O（hk+1）時(shí)（k為正整數(shù)，h為步長(zhǎng)），稱它是一個(gè)k階公式,k越大，則數(shù)值公式的精度越高,歐拉法是一階公式，改進(jìn)的歐拉法是二階公式。 龍格-庫(kù)塔法有二階公式和四階公式。 線性多步法有四階阿達(dá)姆斯外插公式和內(nèi)插公式,三）用Matlab軟件求常微分方程的數(shù)值解,t，x=solver（f,ts,x0,options,1、在解n個(gè)未知函數(shù)的方程組