4. 数值积分与轨迹生成

参数扫描的前提是能稳定地“算出轨迹”。动力系统通常以常微分方程（ODE）形式给出，例如 Duffing 振子：

dx/dt = v
dv/dt = -δ·v - α·x - β·x³ + γ·cos(ω·t)

除极少数线性系统外，这类方程没有解析解，必须用数值积分逐步推进状态。数值积分器决定了轨迹的精度、稳定性与计算成本，也直接决定后续行为识别（周期、混沌、分岔）是否可信。本节围绕四个问题展开：选什么解法、步长取多大、误差如何控制、轨迹数据如何组织。

一、欧拉法与龙格-库塔法：精度阶的差异

最直观的方法是显式欧拉法：用当前点的导数线性外推，x_{n+1} = x_n + h·f(t_n, x_n)。它每步只计算一次导数，实现简单，但只有一阶精度，局部截断误差为 O(h²)。对振荡系统，欧拉法会系统性地“注入能量”或“耗散能量”，长时间积分后振幅漂移明显，因此只适合演示或极短时间窗口。

实践中更常用的是四阶龙格-库塔法（RK4）。它在一步内取四个斜率：起点斜率 k₁、中点斜率 k₂ 与 k₃、终点斜率 k₄，再按 (k₁+2k₂+2k₃+k₄)/6 加权平均。RK4 的局部误差为 O(h⁵)，全局误差为 O(h⁴)，在相同步长下精度远高于欧拉法，代价是每步四次函数求值。对 Duffing 这类含三角函数与三次项的系统，RK4 在 h=0.01 量级即可获得视觉上光滑的相图。

若系统存在刚性（快慢时间尺度差异大，如大阻尼或化学反应动力学），显式方法会因稳定性限制被迫使用极小步长，此时应改用隐式方法（如隐式欧拉、Radau、BDF）。判断依据是：当步长稍大时轨迹发散或出现高频伪振荡，而减小步长后结果突变，往往提示刚性问题。

二、步长选择：精度、稳定性与采样密度的平衡

步长 h 不能凭感觉设定，应结合系统最快时间尺度估计。对含 cos(ωt) 激励的系统，激励周期 T=2π/ω，通常要求每周期至少采样 20～50 步，即 h ≤ T/20。对固有频率 ω₀ 的振子，同理取 h ≤ 2π/(20ω₀)。若系统时间尺度差异大，应按最快尺度取步长。

另一个约束是输出采样与积分步长分离。行为识别常需要等间隔轨迹点（如相图、Poincaré 截面），但积分步长可能更小。推荐做法是：积分器内部用小步长 h_int 保证精度，每隔固定时间 T_s 记录一个输出点，T_s 可大于 h_int。这样既保证数值稳定，又让轨迹数据规整、便于后续分析。

三、误差控制：从固定步长到自适应

固定步长实现简单，但无法兼顾“平缓段浪费、剧烈段不足”。自适应步长方法（如 RK45/Dormand-Prince）用两套不同阶数的解估计局部误差，若误差超过容差 atol+rtol·|x|，则缩小步长重算；若远小于容差，则放大步长。常用容差为 rtol=1e-6、atol=1e-9 量级，具体依问题尺度调整。

误差控制还要关注两个实践细节：其一，瞬态段（初始条件尚未落到吸引子上）常出现剧烈变化，应允许积分器自动缩步，不要强行固定步长；其二，长时间积分会累积舍入误差，对混沌系统，初值微小差异会指数放大，此时应把“轨迹可复现”理解为“统计行为可复现”，而非逐点一致。判断积分是否可信的实用手段是：把步长减半重算，若相图与关键指标（周期、振幅、Lyapunov 指数）变化在可接受范围，则认为收敛。

四、轨迹数据组织：让扫描结果可复用

参数扫描会产生大量轨迹，若数据组织混乱，后续分析将难以进行。推荐把每次积分结果整理为结构化表格，至少包含时间列、状态列、参数组标识和段标记。下面是一个可直接参考的字段设计（示例结构，非实测数据）：

其中“瞬态/稳态”标记非常关键：行为识别应只使用稳态段数据，否则初始条件的影响会污染周期与振幅统计。通常做法是丢弃前若干个激励周期（如 200 个周期）作为瞬态，再取后续数据计算指标。参数组编号则保证扫描结果可追溯，便于在分岔图上定位每个点对应的完整轨迹。

综上，数值积分不是“调个库函数”就结束的环节。解法阶数决定精度基线，步长决定稳定性与采样质量，误差控制决定结果可信度，数据组织决定扫描结果能否被高效分析。把这四件事做扎实，后续的行为识别与分岔分析才有可靠的地基。