SymPy中dynamicsymbols如何设置二阶导数为零计算点P惯性加速度

问题背景 #

N为惯性参考系，伸缩臂沿N.x方向以恒定速率伸出，圆盘以恒定角速度$\dot{\theta}$自旋，需求解点P的惯性加速度表达式。
现有实现的痛点：无法直接给dynamicsymbols添加L、θ二阶导数为0的约束，只能手动删除表达式中对应的二阶导数项再化简，操作繁琐。

初始实现代码 #

import sympy as sym
sym.init_printing(use_latex='mathjax')
from sympy.physics.vector import *
init_vprinting(use_latex='mathjax')

N =  ReferenceFrame('N')
E =  ReferenceFrame('E')

th = dynamicsymbols('theta')
L = dynamicsymbols('L')
t = sym.symbols('t')
r = sym.symbols('r')
p = sym.symbols('p')

E.orient_axis(N, N.z, th)

p = L * N.x + r*E.x
p_vel = p.diff(t,N)
p_acc = p_vel.diff(t,N)
p_acc.simplify()

约束添加方法 #

dynamicsymbols本身没有内置约束定义接口，直接在通用求导结果中用subs()方法替换对应约束项即可，不需要手动修改表达式：

• 恒定伸出速率对应约束：L的二阶导数$\ddot{L}=0$，对应符号为L.diff(t,2)
• 恒定自旋角速度对应约束：θ的二阶导数$\ddot{\theta}=0$，对应符号为th.diff(t,2)

在原有代码末尾添加如下替换和化简代码即可得到正确结果：

p_acc_constrained = p_acc.subs({
    L.diff(t, 2): 0,
    th.diff(t, 2): 0
}).simplify()

最终得到的惯性加速度分量符合物理预期：因为L匀速伸出无轴向加速度，最终P点加速度为圆盘自旋产生的向心加速度，在N系x、y方向的分量分别为$-r\dot{\theta}^{2\cos\theta$、$-r\dot{\theta}}2\sin\theta$。

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