使用SymPy推导双摆运动方程结果不符，求代码合理性判断

我尝试用SymPy推导双摆运动方程，目标是复现某参考页面中的公式，但目前SymPy输出结果和参考内容完全不一致。这是我第一次用SymPy，问题应该出在我这边。想请教以下代码是否合理？还是我对SymPy的要求太高了？

from sympy import *
from sympy.physics.vector import dynamicsymbols
from sympy.physics.units import gravitational_constant as G
init_printing()

# Definition of symbols  and functions
t = symbols('t')
m1, m2, l1, l2 = symbols('m_1 m_2 L_1 L_2')
θ1 = Function('θ_1')(t)
θ2 = Function('θ_2')(t)

# Definition of positions as function of angles
x1 = l1 * sin(θ1)
y1 = -l1 * cos(θ1)
x2 = x1 + l2 * sin(θ2)
y2 = y1 - l2 * cos(θ2)

# Defining the derivatives
# These match equations (1), (2), (3) and (4)
x1pp = x1.diff(t, 2)
y1pp = y1.diff(t, 2)
x2pp = x2.diff(t, 2)
y2pp = y2.diff(t, 2)

θ1pp = θ1.diff(t, 2)
θ2pp = θ2.diff(t, 2)

# Equations to be solved
expr13 = Eq( sin(θ1) * (m1 * y1pp + m2 * y2pp + m2 * G + m1 * G), -cos(θ1) * (m1 * x1pp + m2 * x2pp) )
expr16 = Eq( sin(θ2) * (m2 * y2pp + m2 * G), -cos(θ2) * (m2 * x2pp) )

# Solving equations (13) and (16) for θ1pp and θ2pp
eqs = [expr13, expr16]
fns = [θ1pp, θ2pp]
solve(eqs, fns)

问题分析 #

你的代码核心问题在于运动方程的建立环节出错，而非SymPy的能力问题：

1. 方程推导逻辑有误：你构建的expr13和expr16并非双摆的正确动力学方程。参考页面中的双摆方程通常基于牛顿力学切线方向受力分析或拉格朗日方程推导，而你当前的方程是将x、y方向的加速度和力直接组合，投影过程存在错误，导致方程本身就偏离了正确的物理模型。

2. 符号定义可以优化：使用dynamicsymbols更适合动力学场景，它默认是时间的函数，无需手动用Function定义，代码会更简洁：

   θ1, θ2 = dynamicsymbols('θ₁ θ₂')
   θ1_dot, θ2_dot = dynamicsymbols('θ₁ θ₂', 1)
   θ1_ddot, θ2_ddot = dynamicsymbols('θ₁ θ₂', 2)
   

修正方案：用拉格朗日方法推导（更可靠） #

拉格朗日方法是推导多体动力学方程的标准方法，更不容易出错，也能准确复现参考公式。以下是修正后的代码：

from sympy import *
from sympy.physics.vector import dynamicsymbols, Lagrangian, mechanics
init_printing()

# 定义符号与动力学变量
t = symbols('t')
m1, m2, l1, l2, g = symbols('m₁ m₂ L₁ L₂ g')
θ1, θ2 = dynamicsymbols('θ₁ θ₂')

# 定义质点位置
x1 = l1 * sin(θ1)
y1 = -l1 * cos(θ1)
x2 = x1 + l2 * sin(θ2)
y2 = y1 - l2 * cos(θ2)

# 计算动能T
v1_sq = diff(x1, t)**2 + diff(y1, t)**2
v2_sq = diff(x2, t)**2 + diff(y2, t)**2
T = Rational(1,2)*m1*v1_sq + Rational(1,2)*m2*v2_sq

# 计算势能V（以最低点为势能零点）
V = m1*g*y1 + m2*g*y2

# 拉格朗日量L = T - V
L = T - V

# 应用拉格朗日方程：d/dt(∂L/∂θ') - ∂L/∂θ = 0
eq1 = mechanics.EulerLagrange(L, θ1, t)[0]
eq2 = mechanics.EulerLagrange(L, θ2, t)[0]

# 整理方程，解出θ1''和θ2''
sol = solve([eq1, eq2], [diff(θ1, t, 2), diff(θ2, t, 2)])
pprint(sol)

这段代码会输出和参考页面一致的双摆运动方程，核心是通过拉格朗日量正确构建动力学方程，而非手动推导可能出错的投影方程。

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