張道祥, 孫光訊, 徐明麗, 陳金瓊, 周 文

(安徽師范大學 數(shù)學與統(tǒng)計學院, 安徽 蕪湖 241002)

0 引 言

文獻[1]用反應擴散模型說明了實際應用中的典型斑圖現(xiàn)象. 目前, 關(guān)于時空斑圖的研究已取得許多成果[2-7]： 文獻[2]研究了非線性捕食者收獲效應下時滯對Turing斑圖和螺旋波斑圖形成的影響; 文獻[3]討論了時滯對傳染病模型中Turing斑圖生成的影響; 文獻[4]研究了時滯反應擴散系統(tǒng)中螺旋波的形成和Hopf分支產(chǎn)生周期解的穩(wěn)定性及其方向. 文獻[8]研究了下列帶有Holling-Ⅱ型功能反應和雙曲死亡率的捕食者-食餌模型:

(1)

其中:U(t)和V(t)分別表示在t時刻食餌和捕食者的密度;r,K,a1,n,a2,δ均為正常數(shù), 分別表示食餌的內(nèi)稟增長率、 食餌的環(huán)境容納量、 捕食者對食餌的捕食率、 捕食者對食餌的最大攝入比率、 捕食者的增長率以及死亡率;,η是正常數(shù), 分別表示在浮游生物死亡的背景下由于水和自我保護導致的光衰減系數(shù). 文獻[8]討論了常微分系統(tǒng)的局部穩(wěn)定性和自擴散效應下Hopf分支的方向以及產(chǎn)生周期解的穩(wěn)定性. 本文考慮具有非線性食餌收獲效應的捕食者-食餌的反應擴散系統(tǒng)[9]:

(2)

在自然界中, 時滯會導致系統(tǒng)的周期解由穩(wěn)定變?yōu)椴环€(wěn)定[10]; 在神經(jīng)網(wǎng)絡(luò)中, 時滯會導致網(wǎng)絡(luò)產(chǎn)生振蕩和不穩(wěn)定現(xiàn)象, 從而破壞網(wǎng)絡(luò)的原有性能[11]; 在化學中, 時滯影響控制系統(tǒng)的穩(wěn)定性[12]等.因此, 研究具有時滯效應影響的數(shù)學模型具有一定的應用價值. 為簡化模型, 本文做如下變換:

考慮時滯效應, 引入上述變換可將系統(tǒng)(2)轉(zhuǎn)化為如下時滯偏微分系統(tǒng):

(3)

其中: (X,t)∈Ω×+,Ω?2為邊界光滑的有界區(qū)域;ν為?Ω上單位外法向量;τ是正常數(shù), 表示食餌種群孕育后代引起的時間滯后效應;u(X,t-τ)表示食餌在X處t時刻的密度.

對于系統(tǒng)(3), 文獻[9]研究了在自擴散效應下時滯對Hopf分支產(chǎn)生周期解的穩(wěn)定性及其方向的影響; 文獻[13]討論了當h=0,d12=d21=0時, 時滯效應下系統(tǒng)周期解的存在性. 基于生態(tài)學中物種空間分布的重要性, 本文考慮系統(tǒng)(3)中時滯效應和擴散效應對具有非線性食餌收獲效應系統(tǒng)空間斑圖生成的影響.

1 穩(wěn)定性與分支分析

1.1 線性穩(wěn)定性分析

首先, 考慮系統(tǒng)(3)對應的常微分系統(tǒng):

(4)

基于生態(tài)學意義, 本文只考慮系統(tǒng)(3)和系統(tǒng)(4)的正平衡點. 通過求解易知, 當h<ρ<β成立時, 系統(tǒng)(3)和系統(tǒng)(4)存在唯一的正平衡點E*=(u*,v*), 其中v*=u/(βγ)且u*滿足如下三次方程:

A0u3+A1u2+A2u+A3=0,

(5)

其中:

A0=βγ;A1=β2γ+βγρ+s-βγ;A2=β2γρ+βγh+ρs-β2γ-βγρ;A3=β2γ(h-ρ).

其次, 考慮系統(tǒng)(3)的穩(wěn)定性. 由文獻[14-15]知, 可假設(shè)τ是小量, 將u(X,t-τ)寫成如下形式:

(6)

將式(6)代入系統(tǒng)(3)中, 可得

(7)

其中R(u)=u/(β+u). 將R(u(X,t)-τ(?u(X,t)/?t))在(u,v)處Taylor展開, 且不考慮高階非線性部分, 則系統(tǒng)(7)轉(zhuǎn)化為

(8)

其中Ru(u)=dR(u)/du. 為簡便, 記

(9)

結(jié)合式(8)和式(9), 可得如下系統(tǒng):

(10)

下面對系統(tǒng)(10)的平衡點(u*,v*)進行線性穩(wěn)定性分析. 在系統(tǒng)(10)的平衡點(u*,v*)處做如下微擾:

u(X,t)=u*+δu(X,t),v(X,t)=v*+δv(X,t).

(11)

對系統(tǒng)(10)在平衡點(u*,v*)處進行線性化, 可得

(12)

其中:

在Neumann邊界條件下, 系統(tǒng)(10)的解有如下形式:

δu=δu*eλtei(kxx+kyy),δv=δv*eλtei(kxx+kyy),

(13)

(14)

(15)

將式(13)代入系統(tǒng)(10)中, 可得系統(tǒng)(10)在平衡點(u*,v*)處的特征方程:

λ2-T(k)λ+D(k)=0,

(16)

其中:

1.2 Hopf分支

若存在kH>0, 使得T(kH)=0,D(kH)>0,T(k)≠0,D(k)≠0, ?k≠kH, 則系統(tǒng)(10)發(fā)生Hopf分支. 此時矩陣A的特征值為純虛數(shù), 即tr(A)=T(kH)=0, 其中tr(A)表示A的跡.

假設(shè):

證明: 存在kj∈={0,1,2,…}, 使得當k=kj=kH時, tr(Aj)=0. 因此, 存在一個τ的臨界值τ(i), 對于i≥0, 有

1)q(τ)≠0, 對任意的τ∈I均成立;

綜上所述, 當τ=τc時, 系統(tǒng)(10)在平衡點(u*,v*)處發(fā)生了Hopf分支.

1.3 Turing分支

假設(shè):

由于時滯和擴散可導致系統(tǒng)(10)的Turing不穩(wěn)定. 下面考慮系統(tǒng)(10)發(fā)生Turing不穩(wěn)定的條件, 即Turing分支的條件. 由Turing斑圖形成的必要條件知, 對于系統(tǒng)(10)的常微分形式下平衡點(u*,v*)是穩(wěn)定的, 而對于系統(tǒng)(10)平衡點(u*,v*)是不穩(wěn)定的. 由上述穩(wěn)定性分析和Turing不穩(wěn)定的條件可知, 對于系統(tǒng)(10)的常微分方程在平衡點(u*,v*)的Jacobi矩陣J, 有tr(J)<0和det(J)>0, 且對任意的k>0, 總有tr(A)<0成立. 若要系統(tǒng)(10)在平衡點(u*,v*)處是不穩(wěn)定的, 則存在某個波數(shù)k>0, 使得det(A)<0. 由此可得如下結(jié)論:

記h為分支參數(shù). 通過計算可知Turing分支參數(shù)h的臨界值hT滿足如下Turing分支曲線:

圖1 系統(tǒng)(3)的分支Fig.1 Bifurcation of system (3)

通過上述線性穩(wěn)定性理論和分支理論[16-17], 可得h-τ關(guān)系如圖1所示.圖1中紫色曲線為正解性曲線, 由Turing分支曲線(藍線)和Hopf分支曲線(紅線), 可將整個區(qū)域分成4個參數(shù)空間, 其中: 區(qū)域Ⅰ表示Turing空間, 即在常微分系統(tǒng)中, 平衡點是穩(wěn)定的, 而在偏微分系統(tǒng)中, 平衡點是不穩(wěn)定的; 區(qū)域Ⅱ表示Turing-Hopf空間, 即在常微分系統(tǒng)和偏微分系統(tǒng)下, 系統(tǒng)的平衡點均不穩(wěn)定; 區(qū)域Ⅲ表示由Hopf分支引起的不穩(wěn)定, 即Hopf空間; 區(qū)域Ⅳ表示在常微分系統(tǒng)和偏微分系統(tǒng)下, 平衡點均是穩(wěn)定的.

2 數(shù)值模擬

下面用MATLAB軟件對系統(tǒng)(3)的空間斑圖進行數(shù)值模擬. 所有數(shù)值模擬均在離散格子200×200或400×400內(nèi)進行, 空間間隔Δh=1, 時間間隔Δt=0.1; 所有數(shù)值模擬均采用齊次Neumann邊界條件. 由于食餌和捕食者的空間分布類似, 所以數(shù)值模擬中僅考慮食餌空間斑圖的形成. 其中一些參數(shù)取值如下:

β=0.42,ρ=0.28,γ=1.2,s=0.8,d12=0.05.

(20)

2.1 區(qū)域Ⅰ中Turing斑圖的形成

考慮在圖1區(qū)域Ⅰ中Turing斑圖的形成以及時滯τ對Turing斑圖形成的影響. 初始條件為種群的初始隨機分布, 其他參數(shù)取值為:α=1.5,β=0.8,d11=1,d21=11,d22=11.5.圖2為在不同的食餌捕獲效應參數(shù)h下系統(tǒng)(3)的空間Turing斑圖, 其中(A),(B),(C)和(D),(E),(F)分別表示當τ=0,1.5時的情形.圖2(A)和(D)是當h=0.261時相同迭代步數(shù)下的Turing斑圖, 對比兩圖可見, 斑圖結(jié)構(gòu)一致, 均為點狀, 但圖2(D)相對圖2(A)點狀斑圖形成更快. 當h值增至0.275和0.296時, 同理可見圖2(B)和(E)、
圖2(C)和(F)均為點狀和條狀共存的斑圖, 只是圖2(B)和(E)中點狀多于條狀,圖2(C)和(F)中條狀多于點狀, 且斑圖的結(jié)構(gòu)未發(fā)生改變, 與文獻[2-3]中結(jié)論一致.

2.2 區(qū)域Ⅱ中的空間斑圖

圖3為在常微分系統(tǒng)和偏微分系統(tǒng)下平衡點均不穩(wěn)定時, 食餌u隨時間演化的空間斑圖. 初始條件及各參數(shù)取值為:u(X,t0)=u*-0.000 2×(x-100),v(X,t0)=v*-0.000 3×(y-100),t0∈[-τ,0][7],α=1.5,h=0.1,τc=1.458 5,τ=1.65>τc,d11=0.1,d21=0.8,d22=0.95.圖3(A)為當t=0時的斑圖, 隨著時間的演化, 出現(xiàn)規(guī)則的條狀斑圖和點狀斑圖共存的情形(圖3(B)～(E)), 不同的是圖3(B)～(D)中條狀優(yōu)于點狀, 而圖3(E)中點狀優(yōu)于條狀, 最后點狀斑圖占據(jù)整個區(qū)域(圖3(F)), 系統(tǒng)的動力學行為不再發(fā)生改變.

圖3 當h=0.1, τ=1.65時, 食餌u隨時間演化的空間斑圖Fig.3 Spatial patterns of time evolution of prey u when h=0.1, τ=1.65

2.3 區(qū)域Ⅲ中螺旋波斑圖

考慮圖1的區(qū)域Ⅲ內(nèi)螺旋波的形成. 初始條件[4]為u(X,t0)=u*-0.000 2×(x-200),v(X,t0)=v*-0.000 3×(y-200),t0∈[-τ,0].

2.3.1 時滯對螺旋波斑圖形成的影響 圖4為系統(tǒng)(3)中食餌u的螺旋波斑圖隨時間的演化情況. 各參數(shù)取值為:α=2.15,h=0.05,d11=0.2,d21=0.03,d22=0.6, 其他參數(shù)由式(20)確定, 此時τc=1.524 9.圖4(A)為當t=0時的斑圖, 起初螺旋波斑圖是規(guī)則的(圖4(B)～(D)), 隨著時間的演化, 由螺旋波中心開始, 其規(guī)則結(jié)構(gòu)依次消失, 逐漸演化成不規(guī)則的螺旋波斑圖(圖4(E)～(G)), 最后不規(guī)則的螺旋波斑圖占據(jù)整個區(qū)域(圖4(H)). 由圖4可見, 不規(guī)則的螺旋波斑圖由規(guī)則的螺旋波斑圖隨時間的演化而逐漸形成.

圖4 當τ=1.9時, 食餌u隨時間演化的空間螺旋波斑圖Fig.4 Spatial spiral patterns of time evolution of prey u when τ=1.9

圖5 u(t)的時間演化曲線Fig.5 Time evolution curve of u(t)

圖5為食餌u隨時間演化的密度曲線. 由圖5可見, 隨著時間的變化, 食餌u的密度呈振蕩的變化趨勢. 因此, 在足夠長的時間內(nèi), 食餌u的空間斑圖向不規(guī)則的趨勢變化.圖6為在不同時滯τ下食餌u的空間螺旋波斑圖. 各參數(shù)取值為:α=2.15,h=0.05,d11=0.2,d21=0.4,d22=0.6, 其他參數(shù)由式(20)確定. 由圖6可見: 當τ=2.1時, 規(guī)則的螺旋波斑圖占據(jù)整個區(qū)域(圖6(A)); 隨著τ的增大, 中間部分仍是規(guī)則的, 而邊緣出現(xiàn)不規(guī)則的螺旋波斑圖; 且隨著時滯的增大, 中間部分規(guī)則的螺旋波變得越來越細, 外圍部分越來越混亂(圖6(B)～(D)).

2.3.2 非線性收獲效應對螺旋波斑圖形成的影響

圖7為在同一時刻對系統(tǒng)(3)隨h變化的螺旋波斑圖. 各參數(shù)取值為:α=2.15,d11=0.2,d21=0.15,d22=0.85,τ=1.9, 其他參數(shù)由式(20)確定. 由圖7可見: 當h=0.05時, 系統(tǒng)(3)產(chǎn)生的螺旋波斑圖是規(guī)則的(圖7(A)); 隨著h增至0.08和0.11, 開始出現(xiàn)不規(guī)則的螺旋波(圖7(B),(C)), 且由邊緣處依次向中間延伸破碎; 當h增至0.16時, 完全不規(guī)則的螺旋波斑圖占據(jù)整個區(qū)域(圖7(D)); 隨著收獲效應參數(shù)h的增大, 螺旋波的破碎進程加快. 因此在實際應用中, 人類適當?shù)氖斋@可對物種進行保護[18].

圖6 當時間t=10 000時, 食餌u隨τ變化的空間螺旋波斑圖Fig.6 Spatial spiral patterns of prey u with different τ when t=10 000

圖7 當時間t=10 000時, 食餌u隨h變化的空間螺旋波斑圖Fig.7 Spatial spiral patterns of prey u with different h when t=10 000

圖8 在時間t=10 000時, 食餌u隨d21變化的空間螺旋波斑圖Fig.8 Spatial spiral patterns of prey u with different d21 when t=10 000

2.3.3 交叉擴散效應對螺旋波斑圖形成的影響 圖8為在同一時刻系統(tǒng)(3)隨d21變化的螺旋波斑圖. 各參數(shù)取值為:α=2.15,h=0.05,d11=0.2,d22=0.6,τ=1.9, 其他參數(shù)由式(20)確定. 由圖8可見: 當d21=0.42時, 系統(tǒng)(3)出現(xiàn)了規(guī)則的螺旋波(圖8(A)); 隨著d21減至0.29和0.12, 螺旋波斑圖逐漸變得不規(guī)則(圖8(B),(C)), 且由邊緣處依次向中間延伸破碎; 當d21減至0.03時, 系統(tǒng)(3)出現(xiàn)了完全不規(guī)則的螺旋波斑圖(圖8(D)); 隨著擴散系數(shù)d21的減小, 螺旋波的破碎進程加快.

綜上, 基于時滯對生態(tài)系統(tǒng)動力學行為的影響, 本文從理論和數(shù)值兩方面研究了一類具有時滯反應擴散的非線性收獲效應的捕食者-食餌模型的空間斑圖生成問題, 可得如下結(jié)論: 1) 數(shù)值結(jié)果表明, 時滯未改變系統(tǒng)(3)的Turing斑圖結(jié)構(gòu), 可能會影響Turing斑圖生成的進程; 2) 理論與數(shù)值結(jié)果表明, 時滯和非線性收獲效應以及交叉擴散效應均影響了螺旋波斑圖的生成和破碎; 3) 理論與數(shù)值結(jié)果表明, 系統(tǒng)(3)具有豐富的動力學行為.