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

數(shù)學(xué)建模三峽大學(xué)理學(xué)院俞輝Email:yuhui@精選課件第二局部動態(tài)模型動態(tài)模型介紹動態(tài)模型分析動態(tài)模型模擬精選課件第4章動態(tài)模型介紹常態(tài)分析動力系統(tǒng)離散時間動力系統(tǒng)精選課件4.1常態(tài)分析例4.1在一個未被管理的森林，硬材樹和軟材樹競爭可用的土地和水分。越可用的硬材樹生長的越慢，但越耐久且提供越有價值的木材。軟材樹靠生長快，有效消耗水分和土壤養(yǎng)分與硬材樹競爭。硬材樹靠生長的高度與軟材樹競爭，它們遮擋了小樹的陽光，它們也更耐抗疾病。這兩種樹能否同時在一片森林中共存，或者一種樹會迫使另一種樹滅絕？精選課件五步法H和S分別表示硬材樹和軟材樹種群。生物學(xué)家習(xí)慣使用的計(jì)量單位是每英畝上的木材噸數(shù)。無限制生長(豐富的空間、陽光、水分、土壤養(yǎng)料等〕：rP種群內(nèi)的競爭：-aP2〔小種群的增長率線性依賴于種群的大小，即：-aP)種群生長(率)函數(shù)：g(P)=rP-aP2，(r為內(nèi)稟增長率，a<<r是資源限制強(qiáng)度系數(shù))種群間的競爭：-bSH精選課件變量與假設(shè)變量：H=硬材樹種群〔噸/畝〕S=軟材樹種群〔噸/畝〕gH=硬材樹的生長率〔噸/英畝/年〕gS=軟材樹的生長率〔噸/英畝/年〕cH=與軟材樹競爭的損失〔噸/英畝/年〕cs=與硬材樹競爭的損失〔噸/英畝/年〕精選課件假設(shè)gH=r1H-a1H2gS==r2S-a2S2cH=b1SHcS=b2SHH>=0,S>=0r1,r2,a1,a2,b1,b2是正實(shí)數(shù)目標(biāo)：是否有H->0或S->0。精選課件第二步，選擇建模方法微分方程動力系統(tǒng)理論見書P.90精選課件第三步，構(gòu)造模型公式記x1=H,x2=S為兩個狀態(tài)變量，定義在狀態(tài)空間：{(x1,x2):x1>=0,x2>=0}定常態(tài)方程r1x1-a1x12-b1x1x2=0r2x2-a2x22-b2x1x2=0精選課件第四步，求解模型得四個解：三個解(0，0),(0，r2/a2),(r1/a1,0)在坐標(biāo)軸上第四個解在兩條直線r1-a1x1-b1x2=0r2-a2x2-b2x1=0的交點(diǎn)：精選課件如果兩條直線不相交，那么只存在3個平衡點(diǎn)。在這種情況下，兩個種群不能共存。我們希望知道在什么條件下x1>0且x2>0.假設(shè)ai>bi,2個種群共存的條件：精選課件平衡態(tài)精選課件第五步，答復(fù)以下問題對每種種群存在兩類增長限制。第一種是由于與另一種群的競爭，第二種是由于擁擠造成的同一種群內(nèi)部的競爭。因此，對每一種樹，存在著一點(diǎn)，在這一點(diǎn)數(shù)木由于擁擠會主動停止增長，且存在另一點(diǎn)，在這一點(diǎn)樹木通過競爭阻止另一種群的增長。兩種樹能夠共存的條件是每種樹在到達(dá)限制自己增長的點(diǎn)之前已經(jīng)到達(dá)它限制另一種樹增長的點(diǎn)。精選課件4.2動力系統(tǒng)例4.2藍(lán)鯨和長須鯨是生活在同一海域的相似種群，因此認(rèn)為他們之間存在競爭。藍(lán)鯨的內(nèi)稟增長率每年估計(jì)為5%，長須鯨為每年8%，環(huán)境承載力〔環(huán)境能夠支付的鯨魚的最大數(shù)量〕估計(jì)藍(lán)鯨為150000條，長須鯨為400000條。鯨魚競爭的程度是未知的。在過去的100年劇烈的捕撈已經(jīng)使鯨魚數(shù)量減少，藍(lán)鯨大約為5000條，長須鯨大約為70000條。藍(lán)鯨是否會滅絕？精選課件第一步，提出問題變量：B=藍(lán)鯨的數(shù)量F=長須鯨的數(shù)量gB=藍(lán)鯨種群的增長率〔每年〕gF=長須鯨種群的增長率〔每年〕cB=藍(lán)鯨與長須鯨競爭的影響〔每年的鯨魚數(shù)〕cF=長須鯨與藍(lán)鯨競爭的影響〔每年的鯨魚數(shù)〕精選課件假設(shè)gB=0.05B(1-B/150000)gF=0.08F(1-F/400000)cB=cF=aBFB>=0,F>=0,a是正實(shí)數(shù)目標(biāo)：確定動力系統(tǒng)是否能夠從B=5000,F=70000開始到達(dá)穩(wěn)定的平衡態(tài)。精選課件第二步，選擇建模方法連續(xù)時間動力系統(tǒng)理論見教材p.94.精選課件第三步，構(gòu)造模型公式令x1=B,x2=F,記x’1=f1(x1,x2),x’2=f2(x1,x2)f1(x1,x2)=0.05x1(1-x1/150000)-ax1x2f2(x1,x2)=0.08x2(1-x2/400000)-ax1x2狀態(tài)空間為S={(x1,x2):x1>=0,x2>=0}精選課件第四步，求解模型繪制向量場clearall,closeall,clcsymsx1x2alpha=10^(-7);f1=0.05*x1*(1-x1/150000)-alpha*x1*x2;f2=0.08*x2*(1-x2/400000)-alpha*x1*x2;[x1steady,x2steady]=solve(f1,f2);disp('Theequilibriumpointsare')disp([x1steadyx2steady])精選課件M=10;%numberofsamplespointsx1min=0;x1max=900000;%domainspecificationx2min=0;x2max=600000;[X1,X2]=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max);dX1=0.05*X1.*(1-X1/150000)-alpha*X1.*X2;%x1-componentdX2=0.08*X2.*(1-X2/400000)-alpha*X1.*X2;%x2-componentquiver(X1,X2,dX1,dX2);%matlabroutineaxis([x1minx1maxx2minx2max]);title('Directionfield(thevectorsmayberescaled!)');holdonxlabel('BlueWhales');ylabel('FinWhales');ezplot(f1,[09000000600000]),holdonezplot(f2,[09000000600000])精選課件精選課件四個平衡態(tài)解三個為：〔0，0〕，〔150000，0〕，〔0，400000〕第四個在區(qū)域內(nèi)部，為唯一穩(wěn)定的平衡態(tài)。精選課件第五步，答復(fù)以下問題只要停止捕撈，鯨魚種群將恢復(fù)到原來的水平，生態(tài)系統(tǒng)將處于穩(wěn)定的平衡態(tài)。精選課件靈敏性和穩(wěn)定性對參數(shù)a做靈敏性分析symsalphaf1=.05*x1*(1-x1/150000)-alpha*x1*x2;f2=.08*x2*(1-x2/400000)-alpha*x1*x2;[x1steady,x2steady]=solve(f1,f2)x1alpha=x1steady(4);x2alpha=x2steady(4);pretty(x1alpha),pretty(x2alpha)精選課件解之得x1=150000*(-1+8000000*a)/Dx2=400000*(-1+1875000*a)/D其中，D=-1+15000000000000*a^2alpha1=solve(x1alpha);alpha2=solve(x2alpha);formatshortedouble(alpha1),double(alpha2)可見，對任意的a<1.25*10-7存在一個穩(wěn)定的平衡態(tài)x1>0,x2>0.精選課件figureezplot(x1alpha,[08*10^(-7)]),holdongridonezplot(x2alpha,[08*10^(-7)]),holdontitle('Levelofcoexistingpopulationsvsparameter\alpha');精選課件精選課件a=linspace(0,9*10^(-7));x1a=subs(x1alpha,alpha,a);x2a=subs(x2alpha,alpha,a);ind=find(x1a>0&x2a>0);plot(a(ind),x1a(ind),'bo');holdonplot(a(ind),x2a(ind),'ro');精選課件精選課件formatbankSx1alpha=diff(x1alpha,alpha)*(alpha/x1alpha);Sx2alpha=diff(x2alpha,alpha)*(alpha/x2alpha);Sx1a=subs(Sx1alpha,alpha,10^(-7))Sx2a=subs(Sx2alpha,alpha,10^(-7))精選課件穩(wěn)健性分析穩(wěn)健性分析是考慮上述模型中f1和f2具有更一般的形式。只要向量場具有相同的一般特征，我們的結(jié)論仍然是正確的。精選課件4.3離散時間動力系統(tǒng)例4.3宇航員在訓(xùn)練中要求用手動控制做對接演習(xí)。作為這個演習(xí)的一局部，要求保持一個正在運(yùn)行的太空船與另一個正在運(yùn)行的太空船的相對位置。手控制器提供了不同的加速度和減速度，并且在太空船上有一個裝置測量這兩個飛船的接近速度。建議使用如下的策略進(jìn)行飛船對接。精選課件首先觀察接近速度。如果為零，那么不用再做任何事情。否那么，記住這個接近速度，再看加速度控制器，控制加速度使得它與接近速度相反〔即如果接近速度是正值，那么放慢，如果是負(fù)的，那么加快。〕，且正比于這個差值〔即如果發(fā)現(xiàn)接近速度到達(dá)2倍時，我們將于2倍的速度剎車〕。經(jīng)過一段時間，再觀察接近速度并重復(fù)上面的步驟。在什么環(huán)境下這個策略才是有效的？精選課件第一步，提出問題設(shè)vn表示在時間tn觀測到的接近速度，tn為第n次觀測的時間。太空船接近速度的改變：Δvn=vn+1-vn兩次觀測之間的時間間隔：Δtn=tn+1-tn時間區(qū)間被分成兩局部：Δtn=cn+wncn為調(diào)整控制器的時間，wn為下一次觀測前的等待時間。記an為第n次調(diào)節(jié)后設(shè)定的加速度，那么Δvn=an-1cn+anwn按控制律要求加速度正比于(-vn),因此，an=-kvn精選課件變量：tn=第n次觀測速度的時間(秒)vn=在tn時刻的速度(米/秒)cn=執(zhí)行第n次控制調(diào)節(jié)的時間(秒)an=第n次調(diào)節(jié)后的加速度(米/秒)wn=等待到第n+1次觀測前的等待時間(秒)精選課件假設(shè)：tn+1=tn+Δtn=tn+cn+wnvn+1=vn+Δvn=vn+an-1cn+anwnan=-kvncn>0wn>=0目標(biāo)：確定是否有vn→0.精選課件第二步，選擇建模方法離散時間動力系統(tǒng)理論見教材，p99.精選課件第三步，推導(dǎo)模型公式vn+1-vn=-kvn-1cn-kvnwn為簡化起見，對所有的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精選課件第四步，求解模型平衡態(tài)方程：-kwx1-kcx2=0x1-x2=0平衡點(diǎn)〔0，0〕位于上述兩直線的交點(diǎn)。下面繪制向量場F=(-kwx1-kcx2,x1-x2)精選課件繪制向量場clearall,closeall,clcc=5,w=10,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([x1minx1maxx2minx2max]);title('Directionfield(thevectorsarerescaled!)');holdonxlabel('CurrentSpeed'),ylabel('PreviousSpeed');精選課件精選課件系統(tǒng)演變定義函數(shù)dockfun.mfunctionrhs=dockfun(x,c,w,k);rhs=[-k*w*x(1)-k*c*x(2);x(1)-x(2)];精選課件迭代實(shí)現(xiàn)x=[8;10];N=12;fprintf('nCurrentspeedPrev.speed\n\n')fprintf('%2.0f%5.2f%5.2f\n',0,x(1),x(2))forn=1:Nxnew=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.2f\n',n,x(1),x(2))end精選課件精選課件第五步，答復(fù)以下問題假設(shè)c<<w,有Δx1≈-kwx1當(dāng)kw<2時，我們將得到一個穩(wěn)定的平衡態(tài)。當(dāng)kw<1時，系統(tǒng)將逼近平衡態(tài)，而不會越過它。據(jù)此，可以答復(fù)以下問題如下：只要控制調(diào)節(jié)不是太劇烈，那么控制將起作用。精選課件進(jìn)一步，兩次調(diào)節(jié)之間的間隔時間越長，調(diào)節(jié)幅度必須越小。而且，它們之間呈反比。如果兩次調(diào)節(jié)之間的時間間隔增加兩倍，調(diào)節(jié)幅度可以減半。特別地，如果我們以10秒調(diào)節(jié)一次，那么設(shè)置加速度為1/10，以免超過目標(biāo)速度零。精選課件作業(yè)：習(xí)題4.41，11精選課件第五章動態(tài)模型分析連續(xù)時間系統(tǒng)的特征值方法離散時間系統(tǒng)的特征值方法相圖精選課件連續(xù)時間系統(tǒng)特征值方法例5.1再次考慮例4.1的樹木問題。假設(shè)硬材樹每年增長率為10%，軟材樹每年增長率為25%。一英畝林地可以提供大約10000噸的硬木或6000噸的軟木。競爭的程度還未從數(shù)值上確定。兩種樹能否共存于一個穩(wěn)定的平衡態(tài)？精選課件第一步，提出問題r1=0.10r2=0.25a1=0.10/10000a2=0.25/6000精選課件第二步，選擇建模方法連續(xù)時間動力系統(tǒng)的特征值分析法見教材，p112精選課件第三步，推導(dǎo)模型公式已經(jīng)確定了r1,r2,a1,a2,我們?nèi)约僭O(shè)bi<ai,不妨取bi=ai/2,那么動態(tài)系統(tǒng)方程為x’=F(x),其中F(x)=(f1,f2).f1(x1,x2)=0.10x1-(0.10/10000)(x1)2-(0.05/10000)x1x2f2(x1,x2)=0.25x2-(0.25/6000)(x2)2-(0.125/6000)x1x2精選課件第四步，求解模型區(qū)域內(nèi)部有平衡點(diǎn)：x10≈9333,x20≈1333計(jì)算偏導(dǎo)數(shù)矩陣在平衡點(diǎn)(x10,x20)處的值A(chǔ)。求矩陣A的特征值，得精選課件這兩個特征值均具有負(fù)實(shí)部，所以該平衡點(diǎn)為穩(wěn)定的。上述過程也可以用計(jì)算機(jī)求解。精選課件代碼實(shí)現(xiàn)---求平衡點(diǎn)symsx1x2f1=0.1*x1*(1-(1/10000)*x1-0.5*(1/10000)*x2);f2=0.25*x2*(1-(1/6000)*x2-0.5*(1/6000)*x1);[x1steady,x2steady]=solve(f1,f2);N=length(x1steady);fprintf('Theequilibriumpointsare\n')disp([x1steadyx2steady])精選課件代碼實(shí)現(xiàn)---繪制方向場M=15;x1min=0;x1max=11000;x2min=0;x2max=8000;ezplot(f1,[x1minx1maxx2minx2max]),holdonezplot(f2,[x1minx1maxx2minx2max]),holdon[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([x1minx1maxx2minx2max]);title('Directionfield(thevectorsarerescaled!)');xlabel('Hardwoods'),ylabel('Softwoods'),holdon精選課件精選課件代碼實(shí)現(xiàn)---穩(wěn)定性分析DF=[diff(f1,x1),diff(f1,x2);diff(f2,x1),diff(f2,x2)];fori=1:Nx1num=double(x1steady(i));x2num=double(x2steady(i));A=subs(DF,[x1,x2],[x1num,x2num])lambda=eig(A);fprintf('Theeigenvaluesfortheequilibrium(')fprintf('%1.0f,%1.0f',x1num,x2num);fprintf(')are')fprintf('%1.2f%1.2f\n',lambda(1),lambda(2));plot(x1num,x2num,'ro','MarkerSize',10,'MarkerFaceColor','g');end精選課件精選課件第五步，答復(fù)以下問題我們發(fā)現(xiàn)硬材樹和軟材樹可以共存于一個平衡態(tài)。在一個成熟的穩(wěn)定的樹林中，每英畝大約有9300噸硬材樹和1300噸軟材樹。這個結(jié)論基于對兩類樹種之間競爭程度的近似合理的假設(shè)。精選課件靈敏性分析為做靈敏性分析，這里放松假設(shè)bi=(1/2)ai,而假設(shè)bi=tai.條件:bi<airi/ai<rj/bj隱含著0<t<0.6精選課件定義系統(tǒng)方程symsx1x2tf1=0.1*x1*(1-(1/10000)*x1-t*(1/10000)*x2);f2=0.25*x2*(1-(1/6000)*x2-t*(1/6000)*x1);精選課件求平衡態(tài)[x1steady,x2steady]=solve(f1,f2);N=length(x1steady);fprintf('Theequilibriumpointsare\n')disp([x1steadyx2steady])精選課件每一平衡點(diǎn)對t(0<t<0.6)的靈敏性DF=[diff(f1,x1),diff(f1,x2);diff(f2,x1),diff(f2,x2)]fori=1:NA=subs(DF,[x1,x2],[x1steady(i),x2steady(i)])lambda=eig(A)figurefortt=0:0.01:0.6y=subs(lambda,t,tt);plot(tt,y),holdonendend精選課件平衡態(tài)〔0，0〕精選課件平衡態(tài)〔150000，0〕精選課件平衡態(tài)〔0，400000〕精選課件區(qū)域內(nèi)部的平衡態(tài)精選課件5.3相圖所謂連續(xù)時間動態(tài)系統(tǒng)的相圖就是連續(xù)時間動態(tài)系統(tǒng)有代表性的解曲線在狀態(tài)空間的草