''實驗二：微分方程與差分方程模型Matlab求解一、實驗?zāi)康恼莆战馕觥?shù)值解法，并學(xué)會用圖形觀察解的形態(tài)和進行解的定性分析；精品文檔放心下載熟悉MATLAB軟件關(guān)于微分方程求解的各種命令；感謝閱讀通過范例學(xué)習(xí)建立微分方程方面的數(shù)學(xué)模型以及求解全過程；感謝閱讀熟悉離散Logistic模型的求解與混沌的產(chǎn)生過程。謝謝閱讀二、實驗原理1.微分方程模型與MATLAB求解解析解MATLAB命令dsolve(‘eqn1’,’eqn2’,...)求常微分方程（組）的解析解。其中‘eqni'表示第i個微分方程，Dny表示y的n階導(dǎo)數(shù),默認的自變量為t。謝謝閱讀（1）微分方程例1 求解一階微分方程dydx1y2謝謝閱讀求通解輸入：dsolve('Dy=1+y^2')輸出：ans=tan(t+C1)（2）求特解輸入：dsolve('Dy=1+y^2','y(0)=1','x')謝謝閱讀指定初值為1，自變量為x輸出：ans=tan(x+1/4*pi)''xyxy(x1422例2求解二階微分方程y(/2)2y(/2)2/原方程兩邊都除以2，得1(11)y0xyxy4x2輸入:dsolve('D2y+(1/x)*Dy+(1-1/4/x^2)*y=0','y(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')精品文檔放心下載=exp(3*t)*(C1*sin(4*t)+C2*cos(4*t))感謝閱讀=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')精品文檔放心下載=exp(3*t)*sin(4*t)=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處給出，（3）ts=t0:k:tf,則輸出在區(qū)間[t0,tf]的等分點給出；精品文檔放心下載y0為初值條件；options用于設(shè)定誤差限（缺省是設(shè)定相對誤差是10^(-3)，絕對誤差是10^(-6)）；感謝閱讀ode23是微分方程組數(shù)值解的低階方法，ode45為中階方法，與ode23類似。謝謝閱讀例4 求解一個經(jīng)典的范得波（VanDerpol）微分方程：謝謝閱讀u''(u21)u'u0，u(0)1,u'(0)0謝謝閱讀解形式轉(zhuǎn)化：令yu(t);yu'(t)。則以上方程轉(zhuǎn)化一階微分方程組：12yy;y(1y2)yy。122121編寫M文件如下，必須是M文件表示微分方程組，并保存，一般地，M文件的名字與函數(shù)名相同，保存位置可以為默認的work子目錄，也可以保存在自定義文件夾，這時注意要增加搜索路徑（File\SetPath\AddFolder）謝謝閱讀function dot1=vdpol(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('VanDerPolSolution');精品文檔放心下載xlabel('Time,Second');ylabel('y(1)andy(2)')謝謝閱讀執(zhí)行：VanDerPolSolution321)2d 0a1y-1-2-302468101214161820Time,Second注：VanderPol方程描述具有一個非線性振動項的振動子的運動過程。最初,由于它在非線性電路上的應(yīng)用而引起廣泛興趣。一般形式為精品文檔放心下載u''(u21)u'u0。圖形解無論是解析解還是數(shù)值解，都不如圖形解直觀明了。即使是在得到了解析解或數(shù)值解的情況下，作出解的圖形，仍然是一件深受歡迎的事。這些都可以用Matlab方便地進行。謝謝閱讀(1)圖示解析解如果微分方程（組）的解析解為：y=f(x)，則可以用Matlab函數(shù)fplot作出其圖形：精品文檔放心下載fplot('fun',lims)其中：fun給出函數(shù)表達式；lims=[xminxmaxyminymax]限定坐標(biāo)軸的大小。例如謝謝閱讀fplot('sin(1/x)',[0.010.1-11])謝謝閱讀''10.80.60.40.20-0.2-0.4-0.6-0.8-10.010.020.030.040.050.060.070.080.090.1(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)鯊魚的比例有明顯增加（見下表）。謝謝閱讀年代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、符號說明：① x(t),x(t)分別是食餌、捕食者（鯊魚）在t時刻的數(shù)量；謝謝閱讀1 2② r1,r2是食餌、捕食者的固有增長率；③λ1是捕食者掠取食餌的能力，λ2是食餌對捕食者的供養(yǎng)能力；''2、基本假設(shè)：① 捕食者的存在使食餌的增長率降低，假設(shè)降低的程度與捕食者數(shù)量成正比，謝謝閱讀dxx(rx)1即dt1112② 食餌對捕食者的數(shù)量x2起到增長的作用，其程度與食餌數(shù)量x1成正比，感謝閱讀dx即dt2x2(r22x1)感謝閱讀綜合以上①和②，得到如下模型：模型一：不考慮人工捕獲的情況dxx(rx)1dt1112dx2x(rx)2221該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關(guān)感謝閱讀系，沒有考慮食餌和捕食者自身的阻滯作用，是Volterra提出的最簡單的模型。感謝閱讀給定一組具體數(shù)據(jù)，用matlab軟件求解。食餌：r1=1, λ1=0.1, x10=25;感謝閱讀捕食者(鯊魚)：r2=0.5,λ2=0.02, x20=2;感謝閱讀dxdt1dx2

x(10.1x)12x(0)25,x(0)2x(0.50.02x)1221編制程序如下1、建立m-文件shier.m如下：functiondx=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,[252]);精品文檔放心下載plot(t,x(:,1),'-',t,x(:,2),'*'),grid謝謝閱讀''10090807060504030201000 5 10 15圖中，藍色曲線和綠色曲線分別是食餌和鯊魚數(shù)量隨時間的變化情況，從圖中可以看出它們的數(shù)量都呈現(xiàn)出周期性，而且鯊魚數(shù)量的高峰期稍滯后于食餌數(shù)量的高峰期。謝謝閱讀畫出相軌跡圖：plot(x(:,1),x(:,2))3025201510500 10 20 30 40 50 60 70 80 90 100謝謝閱讀模型二 考慮人工捕獲的情況假設(shè)人工捕獲能力系數(shù)為e，相當(dāng)于食餌的自然增長率由r1降為r1-e，捕食者的死亡率由r2增為r2+e，因此模型一修改為：精品文檔放心下載設(shè)戰(zhàn)前捕獲能力系數(shù)e=0.3,戰(zhàn)爭中降為e=0.1,其它參數(shù)與模型(一)的參感謝閱讀dxx[(re)x]11112dtdxx[(re)x]2dt2221''數(shù)相同。觀察結(jié)果會如何變化？dxx(0.70.1x)1dt12dx2x(0.80.02x)dt21x(0)25,x(0)2121）當(dāng)e=0.3時：2）當(dāng)e=0.1時：dyy(0.90.1y)112dtdyy(0.60.02y)2dt21y1(0)25,y2(0)2分別求出兩種情況下鯊魚在魚類中所占的比例。即計算p(t)x(t);p(t)y(t)221x(t)x(t)2y(t)y(t)1212畫曲線：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',[015],[252]);感謝閱讀[t2,y]=ode45('shier2',[015],[252]);感謝閱讀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.100 5 10 15圖中‘*’曲線為戰(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代，既無世代交疊，則第n+1代蟲口數(shù)為y ry，r＞1時，蟲口會無限制精品文檔放心下載n1

n地增長；r＜1時，蟲口最終會趨于消亡，因此需要對模型進行修正。由于環(huán)境的制約和食精品文檔放心下載1物有限，因爭奪生存空間發(fā)生相互咬斗事件的最大次數(shù)為2yn(yn1)，即制約蟲口數(shù)的因謝謝閱讀素與y2n

成正比，設(shè)咬斗事件的戰(zhàn)死率為β則對蟲口的修正項為謝謝閱讀

y2n

，則有：y ryn1

n

y2.n令xn

r

y

n

，則x rx(1x)

（1）n1

n

n取最大蟲口數(shù)為1，且蟲口數(shù)不能為負，則 ；當(dāng) =0.5時，方程有極大感謝閱讀1n1 4一個抽象的標(biāo)準(zhǔn)蟲口方程（1）。''記映射f為f(x)rx(1x)（2）方程（1）可寫為xf(x)（3）n1n這一迭代關(guān)系通常稱為logistic映射。從[0,1]內(nèi)點x0出發(fā)，由Logistic映射的迭代形成精品文檔放心下載x=fn(x),n=0,1,2,…n0序列{xn}稱為x0的軌道。一個看似簡單的系統(tǒng)，隨著參量的不同會表現(xiàn)出截然不同的行為，當(dāng)r的取值范圍在1~3時，方程（1）有定態(tài)解謝謝閱讀x11r即方程通過多次迭代后趨于一個穩(wěn)定的不動點，此時系統(tǒng)是穩(wěn)定的。謝謝閱讀1x1r為方程f(x)x的解，稱為周期2點。謝謝閱讀當(dāng)r在3~3.448范圍內(nèi)取值時，經(jīng)過多次迭代，方程（1）出現(xiàn)周期2點x和x，即感謝閱讀1 2x和x是方程f2(x)x的解，滿足謝謝閱讀1 2xrx(1x)1 2 2xrx(1x)2 1 13是使解有意義的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ū)域被涂黑了，仔細觀察發(fā)現(xiàn)，混沌區(qū)域并非一片涂斑，而是有粗粗細細的謝謝閱讀白色“豎線”，稱為周期窗口，隨著參量r的增大（如r18）時，混沌突然消失，系統(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ā)（或間歇）道路，即系統(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('Pleaseinputanumber:');fori=1:n精品文檔放心下載x1=r*x*(1-x);x=x1;plot(i,x1,'.');holdon;end下面程序繪制logistic映射r=3,5的軌道，觀察是否有周期點謝謝閱讀clearall;clf;x=0.3;r=3.5;fori=1:100x1=r*x*(1-x);x=x1;

0.510.50.490.480.470.460.450.440.430.420 10 20 30 40 50 60 70 80 90 100感謝閱讀10.90.80.70.60.50.40 10 20 30 40 50 60 70 80 90 100感謝閱讀''plot(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;end

10.90.80.70.60.50.40.30.20.100 10 20 30 40 50 60 70 80 90 100感謝閱讀下面程序繪制在混沌狀態(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;forr=1:0.001:4fori=1:18x1=r*x*(1-x);x=x1;ifi>9plot(r,x);holdon;endendend

10.80.60.40.20-0.2-0.4-0.6-0.8-10 100 200 300 400 500 600 700 800 900 1