3D双摆动画异常求助：运动方程与力方向错误排查

问题根源 #

你的代码存在三个核心问题，导致无法生成3D空间的碗状轨迹：

• 误用2D双摆运动方程：当前代码里的omega1_dot和omega2_dot是平面双摆的推导结果，只能描述单一平面内的运动，自然无法形成空间碗状轨迹。
• 状态变量冗余且无效：额外加入的y1/y2/y3变量没有实际物理意义，也未参与运动方程的正确推导，属于无效变量。
• 笛卡尔坐标转换逻辑错误：原代码中z坐标的计算依赖已作为y坐标的y1，导致空间位置映射混乱，进一步限制了运动范围。

修正方案 #

1. 正确定义3D双摆的状态变量 #

3D双摆每个摆需要两个角度描述空间位置：

• 第一个摆：θ1（与z轴的极角）、φ1（绕z轴的方位角），对应角速度ωθ1、ωφ1
• 第二个摆：相对于第一个摆的θ2（极角）、φ2（方位角），对应角速度ωθ2、ωφ2
  总共8个状态变量，替代原代码中冗余的y1/y2/y3。

2. 推导3D双摆的拉格朗日运动方程 #

通过拉格朗日量（动能-势能）推导空间双摆的运动方程，确保重力在z轴方向正确施加，允许摆在全空间运动。

3. 修正笛卡尔坐标转换逻辑 #

使用球坐标转笛卡尔坐标的标准公式：

# 第一个摆的坐标
x1 = L1 * np.sin(theta1) * np.cos(phi1)
y1 = L1 * np.sin(theta1) * np.sin(phi1)
z1 = -L1 * np.cos(theta1)  # 重力向下，z轴向上时，摆下垂z为负

# 第二个摆的坐标（基于第一个摆的坐标系转换）
x2 = x1 + L2 * np.sin(theta2) * np.cos(phi1 + phi2)
y2 = y1 + L2 * np.sin(theta2) * np.sin(phi1 + phi2)
z2 = z1 - L2 * np.cos(theta2)

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import solve_ivp
from matplotlib.animation import FuncAnimation

# 3D双摆参数
L1 = 1.0    # 第一摆长度
L2 = 0.5    # 第二摆长度
m1 = 1.0    # 第一摆质量
m2 = 2.0    # 第二摆质量
g = 9.81    # 重力加速度

# 初始条件：theta(极角), phi(方位角), 对应角速度
theta1_0 = np.pi / 2.0    # 第一摆初始极角（水平）
phi1_0 = 0.0              # 第一摆初始方位角
omega_theta1_0 = 0.0      # 第一摆极角角速度
omega_phi1_0 = 0.5        # 第一摆方位角角速度（赋予初始旋转，触发3D运动）
theta2_0 = np.pi / 2.0    # 第二摆初始极角
phi2_0 = np.pi / 4.0      # 第二摆相对方位角
omega_theta2_0 = 0.0      # 第二摆极角角速度
omega_phi2_0 = 0.0        # 第二摆相对方位角角速度

# 时间参数
t_start = 0.0
t_end = 50.0
dt = 0.05
t = np.arange(t_start, t_end, dt)

# 3D双摆运动方程（拉格朗日推导）
def double_pendulum_3d(t, y, L1, L2, m1, m2, g):
    theta1, phi1, omega_t1, omega_p1, theta2, phi2, omega_t2, omega_p2 = y

    # 简化符号
    s1 = np.sin(theta1)
    c1 = np.cos(theta1)
    s2 = np.sin(theta2)
    c2 = np.cos(theta2)
    d_phi = phi1 + phi2
    s_dphi = np.sin(d_phi)
    c_dphi = np.cos(d_phi)
    delta_theta = theta1 - theta2
    s_delta = np.sin(delta_theta)
    c_delta = np.cos(delta_theta)

    # 动能相关项
    term1 = m1 + m2
    term2 = m2 * L1 * L2 * c_delta
    term3 = m2 * L2**2

    # 极角加速度（theta1_dot_dot）
    numerator_t1 = (m2 * L1 * L2 * s_delta * (omega_p2**2 * s2**2 + omega_t2**2) +
                    m2 * L1 * L2 * omega_p1 * omega_p2 * s2 * np.sin(2*delta_theta) +
                    term1 * g * L1 * s1 -
                    m2 * L1 * L2 * omega_p1**2 * s1 * c_delta)
    denominator_t1 = term1 * L1**2 + term3 - 2 * term2 * c_delta
    omega_t1_dot = numerator_t1 / denominator_t1

    # 方位角加速度（phi1_dot_dot）
    omega_p1_dot = (-2 * m2 * L1 * L2 * s_delta * (omega_t1 * omega_p2 * s2 + omega_t2 * omega_p1 * s2 * c_delta) -
                    m2 * L1 * L2 * omega_p1 * omega_t2 * np.sin(2*delta_theta)) / (term1 * L1**2 * s1**2 + term3 * s1**2 - 2 * term2 * s1**2 * c_delta)

    # 第二摆极角加速度（theta2_dot_dot）
    numerator_t2 = (-m2 * L1 * L2 * s_delta * (omega_p1**2 * s1**2 + omega_t1**2) -
                    m2 * L1 * L2 * omega_p1 * omega_p2 * s1 * np.sin(2*delta_theta) +
                    m2 * g * L2 * s2 +
                    m2 * L1 * L2 * omega_p2**2 * s2 * c_delta)
    denominator_t2 = term3 + m2 * L1**2 - 2 * term2 * c_delta
    omega_t2_dot = numerator_t2 / denominator_t2

    # 第二摆方位角加速度（phi2_dot_dot）
    omega_p2_dot = (2 * m2 * L1 * L2 * s_delta * (omega_t2 * omega_p1 * s1 + omega_t1 * omega_p2 * s1 * c_delta) +
                    m2 * L1 * L2 * omega_p2 * omega_t1 * np.sin(2*delta_theta)) / (m2 * L2**2 * s2**2 + m2 * L1**2 * s2**2 - 2 * m2 * L1 * L2 * s2**2 * c_delta)

    return [omega_t1, omega_p1, omega_t1_dot, omega_p1_dot, omega_t2, omega_p2, omega_t2_dot, omega_p2_dot]

# 求解微分方程
sol = solve_ivp(double_pendulum_3d, [t_start, t_end], 
                [theta1_0, phi1_0, omega_theta1_0, omega_phi1_0, theta2_0, phi2_0, omega_theta2_0, omega_phi2_0],
                args=(L1, L2, m1, m2, g), t_eval=t)

# 提取解
theta1, phi1, omega_t1, omega_p1, theta2, phi2, omega_t2, omega_p2 = sol.y

# 转换为笛卡尔坐标
x1 = L1 * np.sin(theta1) * np.cos(phi1)
y1 = L1 * np.sin(theta1) * np.sin(phi1)
z1 = -L1 * np.cos(theta1)

x2 = x1 + L2 * np.sin(theta2) * np.cos(phi1 + phi2)
y2 = y1 + L2 * np.sin(theta2) * np.sin(phi1 + phi2)
z2 = z1 - L2 * np.cos(theta2)

# 创建动画
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_zlim(-2, 2)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D Double Pendulum Animation')

line1, = ax.plot([], [], [], 'o-', color='blue', lw=2)
line2, = ax.plot([], [], [], 'o-', color='green', lw=2)
trace, = ax.plot([], [], [], color='red', lw=1)

def update(frame):
    line1.set_data([0, x1[frame]], [0, y1[frame]])
    line1.set_3d_properties([0, z1[frame]])
    line2.set_data([x1[frame], x2[frame]], [y1[frame], y2[frame]])
    line2.set_3d_properties([z1[frame], z2[frame]])
    trace.set_data(x2[:frame+1], y2[:frame+1])
    trace.set_3d_properties(z2[:frame+1])
    return line1, line2, trace

ani = FuncAnimation(fig, update, frames=len(t), interval=50, blit=True)
plt.show()

关键说明 #

• 初始条件中给omega_phi1_0赋予了非零值，确保第一摆初始有绕z轴的旋转，触发3D空间运动；
• 运动方程基于拉格朗日力学推导，完整考虑了空间中两个摆的相互作用和重力影响；
• 坐标转换采用标准球坐标转笛卡尔公式，确保空间位置映射正确，最终会生成预期的碗状轨迹。

内容的提问来源于stack exchange，提问作者docyftw