近年来,随着人类社会高速发展,生态问题日益突出,环境污染、物种灭绝引起了人们的高度重视,进而催生了种群生态学的快速发展.其中运用数学模型分析种群变化规律是生态学中重要的一个方法.人们最关心的问题是如何维持种群数量的合理、保持生态系统的稳定,而在数学上,这类问题就是关于竞争模型解的稳定性问题.本文对三种群的Lotka-Volterra竞争模型进行了分析.

Lotka-Volterra捕食者-食饵模型最初是由Lotka于1910年在自催化反应理论^[1-4]中提出. 1926年数学家和物理学家Volterra提出了同样的方程组^[5-6],并用它解释一战期间亚得里亚海掠食性鱼类数量与渔捞量之间的关系^[7].由于种群竞争模型分别由Lotka和Volterra两人提出,因此后来称之为Lotka-Volterra模型. Lotka-Volterra模型被用来解释捕食者和被捕食者的自然种群的动态,例如哈德逊湾公司的山猫和雪兔的数量^[8]、罗亚岛国家公园里的鹿群和狼群数量^[9].该模型为研究种群竞争关系提供了理论基础,并对现代生态学理论的发展产生了重要影响.基于该模型的基本概念和理论框架,人们进行了一些改进和实践研究^[10-13].徐瑞等人^[14]研究了一类具有时滞和基于比率的三种群Lotka-Volterra模型.证明了该系统在适当条件下的一致持久性,并通过构造李亚普诺夫函数的方法,得到了使该系统正平衡点全局渐近稳定的充分条件.陆忠华等人^[15]研究了既有捕食关系又有竟争关系的三种群Lotka-Volterra模型,从周期系数的观点出发,得到了全局渐近稳定周期解的唯一存在条件.伏升茂等人^[16]研究了具有自扩散和互扩散的三种群Lotka-Volterra模型,从自扩散和交叉扩散的角度讨论了一致有界解的全局存在性,证明了正平衡点的全局渐近稳定性.在竞争模型中,物种的持久性和灭绝与物种的初始数量有关,而分形理论中的Julia集与初始点的轨迹相关,这为我们研究Lotka-Volterra模型提供了一种新思路.据我们所知,目前为止运用分形理论讨论Lotka-Volterra模型尚未少数,因此本文从分形的角度对竞争模型进行了研究.

分形是自然界中常见的现象,它解释了许多以前人们无法解释的非线性现象,比如湍流问题、癌细胞的生长等等,因此成为近几年来热门的前沿课题.其中Julia集是分形理论中十分重要的一支,相应的研究成果在多个领域得到了应用.如Michael Levin指出广义复映射的填充Julia集可以有效地模拟实际生物体或生物器官的形态特征和生长规律^[17]. Sun等人利用M-J集作为加密字典,利用Julia集的初值敏感性大大增加了加密的稳定性和密钥储存空间^[18]. Shudo等人将Julia集应用到量子隧道效应的复半经典理论中,利用分形集合的自相似性从宏观状态研究量子力学中的微观特性^[19].

目前为止运用Julia集研究种群竞争模型主要集中在二维生态系统^[20].但是生物界中大量存在着三维甚至更高维数的种群竞争,如狮子与斑马、水牛,雕与绵羊、鹿.因此本文将Julia集引入三种群捕食-食饵模型进行讨论.在三种群捕食-食饵模型中, $ z_1 $种群以$ x_1 $、$ y_1 $种群为食,食饵种群间又存在资源竞争.

本文首先介绍了种群竞争模型的Julia集,讨论了初始种群数对模型的影响.然后采用辅助参考反馈控制方法对模型的Julia集进行了控制,实现了模型由不稳定向稳定的转换,同时讨论了当一个种群初始数量变化时,另外两个种群相应的变化.最后通过引入耦合控制项,实现两个参数不同的系统之间的同步. MATLAB仿真证明了控制的有效性.

设$ E $是$ {\Bbb R}^{n} $的子集(通常是$ R^{n} $本身), $ f:E\rightarrow E $是一连续映射, $ f^{k} $表示$ f $的$ k $次迭代,则$ E $的子集$ F $被称为$ f $的吸引子,如果$ f $是一个闭集而且是不变的(即$ f(F) = F $),对包含$ F $的一个开集$ V $中的所有$ x $,随着$ k $趋近于无穷大, $ f^{k}(x) $与$ F $之间的距离趋近于零,则称集$ V $为$ F $的吸引域.也就是说,吸引子是一个集,它使附近的所有轨迹都收敛到它.同样的,如果$ F $附近的所有点(但不在$ F $中)经过迭代都远离$ F $,则称闭不变子集$ F $为斥子.

假设$ f:C\rightarrow C $是复平面上的映射.如果$ f(\omega) = \omega $,则称$ \omega $为$ f $的不动点.如果存在大于$ 1 $的整数$ p $使得$ f^{p}(\omega) = \omega $,则称$ \omega $为$ f $的周期点.满足$ f^{p}(\omega) = \omega $的最小的$ p $称为$ \omega $的最小周期,且称$ \omega, f(\omega), \cdot\cdot\cdot, f^{p}(\omega) $为周期$ p $的轨道.设$ \omega $为周期点, $ p $为周期, $ (f^{p})'(\omega) = \lambda $,则当(1) $ \lambda = 0 $, $ \omega $为超吸引的; (2) $ 0 < \mid\lambda\mid < 1 $, $ \omega $为吸引的; (3) $ \mid\lambda\mid = 1 $, $ \omega $为中性的; (4) $ \mid\lambda\mid > 1 $, $ \omega $为排斥的.

定义2.1  设$ f $是复平面$ C $上的复映射,则$ f $的Julia集定义为$ f $的斥性周期点的闭包.

用$ J(f) $表示$ f $的Julia集,具有以下性质

(ⅰ) $ J(f) $是非空有界的;

(ⅱ) $ J(f) $是全不变的,即$ J(f) = f(J(f)) = f^{-1}(J(f)) $;

(ⅲ)对任何正整数$ p $, $ J(f) = J(f^{p}) $成立;

(ⅳ) $ J(f) $是任意吸引不动点$ f $的吸引域的边界,即$ J(f) = \partial A(\omega) $,其中$ \partial A(\omega) $是吸引不动点$ \omega $的吸引域$ A(\omega) $的边界.

本文研究了以下三种物种竞争模型

(3.1) $ \begin{eqnarray} \left\{\begin{array}{l} \frac{{\rm d}x}{{\rm d}t} = x(1-x)-a_{1}xy-b_1xz, \\ \frac{{\rm d}y}{{\rm d}t} = y(1-y)-a_{2}xy-b_2yz, \\ \frac{{\rm d}z}{{\rm d}t} = c_1xz+c_2yz-d_1z-d_2z^2. \end{array} \right. \end{eqnarray} $

根据欧拉离散方法,取步长为1,可得下列离散方程

(3.2) $ \begin{eqnarray} \left\{\begin{array}{l} x_{n+1} = 2x_{n}-x_{n}^2-a_1 x_{n}y_{n}-b_1 x_{n}z_{n}, \\ y_{n+1} = 2y_{n}-y_{n}^2-a_2 x_{n}y_{n}-b_2 y_{n}z_{n}, \\ z_{n+1} = (1-d_1)z_{n}-d_2z_{n}^2+c_1 x_{n}z_{n}+c_2 y_{n}z_{n}. \end{array} \right. \end{eqnarray} $

其中, $ x, y $表示两种食饵种群, $ z $表示捕食种群. $ x_{n} $, $ y_{n} $, $ z_{n} $($ x_{n} > 0, y_{n} > 0, z_{n} > 0 $)为当前的物种密度, $ d_{1} $表示捕食者出生率, $ a_{1}, a_{2}, b_{1}, b_{2}, c_{1}, c_{2}, d_{2} $ (均为正)描述了物种间的相互作用.取$ a_1 = 0.97, \; a_2 = 0.81, \; b_1 = 0.39, \; b_2 = 0.85, \; c_1 = 0.23, \; c_2 = 0.13, \; d_1 = 0.17, d_2 = 0.01 $,则相应的Julia集如图 1所示.为了防止种群数量无限发展,我们需要将种群的初始数量限制在Julia集内,保证种群数量稳定.

设该系统的不动点$ (x^*, y^*, z^*) $,即

(3.3) $ \begin{eqnarray} \left \{\begin{array}{l} 2x^*-(x^*)^2-a_1 x^*y^*-b_1 x^*z^* = x^*, \\ 2y^*-(y^*)^2-a_2 x^*y^*-b_2 y^*z^* = y^*, \\ (1-d_1)z^*-d_2(z^*)^{2}+c_1 x^*z^*+c_2 y^*z^* = z^*. \end{array} \right. \end{eqnarray} $

其解分为以下几种情况

(1) $ (x^*, y^*, z^*) = (0, 0, 0). $

(2) $ x^* = 0, y^*\neq0, z^*\neq0 $,代入(3.3)式得到不动点$ (x^*, y^*, z^*) = (0, \frac{d_{2}+b_{2}d_{1}}{d_{2}+b_{2}c_{2}}, \frac{c_{2}-d_{1}}{d_{2}+b_{2}c_{2}}). $同理, $ x^*\neq0, y^* = 0, z^*\neq0 $时,不动点为$ (x^*, y^*, z^*) = (\frac{d_{2}+b_{1}d_{1}}{d_{2}+b_{2}c_{1}}, 0, \frac{c_{1}-d_{1}}{d_{2}+b_{1}c_{1}}). $$ x^*\neq0, y^*\neq0, $$ z^* = 0 $时,不动点为$ (x^*, y^*, z^*) = (\frac{1-a_{1}}{1-a_{1}a_{2}}, \frac{1-a_{2}}{1-a_{1}a_{2}}, 0). $

(3) $ x^* = 0, y^* = 0, z^*\neq0 $,代入(3.3)式得到不动点$ (x^*, y^*, z^*) = (0, 0, \frac{-d_{1}}{d_{2}}) $同理,当$ x^* = 0, y^*\neq0, z^* = 0 $和$ x^*\neq0, y^* = 0, z^* = 0 $时,不动点分别为$ (x^*, y^*, z^*) = (0, 1, 0). $和$ (x^*, y^*, z^*) = (1, 0, 0). $

(4) $ x^*\neq0, y^*\neq0, z^*\neq0 $,考虑该离散系统的共存不动点$ (x^*, y^*, z^*) = (\frac{D_x}{D}, \; \frac{D_y}{D}, \; \frac{D_z}{D}) $,其中

$ \begin{eqnarray*} &&D = -d_2-b_2c_2+a_1a_2d_2+a_1b_2c_1+a_2b_1c_2-b_1c_1, \\ &&D_x = -d_2-b_2c_2+a_1d_2+a_1d_1b_2+b_1c_2-b_1d_1, \\ &&D_y = -d_2-d_1b_2+a_2d_2+b_2c_1+b_1a_2d_1-b_1c_1, \\ &&D_z = d_1-c_2-a_1a_2d_1+a_1c_1+a_2c_2-c_1. \end{eqnarray*} $

系统(3.2)在$ (x^*, y^*, z^*) $处的Jacobi矩阵为

$ \begin{eqnarray*} J = \left[ \begin{array}{ccc} 1-x^* &\quad -a_1x^* \quad& -b_1x^*\\ -a_2y^* & 1-y^* & -b_2y^*\\ c_1z^* & c_2z^* & -d_2z^*+1\\ \end{array} \right]. \end{eqnarray*} $

令

$ \begin{eqnarray*} \mid\lambda E-J\mid = \left| \begin{array}{ccc} \lambda-1+x^* & a_1x^* & b_1x^*\\ a_2y^* & \quad\lambda-1+y^*\quad & b_2y^*\\ -c_1z^* & -c_2z^* & \lambda+d_2z^*-1\\ \end{array} \right|, \end{eqnarray*} $

其中$ E $为三阶单位矩阵,得到

(3.4) $ \begin{eqnarray} &&\lambda^{3}+(x^*+y^*+d_{2}z^*-3)\lambda^{2}+[(d_{2}+b_{1}c_{1})x^*z^*+(d_{2}+b_{2}c_{2})y^*z^*\\ &&+(1-a_{1}a_{2})x^*y^*-2x^*-2y^*-2d_{2}z^*+3]\lambda+(d_{2}-a_{1}b_{2}c_{1}\\ &&-a_{2}b_{1}c_{2}+b_{1}c_{1}+b_{2}c_{2}-a_{1}a_{2}d_{2})x^*y^*z^*+(a_{1}a_{2}-1)x^*y^*\\ &&-(b_{1}c_{1}+d_{2})x^*z^*-(b_{2}c_{2}+d_{2})y^*z^*+x^*+y^*+d_{2}z^* = 0. \end{eqnarray} $

为了方便表示,令

$ \begin{eqnarray*} m& = &x^*+y^*+d_{2}z^*-3, \\ n& = &(d_{2}+b_{1}c_{1})x^*z^*+(d_{2}+b_{2}c_{2})y^*z^*+(1-a_{1}a_{2})x^*y^*-2x^*-2y^*-2d_{2}z^*+3, \\ w& = &(d_{2}-a_{1}b_{2}c_{1}-a_{2}b_{1}c_{2}+b_{1}c_{1}+b_{2}c_{2}-a_{1}a_{2}d_{2})x^*y^*z^*+(a_{1}a_{2}-1)x^*y^*\\ &&-(b_{1}c_{1}+d_{2})x^*z^*-(b_{2}c_{2}+d_{2})y^*z^*+x^*+y^*+d_{2}z^*. \end{eqnarray*} $

则(3.4)式变为

$ \lambda^{3}+m\lambda^{2}+n\lambda+w = 0. $

可以用Matlab求得上式解,由于解的形式过于复杂,这里就不展开写了,以$ \lambda_{1} $、$ \lambda_{2} $、$ \lambda_{3} $表示上式的三个根.根据不动点理论,当$ \mid\lambda_{1}\mid < 1 $、$ \mid\lambda_{2}\mid < 1 $、$ \mid\lambda_{3}\mid < 1 $时,不动点$ (x^*, y^*, z^*) $是稳定的.即当三个种群的初始数量位于封闭图像中时,种群的数量最终将保持稳定,反之,种群数量将趋向于无限发展.

系统的稳定性在生态系统的分析中起着重要的作用,稳定的生态系统对于保持完整的系统功能是非常重要的.近年来随着工业的发展,生态平衡受到严重破坏,为了保持生态系统稳定,人为地对生态系统进行控制是十分重要的.我们可以根据客观问题的实际需要,通过设计控制系统,消除有害因素,引入有用的控制,有效地约束描述系统初始个数的Julia集,达到控制的效果.

将线性反馈控制器引入到系统(3.2)中,得到被控系统如下

(4.1) $ \begin{eqnarray} \left \{\begin{array}{l} x_{n+1} = 2x_{n}-x_{n}^2-a_1 x_{n}y_{n}-b_1 x_{n}z_{n}+k_{1}(x_{n}-x_{ref}), \\ y_{n+1} = 2y_{n}-y_{n}^2-a_2 x_{n}y_{n}-b_2 y_{n}z_{n}+k_{2}(y_{n}-y_{ref}), \\ z_{n+1} = (1-d_1)z_{n}-d_2z_{n}^2+c_1 x_{n}z_{n}+c_2y_{n}z_{n}+k_{3}(z_{n}-z_{ref}). \end{array} \right. \end{eqnarray} $

其中$ x_{ref} $、$ y_{ref} $、$ z_{ref} $为控制项, $ k_{1} $、$ k_{2} $、$ k_{3} $是反馈控制参数.显然,系统(3.2)的不动点也是系统(4.1)的不动点,这意味着它必须满足以下方程

$ \begin{eqnarray*} \label{system5} \left \{\begin{array}{l} x^* = 2x^*-(x^*)^2-a_1x^*y^*-b_1x^*z^*+k_{1}(x^*-x_{ref}), \\ y^* = 2y^*-(y^*)^2-a_2x^*y^*-b_2y^*z^*+k_{2}(y^*-y_{ref}), \\ z^* = (1-d_1)z^*-d_2(z^*)^2+c_1x^*z^*+c_2y^*z^*+k_{3}(z^*-z_{ref}). \end{array} \right. \end{eqnarray*} $

因此,辅助参考项$ x_{ref} $、$ y_{ref} $、$ z_{ref} $必须被视为系统的不动点$ (x^*, y^*, z^*) $.故竞争控制模型如下

(4.2) $ \begin{eqnarray} \left \{\begin{array}{l} x_{n+1} = 2x_{n}-x_{n}^2-a_1 x_{n}y_{n}-b_1 x_{n}z_{n}+k_{1}(x_{n}-x^*), \\ y_{n+1} = 2y_{n}-y_{n}^2-a_2 x_{n}y_{n}-b_2 y_{n}z_{n}+k_{2}(y_{n}-y^*), \\ z_{n+1} = (1-d_1)z_{n}-d_2z_{n}^2+c_1 x_{n}z_{n}+c_2y_{n}z_{n}+k_{3}(z_{n}-z^*). \end{array} \right. \end{eqnarray} $

系统(4.2)在$ (x^*, y^*, z^*) $处的Jacobi矩阵为

$ \begin{eqnarray*} J = \left[ \begin{array}{ccc} 1-x^*+k_{1} & -a_1x^* & -b_1x^*\\ -a_2y^* &\ 1-y^*+k_{2}\ & -b_2y^*\\ c_1z^* & c_2z^* & -d_2z^*+1+k_{3}\\ \end{array} \right]. \end{eqnarray*} $

令

$ \begin{eqnarray*} &\mid\lambda E-J\mid = \left| \begin{array}{ccc} \lambda-1+x^*-k_1 & a_1x^* & b_1x^*\\ a_2y^* &\ \lambda-1+y^*-k_2\ & b_2y^*\\ -c_1z^* & -c_2z^* & \lambda+d_2z^*-1-k_3\\ \end{array} \right|, \end{eqnarray*} $

其中$ E $为三阶单位矩阵,得到

(4.3) $ \begin{eqnarray} &&\lambda^3+(x^*+y^*+d_2z^*-k_1-k_2-k_3-3)\lambda^2+[(1-a_1a_2)x^*y^*+(b_1c_1+d_2)x^*z^*\\ &&+(b_2c_2+d_2)y^*z^*-(k_2+k_3+2)x^*-(k_1+k_3+2)y^*-(k_1+k_2+2)d_2z^*+ 2(k_1+k_2+k_3)\\ &&+k_1k_2+k_1k_3+k_2k_3+3]\lambda+(b_1c_1+b_2c_2-a_1a_2d_2-a_1b_2c_1-a_2b_1c_2+d_2)x^*y^*z^*\\ &&-(1-a_1a_2)(k_3+1)x^*y^*-(b_1c_1+d_2)(k_2+1)x^*z^*-(b_2c_2+d_2)(k_1+1)y^*z^*\\ &&+(1+k_2+k_3+k_2k_3)x^*+(1+k_1+k_3+k_1k_3)y^*+(1+k_1+k_2+k_1k_2)d_2z^*\\ &&-(k_1+k_2+k_3+k_1k_2+k_1k_3+k_2k_3+k_1k_2k_3+1) = 0. \end{eqnarray} $

为了方便表示,令

$ \begin{eqnarray*} m'& = &x^*+y^*+d_2z^*-k_1-k_2-k_3-3, \\ n'& = &(1-a_1a_2)x^*y^*+(b_1c_1+d_2)x^*z^*+(b_2c_2+d_2)y^*z^*-(k_2+k_3+2)x^* \\ &&-(k_1+k_3+2)y^*-(k_1+k_2+2)d_2z^*+2(k_1+k_2+k_3)+k_1k_2+k_1k_3+k_2k_3+3, \\ w'& = &(b_1c_1+b_2c_2-a_1a_2d_2-a_1b_2c_1-a_2b_1c_2+d_2)x^*y^*z^*-(1-a_1a_2)(k_3+1)x^*y^*\\ &&-(b_1c_1+d_2)(k_2+1)x^*z^*-(b_2c_2+d_2)(k_1+1)y^*z^*+(1+k_2+k_3+k_2k_3)x^*\\ &&+(1+k_1+k_3+k_1k_3)y^*+(1+k_1+k_2+k_1k_2)d_2z^*\\ &&-(k_1+k_2+k_3+k_1k_2+k_1k_3+k_2k_3+k_1k_2k_3+1). \end{eqnarray*} $

则(4.3)式变为

$ \lambda^{3}+m'\lambda^{2}+n'\lambda+w' = 0. $

由于解的形式的复杂性,以$ \lambda_{1}' $、$ \lambda_{2}' $、$ \lambda_{3}' $表示上式的三个根.因此,我们可以通过调整反馈控制参数$ k_{1} $、$ k_{2} $、$ k_{3} $,使$ \mid\lambda_{1}'\mid < 1 $、$ \mid\lambda_{2}'\mid < 1 $、$ \mid\lambda_{3}'\mid < 1 $,从而实现系统的稳定性控制.比如取$ a_1 = 0.97, a_2 = 0.81, b_1 = 0.39, b_2 = 0.85, c_1 = 0.23, c_2 = 0.13, d_1 = 0.17, d_2 = 0.01 $.那么共存不动点为$ (0.4857, 0.461, \; 0.1739) $,取$ k_1 = -0.91, k_2 = -0, 7, k_3 = -0.49 $,则共存不动点对应的特征值为$ -0.7, 0.15, 0.5 $.根据Jury判据,离散系统的闭环特征方程为

$ D(z) = z^3+m'z^2+n'z+w' = 0. $

由于$ n = 3 $,故Jury表为$ 3 $行、$ 4 $列.

其中

$ \begin{eqnarray*} &&a_0 = w' = 0.0528, a_1 = n' = -0.3787, a_2 = m' = 0.0484, a_3 = 1.\\ &&{b_0 = {\left[\begin{array}{cc} a_0\quad&a_3\\ a_3\quad&a_0 \end{array} \right]} = a_{0}^2-a_{3}^2 = {w'}^2-1 = -0.9972, }\\ &&{b_1 = {\left[\begin{array}{cc} a_0\quad &a_2\\ a_3\quad &a_1 \end{array} \right]} = a_{0}a_{1}-a_{2}a_{3} = w'n'-m' = -0.0684, }\\ &&{b_2 = {\left[\begin{array}{cc} a_0\quad &a_1\\ a_3\quad &a_2 \end{array} \right]} = a_{0}a_{2}-a_{1}a_{3} = w'm'-n' = -0.3813.} \end{eqnarray*} $

根据Jury判据

$ \begin{eqnarray*} (1)&&D(1) = 1^3+w'+m'+n' = 0.7225>0, \\ (2)&&D(-1) = (-1)^3+w'+m'-n' = -0.52<0, \\ (3)&&|a_0|<a_3, |b_0|>|b_2|. \end{eqnarray*} $

上述条件满足,故系统是稳定的.

对应的分形图如图 2.通过调整反馈控制参数$ k_i(i = 1, 2, 3) $,可以将原本位于Julia集之外的点引入Julia集,从而扩大了使系统保持稳定的种群初始数量的选择范围.

同时,我们可以控制单一变量,观察某一种群的数量变化对另外两个种群以及整个生态系统的影响.保持$ a_1 = 0.97, a_2 = 0.81, b_1 = 0.39, b_2 = 0.85, c_1 = 0.23, c_2 = 0.13, d_1 = 0.17, $$ d_2 = 0.01 $不变,令$ k_1 = k_2 = -0.0987, k_3 = -0.9423 $及$ k_1 = k_2 = -0.0987, k_3 = -0.0065 $,仿真如图 3(a)–(b),可以看出随着$ k_3 $的增大, Julia集整体呈收缩趋势.也就是说,捕食者数量增大,会导致食饵数量减少,从而加剧了生态系统的不稳定性.同理, $ k_1 = k_3 = -0.9772, k_2 = -0.8931 $及$ k_1 = k_3 = -0.9772, k_2 = -0.0346 $时,仿真如图 4(a)–(b),可以看出随着$ k_2 $的增大, Julia集缩小; $ k_2 = k_3 = -0.0788, k_1 = -0.994 $及$ k_2 = k_3 = -0.0788, k_1 = -0.1273 $时,仿真如图 5(a)–(b),可以看出随着$ k_1 $的增大, Julia集呈收缩趋势.

众所周知,非线性动力系统的同步是一个重要而有趣的研究课题,它广泛应用于力学、通信等领域,人们对此进行研究并取得了很多成果.在这一节中,我们将讨论两种形式相同但参数不同的的Julia集的同步问题.根据Julia集的定义,一旦给出了系统的参数,其相应的Julia集就是确定的和唯一的,并且它们之间没有相关性.通过引入非线性耦合控制项,使得系统(3.2)变成另一个系统,实现两个系统的耦合.

考虑与系统(3.2)形式相同但参数不同的系统

(5.1) $ \begin{eqnarray} \left \{\begin{array}{l} p_{n+1} = 2p_{n}-p_{n}^2-\alpha_1 p_{n}q_{n}-\beta_1 p_{n}w_{n}, \\ q_{n+1} = 2q_{n}-q_{n}^2-\alpha_2 p_{n}q_{n}-\beta_2 q_{n}w_{n}, \\ w_{n+1} = (1-\xi_1)w_{n}-\xi_2w_{n}^2+\gamma_1 p_{n}w_{n}+\gamma_2 q_{n}w_{n}. \end{array} \right. \end{eqnarray} $

为了实现系统(3.2)与系统(5.1)的同步,引入非线性耦合控制项,得到耦合系统如下

$ \begin{eqnarray*} \left \{\begin{array}{rl} p_{n+1} = &2p_{n}-p_{n}^2-\alpha_1 p_{n}q_{n}-\beta_1 p_{n}w_{n}-m_1(p_n-x_n)-h_1(x^2_n-p^2_n)\\ &-h_1(a_1x_ny_n-\alpha_1p_nq_n)-h_1(b_1x_nz_n-\beta_1p_nw_n), \\ q_{n+1} = &2q_{n}-q_{n}^2-\alpha_2 p_{n}q_{n}-\beta_2 q_{n}w_{n}-m_2(q_n-y_n)-h_2(y^2_n-q^2_n)\\ &-h_2(a_2x_ny_n-\alpha_2p_nq_n)-h_2(b_2y_nz_n-\beta_2q_nw_n), \\ w_{n+1} = &(1-\xi_1)w_{n}-\xi_2w_{n}^2+\gamma_1 p_{n}w_{n}+\gamma_2 q_{n}w_{n}-h_3[(1-\xi_1)w_n-(1-d_1)z_n]\\ &-h_3(d_2z^2_n-\xi_2w^2_n)-h_3(\gamma_1p_{n}w_{n}-c_1x_nz_n)-h_3(\gamma_2 q_{n}w_{n}-c_2y_nz_n). \end{array} \right. \end{eqnarray*} $

定理5.1  当$ n\rightarrow\infty, \; m_{1}\rightarrow2, \; h_{1}\rightarrow1 $时, $ |p_{n+1}-x_{n+1}|\rightarrow0 $.

证  由于Julia集以及系统参数的有界性,存在$ M_1 > 0 $满足$ |p^2_n|+|x^2_n| < M_1 $, $ |a_1||x_ny_n|+|\alpha_1p_nq_n| < M_1 $, $ |b_1||x_nz_n|+|\beta_1||p_nw_n| < M_1 $.则

$ \begin{eqnarray*} &&|p_{n+1}-x_{n+1}|\\ &\leq&|2-m_1||p_n-x_n|+|1-h_1|(|x^2_n-p^2_n|+|a_1x_ny_n-\alpha_1 p_nq_n|+|b_1x_nz_n-\beta_1p_nw_n|)\\ &\leq&|2-m_1||p_n-x_n|+|1-h_1|(|x^2_n|+|p^2_n|+|a_1||x_ny_n|+|\alpha_1|| p_nq_n|+|b_1||x_nz_n|+|\beta_1||p_nw_n|)\\ &\leq&|2-m_1||p_n-x_n|+3|1-h_1|M_1\\ &\leq&|2-m_1|^2|p_{n-1}-x_{n-1}|+3|1-h_1|M_1(1+|2-m_1|)\\ &\leq&\cdots\\ &\leq&|2-m_1|^n|p_1-x_1|+3|1-h_1|M_1\frac{1-|2-m_1|^n}{1-|2-m_1|}. \end{eqnarray*} $

当$ m_1\rightarrow2, \; h_1\rightarrow1 $时, $ |p_{n+1}-x_{n+1}|\rightarrow0 $.同理,当$ n\rightarrow+\infty $, $ |q_{n+1}-y_{n+1}|\rightarrow0 $.

证毕.

定理5.2  当$ n\rightarrow\infty, \; h_{3}\rightarrow1 $时, $ |w_{n+1}-z_{n+1}|\rightarrow0 $.

证  由于Julia集以及系统参数的有界性,存在$ M_2 > 0 $满足$ |d_2||z^2_n|+|\xi_2||w^2_n| < M_3 $, $ |\gamma_1||p_nw_n|+|c_1||x_nz_n| < M_3 $, $ |\gamma_2||q_nw_n|+|c_2||y_nz_n| < M_3 $.则

$ \begin{eqnarray*} &&|w_{n+1}-z_{n+1}|\\ &\leq&|1-h_3|(|(1-\xi_1)w_n-(1-d_1)z_n|+|d_2z^2_n-\xi_2w^2_n|+|\gamma_1p_nw_n-c_1x_nz_n|\\ &&+|\gamma_2q_nw_n-c_2y_nz_n|)\\ &\leq&|1-h_3|(|(1-\xi_1)w_n-(1-d_1)z_n|+|d_2||z^2_n|+|\xi_2||w^2_n|+|\gamma_1||p_nw_n|\\ &&+|c_1||x_nz_n|+|\gamma_2||q_nw_n|+|c_2||y_nz_n|)\\ &\leq&|1-h_3|(|(1-\xi_1)w_n-(1-d_1)z_n|+3M_3). \end{eqnarray*} $

易知,当$ h_3\rightarrow1 $, $ |w_{n+1}-z_{n+1}|\rightarrow0 $.

根据Julia集的定义和上述定理,耦合控制系统(5.1)的Julia集变为系统(3.2)的Julia集,从而实现系统(3.2)和系统(5.1)的Julia集的同步.虽然在给定参数时确定了Julia集,但通过采用不同的耦合参数,我们可以看到Julia集的变化过程.图 1是系统(3.2)中参数取$ a_1 = 0.97, \; a_2 = 0.81, \; b_1 = 0.39, \; b_2 = 0.85, \; c_1 = 0.23, \; c_2 = 0.13, \; d_1 = 0.17, d_2 = 0.01 $时的Julia集, 图 6是引入耦合控制项后Julia集的变化.系统(5.1)随着控制参数的变化而变化,很明显,系统(5.1)的Julia集正在变为图 1所示的系统(3.2)的Julia集.

种群竞争模型不仅是理论研究的重要模型,更是实际应用的重要模型,模型的原始数量对种群数量的长期变化具有重要意义,因此可以用于指导生态、研究种群发展.与传统的李雅普诺夫分析方法相对比,分形方法使用更加直观、方便.李雅普诺夫方法是运用李氏判据处理稳定性问题,其关键步骤是构造李雅普诺夫函数$ V(x) $,但李雅普诺夫方法未能给出构造$ V(x) $的一般方法,而且至今仍未形成构造李雅普诺夫函数的通用方法.构造$ V $函数是一种技巧性的问题,多凭借研究工作者本人的经验,而这种经验是在长期工作和实践中积累形成的,这就给理论工作的研究、应用及推广带来了一定困难.分形图像刻画的是收敛域与发散域的边界,在收敛域内的点最终都趋于稳定,而在发散域的点最终都趋于无穷.因此运用分形方法分析时,只需要根据系统方程刻画出分形图像,就可以简洁直观的看出每一点的稳定性,即在分形区域内的点最终都趋于稳定,在实际应用中十分直观、方便.本文在二维生态系统的基础上引入三维生态系统,通过研究竞争模型的Julia集,使得集外的一些初值被引入到集内,这相当于控制一些不具有持久性的种群具有持久性,从而实现了竞争模型的稳定性控制,并分析了一个种群初始数量变化对另外两个种群数量及生态系统的影响.通过引入耦合项,实现了Julia集的同步,使得两种模型具有相同的行为,并且都具有持久性,仿真结果表明了方法的有效性.