實驗二微分方程與差分方程模型Matlab求解

1、實驗二：微分方程與差分方程模型Matlab 求解一、實驗?zāi)康? 掌握解析、數(shù)值解法，并學(xué)會用圖形觀察解的形態(tài)和進行解的定性分析；2 熟悉 MATLAB軟件關(guān)于微分方程求解的各種命令；3 通過范例學(xué)習(xí)建立微分方程方面的數(shù)學(xué)模型以及求解全過程；4 熟悉離散 Logistic 模型的求解與混沌的產(chǎn)生過程。二、實驗原理1. 微分方程模型與 MATLAB 求解解析解用 MATLAB命令 dsolve( eqn1 , eqn2, .)求常微分方程（組）的解析解。其中 eqni' 表示第 i 個微分方程， Dny 表示 y 的 n 階導(dǎo)數(shù) , 默認的自變量為 t 。（ 1） 微分方程例 1 求解一階

2、微分方程 dy 1 y 2 dx(1) 求通解輸入：dsolve('Dy=1+y2')輸出：ans =tan(t+C1)（ 2）求特解輸入：dsolve('Dy=1+y2','y(0)=1','x')指定初值為 1，自變量為 x輸出：ans =tan(x+1/4*pi)x2 yxy (x21 ) y 0例 24求解二階微分方程y(/ 2)2y (/ 2)2 /原方程兩邊都除以 x2，得 y1y (112 ) y 0x4x輸入 :dsolve('D2y+(1/x)*Dy+(1-1/4/x2)*y=0','y(

3、pi/2)=2,Dy(pi/2)=-2/pi','x')ans =- (exp(x*i)*(pi/2)(1/2)*i)/x(1/2) + (exp(x*i)*exp(-x*2*i)*(pi/2)(3/2)*2*i)/(pi*x(1/2)試試能不用用 simplify函數(shù)化簡輸入 : simplify(ans)ans =2(1/2)*pi(1/2)/x(1/2)*sin(x)（ 2）微分方程組例 3 求解df/dx=3f+4g;dg/dx=-4f+3g。（1）通解 :f,g=dsolve('Df=3*f+4*g','Dg=-4*f+3*g'

4、)f =exp(3*t)*(C1*sin(4*t)+C2*cos(4*t)g =exp(3*t)*(C1*cos(4*t)-C2*sin(4*t)特解：f,g=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1')f =exp(3*t)*sin(4*t)g =exp(3*t)*cos(4*t)數(shù)值解在微分方程 ( 組) 難以獲得解析解的情況下，可以用 Matlab 方便地求出數(shù)值解。格式為 :t,y = ode23('F',ts,y0,options)注意：微分方程的形式： y'

5、; = F(t, y)，t 為自變量， y 為因變量（可以是多個，如微分方程組）；t, y 為輸出矩陣，分別表示自變量和因變量的取值；F 代表一階微分方程組的函數(shù)名（ m文件，必須返回一個列向量，每個元素對應(yīng)每個方程的右端）；ts 的取法有幾種，（ 1）ts=t0, tf表示自變量的取值范圍，（2）ts=t0,t1,t2,tf,則輸出在指定時刻t0,t1,t2,tf處給出，（3）ts=t0:k:tf,則輸出在區(qū)間 t0,tf的等分點給出；y0 為初值條件；options用于設(shè)定誤差限（缺省是設(shè)定相對誤差是10(-3) ，絕對誤差是 10(-6) ）；ode23 是微分方程組數(shù)值解的低階方法，o

6、de45 為中階方法，與 ode23 類似。例 4求解一個經(jīng)典的范得波（Van Der pol ）微分方程：u'' (u21)u' u 0，u(0) 1, u' (0) 0解 形式轉(zhuǎn)化：令 y1u(t);y2u'(t ) 。則以上方程轉(zhuǎn)化一階微分方程組：y1 y2 ; y2 (1 y12 ) y2y1 。編寫 M文件如下，必須是 M文件表示微分方程組，并保存，一般地，M文件的名字與函數(shù)名相同，保存位置可以為默認的work 子目錄，也可以保存在自定義文件夾，這時注意要增加搜索路徑（ FileSet PathAdd Folder）functiondot1=v

7、dpol(t,y);dot1=y(2); (1-y(1)2)*y(2)-y(1);在命令窗口寫如下命令 :t,y=ode23('vdpol',0,20,1,0);y1=y(:,1);y2=y(:,2);plot(t,y1,t,y2,'-');title('Van Der Pol Solution ');xlabel('Time,Second');ylabel('y(1)andy(2)')執(zhí)行：Van Der Pol Solution321)2(ydn0a)1(y-1-2-302468101214161820Time

8、,Second注： Van der Pol 方程描述具有一個非線性振動項的振動子的運動過程。最初 , 由于它在非線性電路上的應(yīng)用而引起廣泛興趣。一般形式為u'' (u 21)u' u0 。圖形解無論是解析解還是數(shù)值解， 都不如圖形解直觀明了。 即使是在得到了解析解或數(shù)值解的情況下，作出解的圖形，仍然是一件深受歡迎的事。這些都可以用Matlab 方便地進行。(1) 圖示解析解如果微分方程（組）的解析解為： y=f (x)，則可以用 Matlab 函數(shù) fplot 作出其圖形：fplot('fun',lims)其中： fun給出函數(shù)表達式；lims=xmin

9、 xmax ymin ymax限定坐標(biāo)軸的大小。例如fplot('sin(1/x)', 0.01 0.1 -1 1)10.80.60.40.20-0.2-0.4-0.6-0.8-10.020.030.040.050.060.070.080.090.10.01(2) 圖示數(shù)值解設(shè)想已經(jīng)得到微分方程 （組）的數(shù)值解 (x,y) ?？梢杂?Matlab 函數(shù) plot(x,y)直接作出圖形。其中x 和 y 為向量（或矩陣）。2、Volterra 模型（食餌捕食者模型）意大利生物學(xué)家 Ancona 曾致力于魚類種群相互制約關(guān)系的研究， 他從第一次世界大戰(zhàn)期間 ,地中海各港口捕獲的幾種魚

10、類捕獲量百分比的資料中，發(fā)現(xiàn)鯊魚的比例有明顯增加（見下表） 。年代19141915191619171918百分比11.921.422.121.236.4年代19191920192119221923百分比27.316.015.914.819.7戰(zhàn)爭為什么使鯊魚數(shù)量增加？是什么原因？因為戰(zhàn)爭使捕魚量下降，食用魚增加，顯然鯊魚也隨之增加。但為何鯊魚的比例大幅增加呢？生物學(xué)家 Ancona 無法解釋這個現(xiàn)象， 于是求助于著名的意大利數(shù)學(xué)家 V.Volterra ，希望建立一個 食餌 捕食者系統(tǒng) 的數(shù)學(xué)模型，定量地回答這個問題。1 、符號說明：x1(t), x2(t)分別是食餌、捕食者（鯊魚）在t 時刻

11、的數(shù)量； r1, r2 是食餌、捕食者的固有增長率； 1 是捕食者掠取食餌的能力， 2 是食餌對捕食者的供養(yǎng)能力；2、基本假設(shè)：捕食者的存在使食餌的增長率降低，假設(shè)降低的程度與捕食者數(shù)量成正比，dx1x1( r11 x2 )即 dt食餌對捕食者的數(shù)量x2 起到增長的作用，其程度與食餌數(shù)量 x1 成正比，即dx2x2(r22 x1)dt綜合以上和，得到如下模型：模型一： 不考慮人工捕獲的情況dx1x1(r11 x2 )dtdx2x2( r22 x1 )dt該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關(guān)系，沒有考慮食餌和捕食者自身的阻滯作用，是 Volterra提出的最簡單的模型。

12、給定一組具體數(shù)據(jù)，用matlab 軟件求解。食餌： r1= 1,1= 0.1, x10= 25;捕食者 ( 鯊魚 ) ：r2=0.5, 2=0.02, x20= 2;dx1x1(10.1x2 )dtx1 (0) 25, x2 (0) 2dx2x2(0.5 0.02x1)dt編制程序如下1、建立 m-文件 shier.m如下：function dx=shier(t,x)dx=zeros(2,1); %初始化dx(1)=x(1)*(1-0.1*x(2);dx(2)=x(2)*(-0.5+0.02*x(1);2、在命令窗口執(zhí)行如下程序：t,x=ode45('shier',0:0.1:

13、15,25 2);plot(t,x(:,1),'-',t,x(:,2),'*'),grid1009080706050403020100051015圖中，藍色曲線和綠色曲線分別是食餌和鯊魚數(shù)量隨時間的變化情況，量都呈現(xiàn)出周期性，而且鯊魚數(shù)量的高峰期稍滯后于食餌數(shù)量的高峰期。從圖中可以看出它們的數(shù)畫出相軌跡圖：plot(x(:,1),x(:,2)3025201510500102030405060708090100模型二考慮人工捕獲的情況假設(shè)人工捕獲能力系數(shù)為e，相當(dāng)于食餌的自然增長率由r 1 降為 r 1-e ，捕食者的死亡率由 r 2 增為 r 2+e，因此模型一

14、修改為：設(shè)戰(zhàn)前捕獲能力系數(shù) e=0.3,戰(zhàn)爭中降為 e=0.1,其它參數(shù)與模型 (一 )的參dx1x1( r1e)1x2 dtdx2x2 (r2e)2 x1 dt數(shù)相同。觀察結(jié)果會如何變化？dx1x1(0.70.1x2 )dtdx2x2( 0.80.02x1)dtx1 (0) 25, x2 (0)21）當(dāng) e=0.3 時：2 ）當(dāng) e=0.1 時：dy1y1(0.90.1y2 )dtdy2y2( 0.60.02 y1 )dty1 (0) 25, y2 (0)2分別求出兩種情況下鯊魚在魚類中所占的比例。即計算p1x2 (t);p2y2 (t)(t)(t)x1 (t) x2 (t)y1 (t )

15、y2 (t)畫曲線： plot(t,p1(t),t,p2(t),'*')MATLAB 編程實現(xiàn)建立兩個 M文件functiondx=shier1(t,x)dx=zeros(2,1);dx(1)=x(1)*(0.7-0.1*x(2);dx(2)=x(2)*(-0.8+0.02*x(1);functiondy=shier2(t,y)dy=zeros(2,1);dy(1)=y(1)*(0.9-0.1*y(2);dy(2)=y(2)*(-0.6+0.02*y(1);運行以下程序：t1,x=ode45('shier1',0 15,25 2);t2,y=ode45('

16、;shier2',0 15,25 2);x1=x(:,1);x2=x(:,2);p1=x2./(x1+x2);y1=y(:,1);y2=y(:,2);p2=y2./(y1+y2);plot(t1,p1,'-',t2,p2,'*')0.80.70.60.50.40.30.20.10510150圖中 * 曲線為戰(zhàn)爭中鯊魚所占比例。結(jié)論：戰(zhàn)爭中鯊魚的比例比戰(zhàn)前高。3、 Logistic映射logistic 映射 - 通向混沌的道路混沌系統(tǒng)，由于其行為的復(fù)雜性，往往認為其動態(tài)特性（運動方程）也一定非常復(fù)雜，事實并非如此， 一個參量很少、 動態(tài)特性非常簡單的系統(tǒng)有

17、時也能夠產(chǎn)生混沌現(xiàn)象，以一維蟲口模型為例，假設(shè)某一區(qū)域內(nèi)的現(xiàn)有蟲口數(shù)為y ，昆蟲的繁殖率為r ，且第 n 代昆蟲不能n存活于第 n+1 代，既無世代交疊，則第n+1 代蟲口數(shù)為 yn 1 ryn ，r 1 時，蟲口會無限制地增長； r 1 時，蟲口最終會趨于消亡，因此需要對模型進行修正。由于環(huán)境的制約和食物有限，因爭奪生存空間發(fā)生相互咬斗事件的最大次數(shù)為1 yn ( yn1) ，即制約蟲口數(shù)的因素與 y2 成正比，設(shè)咬斗事件的戰(zhàn)死率為2y2 ，則有： 則對蟲口的修正項為nnyn 1rynyn2.令 xnyn ，則rxn 1 rxn (1 xn )（ 1）取最大蟲口數(shù)為1，且蟲口數(shù)不能為負，則；

18、當(dāng)=0.5 時，方程有極大值， xn 11 r 而又必須小于1，因而 r 4，則參量 r 的取值范圍為1 到 4，這就得到4一個抽象的標(biāo)準(zhǔn)蟲口方程（1）。記映射 f 為f (x) rx (1 x)（ 2）方程（ 1）可寫為xn 1f ( xn )（3）這一迭代關(guān)系通常稱為 logistic映射。從 0,1內(nèi)點 x 出發(fā)，由 Logistic映射的0迭代形成xn= f n( x0),n = 0,1,2,序列 xn 稱為 x0 的軌道。一個看似簡單的系統(tǒng)，隨著參量的不同會表現(xiàn)出截然不同的行為，當(dāng)r 的取值范圍在 13 時，方程（ 1）有定態(tài)解x11r即方程通過多次迭代后趨于一個穩(wěn)定的不動點，此時系

19、統(tǒng)是穩(wěn)定的。x11為方程f ( x)x 的解，稱為周期2 點。r當(dāng) r在33.448范圍內(nèi)取值時，經(jīng)過多次迭代，方程（1）出現(xiàn)周期2 點x1 和x2 ，即x1和 x2 是方程 f2 ( x)x 的解，滿足x1rx2 (1x2 )x2rx1 (1x1 )r 3是使解有意義的 r 最小值。隨著 r 的增大， r=3.449 ； 3.544 ；3.564 依次出現(xiàn)周期4、周期 8、周期 16 的振蕩解， r 的極限值約為 3.569 。這種行為稱為倍周期分岔，直到r 3.5699 時，系統(tǒng)進入了混沌狀態(tài)，如下圖所示，此時系統(tǒng)的狀態(tài)不再具有規(guī)律性，而是發(fā)生隨機的波動，使圖d的右側(cè)的大部分區(qū)域被涂黑了，

20、仔細觀察發(fā)現(xiàn)， 混沌區(qū)域并非一片涂斑， 而是有粗粗細細的白色“豎線”，稱為周期窗口，隨著參量r 的增大（如 r 18 ）時，混沌突然消失，系統(tǒng)出現(xiàn)周期三的穩(wěn)定狀態(tài)，Logistic映射分岔圖接著倍周期分岔以更快的速度進行， 再次進入混沌狀態(tài)。 如果將周期窗口放大， 發(fā)現(xiàn)其結(jié)構(gòu)與分岔圖的整體結(jié)構(gòu)具有相似性，而且是一種無限嵌套的自相似結(jié)構(gòu)?？梢钥闯觯?通過改變系統(tǒng)參量，使系統(tǒng)進入混沌的第一種模式是倍周期分岔，即由不動點周期二周期四 無限倍周期進入混沌狀態(tài)。當(dāng)然通向混沌的道路不只于此，第二種通向的道路是：從平衡態(tài)到周期運動，再到擬周期運動，直到進入混沌狀態(tài)。第三種通向混沌的方式是陣發(fā)（或間歇）道路，

21、即系統(tǒng)在近似周期運動的過程中，通過改變參量，系統(tǒng)會出現(xiàn)陣發(fā)性混沌過程， 隨著參量的調(diào)整， 陣發(fā)性混沌越來越頻繁， 近似的周期運動越來越少，最后進入混沌。Matlab 程序下面程序繪制r=2 ， x0=0.3 的軌道clearall;clf;x=0.3;r=2;n=input( 'Please input a number: ' ); for i=1:nx1=r*x*(1-x);x=x1;plot(i,x1,'.');holdon ;end下面程序繪制 logistic 映射 r=3,5 的軌道，觀察是否有周期點clearall;clf;x=0.3;r=3.5;f

22、ori=1:100x1=r*x*(1-x);x=x1;0.510.50.490.480.470.460.450.440.430.42010203040506070809010010.90.80.70.60.50.40102030405060708090100plot(i,x1,'.');holdon ;end下面程序繪制 logistic 映射 r=4 的軌道，觀察其混沌clearall;clf;x=0.007;fori=1:100x1=4*x*(1-x);x=x1;plot(i,x1,'.');holdon ;end10.90.80.70.60.50.40.3

23、0.20.100102030405060708090100下面程序繪制在混沌狀態(tài)下不同初值x0=0.100001 和 x0=0.1 的軌道的差（對初值的敏感性）clearall;clf;x1=0.100001;x11=0.1;fori=1:1000x2=4*x1*(1-x1);x1=x2;x22=4*x11*(1-x11);x11=x22;plot(i,x11-x1);holdon ;end下面程序繪制logistic 映射的分岔圖clearall;clf;x=0.2;10.80.60.40.20-0.2-0.4-0.6-0.8-10100200300400500600700800900100

24、0forr=1:0.001:4fori=1:18x1=r*x*(1-x);x=x1;ifi>9plot(r,x);holdon ;endendend10.90.80.70.60.50.40.30.20.1011.522.533.54蛛網(wǎng)迭代clearall;clf;x=0.007;y=0;fori=1:200XY=x,y;y=4*x*(1-x);XY1=x,y;line(XY,XY1);holdon ;x=y;XY2=x,x;line(XY1,XY2);holdon ;end10.90.80.70.60.50.40.30.20.1000.10.20.30.40.50.60.70.80.9

25、1三、實驗內(nèi)容1求微分方程的解析解 ,并畫出圖形，y =y+2 x , y(0)1, 0x12求微分方程的數(shù)值解 ,并畫出圖形，yy cos x0, y(0)1, y (0)03兩種相似的群體之間為了爭奪有限的同一種食物來源和生活空間而進行生存競爭時， 往往是競爭力較弱的種群滅亡， 而競爭力較強的種群達到環(huán)境容許的最大數(shù)量。假設(shè)有甲、乙兩個生物種群，當(dāng)它們各自生存于一個自然環(huán)境中，均服從 Logistic 規(guī)律。（1） x1 (t), x2 (t) 是兩個種群的數(shù)量；（2） r1 ,r2 是它們的固有增長率；（3） n1, n2 是它們的最大容量；（4） m2 (m1 ) 為種群乙 ( 甲 ) 占據(jù)甲 ( 乙) 的位置的數(shù)量，并且 m2 x2 ; m1 x1.計算 x1 (t ), x2 (t ) ,畫出圖形及相軌跡圖。解釋其解變化過程。2）設(shè)r1r21, n1n2100, x10x2010 ,=1.5 ，=0.7 ，計算 x1 (t), x2 (t ) ,畫出圖形及相軌跡圖。解釋其解變化過程。4.繪制 Logistic映射的軌道圖 , 分岔圖和蛛網(wǎng)迭代圖。四、實驗心得 其中專業(yè)理論知識內(nèi)容包括：保安理論知識、消防業(yè)務(wù)知識、職業(yè)道德、法律常識、保安禮儀、救護知識。作技