本文是《微分方程完全入门》系列的第四篇,覆盖一阶常微分方程组的相平面分析、平衡点分类和李雅普诺夫稳定性。阅读本文前需掌握特征根的概念——我们将把特征根从”单个方程”推广到”方程组”。


第四篇:微分方程组与相平面——从一维到多维

一、概念引入——为什么需要方程组

一个四旋翼无人机的动力学需要 12 个一阶 ODE 来描述——位置(3)+速度(3)+姿态(3)+角速度(3)。这不是单个方程能处理的。

一阶 ODE 方程组的一般形式:

{dx1dt=f1(x1,x2,,xn)dx2dt=f2(x1,x2,,xn)dxndt=fn(x1,x2,,xn) \begin{cases} \frac{dx_1}{dt} = f_1(x_1, x_2, \dots, x_n) \\ \frac{dx_2}{dt} = f_2(x_1, x_2, \dots, x_n) \\ \vdots \\ \frac{dx_n}{dt} = f_n(x_1, x_2, \dots, x_n) \end{cases}

写成向量形式:

x˙=f(x),xRn \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}), \quad \mathbf{x} \in \mathbb{R}^n

二、线性化——在平衡点附近近似

2.1 平衡点

如果 f(x)=0\mathbf{f}(\mathbf{x}^*) = \mathbf{0},则 x\mathbf{x}^* 是一个平衡点

在平衡点附近,用泰勒展开取一阶近似:

x˙f(x)+fxx(xx)=0+A(xx) \dot{\mathbf{x}} \approx \mathbf{f}(\mathbf{x}^*) + \frac{\partial\mathbf{f}}{\partial\mathbf{x}}\bigg|_{\mathbf{x}^*}(\mathbf{x} - \mathbf{x}^*) = \mathbf{0} + \mathbf{A}(\mathbf{x} - \mathbf{x}^*)

其中 A\mathbf{A}雅可比矩阵

A=[f1x1f1x2f2x1f2x2] \mathbf{A} = \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots \\ \vdots & \vdots & \ddots \end{bmatrix}

线性化系统y˙=Ay\dot{\mathbf{y}} = \mathbf{A}\mathbf{y},其中 y=xx\mathbf{y} = \mathbf{x} - \mathbf{x}^*

2.2 手算例子:单摆

单摆方程:

θ¨+gLsinθ=0 \ddot{\theta} + \frac{g}{L}\sin\theta = 0

写成一阶系统:x1=θx_1 = \theta, x2=θ˙x_2 = \dot{\theta}

{x˙1=x2x˙2=gLsinx1 \begin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = -\frac{g}{L}\sin x_1 \end{cases}

平衡点:x1=0x_1 = 0π\pi(即 θ=0\theta=0θ=π\theta=\pi

雅可比矩阵:

A=[01gLcosx10] \mathbf{A} = \begin{bmatrix} 0 & 1 \\ -\frac{g}{L}\cos x_1 & 0 \end{bmatrix}

θ=0\theta=0(最低点):

A=[01g/L0] \mathbf{A} = \begin{bmatrix} 0 & 1 \\ -g/L & 0 \end{bmatrix}

特征值:λ=±ig/L\lambda = \pm i\sqrt{g/L}——中心点(中性稳定)

θ=π\theta=\pi(倒立摆,最高点):

A=[01g/L0] \mathbf{A} = \begin{bmatrix} 0 & 1 \\ g/L & 0 \end{bmatrix}

特征值:λ=±g/L\lambda = \pm \sqrt{g/L}——鞍点(不稳定!)


三、相平面分析——二维系统的几何

3.1 二维线性系统的四种平衡点

对于 2×22\times 2 系统 y˙=Ay\dot{\mathbf{y}} = \mathbf{A}\mathbf{y},平衡点的类型由特征值决定:

特征值 平衡点类型 图示
两实负 λ1<λ2<0\lambda_1 < \lambda_2 < 0 稳定结点 所有轨道汇入原点
两实正 0<λ1<λ20 < \lambda_1 < \lambda_2 不稳定结点 所有轨道从原点出发
一正一负 λ1<0<λ2\lambda_1 < 0 < \lambda_2 鞍点 沿一个方向稳定,另一方向不稳定
共轭复根 α<0\alpha < 0 稳定焦点 螺旋衰减进入原点
共轭复根 α=0\alpha = 0 中心点 椭圆轨道,不衰减也不发散
共轭复根 α>0\alpha > 0 不稳定焦点 螺旋发散出原点

3.2 与二阶ODE阻尼比的对应

在第二篇中我们学的二阶ODE:

y¨+2ζωny˙+ωn2y=0 \ddot{y} + 2\zeta\omega_n\dot{y} + \omega_n^2 y = 0

写成方程组:

{x˙1=x2x˙2=ωn2x12ζωnx2 \begin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = -\omega_n^2 x_1 - 2\zeta\omega_n x_2 \end{cases}

特征方程:λ2+2ζωnλ+ωn2=0\lambda^2 + 2\zeta\omega_n\lambda + \omega_n^2 = 0

特征值:λ=ζωn±ωnζ21\lambda = -\zeta\omega_n \pm \omega_n\sqrt{\zeta^2 - 1}

ζ\zeta 特征值 平衡点类型 物理含义
ζ>1\zeta > 1 两负实根 稳定结点 过阻尼,缓慢回复
ζ=1\zeta = 1 重负实根 稳定结点(退化) 临界阻尼,最快回复
0<ζ<10 < \zeta < 1 共轭复根(负实部) 稳定焦点 欠阻尼,衰减震荡
ζ=0\zeta = 0 纯虚根 中心点 无阻尼,等幅震荡
ζ<0\zeta < 0 共轭复根(正实部) 不稳定焦点 发散震荡→不稳定

四、李雅普诺夫稳定性

4.1 定义

平衡点 x\mathbf{x}^* 称为:

  • 稳定(Lyapunov stable):从附近出发的解永远停留在附近
  • 渐近稳定:从附近出发的解趋近于平衡点
  • 不稳定:存在任意靠近的初始条件导致的解远离平衡点

4.2 线性化稳定判据

对于线性化系统 y˙=Ay\dot{\mathbf{y}} = \mathbf{A}\mathbf{y}

所有特征值的实部为负 → 渐近稳定
存在特征值的实部为正 → 不稳定

这称为李雅普诺夫间接方法

4.3 在四旋翼中的应用

四旋翼的悬停状态就是一个平衡点:

x=[x0,y0,z0,0,0,0,ϕ0=0,θ0=0,ψ0,0,0,0]T \mathbf{x}^* = [x_0, y_0, z_0, 0, 0, 0, \phi_0=0, \theta_0=0, \psi_0, 0, 0, 0]^T

对 12 维系统做线性化,得到雅可比矩阵 A\mathbf{A} 的特征值——这告诉飞控设计者哪几个模态是稳定的、哪几个需要主动控制。


五、核心公式速查卡

公式 含义
x˙=f(x)\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}) 一阶ODE方程组
A=f/x\mathbf{A} = \partial\mathbf{f}/\partial\mathbf{x} 雅可比矩阵
y˙=Ay\dot{\mathbf{y}} = \mathbf{A}\mathbf{y} 线性化系统
det(AλI)=0\det(\mathbf{A} - \lambda\mathbf{I}) = 0 特征方程
所有 Re(λi)<0\text{Re}(\lambda_i) < 0 → 渐近稳定 李雅普诺夫间接方法

参考文献

  1. Perko, L. (2001). Differential Equations and Dynamical Systems (3rd ed., Chapter 1-2). Springer.
  2. Strogatz, S. H. (2018). Nonlinear Dynamics and Chaos (2nd ed., Chapter 5-6). CRC Press.
  3. Lyapunov, A. M. (1892). The General Problem of the Stability of Motion. — 稳定性理论的奠基之作
  4. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft Control and Simulation (3rd ed.). Wiley.

下一节:数值解法:欧拉法与RK4

大多数微分方程没有解析解——必须依靠数值方法。从最简单的欧拉法到工程标准的RK4,再到现代自适应变步长方法。