1、用用MatlabMatlab求解差分方程問題求解差分方程問題一階線性常系數(shù)差分方程高階線性常系數(shù)差分方程線性常系數(shù)差分方程組差分方程是在離散時(shí)段上描述現(xiàn)實(shí)世界中變化過程的數(shù)學(xué)模型 例1、 某種貨幣1年期存款的年利率是r ，現(xiàn)存入M元，問年后的本金與利息之和是多少？ Xk+1=(1+r)xk , k = 0 , 1 , 2 以k=0時(shí)x0=M代入，遞推n次可得n年后本息為 1nnxrM 污水處理廠每天可將處理池的污水濃度降低一個(gè)固定比例q，問多長時(shí)間才能將污水濃度降低一半？ 記第k天的污水濃度為ck,則第k+1天的污水濃度為 ck+1=(1-q)ck,k=0,1,2, 從k=0開始遞推n次得 以

2、cn=c0/2代入即求解。0(1)nncqc一階線性常系數(shù)差分方程 瀕危物種的自然演變和人工孵化 問題 Florida沙丘鶴屬于瀕危物種，它在較好自然環(huán)境下，年均增長率僅為1.94%，而在中等和較差環(huán)境下年均增長率分別為 -3.24% 和 -3.82%，如果在某自然保護(hù)區(qū)內(nèi)開始有100只鶴，建立描述其數(shù)量變化規(guī)律的模型，并作 數(shù)值計(jì)算。模型建立 記第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年

3、后觀察沙丘鶴的數(shù)量變化情況Matlab實(shí)現(xiàn) 首先建立一個(gè)關(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(

4、k,y1,k,y2,k,y3) 在同一坐標(biāo)系下畫圖 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)在圖上做標(biāo)記。 人工孵化是挽救瀕危物種的措施之一，如果每年孵化5只鶴放入保護(hù)區(qū)，觀察在中等自然條件下沙丘鶴的數(shù)量如何變化Xk+1=aXk +5 ,a=1+r如果我們想考察每年孵化多少只比較合適，可以令Xk+1=aXk +b ,a

5、=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) ; %一個(gè)行向量 y1=(20,-0.0324,5); 也是一個(gè)行向量 round( k , y 1 ) 對k,y1四舍五入，但 是 不改變變量的值 plot( k , y1) k y1 是行向量列向量都可以也可以觀察200年的發(fā)展趨勢，以及在較差條件下的發(fā)展趨勢，也可以考察每年孵化數(shù)量變化的影響。一階線性常系數(shù)差分方程的解、平衡點(diǎn)及其穩(wěn)定性一階線性常系數(shù)差分方程的解、平衡點(diǎn)及其穩(wěn)定性 自然環(huán)境下,b=0 人工孵化條件下 令xk=

6、xk+1=x得 差分方程的平衡點(diǎn) k時(shí)，xkx,稱平衡點(diǎn)是穩(wěn)定的0kkxax10(1)kkkxa xbaa011kkaaxba1bxa1kkxa xb高階線性常系數(shù)差分方程 如果第k+1時(shí)段變量Xk+1不僅取決于第k時(shí)段變量Xk，而且與以前時(shí)段變量有關(guān)，就要用高階差分方程來描述一年生植物的繁殖 一年生植物春季發(fā)芽，夏天開花，秋季產(chǎn)種，沒有腐爛，風(fēng)干，被人為掠取的那些種子可以活過冬天，其中一部分能在第2年春季發(fā)芽，然后開花，產(chǎn)種，其中的另一部分雖未能發(fā)芽，但如又能活過一個(gè)冬天，則其中一部分可在第三年春季發(fā)芽，然后開花，產(chǎn)種，如此繼續(xù)，一年生植物只能活1年，而近似的認(rèn)為，種子最多可以活過兩個(gè)冬天，

7、試建立數(shù)學(xué)模型研究這種植物數(shù)量變化的規(guī)律，及它能一直繁殖下去的條件。模型及其求解 記一棵植物春季產(chǎn)種的平均數(shù)為c,種子能活過一個(gè)冬天的(1歲種子)比例為b,活過一個(gè)冬天沒有發(fā)芽又活過一個(gè)冬天的（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 實(shí)際上

8、，就是Xk= pXk-1 + qXk-2 我們需要知道x0,a1,a2,c, 考察b不同時(shí)，種子繁殖的情況。在這里假設(shè) X0=100,a1=0.5,a2=0.25,c=10,b=0.180.20 這樣可以用matlab計(jì)算了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*(xk-1)+q*(xk-2); end K=(0:20); Y1=zwfz(100,21,0.1

9、8); Y2=zwfz(100,21,0.19); Y3=zwfz(100,21,0,20); Round(k,y1,y2,y3) Plot(k,y1,k,y2,:,k,y3,o), Gtext(b=0.18),gtext(b=0.19),gtext(b=0.20)結(jié)果分析：Xk= pXk-1 + qXk-2 (1) x1+px0=0 (2) 對高階差分方程可以尋求形如的解。代入(1)式得稱為差分方程的特征方程。差分方程的特征根：方程(1)的解可以表為C1,c2 由初始條件x0,x1確定。kkx20pq21 , 242ppq1122kkkxcc1,21,0()kxk 本例中，用待定系數(shù)的方法可

10、以求出b=0.18時(shí),c1=95.64, c2=4.36 ， 這樣實(shí)際上，植物能一直繁殖下去的條件是b0.19112( ,)(0.9430, 0.0430) 95.64(0.9430)4.36( 0.0430)kkkx 1,21,()kxk 1 , 251 02b線性常系數(shù)差分方程組 汽車租賃公司的運(yùn)營一家汽車租賃公司在3個(gè)相鄰的城市運(yùn)營，為方便顧客起見公司承諾，在一個(gè)城市租賃的汽車可以在任意一個(gè)城市歸還。根據(jù)經(jīng)驗(yàn)估計(jì)和市場調(diào)查，一個(gè)租賃期內(nèi)在A市租賃的汽車在A,B,C市歸還的比例分別為0.6,0.3,0.1;在B市租賃的汽車歸還比例0.2,0.7,0.1;C市租賃的歸還比例分別為0.1,0.

11、3,0.6。若公司開業(yè)時(shí)將600輛汽車平均分配到3個(gè)城市，建立運(yùn)營過程中汽車數(shù)量在3個(gè)城市間轉(zhuǎn)移的模型，并討論時(shí)間充分長以后的變化趨勢。0.60.3 A B C A B C A B C假設(shè)在每個(gè)租賃期開始能把汽車都租出去，并都在租賃期末歸還0.10.70.20.10.60.30.1模型及其求解 記第k個(gè)租賃期末公司在ABC市的汽車數(shù)量分別為x1(k),x2(k),x3(k)（也是第k+1個(gè)租賃期開始各個(gè)城市租出去的汽車數(shù)量），很容易寫出第k+1個(gè)租賃期末公司在ABC市的汽車數(shù)量為(k=0,1,2,3)112321233123(1)0.6( )0.2( )0.1( )(1)0.3 ( )0.7(

12、 )0.3( )(1)0.1 ( )0.1( )0.6( )x kx kx kx kx kx kx kx kx kx kx kx k 用矩陣表示用matlab編程，計(jì)算x(k),觀察n年以后的3個(gè)城市的汽車數(shù)量變化情況112233(1)0.60.20.1( )(1)0.30.70.3( )(1)0.10.10.6( )x kx kx kx kx kx k function x=czqc(n) A=0.6,0.2,0.1;0.3,0.7,0.3;0.1,0.1,0.6; x(:,1)=200,200,200; for k=1:n x(:,k+1)=A*x(:,k); end如果直接看10年或者2

13、0年發(fā)展趨勢，可以直接在命令窗口（commond window）作，而不是必須編一個(gè)函數(shù)112233(1)0.60.20.1()(1)0.30.70.3()(1)0.10.10.6()xkxkxkxkxkxk A=0.6,0.2,0.1;0.3,0.7,0.3;0.1,0.1,0.6; n=10; for k=1:n x(:,1)=200,200,200; x(:,k+1)=A*x(:,k); end round(x)作圖觀察數(shù)量變化趨勢012345678910120140160180200220240260280300 x1(k)x2(k)x3(k) k=0:10; plot(k,x) ，g

14、rid gtext(x1(k), gtext(x2(k), gtext(x3(k) 可以看到時(shí)間充分長以后3個(gè)城市汽車數(shù)量趨于180，300，120 可以考察這個(gè)結(jié)果與初始條件是否有關(guān) 若最開始600輛汽車都在A市，可以看到變化時(shí)間充分長以后，各城市汽車數(shù)量趨于穩(wěn)定，與初始值無關(guān)直接輸入x（:,1）的值即可x(:,1)=600,0,0; round(x);plot(k,x),grid0123456789100100200300400500600按年齡分組的種群增長 野生或飼養(yǎng)的動(dòng)物因繁殖而增加，因自然死亡和人為屠殺而減少，不同年齡動(dòng)物的繁殖率，死亡率有較大差別，因此在研究某一種群數(shù)量的變化時(shí)，

15、需要考慮年齡分組的種群增長。 將種群按年齡等間隔的分成若干個(gè)年齡組，時(shí)間也離散化為時(shí)段，給定各年齡組種群的繁殖率和死亡率，建立按年齡分組的種群增長模型，預(yù)測未來各年齡組的種群數(shù)量，并討論時(shí)間充分長以后的變化趨勢。模型及其求解 設(shè)種群按年齡等間隔的分成n個(gè)年齡組，記i=1,2,，n,時(shí)段記作k=0,1,2,且年齡組區(qū)間與時(shí)段長度相等(若5歲為一個(gè)年齡組，則5年為一個(gè)時(shí)段)。以雌性個(gè)體為研究對象 記在時(shí)段k第i年齡組的數(shù)量為xi(k);第i年齡組的繁殖率為bi，表示每個(gè)個(gè)體在一個(gè)時(shí)段內(nèi)繁殖的數(shù)量；第i年齡組死亡率為di，表示一個(gè)時(shí)段內(nèi)死亡數(shù)與總數(shù)的比，si=1-di是存活率。 注意：第k時(shí)段的第i

16、年齡組活過來的，是第k+1時(shí)段的第i+1年齡組 Xi+1(k+1)=sixi(k) i=1，2,n-1, k=0,1, 各年齡組在第k時(shí)段繁殖的數(shù)量和是第k+1時(shí)段的第1年齡組 X1(k+1)= k=0,1, 記在時(shí)段k種群各年齡組的數(shù)量為 X(k)=x1(k),x2(k),xn(k)1()niiibxk121121000000000nnnbbbbsLss 這樣，有x(k+1)=Lx(k),k=0,1, 給定在0時(shí)段，各年齡組的初始數(shù)量x(0) 就可以預(yù)測任意時(shí)段k,各年齡組的數(shù)量 設(shè)一種群分成5個(gè)年齡組， 繁殖率b1=0,b2=0.2,b3=1.8,b4=0.8,b5=0.2 存活率s1=0.5,s2=0.8,s3=0.8,s4=0.1 各年齡組現(xiàn)有數(shù)量都是100只， 用matlab計(jì)算x(k) b=0,0.2,1.8,0.8,0.2; s=diag(0.5,0.8,0.8,0.1); L=b;s,zeros(4,1); x(: