执行此程序之后，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

结果看上去挺复杂，其实只要运用如下的三角公式，就和前面的结果一致了：