Phase Plane Method

Motivation

本文讲述的相平面法[Phase Plane Methods]主要围绕对分析二元微分方程组和一元二次微分方程的向量场，在介绍之前，先引入一个二元微分方程组，The Competitive Lotka-Volterra Equations：

$$\frac{dx}{d\overline{t}} = r_xX(1-\frac{x}{x_c}) - F_X(x,y) \\ \frac{dy}{d\overline{t}} = r_yY(1-\frac{y}{y_c}) - F_Y(x,y) \\$$

其中，我们定义$F_X>0, F_Y >0$，这描述了两个同类动物族群的数量演变，我们并不能直接求出这两个方程的具体解，但是我们可以通过相平面法求出(X, Y)的变化趋势与收敛情况。

---

我们首先来看一个简单的方程组，我们假设在整个过程中，所有函数连续并可导：

$$\frac{dx}{dt} = X(x,y) \ \ \ \ \frac{dy}{dt} = Y(x,y)$$

• Integral Path: 对于一组解: $X_k(t), Y_k(t)$ 满足这一微分方程的随着时间的变化曲线。由于微分方程组有无数多解（没有初始值t=0约束的缘故），但每一组解都是一一对应随时间变化的，我们定义某一组解构成的坐标的点随时间运动的轨迹为Integral Path，当然，一组可解的微分方程往往有很多Integral Path
• Phase portrait：Integral Path的集合，在这是一张x-y坐标系图，我们的目标就是画出这样的平面图。

Equilibrium Point 平衡点

平衡点就是（x, y）不再变化的点，但是不代表（x,y）在实际中静止与此，而仅仅代表：$dX/dt=0, dY/dt=0$，等价于：

$$X(x,y) = Y(x,y) = 0$$

Integral paths cannot cross in the phase plane, except at equilibrium points. （我认为这显而易见，但我无法证明）

• Ordinary point：所有非平衡点的点都是Ordinary point
• Direction of the Integral Path: 向量$\vec d = [X(x,y), Y(x,y)]$ X(x0, y0) , Y (x0, y0) 是integral path 在 ordinary point (x0, y0)的切线方向。

Integral Path 的分类

1. Isocline：当x与y的变化率为常数时：即 $y = kx + C$ 时，或者写作：

   $$\frac{dy}{dx} = \frac{dy}{dt}/\frac{dx}{dt} = \frac{Y(x,y)}{X(x,y)} = Constant$$

2. Nullcline：Y(x, y) 或者 X(x, y) 为0的 Integral Path，Nullcline贯穿所有平衡点

对于Nullcline上的点，可以判断该点的移动趋势，实际上是分析x或者y的变化的正负号：

• 当X=0时：若Y大于0，则代表Y的变化趋势是+，不断增大向上移动

• 当Y=0时：若X大于0，则代表x变大，向右变化

  

  对Nullcline与变化趋势的解释

实际上，对于变化趋势的判断用到了对于特定点周围的线性近似过程，我们也可以使用这个方法来进一步判断平衡点的性质。

平衡点的性质判断

下列式子是对$(x_{eq}, y_{eq})$的线性近似：

$$\frac{d\hat{x}}{dt} = X(x_{eq} + \hat{x}, y_{eq} + \hat{y}) \approx X(x_{eq}, y_{eq}) + \hat{x} \frac{\partial X}{\partial x}(x_{eq}, y_{eq}) + \hat{y} \frac{\partial X}{\partial y}(x_{eq}, y_{eq}) \\ \frac{d\hat{y}}{dt} = Y(x_{eq} + \hat{x}, y_{eq} + \hat{y}) \approx Y(x_{eq}, y_{eq}) + \hat{x} \frac{\partial Y}{\partial x}(x_{eq}, y_{eq}) + \hat{y} \frac{\partial Y}{\partial y}(x_{eq}, y_{eq})$$

由于平衡点的性质，在小范围内，等价于：

$$\frac{d\hat{x}}{dt} = \hat{x} \frac{\partial X}{\partial x}(x_{eq}, y_{eq}) + \hat{y} \frac{\partial X}{\partial y}(x_{eq}, y_{eq}) \\ \frac{d\hat{y}}{dt} = \hat{x} \frac{\partial Y}{\partial x}(x_{eq}, y_{eq}) + \hat{y} \frac{\partial Y}{\partial y}(x_{eq}, y_{eq})$$

写成矩阵的形式，事实上，我们只关心其Jacobian matrix

$$\frac{d\vec y}{dt} = \frac{d}{dt} \begin{pmatrix} \hat{x} \\ \hat{y} \end{pmatrix} = \begin{pmatrix} \frac{\partial X}{\partial x}(x_{eq}, y_{eq}) & \frac{\partial X}{\partial y}(x_{eq}, y_{eq}) \\ \frac{\partial Y}{\partial x}(x_{eq}, y_{eq}) & \frac{\partial Y}{\partial y}(x_{eq}, y_{eq}) \end{pmatrix} \begin{pmatrix} \hat{x} \\ \hat{y} \end{pmatrix}$$

为了方便讨论平衡点的性，我们着重分析 Nullcline，当X为0时，在足够接近平衡点的时候(x,y)可以看做是线性变化，也就意味着对于给定的向量坐标 $\vec v$ ，其变化趋势的方向与自身指向方向一致，即 $\vec x$ 可以看作是Jacobian的特征向量，所以问题的分析变成了对特征值的分析。我们对其求特征向量，有：

$$\textbf{y}_{+} = c_+\textbf{v}_+e^{\lambda_+t} \\ \textbf{y}_{-} = c_-\textbf{v}_-e^{\lambda_-t}$$

Case 1: 特征根为实数

这时候，向量y就有两个方向，约束平衡点周围的点的变化趋势：

• 当 $\vec y$ 对应的特征值大于0时，$|\textbf{y}| \rightarrow \infty \text{ as } t \rightarrow \infty$, 我们说在这一方向上时发散的，或者说Unstable
• 当 $\vec y$ 对应的特征值大于0时，$|\textbf{y}| \rightarrow 0 \text{ as } t \rightarrow \infty$, 我们说在这一方向上时收敛于平衡点，或者说Stable

进一步有如下分类：

1. $\lambda_+ > 0 > \lambda_-$ : 两个方向不一致，为Saddle
2. $\lambda_+ > \lambda_- > 0$：均为Unstable → 平衡点Unstable
3. $0 > \lambda_+ > \lambda_-$: 均为Stable → 平衡点Stable

Case 2: 特征根为虚数

由于特征值共轭，所以在实数部分两者一致，进一步模长随着时间变化的趋势一致，即同时趋向无穷或者0随着时间的变大，在图像上，表达为螺旋状：

左图解趋向于无穷，代表实部为正数，Unstable Spiral；右侧趋向于0，为负数, Stable Spiral

当实部为0时，我们认为这是 linear centre，(x, y) 随着时间变化为周期运动，称之为：conserved quantity

An equilibrium point whose linearised system has no eigenvalues with zero real part is said to be hyperbolic.

我们称没有实部为0的平衡点为：Hyperbolic equilibrium point，除了Linear Center之外都是Hyperbolic Equilibrium Point

An equilibrium point whose linearised system has at least one eigenvalue with zero real part is said to be Non-hyperbolic.

Example: 对Competitive Lotka Volterra Equations分析相平面

在这个例子中，我们进一步化简方程组，我们假设 $F_X 和 F_y$ 与x, y呈线性关系：

$$\frac{dx}{dt} = x(1-x) - \alpha xy \\ \frac{dy}{dt} = r[y(1-y)-\beta xy]$$

Step 1: 求解Jacobian matrix 和平衡点

根据定义有：

$$x(1-x) - \alpha xy =0 \\ r\left\{y(1-y)-\beta xy\right\}=0$$

其中，(0,0),(1,0) and (0,1) 显而易见是方程的平衡点，从现实意义上讲，这都意味着一个种族的人口为0，我们还能发现:

$$x≡x_c = \frac{1-\alpha}{1 - \alpha \beta} \\ y≡y_c = \frac{1-\beta}{1 - \alpha \beta}$$

我们称这组解为coexistence equilibrium state，即两个种族的共生状态。但是这个解并不是一直成立，因为我们定义x, y都是大于等于0的。所以𝛼与𝛽应同时在（0，1）之间或者同时在(1, +∞), 才存在这个平衡点。

Equilibrium Point and Jacobian Matrix

Step 2: Classify the equilibrium points

Invariant set：A set of points $S \in R^2$ is an invariant set if $(x,y) \in S$ at some time t implies that $(x,y) \in S \text{ for } ∀t.$对于这个模型来说，x, y轴就是Invariant Set

Step 3: 相平面绘制

在最后的Case中，我们介绍一个概念：Basin of attraction：

在相平面中，”吸引盆”（Basin of Attraction）是指一个动力系统中，对于一个给定的稳定平衡点或稳定极限环，所有从中起始最终趋向于该平衡点或极限环的初始条件的集合。简单来说，吸引盆是相空间中的一个区域，如果系统的状态从该区域内的任意点开始演化，随着时间的推移，它将逐渐靠近并最终达到稳定状态（如一个吸引子，可以是一个点、闭合轨道或是更复杂的集合）。

Second Order ODE

对于二阶微分方程，我们也可以使用相平面法去探索物体的变化规律，单摆是一个简单的例子：

$$F = -kx$$

我们已知单摆在小角度下成立，将其写成ODE的形式，并将参数定义为ω,有：

$$\frac{d^2\theta}{dt^2} ≡ \frac{d^2x}{dt^2} = -\omega^2x = -\omega^2sin\theta$$

我们令 $q = \theta_t$, q可以理解为一个速度变量，构造二元方程：

$$\frac{d\theta}{dt} = q \\ \frac{dq}{dt} = -\omega^2 sin\theta$$

我们可以发现有平衡点(0, 0)，由于角度的周期性，角度加任意倍数的π都为平衡点，注意，（0， 0）是一个Linear Center，容易理解的，角度加360度不影响性质，所以任意360度的倍数都是Linear Center, 但是对于角度加180度，这并不能构成一个单摆的稳定状态，它将继续运动直到收敛于单摆运动，是saddle的。

Conserved System

在这个场景中，我们认为是速度和位移两个元素构成了一个系统并探究两者的运动关系，在这个系统中，始终遵守：

$$\frac{d\theta}{dq} = \frac{d\theta}{dt} / \frac{dq}{dt} = -\frac{q}{\omega^2sin\theta}$$

那么通过求解微分方程，我们得到：

$$\frac{1}{2}q^2 - \omega^2cos\theta = Constant$$

我们称这个系统是保守系因为两者存在一个定值的约束关系，我们定义这个约束关系为E，称之为Conserved quantity。

根据这个约束条件，我们可以绘制不同定值下两者的变化关系。

• An integral path that connects one saddle point to another is called a heteroclinic connection or heteroclinic orbit串联所有的Saddle Point

• Homoclinic path ：过所有的平衡点的线

  

根据定义，我们可以知道红色轨迹是heteroclinic路径，因为这条路径经过所有的Saddle Point

当角度初始为小角度时，轨迹围绕Linear Center平衡点，在Limit Cycle中运动，即有周期性的简谐运动。当角度过大，或者q过大，则会做圆周运动，如蓝线所示。

Limit Cycle

在保守系统中，我们定义：极限环（Limit Cycle）主要用于描述动态系统的行为。极限环是一个封闭的轨道，这意味着系统的状态（通常表现为某些变量的集合）随时间演化后，将返回到起始位置，并重复这一周期性的行为。接下来我们将讨论如何判断一个封闭空间有没有Limit Cycle

Method 1: Poincaré index

对于一个封闭曲线，我们为这个曲线计算一个参数，叫Poincare Index，有以下准则：

• 如果这个曲线围城的区域中不含有任何平衡点，则其Index=0
• 若$L = L_1 +L_2$,则 $Index(R) = Index(R_1) + Index(R_2)$
• 若这个曲线是Limit Cycle，Index = 1
• 曲线每包含一个：
  • Saddle Point: -1
  • Center: +1 包括 Nodes， spirals

Theorem：若曲线围的所有平衡点的Index和不为1，曲线围的区域不存在Limit Center

Method 2: Bendixson’s negative criterion

对于一个简单的封闭区域R，若对于R中的点(x,y):

$$D = \frac{\partial X}{\partial x} + \frac{\partial Y}{\partial y}$$

若D在遍历所有点时符号无改变，则不存在Limit Cycle在R中。但这个式子并不是非常有用，进一步，可以结合Dulac’s extension

构造在R中连续函数 $\rho(x, y)$ ，若

$$D = \frac{\partial \rho X}{\partial x} + \frac{\partial \rho Y}{\partial y}$$

在R中所有点仍不变号，则无Limit Cycle

我们可以用这个方法检测Lotka-Volterra system

$$\frac{dx}{dt} = x(1-x) - \alpha xy \\ \frac{dy}{dt} = r[y(1-y)-\beta xy]$$

在不结合Dulac’s extension, $X_x = 1 - 2x - \alpha y, Y_y = r - 2r - r\beta x$

$$D= 1 − 2x − αy + r(1 − 2y) − rβx.$$

这看起来并并不直观，但是当结合了Dulac’s extension，我们定义 $\rho(x,y) = {1\over xy}$:

$$D = {\partial \over \partial x} {1\over y}[(1-x) - \alpha]+ {\partial \over \partial x} {r\over x}[(1-y) - \beta] = -{1\over y } - {r \over x}$$

由于(x,y)总在第一象限，所以D恒小于0，则在第一象限内不存在Limit Cycle.

Poincaré-Bendixson theorem

把这两个人的名字和在一起，我们有这一定理：

Theorem：在二维平面中，若一个区域不存在平衡点，而Integral Path从所有边缘方向进入这个区域，那么这个区域一定存在Limit Cycle。

在这个定理的正式定义前，需要解释几个概念：

• Positively-invariant set：t=0 随着时间变大，Solution一直在这个区域内，如稳定螺旋
• Negatively-invariant set：t=0 随着时间变小（小于0），Solution一直在这个区域内，如不稳定螺旋
• 若一个Integral Path不仅是Positively-invariant set还是Negatively-invariant set，那么这就是Limit Cycle

Poincaré-Bendixson: If there exists a bounded, invariant region I in the phase plane, and I contains no equilibrium points, then I contains at least one limit cycle.