一类含参数的二维非线性系统的全局稳定性研究

Study on Global Stability of a Class of Two-Dimensional Nonlinear Systems with Parameters

摘要: 应用常微分定性理论，研究一类含参数的二维非线性系统的全局稳定性。首先分析系统随参数变化，其平衡点的个数和类型以及在平衡点处的稳定性发生什么样的改变，再研究该系统极限环的存在性、个数问题、稳定性及其位置，其次讨论了系统在无穷远点处的轨线情况，最后得到系统的全局相图。

Abstract: The global stability of a class of two-dimensional nonlinear systems with parameters is studied by using the ordinary differential qualitative theory. Firstly, the variation of the system with param-eters is analyzed, the number and type of the equilibrium point and the stability at the equilibrium point of the system are analyzed. Then the existence, number, stability and position of limit cycles of the system are studied. Secondly, the trajectory of the system at infinity is discussed and the global phase diagram is obtained.

1. 引言

在微分方程相关的一些课程中，我们发现求解一些微分方程比较简单，但一些微分方程的解你可能无法求出来，我们就需要对该系统进行一些定性理论分析，本文则是对一类含参数系统的解进行研究 [1] ，文献 [2] 中只对下面系统

$\{\begin{array}{l}\stackrel{\dot{}}{x}=\alpha x-y+x\left(x^2+y^2\right)=P\left(x,y\right)\\ \stackrel{\dot{}}{y}=x+\alpha y+y\left(x^2+y^2\right)=Q\left(x,y\right)\end{array}$ (1)

平衡点 [3] [4] 的稳定性 [5] [6] 和极限环 [7] 的存在进行了直接说明，缺少证明过程。但本文讨论了因参数 $\alpha$ 改变，该系统的平衡点发生的分叉 [8] 情况，以及对极限环的位置和稳定性进行了详细的说明，也给出了相应的证明过程，然后又讨论了系统在无穷远处 [1] 的情况，最后在参数的不同范围内分别取特定的数值，用相图 [9] 来描述不同情况下的轨线 [10] 走向情况，由此来验证前面理论部分的准确性。

2. 系统(1)的平衡点分析

定理1 当 $\alpha <0$ 时， $\left(0,0\right)$ 是系统(1)稳定的焦点；当 $\alpha \ge 0$ 时， $\left(0,0\right)$ 是系统(1)不稳定的焦点。

证明：

由平衡点的定义：

令 $\{\begin{array}{l}\alpha x-y+x\left(x^2+y^2\right)=0\\ x+\alpha y+y\left(x^2+y^2\right)=0\end{array}$，得 $\left(0,0\right)$ 是非线性系统(1)的平衡点。舍去方程(1)中非线性项，得到一个常系数线性方程：

$\{\begin{array}{l}\stackrel{\dot{}}{x}=\alpha x-y\\ \stackrel{\dot{}}{y}=x+\alpha y\end{array}$ (2)

其系数矩阵为 $A=\left(\begin{array}{cc}\alpha & -1\\ 1& \alpha \end{array}\right)$，其特征方程为 ${\lambda }^2-2\alpha \lambda +{\alpha }^2+1=0$，解得其特征值为 ${\lambda }_1=\alpha +i$， ${\lambda }_2=\alpha -i$。

1) 当 $\alpha <0$ 时， ${\lambda }_1,{\lambda }_2$ 为有负实部的虚数，所以平衡点 $\left(0,0\right)$ 是系统(1)稳定的焦点，如图1所示：

Figure 1. Track diagram of system (1) when α < 0

图1. α < 0时系统(1)的轨线图

2) 当 $\alpha >0$ 时， ${\lambda }_1,{\lambda }_2$ 为有正实部的虚数，所以平衡点 $\left(0,0\right)$ 是系统(1)不稳定的焦点。

3) 当 $\alpha =0$ 时， ${\lambda }_1=i,{\lambda }_2=-i$ 为纯虚数系统(2)以 $\left(0,0\right)$ 为中心，此时(1)就变

$\{\begin{array}{l}\stackrel{\dot{}}{x}=-y+x\left(x^2+y^2\right)\\ \stackrel{\dot{}}{y}=x+y\left(x^2+y^2\right)\end{array}$ (3)

用后继函数法研究此非线性方程组在 $\left(0,0\right)$ 的稳定性：在方程组(3)中，令

$x=r\cos \theta ,y=r\sin \theta$

有 $\{\begin{array}{l}\stackrel{\dot{}}{x}=\frac{\text{d}r}{\text{d}t}\cos \theta -r\sin \theta \frac{\text{d}\theta }{\text{d}t}=-r\sin \theta +r^3\cos \theta \\ \stackrel{\dot{}}{y}=\frac{\text{d}r}{\text{d}t}\sin \theta +r\cos \theta \frac{\text{d}\theta }{\text{d}t}=r\cos \theta +r^3\sin \theta \end{array}$

得： $\{\begin{array}{l}\frac{\text{d}r}{\text{d}t}=r^3\\ \frac{\text{d}\theta }{\text{d}t}=1\end{array}$ 消去t，得到

$\frac{\text{d}r}{\text{d}\theta }=r^3$ (4)

由于对初条件 $\theta =0$ 时， $r=0$，有整体解 $r\left(\theta \right)\equiv 0\left(-\infty <\theta <\infty \right)$，由文献 [2] 中第一章的定理4可知，对充分小的c， $\theta =0$ 时 $r=c$ 的解 $r\left(\theta ,c\right)$ 在 $\theta \in \left[-4\pi ,4\pi \right]$ 上有意义且 $r\in \left[0,r_1\right)$。

由于 $r\left(\theta ,c\right)$ 对c是解析的，所以可以展开成c的幂级数：

$r\left(\theta ,c\right)=r_1\left(\theta \right)c+r_2\left(\theta \right)c^2+\cdots$

并且由初始条件 $r\left(0,c\right)=c$ 得：

$r_1\left(0\right)=1$， $r_2\left(0\right)=r_3\left(0\right)=\cdots =0$

所以设解为 $r=c+r_2{c}^2+r_3{c}^3+\cdots$，其中 $r_2\left(0\right)=r_3\left(0\right)=\cdots =0$

将此解带入方程(4)，比较 $c^2$ 的系数得： $\frac{\text{d}{r}_2}{\text{d}\theta }=0$ 且 $r_2\left(0\right)=0$，所以 $r_2=0$

$r_2$ 是周期函数，再比较 $c^3$ 的系数得： $\frac{\text{d}{r}_3}{\text{d}\theta }=1$ 且 $r_3\left(0\right)=0$，所以 $r_3\left(\theta \right)=\theta$

因此 $r_3$ 不是周期函数

所以 $r\left(\theta ,c\right)=c+\theta c^3+r_4\left(\theta \right)c^4+\cdots$

$r\left(2\pi ,c\right)-r\left(0,c\right)=c+2\pi c^3+r_4\left(2\pi \right)c^4+\cdots -c-0\cdot c^3-r_4\left(0\right)c^4+\cdots =2\pi c^3+\cdots$

当c充分小时， $r\left(2\pi ,c\right)-r\left(0,c\right)>0\Rightarrow r\left(2\pi ,c\right)>r\left(0,c\right)$ 所以 $\left(0,0\right)$ 是系统(1)不稳定的焦点，如图2所示：

Figure 2. Track diagram of system (1) when α ≥ 0

图2. α ≥ 0时系统(1)的轨线图

证毕。

显然，当参数 $\alpha$ 由负变到正时， $\lambda$ 沿实轴上方或下方穿过虚轴，平衡点由稳定的焦点 $\left(\alpha <0\right)$ 变为不稳定的焦点 $\left(\alpha >0\right)$，这个系统在 $\alpha =0$ 处发生了Hopf分叉，如图3所示：

Figure 3. Hopf bifurcation

图3. Hopf分叉

该系统的分叉图如图4所示：

Figure 4. Subcritical bifurcation

图4. 亚临界分岔

从这个分叉图可以看出：当 $\alpha <0$，系统有一个稳定奇点和一个闭轨，而通过分叉点 $\alpha =0$ 后，奇点又变的不稳定了，则系统(1)的零解 $\left(x,y\right)=\left(0,0\right)$ 在 $\alpha =0$ 处发生亚临界分叉。

3. 系统(1)存在一个不稳定的极限环情况

定理2 当 $\alpha <0$ 时，系统存在一个不稳定的极限环

证明：

先讨论系统(1)在全平面上的极限环存在情况：

用Bendixson判断 [1] ，因为

$\frac{\partial P}{\partial x}+\frac{\partial Q}{\partial y}=2\alpha +4x^2+4y^2$

所以当参数 $\alpha \ge 0$，系统(1)在全平面上无闭轨。

下面考虑参数 $\alpha <0$ 的情况 [10]：

取 $D\left(x,y\right)=A{x}^m+B{y}^n+C{x}^p{y}^q{\text{e}}^{lx+gy}$

$\begin{array}{c}DP=A\alpha x^{m+1}-A{x}^{m}y+A{x}^{m+3}+A{x}^{m+1}{y}^2+B\alpha x{y}^n-B{y}^{n+1}+B{x}^3{y}^n+Bx{y}^{n+2}\\ +\alpha C{x}^{p+1}{y}^q{\text{e}}^{lx+gy}-C{x}^p{y}^{q+1}{\text{e}}^{lx+gy}+C{x}^{p+3}{y}^q{\text{e}}^{lx+gy}+C{x}^{p+1}{y}^{q+2}{\text{e}}^{lx+gy}\end{array}$

$\begin{array}{c}DQ=A{x}^{m+1}+A\alpha x^{m}y+A{x}^{m+2}y+A{x}^m{y}^3+Bx{y}^n+B\alpha y^{n+1}+B{x}^2{y}^{n+1}+B{y}^{n+3}\\ +C{x}^{p+1}{y}^q{\text{e}}^{lx+gy}+\alpha C{x}^p{y}^{q+1}{\text{e}}^{lx+gy}+C{x}^{p+2}{y}^{q+1}{\text{e}}^{lx+gy}+C{x}^p{y}^{q+3}{\text{e}}^{lx+gy}\end{array}$

$\begin{array}{c}\frac{\partial \left(DP\right)}{\partial x}=A\alpha \left(m+1\right)x^m-Am{x}^{m-1}y+A\left(m+3\right)x^{m+2}+A\left(m+1\right)x^m{y}^2+B\alpha y^n+3B{x}^2{y}^n+B{y}^{n+2}\\ +\alpha C\left(p+1\right)x^p{y}^q{\text{e}}^{lx+gy}-Cp{x}^{p-1}{y}^{q+1}{\text{e}}^{lx+gy}+C\left(p+3\right)x^{p+2}{y}^q{\text{e}}^{lx+gy}+C\left(p+1\right)x^p{y}^{q+2}{\text{e}}^{lx+gy}\\ +\alpha lC{x}^{p+1}{y}^q{\text{e}}^{lx+gy}-Cl{x}^p{y}^{q+1}{\text{e}}^{lx+gy}+Cl{x}^{p+3}{y}^q{\text{e}}^{lx+gy}+Cl{x}^{p+1}{y}^{q+2}{\text{e}}^{lx+gy}\end{array}$

$\begin{array}{c}\frac{\partial \left(DQ\right)}{\partial y}=A\alpha x^m+A{x}^{m+2}+3A{x}^m{y}^2+Bnx{y}^{n-1}+B\left(n+1\right)\alpha y^n+B\left(n+1\right)x^2{y}^n+B\left(n+3\right)y^{n+2}\\ +Cq{x}^{p+1}{y}^{q-1}{\text{e}}^{lx+gy}+\alpha C\left(q+1\right)x^p{y}^q{\text{e}}^{lx+gy}+C\left(q+1\right)x^{p+2}{y}^q{\text{e}}^{lx+gy}+C\left(q+3\right)x^p{y}^{q+2}{\text{e}}^{lx+gy}\\ +Cg{x}^{p+1}{y}^q{\text{e}}^{lx+gy}+\alpha gC{x}^p{y}^{q+1}{\text{e}}^{lx+gy}+gC{x}^{p+2}{y}^{q+1}{\text{e}}^{lx+gy}+gC{x}^p{y}^{q+3}{\text{e}}^{lx+gy}\end{array}$

1) 令 $A=C=0,B>0$，则

$\begin{array}{c}\frac{\partial \left(DP\right)}{\partial x}+\frac{\partial \left(DQ\right)}{\partial y}=B\alpha y^n+3B{x}^2{y}^n+B{y}^{n+2}+Bnx{y}^{n-1}+B\left(n+1\right)\alpha y^n+B\left(n+1\right)x^2{y}^n+B\left(n+3\right)y^{n+2}\\ =B\left(n+2\right)\alpha y^n+B\left(n+4\right)y^{n+2}+B\left(n+4\right)x^2{y}^n+Bnx{y}^{n-1}\end{array}$

取 $n=-2,B=1,D\left(x,y\right)=y^{-2}$

$\frac{\partial \left(DP\right)}{\partial x}+\frac{\partial \left(DQ\right)}{\partial y}=2+2x^2{y}^{-2}-2x{y}^{-3}$

所以系统(1)在取第二、四象限的点时，

$\frac{\partial \left(DP\right)}{\partial x}+\frac{\partial \left(DQ\right)}{\partial y}>0$

由Dulac判别法 [11] 可知，系统(1)在二、四象限不存在闭轨。

2) 令 $B=C=0,A>0$，则

$\begin{array}{c}\frac{\partial \left(DP\right)}{\partial x}+\frac{\partial \left(DQ\right)}{\partial y}=A\alpha \left(m+1\right)x^m-Am{x}^{m-1}y+A\left(m+3\right)x^{m+2}+A\left(m+1\right)x^m{y}^2+A\alpha x^m+A{x}^{m+2}+3A{x}^m{y}^2\\ =A\alpha \left(m+2\right)x^m+A\left(m+4\right)x^{m+2}+A\left(m+4\right)x^m{y}^2-Am{x}^{m-1}y\end{array}$

取 $m=-2,A=1,D\left(x,y\right)=x^{-2}$

$\frac{\partial \left(DP\right)}{\partial x}+\frac{\partial \left(DQ\right)}{\partial y}=2+2x^{-2}{y}^2+x^{-3}y$

所以系统(1)在取第一、三象限的

$\frac{\partial \left(DP\right)}{\partial x}+\frac{\partial \left(DQ\right)}{\partial y}>0$

同样的由Dulac判别法可知，系统(1)在一、三象限不存在闭轨。

又因为显然 $x=0$ 与 $y=0$ 不是系统的轨线，则当 $\alpha <0$ 时，系统(1)可能存在与x轴，y轴相交的轨线。

下面来证明当 $\alpha <0$ 时，系统(1)极限环的存在性：

为了解系统(1)，做极坐标变换 $x=r\cos \theta$， $y=r\sin \theta$，则有

$\{\begin{array}{l}\frac{\text{d}r}{\text{d}t}=r\left(r^2+\alpha \right)\\ \frac{\text{d}\theta }{\text{d}t}=1\end{array}$

容易看出当 $\alpha <0$ 时，有一个定常解 $r^2=-\alpha$，用直角坐标x与y表示圆周 $x^2+y^2=-\alpha$。

对于其他轨线有两种可能性：

1) 从圆周内部的点出发的轨线：当 $r<\sqrt{-\alpha }$ 时，有 $\frac{\text{d}r}{\text{d}t}<0$，因此，当t增加时，从圆周内部出发的轨线越来越远离该圆周；

2) 从圆周外部的点出发的轨线：当 $r>\sqrt{-\alpha }$ 时，有 $\frac{\text{d}r}{\text{d}t}>0$，因此，当t增加时，从圆周外部的点出发的轨线越来越远离该圆周。

因此该圆周是一个孤立的不稳定的闭轨——不稳定的极限环，如图5所示。

Figure 5. Limit cycle of the system

图5. 系统的极限环

下面讨论极限环的唯一性：

用文献 [7] 中的定理6.1和定理6.8来分别证明系统(1)只存在一个极限环。

由前面可知系统(1)的极坐标形式为：

$\frac{\text{d}r}{\text{d}\theta }=r\left(r^2+\alpha \right)=\phi \left(r,\theta \right)$

取

$\psi \left(r\right)=\frac{1}{r}$ 且 $0<r^{\prime }<\sqrt{-\alpha }<r^{″}$

显然在圆周 $C_{r^{\prime }}:x^2+y^2={r^{\prime }}^2$ 和圆周 $C_{r^{″}}:x^2+y^2={r^{″}}^2$ 之间有：

$\frac{\text{d}\psi }{\text{d}r}=-\frac{1}{r^2}\ne \infty$, $\phi \left(r,\theta \right)\ne \infty$, $\frac{\partial \phi }{\partial r}=3r^2+\alpha \ne \infty$, $\frac{{\text{d}}^2\psi }{\text{d}{r}^2}=\frac{2}{r^3}\ne \infty$

最后有 $\frac{\partial }{\partial r}\left(\phi \frac{\text{d}\psi }{\text{d}r}\right)=\frac{\partial }{\partial r}\left[-\frac{1}{r}\left(r^2+\alpha \right)\right]=\frac{\alpha }{r^2}-1\ne 0$ 因为 $\alpha <0$，故由定理6.1可知系统(1)的极限环唯一。

对于系统(1)

$\{\begin{array}{l}\stackrel{\dot{}}{x}=\alpha x-y+x\left(x^2+y^2\right)=P\left(x,y\right)\\ \stackrel{\dot{}}{y}=x+\alpha y+y\left(x^2+y^2\right)=Q\left(x,y\right)\end{array}$

有

$\begin{array}{l}P\left(x,y\right)Q\left(\lambda x,\lambda y\right)-P\left(\lambda x,\lambda y\right)Q\left(x,y\right)\\ =\left[\alpha x-y+x\left(x^2+y^2\right)\right]\left[\lambda x+\lambda \alpha y+{\lambda }^{3}y\left(x^2+y^2\right)\right]\\ -\left[\lambda \alpha x-\lambda y+{\lambda }^{3}x\left(x^2+y^2\right)\right]\left[x+\alpha y+y\left(x^2+y^2\right)\right]\\ =\lambda \left(1-{\lambda }^2\right){\left(x^2+y^2\right)}^2\end{array}$

对于任意的 $\lambda >1$，恒有 $P\left(x,y\right)Q\left(\lambda x,\lambda y\right)-P\left(\lambda x,\lambda y\right)Q\left(x,y\right)\le 0$ 并且只有取 $\left(0,0\right)$ 点时等号成立，故由定理6.8可知系统(1)的极限环唯一。

证毕。

4. 系统(1)在无穷远处的轨线情况

下面考虑系统(1)在无穷远的情况：

分以下几步进行：

1) 作Poincare’变换

$x=\frac{1}{z}$, $y=\frac{u}{z}$

即：

$u=\frac{y}{x}$, $z=\frac{1}{x}$

把系统(1)变成系统：

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=\frac{u^{2}z+z}{z}\\ \frac{\text{d}z}{\text{d}t}=\frac{u{z}^2-\alpha z^2-u^2-1}{z}\end{array}$ (5)

再令 $\text{d}\tau =\frac{\text{d}t}{z}$，系统(5)变成系统：

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}\tau }=u^{2}z+z=P_1\left(u,z\right)\\ \frac{\text{d}z}{\text{d}\tau }=u{z}^2-\alpha z^2-u^2-1=Q_1\left(u,z\right)\end{array}$ (6)

2) 求出系统(6)在u轴 $\left(z=0\right)$ 上的平衡点，并判断其稳定性：

令 $z=0$，则有

$\{\begin{array}{l}{P}_1\left(u,0\right)=0\\ Q_1\left(u,0\right)=-u^2-1=0\end{array}$

所以方程无实数解，则赤道上无其他奇点。

3) 作Poincare’变换：

$x=\frac{v}{z}$, $y=\frac{1}{z}$

即：

$v=\frac{x}{y}$, $z=\frac{1}{y}$

系统(1)化为：

$\{\begin{array}{l}\frac{\text{d}v}{\text{d}t}=\frac{-v^{2}z-z}{z}\\ \frac{\text{d}z}{\text{d}t}=\frac{-v{z}^2-\alpha z^2-v^2-1}{z}\end{array}$ (7)

再令 $\text{d}\tau =\frac{\text{d}t}{z}$，系统(7)化为：

$\{\begin{array}{l}\frac{\text{d}v}{\text{d}\tau }=-v^{2}z-z=P_2\left(v,z\right)\\ \frac{\text{d}z}{\text{d}\tau }=-v{z}^2-\alpha z^2-v^2-1=Q_2\left(v,z\right)\end{array}$ (8)

求出系统(8)在v轴 $\left(z=0\right)$ 上的平衡点，令 $z=0$，有：

$\{\begin{array}{l}{P}_2\left(v,0\right)=0\\ Q_2\left(v,0\right)=-v^2-1=0\end{array}$

则方程也无实数解。

综上：系统(1)没有无穷远点。

5. 系统的全局相图

由前面的推证，可大致画出系统的全局相图，如图6所示。

Figure 6. Global phase diagram

图6. 全局相图

6. Matlab作图

为验证前面理论部分的准确性，分别取参数α = −1, 0, 1，利用Matlab来做出系统的全局相图，如图7~9所示。

Figure 7. Global phase diagram of system (1) when α = −1

图7. α = −1时系统(1)的全局相图

Figure 8. Global phase diagram of system (1) when α = 0

图8. α = 0时系统(1)的全局相图

Figure 9. Global phase diagram of system (1) when α = 1

图9. α = 1时系统(1)的全局相图

7. 结论

通过对系统的定性分析以及相图中对轨线走向的描述，我们能够清晰地了解到系统(1)的平衡点随参数的改变而发生的分岔情况，当参数 $\alpha <0$ 时，平衡点 $\left(0,0\right)$ 是系统(1)稳定的焦点；当 $\alpha \ge 0$ 时，平衡点 $\left(0,0\right)$ 是系统(1)不稳定的焦点；只有当参数 $\alpha <0$ 时，系统存在一个不稳定的极限环；并且无论参数取多少，系统(1)都没有无穷远点。

但本文只是对这类含一个参数的系统进行研究，研究类型单一，后续还会尝试对含两个参数的系统进行研究，使得研究结果更为完善。

NOTES

^*第一作者。

^#通讯作者。