王立葉, 楊瑞智

(東北林業(yè)大學(xué) 理學(xué)院, 黑龍江 哈爾濱 150040)

水體富有營(yíng)養(yǎng)化引起藻類(lèi)和浮游植物大量繁殖,形成水華,造成水質(zhì)污染,控制甚至解決水華問(wèn)題不容忽視。除了生物法以及藻華防控技術(shù)等方法[1-4],學(xué)者們還構(gòu)建并研究了許多浮游植物模型并取得了豐富的研究成果。文獻(xiàn)[5]研究了物質(zhì)循環(huán)和狀態(tài)反饋控制下的營(yíng)養(yǎng)-浮游植物模型,文獻(xiàn)[6]分析了具趨化的浮游動(dòng)植物模型,文獻(xiàn)[7]探究了非退化擴(kuò)散和退化擴(kuò)散下的隨機(jī)浮游動(dòng)植物模型。本文考慮到由于營(yíng)養(yǎng)物質(zhì)和浮游植物在空間分布是不均勻的,并且會(huì)產(chǎn)生擴(kuò)散現(xiàn)象,故在文獻(xiàn)[8]已有的浮游植物模型的基礎(chǔ)上加入擴(kuò)散項(xiàng),研究了帶擴(kuò)散的浮游植物模型,其數(shù)學(xué)模型如下所示:

(1)

式中:x是空間變量;t是時(shí)間變量;d1和d2分別表示營(yíng)養(yǎng)物質(zhì)和浮游植物的擴(kuò)散系數(shù);Nx和Px分別表示N和P對(duì)x的一階導(dǎo)數(shù);Nx,x和Px,x分別表示N和P對(duì)x的二階導(dǎo)數(shù);N表示營(yíng)養(yǎng)物質(zhì)的濃度;P表示浮游植物的數(shù)量;a表示持續(xù)的外在營(yíng)養(yǎng)物質(zhì)的注入;b表示浮游植物對(duì)營(yíng)養(yǎng)物質(zhì)的最大吸收速率;c表示營(yíng)養(yǎng)物質(zhì)到浮游植物的最大轉(zhuǎn)變速率;d表示浮游植物的平均死亡速度;e表示營(yíng)養(yǎng)物質(zhì)的平均損失率;θ表示產(chǎn)毒浮游植物釋放有毒的化學(xué)物質(zhì)的速率;μ表示半飽和常量。

1 穩(wěn)定性分析

(2)

hαP3+(hβ-(1-sα))P2+(sβ+hμ2α)P-μ2(1-sα)=0。

(3)

P*為式(3)的解。

定理1 若1-sα≤0,則系統(tǒng)(2)沒(méi)有正平衡解;若0<1-sα≤hβ,則系統(tǒng)只有1個(gè)正平衡解;若hβ<1-sα,則系統(tǒng)有1個(gè)或3個(gè)正平衡解。

系統(tǒng)(2)在正平衡解E*處的雅可比矩陣為

λ2-Tr(J)λ+Det(J)=0。

(4)

式中:Tr(J)=b2-a1;Det(J)=a2b1-a1b2。則當(dāng)條件H1(b2-a1<0,a2b1-a1b2>0)滿(mǎn)足時(shí),式(4)有2個(gè)具有負(fù)實(shí)部的根,即系統(tǒng)(2)的正平衡解E*局部漸近穩(wěn)定。

2 Turing不穩(wěn)定性分析

利用反應(yīng)擴(kuò)散方程理論,可以得到反應(yīng)擴(kuò)散系統(tǒng)(2)在正平衡解E*處的特征方程為

令z=(n/l)2,可化為

λ2-Tnλ+Dn=0,(n=0,1,2,…)。

(5)

式中:Tn=-(d1+d2)z-a1+b2;Dn=d1d2z2+(a1d2-b2d1)z+a2b1-a1b2。

根據(jù)文獻(xiàn)[9],若想使E*不穩(wěn)定,式(5)至少有1個(gè)實(shí)部為正的特征根。

證明 令D(z)=0有2個(gè)根分別為z1,2且z1

3 Hopf分支分析

在這一部分選取a為分支參數(shù)。部分公式來(lái)源于文獻(xiàn)[10]。由Tn=0,得到關(guān)于參數(shù)a的一系列臨界值

取條件H3(①b2d1-a1d2<0;②b2d1-a1d2>0且Δ<0;③b2d1-a1d2>0,Δ>0且z?(z1,z2))。

引理1 當(dāng)a=a(n),(n=0,1,2,…)時(shí),Tn=0,若滿(mǎn)足條件H3中的任意一條時(shí),都有Dn>0,此時(shí)式(4)有1對(duì)純虛根。

定理3 若條件H3中任意一條成立,當(dāng)a=a(n)(n=0,1,2,…)時(shí),系統(tǒng)(2)在E*附近產(chǎn)生Hopf分支。

接下來(lái),將對(duì)Hopf分支的性質(zhì)進(jìn)行分析,為了方便,只給出了在a=a(0)時(shí)決定Hopf分支的性質(zhì)的參數(shù)計(jì)算公式。

系統(tǒng)(2)在E*附近的特征方程記為L(zhǎng)(a),L*(a)為算子L(a)的共軛算子。

可得

當(dāng)a=a(0)時(shí),可得

其中:

當(dāng)Re(c1(a(0)))<0(>0)時(shí),Hopf分支是向后的(向前的),并且分支周期解是局部漸近穩(wěn)定的(不穩(wěn)定的)。

4 數(shù)值模擬

令Ω=(0,4.5π),h=0.5,s=0.1,α=0.247,β=2.0,μ=0.5,c=1.0。計(jì)算得系統(tǒng)的正平衡解為(N*,P*)≈(2.084,0.76)。 在常微分方程系統(tǒng)中d1=d2=0,取a=2.098,其余參數(shù)保持一致,則有a1>b2,a2b1-a1b2=0.929>0,滿(mǎn)足條件H1,所以系統(tǒng)(2)的正平衡解為局部漸近穩(wěn)定,如圖1所示。選取參數(shù)d1=0.004,d2=0.002,偏微分狀態(tài)下系統(tǒng)(2)的正平衡解是局部漸近穩(wěn)定的,見(jiàn)圖2所示。

圖1 系統(tǒng)正平衡解局部漸近穩(wěn)定Fig.1 Locally asymptotically stable of the positive equilibrium solution of the system

(a)N(x,t)局部漸近穩(wěn)定(b) P(x,t)局部漸近穩(wěn)定

重新選取參數(shù)d1=3.0,d2=0.3,解得z1=0.805,z2=1.283。取n=5時(shí)有z∈(0.805,1.283)使Dn<0。即滿(mǎn)足產(chǎn)生Turing不穩(wěn)定的條件,所以系統(tǒng)(2)在E*處產(chǎn)生Turing不穩(wěn)定,如圖3所示。

(a) N(x,t)不穩(wěn)定圖像(b) P(x,t)不穩(wěn)定圖像

選取參數(shù)d1=0.001,d2=0.004,a=0.8,其余參數(shù)保持一致。計(jì)算可得Re(c1(a(0)))=-0.052 1<0,因此系統(tǒng)(2)在E*處有穩(wěn)定的周期解,如圖4所示。

(a) N(x,t)周期解圖像(b) P(x,t)周期解圖像

5 結(jié) 論

本文研究了模型的穩(wěn)定點(diǎn)存在的條件以及穩(wěn)定點(diǎn)穩(wěn)定性的分析。得出:當(dāng)條件H1成立時(shí),正平衡解是穩(wěn)定的;當(dāng)條件H1、H2同時(shí)成立時(shí),正平衡解是Turing不穩(wěn)定的;若條件H3中任意一條成立,當(dāng)a=a(n)(n=0,1,2,…)時(shí),系統(tǒng)(2)在E*附近產(chǎn)生Hopf分支。