执行此程序之后,eq1对应于的拉格朗日方程,eq2对应于的拉格朗日方程。
由于SymPy只能对符号变量求导数,即只能计算D(L, t),而不能计算D(f, v(t))。 ?因此在求偏导数之前,将偏导数变量置换为一个tmp变量,然后对tmp变量求导数,例如下面的程序对D(v(t), t)求偏导数,即计算:
L.subs(D(v(t), t), tmp).diff(tmp).subs(tmp, D(v(t), t))
而在计算时,需要将v(t)替换为v之后再进行微分计算。由于将v(t)替换为v的同时,会将D(v(t), t)中的也进行替换,这不是我们希望的结果,因此先将D(v(t), t)替换为tmp,微分计算完毕之后再替换回去:
L.subs(D(v(t), t), tmp).subs(v(t), v).diff(v).subs(v, v(t)).subs(tmp, D(v(t), t))
最后得到eq1和eq2的值为:
>>> eq1
ddth1*(m1*l1**2 + m2*l1**2) +
ddth2*(l1*l2*m2*cos(th1)*cos(th2) + l1*l2*m2*sin(th1)*sin(th2)) +
dth2**2*(l1*l2*m2*cos(th2)*sin(th1) - l1*l2*m2*cos(th1)*sin(th2)) +
g*l1*m1*sin(th1) + g*l1*m2*sin(th1)
>>> eq2
ddth1*(l1*l2*m2*cos(th1)*cos(th2) + l1*l2*m2*sin(th1)*sin(th2)) +
dth1**2*(l1*l2*m2*cos(th1)*sin(th2) - l1*l2*m2*cos(th2)*sin(th1)) +
g*l2*m2*sin(th2) + ddth2*m2*l2**2
结果看上去挺复杂,其实只要运用如下的三角公式,就和前面的结果一致了: