3. 平衡点稳定性与线性化分析

对非线性系统做稳定性分析，最常用的工程手段是“先求平衡点，再在平衡点处线性化，用雅可比矩阵特征值判断局部稳定性”。这一节我们用一个具体的二阶系统把整条链路走通：从平衡点方程、雅可比矩阵推导，到特征值判据，再到数值仿真验证线性化预测的适用范围与失效边界。

一、系统与平衡点

考虑带阻尼的单摆方程（取摆长与重力加速度归一化，状态为角度与角速度）：

θ' = ω
ω' = -sin(θ) - 0.3·ω

平衡点要求 θ' = 0、ω' = 0，即 ω = 0 且 sin(θ) = 0。在 [0, 2π) 内得到两个平衡点：P₁ = (0, 0)（下垂位置）与 P₂ = (π, 0)（倒立位置）。注意平衡点是“状态空间中的点”，不是轨迹，求平衡点本质上是解代数方程组，与积分无关。

二、雅可比矩阵与特征值判据

记 f₁ = ω，f₂ = -sin(θ) - 0.3ω，则雅可比矩阵为

J(θ, ω) = [ ∂f₁/∂θ  ∂f₁/∂ω ] = [ 0        1    ]
          [ ∂f₂/∂θ  ∂f₂/∂ω ]   [ -cos(θ) -0.3 ]

在 P₁ 处 cos(0) = 1，J₁ = [[0, 1], [-1, -0.3]]，特征方程 λ² + 0.3λ + 1 = 0，特征值 λ = -0.15 ± 0.9887i。实部为负、虚部非零，故 P₁ 是稳定焦点：小扰动下系统以衰减振荡回到下垂位置。

在 P₂ 处 cos(π) = -1，J₂ = [[0, 1], [1, -0.3]]，特征方程 λ² + 0.3λ - 1 = 0，特征值 λ₁ ≈ 0.861、λ₂ ≈ -1.161。一正一负，故 P₂ 是鞍点：稳定流形与不稳定流形在此交叉，倒立位置对扰动极其敏感。

三、数值仿真验证

用四阶龙格–库塔法（步长 h = 0.01，仿真时长 20 s）分别从两个平衡点附近施加小扰动，观察轨迹行为。下面给出可直接运行的 Python 片段：

import numpy as np

def f(s):
    th, w = s
    return np.array([w, -np.sin(th) - 0.3*w])

def rk4(s0, T=20.0, h=0.01):
    n = int(T/h); s = np.array(s0, float); traj = [s.copy()]
    for _ in range(n):
        k1 = f(s); k2 = f(s + h/2*k1)
        k3 = f(s + h/2*k2); k4 = f(s + h*k3)
        s = s + h/6*(k1 + 2*k2 + 2*k3 + k4)
        traj.append(s.copy())
    return np.array(traj)

# 下垂点附近：初值 (0.1, 0)
t1 = rk4([0.1, 0.0])
# 倒立点附近：初值 (pi + 0.01, 0)
t2 = rk4([np.pi + 0.01, 0.0])

预期结果（示例性数值，供对照）：从 (0.1, 0) 出发，θ 以约 6.3 s 的周期振荡并逐步衰减，20 s 后 |θ| 已降至 10⁻³ 量级，与“稳定焦点”预测一致；从 (π + 0.01, 0) 出发，θ 在约 4 s 内迅速偏离 π，最终翻越到 0 附近并进入下垂点的小振荡，说明倒立点附近的微小扰动被指数放大——这正是鞍点不稳定方向的体现。

四、线性化预测的适用范围与局限

线性化结论是局部的：它只保证在平衡点充分小的邻域内成立，且要求平衡点是双曲的。本例中，若把下垂点的初始角位移从 0.1 rad 增大到 1.5 rad，线性模型预测的衰减速率与振荡频率会明显偏离真实轨迹，因为 sin(θ) ≈ θ 的近似误差随幅度增大而累积。更关键的是，线性化无法回答“吸引域有多大”——即多大的初始扰动仍能被拉回 P₁。这需要借助能量函数或数值打靶来估计。

另一类局限出现在非双曲平衡点。若把阻尼项改为 0（无阻尼单摆），P₁ 处特征值为 ±i，实部为零，线性化只能给出“中心”这一临界判断，而真实系统是否为周期轨道、是否被非线性项拉向稳定，必须用中心流形约化或守恒量分析。这说明：特征值判据是充分条件而非必要条件，遇到零实部必须换工具。

五、实操要点与常见误区

第一，求平衡点时要考虑周期性。单摆的 θ 具有 2π 周期，平衡点应写成 (2kπ, 0) 与 ((2k+1)π, 0)，只取一个周期内的代表点做分析即可，但仿真时若初值跨越周期边界，需注意角度回绕处理。第二，雅可比矩阵必须在平衡点处求值，不能在任意点求值后套用判据。第三，特征值计算建议用数值库（如 numpy.linalg.eigvals）复核手算结果，尤其是高阶系统。第四，仿真验证时应同时记录状态轨迹与“线性化预测轨迹”，把两者叠在同一张图上对比，才能直观看出近似何时开始失效。

下表汇总本例的关键判断依据，可作为同类二阶系统的分析模板：

最后强调：线性化分析给出的是“局部定性结论”，它回答“附近扰动会不会回来”，但不回答“回来得多快、能容忍多大扰动、全局是否存在其他吸引子”。工程实践中应把特征值判据与数值仿真、能量函数或分岔分析结合使用，才能对非线性系统的稳定性形成完整判断。