引言

19 96年，Kopel提出了Kopel系统^[1]，它是一个二维离散古诺双寡头模型，其表达式为

其中，x_n，y_n分别为公司x和y在时间t内生产的物资数量，ρ与μ均为正参数，ρ可以理解为耦合参数，限制0＜ρ＜1.

现回顾Kopel系统的整个研究过程.Kopel系统最早于1996年被提出^[1].1999年，Agiza分析了其不动点的稳定性、混沌相图、分岔的数值模拟及混沌控制^[2].2005年，Anderson.等分析了该系统的吸引盆问题^[3].2006年，于晋臣研究了Kopel系统的跨临界分岔、倍周期分岔、叉式分岔及N-S分岔的数值模拟^[4].2008年，Govaerts.等研究了其不动点及周期轨的稳定性^[5].2009年，于晋臣等分析了系统的分岔及数值模拟^[6]，同年，于晋臣等^[7]给出其倍周期分岔与混沌.2010年，吴文娟等^[8，9]给出了系统的马蹄混沌、间歇混沌及混沌控制.2020年，李波等^[10] 研究了其双参数分支情形-1:4共振.Gao等^[11-18]分析了基于Kopel系统的推广模型的稳定性、分岔与混沌现象，这些分岔研究主要基于文献^[23-27] 的方法.

Murakami提出分析N-S分岔的新方法，并成功得到验证^[19，22].现基于此新方法，再次给出Kopel系统（1）的N-S分岔的证明，得到分岔方向及稳定性结果，并计算出N-S分岔所产生的不变曲线的近似表达式，同时，对N-S分岔的扰动参数与不变曲线中的参数之间的关系进行了分析.这些是对Kopel系统（1）的内在复杂动力学研究的进一步补充完善.

1 Kopel系统的N-S分岔及不变曲线

由离散系统的不动点定义，令

$\left\{\begin{array}{c}(1-\rho )x_n+\rho \mu y_n\left(1-y_n\right)=x_n,\\ (1-\rho )y_n+\rho \mu x_n\left(1-x_n\right)=y_n,\end{array}\right.$

得Kopel系统的四个非负不动点^[1-8]，分别记为

$\begin{array}{c}{P}_1(x*,y*)=\left(\mathrm{0,\; 0}\right),\\ P_2(x*,y*)=\left(\frac{\mu -1}{\mu },\frac{\mu -1}{\mu }\right),\mu >1,\\ P_3(x*,y*)=\left(P_{31},P_{32}\right),\end{array}$

其中，

$\begin{array}{c}{P}_{31}=\frac{\mu +1+\sqrt{(\mu +1)(\mu -3)}}{2\mu },\\ P_{32}=\frac{\mu +1-\sqrt{(\mu +1)(\mu -3)}}{2\mu },\end{array}$

μ≥3，不动点P₄与P₃关于直线x=y对称.文献^[4，6]取定一组特殊参数值，给出了P₃不动点处的N-S分岔数值模拟，但对此系统N-S分岔的不变曲线的研究仍未涉及.

系统（1）在不动点处的雅可比矩阵为

其特征方程为

$\begin{array}{c}{\lambda }^2-2(1-\rho )\lambda +(1-\rho {)}^2-\\ (\rho \mu -2\rho \mu x*)(\rho \mu -2\rho \mu y*)=0\end{array}$

则特征值

$\begin{array}{c}{\left.{\lambda }_{\mathrm{1,\; 2}}\right|}_{(x*,y*)}=(1-\rho )\pm \\ \rho \mu \sqrt{(1-2x*)(1-2y*)}.\end{array}$

现在分析P₃处的N-S分岔及分岔所产生的不变曲线的近似表达式.

引理1^[22] 给定二维方阵A，如果迹trA及行列式det A满足|trA|-1≤detA =1，则A的特征值为λ=e^±iω，其中ω=cos^-1（trA/2），i为虚数单位.

引理2 如果0＜ρ＜1，μ＞1+ $\sqrt{6}$，则当ρ=2/（μ²-2μ-3）时，系统（1）在不动点P₃处的雅可比矩阵有一对共轭纯虚数特征值λ=e^±iω，其中，ω=cos^-1（1-ρ），并且λ满足

$\text{(}\text{i}\text{)}{\lambda }^j\ne 1(j=\mathrm{1,\; 2},\mathrm{3,\; 4})$

$\text{(ii)}{\left.\frac{\mathrm{d}\left|\lambda \right|}{\mathrm{d}\rho }\right|}_{\rho =\frac{2}{\mu 2-2\mu -3}}=1>0\text{.}$

证明:（i）Kopel系统（1）在不动点P₃处的雅可比矩阵为

其中

$\begin{array}{c}{A}_{12}=-\rho [1-\sqrt{(\mu +1)(\mu -3)}],\\ A_{21}=-\rho [1+\sqrt{(\mu +1)(\mu -3)}],\end{array}$

则$\mathrm{t}\mathrm{r}\mathit{A}=2（1-\rho ），\mathrm{d}\mathrm{e}\mathrm{t}\mathit{A}=\left({\mu }^2-2\mu -3\right){\rho }^2-2\rho +1$，由0＜ρ＜1及$\rho =2/\left({\mu }^2-2\mu -3\right)$，故|trA|-1≤1，且detA =1，又$\mu >1+\sqrt{6}$，也保证了0＜ρ＜1，即此引理的三个条件是不冲突的，是相容的.因此，由引理1，系统（1）在P₃处的雅可比矩阵有一对共轭纯虚数特征值λ=e^±iω，其中ω=cos^-1（1-ρ），又${\lambda }^j={\left({\mathrm{e}}^{\pm \mathrm{i}\omega }\right)}^j=\mathrm{c}\mathrm{o}\mathrm{s}（j\omega ）\pm \mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}（j\omega ）$，因 0＜cos ω=1-ρ＜1，故$\omega \ne 0，\frac{\pi }{2}，\frac{2\pi }{3}，\pi$，从而cos（jω）≠1，得λ^j≠1（j=1，2，3，4）.

（ii）λ=e^±iω是一对共轭纯虚数，故$|\lambda {|}^2=\lambda \stackrel{-}{\lambda }=\mathrm{d}\mathrm{e}\mathrm{t}\mathit{A}$，得$\left|\lambda \right|=（\mathrm{d}\mathrm{e}\mathrm{t}\mathit{A}{）}^{\frac{1}{2}}$，又当$\rho =\frac{2}{{\mu }^2-2\mu -3}$时，detA =1.又

$\mathrm{d}\mathrm{e}\mathrm{t}\mathit{A}=\left({\mu }^2-2\mu -3\right){\rho }^2-2\rho +1$

故$\frac{\mathrm{d}}{\mathrm{d}\rho }（\mathrm{d}\mathrm{e}\mathrm{t}\mathit{A}）=2\left({\mu }^2-2\mu -3\right)\rho -2$，从而

证毕.

由文献^[24]的定理3.5.2可知，当0＜ρ＜1，$\mu >1+\sqrt{6}$及$\rho =\frac{2}{{\mu }^2-2\mu -3}$时，Kopel系统（1）在不动点P₃处已经满足了文献^[24]的定理3.5.2 的（SH1）与（SH2）条件，现Kopel系统（1）可能会发生N-S分岔.同时，存在光滑的坐标变换h，采用文献^[24]的记号，使得系统（1）可转化为极坐标形式$\left(r+d\left(\mu -{\mu }_0\right)r+a{r}^3，\theta +c+b{r}^2\right)$ +高阶项.

还需验证条件（SH3）:r³的系数a≠0.当0＜ρ＜1，$\mu >1+\sqrt{6}$及$\rho =\frac{2}{{\mu }^2-2\mu -3}$时，对μ与ρ取一般值，由文献^[23，24]的方法或文献^[19，22]的方法计算a，即使利用Mathematica软件计算r³的系数都非常困难.下面将在一组特殊参数下，利用文献^[19，22]新方法分析N-S分岔并给出不变曲线的近似表达式，有别于文献^[23，24]的方法.

引理3 如果$\rho =\frac{2}{{\mu }^2-2\mu -3}$，则系统（1）在不动点P₃处的雅可比矩阵的特征值λ=e^±iω所对应的特征向量为

$\mathit{q}={\left(-\frac{\mathrm{i}\sqrt{{\mu }^2-2\mu -4}}{1+\sqrt{(\mu +1)(\mu -3)}},1\right)}^{\mathrm{T}}$

且q与p满足$\mathit{A}q=\lambda q，\mathit{p}\mathit{A}=\lambda p\text{与}pq=2，\omega ={\mathrm{c}\mathrm{o}\mathrm{s}}^{-1}（1-\rho ）$.

证明:矩阵A的特征方程

${\lambda }^2-2（1-\rho ）\lambda +\left({\mu }^2-2\mu -3\right){\rho }^2-2\rho +1=0$，其中

$\begin{array}{c}\lambda ={\mathrm{e}}^{\mathrm{i}\omega }=cos\left(\omega \right)+isin\left(\omega \right)=\\ (1-\rho )+i\sqrt{1-(1-\rho {)}^2}=\\ \frac{{\mu }^2-2\mu -5+2i\sqrt{{\mu }^2-2\mu -4}}{{\mu }^2-2\mu -3},\\ \end{array}$

易验证Aq=λq、pA=λp且ρ.此处ρ与ρ分别是A的右特征向量与左特征向量.证毕.

现基于文献 ^[19，22]的新方法，选取μ=5，$\rho =\frac{2}{{\mu }^2-2\mu -3}=\frac{1}{6}$的特殊情形进行分析.此组参数情形下，引理2与引理3仍是成立的.

定理1 当ρ$=\frac{1}{6}$且μ=5时，系统（1）在不动点P₃处发生N-S分岔.对参数ρ进行扰动ρ+ε₁，ε₁是充分小的扰动参数，那么，当ε₁＞0时，存在一个稳定的闭不变曲线，且其近似表达式为

其中

$\begin{array}{c}\left(\begin{array}{c}{x}^{*}\\ y^{*}\end{array}\right)=\left(\begin{array}{c}\frac{3+\sqrt{3}}{5}\\ \frac{3-\sqrt{3}}{5}\end{array}\right),\\ q=\left(\begin{array}{c}-\frac{\sqrt{11}}{1+2\sqrt{3}}i\\ 1\end{array}\right),\\ r_0=\sqrt{\frac{\epsilon }{7.1851}},\end{array}$

${\mathit{K}}_{20}=\left(\begin{array}{c}\frac{5\left(171+23\sqrt{3}+(45\sqrt{11}+35\sqrt{33}{)}_{\mathrm{i}}\right)}{48(13+4\sqrt{3})}\\ \frac{5\left(-177-35\sqrt{3}+(25\sqrt{33}-15\sqrt{11}{)}_{\mathrm{i}}\right)}{48(13+4\sqrt{3})}\end{array}\right),$

${\mathit{K}}_{11}=\left(\begin{array}{c}\frac{5+65\sqrt{3}}{-39-12\sqrt{3}}\\ \frac{65+75\sqrt{3}}{39+12\sqrt{3}}\end{array}\right),$

ε是一个适当小的正数，θ为旋转角度.

证明:由于ρ$=\frac{1}{6}$且μ=5，因此，引理2及引理3的结论都成立.现在进一步分析N-S分岔及分岔所产生的不变曲线的近似表达式.

首先，把不动点P₃平移到原点，作变换$u_n=x_n-x^{*}，v_n=y_n-y^{*}$，且

$\left\{\begin{array}{c}(1-\rho )x^{*}+\rho \mu y^{*}\left(1-y^{*}\right)=x^{*},\\ (1-\rho )y^{*}+\rho \mu x^{*}\left(1-x^{*}\right)=y^{*}.\end{array}\right.$

得

由引理2及分岔理论，系统（7）可转化为规范型

在极坐标下，式（8）可写成

其中

$a\left(\epsilon \right)=\mathrm{R}\mathrm{e}\left(\frac{c\left(\epsilon \right)}{a\left(\epsilon \right)}\right),b\left(\epsilon \right)=\mathrm{I}\mathrm{m}\left(\frac{c\left(\epsilon \right)}{a\left(\epsilon \right)}\right).$

对（9）式的第一个方程的系数在ε=0处进行Taylor展开，得

下面计算立方项的系数a（0），判定其正负.记（7）式的非线性项为

则（7）式可简记为

其中$\mathit{A}=\left(\begin{array}{cc}1-\rho & \rho \mu \left(1-2y^{*}\right)\\ \rho \mu \left(1-2x^{*}\right)& 1-\rho \end{array}\right)$.

定义矩阵$\mathbf{\Phi }=（q，\stackrel{-}{q}）$，作变换

其中

${\mathrm{\Phi }}_{11}=\frac{-\mathrm{i}z\sqrt{{\mu }^2-2\mu -4}+\mathrm{i}\stackrel{-}{z}\sqrt{{\mu }^2-2\mu -4}}{1+\sqrt{(\mu +1)(\mu -3)}}.$

把（13）式代入（11）式，当z=0时，也有

$\mathit{F}\left(\mathbf{\Phi }\left(\begin{array}{c}z\\ z\end{array}\right)\right)=\left(\begin{array}{c}-\rho \mu \left[z^2+2z\stackrel{-}{z}+{\stackrel{-}{z}}^2\right]\\ F_{21}\end{array}\right)$

得

其中

$F_{21}=\rho \mu \frac{\left({\mu }^2-2\mu -4\right)\left[z^2-2z\stackrel{-}{z}+{\stackrel{-}{z}}^2\right]}{[1+\sqrt{(\mu +1)(\mu -3)}{]}^2}$

因此

其中

$\begin{array}{c}{f}_{201}=\\ \frac{4\mu \left(-4-2\mu +{\mu }^2\right)}{(-3+\mu )(1+\mu )\left(1+\sqrt{(-3+\mu )(1+\mu ){)}^2}\right.}\end{array}$

又

其中

$\begin{array}{c}{f}_{112}=\\ \frac{-4\mu \left(-4-2\mu +{\mu }^2\right)}{(-3+\mu )(1+\mu )\left(1+\sqrt{(-3+\mu )(1+\mu ){)}^2}\right.}\end{array}$

且

$\begin{array}{c}{\mathit{f}}_{21}=\\ {\left.\frac{{\partial }^3}{\partial z^2{\partial }_z^{-}}\mathit{F}\left(\mathbf{\Phi }\left(\begin{array}{c}z\\ z\end{array}\right)+\frac{K_{20}{z}^2}{2}+K_{11}z\stackrel{-}{z}+\frac{K_{02}{\stackrel{-}{z}}^2}{2}\right)\right|}_{z=0}\end{array}$

现可计算

得

由于f₂₁/2是（15）式

$\mathit{F}\left(\mathbf{\Phi }\left(\begin{array}{c}z\\ \stackrel{-}{z}\end{array}\right)+\frac{{\mathit{K}}_{20}{z}^2}{2}+{\mathit{K}}_{11}z\stackrel{-}{z}+\frac{{\mathit{K}}_{02}{\stackrel{-}{z}}^2}{2}\right)$

在（z，$\stackrel{-}{z}$）=（0，0）处展开式z²$\stackrel{-}{z}$的系数，因此实际计算过程中，取

把μ=5，ρ=16代入（7）式、（8）式及（11）~（19）式，得

$\begin{array}{c}{f}_{21}=\\ \left(\begin{array}{c}\frac{25(-239-445\sqrt{3}+(15\sqrt{11}-25\sqrt{33}\left)i\right)}{144(13+4\sqrt{3})}\\ \frac{25(495+385\sqrt{3}-(203\sqrt{11}+439\sqrt{33}\left)i\right)}{144(37+30\sqrt{3})}\end{array}\right)\end{array}$

因此有

$\begin{array}{c}c\left(0\right)=\frac{1}{2}\frac{1}{pq}p{\mathit{f}}_{21}=\frac{1}{2}\times \frac{1}{2}p{\mathit{f}}_{21}=\\ \frac{125(90+37\sqrt{3})(55-89\sqrt{11}\mathrm{i})}{792(217+104\sqrt{3})}\end{array}$

得

因此，由引理2及（20）式，结合文献^[19，22]的定理，Kopel系统（1）在不动点P₃处发生N-S分岔.对充分小的扰动ε₁＞0，存在稳定的不变曲线，且不动点P₃是渐近稳定的.

基于以上分析及文献^[19，22]的结论，得N-S分岔的不变曲线的近似表达式为

当取参数μ=5，ρ=1/6时，计算可得

$\begin{array}{c}{r}_0=\sqrt{-\frac{d\epsilon }{\alpha \left(0\right)}}=\sqrt{\frac{\epsilon }{7.1851}},\\ q=\left(\begin{array}{c}-\frac{\sqrt{11}}{1+2\sqrt{3}}i\\ 1\end{array}\right),\end{array}$

$\begin{array}{c}{\mathit{K}}_{20}=\\ \left(\begin{array}{c}\frac{5(171+23\sqrt{3}+(45\sqrt{11}+35\sqrt{33}\left)i\right)}{48(13+4\sqrt{3})}\\ \frac{5(-177-35\sqrt{3}+(25\sqrt{33}-15\sqrt{11}\left)i\right)}{48(13+4\sqrt{3})}\end{array}\right),\end{array}$

${\mathit{K}}_{11}=\left(\begin{array}{c}\frac{5+65\sqrt{3}}{-39-12\sqrt{3}}\\ \frac{65+75\sqrt{3}}{39+12\sqrt{3}}\end{array}\right)$，θ为旋转角度，存在一个适当小的正数ε，从而得到N-S分岔的不变曲线的具体近似表达式（6）.证毕.

注:闭不变曲线的近似表达式（6）中的ε是一个适当小的正数，它不一定等于扰动参数ε₁.

2 数值模拟

现对N-S分岔进行数值模拟，当取μ=5，ε₁=0.001，则扰动后ρ=1/6+ε₁，（6）式中的ε=ε₁=0.001，初始值取ε=ε₁=0.001，初始值的选取不能远离不动点P₃，迭代次数2000次，去掉部分暂态数据，得N-S分岔图及近似不变曲线，如图1所示，其中心的星号表示P₃不动点. 近似不变曲线为闭虚线，有很好的近似效果.

当取μ=5，ε₁=0.0001，则扰动后ρ=1/6+ε₁，（6）式中的ε=ε₁=0.0001时，初始值取（0.3，0.2），得到数值模拟结果如图2所示，此时ε=0.0001所得的近似不变曲线已经不能近似N-S分岔的结果.而若选取ε₁=0.001≠ε时，才能得到图1类似的近似效果.

图1与图2两种情形也说明了当系统的参数发生小扰动ε₁时，会发生N-S分岔，则存在某个ε，可构造N-S分岔所产生的不变曲线近似方程（6），但此处的ε不一定是与小扰动参数ε₁取同一个值.

现结合文献^[22]的定理4.4及其第916页的数值模拟例子印证以上观点.文献^[22]中所研究的系统为

其不动点P₃（β，α（1-β）-1）.当$\beta =\frac{1}{2}\left(1-\frac{1}{\alpha }\right)$时，系统（22）将发生N-S分岔，分岔所产生的近似不变曲线的近似表达式为

其中${\rho }_0=\frac{\sqrt{\mu }}{{\alpha }_0^2}，\mathit{q}、{\mathit{K}}_{20}、{\mathit{K}}_{11}$等取值情况见文献^[22]的定理4.4.

当α=2.5+μ，即α₀=2.5，扰动参数μ=0.1，β=0.3，取初始值（0.6，0.6）. 类似前面的处理技巧，把（23）式的参数${\rho }_0=\sqrt{\mu }/{\alpha }_0^2$中的μ重新记成μ₁，以示区分. 当μ=μ₁=0.1时，即为文献^[22]的定理4.4的数值模拟的例子，所得结果如图3所示，（23）式能很好地近似N-S分岔所产生的不变曲线.

当μ=μ₁=0.001时，其他参数不变，去掉部分暂态数据，此时所得近似效果已较差.而此时若选择某个恰当的μ₁=0.0019，则可取得较好的近似效果.

通过以上算例，文献^[19]的推论2.1与文献^[22]的定理4.4中，不变曲线的近似表达式方程为

如果在N-S分岔中，对系统的参数扰动是μ，那么近似表达方程式（24）中${\rho }_0=\sqrt{-d/a\mu }$处的μ不一定取参数扰动μ相等的值，而只是存在某个μ₁，使得当${\rho }_0=\sqrt{-d/a{\mu }_1}$时，有较好的近似效果，否则$\underset{\mu \to 0}{lim}{\rho }_0=\underset{\mu \to 0}{lim}\sqrt{-d/a{\mu }_1}=0$，其中a、d是与μ无关的固定值.因此，由（24）式，当扰动参数μ趋于0时，近似表达式会缩小至不动点x*，而不是N-S分岔所产生的不变曲线的近似，只有找到恰当的μ₁才能实现较好的近似.

3 结论

首先阐述了Kopel系统的分岔、混沌与稳定性等方面已有的研究成果.在分岔方面，由于推导计算的复杂性，对N-S分岔的方向与稳定性的一般性证明仍是公开的问题.对N-S分岔所产生的不变曲线的研究也没有涉及.因此，本文采用Murakami K提出的新方法，分析了Kopel系统的N-S分岔，给出了部分一般性证明，即引理2与引理3，再分析了一组特殊参数（ρ=1/6，μ=5）情形下 Kopel系统的N-S分岔.同时，给出了N-S分岔所产生的不变曲线的近似表达式，讨论了扰动参数ε₁与近似表达式（6）中ε的选取问题.这些研究是对Kopel系统的复杂动力学的一个完善.但对其N-S分岔的一般性分析，即ρ=2/（μ²-2μ-3）与$\mu >1+\sqrt{6}$时，要克服（17）式、（19）式与（20）式的推导计算的困难，仍是将来需要解决的公开问题.Kopel系统的高余维分岔也是可研究问题，这些分岔与混沌的研究有助于理解Kopel系统的内在复杂动力学与经济现象.

参考文献