Logo

离散动力系统

1. 理论基础与应用背景

特征值和特征向量是分析动力系统长期行为最为经典且广泛的工具,它们在描述系统的稳定性、周期性和长期趋势方面具有重要意义。本节我们主要借助差分方程来分析离散动力系统,探讨系统在离散时间点上的状态演化。我们在前面介绍过特征值在离散动力系统中的应用,下面我们将其模型进行一般化处理。

2. 离散动力系统的数学模型

我们先假设矩阵 A ~\mathbf{A}~是可对角化的,并且具有 n ~n~个线性无关的特征向量 {v1,v2,,vn} ~\{\mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_n\}~,对应的特征值为 {λ1,λ2,,λn} ~\{\lambda_1, \lambda_2, \ldots, \lambda_n\}~。为方便分析,特征向量按特征值的模大小排序,使得 λ1λ2λn ~|\lambda_1| \geq |\lambda_2| \geq \cdots \geq |\lambda_n|~。基于此,对模型的分析和推广步骤如下:

任何初始状态向量 x0\mathbf{x}_0 都可以用特征向量的线性组合表示:
x0=c1v1+c2v2++cnvn\mathbf{x}_0 = c_1\mathbf{v}_1 + c_2\mathbf{v}_2 + \cdots + c_n\mathbf{v}_n
其中,系数 cic_i 是唯一确定的。

通过差分方程 xk+1=Axk\mathbf{x}_{k+1} = A\mathbf{x}_k,可以推导出系统在任意时间步的状态。对于第一步 k=1k = 1,状态可以写为:
x1=Ax0=c1Av1+c2Av2++cnAvn\mathbf{x}_1 = A\mathbf{x}_0 = c_1A\mathbf{v}_1 + c_2A\mathbf{v}_2 + \cdots + c_nA\mathbf{v}_n
由于 vi\mathbf{v}_iAA 的特征向量,且满足 Avi=λiviA\mathbf{v}_i = \lambda_i\mathbf{v}_i,我们可以得到:
x1=c1λ1v1+c2λ2v2++cnλnvn\mathbf{x}_1 = c_1\lambda_1\mathbf{v}_1 + c_2\lambda_2\mathbf{v}_2 + \cdots + c_n\lambda_n\mathbf{v}_n

按照同样的逻辑,对于第二步 k=2k = 2,状态进一步演化为:
x2=Ax1=c1λ1Av1+c2λ2Av2++cnλnAvn\mathbf{x}_2 = A\mathbf{x}_1 = c_1\lambda_1A\mathbf{v}_1 + c_2\lambda_2A\mathbf{v}_2 + \cdots + c_n\lambda_nA\mathbf{v}_n
利用特征向量的性质,状态可以写成:
x2=c1(λ1)2v1+c2(λ2)2v2++cn(λn)2vn\mathbf{x}_2 = c_1(\lambda_1)^2\mathbf{v}_1 + c_2(\lambda_2)^2\mathbf{v}_2 + \cdots + c_n(\lambda_n)^2\mathbf{v}_n
从这里可以看出,系统的状态在每一步都会受到特征值的幂次放大的影响。

以此类推,系统在任意状态可以表示为:
xk=c1(λ1)kv1++cn(λn)kvn(k=0,1,2,)(1)\colorbox{#F0F8FF}{$\mathbf{x}_k = c_1(\lambda_1)^k\mathbf{v}_1 + \cdots + c_n(\lambda_n)^k\mathbf{v}_n \quad (k = 0, 1, 2, \ldots)$} \tag{1}

由递推公式可知,动力系统的长期行为主要由模最大的特征值λ1\lambda_1 决定。随着kk \to \infty,那些模较小的特征值对应的项ci(λi)kvic_i(\lambda_i)^k\mathbf{v}_i会快速衰减,仅剩下主特征值 λ1\lambda_1及其对应的特征向量 v1\mathbf{v}_1对系统的长期状态产生显著影响。

3. 捕食者与猎物模型的案例分析

开通会员解锁全部动画

捕食者-猎物模型用于研究猫头鹰(捕食者)和老鼠(猎物)种群的动态变化,描述捕食行为和种群繁殖对数量的影响。作为一个简单的离散动力系统,该模型可通过矩阵的特征值和特征向量进行分析。以下是详细的分析过程:

  • 不同时间 k ~k~时刻的种群状态用一个向量表示:
    xk=[OkRk]\mathbf{x}_k = \begin{bmatrix} O_k \\ R_k \end{bmatrix}
    其中:
    •  Ok ~O_k~:第 k ~k~月猫头鹰的数量。

    •  Rk ~R_k~:第 k ~k~月老鼠的数量(以千为单位)。

  • 种群随时间变化由以下公式给出:
    Ok+1=0.5Ok+0.4RkRk+1=pOk+1.1Rk\begin{aligned} O_{k+1} &= 0.5O_k + 0.4R_k \\[1ex] R_{k+1} &= -p \cdot O_k + 1.1R_k \end{aligned}
    公式的含义:
    • 0.5Ok 0.5O_k~:如果没有老鼠( Rk=0 ~R_k=0~),每个月只有 50% ~50\%~的猫头鹰能存活下来。

    • 0.4Rk 0.4R_k~:老鼠的数量和猫头鹰正相关,老鼠越多猫头鹰越多。

    • 1.1Rk 1.1R_k~:在没有猫头鹰的情况( Ok=0 ~O_k=0~),老鼠每个月的自然增长率是 10~10%~

    • pRk -p\cdot R_k~:猫头鹰的捕食会造成老鼠数量减少, p ~p~的值一般由统计数据给出。

上述公式可以用矩阵形式表示为:
xk+1=Axk,A=[0.50.4p1.1]\mathbf{x_{k+1}} = \mathbf{A} \mathbf{x_k}, \quad \mathbf{A} = \begin{bmatrix} 0.5 & 0.4 \\ -p & 1.1 \end{bmatrix}

 p=0.104 ~p=0.104~时,矩阵 A ~\mathbf{A}~的特征值为:
λ1=1.02,λ2=0.58\lambda_1=1.02,\quad \lambda_2 = 0.58
对应的特征向量为:
v1=[1013],v2=[51]\mathbf{v_1} = \begin{bmatrix} 10 \\ 13 \end{bmatrix}, \quad \mathbf{v_2} = \begin{bmatrix} 5 \\ 1 \end{bmatrix}
这些特征值和特征向量描述了系统的长期行为。

初始状态 x0 ~\mathbf{x}_0~可以表示为特征向量的线性组合:
x0=c1v1+c2v2\mathbf{x}_0=c_1\mathbf{v}_1 + c_2\mathbf{v}_2
随时间的演化公式为:
xk=c1λ1kv1+c2λ2kv2\mathbf{x_k} = c_1 \lambda_1^k \mathbf{v_1} + c_2 \lambda_2^k \mathbf{v_2}
其中,
  •  λ1k ~\lambda_1^k~:它主导系统的长期行为(λ1>λ2\lambda_1 > \lambda_2)。

  •  λ2k ~\lambda_2^k~:随时间迅速衰减(因为0.58k0, 当 k0.58^k \to 0,~ \small{当}~ k \to \infty)。

 k ~k~足够大时, λ2k ~\lambda_2^k~项可以忽略,因此:
xkc1λ1kv1=c1(1.02)k[1013]\mathbf{x}_k \approx c_1 \lambda_1^k \mathbf{v}_1 = c_1 (1.02)^k \begin{bmatrix} 10 \\ 13 \end{bmatrix}
这说明:
  • 种群每月以大约 2% ~2\%~的速度增长

  • 无论初始状态如何,猫头鹰与老鼠的数量比接近于 10:13 ~10:13~(老鼠的数量单位是:千)。

捕食者-猎物动力系统的长期行为由主特征值 λ1 ~\lambda_1~和对应的特征向量 v1 ~\mathbf{v}_1~决定。对于足够大的 k ~k~,系统满足以下近似:
xk+1λ1xk,xkc1(λ1)kv1\mathbf{x_{k+1}} \approx \lambda_1 \mathbf{x_k}, \quad \mathbf{x_k} \approx c_1 (\lambda_1)^k \mathbf{v_1}
这表明系统的增长率由 λ1 ~\lambda_1~确定,状态向量的方向趋于 v1 ~\mathbf{v}_1~

4. 解的几何意义(对角矩阵)

我们还是先将离散动力系统的研究范围限制在对角矩阵的情况。当 A ~\mathbf{A}~ 2×2 ~2\times 2~的矩阵时,我们可以通过几何形式来描述系统的演化过程。方程:
xk+1=Axk\mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k
可以看作是描述初始点 x0 ~\mathbf{x}_0~在二维空间 R2 ~\mathbb{R^2}~中,通过映射xAx\mathbf{x}\mapsto \mathbf{A}\mathbf{x}的不断变换后发生的变化。连接点x0,x1,\mathbf{x}_0,\mathbf{x}_1,\dots的图形称为动力系统的轨迹trajectory\textbf{trajectory} )。根据特征值的情况不同,离散动力系统的轨迹主要可分为吸引点(Attractor\textbf{Attractor})、排斥点(Repeller\textbf{Repeller})、鞍点(Saddle Point\textbf{Saddle Point})。下面通过示例来逐个介绍。

3.1 吸引点

吸引点发生在所有特征值的绝对值都小于 1 ~\mathbf{1}~的情况。此时,系统的状态向量 xk ~\mathbf{x}_k~随时间趋向于原点,轨迹逐渐衰弱。状态分量沿特征值绝对值较大的方向收敛更快,而特征值绝对值较小的方向对长期轨迹的影响逐渐减弱。

设系数矩阵为:
A=[0.8000.64]\mathbf{A} = \begin{bmatrix} 0.8 & 0 \\ 0 & 0.64 \end{bmatrix}
这是一个对角矩阵,其特征值为 λ1=0.8,λ2=0.64 ~\lambda_1=0.8,\,\lambda_2=0.64~。对应的特征向量为:
v1=[10],v2=[01]\mathbf{v}_1 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \quad \mathbf{v}_2 = \begin{bmatrix} 0 \\ 1 \end{bmatrix}
任意初始状态 x0\mathbf{x}_0可以表示为特征向量的线性组合:
x0=c1v1+c2v2\mathbf{x}_0 = c_1 \mathbf{v}_1 + c_2 \mathbf{v}_2
在经过 kk 次迭代后,系统的状态向量变为:
xk=c1(0.8)kv1+c2(0.64)kv2\mathbf{x}_k = c_1 (0.8)^k \mathbf{v}_1 + c_2 (0.64)^k \mathbf{v}_2
随着 kk \to \infty,两个特征值的幂次 (0.8)k(0.8)^k(0.64)k(0.64)^k 均趋于 0 ~0~,因此xk0\mathbf{x}_k \to \mathbf{0}

开通会员解锁全部动画

3.2 排斥点

排斥点发生在所有特征值的绝对值均大于 1 ~\mathbf{1}~的情况下。此时,系统的状态向量xk\mathbf{x}_k随时间远离原点,轨迹逐渐发散。其中,沿着特征值绝对值较大的方向扩散更快。

设系数矩阵为:
A=[1.44001.2]\mathbf{A} = \begin{bmatrix} 1.44 & 0 \\ 0 & 1.2 \end{bmatrix}
这也是一个对角矩阵,其特征值为 λ1=1.44\lambda_1 = 1.44λ2=1.2\lambda_2 = 1.2。对应的特征向量分别为:
v1=[10],v2=[01]\mathbf{v}_1 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \quad \mathbf{v}_2 = \begin{bmatrix} 0 \\ 1 \end{bmatrix}
任意初始状态 x0\mathbf{x}_0可以表示为特征向量的线性组合:
x0=c1v1+c2v2\mathbf{x}_0 = c_1 \mathbf{v}_1 + c_2 \mathbf{v}_2
经过 kk 次迭代后,系统的状态向量变为:
xk=c1(1.44)kv1+c2(1.2)kv2\mathbf{x}_k = c_1 (1.44)^k \mathbf{v}_1 + c_2 (1.2)^k \mathbf{v}_2
由于 (1.44)k(1.44)^k(1.2)k(1.2)^k随着 kk \to \infty 快速增长,轨迹逐渐远离原点。

开通会员解锁全部动画

3.3 鞍点

鞍点发生在一个特征值的绝对值大于 1 ~\mathbf{1}~另一个特征值的绝对值小于 1 ~\mathbf{1}~的情况下。此时,系统的状态向量 xk ~\mathbf{x}_k~表现为:沿着特征值绝对值小于 1 ~1~的特征向量方向,轨迹逐渐靠近原点(吸引方向)、沿着特征值绝对值大于 1 ~1~的特征向量方向,轨迹远离原点(排斥方向)。

设系数矩阵为:
A=[2.0000.5]\mathbf{A} = \begin{bmatrix} 2.0 & 0 \\ 0 & 0.5 \end{bmatrix}
这是一个对角矩阵,其特征值为 λ1=2.0\lambda_1 = 2.0λ2=0.5\lambda_2 = 0.5。对应的特征向量分别为:
v1=[10],v2=[01]\mathbf{v}_1 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \quad \mathbf{v}_2 = \begin{bmatrix} 0 \\ 1 \end{bmatrix}
任意初始状态 x0\mathbf{x}_0可以表示为特征向量的线性组合:
x0=c1v1+c2v2\mathbf{x}_0 = c_1 \mathbf{v}_1 + c_2 \mathbf{v}_2
经过 kk 次迭代后,系统的状态向量变为:
xk=c1(2.0)kv1+c2(0.5)kv2\mathbf{x}_k = c_1 (2.0)^k \mathbf{v}_1 + c_2 (0.5)^k \mathbf{v}_2
动态行为的解读
  1. 吸引方向:

    • 对应特征值 λ2=0.5 ~\lambda_2=0.5~,该方向上的分量(0.5)k(0.5)^k会随着kk \rightarrow \infty指数衰减,轨迹逐渐靠近原点。

    • 吸引速度由特征值的绝对值决定,绝对值越小,收敛速度越快。

  2. 排斥方向:

    • 对应特征值 λ1=2.0 ~\lambda_1=2.0~,该方向上的分量(2.0)k(2.0)^k会随着kk \rightarrow \infty指数增长,轨迹逐渐远离原点。

    • 轨迹扩散速度由特征值的大小决定,绝对值越大,扩散越快。

  3. 混合行为:

    • 在初期,轨迹可能受到两个特征值的共同作用,呈现弯曲的轨迹。

    • 长期来看,轨迹表现为沿排斥方向的快速发散,但在吸引方向上逐渐靠近原点。

开通会员解锁全部动画

5. 非对角矩阵的变量替换

接下来我们将讨论的对象从前面的对角矩阵扩展至更一般的可对角化非对角矩阵。对于给定的动力系统:
xk+1=Axk\mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k
其中 A ~\mathbf{A}~是一个 n×n ~n\times n~矩阵。当矩阵 A ~\mathbf{A}~是非对角矩阵时,我们的处理方法是首先对其进行对角化(即A=PDP1\mathbf{A}=\mathbf{P}\mathbf{D}\mathbf{P}^{-1})。通过这一过程,将动力系统中的变量 xk ~\mathbf{x}_k~转换为新的变量 yk ~\mathbf{y}_k~,从而实现系统的解耦。以下是变量替换的过程:

给定动力系统:
xk+1=Axk\mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k
其中 A ~\mathbf{A}~是一个 n×n ~n\times n~矩阵。

为了简化分析,引入新的变量 yk\mathbf{y}_k,定义为:
yk=P1xk\mathbf{y}_k = \mathbf{P}^{-1} \mathbf{x}_k
其中 P\mathbf{P} 是由矩阵A\mathbf{A} 的特征向量组成的矩阵:
P=[v1v2vn]\mathbf{P} = \begin{bmatrix} \mathbf{v}_1 & \mathbf{v}_2 & \cdots & \mathbf{v}_n \end{bmatrix}
通过这一替换,xk\mathbf{x}_k 可以重新表示为:
xk=Pyk\mathbf{x}_k = \mathbf{P} \mathbf{y}_k

将变量替换代入原始动力系统:
xk+1=Axk\mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_k
结合变量替换xk=Pyk\mathbf{x}_k = \mathbf{P} \mathbf{y}_k,得到:
Pyk+1=APyk\mathbf{P} \mathbf{y}_{k+1} = \mathbf{A} \mathbf{P} \mathbf{y}_k
由于A=PDP1\mathbf{A} = \mathbf{P} \mathbf{D} \mathbf{P}^{-1}(其中 D\mathbf{D} 是由A\mathbf{A} 的特征值构成的对角矩阵),代入后有:
Pyk+1=PDyk\mathbf{P} \mathbf{y}_{k+1} = \mathbf{P} \mathbf{D} \mathbf{y}_k
左右同乘以 P1\mathbf{P}^{-1},得到:
yk+1=Dyk\mathbf{y}_{k+1} = \mathbf{D} \mathbf{y}_k

通过变量替换,原始系统的递推公式被解耦为:
[y1(k+1)y2(k+1)yn(k+1)]=[λ1000λ2000λn][y1(k)y2(k)yn(k)]\begin{bmatrix} y_1(k+1) \\ y_2(k+1) \\ \vdots \\ y_n(k+1) \end{bmatrix} = \begin{bmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_n \end{bmatrix} \begin{bmatrix} y_1(k) \\ y_2(k) \\ \vdots \\ y_n(k) \end{bmatrix}
其中,λ1,λ2,,λn\lambda_1, \lambda_2, \dots, \lambda_n 是矩阵A\mathbf{A} 的特征值。

解耦后的系统形式表明每个分量的变化是独立的。例如:
y1(k+1)=λ1y1(k),y2(k+1)=λ2y2(k),,yn(k+1)=λnyn(k)y_1(k+1) = \lambda_1 \cdot y_1(k), \quad y_2(k+1) = \lambda_2 \cdot y_2(k), \quad \dots, \quad y_n(k+1) = \lambda_n \cdot y_n(k)
这说明系统的每个分量仅由对应的特征值 λi\lambda_i 决定。
下面看一个具体的示例,对于离散动态系xk+1=xk\mathbf{x}_{k+1} = \mathbf{x}_k,其中:
A=[1.250.750.751.25]\mathbf{A} = \begin{bmatrix} 1.25 & -0.75 \\ -0.75 & 1.25 \end{bmatrix}
 A ~\mathbf{A}~的特征值为 λ1=2,λ2=0.5 ~\lambda_1=2,\,\lambda_2=0.5~,对应的特征向量为:
v1=[11],v2=[11]\mathbf{v}_1 = \begin{bmatrix} 1 \\ -1 \end{bmatrix}, \quad \mathbf{v}_2 = \begin{bmatrix} 1 \\ 1 \end{bmatrix}
由于λ1>1 且 λ2<1|\lambda_1| > 1~\text{且}~|\lambda_2| < 1,因此动力系统的原点是一个鞍点,经过 k ~k~次迭代后,系统的状态向量变为:
xk=c1(2.0)kv1+c2(0.5)kv2\mathbf{x}_k = c_1 (2.0)^k \mathbf{v}_1 + c_2 (0.5)^k \mathbf{v}_2
该公式与前面示例3中的形式类似,但有关键区别:示例3使用的是标准基 E={e1,e2} ~\mathcal{E}=\{\mathbf{e}_1,\mathbf{e}_2\}~换成了新基向量B={v1,v2}\mathcal{B} = \{\mathbf{v}_1,\mathbf{v}_2\}, 其轨迹如下:

开通会员解锁全部动画

6. 复特征值的情况

当一个实数 2×2 ~2\times 2~矩阵 A ~\mathbf{A}~具有复特征值时(即不可对角化的情况),动力系统的行为取决于复特征值的模长(r=a2+b2|r| = \sqrt{a^2 + b^2})。如果复特征值的模长大于 1 ~1~,系统轨迹会围绕原点螺旋向外,原点充当排斥点;如果模长小于 1 ~1~,轨迹则会螺旋向原点收敛,原点充当吸引点。复特征值的这种性质使得系统的动态表现可以是螺旋状的收敛或发散。例如下面的矩阵:
A=[0.80.50.11.0]\mathbf{A} = \begin{bmatrix}0.8 & 0.5 \\ -0.1 & 1.0\end{bmatrix}
A \mathbf{A}~的特征值为λ=0.9±0.2i\lambda=0.9\pm 0.2i,对应的特征向量为:[12i1]T\begin{bmatrix} 1\mp 2i & 1\end{bmatrix}^T,复特征值的模长r=0.92+0.220.92|r| = \sqrt{0.9^2 + 0.2^2} \approx 0.92,那么由于这个模长小于 1 ~1~,因此系统的轨迹应该螺旋向原点收敛,即原点是一个吸引点。下面动画演示示了三个不同初始向量下的轨迹:

开通会员解锁全部动画