2. 双摆系统的混沌运动仿真

双摆的建模思路与状态变量

双摆由两根无质量刚性杆和两个集中质量组成，杆长分别为 L1、L2，质量分别为 m1、m2。取两杆与竖直向下方向的夹角 θ1、θ2 为广义坐标，角速度 ω1、ω2 为广义速度，则系统状态向量为 y = [θ1, ω1, θ2, ω2]，共 4 个状态量，构成 4 维一阶常微分方程组。

双摆的动力学方程可由拉格朗日方程推导。设重力加速度为 g，忽略摩擦与杆的质量，角加速度满足如下耦合关系：

dω1/dt = [ -g(2m1+m2)sinθ1 - m2·g·sin(θ1-2θ2)
           - 2·sin(θ1-θ2)·m2·(ω2²·L2 + ω1²·L1·cos(θ1-θ2)) ]
         / [ L1·(2m1 + m2 - m2·cos(2θ1-2θ2)) ]

dω2/dt = [ 2·sin(θ1-θ2)·( ω1²·L1·(m1+m2)
           + g·(m1+m2)·cosθ1 + ω2²·L2·m2·cos(θ1-θ2) ) ]
         / [ L2·(2m1 + m2 - m2·cos(2θ1-2θ2)) ]

该方程组的显著特点是两个自由度之间存在非线性耦合项，且分母含有 cos(2θ1-2θ2)，这正是混沌产生的数学根源。下面用四阶龙格-库塔法（RK4）对状态方程进行数值积分。

数值积分实现

采用固定步长 RK4，步长取 dt = 0.001 s，仿真时长 T = 20 s。参数取 L1 = L2 = 1.0 m，m1 = m2 = 1.0 kg，g = 9.81 m/s²。初始条件设置两组仅相差 1e-3 rad 的微小扰动，用于观察轨迹分离。

import numpy as np

def derivs(state, L1, L2, m1, m2, g):
    th1, w1, th2, w2 = state
    d = 2*m1 + m2 - m2*np.cos(2*th1 - 2*th2)
    a1 = (-g*(2*m1+m2)*np.sin(th1)
          - m2*g*np.sin(th1-2*th2)
          - 2*np.sin(th1-th2)*m2*(w2**2*L2 + w1**2*L1*np.cos(th1-th2))) / (L1*d)
    a2 = (2*np.sin(th1-th2)*(w1**2*L1*(m1+m2)
          + g*(m1+m2)*np.cos(th1)
          + w2**2*L2*m2*np.cos(th1-th2))) / (L2*d)
    return np.array([w1, a1, w2, a2])

def rk4_step(state, dt, L1, L2, m1, m2, g):
    k1 = derivs(state, L1, L2, m1, m2, g)
    k2 = derivs(state + 0.5*dt*k1, L1, L2, m1, m2, g)
    k3 = derivs(state + 0.5*dt*k2, L1, L2, m1, m2, g)
    k4 = derivs(state + dt*k3, L1, L2, m1, m2, g)
    return state + dt/6.0*(k1 + 2*k2 + 2*k3 + k4)

def simulate(y0, T=20.0, dt=0.001):
    n = int(T/dt)
    traj = np.zeros((n+1, 4))
    traj[0] = y0
    for i in range(n):
        traj[i+1] = rk4_step(traj[i], dt, 1.0, 1.0, 1.0, 1.0, 9.81)
    return traj

y0_a = np.array([np.pi/2, 0.0, np.pi/2, 0.0])
y0_b = np.array([np.pi/2 + 1e-3, 0.0, np.pi/2, 0.0])
traj_a = simulate(y0_a)
traj_b = simulate(y0_b)
sep = np.abs(traj_a[:, 0] - traj_b[:, 0])
print("t=5s  θ1偏差:", sep[5000])
print("t=10s θ1偏差:", sep[10000])
print("t=20s θ1偏差:", sep[20000])

运行上述代码可得到两组轨迹的 θ1 偏差随时间的变化。预期现象：初始偏差仅 1e-3 rad，在 t ≈ 5 s 时偏差仍较小，t ≈ 10 s 时明显增大，t ≈ 20 s 时两组轨迹已完全分离。这正是混沌系统对初值敏感依赖的典型表现。

混沌现象的可视化与判读

下图给出双摆的几何结构与两组初值下 θ1(t) 的分离过程示意。左侧为双摆的杆-质量拓扑，右侧为偏差随时间的指数增长趋势。判读要点：偏差曲线在早期近似指数增长，其斜率对应最大李雅普诺夫指数；当偏差增长到与系统尺度相当时，两条轨迹不再相关。

结果判读与验证方法

判断混沌是否出现的核心依据是偏差曲线的形态：若偏差在早期呈近似指数增长，且增长速率对初值扰动不敏感（即不同扰动幅度下斜率一致），则可认为系统处于混沌状态。验证时建议做以下检查：

• 将步长减半（dt = 0.0005 s），若偏差曲线形态基本不变，说明数值积分稳定，混沌现象不是数值误差造成的。
• 将初始扰动改为 1e-4 rad 和 1e-2 rad，观察偏差增长斜率是否一致；若一致，则最大李雅普诺夫指数为正。
• 绘制相图 (θ1, ω1)，混沌状态下相轨迹不会闭合，而是填充一定区域。