4. 线性稳定性与变分方程

要计算混沌系统的Lyapunov指数，第一步不是直接写代码跑数值积分，而是先回答一个更基础的问题：系统在平衡点附近到底“稳不稳”？线性稳定性分析给出了这个答案，而变分方程则把这种局部线性化的思想推广到整条轨迹上，最终导出Lyapunov指数的连续系统定义。本节沿着“平衡点→线性化→雅可比矩阵→变分方程→Lyapunov指数”这条主线展开，每一步都给出可操作的计算细节。

一、平衡点与线性稳定性

考虑连续自治系统 ẋ = f(x)，x ∈ ℝⁿ。平衡点 x* 满足 f(x*) = 0，即系统在该点处“静止”。但静止不等于稳定：若给 x* 一个微小扰动 δx，系统是把它拉回来、推出去，还是让它绕圈，取决于 f 在 x* 附近的局部行为。

设 x(t) = x* + δx(t)，对 f 在 x* 处做一阶泰勒展开：

f(x* + δx) ≈ f(x*) + J(x*)·δx = J(x*)·δx

其中 J(x*) 是 f 在 x* 处的雅可比矩阵，元素为 J_ij = ∂f_i/∂x_j。于是扰动满足线性方程 δẋ = J(x*)·δx，其解为 δx(t) = e^{J(x*)t}·δx(0)。稳定性完全由 J(x*) 的特征值 λ_i 决定：

• 所有 Re(λ_i) < 0：扰动指数衰减，平衡点渐近稳定（吸引子）；
• 存在 Re(λ_i) > 0：扰动沿该特征方向指数增长，平衡点不稳定；
• 存在 Re(λ_i) = 0 且其余为负：临界情形，线性分析失效，需借助中心流形或非线性项判断。

以经典Lorenz系统为例，参数取 σ = 10、ρ = 28、β = 8/3。原点 (0,0,0) 处雅可比矩阵为 diag(-10, -1, -8/3)，三个特征值全为负，原点稳定。另外两个平衡点 C± = (±√(β(ρ−1)), ±√(β(ρ−1)), ρ−1) ≈ (±8.485, ±8.485, 27)，其雅可比矩阵的特征值为一个负实数和一对共轭复根，实部约为 −13.85 与 +0.094。正实部的出现说明 C± 是不稳定的鞍焦点——这正是Lorenz系统产生混沌的必要条件之一。

二、从平衡点到轨迹：变分方程

平衡点分析只回答“静止点附近如何”，而混沌吸引子上的轨迹从不停留。要刻画沿任意轨迹 x(t) 的扰动演化，需把线性化推广到整条轨迹。设 x(t) 是系统的一条解，在其上叠加小扰动 δx(t)，代入 ẋ + δẋ = f(x + δx) 并保留一阶项，得到变分方程：

δẋ = J(x(t))·δx

与平衡点情形不同，这里的雅可比矩阵 J(x(t)) 沿轨迹随时间变化，因此不能简单写成矩阵指数。变分方程是一个线性时变系统，它的解可形式地写为 δx(t) = Φ(t, t₀)·δx(t₀)，其中 Φ 称为状态转移矩阵，满足 Φ̇ = J(x(t))·Φ，Φ(t₀, t₀) = I。

实际计算中，我们并不显式构造 Φ，而是把变分方程与原系统联立，作为 2n 维系统一起数值积分：

ẋ = f(x)
Φ̇ = J(x)·Φ,   Φ(0) = I

以Lorenz系统为例，雅可比矩阵为：

J(x,y,z) = [[-σ,   σ,   0],
            [ρ-z, -1,  -x],
            [y,    x,  -β]]

在每一步积分中，用当前 (x, y, z) 实时计算 J，再推进 Φ 的一列或整个矩阵。这就是后续章节数值计算Lyapunov指数的核心引擎。

三、Lyapunov指数的连续系统定义

有了状态转移矩阵，就可以定义Lyapunov指数。取初始扰动方向为单位向量 v₀，考察扰动范数随时间的增长：

λ(v₀) = lim_{t→∞} (1/t) · ln( ‖δx(t)‖ / ‖δx(0)‖ )

由于 δx(t) = Φ(t,0)·δx(0)，上式等价于 λ(v₀) = lim (1/t)·ln‖Φ(t,0)·v₀‖。对 n 维系统，存在 n 个Lyapunov指数 λ₁ ≥ λ₂ ≥ … ≥ λₙ，它们由 ΦᵀΦ 的特征值 μ_i(t) 的对数增长率给出：λ_i = lim (1/2t)·ln μ_i(t)。最大Lyapunov指数 λ₁ > 0 是判断混沌的定量标志；对Lorenz系统，λ₁ ≈ 0.906，λ₂ = 0，λ₃ ≈ −14.57，三者之和等于散度 ∇·f = −σ − 1 − β ≈ −13.67，与理论一致。

四、实操要点：把定义变成可运行的计算

把上述定义落到代码层面，关键有三点。第一，雅可比矩阵要解析给出，而不是数值差分，否则误差会污染指数估计；对Lorenz系统直接按上面的 3×3 公式写函数即可。第二，变分方程与原系统同步积分，推荐用四阶Runge–Kutta，步长 h = 0.01，积分时长先取 T = 1000 个时间单位，前 100 个单位作为暂态丢弃。第三，必须做重正交化：Φ 的列向量会迅速趋于同一方向（数值上全部对齐到最大增长方向），若不处理，λ₂、λ₃ 会被 λ₁ 淹没。标准做法是每隔 Δt = 1 个时间单位对 Φ 做QR分解，把 R 的对角元累加进指数估计：

for each step:
    x, Phi = rk4_step(f, J, x, Phi, h)
    if t % 1 == 0:
        Q, R = qr(Phi)
        log_sum[i] += log(abs(R[i,i]))
        Phi = Q
lambda_i = log_sum / T

判断结果是否可信，可做两项自检：其一，三个指数之和应接近散度 −σ−1−β ≈ −13.67，偏差超过 5% 说明步长过大或积分时长不足；其二，λ₂ 应接近 0（沿轨迹方向的扰动既不增长也不衰减），若明显偏离 0，通常是重正交化间隔太长。这两条检查在后续章节的数值实验中会反复用到。

至此，我们完成了从“平衡点稳不稳”到“轨迹上扰动平均增长多快”的完整推导链：平衡点给出局部判据，雅可比矩阵提供线性化工具，变分方程把线性化沿轨迹延拓，Lyapunov指数则是对状态转移矩阵长期增长率的极限刻画。这条链条既是理论骨架，也是后续数值实现的直接依据。