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

.10/10實(shí)驗(yàn)二：微分方程與差分方程模型Matlab求解一、實(shí)驗(yàn)?zāi)康腫1]掌握解析、數(shù)值解法,并學(xué)會(huì)用圖形觀察解的形態(tài)和進(jìn)行解的定性分析；[2]熟悉MATLAB軟件關(guān)于微分方程求解的各種命令；[3]通過范例學(xué)習(xí)建立微分方程方面的數(shù)學(xué)模型以及求解全過程；[4]熟悉離散Logistic模型的求解與混沌的產(chǎn)生過程。二、實(shí)驗(yàn)原理1.微分方程模型與MATLAB求解解析解用MATLAB命令dsolve<‘eqn1’,’eqn2’,...>求常微分方程〔組的解析解。其中〔1微分方程例1求解一階微分方程<1>求通解輸入：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>例2求解二階微分方程原方程兩邊都除以,得輸入: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ù)化簡(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〔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*t>g=exp<3*t>*cos<4*t>數(shù)值解在微分方程<組>難以獲得解析解的情況下,可以用Matlab方便地求出數(shù)值解。格式為:[t,y]=ode23<'F',ts,y0,options>注意：微分方程的形式：y'=F<t,y>,t為自變量,y為因變量〔可以是多個(gè),如微分方程組；[t,y]為輸出矩陣,分別表示自變量和因變量的取值；F代表一階微分方程組的函數(shù)名〔m文件,必須返回一個(gè)列向量,每個(gè)元素對(duì)應(yīng)每個(gè)方程的右端；ts的取法有幾種,〔1ts=[t0,tf]表示自變量的取值范圍,〔2ts=[t0,t1,t2,…,tf],則輸出在指定時(shí)刻t0,t1,t2,…,tf處給出,〔3ts=t0:k:tf,則輸出在區(qū)間[t0,tf]的等分點(diǎn)給出；y0為初值條件；options用于設(shè)定誤差限〔缺省是設(shè)定相對(duì)誤差是10^<-3>,絕對(duì)誤差是10^<-6>；ode23是微分方程組數(shù)值解的低階方法,ode45為中階方法,與ode23類似。例4求解一個(gè)經(jīng)典的范得波〔VanDerpol微分方程：解形式轉(zhuǎn)化：令。則以上方程轉(zhuǎn)化一階微分方程組：。編寫M文件如下,必須是M文件表示微分方程組,并保存,一般地,M文件的名字與函數(shù)名相同,保存位置可以為默認(rèn)的work子目錄,也可以保存在自定義文件夾,這時(shí)注意要增加搜索路徑〔File\SetPath\AddFolderfunctiondot1=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í)行：注：VanderPol方程描述具有一個(gè)非線性振動(dòng)項(xiàng)的振動(dòng)子的運(yùn)動(dòng)過程。最初,由于它在非線性電路上的應(yīng)用而引起廣泛興趣。一般形式為。圖形解無論是解析解還是數(shù)值解,都不如圖形解直觀明了。即使是在得到了解析解或數(shù)值解的情況下,作出解的圖形,仍然是一件深受歡迎的事。這些都可以用Matlab方便地進(jìn)行。<1>圖示解析解如果微分方程〔組的解析解為：y=f<x>,則可以用Matlab函數(shù)fplot作出其圖形：fplot<'fun',lims>其中：fun給出函數(shù)表達(dá)式；lims=[xminxmaxyminymax]限定坐標(biāo)軸的大小。例如fplot<'sin<1/x>',[0.010.1-11]><2>圖示數(shù)值解設(shè)想已經(jīng)得到微分方程〔組的數(shù)值解<x,y>?？梢杂肕atlab函數(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)爭(zhēng)為什么使鯊魚數(shù)量增加？是什么原因？因?yàn)閼?zhàn)爭(zhēng)使捕魚量下降,食用魚增加,顯然鯊魚也隨之增加。但為何鯊魚的比例大幅增加呢？生物學(xué)家Ancona無法解釋這個(gè)現(xiàn)象,于是求助于著名的意大利數(shù)學(xué)家V.Volterra,希望建立一個(gè)食餌—捕食者系統(tǒng)的數(shù)學(xué)模型,定量地回答這個(gè)問題。1、符號(hào)說明：①x1<t>,x2<t>分別是食餌、捕食者〔鯊魚在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成正比,即綜合以上①和②,得到如下模型：模型一：不考慮人工捕獲的情況該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關(guān)系,沒有考慮食餌和捕食者自身的阻滯作用,是Volterra提出的最簡(jiǎn)單的模型。給定一組具體數(shù)據(jù),用matlab軟件求解。食餌：r1=1,λ1=0.1,x10=25;捕食者<鯊魚>：r2=0.5,λ2=0.02,x20=2;編制程序如下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圖中,藍(lán)色曲線和綠色曲線分別是食餌和鯊魚數(shù)量隨時(shí)間的變化情況,從圖中可以看出它們的數(shù)量都呈現(xiàn)出周期性,而且鯊魚數(shù)量的高峰期稍滯后于食餌數(shù)量的高峰期。畫出相軌跡圖：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ù)相同。觀察結(jié)果會(huì)如何變化？1當(dāng)e=0.3時(shí)：2當(dāng)e=0.1時(shí)：分別求出兩種情況下鯊魚在魚類中所占的比例。即計(jì)算畫曲線：plot<t,p1<t>,t,p2<t>,'*'>MATLAB編程實(shí)現(xiàn)建立兩個(gè)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>>;運(yùn)行以下程序：[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,'*'>圖中‘*’曲線為戰(zhàn)爭(zhēng)中鯊魚所占比例。結(jié)論：戰(zhàn)爭(zhēng)中鯊魚的比例比戰(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)象,以一維蟲口模型為例,假設(shè)某一區(qū)域內(nèi)的現(xiàn)有蟲口數(shù)為yn,昆蟲的繁殖率為r,且第n代昆蟲不能存活于第n+1代,既無世代交疊,則第n+1代蟲口數(shù)為,r＞1時(shí),蟲口會(huì)無限制地增長(zhǎng)；r＜1時(shí),蟲口最終會(huì)趨于消亡,因此需要對(duì)模型進(jìn)行修正。由于環(huán)境的制約和食物有限,因爭(zhēng)奪生存空間發(fā)生相互咬斗事件的最大次數(shù)為,即制約蟲口數(shù)的因素與成正比,設(shè)咬斗事件的戰(zhàn)死率為β則對(duì)蟲口的修正項(xiàng)為,則有：.令,則〔1取最大蟲口數(shù)為1,且蟲口數(shù)不能為負(fù),則；當(dāng)=0.5時(shí),方程有極大值,而又必須小于1,因而r＜4,則參量r的取值范圍為1到4,這就得到一個(gè)抽象的標(biāo)準(zhǔn)蟲口方程〔1。記映射為〔2方程〔1可寫為〔3這一迭代關(guān)系通常稱為logistic映射。從[0,1]內(nèi)點(diǎn)x0出發(fā),由Logistic映射的迭代形成xn=fn<x0>,n=0,1,2,…序列{xn}稱為x0的軌道。一個(gè)看似簡(jiǎn)單的系統(tǒng),隨著參量的不同會(huì)表現(xiàn)出截然不同的行為,當(dāng)r的取值范圍在1~3時(shí),方程〔1有定態(tài)解即方程通過多次迭代后趨于一個(gè)穩(wěn)定的不動(dòng)點(diǎn),此時(shí)系統(tǒng)是穩(wěn)定的。為方程的解,稱為周期2點(diǎn)。當(dāng)r在3~3.448范圍內(nèi)取值時(shí),經(jīng)過多次迭代,方程〔1出現(xiàn)周期2點(diǎn)和,即和是方程的解,滿足是使解有意義的r最小值。隨著r的增大,r=3.449；3.544；3.564…依次出現(xiàn)周期4、周期8、周期16…的振蕩解,r的極限值約為3.569。這種行為稱為倍周期分岔,直到r＞3.5699時(shí),系統(tǒng)進(jìn)入了混沌狀態(tài),如下圖所示,此時(shí)系統(tǒng)的狀態(tài)不再具有規(guī)律性,而是發(fā)生隨機(jī)的波動(dòng),使圖d的右側(cè)的大部分區(qū)域被涂黑了,仔細(xì)觀察發(fā)現(xiàn),混沌區(qū)域并非一片涂斑,而是有粗粗細(xì)細(xì)的白色"豎線",稱為周期窗口,隨著參量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)具有相似性,而且是一種無限嵌套的自相似結(jié)構(gòu)?？梢钥闯?通過改變系統(tǒng)參量,使系統(tǒng)進(jìn)入混沌的第一種模式是倍周期分岔,即由不動(dòng)點(diǎn)→周期二→周期四→…無限倍周期→進(jìn)入混沌狀態(tài)。當(dāng)然通向混沌的道路不只于此,第二種通向的道路是：從平衡態(tài)到周期運(yùn)動(dòng),再到擬周期運(yùn)動(dòng),直到進(jìn)入混沌狀態(tài)。第三種通向混沌的方式是陣發(fā)〔或間歇道路,即系統(tǒng)在近似周期運(yùn)動(dòng)的過程中,通過改變參量,系統(tǒng)會(huì)出現(xiàn)陣發(fā)性混沌過程,隨著參量的調(diào)整,陣發(fā)性混沌越來越頻繁,近似的周期運(yùn)動(dòng)越來越少,最后進(jìn)入混沌。Matlab程序下面程序繪制r=2,x0=0.3的軌道clearall;clf;x=0.3;r=2;n=input<'Pleaseinputanumber:'>;fori=1:nx1=r*x*<1-x>;x=x1;plot<i,x1,'.'>;holdon;end下面程序繪制logistic映射r=3,5的軌道,觀察是否有周期點(diǎn)clearall;clf;x=0.3;r=3.5;fori=1:100x1=r*x*<1-x>;x=x1;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下面程序繪制在混沌狀態(tài)下不同初值x0=0.100001和x0=0.1的軌道的差〔對(duì)初值的敏感性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>;