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

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

2、ns =tan(t+C1)(2)求特解輸入：dsolve(Dy=1+yA2,y(0)=1,x)指定初值為1,自變量為x輸出：ans =tan(x+1/4*pi)2, 21、x y xy (x )y 0 4例2求解二階微分方程y( /2) 2y( /2)2/原方程兩邊都除以x2,得y y (1y)y 0 x 4x輸入：dsolve(D2y+(1/x)*Dy+(1-1/4/xA2)*y=0,y(pi/2)=2,Dy(pi/2)=-2/pi,x)ans =-(exp(x*i)*(pi/2)A(1/2)*i)/xA(1 +(exp(x*i)*exp(-x*2*i)*(pi/2)A(3/2)*2*i)/

3、(pi*xA(1/2)試試能不用用simplify 函數(shù)化簡輸入:simplify(ans)ans =2A(1/2)*piA(1/2)/xA(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*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、(4*t)g =exp(3*t)*cos(4*t)數(shù)值解在微分方程(組)難以獲得解析解的情況下,可以用Matlab方便地求出數(shù)值 解.格式為：t,y = ode23(F,ts,y0,options)注意：? 微分方程的形式：y= 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處給出,(

5、3)ts=t0:k:tf, 那么輸出在區(qū)間t0,tf的等分點(diǎn)給出；? y0為初值條件；? options用于設(shè)定誤差限(缺省就是設(shè)定相對誤差就是10A(-3),絕對誤 差就是10A(-6);ode23就是微分方程組數(shù)值解的低階方法,ode45為中階方法,與ode23類似.例4求解一個經(jīng)典的范得波(Van Der pol)微分方程：2u (u 1)u u0, u(0) 1, u(0) 0解 形式轉(zhuǎn)化：令y u(t)； y u(t).那么以上方程轉(zhuǎn)化一階微分方程組：yy2；y2(1v；)yy1.編寫M文件如下,必須就是M文件表示微分方程組,并保存,一般地,M文件的 名字與函數(shù)名相同,保存位置可以為

6、默認(rèn)的work子目錄,也可以保存在自定義文件夾,這時注意要增加搜索路徑(PathAdd Folder)function dot1=vdpol(t,y)；dot1=y(2)； (1-y(1)A2)*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 Solution3 1ri9ncr:2 .,./ f i f / /

7、,.1 1T./,/ I 1 JBya 0 - 1i :!iy- 1 - iy、i /i J /- 2 -. / KZJI- 302468101214161820Time,Second注:Van der Pol方程描述具有一個非線性振動項(xiàng)的振動子的運(yùn)動過程.最 初,由于它在非線性電路上的應(yīng)用而引起廣泛興趣.一般形式為u (u2 1)u u 0.圖形解無論就是解析解還就是數(shù)值解,都不如圖形解直觀明了.即使就是在得到了 解析解或數(shù)值解的情況下,作出解的圖形,仍然就是一件深受歡送的事.這些都可 以用Matlab方便地進(jìn)行.(1)圖示解析解如果微分方程(組)的解析解為：y=f (x),那么可以用Mat

8、lab函數(shù)fplot作出其 圖形：fplot(fun,lims)其中:fun給出函數(shù)表達(dá)式;lims=xmin xmax ymin ymax限定坐標(biāo)軸的大小. 例如fplot(sin(1/x), 001 0 1 -11)(2)圖示數(shù)值解設(shè)想已經(jīng)得到微分方程(組)的數(shù)值解(x,y).可以用Matlab函數(shù)plot(x,y) 直接作出圖形.其中x與y為向量(或矩陣).2、Volterra模型(食餌捕食者模型)意大利生物學(xué)家Ancona曾致力于魚類種群相互制約關(guān)系的研究,她從第一 次世界大戰(zhàn)期間,地中海各港口捕獲的幾種魚類捕獲量百分比的資料中 ,發(fā)現(xiàn)鯊魚 的比例有明顯增加(見下表).年代191419

9、15191619171918百分比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é)模型,定量地答復(fù)這個問題.1、符號說明：x1(t), x2(t)分別就是食餌、捕食者(鯊魚)在t時刻的數(shù)量；r1, r2就是食餌、捕食者的固有增長率；入1就是捕食者掠取食餌的水平,入2就是

10、食餌對捕食者的供養(yǎng)水平；2、根本假設(shè)：捕食者的存在使食餌的增長率降低,假設(shè)降低的程度與捕食者數(shù)量成正比,即dx1dtXi(r11 X2)食餌對捕食者的數(shù)量x2起到增長的作用,其程度與食餌數(shù)量X1成正比,即dx2dtX2( J2X1)綜合以上與,得到如下模型：模型一:不考慮人工捕獲的情況1X2)dXiXi(ri dtdX2dTX222X1)dX1dtdx2dtX1(0) 25,X2(0) 2dx(1)=x(1)*(1-0 dx(2)=x(2)*(-0、1*x(2);、5+0、02*x(1);該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關(guān)系沒有考慮食餌與捕食者自身的阻滯作用,就是V

11、olterra提出的最簡單的模型. 給定一組具體數(shù)據(jù),用matlab軟件求解.食餌：“=1, 入 1= 0、1,X10= 25;捕食者鯊魚：r2=0、5,入 2=0、02, x20= 2;X1(1 0.1X2)X2( 0.5 0.02x1)編制程序如下1、建立m-文件shier、m如下:function dx=shier(t,x)dx=zeros(2,1); % 初始化、1:15,25 2);2、在命令窗口執(zhí)行如下程序t,x=ode45(shier,0:0plot(t,x(:,1),-,t,x(:,2),*),grid10090f 1807060J i50401I13020c*10/ M- -

12、14 | lyjMjjfri.0051015圖中,藍(lán)色曲線與綠色曲線分別就是食餌與鯊魚數(shù)量隨時間的變化情況,從圖中可以瞧出它們的數(shù)量都呈現(xiàn)出周期性,而且鯊魚數(shù)量的頂峰期稍滯后于食餌數(shù)量的頂峰期.畫出相軌跡圖：plot(x(:,1),x(:,2)30模型二考慮人工捕獲的情況假設(shè)人工捕獲水平系數(shù)為e,相當(dāng)于食餌的自然增長率由ri降為ri-e,捕食者的死亡率由2增為r 2+e,因此模型一修改為：設(shè)戰(zhàn)前捕獲水平系數(shù)e=0、3,戰(zhàn)爭中降為e=0、1,其它參數(shù)與模型(一)的dx1dtdx2dtxiie1X2X2 2e2X1參數(shù)相同.觀察結(jié)果會如何變化?dx1dtx(0.7 0.IX2)dx2dtx2( 0

13、.8 0.02x1)Xi(0)25,X2(0) 2i)當(dāng) e=0、3 時:2)當(dāng) e=0、i 時:dyidt dy2yi(0.9 0.1y2)dtyi(0)y2( 0.6 0.02yi)25*(0)分別求出兩種情況下鯊魚在魚類中所占的比例. 即計(jì)算Pi (t)X2 (t)xi(t)x2(t)P2(t)y2(t)yi(t) y2 (t)畫曲線:plot(t,pi(t),t,p2(t),*)MATLAB編程實(shí)現(xiàn)建立兩個M文件function dx=shieri(t,x)dx=zeros(2,1); dx(i)=x(i)*(0 dx(2)=x(2)*(-0、7-0、8+0i*x(2);、02*x(i

14、);function dy=shier2(t,y)dy=zeros(2,i);dy(i)=y(i)*(0 dy(2)=y(2)*(-0 運(yùn)行以下程序:、9-0、6+0i*y(2);、02*y(i);ti,x=ode45(shieri,0 i5,25 2);t2,y=ode45(shier2,0 i5,25 2);xi=x(:,i);x2=x(:,2);pi=x2、/(xi+x2);yi=y(:,i);y2=y(:,2);p2=y2、/(yi+y2);plot(ti,pi,-,t2,p2,*)圖中*曲線為戰(zhàn)爭中鯊魚所占比例 結(jié)論：戰(zhàn)爭中鯊魚的比例比戰(zhàn)前高3、 Logistic 映射logisti

15、c映射-通向混沌的道路混沌系統(tǒng),由于其行為的復(fù)雜性,往往認(rèn)為其動態(tài)特性(運(yùn)動方程)也一定非常復(fù)雜,事實(shí) 并非如此,一個參量很少、動態(tài)特性非常簡單的系統(tǒng)有時也能夠產(chǎn)生混沌現(xiàn)象,以一維蟲口模型為例,假設(shè)某一區(qū)域內(nèi)的現(xiàn)有蟲口數(shù)為yn,昆蟲的繁殖率為r,且第n代昆蟲不能存活于第n+1代,既無世代交疊,那么第n+1代蟲口數(shù)為yn 1 ryn ,r 1時,蟲口會無限制地增長;r 1 時,蟲口最終會趨于消亡,因此需要對模型進(jìn)行修正.由于環(huán)境的制約與食物有限,因爭奪生1 2存空間發(fā)生相互咬斗事件的最大次數(shù)為一yn(yn1),即制約蟲口數(shù)的因素與 yn成正比,設(shè)2咬斗事件的戰(zhàn)死率為3那么對蟲口的修正項(xiàng)為y2 ,

16、那么有：2yn 1 ryn yn、令Xn - Yn,那么rxn 1rx n(1 xn )(1)取最大蟲口數(shù)為1,且蟲口數(shù)不能為負(fù),那么/巨10口 ;當(dāng)/ =0、5時,方程有極大值,xn 11而 /十1又必須小于1,因而r3、5699時,系統(tǒng)進(jìn)入了混 沌狀態(tài),如下列圖所示,此時系統(tǒng)的狀態(tài)不再具有規(guī)律性,而就是發(fā)生隨機(jī)的波動,使圖d的右側(cè)的大局部區(qū)域被涂黑了,仔細(xì)觀察發(fā)現(xiàn),混沌區(qū)域并非一片涂斑,而就是有粗粗細(xì)細(xì)的白色“豎線,稱為周期窗口 ,隨著參量r的增大(如r 1 J8)時,混沌忽然消失,系統(tǒng)出現(xiàn)周 期三的穩(wěn)定狀態(tài),Logistic 映射分岔圖接著倍周期分岔以更快的速度進(jìn)行,再次進(jìn)入混沌狀態(tài).如

17、果將周期窗口放大,發(fā)現(xiàn)其結(jié) 構(gòu)與分岔圖的整體結(jié)構(gòu)具有相似性,而且就是一種無限嵌套的自相似結(jié)構(gòu).可以瞧出,通過改變系統(tǒng)參量,使系統(tǒng)進(jìn)入混沌的第一種模式就是倍周期分岔,即由不動點(diǎn)一周期二一周期四一無限倍周期一進(jìn)入混沌狀態(tài).當(dāng)然通向混沌的道路不只于此,第二種通向的道路就是：從平衡態(tài)到周期運(yùn)動,再到擬周期運(yùn)動,直到進(jìn)入混沌狀態(tài).第三種通向 混沌的方式就是陣發(fā)或間歇道路,即系統(tǒng)在近似周期運(yùn)動的過程中,通過改變參量,系統(tǒng)會出現(xiàn)陣發(fā)性混沌過程,隨著參量的調(diào)整,陣發(fā)性混沌越來越頻繁,近似的周期運(yùn)動越來越少, 最后進(jìn)入混沌.Matlab程序下面程序繪制r=2,X0=0、3的軌道clear all ;clf;x=

18、0、3;r=2;n=input( Please input a number:);for i=1:nx1=r*x*(1-x);x=x1;plot(i,x1,、)；hold on;end下面程序繪制logistic映射r=3,5的軌道,觀察就是否有周期點(diǎn)clear all ;clf;x=0、3;r=3 、5;for i=1:100x1=r*x*(1-x);x=x1;plot(i,x1,、)；hold on;end下面程序繪制logistic映射r=4的軌道,觀察其 混沌clear all ;clf;x=0 、 007;for i=1:100x1=4*x*(1-x);x=x1;plot(i,x1,

19、、)；hold on;end下面程序繪制在混沌狀態(tài)下不同初值刈=0、100001與x0=0、1的軌道的差(對初值的敏感0.510.50.490.480.470.460.450.440.431000.90.80.70.60.50.4800.90.80.70.60.50.40.41040-0.2-0.4-0.6-0.80.42 160 7 -70, -80- 90 , , 100 .0.2 (010020030040050060070080090010000.2性)clear all ;clf;x1=0、100001;x11=0、1;for i=1:1000x2=4*x1*(1-x1);x1=x2;x22=4*x11*(1-x11);x11=x22;plot(i,x11-x1);hold on;end下面程序繪制logistic映射的分岔圖clear all ;clf;x=0、2;for r=1:0、001:4for i=1:18x1=r*x*(1-x);x=x1;if i9plot(r,x);hold on;endendend10.90.80.70.60.50.40.30.20.101.52.53.5蛛網(wǎng)迭代clear all ;clf;x=0、