數(shù)學(xué)建模 第二部分 動(dòng)態(tài)模型

1、編輯課件1數(shù)學(xué)建模三峽大學(xué) 理學(xué)院 俞 輝Email:編輯課件2第二部分 動(dòng)態(tài)模型 動(dòng)態(tài)模型介紹 動(dòng)態(tài)模型分析 動(dòng)態(tài)模型模擬編輯課件3第4章 動(dòng)態(tài)模型介紹 常態(tài)分析 動(dòng)力系統(tǒng) 離散時(shí)間動(dòng)力系統(tǒng)編輯課件44.1 常態(tài)分析 例4.1 在一個(gè)未被管理的森林，硬材樹和軟材樹競(jìng)爭(zhēng)可用的土地和水分。越可用的硬材樹生長(zhǎng)的越慢，但越耐久且提供越有價(jià)值的木材。軟材樹靠生長(zhǎng)快，有效消耗水分和土壤養(yǎng)分與硬材樹競(jìng)爭(zhēng)。硬材樹靠生長(zhǎng)的高度與軟材樹競(jìng)爭(zhēng)，它們遮擋了小樹的陽光，它們也更耐抗疾病。這兩種樹能否同時(shí)在一片森林中共存，或者一種樹會(huì)迫使另一種樹滅絕？編輯課件5五步法 H和S分別表示硬材樹和軟材樹種群。生物學(xué)家習(xí)慣使用

2、的計(jì)量單位是每英畝上的木材噸數(shù)。 無限制生長(zhǎng)(豐富的空間、陽光、水分、土壤養(yǎng)料等）：rP 種群內(nèi)的競(jìng)爭(zhēng)：-aP2（小種群的增長(zhǎng)率線性依賴于種群的大小 ，即：-aP) 種群生長(zhǎng)(率)函數(shù)：g(P)= rP -aP2，(r為內(nèi)稟增長(zhǎng)率，a=0,S=0 r1,r2,a1,a2,b1,b2是正實(shí)數(shù) 目標(biāo)：是否有H-0或S-0。編輯課件8第二步，選擇建模方法 微分方程動(dòng)力系統(tǒng)理論 見書P.90編輯課件9第三步，構(gòu)造模型公式 記x1=H,x2=S為兩個(gè)狀態(tài)變量，定義在狀態(tài)空間：(x1,x2):x1=0,x2=0 定常態(tài)方程r1x1-a1x12-b1x1x2=0r2x2-a2x22-b2x1x2=0編輯課件

3、10第四步，求解模型 得四個(gè)解： 三個(gè)解(0，0),(0，r2/a2),(r1/a1,0)在坐標(biāo)軸上 第四個(gè)解在兩條直線 r1-a1x1-b1x2=0 r2-a2x2-b2x1=0 的交點(diǎn)：122 11 22 112121 2121 2,rar brb rxxaabbbbaaa編輯課件11 如果兩條直線不相交，則只存在3個(gè)平衡點(diǎn)。在這種情況下，兩個(gè)種群不能共存。 我們希望知道在什么條件下x10且x20. 假設(shè)aibi,2個(gè)種群共存的條件：21122112rrrrabab且編輯課件12平衡態(tài)1x2x11rb22rb22ra11ra編輯課件13第五步，回答問題 對(duì)每種種群存在兩類增長(zhǎng)限制。第一種是

4、由于與另一種群的競(jìng)爭(zhēng)，第二種是由于擁擠造成的同一種群內(nèi)部的競(jìng)爭(zhēng)。因此，對(duì)每一種樹，存在著一點(diǎn)，在這一點(diǎn)數(shù)木由于擁擠會(huì)主動(dòng)停止增長(zhǎng)，且存在另一點(diǎn)，在這一點(diǎn)樹木通過競(jìng)爭(zhēng)阻止另一種群的增長(zhǎng)。兩種樹能夠共存的條件是每種樹在達(dá)到限制自己增長(zhǎng)的點(diǎn)之前已經(jīng)達(dá)到它限制另一種樹增長(zhǎng)的點(diǎn)。編輯課件144.2 動(dòng)力系統(tǒng) 例 4.2 藍(lán)鯨和長(zhǎng)須鯨是生活在同一海域的相似種群，因此認(rèn)為他們之間存在競(jìng)爭(zhēng)。藍(lán)鯨的內(nèi)稟增長(zhǎng)率每年估計(jì)為5%，長(zhǎng)須鯨為每年8%，環(huán)境承載力（環(huán)境能夠支付的鯨魚的最大數(shù)量）估計(jì)藍(lán)鯨為150000條，長(zhǎng)須鯨為400000條。鯨魚競(jìng)爭(zhēng)的程度是未知的。在過去的100年劇烈的捕撈已經(jīng)使鯨魚數(shù)量減少，藍(lán)鯨大約為

5、5000條，長(zhǎng)須鯨大約為70000條。藍(lán)鯨是否會(huì)滅絕？編輯課件15第一步，提出問題 變量： B=藍(lán)鯨的數(shù)量 F=長(zhǎng)須鯨的數(shù)量 gB=藍(lán)鯨種群的增長(zhǎng)率（每年） gF=長(zhǎng)須鯨種群的增長(zhǎng)率（每年） cB=藍(lán)鯨與長(zhǎng)須鯨競(jìng)爭(zhēng)的影響（每年的鯨魚數(shù)） cF=長(zhǎng)須鯨與藍(lán)鯨競(jìng)爭(zhēng)的影響（每年的鯨魚數(shù)）編輯課件16 假設(shè) gB=0.05B(1-B/150000) gF=0.08F(1-F/400000) cB=cF=aBF B=0,F=0, a是正實(shí)數(shù) 目標(biāo)：確定動(dòng)力系統(tǒng)是否能夠從B=5000, F=70000開始達(dá)到穩(wěn)定的平衡態(tài)。編輯課件17第二步，選擇建模方法 連續(xù)時(shí)間動(dòng)力系統(tǒng)理論 見教材p.94.編輯課件18

6、第三步，構(gòu)造模型公式 令x1=B,x2=F,記x1=f1(x1,x2),x2=f2(x1,x2) f1(x1,x2)=0.05x1(1-x1/150000)-ax1x2 f2(x1,x2)=0.08x2(1-x2/400000)-ax1x2 狀態(tài)空間為 S=(x1,x2):x1=0,x2=0編輯課件19第四步，求解模型 繪制向量場(chǎng)clear all, close all, clc syms x1 x2 alpha=10(-7); f1=0.05*x1*(1-x1/150000) - alpha*x1*x2;f2=0.08*x2*(1-x2/400000) - alpha*x1*x2; x1st

7、eady,x2steady=solve(f1,f2); disp(The equilibrium points are) disp(x1steady x2steady) 編輯課件20M=10; % number of samples points x1min=0; x1max=900000; % domain specification x2min=0; x2max=600000; X1,X2=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max); dX1=0.05*X1.*(1-X1/150000) - alpha

8、*X1.*X2; % x1-componentdX2=0.08*X2.*(1-X2/400000) - alpha*X1.*X2; % x2-component quiver(X1,X2,dX1,dX2); % matlab routine axis(x1min x1max x2min x2max); title(Direction field (the vectors may be rescaled!); hold on xlabel(Blue Whales); ylabel(Fin Whales); ezplot(f1,0 900000 0 600000), hold on ezplot(

9、f2,0 900000 0 600000) 編輯課件210123456789x 1050123456x 1052/25 x2 (1-1/400000 x2)-1/10000000 x1 x2 = 0 x1x2編輯課件22四個(gè)平衡態(tài)解 三個(gè)為： （0，0），（150000，0），（0，400000） 第四個(gè)在區(qū)域內(nèi)部，為唯一穩(wěn)定的平衡態(tài)。編輯課件23第五步，回答問題 只要停止捕撈，鯨魚種群將恢復(fù)到原來的水平，生態(tài)系統(tǒng)將處于穩(wěn)定的平衡態(tài)。編輯課件24靈敏性和穩(wěn)定性 對(duì)參數(shù)a做靈敏性分析syms alpha f1=.05*x1*(1-x1/150000)-alpha*x1*x2; f2=.08*x2

10、*(1-x2/400000) - alpha*x1*x2; x1steady,x2steady = solve(f1,f2) x1alpha=x1steady(4); x2alpha=x2steady(4); pretty(x1alpha), pretty(x2alpha) 編輯課件25 解之得x1=150000*(-1+8000000*a)/Dx2=400000*(-1+1875000*a)/D其中，D=-1+15000000000000*a2alpha1 = solve(x1alpha); alpha2 = solve(x2alpha); format short e double(alp

11、ha1), double(alpha2) 可見，對(duì)任意的a0,x20.編輯課件26figure ezplot(x1alpha,0 8*10(-7), hold on grid on ezplot(x2alpha,0 8*10(-7), hold on title(Level of coexisting populations vs parameter alpha); 編輯課件27012345678x 10-7-6-4-20246810 x 105Level of coexisting populations vs parameter 編輯課件28a = linspace(0,9*10(-7);

12、 x1a=subs(x1alpha,alpha,a); x2a=subs(x2alpha,alpha,a); ind=find(x1a0 & x2a0); plot(a(ind),x1a(ind),bo); hold on plot(a(ind), x2a(ind),ro); 編輯課件29012345678x 10-7-6-4-20246810 x 105Level of coexisting populations vs parameter 編輯課件30format bank Sx1alpha = diff(x1alpha, alpha)*(alpha/x1alpha); Sx2al

13、pha = diff(x2alpha, alpha)*(alpha/x2alpha); Sx1a = subs(Sx1alpha, alpha, 10(-7) Sx2a = subs(Sx2alpha, alpha, 10(-7) 編輯課件31穩(wěn)健性分析 穩(wěn)健性分析是考慮上述模型中f1和f2具有更一般的形式。 只要向量場(chǎng)具有相同的一般特征，我們的結(jié)論仍然是正確的。編輯課件324.3 離散時(shí)間動(dòng)力系統(tǒng) 例4.3 宇航員在訓(xùn)練中要求用手動(dòng)控制做對(duì)接演習(xí)。作為這個(gè)演習(xí)的一部分，要求保持一個(gè)正在運(yùn)行的太空船與另一個(gè)正在運(yùn)行的太空船的相對(duì)位置。手控制器提供了不同的加速度和減速度，并且在太空船上有一個(gè)裝置

14、測(cè)量這兩個(gè)飛船的接近速度。建議使用如下的策略進(jìn)行飛船對(duì)接。編輯課件33 首先觀察接近速度。如果為零，則不用再做任何事情。否則，記住這個(gè)接近速度，再看加速度控制器，控制加速度使得它與接近速度相反（即如果接近速度是正值，則放慢，如果是負(fù)的，則加快。），且正比于這個(gè)差值（即如果發(fā)現(xiàn)接近速度達(dá)到2倍時(shí)，我們將于2倍的速度剎車）。經(jīng)過一段時(shí)間，再觀察接近速度并重復(fù)上面的步驟。在什么環(huán)境下這個(gè)策略才是有效的？編輯課件34第一步，提出問題 設(shè)vn表示在時(shí)間tn觀測(cè)到的接近速度，tn為第n次觀測(cè)的時(shí)間。 太空船接近速度的改變：vn =vn+1-vn 兩次觀測(cè)之間的時(shí)間間隔：tn=tn+1-tn 時(shí)間區(qū)間被分成

15、兩部分： tn=cn+wn cn為調(diào)整控制器的時(shí)間，wn為下一次觀測(cè)前的等待時(shí)間。 記an為第n次調(diào)節(jié)后設(shè)定的加速度，則vn =an-1cn+anwn 按控制律要求加速度正比于(-vn),因此，an=-kvn編輯課件35變量： tn=第n次觀測(cè)速度的時(shí)間(秒) vn=在tn時(shí)刻的速度(米/秒) cn=執(zhí)行第n次控制調(diào)節(jié)的時(shí)間(秒) an=第n次調(diào)節(jié)后的加速度(米/秒) wn=等待到第n+1次觀測(cè)前的等待時(shí)間(秒)編輯課件36假設(shè)： tn+1=tn+tn=tn+cn+wn vn+1=vn+vn=vn+an-1cn+anwn an=-kvn cn0 wn=0 目標(biāo)：確定是否有vn0.編輯課件37第

16、二步，選擇建模方法 離散時(shí)間動(dòng)力系統(tǒng)理論 見教材，p99.編輯課件38第三步，推導(dǎo)模型公式 vn+1-vn=-kvn-1cn-kvnwn 為簡(jiǎn)化起見，對(duì)所有的n，設(shè)cn=c,wn=w 有vn+1-vn=-kwvn-kcvn-1 設(shè)x1(n)=vn,x2(n)=vn-1,則 x1 =-kwx1-kcx2 x2=x1-x2編輯課件39第四步，求解模型 平衡態(tài)方程： -kwx1-kcx2=0 x1-x2=0 平衡點(diǎn)（0，0）位于上述兩直線的交點(diǎn)。 下面繪制向量場(chǎng)F=(-kwx1-kcx2, x1-x2)編輯課件40繪制向量場(chǎng)clear all, close all, clc c = 5,w = 10

17、,k = 0.1; M=10,x1min=-15,x1max=15,x2min=-15,x2max=15;X1,X2=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max);dX1=-k*w*X1-k*c*X2; dX2= X1-X2; quiver(X1,X2,dX1,dX2); axis(x1min x1max x2min x2max); title(Direction field (the vectors are rescaled!); hold on xlabel(Current Speed), ylabel(

18、Previous Speed); 編輯課件41-15-10-5051015-15-10-5051015Direction field (the vectors are rescaled!)Current SpeedPrevious Speed編輯課件42系統(tǒng)演變 定義函數(shù)dockfun.mfunction rhs = dockfun(x,c,w,k); rhs = -k*w*x(1)-k*c*x(2); x(1)-x(2);編輯課件43迭代實(shí)現(xiàn)x=8;10; N = 12; fprintf(n Current speed Prev. speednn) fprintf(%2.0f %5.2f %

19、5.2fn, 0, x(1), x(2)for n=1:N xnew = x + dockfun(x,c,w,k); plot(x(1),xnew(1), x(2),xnew(2),-ro,. MarkerFaceColor,k,MarkerSize,2) x = xnew; fprintf(%2.0f %5.2f %5.2fn, n, x(1), x(2) end 編輯課件44-15-10-5051015-15-10-5051015Direction field (the vectors are rescaled!)Current SpeedPrevious Speed編輯課件45第五步，回

20、答問題 假設(shè)cw,有x1 -kwx1 當(dāng)kw2時(shí)，我們將得到一個(gè)穩(wěn)定的平衡態(tài)。 當(dāng)kw1時(shí)，系統(tǒng)將逼近平衡態(tài)，而不會(huì)越過它。 據(jù)此，可以回答問題如下： 只要控制調(diào)節(jié)不是太劇烈，則控制將起作用。編輯課件46 進(jìn)一步，兩次調(diào)節(jié)之間的間隔時(shí)間越長(zhǎng)，調(diào)節(jié)幅度必須越小。而且，它們之間呈反比。如果兩次調(diào)節(jié)之間的時(shí)間間隔增加兩倍，調(diào)節(jié)幅度可以減半。 特別地，如果我們以10秒調(diào)節(jié)一次，則設(shè)置加速度為1/10，以免超過目標(biāo)速度零。編輯課件47 作業(yè)： 習(xí)題4.4 1，11編輯課件48第五章 動(dòng)態(tài)模型分析 連續(xù)時(shí)間系統(tǒng)的特征值方法 離散時(shí)間系統(tǒng)的特征值方法 相圖編輯課件49連續(xù)時(shí)間系統(tǒng)特征值方法 例5.1 再次考

21、慮例4.1的樹木問題。假設(shè)硬材樹每年增長(zhǎng)率為10%，軟材樹每年增長(zhǎng)率為25%。一英畝林地可以提供大約10000噸的硬木或6000噸的軟木。競(jìng)爭(zhēng)的程度還未從數(shù)值上確定。兩種樹能否共存于一個(gè)穩(wěn)定的平衡態(tài)？編輯課件50第一步，提出問題 r1=0.10 r2=0.25 a1=0.10/10000 a2=0.25/6000編輯課件51第二步，選擇建模方法 連續(xù)時(shí)間動(dòng)力系統(tǒng)的特征值分析法 見教材，p112編輯課件52第三步，推導(dǎo)模型公式 已經(jīng)確定了r1,r2,a1,a2,我們?nèi)约僭O(shè)biai,不妨取bi=ai/2,則 動(dòng)態(tài)系統(tǒng)方程為x=F(x),其中F(x)=(f1,f2). f1(x1,x2)=0.10

22、x1-(0.10/10000)(x1)2-(0.05/10000)x1x2 f2(x1,x2)=0.25x2-(0.25/6000)(x2)2-(0.125/6000)x1x2編輯課件53第四步，求解模型 區(qū)域內(nèi)部有平衡點(diǎn)： x109333,x201333 計(jì)算偏導(dǎo)數(shù)矩陣在平衡點(diǎn)(x10,x20)處的值A(chǔ)。 求矩陣A的特征值，得671339900編輯課件54 這兩個(gè)特征值均具有負(fù)實(shí)部，所以該平衡點(diǎn)為穩(wěn)定的。 上述過程也可以用計(jì)算機(jī)求解。編輯課件55代碼實(shí)現(xiàn)-求平衡點(diǎn)syms x1 x2 f1 = 0.1*x1*(1-(1/10000)*x1- 0.5*(1/10000)*x2); f2 = 0

23、.25*x2*(1-(1/6000)*x2- 0.5*(1/6000)*x1); x1steady,x2steady = solve(f1,f2); N=length(x1steady); fprintf(The equilibrium points aren) disp(x1steady x2steady) 編輯課件56代碼實(shí)現(xiàn)-繪制方向場(chǎng)M=15; x1min=0; x1max=11000; x2min=0; x2max=8000; ezplot(f1,x1min x1max x2min x2max),hold on ezplot(f2,x1min x1max x2min x2max),h

24、old on X1,X2=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max); dX1 = 0.1*X1.*(1-(1/10000)*X1-0.5*(1/10000)*X2); dX2 = 0.25*X2.*(1-(1/6000)*X2-0.5*(1/6000)*X1); quiver(X1,X2,dX1,dX2),axis(x1min x1max x2min x2max); title(Direction field (the vectors are rescaled!); xlabel(Hardwoods),y

25、label(Softwoods),hold on 編輯課件5701000 20003000 40005000 60007000 80009000 10000 11000010002000300040005000600070008000HardwoodsSoftwoodsDirection field (the vectors are rescaled!)編輯課件58代碼實(shí)現(xiàn)-穩(wěn)定性分析DF = diff(f1,x1), diff(f1,x2); diff(f2,x1), diff(f2,x2); for i=1:N x1num = double(x1steady(i); x2num = dou

26、ble(x2steady(i); A=subs(DF,x1,x2,x1num,x2num) lambda = eig(A); fprintf(The eigenvalues for the equilibrium () fprintf(%1.0f, %1.0f,x1num, x2num); fprintf() are ) fprintf(%1.2f %1.2fn,lambda(1), lambda(2); plot(x1num,x2num,ro,MarkerSize,10, MarkerFaceColor,g); end 編輯課件5901000 20003000 40005000 600070

27、00 80009000 10000 11000010002000300040005000600070008000HardwoodsSoftwoodsDirection field (the vectors are rescaled!)編輯課件60第五步，回答問題 我們發(fā)現(xiàn)硬材樹和軟材樹可以共存于一個(gè)平衡態(tài)。在一個(gè)成熟的穩(wěn)定的樹林中，每英畝大約有9300噸硬材樹和1300噸軟材樹。這個(gè)結(jié)論基于對(duì)兩類樹種之間競(jìng)爭(zhēng)程度的近似合理的假設(shè)。編輯課件61靈敏性分析 為做靈敏性分析，這里放松假設(shè)bi=(1/2)ai,而假設(shè)bi=tai. 條件: biai ri/airj/bj 隱含著0t0.6編輯課件62定

28、義系統(tǒng)方程syms x1 x2 tf1 = 0.1*x1*(1-(1/10000)*x1- t*(1/10000)*x2); f2 = 0.25*x2*(1-(1/6000)*x2- t*(1/6000)*x1); 編輯課件63求平衡態(tài)x1steady,x2steady = solve(f1,f2); N=length(x1steady); fprintf(The equilibrium points aren) disp(x1steady x2steady) 編輯課件64每一平衡點(diǎn)對(duì)t(0t0.6)的靈敏性DF = diff(f1,x1), diff(f1,x2); diff(f2,x1),

29、 diff(f2,x2) for i=1:N A=subs(DF,x1,x2,x1steady(i) , x2steady(i) lambda = eig(A) figurefor tt=0:0.01:0.6y=subs(lambda,t,tt);plot(tt,y),hold onend end 編輯課件65平衡態(tài)（0，0）00.10.20.30.40.50.60.70.10.150.20.25編輯課件66平衡態(tài)（150000，0）00.10.20.30.40.50.60.7-0.1-0.0500.050.10.150.20.25編輯課件67平衡態(tài)（0，400000）00.10.20.30.

30、40.50.60.7-0.25-0.2-0.15-0.1-0.0500.050.1編輯課件68區(qū)域內(nèi)部的平衡態(tài)00.10.20.30.40.50.60.7-0.25-0.2-0.15-0.1-0.050編輯課件695.3 相圖 所謂連續(xù)時(shí)間動(dòng)態(tài)系統(tǒng)的相圖就是連續(xù)時(shí)間動(dòng)態(tài)系統(tǒng)有代表性的解曲線在狀態(tài)空間的草圖。 通常與向量場(chǎng)的勾畫一起運(yùn)用，從而獲得動(dòng)態(tài)系統(tǒng)在狀態(tài)空間動(dòng)態(tài)行為的圖形描述。編輯課件70例5.3 考慮如下圖所示的電路圖。電路由一個(gè)電容，一個(gè)電阻和一個(gè)電感器構(gòu)成一個(gè)簡(jiǎn)單閉路。電路中每個(gè)元件的作用由在這個(gè)回路中的電流和電壓之間的關(guān)系表示。一個(gè)理想的物理模型給出這個(gè)關(guān)系：( )CCRRLLdv

31、Cidtvf idivdt（電容）（電阻）L（電感）編輯課件71RCL編輯課件72 稱函數(shù)f(x)為電阻的v-i特征。在古典的RLC電路理論中我們假設(shè)f(x)=Rx,其中R表示電阻?；鶢柣舴螂娏鞫ɡ碚f明：進(jìn)入一個(gè)節(jié)點(diǎn)的電流之和等于流出電流之和?；鶢柣舴螂妷憾ɡ碚f明：閉路上所有電壓差之和為零。 對(duì)情形L=1,C=1/3和f(x)=x3+4x確定這個(gè)電路隨時(shí)間變化的行為。編輯課件73第一步，提出問題 變量：vC,iC,vR,iR,vL,iL 假設(shè) CdvC/dt=iC vR=f(iR) LdiL/dt=vL iC=iR=iL,vC+vR+vL=0 L=1,C=1/3,f(x)=x3+4x 目標(biāo)：確定六個(gè)變量隨時(shí)間變化的性質(zhì)。編輯課件74第二步，選擇建模方法 建立連