關(guān)于競(jìng)爭(zhēng)模型matlab求解例題

《關(guān)于競(jìng)爭(zhēng)模型matlab求解例題》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在工程資料-天天文庫(kù)。

1、關(guān)于兩個(gè)生物種群的競(jìng)爭(zhēng)模型求解題：，=0.5，=1.6，=2.5，=1.8，=1.6，=1，初始條件1），初始條件2）問(wèn)：（1）.用數(shù)值方法畫出的曲線圖；（2）.畫出相軌線；（3）.從的圖中觀察時(shí)，的變化趨勢(shì)，試圖用穩(wěn)定性的理論進(jìn)行分析。答：(1)先依據(jù)微分方程建立函數(shù)文件xt.m。functiondy=xt(t,y)p1=0.5;p2=1.6;r1=2.5;r2=1.8;N1=1.6;N2=1;dy=zeros(2,1);dy(1)=r1*y(1).*(1-y(1)./N1-p1*y(2)./N2);dy(2)=r2*y(

2、2).*(1-y(2)./N2-p2*y(1)./N1);利用以下程序畫出的曲線圖。t0=0;tf=10;[t,y]=ode45('xt',[t0tf],[0.10.1]);%取初始條件一時(shí)，求微分方程數(shù)值解subplot(1,2,1);plot(t,y(:,1),t,y(:,2),'r');%畫出x1(t),x2(t)曲線圖xlabel('時(shí)間t');ylabel('種群密度x');gtext('x1(t)');gtext('x2(t)');%作標(biāo)記title('初值分別為0.1,0.1時(shí)兩種群密度與時(shí)間關(guān)系');grid

3、on;[t2,y2]=ode45('xt',[t0tf],[12]);%取初始條件二時(shí)，求微分方程數(shù)值解subplot(1,2,2);plot(t2,y2(:,1),t2,y2(:,2),'r');%畫出x1(t),x2(t)曲線圖xlabel('時(shí)間t');ylabel('種群密度x');gtext('x1(t)');gtext('x2(t)');%作標(biāo)記title('初值分別為1,2時(shí)兩種群密度與時(shí)間關(guān)系');gridon;當(dāng)取初始條件一時(shí),的曲線圖如圖(一)：圖(一)當(dāng)取初始條件二時(shí),的曲線圖如圖(二)：圖(二)(2)

4、畫出相軌線程序如下：t0=0;tf=10;[t,y]=ode45('xt',[t0tf],[0.10.1]);plot(y(:,1),y(:,2),'b');holdon;plot(y(1,1),y(1,2),'r+');xlabel('甲種群密度x1');ylabel('乙種群密度x2');%作標(biāo)記[t2,y2]=ode45('xt',[t0tf],[12]);holdon;plot(y2(:,1),y2(:,2),'g');holdon;plot(y2(1,1),y2(1,2),'r*');gridon;title('甲

5、乙種群相軌線');作出圖像如圖(三),其中在初始條件一、二下的開始位置分別用“+”和“*”標(biāo)記。圖(三)(3)從圖(一)和圖(二)我們可以看出時(shí)，總趨向其最大容量1.6，而總趨向于0(趨向消亡)。我們可從穩(wěn)定性的理論進(jìn)行分析。求平衡點(diǎn)：令f(x1,x2)=g(x1,x2)=0，得到四個(gè)平衡點(diǎn)：P1(1.6,0)，P2(0,1)，P3(0,0)，P4(4,-3)。令穩(wěn)定條件為p>0而且q>0，作表(一)：pq穩(wěn)定條件P1r1-r2(1-s2)-r1r2(1-s2)?s2>1(s1<1)P2-r1(1-s1)+r2-r1r2(1

6、-s1)s1>1(s2<1)P3-(r1+r2)r1r2?不穩(wěn)定P4[r1(1-s1)+r2(1-s2)](1-s1s2)-1r1r2(1-s1)(1-s2)(1-s1s2)-1s1<1s2<1表(一)又因?yàn)?0.5，=1.6，滿足s2>1且s1<1。依表(一)知P1是唯一穩(wěn)定的平衡點(diǎn)。相平面上，即在x1-x2平面上研究方程解沿著t增加所表現(xiàn)出的趨勢(shì)?？芍獂1(t)和x2(t)的增減由j=1-x1/N1-s1x2/N2和y=1-s2x1/N1-x2/N2決定。在s1<1，s2>1情況下，j=0和y=0將相平面分為三個(gè)區(qū)域：S

7、1：x’1>0,x’2>0;S2：x’1>0,x’2<0;S3：x’1<0,x’2<0.用程序作圖：t0=0;tf=10;x1=0:0.1:1.6;x11=0:0.1:1;p1=0.5;p2=1.6;N1=1.6;N2=1;r1=2.5;r2=1.8;x21=N2*(1-x1./N1)/p1;x22=N2*(1-(p2/N1).*x11);plot(x1,x21,'r',x11,x22,'r')%作直線[t,y]=ode45('xt',[t0tf],[0.10.1]);[t1,y1]=ode45('xt',[t0tf],[1

8、2]);[t2,y2]=ode45('xt',[t0tf],[0.41.2]);holdon;plot(y(:,1),y(:,2),y(1,1),y(1,2),'r*');holdon;plot(y1(:,1),y1(:,2),y1(1,1),y1(1,2),'r*');holdon;pl