一类具有双时滞效应的氮–磷–浮游植物模型的动力学研究

Dynamic Analysis of a Nitrogen-Phosphorus-Plankton Model with Two Time Delays

摘要: 在水体富营养化的预防和控制的研究背景下，本论文构建了一类具有双时滞效应的氮–磷–浮游植物动力学模型，并对其相关动力学性质进行理论分析和数值仿真，解析出模型正平衡点具有局部渐近稳定性和发生Hopf分支的阈值条件，模拟出营养及浮游植物密度和Hopf分岔的动态演化过程。本论文的研究结果有利于从动力学的角度揭示营养与浮游植物之间的互作机制，有助于深入分析时滞效应对海洋生态系统的影响，为进一步预防和控制水体富营养化提供一定的理论支撑。

Abstract: Under the framework of water eutrophication prevention and control research, a nitro-gen-phosphorus-phytoplankton model with two delays will be investigated mathematically and numerically. Mathematical works have comprehensively explored the locally asymptotic stability of positive equilibrium and the threshold conditions of occurring Hopf bifurcation. Numerical si-mulation works have depicted the dynamic evolution process of nutrient and phytoplankton den-sity and Hopf bifurcation. These results provide a great help for the interaction between nutrition and phytoplankton in dynamics, and are helpful to deeply understand how the delay affects the dynamic trend of the marine ecosystem. Furthermore, they can provide certain theoretical support for further prevention and control of eutrophication of water bodies as well.

1. 引言

富营养化是指生物所需的氮、磷等无机营养物质大量进入湖泊、河口、海湾等相对封闭、水流缓慢的水体，在适宜的外界环境因素综合作用下，引起藻类及其它浮游生物迅速繁殖，水体溶解氧量下降，水质恶化，鱼类及其它水生生物大量死亡的现象 [1] [2] [3] 。此外，人类活动的影响可加剧这一过程，特别是在现代生产和生活中，人类对环境资源的开发利用日益频繁，工农业发展迅速，大量的营养物质进入水体并在其中积累，导致富营养化在短期内出现。目前，人类活动导致的水体富营养化已成为当今世界水污染治理的难题，是全球范围内普遍存在的环境问题 [4] 。

如今，许多实验生态学家和数学生态学家一直在关注着浮游植物生长的动态行为，以便找出如何控制富营养化，如何预测藻类爆发，以及如何模拟藻类传播趋势等等 [5] [6] [7] 。他们的研究表明，光照，温度，矿物元素等都是影响浮游植物生长，生物量和群落分布的关键因素，这一工作为准确预防藻类爆发提供了良好的理论依据，对探索自然生态系统中的富营养化问题至关重要。其中有些学者也尝试着运用数学模型来研究富营养化问题 [8] [9] [10] 。他们的工作表明，这些数学模型不但可以近似地刻画生物种群的动态演化进程，还可以解释种群之间的增长和相互作用等。近年来越来越多的学者开始对生物种群动态的数学模型进行研究 [11] [12] [13] 。Mei and Zhang [11] 研究了一种非局部的反应扩散–平流系统，模拟了多种竞争性浮游植物物种的生长。Wang et al. [12] 考虑到收获和毒素释放对浮游植物产生的影响，建立了一类带有时滞效应和选择性收获的产毒浮游植物–浮游动物模型。Wang et al. [13] 构建了一类具有时滞效应的营养–浮游生物模型。他们的研究成果都具有重要意义，为解决水体富营养化问题提供了不同的视角。文献 [14] [15] [16] 也详细研究了带有时滞效应的生物种群动态的数学模型，它们的研究成果为我们构建模型提供了一定的理论依据，从而保证我们所构建的模型具有一定的实际意义。

Li and Liu [17] 研究了一类具有时滞效应的浮游生物模型，解析出系统正平衡点全局稳定性和发生Hopf分岔的阈值条件，这一方法很值得我们借鉴。Martin and Ruan [18] 考虑了收获率和时间延迟的综合影响，构建了一类带有时滞效应和选择性收获的浮游植物–浮游动物模型，他们发现时滞因素可能会导致系统发生失稳，但是当收获值处于临界收获水平时，收获率对系统的平衡有稳定作用，这项工作可以为研究收获率和时间延迟之间的动态机制提供一个良好的框架。Shi et al. [19] 研究了一类具有分布时滞和Dirichlet边界条件的反应–扩散模型，他们的研究成果具有非常重要的意义，为我们研究营养循环如何影响浮游植物生长的问题提供了一套系统的研究方法。Shi and Jiang [20] 研究了一类具有双时滞效应的浮游动物–浮游植物模型，这一研究为我们以后研究多时滞模型提供了一套完整的理论和模拟方法等。他们的研究成果对预防和控制水体的富营养化具有一定的指导意义。

基于以上讨论和Deng et al. [21] 研究的相关结果，本文构建了一类具有时滞效应的氮–磷–浮游植物模型，其可以表示如下：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=I_1-b_{1}x-\frac{a_{1}xz}{c_1+x}\\ \frac{\text{d}y}{\text{d}t}=I_2-b_{2}y-\frac{a_{2}yz}{c_2+y}\\ \frac{\text{d}z}{\text{d}t}=\frac{e_1{a}_{1}x\left(t-{\tau }_1\right)z}{c_1+x\left(t-{\tau }_1\right)}+\frac{e_2{a}_{2}y\left(t-{\tau }_2\right)z}{c_2+y\left(t-{\tau }_2\right)}-mz\end{array}$ , (1)

其中x，y分别是总氮(mg/L)和总磷(mg/L)在t时刻的密度，z是浮游植物种群在t时刻的生物量(mg/L)， ${\tau }_1$ 是浮游植物吸收氮引起的时滞项， ${\tau }_2$ 是浮游植物吸收磷引起的时滞项。系统(1)中的其他参数与 [21] 中的参数相同。从生物学的角度来看，我们只关注闭合第一象限系统(1)的相关动力学行为，因此系统(1)满足以下初始条件：

$\begin{array}{l}x\left(t\right)={\varphi }_1\left(\theta \right)\ge 0\\ y\left(t\right)={\varphi }_2\left(\theta \right)\ge 0\\ z\left(t\right)={\varphi }_3\left(\theta \right)\ge 0\end{array}\}\theta \in \left[-\tau ,0\right]$ , (2)

其中 ${\varphi }_r\in C\left(\left[-\tau ,0\right]\to R_{+}^3\right)$ ， $\left(r=1,2,3\right)$ ， $R_{+}^3=\left\{\left(x,y,z\right)|x>0,y>0,z>0\right\}$ 。

2. 局部稳定性和Hopf分支

现在主要讨论系统(1)正平衡点的稳定性和Hopf分支的存在性。

首先对系统(1)进行换元变形： $u\left(t\right)=x\left(t\right)-x^{\text{*}}$ ， $v\left(t\right)=y\left(t\right)-y^{\text{*}}$ ， $w\left(t\right)=z\left(t\right)-z^{*}$ ，将正平衡点 $E_2=\left(x^{\text{*}},y^{*},z^{*}\right)$ 转化到原点，然后将系统(1)线性化处理，则系统(1)可以转化为：

$\frac{\text{d}U}{\text{d}t}=MU\left(t\right)+NU\left(t-{\tau }_1\right)+OU\left(t-{\tau }_2\right)$ , (3)

其中：

$U\left(t\right)={\left[u,v,w\right]}^{\text{T}}$ , $M={\left(m_{ij}\right)}_{3\times 3}$ , $N={\left(n_{ij}\right)}_{3\times 3}$ , $O={\left(o_{ij}\right)}_{3\times 3}$ ,

$m_{11}=-b_1-\frac{a_1{c}_1{z}^{*}}{{\left(c_1+x^{*}\right)}^2}$ , $m_{13}=-\frac{a_1{x}^{*}}{c_1+x^{*}}$ ,

$m_{22}=-b_2-\frac{a_2{c}_2{z}^{*}}{{\left(c_2+y^{*}\right)}^2}$ , $m_{23}=-\frac{a_2{y}^{*}}{c_2+y^{*}}$ ,

$n_{31}=\frac{e_1{a}_1{c}_1{z}^{*}}{{\left(c_1+x^{*}\right)}^2}$ , $n_{32}=\frac{e_2{a}_2{c}_2{z}^{*}}{{\left(c_2+y^{*}\right)}^2}$ ,

其他 $m_{ij}=0$ ， $n_{ij}=0$ ， $o_{ij}=0$ 。

由于 $\left|\lambda E-A\right|=0$ ，可得：

$\left|\begin{array}{ccc}\lambda -a_{11}& 0& -a_{13}\\ 0& \lambda -a_{22}& -a_{23}\\ -a_{31}\cdot {\text{e}}^{-\lambda {\tau }_1}& -a_{32}\cdot {\text{e}}^{-\lambda {\tau }_2}& \lambda \end{array}\right|=0$ , (4)

其中：

$a_{11}=-b_1-\frac{a_1{c}_1{z}^{*}}{{\left(c_1+x^{*}\right)}^2}$ , $a_{13}=-\frac{a_1{x}^{*}}{c_1+x^{*}}$ ,

$a_{22}=-b_2-\frac{a_2{c}_2{z}^{*}}{{\left(c_2+y^{*}\right)}^2}$ , $a_{23}=-\frac{a_2{y}^{*}}{c_2+y^{*}}$ ,

$a_{31}=\frac{e_1{a}_1{c}_1{z}^{*}}{{\left(c_1+x^{*}\right)}^2}$ , $a_{32}=\frac{e_2{a}_2{c}_2{z}^{*}}{{\left(c_2+y^{*}\right)}^2}$ ,

系统(3)对应的特征方程如下：

${\lambda }^3+A{\lambda }^2+B\lambda +\left(C+D\lambda \right){\text{e}}^{-\lambda {\tau }_1}+\left(E+F\lambda \right){\text{e}}^{-\lambda {\tau }_2}=0$ , (5)

其中：

$A=-\left(a_{11}+a_{22}\right)$ , $B=a_{11}{a}_{22}$ , $B=a_{11}{a}_{22}$ ,

$D=-a_{13}{a}_{31}$ , $D=-a_{13}{a}_{31}$ , $D=-a_{13}{a}_{31}$ .

下面分五种情况进行讨论：

第一种情况： ${\tau }_1={\tau }_2=0$ 。

定理2.1： ${\tau }_1={\tau }_2=0$ 时，满足条件 $\left(H_1\right)$ ： $A_{11}>0$ ， $A_{11}{A}_{12}>A_{13}$ 时，系统正平衡点 $E_2=\left(x^{\text{*}},y^{*},z^{*}\right)$ 是局部稳定的。

证明：当 ${\tau }_1={\tau }_2=0$ 时，方程(5)可以转化为：

${\lambda }^3+A_{11}{\lambda }^2+A_{12}\lambda +A_{13}=0$ , (6)

其中：

$A_{11}=A$ , $A_{12}=B+D+F$ , $A_{13}=C+E$ .

根据Routh-Hurwitz判据，满足条件 $\left(H_1\right)$ ： $A_{11}>0$ ， $A_{11}{A}_{12}>A_{13}$ 时， $E_2=\left(x^{\text{*}},y^{*},z^{*}\right)$ 是局部渐近稳定的。

第二种情况： ${\tau }_1>0,{\tau }_2=0$ 。

定理2.2： ${\tau }_1>0,{\tau }_2=0$ 时，满足条件：

1) $\left(H_{21}\right)$ ：文献 [22] 中定理(2.1)情况(i)，(iii)成立；

2) $\left(H_{22}\right)$ ： ${f^{\prime }}_1\left({\omega }_{10}^2\right)\ne 0$ ， $A_{21}^2-2A_{22}>0$ ， $A_{22}^2-2A_{21}{A}_{23}>0$ 。

系统正平衡点 $E_2=\left(x^{\text{*}},y^{*},z^{*}\right)$ 在 ${\tau }_1\in \left[0,{\tau }_{10}\right)$ 是局部渐近稳定的，并且当 ${\tau }_1={\tau }_{10}$ 时出现Hopf分支。

证明：当 ${\tau }_1>0,{\tau }_2=0$ 时，方程(5)可以转化为：

${\lambda }^3+A_{21}{\lambda }^2+A_{22}\lambda +A_{23}+\left(B_{21}\lambda +B_{22}\right){\text{e}}^{-\lambda {\tau }_1}=0$ , (7)

其中：

$A_{21}=-\left(a_{11}+a_{22}\right)$ , $A_{22}=a_{11}{a}_{22}-a_{23}{a}_{32}$ , $A_{23}=a_{23}{a}_{32}{a}_{11}$ ,

$B_{21}=-a_{13}{a}_{31}$ , $B_{22}=a_{13}{a}_{31}{a}_{22}$ ,

假设 $\lambda =i{\omega }_1\left({\omega }_1>0\right)$ 是方程(7)的纯虚根，将其代入并且分离实部和虚部可得：

$\{\begin{array}{l}-{\omega }_1^3+A_{22}{\omega }_1+B_{21}{\omega }_1\cos \left({\omega }_1{\tau }_1\right)-B_{22}\sin \left({\omega }_1{\tau }_1\right)=0\\ -A_{21}{\omega }_1^2+A_{23}+B_{21}{\omega }_1\sin \left({\omega }_1{\tau }_1\right)+B_{22}\cos \left({\omega }_1{\tau }_1\right)=0\end{array}$ , (8)

将(8)中两个方程平方然后相加，可得：

${\omega }_1^6+E_{21}{\omega }_1^4+E_{22}{\omega }_1^2+E_{23}=0$ ,(9)

其中：

$E_{21}=A_{21}^2-2A_{22}$ , $E_{22}=A_{22}^2-B_{21}^2-2A_{21}{A}_{23}$ , $E_{23}=A_{23}^2-B_{22}^2$ .

假设 $v_1={\omega }_1^2$ ，则方程(9)可以转化为：

$v_1^3+E_{21}{v}_1^2+E_{22}{v}_1+E_{23}=0$ . (10)

接下来定义一个函数：

$f_1\left(v_1\right)=v_1^3+E_{21}{v}_1^2+E_{22}{v}_1+E_{23}$ . (11)

假设 $\left(H_{21}\right)$ ：文献 [22] 中定理(2.1)情况(i)，(iii)成立，则 $f_1\left(v_1\right)$ 至少有一个正根，不失一般性，假设方程(11)的根为 $v_{11},v_{12},v_{13}$ ，则有 ${\omega }_{1k}=\sqrt{v_{1k}},\left(k=1,2,3\right)$

由(8)可得：

$\cos \left({\omega }_1{\tau }_1\right)=\frac{B_{21}{\omega }_1^4+\left(B_{22}{A}_{21}-B_{21}{A}_{22}\right){\omega }_1^2-B_{22}{A}_{23}}{B_{21}^2{\omega }_1^2+B_{22}^2}$ , (12)

$\sin \left({\omega }_1{\tau }_1\right)=\frac{\left(B_{21}{A}_{21}-B_{22}\right){\omega }_1^3+\left(B_{22}{A}_{22}-B_{21}{A}_{23}\right){\omega }_1}{B_{21}^2{\omega }_1^2+B_{22}^2}$ , (13)

由(12)和(13) 可得：

${\phi }_1=\{\begin{array}{l}2\text{π}-\arccos \left[\frac{B_{21}{\omega }_1^4+\left(B_{22}{A}_{21}-B_{21}{A}_{22}\right){\omega }_1^2-B_{22}{A}_{23}}{B_{21}^2{\omega }_1^2+B_{22}^2}\right],\text{if}\left(13\right)<0\\ \cos \left[\frac{B_{21}{\omega }_1^4+\left(B_{22}{A}_{21}-B_{21}{A}_{22}\right){\omega }_1^2-B_{22}{A}_{23}}{B_{21}^2{\omega }_1^2+B_{22}^2}\right],\text{if}\left(13\right)\ge 0\end{array}$ , (14)

${\tau }_{1k}^j=\frac{1}{{\omega }_{1k}}\left({\phi }_1+2j\text{π}\right),\left(j=0,1,2,\cdots ,k=1,2,3\right)$ ,(15)

定义 ${\tau }_{10}=\min \left({\tau }_{1k}^{\left(0\right)}|k=1,2,3\right)$ ，令 $\lambda \left({\tau }_1\right)=\alpha \left({\tau }_1\right)+i{\omega }_1\left({\tau }_1\right)$ 是特征方程(7)的一个根，且近于 ${\tau }_1={\tau }_{10}$ ，并且满足 $\alpha \left({\tau }_{10}\right)=0$ ， ${\omega }_1\left({\tau }_{10}\right)={\omega }_{10}$ 。对关于 ${\tau }_1$ 的特征方程(7)两边进行求导，得到：

${\left(\frac{\text{d}\lambda }{\text{d}{\tau }_1}\right)}^{-1}=\frac{3{\lambda }^2+2A_{21}\lambda +A_{22}}{-\lambda \left({\lambda }^3+A_{21}{\lambda }^2+A_{22}\lambda +A_{23}\right)}+\frac{B_{21}}{B_{21}{\lambda }^2+B_{22}\lambda }-\frac{{\tau }_1}{\lambda }$ , (16)

则

$\begin{array}{c}{\left(\frac{\text{d}\mathrm{Re}\left(\lambda \right)}{\text{d}{\tau }_1}\right)}_{\lambda =i{\omega }_{10}}^{-1}=\frac{3{\omega }_{10}^6+2\left(A_{21}^2-2A_{22}\right){\omega }_{10}^4+\left(A_{22}^2-B_{21}^2-2A_{21}{A}_{23}\right){\omega }_{10}^2}{{\omega }_{10}^8+\left(A_{21}^2-2A_{22}\right){\omega }_{10}^6+\left(A_{22}^2-2A_{21}{A}_{23}\right){\omega }_{10}^4+A_{23}^2{\omega }_{10}^2}\\ =\frac{1}{{\omega }_{10}^2}\left(\frac{{f^{\prime }}_1\left({\omega }_{10}^2\right)}{{\omega }_{10}^6+\left(A_{21}^2-2A_{22}\right){\omega }_{10}^4+\left(A_{22}^2-2A_{21}{A}_{23}\right){\omega }_{10}^2+A_{23}^2}\right),\end{array}$ (17)

得知当 $\left(H_{22}\right)$ ： ${f^{\prime }}_1\left({\omega }_{10}^2\right)\ne 0$ ， $A_{21}^2-2A_{22}>0$ ， $A_{22}^2-2A_{21}{A}_{23}>0$ 成立时， $\mathrm{Re}{\left(\frac{\text{d}\lambda }{\text{d}{\tau }_1}\right)}_{\lambda =\text{i}{\omega }_{10}}^{-1}\ne 0$ 。证毕。

第三种情况： ${\tau }_2>0,{\tau }_1=0$ 。

定理2.3： ${\tau }_2>0,{\tau }_1=0$ 时，满足条件：

1) $\left(H_{21}\right)$ ：同定理2.2；

2) $\left(H_{32}\right)$ ： ${f^{\prime }}_2\left({\omega }_{20}^2\right)\ne 0$ ， $A_{31}^2-2A_{32}>0$ ， $A_{32}^2-2A_{31}{A}_{33}>0$ 。

系统正平衡点 $E_2=\left(x^{\text{*}},y^{*},z^{*}\right)$ 在 ${\tau }_2\in \left[0,{\tau }_{20}\right)$ 是局部渐近稳定的，并且当 ${\tau }_2={\tau }_{20}$ 时出现Hopf分支。

证明：当 ${\tau }_2>0,{\tau }_1=0$ 时，方程(5)可以转化为：

${\lambda }^3+A_{31}{\lambda }^2+A_{32}\lambda +A_{33}+\left(B_{31}\lambda +B_{32}\right){\text{e}}^{-\lambda {\tau }_2}=0$ , (18)

其中：

$A_{31}=-\left(a_{11}+a_{22}\right)$ , $A_{32}=a_{11}{a}_{22}-a_{13}{a}_{31}$ , $A_{33}=a_{13}{a}_{31}{a}_{22}$ , $B_{31}=-a_{23}{a}_{32}$ , $B_{32}=a_{23}{a}_{32}{a}_{11}$ ,

假设 $\lambda =i{\omega }_2\left({\omega }_2>0\right)$ 是方程(18)的根，将其代入并且分离实部和虚部可得：

$\{\begin{array}{c}-{\omega }_2^3+A_{32}{\omega }_2+B_{31}{\omega }_2\cos \left({\omega }_2{\tau }_2\right)-B_{32}\sin \left({\omega }_2{\tau }_2\right)=0\\ -A_{31}{\omega }_2^2+A_{33}+B_{31}{\omega }_2\sin \left({\omega }_2{\tau }_2\right)+B_{32}\cos \left({\omega }_2{\tau }_2\right)=0\end{array}$ , (19)

将(19)中两个方程平方然后相加，可得：

${\omega }_2^6+E_{31}{\omega }_2^4+E_{32}{\omega }_2^2+E_{33}=0$ , (20)

其中：

$E_{31}=A_{31}^2-2A_{32}$ , $E_{32}=A_{32}^2-B_{31}^2-2A_{31}{A}_{33}$ , $E_{33}=A_{33}^2-B_{32}^2$ .

假设 $v_2={\omega }_2^2$ ，则方程(20)可以转化为：

$v_2^3+E_{31}{v}_2^2+E_{32}{v}_2+E_{33}=0$ . (21)

定义一个函数：

$f_2\left(v_2\right)=v_2^3+E_{31}{v}_2^2+E_{32}{v}_2+E_{33}$ . (22)

假设 $\left(H_{21}\right)$ 成立，则方程(22)至少有一个正根，不失一般性，将方程(22)的正根记为 $v_{21},v_{22},v_{23}$ ，则有 ${\omega }_{2k}=\sqrt{v_{2k}},\left(k=1,2,3\right)$ 。根据(19)可得：

$\cos \left({\omega }_2{\tau }_2\right)=\frac{B_{31}{\omega }_2^4+\left(B_{32}{A}_{31}-B_{31}{A}_{32}\right){\omega }_2^2-B_{32}{A}_{33}}{B_{31}^2{\omega }_2^2+B_{32}^2}$ , (23)

$\sin \left({\omega }_2{\tau }_2\right)=\frac{\left(B_{31}{A}_{31}-B_{32}\right){\omega }_2^3+\left(B_{32}{A}_{32}-B_{31}{A}_{33}\right){\omega }_2}{B_{31}^2{\omega }_2^2+B_{32}^2}$ , (24)

由(23)和(24)可得：

${\phi }_2=\{\begin{array}{l}2\text{π}-\arccos \left[\frac{B_{31}{\omega }_2^4+\left(B_{32}{A}_{31}-B_{31}{A}_{32}\right){\omega }_2^2-B_{32}{A}_{33}}{B_{31}^2{\omega }_2^2+B_{32}^2}\right],\text{if}\left(24\right)<0\\ \cos \left[\frac{B_{31}{\omega }_2^4+\left(B_{32}{A}_{31}-B_{31}{A}_{32}\right){\omega }_2^2-B_{32}{A}_{33}}{B_{31}^2{\omega }_2^2+B_{32}^2}\right],\text{if}\left(24\right)\ge 0\end{array}$ ,

${\tau }_{2k}^j=\frac{1}{{\omega }_{2k}}\left({\phi }_2+2j\text{π}\right),\left(j=0,1,2,\cdots ,k=1,2,3\right)$ . (25)

定义 ${\tau }_{20}=\min \left({\tau }_{2k}^{\left(0\right)}|k=1,2,3\right)$ ，令 $\lambda \left({\tau }_2\right)=\alpha \left({\tau }_2\right)+i{\omega }_2\left({\tau }_2\right)$ 是特征方程(18)的一个根，且近于 ${\tau }_2={\tau }_{20}$ ，并且满足 $\alpha \left({\tau }_{20}\right)=0$ ， ${\omega }_2\left({\tau }_{20}\right)={\omega }_{20}$ 。对关于 ${\tau }_2$ 的特征方程(18)两边进行求导，得到：

${\left(\frac{\text{d}\lambda }{\text{d}{\tau }_2}\right)}^{-1}=\frac{3{\lambda }^2+2A_{31}\lambda +A_{32}}{-\lambda \left({\lambda }^3+A_{31}{\lambda }^2+A_{32}\lambda +A_{33}\right)}+\frac{B_{31}}{B_{31}{\lambda }^2+B_{32}\lambda }-\frac{{\tau }_2}{\lambda }$ , (26)

故

${\left(\frac{\text{d}\mathrm{Re}\left(\lambda \right)}{\text{d}{\tau }_2}\right)}_{\lambda =i{\omega }_{20}}^{-1}=\frac{1}{{\omega }_{20}^2}\left(\frac{{f^{\prime }}_2\left({\omega }_{20}^2\right)}{{\omega }_{20}^6+\left(A_{31}^2-2A_{32}\right){\omega }_{20}^4+\left(A_{32}^2-2A_{31}{A}_{33}\right){\omega }_{20}^2+A_{33}^2}\right)$ , (27)

得知当 $\left(H_{32}\right)$ ： ${f^{\prime }}_2\left({\omega }_{20}^2\right)\ne 0$ ， $A_{31}^2-2A_{32}>0$ ， $A_{32}^2-2A_{31}{A}_{33}>0$ ，则 $\mathrm{Re}{\left(\frac{\text{d}\lambda }{\text{d}{\tau }_2}\right)}_{\lambda =\text{i}{\omega }_{20}}^{-1}\ne 0$ 。证毕。

第四种情况： ${\tau }_1\in \left(0,{\tau }_{10}\right),{\tau }_2>0$ 。

定理2.4： ${\tau }_1\in \left(0,{\tau }_{10}\right),{\tau }_2>0$ 时，满足条件：

1) $\left(H_{41}\right)$ ：(29)存在有限个正根；

2) $\left(H_{42}\right)$ ： $\left(A_1{C}_1+B_1{D}_1\right)\left({\omega }_{20}^{\ast }{}^2+a_{11}^2\right)-\left(A_1^2+B_1^2\right)\ne 0$ 。

系统正平衡点 $E_2=\left(x^{\text{*}},y^{*},z^{*}\right)$ 在 ${\tau }_2\in \left[0,{\tau }_{20}^{*}\right)$ 是局部渐近稳定的，并且当 ${\tau }_2={\tau }_{20}^{*}$ 时出现Hopf分支。

证明：当 ${\tau }_1\in \left(0,{\tau }_{10}\right),{\tau }_2>0$ 时，将 ${\tau }_2$ 作为参数，假设 $\lambda =i{\omega }_2^{*}\left({\omega }_2^{*}>0\right)$ 是方程(5)的根，将其代入并且分离实部和虚部可得：

$\{\begin{array}{l}-{\omega }_2^{\ast }{}^3+a_{11}{a}_{22}{\omega }_2^{\ast }-\left[a_{13}{a}_{31}{a}_{22}\sin \left({\omega }_2^{\ast }{\tau }_1\right)+a_{13}{a}_{31}{\omega }_2^{\ast }\cos \left({\omega }_2^{\ast }{\tau }_1\right)\right]=a_{23}{a}_{32}{a}_{11}\sin \left({\omega }_2^{\ast }{\tau }_2\right)+a_{23}{a}_{32}{\omega }_2^{\ast }\cos \left({\omega }_2^{\ast }{\tau }_2\right)\\ \left(a_{11}+a_{22}\right){\omega }_2^{\ast }{}^2++a_{13}{a}_{31}{a}_{22}\cos \left({\omega }_2^{\ast }{\tau }_1\right)-a_{13}{a}_{31}{\omega }_2^{\ast }\sin \left({\omega }_2^{\ast }{\tau }_1\right)=-a_{23}{a}_{32}{a}_{11}\cos \left({\omega }_2^{\ast }{\tau }_2\right)+a_{23}{a}_{32}{\omega }_2^{\ast }\sin \left({\omega }_2^{\ast }{\tau }_2\right)\end{array}$ , (28)

将(28)两个方程平方然后相加，可得：

${\omega }_2^{\ast }{}^6+E_{41}{\omega }_2^{\ast }{}^4+E_{42}{\omega }_2^{\ast }{}^2+E_{43}+E_{44}\sin \left({\omega }_2^{\ast }{\tau }_1\right)+E_{45}\cos \left({\omega }_2^{\ast }{\tau }_1\right)=0$ , (29)

其中：

$E_{41}=a_{11}^2+a_{22}^2$ , $E_{42}=a_{13}^2{a}_{31}^2-a_{23}^2{a}_{32}^2+a_{11}^2{a}_{22}^2$ , $E_{43}=a_{13}^2{a}_{31}^2{a}_{22}^2-a_{23}^2{a}_{32}^2{a}_{11}^2$ ,

$E_{44}=-2a_{13}{a}_{31}{a}_{11}\left({\omega }_2^{\ast }{}^3+a_{22}^2{\omega }_2^{\ast }\right)$ , $E_{45}=2a_{13}{a}_{31}\left({\omega }_2^{\ast }{}^4+a_{22}^2{\omega }_2^{\ast }{}^2\right)$ .

假设 $\left(H_{41}\right)$ ：(29)存在有限个正根成立，且记为。

根据(28)可得：

,

, (30)

定义，令是特征方程(5)的一个根，且近于，并且满足，。对关于的特征方程(5)两边进行求导，得到：

,

其中：

,

,

,

.

故

, (31)

假设当：成立时，则有。证毕。

第五种情况：。

定理2.5：时，满足

1)：(33)存在有限个正根；

2)：。

系统正平衡点在是局部渐近稳定的，并且当时出现Hopf分支。

证明：当时，将作为参数，假设是方程(5)的根，将其代入并且分离实部和虚部可得：

, (32)

将(32)式中两个方程平方然后相加，可得：

, (33)

其中：

, , ,

,.

假设：(33)存在有限个正根成立，且记为。

根据(32)可得：

,

, (34)

定义，令是特征方程(5)的一个根，且近于，并且满足，。对关于的特征方程(5)两边进行求导，得到：

,

其中：

,

,

,

.

故

, (35)

假设当：成立时，则有。证毕。

3. 数值仿真

前面已经对系统正平衡点的稳定性和Hopf分支做了理论分析，为了进一步探讨系统(1)的复杂的动力学行为，并核实理论结果的有效性和可行性，对系统(1)进行相关动力学模拟试验。对系统的部分参数作如下处理：

, , , , , , , , , ,.

从图1~图4的数值模拟可知，当时滞参数值小于临界值时，系统的正平衡点是局部渐近稳定的，当时滞参数值不断增大时，系统的稳定性会发生变化，而且当时滞参数值大于临界值时系统会出现稳定的周期解。从图5~图8的数值模拟可知，我们同样可以得到相同的结论，当固定一个时滞参数在系统稳定的区间内，当另一个时滞参数值小于临界值时，系统的正平衡点是局部渐近稳定的，而当另一个时滞参数值不断增大时，系统的稳定性会发生变化，当它大于临界值时系统仍会出现稳定的周期解。

(a) (b)

Figure 1. (a) Time series; (b) Phase portrait, the initial value is (3, 3, 10), when,

图1. (a) 当，，初始值为(3, 3, 10)时，系统(1)的时间序列图；(b) 对应系统(1)相图

(a) (b)

Figure 2. (a) Time series; (b) Phase portrait, the initial value is (3, 3, 10), when,

图2. (a) 当，，初始值为(3, 3, 10)时，系统(1)的时间序列图；(b) 对应系统(1)相图

(a) (b)

Figure 3. (a) Time series; (b) Phase portrait, the initial value is (3, 3, 5), when,

图3. (a) 当，，初始值为(3, 3, 5)时，系统(1)的时间序列图；(b) 对应系统(1)相图

(a) (b)

Figure 4. (a) Time series; (b) Phase portrait, the initial value is (3, 3, 5), when,

图4. (a) 当，，初始值为(3, 3, 5)时，系统(1)的时间序列图；(b) 对应系统(1)相图

(a) (b)

Figure 5. (a) Time series; (b) Phase portrait, the initial value is (3, 3, 10), when,

图5. (a) 当，，初始值为(3, 3, 10)时，系统(1)的时间序列图；(b)对应系统(1)相图

(a) (b)

Figure 6. (a) Time series; (b) Phase portrait, the initial value is (3, 3, 10), when,

图6. (a) 当，，初始值为(3, 3, 10)时，系统(1)的时间序列图；(b)对应系统(1)相图

(a) (b)

Figure 7. (a) Time series; (b) Phase portrait, the initial value is (3, 3, 5), when,

图7. (a) 当，，初始值为(3, 3, 5)时，系统(1)的时间序列图；(b) 对应系统(1)相图

(a)(b)

Figure 8. (a) Time series; (b) Phase portrait, the initial value is (3, 3, 5), when,

图8. (a) 当，，初始值为(3, 3, 5)时，系统(1)的时间序列图；(b) 对应系统(1)相图

(a) (b)

Figure 9. (a) Bifurcation diagram with and z, the initial value is (3, 3, 5), when; (b) Bifurcation diagram with and z, the initial value is (3, 3, 10), when

图9. (a) 将作为一个参数，，初始值为(3, 3, 5)时，模型(1)的Hopf分支图；(b) 当作为一个参数，，初始值为(3, 3, 10)时，模型(1)的Hopf分支图

从图9的数值模拟可知，当固定一个时滞参数在系统的稳定区间内，改变另一个时滞参数的值时，系统的稳定性会发生改变，而且当参数值与临界值相等时系统会出现Hopf分支。

4. 结论

本文建立了一类具有双时滞效应的氮–磷–浮游植物动力学新模型，以此为基础对模型的相关动力学性质进行了数学分析与理论推导，解析出模型正平衡点具有局部渐近稳定性和发生Hopf分支的阈值条件，给出一些理论结果，同时模拟仿真出营养及浮游植物密度和Hopf分岔的动态演化过程。本论文的研究结果有助于深入理解时滞效应对海洋生态系统的影响，对预防和控制水体的富营养化具有一定的指导意义。

基金项目

国家自然科学基金面上项目(31570364)。