Sympy lambdify报'Derivative未定义'，3D双摆数值求解受阻

问题概述 #

在Jupyter Notebook开发3D双摆计算与动画程序时，先通过分解拉格朗日系统得到各角加速度的符号解并合并，随后尝试用sympy.lambdify将符号解转为数值解时，触发“Derivative is not defined”错误。尽管尝试了导数项替换操作，表达式中仍残留导数，无法使用odeint进行数值求解。

符号解求解代码及输出 #

以下是求解各角加速度符号解的核心代码：

sol_the1_dd = smp.solve(l1_O, O1_dd, dict=True, simplify=False, rational=False)

# 代入the1_dd的解到LE2，求解the2_dd
LE2_substituted = l2_O.subs(O1_dd, sol_the1_dd[0][O1_dd])
sol_the2_dd = smp.solve(LE2_substituted, O2_dd, dict=True, simplify=False, rational=False)

# 代入the2_dd的解到LE3，求解phi1_dd
LE3_substituted = l1_f.subs(O2_dd, sol_the2_dd[0][O2_dd])
sol_phi1_dd = smp.solve(LE3_substituted, f1_dd, dict=True, simplify=False, rational=False)

# 代入phi1_dd的解到LE4，求解phi2_dd
LE4_substituted = l2_f.subs(f1_dd, sol_phi1_dd[0][f1_dd])
sol_phi2_dd = smp.solve(LE4_substituted, f2_dd, dict=True, simplify=False, rational=False)

# 合并所有符号解
combined_solutions = {**sol_the1_dd[0], **sol_the2_dd[0], **sol_phi1_dd[0], **sol_phi2_dd[0]}

# 输出各角加速度表达式
print("Espressione per the1_dd:", combined_solutions[O1_dd])
print("Espressione per the2_dd:", combined_solutions[O2_dd])
print("Espressione per phi1_dd:", combined_solutions[f1_dd])
print("Espressione per phi2_dd:", combined_solutions[f2_dd])

输出的the1_dd表达式（可见残留导数项）：

Espressione per the1_dd: -gm1sin(\theta_1(t))/(m1r1 + m2r1) - gm2sin(\theta_1(t))/(m1r1 + m2r1) + m1r1sin(\theta_1(t))*cos(\theta_1(t))*Derivative(\phi_1(t), t)**2

未生效的导数替换代码 #

尝试通过以下代码替换导数项，但操作后表达式仍存在导数：

from sympy import Derivative, Function,symbols
# 定义符号函数变量
t = symbols('t')
phi1 = Function('phi_1')(t)
phi2 = Function('phi_2')(t)
theta1 = Function('theta_1')(t)
theta2 = Function('theta_2')(t)

# 定义导数符号
phi1_d = Derivative(phi1, t)
phi2_dd = Derivative(phi2, t, t)
theta2_d = Derivative(theta2, t)
theta2_dd = Derivative(theta2, t, t)

# 第一次替换导数项
simplified_the1_dd = combined_solutions[O1_dd].subs({
    phi1_d: f1_d, 
    phi2_dd: f2_dd,
    theta2_d: O2_d,
    theta2_dd: O2_dd
})

# 简化表达式
simplified_the1_dd = smp.simplify(simplified_the1_dd)
# 第二次替换导数项
simplified_the1_dd_substituted = simplified_the1_dd.subs({
    Derivative(phi1, t): f1_d,
    Derivative(phi2, t, t): f2_dd,
    Derivative(theta2, t): O2_d,
    Derivative(theta2, t, t): O2_dd
})

# 后续简化操作

预期目标 #

将符号表达式中的导数项完全替换为自定义变量，通过sympy.lambdify转换为可用于odeint的数值求解函数，完成3D双摆的运动数值模拟。

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