4. 单摆平衡点线性化与局部稳定性判定

问题设定与线性化对象

无阻尼单摆方程取 θ'' + (g/L)·sinθ = 0，状态变量 x = [θ, ω]ᵀ，ω = θ'。写成状态方程：

θ' = ω
ω' = -(g/L)·sinθ

取 g = 9.81 m/s²、L = 1.0 m，则 ω₀² = g/L = 9.81 s⁻²。平衡点满足 ω = 0 且 sinθ = 0，在 [0, 2π) 内有两个：θ = 0（下垂，稳定）与 θ = π（倒立，不稳定）。

雅可比矩阵为

J(θ,ω) = [ 0        1        ]
         [ -ω₀²cosθ  0        ]

在 θ = 0 处 cosθ = 1，J 的特征值为 λ = ±i·ω₀ ≈ ±3.132i，纯虚根，线性化给出中心型结论；在 θ = π 处 cosθ = -1，特征值为 λ = ±ω₀ ≈ ±3.132，一正一负，鞍点，不稳定。注意：纯虚特征值属于线性化失效的临界情形，稳定性需借助能量函数（守恒量 E = ½ω² + ω₀²(1-cosθ)）判定，结论是稳定但非渐近稳定。

数值验证：小扰动下的局部行为

用四阶龙格–库塔积分两个平衡点附近的小扰动初值，观察 θ(t) 的演化。稳定点取 θ₀ = 0.1 rad，不稳定点取 θ₀ = π − 0.01 rad，ω₀ 均为 0，步长 h = 0.001 s，仿真 10 s。

import numpy as np

g, L = 9.81, 1.0
w2 = g / L

def f(x):
    theta, omega = x
    return np.array([omega, -w2 * np.sin(theta)])

def rk4(x, h):
    k1 = f(x)
    k2 = f(x + 0.5*h*k1)
    k3 = f(x + 0.5*h*k2)
    k4 = f(x + h*k3)
    return x + h*(k1 + 2*k2 + 2*k3 + k4)/6

def simulate(theta0, T=10.0, h=1e-3):
    x = np.array([theta0, 0.0])
    n = int(T/h)
    th = np.empty(n+1)
    th[0] = x[0]
    for i in range(n):
        x = rk4(x, h)
        th[i+1] = x[0]
    return np.linspace(0, T, n+1), th

t1, th1 = simulate(0.1)              # 稳定点附近
t2, th2 = simulate(np.pi - 0.01)     # 不稳定点附近

# 稳定点：检查振幅是否保持（能量守恒）
print("稳定点振幅范围:", th1.min(), th1.max())
# 不稳定点：检查偏离量是否指数增长
dev = np.abs(th2 - np.pi)
print("t=0 偏离:", dev[0], " t=3s 偏离:", dev[3000])

预期现象（由方程直接推出，非实测）：稳定点 θ(t) 在 ±0.1 rad 内等幅振荡，周期约 2π/ω₀ ≈ 2.006 s，与线性化预测一致；不稳定点偏离量按 e^{ω₀t} 增长，t = 3 s 时约 0.01·e^{9.4} ≈ 1.2×10³ rad，实际会因 sinθ 非线性在翻落过程中偏离指数律，但初始阶段与 λ = +3.132 的增长率吻合。

线性化有效性边界

线性化只在平衡点邻域成立。对 θ = 0，线性化周期 2.006 s 与真实周期 T = 4K(sin(θ₀/2))·√(L/g) 的偏差随振幅增大：θ₀ = 0.1 rad 时偏差约 0.06%，θ₀ = 0.5 rad 时约 1.6%，θ₀ = 1.0 rad 时约 7%。因此用线性化判定稳定性是可靠的，但用线性化预测大幅振荡的周期或波形则不可靠。对 θ = π，只要扰动足够小，指数增长率 ω₀ 在初始阶段始终成立，这是鞍点型不稳定的典型特征。