張佩雪,牛利娟,龐茹一
(西安工程大學(xué)理學(xué)院,陜西 西安 710048)
捕食者和食餌之間動態(tài)關(guān)系的普遍性和重要性,使得捕食者-食餌模型長期以來一直是生態(tài)學(xué)和數(shù)學(xué)生態(tài)學(xué)的主要課題之一[1].如果捕食者或者天敵的數(shù)量過多,食餌有可能滅絕.避難所對捕食者-食餌的相互作用具有穩(wěn)定作用.Feller[2]建立了具有獵物避難所的捕食者-食餌模型.Bick 等[3]研究發(fā)現(xiàn)避難所的存在可以保護(hù)食餌不被滅絕.近年來,離散捕食者-食餌模型受到許多學(xué)者的關(guān)注[4].物種間相互作用的微分方程模型是數(shù)學(xué)在生物學(xué)中的經(jīng)典應(yīng)用之一.馬兆芝等[5]研究發(fā)現(xiàn),避難所大小是影響動力學(xué)行為的一個(gè)關(guān)鍵性因素.Rodriguez等[6]研究了在避難區(qū),捕食者無法捕食獵物.Teng[7]研究了捕食者-食餌模型的復(fù)雜動力學(xué)行為.Fu 等[8]討論了Holling II 型捕食者-食餌模型的動力學(xué)性質(zhì).Hu 等[9]研究了捕食者-食餌模型的Flip 分岔和Neimark-Sacker 分岔.程利芳[10]討論了分岔理論在生態(tài)模型中的應(yīng)用.Khan 等[11]借助中心流形定理和分岔理論研究了離散捕食者-食餌模型的分岔.
介紹離散模型之前,我們引入具有獵物避難所和Holling II 型功能反應(yīng)函數(shù)的連續(xù)捕食者-食餌模型[12]:
其中:x,y分別表示在t時(shí)刻的食餌和捕食者種群密度,a,k,α,β,γ,c都是正常數(shù).這里α代表食餌的內(nèi)在生長速率;k代表食餌的攜帶能力;γ是捕食者的死亡率;β/α是單位時(shí)間內(nèi)每個(gè)捕食者可吃掉的最大食餌數(shù)量; 1/a是達(dá)到該比率一半所需的食餌密度;c是表示每個(gè)捕獲的食餌的新生捕食者數(shù)量的轉(zhuǎn)換因子;m表示食餌避難率,βx/(1+ax)表示Holling II 型[13]功能反應(yīng)函數(shù).
下面我們用歐拉向前方法離散化模型(1):
本文討論了模型(2)平衡點(diǎn)的存在性,利用中心流形定理和分岔理論給出模型(2)的Flip 分岔和Neimark-Sacker 分岔的存在條件.通過數(shù)值模擬發(fā)現(xiàn)模型(2)出現(xiàn)周期解、極限環(huán)、穩(wěn)定周期解、穩(wěn)定極限環(huán)、不穩(wěn)定極限環(huán).當(dāng)參數(shù)在某個(gè)特定值時(shí),模型(2)的解會出現(xiàn)混沌現(xiàn)象.螨類捕食者-食餌的相互作用通常表現(xiàn)出空間避難所,這為食餌提供了一定程度的保護(hù),使其免受捕食,減少了因捕食而滅絕的機(jī)會,此模型對實(shí)際生活和理論指導(dǎo)有一定的參考價(jià)值.
很容易看出E0(0,0)和E1(k,0)是模型(2)的兩個(gè)平衡點(diǎn).
此外,當(dāng)滿足條件cβ-γa>0 和0≤m≤1-γ/(k(cβ-γa))時(shí),E2(x2,y2)存在正平衡點(diǎn),x2=γ/k(cβ-γa)(1-m),y2=ac(k(cβ-γa)(1-m))/k(cβ-γa)(1-m)2.
下面討論正平衡點(diǎn)E2(x2,y2)是否存在Flip 分岔和Neimark-Sacker 分岔.記在E2(x2,y2)的雅可比矩陣為
J(E2)的相應(yīng)特征方程可以寫成F(λ)=λ2-trJ(E2)λ+detJ(E2).其中,
顯然F(1)>0,當(dāng)α=α1=[4k(1+a(1-m)x2)2+2kaβ(1-m)2x2y2+kγβ(1-m)y2]/2x2(1+a(1-m)x2)2,此時(shí)F(-1)=0,trJ(E2)≠0.當(dāng)trJ(E2)=0 時(shí),有α=α*=[4k(1+a(1-m)x2)2+kaβ(1-m)2x2y2]/2x2(1+a(1-m)x2)2.則在E2可能存在Flip 分岔,定義M={(a,m,k,α,β) :α=α1,α≠α*,a>0,m>0,k>0,γ>0,β>0}.當(dāng)detJ(E2)=1 時(shí),存在共軛復(fù)根,且|λ1|=|λ2|=1.c=c1為(cβ-γa)2/c(k(cβ-γa)(1-m)-γ)-[x2aβ(1-m)+γβ]/x2(1-m)(1+a(1-m)x2)2=0 的解.則在E2可能存在Neimark-Sacker 分岔,定義N={(a,m,k,α,β):c=c1,a>0,m>0,k>0,γ>0,β>0}.
本節(jié)應(yīng)用中心流形定理[14]研究了Flip 分岔的存在性.選擇α作為分岔參數(shù),給參數(shù)α一個(gè)小擾動,并且令un=xn-x2,vn=yn-y2,模型(2)變?yōu)?/p>
將模型(3)在(un,vn,)=(0,0,0)時(shí)進(jìn)行泰勒級數(shù)展開至二階
令
由于|λ1|=-1,|λ2|≠-1,通過計(jì)算,得到矩陣J(E2)的特征值分別為λ1=-1,λ2=3-αx2/k+αβ(1-m)2x2y2/(1+a(1-m)x2)2.
設(shè)矩陣
使用翻轉(zhuǎn)
則模型(3)變成以下形式
對模型(5)應(yīng)用中心流形定理.假設(shè)Wc(0,0,0)表示模型(5)在=0 的小鄰域內(nèi)的中心流形,則Wc(0,0,0)計(jì)算為
此外,對于限制于中心流形Wc(0,0,0)的映射,我們考慮映射G*
在(un,vn,)=(0,0,0)時(shí)定義兩個(gè)判別量l1≠0 和l2≠0.
因此,根據(jù)上述分析和文獻(xiàn)[14]中的定理3.1,我們給出了在正平衡點(diǎn)Flip 分岔存在的條件.
定理1如果l2≠0,模型(2)在E2(x2,y2)處存在Flip 分岔.如果l2>0,則從E2(x2,y2)分岔的周期2 點(diǎn)是穩(wěn)定的; 如果l2<0,則從E2(x2,y2) 分岔的周期2 點(diǎn)是不穩(wěn)定的.
本節(jié)選擇了c作為分岔參數(shù)來研究在E2(x2,y2)的Neimark-Sacker 分岔.令c=c1+,un=xn-x2,vn=yn-y2,其中||?1,模型(2)變?yōu)?/p>
將模型(7)在(un,vn)=(0,0)時(shí)進(jìn)行泰勒級數(shù)展開至三階
其中c11,c12,···,c15,c21,c22,···,c25同模型(4)中的a11,a12,···,a15,a21,a22,···,a25,
模型(8)在原點(diǎn)(0,0)處線性化的特征方程為λ2-P()λ+Q(),其中,
使用翻轉(zhuǎn)
則模型(7)變成以下形式
其中,
定義第一Lyapunov 指數(shù)如下所示[14]
其中,
因此,根據(jù)上述分析和文獻(xiàn)[15]中的Neimark-Sacker 分岔存在理論,我們得到定理2,表明Neimark-Sacker分岔存在和方向的參數(shù)條件.
定理2如果滿足非退化線性條件且L≠0,模型(2)在E2(x2,y2)處存在Neimark-Sacker 分岔.若L<0(或者L>0),則當(dāng)>0(或者<0) 時(shí),吸引(或者排斥) 不變閉合曲線從E2(x2,y2)分岔.
例1選擇參數(shù)a=0.1,c=4,k=0.4,m=0.5,β=0.3,γ=0.2,通過數(shù)值計(jì)算得平衡點(diǎn)E2(0.338 983,1.269 07)和α=2.454 24,λ1=-1,λ2=0.926 375.模型(2)具有平衡點(diǎn)E2(x2,y2)和(a,m,k,α,β)∈M.我們只改變參數(shù)α以查看模型(2)的動力學(xué)行為變化.其中α∈[2.2,3.4].如圖1(a)所示,對于平衡點(diǎn)E2(x2,y2),在α<2.454 24 時(shí)是穩(wěn)定的,當(dāng)α=2.454 24 時(shí)失去了穩(wěn)定性,當(dāng)α>2.454 24 時(shí),存在Flip 分岔.此外,隨著α的增大,出現(xiàn)了一個(gè)混沌集合.圖2(a)是參數(shù)α∈[2.2,3.4] 時(shí)模型(2)的相圖.

圖1 模型(2)的分岔圖

圖2 模型(2)不同α 值、c 值的相圖
例2選擇參數(shù)a=5,k=0.4,α=3,m=0.1,β=0.4,γ=1.5,通過數(shù)值計(jì)算得出平衡點(diǎn)E2(0.222 222,0.740 741)和c=3.75,λ1,2=0.499 999±0.866 026i.模型(2)具有平衡點(diǎn)E2(x2,y2)和(a,m,k,α,β)∈N.我們只改變參數(shù)c以查看模型(2)的動力學(xué)行為變化.其中c∈[3.5,4.2].如圖1(b)所示,對于平衡點(diǎn)E2(x2,y2),在c<3.75 時(shí)是穩(wěn)定的,當(dāng)c=3.75 時(shí)失去了穩(wěn)定性,當(dāng)c>3.75 時(shí),存在Neimark-Sacker 分岔.此外,隨著c的增大,出現(xiàn)了一個(gè)混沌集合.圖2(b)是參數(shù)c∈[3.5,4.2]時(shí)模型(2)的相圖,有周期解和極限環(huán)出現(xiàn).