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

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

2、 (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) 例2 求解二階微分方程 原方程兩邊都除以，得輸入:dsolve('D2y+(1/x)*Dy+(1-1/4/x2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x') ans = - (exp(x*i)*(pi/2)(1/2)

3、*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') f =exp(3*t)*(C1*sin(

4、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' = F(t, y)，t為自變量，

5、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ù)值解的低階

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

7、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方程描述具有一個非線性振動項的振動子的運動過程。最初,由于它在非線性電路上的應(yīng)用而引起廣泛興趣。一般形式為。圖形解無論是解析解還是數(shù)值解，都不如圖形解直觀明了。即使是在得到了解析解或數(shù)值解的情況下，作出

8、解的圖形，仍然是一件深受歡迎的事。這些都可以用Matlab方便地進行。(1)圖示解析解如果微分方程（組）的解析解為：y=f (x)，則可以用Matlab函數(shù)fplot作出其圖形：fplot('fun',lims)其中：fun給出函數(shù)表達式；lims=xmin xmax ymin ymax限定坐標軸的大小。例如fplot('sin(1/x)', 0.01 0.1 -1 1) (2)圖示數(shù)值解設(shè)想已經(jīng)得到微分方程（組）的數(shù)值解(x,y)。可以用Matlab函數(shù)plot(x,y)直接作出圖形。其中x和y為向量（或矩陣）。2、Volterra模型（食餌捕食者模型）意大利

9、生物學家Ancona曾致力于魚類種群相互制約關(guān)系的研究，他從第一次世界大戰(zhàn)期間,地中海各港口捕獲的幾種魚類捕獲量百分比的資料中，發(fā)現(xiàn)鯊魚的比例有明顯增加（見下表）。年代19141915191619171918百分比11.921.422.121.236.4年代19191920192119221923百分比27.316.015.914.819.7戰(zhàn)爭為什么使鯊魚數(shù)量增加？是什么原因？因為戰(zhàn)爭使捕魚量下降，食用魚增加，顯然鯊魚也隨之增加。 但為何鯊魚的比例大幅增加呢？生物學家Ancona無法解釋這個現(xiàn)象，于是求助于著名的意大利數(shù)學家V.Volterra，希望建立一個食餌捕食者系統(tǒng)的數(shù)學模型，定量地回

10、答這個問題。 1、符號說明：x1(t), x2(t)分別是食餌、捕食者（鯊魚）在t時刻的數(shù)量； r1, r2是食餌、捕食者的固有增長率；1是捕食者掠取食餌的能力， 2是食餌對捕食者的供養(yǎng)能力；2、基本假設(shè)： 捕食者的存在使食餌的增長率降低，假設(shè)降低的程度與捕食者數(shù)量成正比，即食餌對捕食者的數(shù)量x2起到增長的作用， 其程度與食餌數(shù)量x1成正比，即綜合以上和，得到如下模型：模型一：不考慮人工捕獲的情況 該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關(guān)系，沒有考慮食餌和捕食者自身的阻滯作用，是Volterra提出的最簡單的模型。給定一組具體數(shù)據(jù)，用matlab軟件求解。 食餌： r1=

11、 1, 1= 0.1, x10= 25; 捕食者(鯊魚)：r2=0.5, 2=0.02, x20= 2;編制程序如下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:15,25 2); plot(t,x(:,1),'-',t,x(:,2),'*'),grid 圖中，藍色曲線和綠色曲線分別是食餌和鯊

12、魚數(shù)量隨時間的變化情況，從圖中可以看出它們的數(shù)量都呈現(xiàn)出周期性，而且鯊魚數(shù)量的高峰期稍滯后于食餌數(shù)量的高峰期。畫出相軌跡圖：plot(x(:,1),x(:,2) 模型二 考慮人工捕獲的情況假設(shè)人工捕獲能力系數(shù)為e，相當于食餌的自然增長率由r1 降為r1-e，捕食者的死亡率由r2 增為 r2+e，因此模型一修改為：設(shè)戰(zhàn)前捕獲能力系數(shù)e=0.3, 戰(zhàn)爭中降為e=0.1, 其它參數(shù)與模型(一)的參數(shù)相同。觀察結(jié)果會如何變化？1）當e=0.3時：2）當e=0.1時：分別求出兩種情況下鯊魚在魚類中所占的比例。即計算畫曲線：plot(t,p1(t),t,p2(t),'*')MATLAB編程

13、實現(xiàn)建立兩個M文件function dx=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); function dy=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('shier2',0 15,25 2); x1=x(:,1);x2=x(:,2)

14、; p1=x2./(x1+x2); y1=y(:,1);y2=y(:,2); p2=y2./(y1+y2); plot(t1,p1,'-',t2,p2,'*') 圖中*曲線為戰(zhàn)爭中鯊魚所占比例。結(jié)論：戰(zhàn)爭中鯊魚的比例比戰(zhàn)前高。  3、 Logistic映射logistic映射-通向混沌的道路 混沌系統(tǒng)，由于其行為的復(fù)雜性，往往認為其動態(tài)特性（運動方程）也一定非常復(fù)雜，事實并非如此，一個參量很少、動態(tài)特性非常簡單的系統(tǒng)有時也能夠產(chǎn)生混沌現(xiàn)象，以一維蟲口模型為例，假設(shè)某一區(qū)域內(nèi)的現(xiàn)有蟲口數(shù)為yn，昆蟲的繁殖率為r，且第n代昆蟲不能存活于第n+1代，既無世代

15、交疊，則第n+1代蟲口數(shù)為，r1時，蟲口會無限制地增長；r1時，蟲口最終會趨于消亡，因此需要對模型進行修正。由于環(huán)境的制約和食物有限，因爭奪生存空間發(fā)生相互咬斗事件的最大次數(shù)為，即制約蟲口數(shù)的因素與 成正比，設(shè)咬斗事件的戰(zhàn)死率為則對蟲口的修正項為 ，則有： .令，則 （1）取最大蟲口數(shù)為1，且蟲口數(shù)不能為負，則 ；當 =0.5時，方程有極大值，而 又必須小于1，因而r4，則參量r的取值范圍為1到4，這就得到一個抽象的標準蟲口方程（1）。記映射為 （2）方程（1）可寫為 （3） 這一迭代關(guān)系通常稱為logistic映射。從0,1內(nèi)點x0出發(fā)，由Logistic映射的迭代形成xn= f n(x0)

16、, n = 0,1,2,序列xn稱為x0的軌道。一個看似簡單的系統(tǒng)，隨著參量的不同會表現(xiàn)出截然不同的行為，當r的取值范圍在13時，方程（1）有定態(tài)解 即方程通過多次迭代后趨于一個穩(wěn)定的不動點，此時系統(tǒng)是穩(wěn)定的。為方程的解，稱為周期2點。當r在33.448范圍內(nèi)取值時，經(jīng)過多次迭代，方程（1）出現(xiàn)周期2點和，即和是方程的解，滿足是使解有意義的r最小值。隨著r的增大，r=3.449；3.544；3.564依次出現(xiàn)周期4、周期8、周期16的振蕩解，r的極限值約為3.569。這種行為稱為倍周期分岔，直到r3.5699時，系統(tǒng)進入了混沌狀態(tài)，如下圖所示，此時系統(tǒng)的狀態(tài)不再具有規(guī)律性，而是發(fā)生隨機的波動，

17、使圖d的右側(cè)的大部分區(qū)域被涂黑了，仔細觀察發(fā)現(xiàn)，混沌區(qū)域并非一片涂斑，而是有粗粗細細的白色“豎線”，稱為周期窗口，隨著參量r的增大（如 ）時，混沌突然消失，系統(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)。當然通向混沌的道路不只于此，第二種通向的道路是：從平衡態(tài)到周期運動，再到擬周期運動，直到進入混沌狀態(tài)。第三種通向混沌的方式是陣發(fā)（或

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

19、r i=1:100 x1=r*x*(1-x); x=x1; plot(i,x1,'.'); hold on;end下面程序繪制logistic映射r=4的軌道，觀察其混沌clear all;clf;x=0.007;for i=1:100 x1=4*x*(1-x); x=x1; plot(i,x1,'.'); hold on;end下面程序繪制在混沌狀態(tài)下不同初值x0=0.100001和x0=0.1的軌道的差（對初值的敏感性）clear all;clf;x1=0.100001;x11=0.1;for i=1:1000 x2=4*x1*(1-x1); x1=x2; x22=4*x11*(1-x11); x11=x22; plot(i,x11-x1); hold on;