具有獵物避難所和Holling II 型功能反應函數的離散捕食者-食餌模型的分岔分析*

張佩雪，牛利娟，龐茹一

（西安工程大學理學院，陜西 西安 710048）

0 引言

捕食者和食餌之間動態關系的普遍性和重要性,使得捕食者-食餌模型長期以來一直是生態學和數學生態學的主要課題之一[1].如果捕食者或者天敵的數量過多,食餌有可能滅絕.避難所對捕食者-食餌的相互作用具有穩定作用.Feller[2]建立了具有獵物避難所的捕食者-食餌模型.Bick 等[3]研究發現避難所的存在可以保護食餌不被滅絕.近年來,離散捕食者-食餌模型受到許多學者的關注[4].物種間相互作用的微分方程模型是數學在生物學中的經典應用之一.馬兆芝等[5]研究發現,避難所大小是影響動力學行為的一個關鍵性因素.Rodriguez等[6]研究了在避難區,捕食者無法捕食獵物.Teng[7]研究了捕食者-食餌模型的復雜動力學行為.Fu 等[8]討論了Holling II 型捕食者-食餌模型的動力學性質.Hu 等[9]研究了捕食者-食餌模型的Flip 分岔和Neimark-Sacker 分岔.程利芳[10]討論了分岔理論在生態模型中的應用.Khan 等[11]借助中心流形定理和分岔理論研究了離散捕食者-食餌模型的分岔.

介紹離散模型之前,我們引入具有獵物避難所和Holling II 型功能反應函數的連續捕食者-食餌模型[12]:

其中:x,y分別表示在t時刻的食餌和捕食者種群密度,a,k,α,β,γ,c都是正常數.這里α代表食餌的內在生長速率;k代表食餌的攜帶能力;γ是捕食者的死亡率;β/α是單位時間內每個捕食者可吃掉的最大食餌數量; 1/a是達到該比率一半所需的食餌密度;c是表示每個捕獲的食餌的新生捕食者數量的轉換因子;m表示食餌避難率,βx/(1+ax)表示Holling II 型[13]功能反應函數.

下面我們用歐拉向前方法離散化模型(1):

本文討論了模型(2)平衡點的存在性,利用中心流形定理和分岔理論給出模型(2)的Flip 分岔和Neimark-Sacker 分岔的存在條件.通過數值模擬發現模型(2)出現周期解、極限環、穩定周期解、穩定極限環、不穩定極限環.當參數在某個特定值時,模型(2)的解會出現混沌現象.螨類捕食者-食餌的相互作用通常表現出空間避難所,這為食餌提供了一定程度的保護,使其免受捕食,減少了因捕食而滅絕的機會,此模型對實際生活和理論指導有一定的參考價值.

1 平衡點分析

很容易看出E0(0,0)和E1(k,0)是模型(2)的兩個平衡點.

此外,當滿足條件cβ-γa＞0 和0≤m≤1-γ/(k(cβ-γa))時,E2(x2,y2)存在正平衡點,x2=γ/k(cβ-γa)(1-m),y2=ac(k(cβ-γa)(1-m))/k(cβ-γa)(1-m)2.

下面討論正平衡點E2(x2,y2)是否存在Flip 分岔和Neimark-Sacker 分岔.記在E2(x2,y2)的雅可比矩陣為

J(E2)的相應特征方程可以寫成F(λ)=λ2-trJ(E2)λ+detJ(E2).其中,

顯然F(1)＞0,當α=α1=[4k(1+a(1-m)x2)2+2kaβ(1-m)2x2y2+kγβ(1-m)y2]/2x2(1+a(1-m)x2)2,此時F(-1)=0,trJ(E2)≠0.當trJ(E2)=0 時,有α=α*=[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}.當detJ(E2)=1 時,存在共軛復根,且|λ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}.

2 Flip 分岔

本節應用中心流形定理[14]研究了Flip 分岔的存在性.選擇α作為分岔參數,給參數α一個小擾動，并且令un=xn-x2,vn=yn-y2,模型(2)變為

將模型(3)在(un,vn,)=(0,0,0)時進行泰勒級數展開至二階

令

由于|λ1|=-1,|λ2|≠-1,通過計算,得到矩陣J(E2)的特征值分別為λ1=-1,λ2=3-αx2/k+αβ(1-m)2x2y2/(1+a(1-m)x2)2.

設矩陣

使用翻轉

則模型(3)變成以下形式

對模型(5)應用中心流形定理.假設Wc(0,0,0)表示模型(5)在=0 的小鄰域內的中心流形,則Wc(0,0,0)計算為

此外,對于限制于中心流形Wc(0,0,0)的映射,我們考慮映射G*

在(un,vn,)=(0,0,0)時定義兩個判別量l1≠0 和l2≠0.

因此,根據上述分析和文獻[14]中的定理3.1,我們給出了在正平衡點Flip 分岔存在的條件.

定理1如果l2≠0,模型(2)在E2(x2,y2)處存在Flip 分岔.如果l2＞0,則從E2(x2,y2)分岔的周期2 點是穩定的; 如果l2＜0,則從E2(x2,y2) 分岔的周期2 點是不穩定的.

3 Neimark-Sacker 分岔

本節選擇了c作為分岔參數來研究在E2(x2,y2)的Neimark-Sacker 分岔.令c=c1+,un=xn-x2,vn=yn-y2,其中||?1,模型(2)變為

將模型(7)在(un,vn)=(0,0)時進行泰勒級數展開至三階

其中c11,c12,···,c15,c21,c22,···,c25同模型(4)中的a11,a12,···,a15,a21,a22,···,a25,

模型(8)在原點(0,0)處線性化的特征方程為λ2-P()λ+Q(),其中,

使用翻轉

則模型(7)變成以下形式

其中,

定義第一Lyapunov 指數如下所示[14]

其中,

因此,根據上述分析和文獻[15]中的Neimark-Sacker 分岔存在理論,我們得到定理2,表明Neimark-Sacker分岔存在和方向的參數條件.

定理2如果滿足非退化線性條件且L≠0,模型(2)在E2(x2,y2)處存在Neimark-Sacker 分岔.若L＜0(或者L＞0),則當＞0(或者＜0) 時,吸引(或者排斥) 不變閉合曲線從E2(x2,y2)分岔.

4 數值模擬

例1選擇參數a=0.1,c=4,k=0.4,m=0.5,β=0.3,γ=0.2,通過數值計算得平衡點E2(0.338 983,1.269 07)和α=2.454 24,λ1=-1,λ2=0.926 375.模型(2)具有平衡點E2(x2,y2)和(a,m,k,α,β)∈M.我們只改變參數α以查看模型(2)的動力學行為變化.其中α∈[2.2,3.4].如圖1(a)所示,對于平衡點E2(x2,y2)，在α＜2.454 24 時是穩定的,當α=2.454 24 時失去了穩定性,當α＞2.454 24 時,存在Flip 分岔.此外,隨著α的增大,出現了一個混沌集合.圖2(a)是參數α∈[2.2,3.4] 時模型(2)的相圖.

圖1 模型(2)的分岔圖

圖2 模型(2)不同α 值、c 值的相圖

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