《差分方程模型》PPT課件.ppt

第五講 差分方程模型,一 差分方程模型 二 差分方程解法 三 差分方程的平衡點及穩(wěn)定性 四 建模案例 五 用Matlab求解差分方程問題,對一數(shù)列an，把數(shù)列中的an和前面的ai(0=in)關(guān)聯(lián)起來的方程叫差分方程，也叫遞推關(guān)系。 例：設(shè)第一月初有雌雄各一的一對小兔。假定兩月后長成成兔，同時從第三個月開始每月初產(chǎn)雌雄各一對一對小兔，新增小兔也按次規(guī)律繁殖。設(shè)第n月末共有Fn對兔子，試建立關(guān)于Fn的差分方程。 解：因為第n月末的兔子包括兩部分，一部分為上月留下的，另一部分為當月新生的，而由題設(shè)當月生的小兔等于前月末的兔子數(shù)，所以有 Fn=Fn-1+ Fn-2 ,F1=F2=1.返回,一 差分方程模型,1 常系數(shù)線性齊次差分方程的解法 形如an+b1an-1+b2an-2+bkan-k=0（1）(其中bi為常數(shù)，bk0,n=k.)的差分方程，稱為an的k階常系數(shù)線性齊次差分方程。 Xk+b1xk-1+bk=0為上述差分方程的特征方程，其根稱為特征根。 解分為三種情況： （1） 單根 若差分方程（1）的特征方程Xk+b1xk-1+bk=0有k個相異的特征根x1,x2,xk，則an=c1x1n+c2x2n+ckxkn是一個通解，其中ci為常數(shù)，由初始條件a0=u0,a1=u1,ak-1=uk-1可確定一個滿足初始條件的特解。,二 差分方程解法,（2） 重根 若差分方程（1）的特征方程Xk+b1xk-1+bk=0的相異特征根x1,x2,xt ，重數(shù)依次為m1,m2,mt, m1+m2+mt=k，則差分方程的通解為,（3） 共軛復(fù)根 若差分方程（1）的特征方程Xk+b1xk-1+bk=0有一對共軛復(fù)根 , 和相異的k-2個實根x3,xk，則差分方程的通解為，,其中,常系數(shù)線性非齊次差分方程的解法 形如an+b1an-1+b2an-2+bkan-k=f(n)（1）(其中bi為常數(shù)，bk0,n=k,f(n) 0)的差分方程，稱為an的k階常系數(shù)線性非齊次差分方程。 非齊次差分方程的通解等于對應(yīng)的齊次差分方程的通解加上非齊次差分方程的一個特解。,2.,返回,三、差分方程的平衡點及穩(wěn)定性 1 一階線性常系數(shù)差分方程的平衡點及穩(wěn)定性 一階線性常系數(shù)差分方程 xk+1+axk=b,k=0,1,2,（1） 的平衡點由x+ax=b解得，為 ，當 時，若xkx*，則x*是穩(wěn)定的。 方程（1）的平衡點的穩(wěn)定性問題可以通過變量代換轉(zhuǎn)換為齊次方程 xk+1+axk=0,k=0,1,2 (2),的平衡點x*=0的穩(wěn)定性問題。而對于方程（2），其解可以表示為 xk=(-a)kx0, k=1,2, (3) 所以當且僅當|a|1時，方程（2）（從而方程（1）的平衡點是穩(wěn)定的。,對于n維向量x(k)和nn常數(shù)矩陣A構(gòu)成的方程組 x(k+1)+Ax(k)=0 其平衡點穩(wěn)定的條件是A的特征根i,I=1,2,，均有|i|1。,2 二階線性差分方程的平衡點及穩(wěn)定性 考察二階線性差分方程xk+a1xk+1+a2xk+2=0 (4) 在平衡點x*=0的穩(wěn)定性。為求（4）的通解，先寫出他的特征方程 記它的根為1，2，則（4）的通解可以表示為 ，其中常數(shù)c1,c2由初始條件x0,x1確定，從而可知，當且僅當|1|1, |2|1時方程（4）的平衡點是穩(wěn)定的。,3 一階非線性差分方程的平衡點及穩(wěn)定性 考察一階非線性差分方程xk+1=f(xk) (7) 的平衡點的穩(wěn)定性。其平衡點x*由x=f(x)解出。將（7）的右端在x*點做泰勒展開，只取一次項，則（7）可以近似為： （8） x*也是（8）的平衡點。線性方程（8） 的平衡點的穩(wěn)定性討論同（1），而當|f(x*)|1時（7）與（8）的平衡點的穩(wěn)定性相同。從而有： 當|f(x*)|1時，方程（7）的平衡點是不穩(wěn)定的。,返回,四、建模案例-最優(yōu)捕魚策略 問題簡介 生態(tài)學(xué)原理：對可再生資源的開發(fā)策略應(yīng)為在可持續(xù)收獲的前提下追求最大經(jīng)濟效益。 考慮4個年齡組：1齡魚，2齡魚，3齡魚，4齡魚的某魚類。該魚類在每年后4個月產(chǎn)卵繁殖。因而捕撈只能在前8個月進行。每年投入的捕撈能力不變，單位時間捕撈量與各年齡組魚群條數(shù)的比例稱為捕撈強度系數(shù)。且只能捕撈3、4齡魚，兩個捕撈強度系數(shù)比為0.42：1。即為固定努力量捕撈。,魚群數(shù)據(jù)為： (1) 各年齡組魚的自然死亡率為0.8（1/年），其平均重量分別為5.07，11.55，17.86，22.99（g）; (2) 1齡魚和2齡魚不產(chǎn)卵，產(chǎn)卵期間，平均每條4齡魚產(chǎn)卵量為1.109105（個），3齡魚為其一半； (3) 卵孵化的成活率為1.221011/（1.221011+n）(n為產(chǎn)卵總量)；,問題描述如下： 如何實現(xiàn)可持續(xù)捕獲（即每年開始捕撈時漁場中各年齡組魚群不變），并在此前提下得到最高收獲量； 合同要求某漁業(yè)公司在5年合同期滿后魚群的生產(chǎn)能力不能受到太大破壞，承包時各年齡組魚群數(shù)量為122，29.7，10.1，3.29（109條）。在固定努力量的捕撈方式下，問該公司應(yīng)采取怎樣的捕撈策略，才能使總收獲量最高。,(1) 引入變量 X=(X1,X2,X3,X4)T 為魚群數(shù)向量； -單位時間的自然死亡率； c-年存活率，c=1-0.8=0.2; k-單位時間4齡魚的捕撈強度系數(shù)； -孵化卵成活率，=1.221011/（1.221011+n）； m-4齡魚的平均產(chǎn)卵量，m為1.109105（個），3齡魚為其一半。,(2)模型的建立 以1年為一個離散化的時間單位。 記年初魚群為X(t)=(X1(t), X2(t), X3(t), X4(t)T, 下一年的魚群數(shù)為X(t+1)=(X1(t+1), X2(t+1), X3(t+1), X4(t+1)T。顯然，Xi(t+1)是Xi-1(t+1)到年底存活下來的魚群數(shù)（i=1,2,3,i=4時X4(t+1)中還包括X4(t)中的存活數(shù)。X0(t)是指上一年由卵孵化而得到的1齡魚），據(jù)此可建立如下差分方程： X2(t+1)=c X1(t); X3(t+1)= c X2(t); X4(t+1)=(c-k3)X3(t)+(c-k4)X4(t);,因為3、4齡魚的捕撈強度系數(shù)比為0.42：1，所以有k3=0.42k4=0.42k,寫成矩陣形式有： X(t+1)=PX(t); 其中 當4齡魚的捕撈強度系數(shù)kc/0.42時，不論上一年魚群數(shù)目如何，下一年魚群將出現(xiàn)負數(shù)。說明模型存在問題，原因是離散化程度不夠精細。,假設(shè)單位時間為一個月，定義月死亡率為，月存活率為（1-）， 月捕撈系數(shù)為k,則年存活率為（1-）12=c=0.2,從而=0.1255。 考慮一年中各月魚群數(shù)目的分布，則有： 一個月的實際存活率：（1-k）； 兩個月的實際存活率：（1-k）2； 三個月的實際存活率：（1-k）3； 。 八個月的實際存活率：（1-k）8； 九個月的實際存活率：（1-k）8（1-）； 。 一年后實際存活率：（1-k）8（1-）4。 同理可得第i月的捕撈率: （1-k）i-1k,i=1,2,8.,從而有： 一年后3齡魚實際存活數(shù)：（1-k3）8（1-）4X3; 一年后4齡魚實際存活數(shù)：（1-k4）8（1-）4X4; 該年3齡魚總捕撈量： ， 該年4齡魚總捕撈量： ， 該年3齡魚產(chǎn)卵總量： ； 該年4齡魚產(chǎn)卵總量： ；,因此矩陣應(yīng)修正為： 關(guān)于魚群的差分方程為：X(t+1)=PX(t) （1） 為實現(xiàn)持續(xù)捕獲，（1）式必須存在穩(wěn)定解：X(t)=PX(t)。 由差分方程穩(wěn)定性理論知其充要條件為：對P的所有特征根i,均有|i|1。由此可求得最佳策略。,五 用Matlab求解差分方程問題,1、一階線性常系數(shù)差分方程 2、高階線性常系數(shù)差分方程,1、一階線性常系數(shù)差分方程,瀕危物種的自然演變和人工孵化 問題： Florida沙丘鶴屬于瀕危物種，它在較好自然環(huán)境下，年均增長率僅為1.94%，而在中等和較差環(huán)境下年均增長率分別為 -3.24% 和 -3.82%，如果在某自然保護區(qū)內(nèi)開始有100只鶴，建立描述其數(shù)量變化規(guī)律的模型，并作數(shù)值計算。,模型建立,記第k年沙丘鶴的數(shù)量為xk,年均增長率為r，則第k+1年鶴的數(shù)量為 xk+1=(1+r)xk k=0,1,2 已知x0=100, 在較好，中等和較差的自然環(huán)境下 r=0.0194, -0.0324,和-0.0382 我們利用Matlab編程，遞推20年后觀察沙丘鶴的數(shù)量變化情況,Matlab實現(xiàn),首先建立一個關(guān)于變量n ，r的函數(shù) function x=sqh(n,r) a=1+r; x=100; for k=1:n x(k+1)=a*x(k); end,在command窗口里調(diào)用sqh函數(shù),k=(0:20); y1=sqh(20,0.0194); y2=sqh(20,-0.0324); y3=sqh(20,-0.0382); round(k,y1,y2,y3),利用plot 繪圖觀察數(shù)量變化趨勢,可以用不同線型和顏色繪圖 r g b c m y k w 分別表示 紅綠蘭蘭綠洋紅黃黑白色 ： + o * . X s d 表示不同的線型,plot(k,y1,k,y2,k,y3) 在同一坐標系下畫圖,plot(k,y2,:) plot(k,y2,-) plot(k,y2,r) plot(k,y2,y) plot(k,y2,y,k,y1,:) plot(k,y2,k,y1,:) plot(k,y2,oy,k,y1,:) 用gtext(r=0.0194),gtext(r=-0.0324),gtext(r=-0.0382)在圖上做標記。,人工孵化是挽救瀕危物種的措施之一，如果每年孵化5只鶴放入保護區(qū)，觀察在中等自然條件下沙丘鶴的數(shù)量如何變化 Xk+1=aXk +5 ,a=1+r 如果我們想考察每年孵化多少只比較合適，可以令 Xk+1=aXk +b ,a=1+r,function x=fhsqh(n,r,b) a=1+r; x=100; for k=1:n x(k+1)=a*x(k)+b; end,k=(0:20) ; %一個行向量 y1=(20,-0.0324,5); 也是一個行向量 round( k , y 1 ) 對k,y1四舍五入，但 是 不改變變量的值 plot( k , y1) k y1 是行向量列向量都可以 也可以觀察200年的發(fā)展趨勢，以及在較差條件下的發(fā)展趨勢，也可以考察每年孵化數(shù)量變化的影響。,2、高階線性常系數(shù)差分方程,如果第k+1時段變量Xk+1不僅取決于第k時段變量Xk，而且與以前時段變量有關(guān)，就要用高階差分方程來描述,一年生植物的繁殖,一年生植物春季發(fā)芽，夏天開花，秋季產(chǎn)種，沒有腐爛，風干，被人為掠取的那些種子可以活過冬天，其中一部分能在第2年春季發(fā)芽，然后開花，產(chǎn)種，其中的另一部分雖未能發(fā)芽，但如又能活過一個冬天，則其中一部分可在第三年春季發(fā)芽，然后開花，產(chǎn)種，如此繼續(xù)，一年生植物只能活1年，而近似的認為，種子最多可以活過兩個冬天，試建立數(shù)學(xué)模型研究這種植物數(shù)量變化的規(guī)律，及它能一直繁殖下去的條件。,模型及其求解,記一棵植物春季產(chǎn)種的平均數(shù)為c,種子能活過一個冬天的(1歲種子)比例為b,活過一個冬天沒有發(fā)芽又活過一個冬天的（2歲種子）比例仍為b,1歲種子發(fā)芽率a1，2歲種子發(fā)芽率a2。 設(shè)c,a1,a2固定，b是變量，考察能一直繁殖的條件 記第k年植物數(shù)量為Xk，顯然Xk與Xk-1,Xk-2有關(guān)，由 Xk-1決定的部分是 a1bcXk-1,由Xk-2決定的部分是 a2b(1-a1)bcXk-2,Xk= a1bcXk-1 + a2b(1-a1)bcXk-2,Xk= a1bcXk-1 + a2b(1-a1)bcXk-2,實際上，就是Xk= pXk-1 + qXk-2 我們需要知道x0,a1,a2,c, 考察b不同時，種子繁殖的情況。在這里假設(shè) X0=100,a1=0.5,a2=0.25,c=10,b=0.180.20 這樣可以用matlab計算了,Xk= a1bcXk-1 + a2b(1-a1)bcXk-2,function x=zwfz(x0,n,b) c=10;a1=0.5;a2=0.25; p=a1*b*c;q=a2*b*(1-a1)*b*c; x1=x0; x2=p*x1; for k=3:n x(k)=p*x(k-1)+q*x(k-2); end,K=(3:20); y1=zwfz(100,21,0.18); y2=zwfz(100,21,0.19); y3=zwfz(100,21,0,20); round