2. 动力系统基础回顾

在进入周期轨道识别之前，需要先把“系统如何随时间演化”这件事用统一的语言描述清楚。周期轨道识别与数值延拓的全部算法，最终都建立在这一节的几个概念之上：连续系统用流描述、离散系统用映射描述、平衡点是流的静止点、周期轨道是流的闭合回路，而庞加莱映射则把“找闭合回路”转化为“找映射的不动点”。理解这条转换链，后面的打靶法、延拓法才不会变成黑箱操作。

连续系统与离散系统：两种时间观

连续系统用常微分方程描述，状态 x ∈ R^n 随时间连续变化：

dx/dt = f(x, μ),   x(0) = x₀

其中 μ 是参数。给定初值 x₀，解是一条随时间演化的曲线，称为轨线；把所有初值的轨线放在一起看，就构成系统的流 φ_t(x₀)，含义是“从 x₀ 出发、经过时间 t 后到达的位置”。流满足两条基本性质：φ_0 是恒等映射，且 φ_{t+s} = φ_t ∘ φ_s。这两条性质正是“时间可以叠加”的数学表达，也是后续用流构造庞加莱映射的基础。

离散系统则用迭代映射描述，时间只取整数步：

x_{k+1} = g(x_k, μ)

它的解是一个点列 x₀, x₁, x₂, …。连续与离散并非互不相干：对连续系统以固定步长 h 采样，就得到离散映射 x_{k+1} = φ_h(x_k)，称为时间-h 映射或频闪映射。数值计算中我们几乎总是通过这种采样把连续问题转成离散问题来处理，因此“流”和“映射”两套语言必须能自由切换。

平衡点与周期轨道：两类特殊解

平衡点满足 f(x*) = 0，轨线停在那里不动。它的稳定性由雅可比矩阵 Df(x*) 的特征值判定：全部特征值实部为负则渐近稳定，存在实部为正则不稳定，出现实部为零的特征值则处于临界情形，线性化失效，需要非线性分析。

周期轨道是满足 φ_T(x*) = x* 且 T > 0 的最小正数 T 的解，轨线绕行一圈后精确回到起点，形成闭合回路。它与平衡点的关键区别在于：平衡点是“点”，周期轨道是“环”，环上每一点都不静止，但整体结构重复。周期轨道同样有稳定性概念，但判定不能只看单点雅可比，而要看一个周期内的累积效应，即单值矩阵（monodromy matrix）M = Φ_T(x*)，它是变分方程沿周期轨道的状态转移矩阵。对平面系统，M 的两个特征值为 1 和 Λ，其中 Λ 就是弗洛凯乘子；|Λ| < 1 对应稳定极限环，|Λ| > 1 对应不稳定极限环。特征值 1 对应沿轨道切向的扰动，它不改变轨道本身，因此总会出现。

庞加莱映射：把环变成点

直接在高维相空间里搜索闭合轨道代价很高，因为周期轨道是“一维对象”，用打靶法求解时方程数与未知数容易不匹配。庞加莱映射的思路是：在相空间中取一个横截面 Σ（通常选在周期轨道附近，要求流与截面横截相交），记录轨线每次从同一方向穿过 Σ 的落点，就得到一个 n-1 维的离散映射 P: Σ → Σ。此时，原系统的周期轨道恰好对应 P 的不动点 P(x*) = x*，而周期轨道的稳定性由 DP(x*) 的特征值决定，这些特征值与单值矩阵除特征值 1 之外的部分一致。

这一步转换的价值在于降维：原本要在 n 维相空间里找闭合曲线，现在只需在 n-1 维截面上解一个不动点方程。以平面系统为例，截面是一维的，问题退化为求一个标量函数的零点，用牛顿法即可高效求解。实际计算中，P(x) 没有解析表达式，只能通过“从 x 出发积分到下一次穿截面”来数值求值，因此每次迭代都要付出一次数值积分的代价，这也是周期轨道延拓计算量偏大的根源。

一个可动手的示例：用庞加莱映射求极限环

考虑平面系统（示例参数，用于说明流程，非实测数据）：

dx/dt = x - y - x(x² + y²)
dy/dt = x + y - y(x² + y²)

该系统在极坐标下化为 dr/dt = r(1 - r²)、dθ/dt = 1，存在单位圆 r = 1 的极限环，且 r < 1 时 dr/dt > 0、r > 1 时 dr/dt < 0，故该环稳定。取截面 Σ = { (x, y) : y = 0, x > 0 }，用 x 坐标作为截面坐标 s。从 (s, 0) 出发积分，直到 y 再次由负变正穿过 0，得到落点 P(s)。不动点方程 P(s) - s = 0 用牛顿法求解，导数 dP/ds 可用中心差分近似：

s = 1.2                      # 初始猜测
for k in range(20):
    P  = poincare(s)         # 数值积分到下一次穿截面
    Pp = poincare(s + 1e-6)
    Pm = poincare(s - 1e-6)
    dP = (Pp - Pm) / 2e-6    # 差分近似 dP/ds
    s  = s - (P - s) / (dP - 1)
    if abs(P - s) < 1e-10: break

收敛后应得到 s ≈ 1.0，与解析解一致，可作为代码正确性的校验。同时 dP/ds 在不动点处的值就是弗洛凯乘子 Λ，本例中应接近 e^{-4π} ≈ 3.5×10⁻⁶，远小于 1，说明极限环强稳定。这个“解析解校验 + 乘子判稳”的组合，是实际工作中验证周期轨道求解器是否写对的标准手段。

几个容易踩的坑

• 截面必须横截：若截面与轨道相切，穿截面次数可能为零或无穷，映射无定义。选取时应让流方向与截面法向的点积远离零。
• 周期未知时不能直接打靶：φ_T(x) = x 中 T 也未知，方程数与未知数不匹配。标准做法是补充相位条件（如固定截面）或把 T 作为未知量并加一个归一化方程。
• 特征值 1 不是不稳定：单值矩阵总有特征值 1，它对应沿轨道的相位漂移，判断稳定性时必须排除它，否则会误判所有周期轨道为临界。
• 数值积分精度影响映射精度：P(s) 的误差会直接进入牛顿迭代，积分容差应比求解容差至少紧一到两个量级。

把这几件事理清后，周期轨道问题就统一成了“求映射不动点”的问题，后续的延拓算法只需在这个框架上增加参数维度和预测-校正步骤即可。