MATLAB求差分方程模型課件

用Matlab求解差分方程問題一階線性常系數(shù)差分方程高階線性常系數(shù)差分方程線性常系數(shù)差分方程組一階線性常系數(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······得Xn＝(1+r)^nX0

已知X0=100,在較好，中等和較差的自然環(huán)境下r=0.0194,-0.0324和-0.0382我們利用Matlab編程，遞推20年后觀察沙丘鶴的數(shù)量變化情況Matlab實現(xiàn)Xn＝(1+r)^nX01.在指令窗口輸入：>>X0=100;>>r=0.0194;>>n=20;>>Xn=(1+r)^n*X0執(zhí)行后得到：Xn=146.85633.建立M函數(shù)文件格式：function輸出變量＝函數(shù)名稱（輸入變量）例如：輸入functionXn=sqh(n,r)X0=100;Xn=(1+r)^n*X0;end單擊工具欄中的保存按鈕，文件名sqh2.m保存以后在指令窗口調(diào)用sqh2函數(shù),輸入：>>Xn＝sqh2(20,0.0194)輸出同上。functionXn=sqh(n,r,X0)Xn=(1+r)^n*X0;end單擊工具欄中的保存按鈕，文件名sqh3.m保存以后在指令窗口調(diào)用sqh3函數(shù),輸入：>>Xn＝sqh3(20,0.0194,100)輸出同上。Matlab實現(xiàn)Xk+1＝(1+r)Xk

functionx=sqh(n,r)x(1)=100;fork=1:nx(k+1)=(1+r)*x(k);end單擊工具欄中的保存按鈕，文件名sqh4.m保存以后在指令窗口調(diào)用sqh4函數(shù),輸入：>>xn＝sqh4(20,0.0194)輸出同上。xn=Columns1through6100.0000101.9400103.9176105.9336107.9888110.0837Columns7through12112.2194114.3964116.6157118.8780121.1843123.5353Columns13through18125.9318128.3749130.8654133.4042135.9922138.6305Columns19through21141.3199144.0615146.8563在同一坐標系下畫圖

>>k=0:20;%一個行向量

>>y1=sqh4(20,0.0194);>>y2=sqh4(20,-0.0324);>>y3=sqh4(20,-0.0382);>>plot(k,y1,k,y2,':',k,y3,'r')人工孵化是挽救瀕危物種的措施之一，如果每年孵化a只鶴放入保護區(qū)，觀察在自然條件下沙丘鶴的數(shù)量如何變化Xk+1=（1+r）Xk+a在M文件中輸入：functionX=rgsqh(n,r,a)X(1)=100;fork=1:nX(k+1)=（1＋r)*X(k)+a;end設(shè)在中等自然條件下每年孵化5只鶴放入保護區(qū).>>k=0:20;%一個行向量>>y=rgsqh(20,-0.0324,5);>>plot(k,y)高階線性常系數(shù)差分方程

——一年生植物的繁殖一年生植物春季發(fā)芽，夏天開花，秋季產(chǎn)種，沒有腐爛，風干，被人為掠取的那些種子可以活過冬天，其中一部分能在第2年春季發(fā)芽，然后開花，產(chǎn)種，其中的另一部分雖未能發(fā)芽，但如又能活過一個冬天，則其中一部分可在第三年春季發(fā)芽，然后開花，產(chǎn)種，如此繼續(xù)，一年生植物只能活1年，而近似的認為，種子最多可以活過兩個冬天，試建立數(shù)學模型研究這種植物數(shù)量變化的規(guī)律，及它能一直繁殖下去的條件。實際上，就是Xk=pXk-1+qXk-2我們需要知道x0,a1,a2,c,考察b不同時，種子繁殖的情況。在這里假設(shè)X0=100,a1=0.5,a2=0.25,c=10,b=0.18~0.20這樣可以用matlab計算了Xk=a1bcXk-1+a2b(1-a1)bcXk-2Xk=a1bcXk-1+a2b(1-a1)bcXk-2functionx=zwfz(x0,n,b)c=10;a1=0.5;a2=0.25;p=a1*b*c;q=a2*b*(1-a1)*b*c;x(1)=x0;x(2)=p*x(1);fork=3:nx(k)=p*x(k-1)+q*x(k-2);end假設(shè)X0=100,a1=0.5,a2=0.25,c=10,b=0.18~0.20簡化為Xk=pXk-1+qXk-2>>k=0:20;>>y1=zwfz(100,21,0.18);>>y2=zwfz(100,21,0.19);>>y3=zwfz(100,21,0.20);>>plot(k,y1,k,y2,':',k,y3,'r')結(jié)果分析：Xk＋pXk-1+qXk-2＝0(1)

x1+px0=0(2)

對高階差分方程可以尋求形如的解。代入(1)式得稱為差分方程的特征方程。差分方程的特征根：方程(1)的解可以表為C1,C2由初始條件x0,x1確定。本例中，用待定系數(shù)的方法可以求出b=0.18時,c1=95.64,c2=4.36，這樣實際上，植物能一直繁殖下去的條件是b>0.191插入的數(shù)的大小，其必須滿足一個遞減趨勢，增長值為負模型及其求解記第k個租賃期末公司在ABC市的汽車數(shù)量分別為x1(k),x2(k),x3(k)（也是第k+1個租賃期開始各個城市租出去的汽車數(shù)量），很容易寫出第k+1個租賃期末公司在ABC市的汽車數(shù)量為(k=0,1,2,3···)用矩陣表示觀察n年以后的3個城市的汽車數(shù)量變化情況1.直接在命令窗口（commondwindow）作>>A=[0.6,0.2,0.1;0.3,0.7,0.3;0.1,0.1,0.6];>>n=10;>>fork=1:nx(:,1)=[200,200,200]';x(:,k+1)=A*x(:,k);end>>x(:,k+1)輸出：ans=179.9324299.9895120.07812.建立M文件functionx=qczl(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]';fork=1:nx(:,k+1)=A*x(:,k);end作圖觀察5年以后數(shù)量的變化趨勢：在指令窗口輸入>>y1=qczl(5);>>k=0:5;>>plot(k,y1)練習：若最初600輛汽車平均分配在三個城市，作圖觀察30年以后數(shù)量的變化趨勢若最開始600輛汽車都在A市，變化時間充分長以后，各城市汽車數(shù)量的變化情況如何？可以看到時間充分長以后3個城市汽車數(shù)量趨于180，300，120若最開始600輛汽車都在A市，可以看到變化時間充分長以后，各城市汽車數(shù)量趨于穩(wěn)定，與初始值無關(guān)利用線性規(guī)劃知識；找到凸集，因為所謂的最優(yōu)解，一般都在凸集按年齡分組的種群增長野生或飼養(yǎng)的動物因繁殖而增加，因自然死亡和人為屠殺而減少，不同年齡動物的繁殖率，死亡率有較大差別，因此在研究某一種群數(shù)量的變化時，需要考慮年齡分組的種群增長。將種群按年齡等間隔的分成若干個年齡組，時間也離散化為時段，給定各年齡組種群的繁殖率和死亡率，建立按年齡分組的種群增長模型，預(yù)測未來各年齡組的種群數(shù)量，并討論時間充分長以后的變化趨勢。模型及其求解設(shè)種群按年齡等間隔的分成n個年齡組，記i=1,2,···，n,時段記作k=0,1,2···,且年齡組區(qū)間與時段長度相等(若5歲為一個年齡組，則5年為一個時段)。以雌性個體為研究對象.記在時段k第i年齡組的數(shù)量為Xi(k);第i年齡組的繁殖率為bi，表示每個個體在一個時段內(nèi)繁殖的數(shù)量；第i年齡組死亡率為Si，表示一個時段內(nèi)死亡數(shù)與總數(shù)的比，Ei=1-Si是存活率。注意：第k時段的第i年齡組活過來的數(shù)量，是第k+1時段第i+1年齡組的數(shù)量

Xi+1(k+1)=EiXi(k)i=1，2,···,n-1,k=0,1,····各年齡組在第k時段繁殖的數(shù)量和是第k+1時段的第1年齡組的數(shù)量

X1(k+1)=k=0,1,····記在時段k種群各年齡組的數(shù)量為

X(k)=[X1(k),X2(k),····,Xn(k)]’得X(k+1)=LX(k),k=0,1,····設(shè)一種群分成5個年齡組，繁殖率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計算X(k)提示：L矩陣的輸入>>b=[0,0.2,1.8,0.8,0.2];>>E=diag([0.5,0.8,0.8,0.1]);>>L=[b;E,zeros(4,1)];>>b=[0,0.2,1.8,0.8,0.2];>>E=diag([0.5,0.8,0.8,0.1