4. 相轨迹数值计算与可视化
相平面分析的核心困难在于：绝大多数非线性系统的相轨迹无法用解析式写出。我们通常只能拿到状态方程本身，例如二阶自治系统
dx/dt = f(x, y)
dy/dt = g(x, y)其中 x 为位移类状态，y 为速度类状态。相轨迹就是该向量场在 (x, y) 平面上的积分曲线。既然解析解不可得，就改用数值积分：从某个初始点出发，沿向量场方向一小步一小步地“走”，把走过的点连起来，就得到一条近似相轨迹。本节要解决两件事——怎样把这条轨迹算准，以及怎样把它和方向场一起画出来，让相图的整体结构一眼可见。
一、数值积分：把向量场离散成折线
先明确符号约定：本文所有相图均以 x 为横轴、y 为纵轴，时间 t 增大方向即向量场箭头方向；轨迹上的箭头表示 t 增加的方向，与向量场箭头一致。
最直接的方法是显式欧拉法。设步长为 h，第 k 步的状态为 (xk, yk)，则
x_{k+1} = x_k + h * f(x_k, y_k)
y_{k+1} = y_k + h * g(x_k, y_k)欧拉法实现简单，但只有一阶精度，且对刚性或快速旋转的轨迹（如无阻尼振子的圆轨迹）会明显“外扩”或“内缩”，长时间积分后轨迹形状失真。因此实际绘图更常用四阶龙格-库塔法（RK4）：
k1 = F(x_k, y_k)
k2 = F(x_k + h/2 * k1x, y_k + h/2 * k1y)
k3 = F(x_k + h/2 * k2x, y_k + h/2 * k2y)
k4 = F(x_k + h * k3x, y_k + h * k3y)
x_{k+1} = x_k + h/6 * (k1x + 2k2x + 2k3x + k4x)
y_{k+1} = y_k + h/6 * (k1y + 2k2y + 2k3y + k4y)其中 F 即向量场 (f, g)。RK4 每步计算四次向量场，代价约为欧拉法的四倍，但精度提升到四阶，在相同步长下轨迹闭合性明显更好。工程上常用折中：先用 RK4 保证形状，再用较小步长（如 h = 0.01～0.05 个时间单位）控制累积误差。
步长选择有一个实用判据：单步内状态变化不应超过相图特征尺度的百分之几。若系统特征频率为 ω，则建议 h < 0.1/ω。步长过大，轨迹会在平衡点附近“跳过”细节甚至发散；步长过小，计算量与绘图点数激增。可以先用粗步长试算，观察轨迹是否光滑闭合，再逐步减半验证结果稳定。
二、方向场：让结构先于轨迹显现
只画几条轨迹，读者很难判断整体走向。方向场（斜率场）在网格点上计算 (f, g) 的方向，用短线段或箭头表示，能在不积分的情况下展示向量场的骨架。绘制时通常对每个网格点做归一化，只保留方向、统一长度，否则大速度区域箭头过长会互相遮挡。
方向场与轨迹叠加后，相图的信息层次就清楚了：方向场给出“处处往哪走”，轨迹给出“从某点出发具体怎么走”，平衡点则是方向场中箭头长度趋零的位置。判断平衡点类型时，可结合该点附近方向场的汇聚、发散或旋转形态，与雅可比矩阵特征值结论互相印证。
三、一个可复现的示例
以无阻尼单摆的小角度近似为例（示例参数，非实测）：取方程
dx/dt = y
dy/dt = -sin(x)初始点选 (x₀, y₀) = (2.0, 0)，步长 h = 0.02，用 RK4 积分约 3000 步。由于系统保守，轨迹应是一条闭合曲线，围绕原点（稳定平衡点）旋转；若改用欧拉法同样步长，轨迹会逐渐向外螺旋，这就是数值耗散/增益造成的假象，可作为方法对比的直观证据。
绘图流程可归纳为：① 确定绘图窗口，如 x ∈ [-4, 4]、y ∈ [-3, 3]；② 在窗口内均匀取 20×20 网格计算方向场；③ 对若干初始点分别积分，得到多条轨迹；④ 将方向场、轨迹、平衡点分层叠加，轨迹用不同颜色区分初始条件；⑤ 在轨迹上按固定步数间隔加箭头，标明时间方向。
四、结果判读与常见问题
叠加图完成后，判读顺序建议是：先看方向场中箭头汇聚、发散或旋转的区域，定位平衡点；再看轨迹是否闭合、是否被吸引到某点或极限环；最后检查轨迹与方向场是否一致——若某条轨迹明显逆着方向场箭头走，多半是积分步长过大或符号写反。
几个高频问题值得提前防范。其一，坐标轴比例不一致会把圆轨迹画成椭圆，误判为极限环，绘图时应保持等比例或明确标注。其二，轨迹在平衡点附近速度趋零，等步长积分会“走不动”，可改用自适应步长或按弧长参数化。其三，多条轨迹叠加时颜色过多反而混乱，建议主轨迹用饱和色、辅助轨迹用低透明度，方向场用浅灰，形成清晰的信息层级。
至此，数值计算与可视化形成了闭环：方程给出向量场，RK4 把向量场积成轨迹，方向场提供全局参照，叠加绘图把两者合成一张可判读的相图。后续分析平衡点稳定性、判断极限环存在性时，这张图就是最直接的出发点。