一类生物入侵模型的解析极限环解

Analytical Limit Cycle Solutions of Abiological Invasion Mode

摘要: 本文利用平方广义谐波函数摄动法研究了一类具有对流–反应–扩散类型的单种群生物入侵模型，求得了该系统极限环解的近似表达式，并得到了极限环幅值与控制参数之间的关系。利用此关系预测了系统在不同参数下的极限环振幅，通过将本方法计算所得结果与数值结果进行比较，验证了所得结果的有效性和可靠性。该结果表明平方广义谐波函数摄动法亦可推广至一类无限维动力系统Hopf分岔的解析定量分析，为无限维动力系统的研究相提供了新的思路和参考方法。

Abstract: In this paper, a single species biological invasion model with convection diffusion reaction type is studied by using the square generalized harmonic function perturbation method. The approximate expression of the limit cycle solution of the system is obtained, and the relationship between the limit cycle amplitudes and the control parameters is also obtained. Based on this relationship, the limit cycle amplitudes of the system under different parameters are predicted. By comparing the results calculated by this method with the numerical results, the effectiveness and reliability of the results are verified. The medium and long-term prediction of the population density development of a class of species and the effective control of Hopf bifurcation behavior are realized, which provides a new idea and reference method for the analytical and quantitative study of Hopf bifurcation of a class of infinite dimensional dynamic systems.

1. 引言

生物数学作为数学与生命科学的交叉学科，对于揭示纷繁复杂的自然现象背后所隐藏的奥妙起着重要意义。目前所发现史料中有关种群动态学模型的最早记载可以追溯到我国宋朝时期。1088年，沈括曾在《梦溪笔谈》中用“纳甲法”推出“胎育之理”[1]，其实质就是一种生物学模型。之后，生物数学不断发展完善，并于1974年被列为一门独立学科[2]。

种群动力学是生物数学一个非常重要的分支领域[3]，主要研究种群与环境、种群与种群之间的相互影响与相互作用，并通过建立动力学模型来分析一些生态现象，从而达到预测、调节和控制物种的发展过程和发展趋势的目的。因此，种群动力学研究引起了国内外众多学者的关注。2010年，Spiegler等通过余维一分岔分析分别研究了输入和树突时间常数对神经群极限环振荡特性的影响规律[4]；2010年，Basu等描述了一类神经元的非线性动力学现象，通过分岔分析研究了当模型接近Hopf分岔和鞍结分岔时，简单和复杂模型动力学行为的相似性[5]；2011年，Touboul等通过神经群的突触连接总数目与外界输入的余维二分岔分析研究了二者的相互作用对神经群动力学特征性的影响[6]；2013年，唐文艳等对具有脉冲扩散效应的Gomportz种群动力学模型研究，得到系统存在周期解的结论，该结论表明通过种群扩散可以使得生物物种持续生存，为现实的生物资源管理提供了可靠的策略依据[7]；2016年，任善静等研究了具脉冲与基因单点突变效应的种群动力学模型，利用离散动力系统频闪映射理论，得到了控制脉冲出生、脉冲收获及脉冲剔除在不同时刻的基因单点突变单种群动力系统的充分条件[8]；2018年，焦建军等利用常微分方程差分分析方法，研究了毒素脉冲输入与脉冲出生切换阶段结构单种群动力学模型和污染环境下具瞬时与非瞬时脉冲出生单种群动力学模型，获得了上述系统种群灭绝和持久生存的控制条件结果，为污染环境中的生物资源管理提供了可靠的管理策略[9] [10]。2019年，Duan等研究了一类具有中性泛函微分方程形式的年龄结构种群模型。通过对特征方程的分析，讨论了正平衡的稳定性，选择成熟时延作为分岔参数，得到了局部Hopf分岔结果[11]。2019年，曹建智等研究了一类具有时滞的云杉蚜虫种群阶段结构模型的动力学行为，讨论了模型正平衡点的存在唯一性并分析了该平衡点的局部稳定性和出现Hopf分岔的充分条件[12]。可见，对种群动力学模型进行非线性分析是当前生物数学领域的热点问题之一。由于大多数种群动力学模型都是非线性系统，不可避免的存在复杂的分岔现象。其中，Hopf分岔是非线性动力系统典型的全局动力学行为，主要研究从奇点邻域内产生极限环的问题。当系统在某奇点的线性化矩阵有一对纯虚根时，只要参数变化引起奇点稳定性的改变，并超过某一临界值时，系统将从奇点分出一族周期轨道而产生分岔现象。因此，Hopf分岔问题对研究种群动力学模型长时间行为有重要意义。

Petrovskii S等建立的单物种种群密度的动力学模型考虑了种群迁移和Allee效应的影响[13]，可用如下方程描述：

$\frac{\partial U\left(X,T\right)}{\partial T}+\left(A_0+2A_{1}U\right)\frac{\partial U}{\partial X}=D\frac{{\partial }^{2}U}{\partial X^2}+\alpha U\left(U-U_0\right)\left(K-U\right)$ (1)

其中$U$ 表示种群密度，$K$ 表示种群承载能力，$U_0$ 为Allee效应的度量，$\alpha$ 为系数，$A_0$ 表示在风速和水流影响下的对流速度，$A_1$ 表示迁徙速度。将上式进行如下的无量纲化变换：

$u=\frac{U}{K},t=T\alpha K^2,x=X\sqrt{\frac{\alpha K^2}{D}}$ (2)

可化为：

$u_t+\left(a_0+a_{1}u\right)u_x=u_{xx}-\beta u+\left(1+\beta \right)u^2-u^3$ (3)

其中$\beta =U_0/K$ ，$a_0=A_0{K}^{-1}{\left(\alpha D\right)}^{-1/2}$ ，$a_1=2A_1{\left(\alpha D\right)}^{-1/2}$ 为无量纲参数，下标x和t表示无量纲的空间和时间。Sheratt等在文献[14] [15]中，利用数值延拓法研究了该类系统的Hopf分岔和周期行波解。文献[16]利用文献[17] [18]中介绍的方法，计算了系统(3)平衡点的奇点量，得到了该系统存在Hopf分岔的条件。文献[13]通过一个特殊的变量变换，求得了系统(3)在特定参数条件下的一个精确解。本文则利用平方广义谐波函数摄动法[19] [20]，对该生物入侵模型进行了定量研究，求得了该系统极限环解的近似表达式，并得到了极限环初值与控制参数之间的关系。并利用该关系，对不同参数下系统的极限环初值进行预测。通过将预测结果与数值结果进行比较，验证该方法的有效性和可靠性。

2. 平方广义谐波函数摄动法

本节主要介绍利用平方广义谐波函数摄动法求解无限维动力系统解析近似解的主要步骤。假设一类具有对流–反应–扩散(advection-diffusion-reaction)类型的无限维动力系统，可化为如下常微分方程：

$u^{″}\left(\xi \right)+p\left(u\right)=\delta f\left(\nu ,u,u_{\xi }\right)$ (4)

其中$\delta >0$ ，$\nu$ 为控制参数。一般而言，当$\delta \ne 0$ 时，不易求出上式的精确值，但当$\delta$ 不是十分大时，采用平方广义谐波函数摄动法，可求解该系统的近似解。下面对系统(4)进行如下非线性变量变换$\phi \left(\xi \right)$ ：

$\frac{\text{d}\phi }{\text{d}\xi }=\Phi \left(\phi \right),\Phi \left(\phi +\pi \right)=\Phi \left(\phi \right)$ (5)

并设系统(4)的解具有如下形式：

$u\left(\phi \left(\xi \right)\right)=b{\cos }^2\phi +d$ (6)

其中$d\in R$ ，$b$ 和$\phi$ 均可展成如下形式的级数：

$b=b_0+\delta b_1+\cdots$ (7)

$\frac{\text{d}\phi }{\text{d}\xi }=\Phi \left(\phi \right)={\Phi }_0+\delta {\Phi }_1+\cdots$ (8)

则解(6)可以写成：

$u=u_0+\epsilon u_1+\cdots$ (9)

其中：

$\{\begin{array}{l}{u}_0=b_0{\cos }^2\phi +d\\ u_n=b_n{\cos }^2\phi \left(n=1,2,\cdots \right)\end{array}$ (10)

将式(8)和式(9)代入式(4)，并比较方法两边同阶$\epsilon$ 的系数可得：

${\delta }^0:{\Phi }_0\frac{\text{d}}{\text{d}\phi }\left({\Phi }_0{u^{\prime }}_0\right)+p\left(u_0\right)=0$ (11)

$\begin{array}{l}{\delta }^1:{\Phi }_0\frac{\text{d}}{\text{d}\phi }\left({\Phi }_1{u^{\prime }}_0\right)+{\Phi }_1\frac{\text{d}}{\text{d}\phi }\left({\Phi }_0{u^{\prime }}_0\right)+{\Phi }_0\frac{\text{d}}{\text{d}\phi }\left({\Phi }_0{u^{\prime }}_1\right)\\ +p^{\prime }\left(u_0\right)u_1=f\left(\nu ,u_0,{\Phi }_0{u^{\prime }}_0\right)\end{array}$ (12)

$\begin{array}{l}{\delta }^2:{\Phi }_0\frac{\text{d}}{\text{d}\phi }\left({\Phi }_2{u^{\prime }}_0\right)+{\Phi }_2\frac{\text{d}}{\text{d}\phi }\left({\Phi }_0{u^{\prime }}_0\right)+{\Phi }_0\frac{\text{d}}{\text{d}\phi }\left({\Phi }_0{u^{\prime }}_2\right)+{\Phi }_1\frac{\text{d}}{\text{d}\phi }\left({\Phi }_0{u^{\prime }}_1\right)\\ +{\Phi }_0\frac{\text{d}}{\text{d}\phi }\left({\Phi }_1{u^{\prime }}_1\right)+{\Phi }_1\frac{\text{d}}{\text{d}\phi }\left({\Phi }_1{u^{\prime }}_0\right)+g^{\prime }\left(u_0\right)u_2+\frac{1}{2}{g}^{″}\left(u_0\right)u_1^2\\ ={f^{\prime }}_u\left(\nu ,u_0,{\Phi }_0{u^{\prime }}_0\right)u_1+{f^{\prime }}_{u^{\prime }}\left(\nu ,u_0,{\Phi }_0{u^{\prime }}_0\right)\left({\Phi }_0{u^{\prime }}_1+{\Phi }_1{u^{\prime }}_0\right)\end{array}$ (13)

$⋮$

$⋮$

其中$u^{\prime }=\text{d}u/\text{d}\phi$ ，${f^{\prime }}_u=\partial f/\partial u$ ，${f^{\prime }}_{u^{\prime }}=\partial f/\partial u^{\prime }$ 。

第一步，求解零阶摄动方程。首先，构造一类新的广义谐波函数，即平方广义谐波函数，来描述非线性振子的解。考虑如下形式的非线性振子：

$u^{″}+\text{g}\left(u\right)=0$ (14)

其初值为：

$u\left(0\right)=e_1,u^{\prime }\left(0\right)=e_2$ (15)

其中$e_1,e_2\in R$ 。引入如下形式的非线性时间变换$\phi \left(t\right)$ ：

$\frac{\text{d}\phi }{\text{d}t}=\Phi \left(\phi \right),\Phi \left(\phi +\pi \right)=\Phi \left(\phi \right)$ (16)

将上述变换代入方程(14)，可将该方程变换为如下形式：

$\Phi \frac{\text{d}}{\text{d}\phi }\left(\Phi u^{\prime }\right)+g\left(u\right)=0$ (17)

其中求导运算表示对$\phi$ 求导，$g\left(u\right)$ 为非线性函数。本文构造如下形式的广义谐波函数：

$u\left(\xi \right)=a{\cos }^2\phi \left(\xi \right)+b$

其中$a,b\in R$ 。将(14)式两边同时式两边同时乘以$u^{\prime }\left(\phi \right)=-2a\cos \phi \sin \phi$ 并关于变量$\phi$ 积分，可以得到：

$\frac{1}{2}{\left(\Phi {u^{\prime }}_{}\right)}^2=-{\displaystyle {\int }_{u(0)}^{u(\phi )}g\left(u\right)\text{d}u}$ (18)

令上式中$\phi =\pi /2$ ，可得：

${\displaystyle {\int }_{a+b}^{b}g\left(u\right)\text{d}u}=0$ (19)

其中，$a+b$ 等于初值$e_1$ 。令$e_2=0$ ，则$a$ 和$b$ 可由(20)式求得。由(19)式可得：

$\Phi \left(\phi \right)=\sqrt{\frac{2\left[V\left(a+b\right)-V\left(a{\cos }^2\phi +b\right)\right]}{{\left({u^{\prime }}_{}\left(\phi \right)\right)}^2}},\phi \ne \frac{n}{2}\pi \left(n=0,\pm 1,\pm 2,\cdots \right)$ (20)

其中$V\left(u\right)={\displaystyle {\int }_0^{x}g\left(u\right)\text{d}u}$ ，$\Phi \left(0\right)=\sqrt{\left|g\left(a+b\right)/2a\right|}$ ，$\Phi \left(\pi /2\right)=\sqrt{\left|g\left(b\right)/2a\right|}$ 。

基于上述步骤，即可求出方程(11)的周期解，即零阶摄动解：

$u_0=b_0{\cos }^2\phi +d$ (21)

其中：

${\displaystyle {\int }_{b_0+d}^{d}p\left(u\right)\text{d}u}=0$ (22)

$\{\begin{array}{l}{\Phi }_0\left(\phi \right)=\sqrt{\frac{2{{\displaystyle \int }}_{b_0{\cos }^2\phi +d}^{b_0+d}p\left(u\right)\text{d}u}{{\left({u^{\prime }}_0\left(\phi \right)\right)}^2}},\phi \ne \frac{n}{2}\pi \left(n=0,\pm 1,\pm 2,\cdots \right)\\ \Phi (0)=\sqrt{\left|p\left(b_0+d\right)/2a\right|},\Phi \left(\pi /2\right)=\sqrt{\left|p\left(d\right)/2a\right|}\end{array}$ (23)

第二步，求一阶摄动解。将式(12)两边同时乘以${u^{\prime }}_0$ 并关于$\phi$ 积分，可得：

${{\Phi }_0{\Phi }_1{u^{\prime }}_0^2|}_{{\phi }_0}^{\phi }={\displaystyle {\int }_{{\phi }_0}^{\phi }f\left(\mu ,x_0,{\Phi }_0{u^{\prime }}_0\right)}{u^{\prime }}_0\text{d}\phi -{\frac{b_1}{2b_0}{\left({\Phi }_0{u^{\prime }}_0\right)}^2|}_{{\phi }_0}^{\phi }-{u_{1}p\left(u_0\right)|}_{{\phi }_0}^{\phi }$ (24)

令式(25)中的${\phi }_0=0$ ，$\phi =\pi$ ，可得：

${\displaystyle {\int }_0^{\pi }f\left(\mu ,u_0,{\Phi }_0{u^{\prime }}_0\right)}{u^{\prime }}_0\text{d}\phi =0$ (25)

令式(25)中的${\phi }_0=0$ ，$\phi =\pi /2$ ，可得：

${\displaystyle {\int }_0^{\frac{\pi }{2}}f\left(\mu ,u_0,{\Phi }_0{u^{\prime }}_0\right)u_0\text{d}\phi }+b_{1}p\left(b_0+d\right)=0$ (26)

再令式(25)中的${\phi }_0=0$ ，则有：

${\displaystyle {\int }_0^{\phi }f\left(\mu ,u_0,{\Phi }_0{u^{\prime }}_0\right)}{u^{\prime }}_0\text{d}\phi -{\Phi }_0{\Phi }_1{\left({u^{\prime }}_0\right)}^2-\frac{b_1}{2b_0^{}}{\left({\Phi }_0{u^{\prime }}_0\right)}^2-u_{1}p\left(u_0\right)+b_{1}p\left(b_0+d\right)=0$ (27)

在极限环的求解中，$b_0$ ，$b_1$ ，$d$ ，${\Phi }_0$ 和${\Phi }_1$ 为未知数。通过式(23)至(28)，可求出上述未知数。由于${\Phi }_0\left(\phi \right)$ 是周期函数，为了便于后续的计算，将${\Phi }_0\left(\phi \right)$ 展成如下傅里叶级数：

${\Phi }_0\left(\phi \right)=\sqrt{\frac{2{{\displaystyle \int }}_{b_0{\cos }^2\phi +d}^{b_0+d}p\left(u\right)\text{d}u}{{\left({u^{\prime }}_0\left(\phi \right)\right)}^2}}={\displaystyle \sum _{n=0}^M\left(p_{2i}\cos 2n\phi +q_{2i}\sin 2n\phi \right)}$ (28)

其中$M\in N^{*}$ 。令：

${\Phi }_0\left(\phi \right)=\sqrt{\frac{2{{\displaystyle \int }}_{b_0{\cos }^2\phi +d}^{b_0+d}p\left(u\right)\text{d}u}{{\left({u^{\prime }}_0\left(\phi \right)\right)}^2}}={\displaystyle \sum _{n=0}^M\left(p_{2i}\cos 2n\phi +q_{2i}\sin 2n\phi \right)}$ (29)

其中$u_0$ 为零阶摄动解。将式(22)，(24)代入式(30)后，得到一个含有未知参数$b_0$ 和$d$ 的表达式，这意味着式(30)可以简记为$I\left(\phi ;b_0,d\right)$ 。因此，式(26)可以改写为：

${I\left(\phi ;b_0,d\right)|}_0^{\pi }=0$ (30)

将上述方程与式(23)联立，即可求出$b_0$ 和$d$ 。将式(30)代入式(27)，可得：

$b_1=-\frac{{I\left(\phi \right)|}_0^{\pi /2}}{p\left(b_0+d\right)}$ (31)

通过将式(30)代入式(28)，可得：

${\Phi }_1=\frac{I\left(\phi \right)-\frac{b_1}{2b_0^{}}{\left({\Phi }_0{u^{\prime }}_0\right)}^2-b_{1}p\left(u_0\right)\cos \left(\phi \right)+b_{1}p\left(b_0+d\right)}{{\Phi }_0{u^{\prime }}_0^2}$ (32)

最后，方程(4)的一阶近似极限环解可表达为：

$u\left(\xi \right)=b_0{\cos }^2\phi \left(\xi \right)+b_1{\cos }^2\phi \left(\xi \right)+d+O\left({\delta }^2\right)$ (33)

$u^{\prime }\left(\xi \right)=-2b_0\left[{\Phi }_0+\delta {\Phi }_1\right]\cos \phi \left(\xi \right)\sin \phi \left(\xi \right)-2b_1\left[{\Phi }_0+\delta {\Phi }_1\right]\cos \phi \left(\xi \right))\sin \phi \left(\xi \right)+O\left({\delta }^2\right)$ (34)

继续按照类似的步骤，可求解二阶近似解$u_2$ 和${\Phi }_2$ ，但是随着求解阶数的增加，计算会越来越复杂，更重要的是，在$\epsilon$ 取1或2时，利用该方法求得的一阶近似解已经具有较理想的精度。下一节将通过一个实例来表明上述方法的有效性。

3. 一个单物种种群密度动力学模型的极限环

首先对系统(3)进行如下的行波变换：

$u\left(x,t\right)=u\left(\xi \right),\xi \to x-ct$ (36)

其中$c$ 为行波波速。基于该变换，系统(1.3)可化为如下常微分系统：

$u^{″}\left(\xi \right)-\beta u\left(\xi \right)+\left(1+\beta \right)u^2\left(\xi \right)-u^3\left(\xi \right)=\left(a_0-c+a_{1}u\left(\xi \right)\right)u^{\prime }\left(\xi \right)$ (35)

将上式右端视为扰动项，并将上式改写为

$u^{″}-\beta u+\left(1+\beta \right)u^2-u^3=\delta \left({\nu }_c+a_{1}u\right)u^{\prime }$ (36)

其中求导运算表示对$\xi$ 求导，参数$\delta >0$ ，${\nu }_c=a_0-c$ 为控制参数。该方程有三个平衡点，其坐标为$\left(u,u^{\prime }\right)=\left(0,0\right),\left(1,0\right),\left(\beta ,0\right)$ 。由于系统的应用背景，本文仅讨论$u>0$ 时的情况，因此仅讨论平衡点$\left(1,0\right)$ 和$\left(\beta ,0\right)$ 。通过分析可知，当参数$\delta =0$ 时，若$\beta >1$ ，则点$\left(1,0\right)$ 为中心，$\left(\beta ,0\right)$ 为鞍点；若$0<\beta <1$ 时，点$\left(\beta ,0\right)$ 为中心，$\left(1,0\right)$ 为鞍点。本节将要讨论的是当$\delta \ne 0$ 时，系统(3.3)存在极限环的条件，以及极限环振幅与参数之间的关系。将

$p\left(u\right)=-\beta u+\left(1+\beta \right)u^2-u^3$ (37)

$f\left(\nu ,u,u_{\xi }\right)=\left({\nu }_c+a_{1}u\right)u^{\prime }$ (38)

代入式(23)和式(31)，并令式(29)中的正整数$M=10$ ，即可求出$b_0$ 的表达式为：

$b_0=-\frac{4\left(2d{a}_1{p}_0-d{a}_1{p}_4)+4{\nu }_c(2p_0-p_4\right)}{a_1\left(4p_0+p_2-2p_4-p_6\right)}$ (39)

其中$p_i$ 为${\Phi }_0$ 的傅里叶展开系数，${\Phi }_0$ 由式(24)给出，$a_1$ ，${\nu }_c$ 为系统参数。上式与式(23)联立，即可求出$b_0$ 和$d$ 。将$b_0$ ，$d$ 和${\Phi }_0$ 代入式(32)和式(33)可得：

$b_1=0$ (40)

$\{\begin{array}{l}{\Phi }_1=\frac{I\left(\phi \right)}{{\Phi }_0{u^{\prime }}_0^2},\phi \ne \frac{n}{2}\pi \left(n=0,\pm 1,\pm 2,\cdots \right)\\ {\Phi }_1\left(0\right)={\Phi }_1\left(\pi \right)={I\left(\phi \right)|}_{\phi =0or\pi }^{},{\Phi }_1\left(\pi /2\right)={I\left(\phi \right)|}_{\phi =\pi /2}^{}\end{array}$ (41)

为了便于运算和表达，将${\Phi }_1\left(\phi \right)$ 也展成如下傅里级数：

${\Phi }_1\left(\phi \right)={\displaystyle \sum _{n=0}^{\tilde{M}}\left(r_{2i}\cos 2n\phi +w_{2i}\sin 2n\phi \right)}$(42)

其中$\tilde{M}$ 为正整数。将求得的$b_0,d,{\Phi }_0,b_1,{\Phi }_1$ 代入式(34)和(35)即可求得系统(3)的极限环解。由于$b_1=0$ ，因此极限环的初始值A可由下式近似表示：

$A=b_0+d=-\frac{4\left(2d{a}_1{p}_0-d{a}_1{p}_4)+4{\nu }_c(2p_0-p_4\right)}{a_1\left(4p_0+p_2-2p_4-p_6\right)}+d$ (43)

由此也可得到系统(3)存在极限环的一个必要条件，那就是将系统参数代入式(23)，(31)，(29)和式(41)后，由式(45)给出的A为非零实数。

Figure 1. The phase portrait of system (3)

图1. 系统(3)的大致相图

下面根据一组具体的参数，来说明该方法的有效性。首先，讨论(1，0)附近的极限环，令$\beta =1.5$ ，${\nu }_c=1.01$ ，$a_1=-1$ 。此时，系统(1.3)的大致相图如图1所示。由图可见，平衡点$\left(1.5,0\right)$ 的左侧存在一条同宿轨，在该同宿轨内存在极限环。将上述参数代入(23)、(24)、(41)、(42)、(43)，即可求得极限环的解为：

$x=0.3866{\cos }^2\phi +0.8204+O\left({\epsilon }^2\right)$ (44)

$\frac{\text{d}\phi }{\text{d}t}={\Phi }_0\left(\phi \right)+\delta {\Phi }_1\left(\phi \right)+O\left({\delta }^2\right)$ (45)

(46)

${\Phi }_1=-0.0340\sin \left(2\phi \right)-0.0003\sin \left(4\phi \right)-0.0001\sin \left(6\phi \right)-0.00002\sin \left(8\phi \right)+HH$ (47)

其中HH为高阶谐波项。为了书面表达的整洁，这里只给出了傅里叶展开的前四项。$\delta =1$ 时，解(46)的相图如图2(d)所示。从图2(d)中的比较不难看出，所求得的解与数值解吻合的较好。该图中，$u$ 表示无量纲的种群密度；$u^{\prime }$ 表示无量纲的种群密度演化速度；$\xi$ 表示无量纲的时空变量。

下面继续令$\beta =1.5$ ，$a_1=-1$ ，$\delta =1$ ，并通过式(45)计算了${\nu }_c$ 在不同值时，系统(3)极限环的初始值，如表1所示。其中$A$ 表示本文算得的极限环初值，$A_N$ 表示利用数值方法算得的极限环初值。从该表可以看出，式(45)预测的极限环初值有着较好的精度，从而进一步表明了该方法的有效性和可靠性。

Figure 2. The limit cycle for system (3). (a) The evolution of limit cycle on $u-u^{\prime }$ plane; (b) (c) The evolution of limit cycle on $\xi -u$ plane; (d) The phase portrait of limit cycle solution, where — denotes the solution obtained by Runge-Kutta method; $\cdot \cdot \cdot \cdot \cdot$ denotes the present solution (46)

图2. 系统(3)的极限环。(a) $u-u^{\prime }$ 平面上的极限环演化图；(b) (c) $\xi -u$ 上的极限环演化图；(d) 极限环解的相图，其中—表示Runge-Kutta所得之解；$\cdot \cdot \cdot \cdot \cdot$ 表示本文所得之解(46)

Table 1. The comparisons of amplitude of limit cycles with different control parameters for system (3)

表1. 系统(3)在不同控制参数下的极限环初始值比较

${\nu }_c$	$A$	$A_N$	相对误差
1.00238	1.10000	1.10000	0%
1.00529	1.15000	1.15030	0.02608%
1.00934	1.20000	1.20070	0.05683%
1.01456	1.25000	1.25130	0.10400%
1.02098	1.30000	1.30240	0.18461%
1.02857	1.35000	1.35400	0.29629%
1.03715	1.40000	1.40660	0.47143%
1.04598	1.45000	1.46190	0.82069%

本节讨论了系统(3)在平衡点(1，0)附近的极限环，预测了在不同控制参数下的极限环幅值，通过与数值结果进行比较，验证了本文所得之结果具有较好的精度。同理，平衡点$\left(\beta ,0\right)$ 附近的极限环可用相同的方法讨论。

4. 结论

本文主要研究了平方广义谐波函数摄动法在一类无限维动力系统极限环求解中的应用。利用平方广义谐波函数摄动法对一类生物入侵模型展开了定量研究，不仅得到了该系统极限环的解析近似表达式，还得到了极限环初值与控制参数之间的关系式。利用该关系式对该模型不同控制参数下的极限环幅值进行了预测，通过将预测结果与数值结果进行比较，证明了该方法所得结果有着较好的精度。同时也表明，平方广义谐波函数摄动法在无限维动力系统全局分岔的研究中依然有效。对该方法继续展开深入研究，有望形成针对一类无限维动力系统Hopf分岔进行解析定量研究的新方法。为无限维动力系统的定量研究提供新的思路和参考方法。

致 谢

感谢湖南省大学生创新创业训练项目对本文的支持。

基金项目

湖南省大学生创新创业训练项目(编号：S202210555196, X202310555175)资助的课题。

NOTES

^*通讯作者。