實(shí)驗(yàn)二微分方程與差分方程模型Matlab求解

1、實(shí)驗(yàn)二： 微分方程與差分方程模型Matlab求解一、實(shí)驗(yàn)?zāi)康? 掌握解析、數(shù)值解法，并學(xué)會(huì)用圖形觀(guān)察解的形態(tài)和進(jìn)行解的定性分析；2 熟悉MATLAB軟件關(guān)于微分方程求解的各種命令；3 通過(guò)范例學(xué)習(xí)建立微分方程方面的數(shù)學(xué)模型以及求解全過(guò)程；4 熟悉離散 Logistic模型的求解與混沌的產(chǎn)生過(guò)程。 二、實(shí)驗(yàn)原理1. 微分方程模型與MATLAB求解解析解用MATLAB命令dsolve(eqn1,eqn2, .) 求常微分方程（組）的解析解。其中eqni'表示第i個(gè)微分方程，Dny表示y的n階導(dǎo)數(shù),默認(rèn)的自變量為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ù)化簡(jiǎn)輸入: 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為因變量（可以是多個(gè)，如微分方程組）；Ø t, y為輸出矩陣，分別表示自變量和因變量的取值；Ø F代表一階微分方程組的函數(shù)名（m文件，必須返回一個(gè)列向量，每個(gè)元素對(duì)應(yīng)每個(gè)方程的右端）；Ø ts的取法有幾種，（1）ts=t0, tf 表示自變量的取值范圍，（2）ts=t0,t1,t2,tf,則輸出在指定時(shí)刻t0,t1,t2,tf處給出，（3）ts=t0:k:tf,則輸出在區(qū)間t0,tf的等分點(diǎn)給出；Ø y0為初值條件；Ø options用于設(shè)定誤差限（缺省是設(shè)定相對(duì)誤差是10(-3)，絕對(duì)誤差是10(-6)）；ode23是微分方程組數(shù)值解的低階

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

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

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

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

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

11、 1, 1= 0.1, x10= 25; 捕食者(鯊魚(yú))：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 圖中，藍(lán)色曲線(xiàn)和綠色曲線(xiàn)分別是食餌和鯊

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

13、實(shí)現(xiàn)建立兩個(gè)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);運(yùn)行以下程序：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,'*') 圖中*曲線(xiàn)為戰(zhàn)爭(zhēng)中鯊魚(yú)所占比例。結(jié)論：戰(zhàn)爭(zhēng)中鯊魚(yú)的比例比戰(zhàn)前高。  3、 Logistic映射logistic映射-通向混沌的道路 混沌系統(tǒng)，由于其行為的復(fù)雜性，往往認(rèn)為其動(dòng)態(tài)特性（運(yùn)動(dòng)方程）也一定非常復(fù)雜，事實(shí)并非如此，一個(gè)參量很少、動(dòng)態(tài)特性非常簡(jiǎn)單的系統(tǒng)有時(shí)也能夠產(chǎn)生混沌現(xiàn)象，以一維蟲(chóng)口模型為例，假設(shè)某一區(qū)域內(nèi)的現(xiàn)有蟲(chóng)口數(shù)為yn，昆蟲(chóng)的繁殖率為r，且第n代昆蟲(chóng)不能存活于第n+1代，既無(wú)世代

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

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

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

18、間歇）道路，即系統(tǒng)在近似周期運(yùn)動(dòng)的過(guò)程中，通過(guò)改變參量，系統(tǒng)會(huì)出現(xiàn)陣發(fā)性混沌過(guò)程，隨著參量的調(diào)整，陣發(fā)性混沌越來(lái)越頻繁，近似的周期運(yùn)動(dòng)越來(lái)越少，最后進(jìn)入混沌。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的軌道，觀(guān)察是否有周期點(diǎn)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的軌道，觀(guān)察其混沌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的軌道的差（對(duì)初值的敏感性）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;