常微分方程模型

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

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

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

4、.快速靜脈注射,t=0 瞬時注射劑量d 的藥物進入中心室,血藥濃度立即為 d/V1,給藥速率 f (t) 和初始條件,2.恒速靜脈滴注,t T時, c1(t)和 c2(t)按指數規(guī)律衰減趨于零,3.口服或肌肉注射,相當于藥物( 劑量d)先進入吸收室，吸收后再進入中心室,吸收室藥量x0(t,參數估計,各種給藥方式下的 c1(t), c2(t) 取決于參數k12, k21, k13, V1,V2,以快速靜脈注射為例 ,在ti(i=1,2,n)測得c1(ti,由較大的 用最小二乘法定A,由較小的 用最小二乘法定B,參數估計法一,進入中心室的藥物全部排除,參數 估計法二, 構造非線性擬合函數 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 傳染病模型,問題,描述傳染病的傳播過程,分析受感染人數的變化規(guī)律,預報傳染病高潮到來的時刻,預防傳染病蔓延的手段,

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

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

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

9、的資料中，發(fā)現鯊魚等的比例有明顯增加（見下表），而供其捕食的食用魚的百分比卻明顯下降.顯然戰(zhàn)爭使捕魚量下降，食用魚增加，鯊魚等也隨之增加，但為何鯊魚的比例大幅增加呢,他無法解釋這個現象，于是求助于其岳父，著名的意大利數學家V.Volterra，希望建立一個食餌捕食系統(tǒng)的數學模型，定性或定量地回答這個問題,該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關系，沒有考慮食餌和捕食者自身的阻滯作用，是Volterra提出的最簡單的模型,模型（一）不考慮捕獲,定理 Volterra微分方程組對應初值問題,的解,是周期函數，且解的周期平均值為,首先，建立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,求解結果,左圖反映了x1（t）與x2（t）的關系。 可以猜測： x1（t）與x2（t）都是周期函數,模型（二） 考慮人工捕獲,設表示捕獲能力的系數為e，相當于食餌的自然增長率由a降為a-e，捕食者的死亡率由c增為

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

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

13、種群的環(huán)境容納量,模型4,考慮雙方內部都存在生存競爭，且捕食者另有食物來源 E.C.Pielou,平衡點,模型,兩種群競爭模型,競爭模型,競爭排斥原理（Competition Exclution Law,多個種群依靠同一個生存資源而生活，如果生活在同一個地理空間，獵取相同食物或營養(yǎng)物。在有限的相同生存資源條件下，如果存在競爭關系，它們必然相互排斥，展開激烈的生存競爭。結局是競爭力較弱的種群滅絕，競爭力最強的種群達到其環(huán)境容納量,模型討論,競爭排斥原理 是針對前三種情形得出的結論， 第四種情況極為罕見, M 函數 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,模型,兩種群互惠模型,互惠模型,研究多個種群之間相互依賴、共生現象,模型討論,模型有三個平衡點，分別為,兩種群最終達到穩(wěn)定平衡態(tài),兩種群共生,P3

15、為正平衡點,P3 穩(wěn)定,兩種群最終達不到穩(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常微分方程的數值解及實驗,一）常微分方程數值解的定義,在生產和科研中所處理的微分方程往往很復雜且大多得不出一般解。而在實際上對初值問題，一般是要求得到解在若干個點上滿足規(guī)定精確度的近似值，或者得到一個滿足精確度要求的便于計算的表達式,因此，研究常微分方程的數值解法是十分必要的,歐拉公式,1、用差商代替導數,若步長h較小，則有,故有公式,此即歐拉法（向前歐拉法,2、使用數值積分,對方程y=f(x,y), 兩邊由xi到xi+1積分，并利用梯形公式，有,實際應用時，與歐拉公式結合使用,

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