1. 单摆多体动力学建模与仿真

单摆的多体动力学建模思路

单摆是最简单的多体系统之一：一个无质量刚性杆连接固定铰点与集中质量。虽然可以用牛顿第二定律直接写出运动方程，但用多体动力学方法建模时，需要明确广义坐标、约束方程、质量矩阵和广义力，这套流程可以自然推广到双摆、连杆等更复杂系统。

取摆角 θ 为广义坐标（从竖直向下方向量起，逆时针为正）。铰点 O 固定，摆长为 L，质量 m 位于杆端。约束条件为质量点到铰点距离恒为 L，系统只有 1 个自由度。

对单摆，拉格朗日方程给出运动方程：

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

其中 g 为重力加速度，c 为阻尼系数（无阻尼时取 0）。这是一个二阶非线性常微分方程，除小角度近似外没有初等解析解，需要用数值积分求解。

数值积分实现

将二阶方程降为一阶方程组，令 y₁ = θ，y₂ = θ'：

y₁' = y₂
y₂' = -(g/L)·sin(y₁) - c·y₂

使用四阶龙格-库塔法（RK4）积分。以下代码完整可运行，输出摆角随时间的变化数据：

import numpy as np

def pendulum_deriv(state, g, L, c):
    theta, omega = state
    dtheta = omega
    domega = -(g / L) * np.sin(theta) - c * omega
    return np.array([dtheta, domega])

def rk4_step(state, dt, g, L, c):
    k1 = pendulum_deriv(state, g, L, c)
    k2 = pendulum_deriv(state + 0.5 * dt * k1, g, L, c)
    k3 = pendulum_deriv(state + 0.5 * dt * k2, g, L, c)
    k4 = pendulum_deriv(state + dt * k3, g, L, c)
    return state + (dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4)

# 参数设置
g = 9.81      # 重力加速度 m/s^2
L = 1.0       # 摆长 m
c = 0.0       # 无阻尼
theta0 = np.radians(30.0)  # 初始摆角 30°
omega0 = 0.0               # 初始角速度

dt = 0.001    # 积分步长 s
T_total = 5.0 # 仿真总时长 s
n_steps = int(T_total / dt)

state = np.array([theta0, omega0])
t_hist = np.zeros(n_steps + 1)
theta_hist = np.zeros(n_steps + 1)
t_hist[0] = 0.0
theta_hist[0] = theta0

for i in range(n_steps):
    state = rk4_step(state, dt, g, L, c)
    t_hist[i + 1] = (i + 1) * dt
    theta_hist[i + 1] = state[0]

# 输出部分结果
for i in range(0, n_steps + 1, 500):
    print(f"t={t_hist[i]:.3f}s  theta={np.degrees(theta_hist[i]):.4f} deg")

结果分析与验证

运行上述代码，可以观察以下规律：

• 周期验证：小角度下单摆周期 T ≈ 2π√(L/g) ≈ 2.006 s。初始 30° 时周期略大于该值，这是非线性效应。
• 能量守恒：无阻尼时总能量 E = ½mL²θ'² + mgL(1−cosθ) 应保持恒定。可用它检验积分精度。
• 阻尼影响：将 c 设为 0.5，摆角幅值会指数衰减，最终趋于零。

能量守恒检验代码：

def total_energy(theta, omega, g, L, m=1.0):
    kinetic = 0.5 * m * L**2 * omega**2
    potential = m * g * L * (1 - np.cos(theta))
    return kinetic + potential

# 在积分循环中记录能量
E0 = total_energy(theta0, omega0, g, L)
# 每步计算 E，观察 |E - E0| / E0 是否在 1e-6 量级

关键参数与结果判读