1. 单摆自由振动仿真与相图绘制

无阻尼单摆的建模与数值积分

取摆长 L = 1.0 m、重力加速度 g = 9.81 m/s²，忽略阻尼与摆杆质量。以摆角 θ（rad）为广义坐标，从最低点逆时针为正，运动方程为：

θ'' = -(g/L)·sin(θ)

令状态向量 y = [θ, ω]，ω = θ'，则化为一阶方程组：

dθ/dt = ω
dω/dt = -(g/L)·sin(θ)

采用四阶龙格–库塔法（RK4）积分，步长 h = 0.005 s，积分时长 20 s。初始条件取 θ₀ = 0.3 rad、ω₀ = 0（小振幅）与 θ₀ = 2.5 rad、ω₀ = 0（大振幅）两组，用于对比周期与相轨迹形态。

import numpy as np

L, g = 1.0, 9.81

def deriv(y):
    theta, omega = y
    return np.array([omega, -(g/L)*np.sin(theta)])

def rk4(y, h):
    k1 = deriv(y)
    k2 = deriv(y + 0.5*h*k1)
    k3 = deriv(y + 0.5*h*k2)
    k4 = deriv(y + h*k3)
    return y + (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)

def simulate(theta0, omega0, T=20.0, h=0.005):
    n = int(T/h)
    y = np.array([theta0, omega0], dtype=float)
    ts = np.zeros(n+1); ys = np.zeros((n+1, 2))
    ts[0] = 0.0; ys[0] = y
    for i in range(n):
        y = rk4(y, h)
        ts[i+1] = (i+1)*h
        ys[i+1] = y
    return ts, ys

ts_s, ys_s = simulate(0.3, 0.0)   # 小振幅
ts_l, ys_l = simulate(2.5, 0.0)   # 大振幅

能量守恒校验与周期估计

单摆总能量（单位质量）为 E = 0.5·L²·ω² + g·L·(1 - cos θ)。无阻尼时 E 应保持常数，可用相对漂移 |E(t) - E(0)| / E(0) 判断积分精度。周期则通过检测 θ 由正变负的过零点间隔来估计。

def energy(theta, omega):
    return 0.5*(L**2)*omega**2 + g*L*(1 - np.cos(theta))

def period_estimate(ts, ys):
    theta = ys[:, 0]
    crossings = []
    for i in range(1, len(theta)):
        if theta[i-1] > 0 and theta[i] <= 0:
            # 线性插值求过零时刻
            t0, t1 = ts[i-1], ts[i]
            th0, th1 = theta[i-1], theta[i]
            crossings.append(t0 + (t1-t0)*th0/(th0-th1))
    if len(crossings) < 2:
        return None
    return np.mean(np.diff(crossings))

for name, ts, ys in [("小振幅 θ₀=0.3", ts_s, ys_s),
                     ("大振幅 θ₀=2.5", ts_l, ys_l)]:
    E0 = energy(ys[0,0], ys[0,1])
    drift = np.max(np.abs(energy(ys[:,0], ys[:,1]) - E0)) / E0
    T = period_estimate(ts, ys)
    print(f"{name}: 周期≈{T:.4f} s, 能量最大相对漂移={drift:.2e}")

预期输出（示例运行结果）：小振幅周期约 2.006 s，接近线性近似 2π√(L/g) ≈ 2.006 s；大振幅周期明显增大（约 2.5 s 量级），体现非线性带来的周期–振幅依赖。能量相对漂移应在 1e-8 量级，说明 RK4 在 h=0.005 s 下守恒性良好。

相平面轨迹与周期–振幅关系

以 θ 为横轴、ω 为纵轴绘制相图：小振幅轨迹接近椭圆，大振幅轨迹呈明显的非椭圆闭合曲线，且两条轨迹都闭合（无阻尼、无耗散）。周期随振幅增大而增大，可用椭圆积分 T = 4√(L/g)·K(sin(θ₀/2)) 校核。

import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
ax[0].plot(ys_s[:,0], ys_s[:,1], label="θ₀=0.3")
ax[0].plot(ys_l[:,0], ys_l[:,1], label="θ₀=2.5")
ax[0].set_xlabel("θ (rad)"); ax[0].set_ylabel("ω (rad/s)")
ax[0].set_title("相平面轨迹"); ax[0].legend(); ax[0].grid(alpha=0.3)

ax[1].plot(ts_s, ys_s[:,0], label="θ₀=0.3")
ax[1].plot(ts_l, ys_l[:,0], label="θ₀=2.5")
ax[1].set_xlabel("t (s)"); ax[1].set_ylabel("θ (rad)")
ax[1].set_title("时间序列"); ax[1].legend(); ax[1].grid(alpha=0.3)
plt.tight_layout(); plt.show()

结果判读要点

• 相轨迹闭合：无阻尼单摆是保守系统，相轨迹为闭合曲线，平衡点 (0,0) 是中心型奇点，不收敛也不发散。
• 能量守恒：RK4 下能量相对漂移应远小于 1e-5；若明显增大，说明步长过大或方程符号写错。
• 周期–振幅依赖：小振幅时 T 接近线性值 2.006 s；振幅增大后 T 单调增大，这是 sin θ 非线性项的直接体现，不能用线性模型解释。
• 大振幅轨迹形态：θ₀ 接近 π 时轨迹趋于分界线（separatrix），周期趋于无穷，仿真中应避免取 θ₀ 过于接近 π 导致数值困难。