拉格朗日方法求解转盘

本文是对官方教程拉格朗日方法解转盘中代码的解读。

创建坐标系

from sympy import symbols, cos, sin, print_latex
from sympy.physics.mechanics import *
mechanics_printing(pretty_print=False)
q1, q2, q3 = dynamicsymbols('q1 q2 q3')
q1d, q2d, q3d = dynamicsymbols('q1 q2 q3', 1)
r, m, g = symbols('r m g')

【dynamicsymbols】自动生成的符号是时间的函数，其第二个参数为对时间的导数阶数。故而上述代码中，$q_i$表示位置，$q_i^d$可表示速度。

N = ReferenceFrame('N')
Y = N.orientnew('Y', 'Axis', [q1, N.z])
L = Y.orientnew('L', 'Axis', [q2, Y.x])
R = L.orientnew('R', 'Axis', [q3, L.y])

其中，$N$是参考坐标系，随后用【orientnew】方法创建了三个关联坐标系，其中$Y$绕$N$的$z$轴旋转角度$q_1$；$L$绕$Y$的$x$轴再旋转$q_2$角度，将用于描述刚体的惯性并矢；$R$绕$L$的$y$轴旋转角度$q_3$，将作为刚体内部的坐标框架。

创建刚体

C = Point('C')
C.set_vel(N, 0)
D = C.locatenew('D', r * L.z)
D.v2pt_theory(C, N, R)
I = inertia(L, m / 4 * r**2, m / 2 * r**2, m / 4 * r**2)
print_latex(I)
BodyD = RigidBody('BodyD', D, R, m, (I, D))
BodyD.potential_energy = - m * g * r * cos(q2)
Lag = Lagrangian(N, BodyD)

$C$是一个点，在坐标系$N$中的速度为0。

以$C$开始，沿着$L$坐标系的$z$轴方向移动距离$r$，定义了点$D$。接下来，通过【v2pt_theory】方法，计算了$D$在$N$中的速度。

【itertial】用于创建惯性并矢，其值为$\frac{m{r}^2}{4}{\hat{\mathbf{l}}}_{\mathbf{x}}\otimes {\hat{\mathbf{l}}}_{\mathbf{x}}+\frac{m{r}^2}{2}{\hat{\mathbf{l}}}_{\mathbf{y}}\otimes {\hat{\mathbf{l}}}_{\mathbf{y}}+\frac{m{r}^2}{4}{\hat{\mathbf{l}}}_{\mathbf{z}}\otimes {\hat{\mathbf{l}}}_{\mathbf{z}}$。

【bodyD】即刚体，位于坐标系$R$中，与$N$在三个坐标轴方向分别有$q_1,q_2,q_3$的夹角。势能为$-mgr\cos (q_2)$。【Lag】为其拉格朗日量。

拉格朗日方法

q = [q1, q2, q3]
l = LagrangesMethod(Lag, q)
le = l.form_lagranges_equations()
le.simplify()
print_latex(le)
lrhs = l.rhs()
lrhs.simplify()
print_latex(lrhs)

【LagrangesMethod】可以将拉格朗日量和广义坐标为参数，构建一个拉格朗日方法对象。该对象的【form_lagranges_equations】方法，可构造出一个规范的拉格朗日方程。

$$\left[\begin{array}{c}\frac{m{r}^2\left(12\sin \left(q_2\left(t\right)\right)\frac{d^2}{d{t}^2}{q}_3\left(t\right)+10\sin \left(2q_2\left(t\right)\right)\frac{d}{dt}{q}_1\left(t\right)\frac{d}{dt}{q}_2\left(t\right)+12\cos \left(q_2\left(t\right)\right)\frac{d}{dt}{q}_2\left(t\right)\frac{d}{dt}{q}_3\left(t\right)-5\cos \left(2q_2\left(t\right)\right)\frac{d^2}{d{t}^2}{q}_1\left(t\right)+7\frac{d^2}{d{t}^2}{q}_1\left(t\right)\right)}{8}\\ \frac{mr\left(8g\sin \left(q_2\left(t\right)\right)-5r\sin \left(2q_2\left(t\right)\right){\left(\frac{d}{dt}{q}_1\left(t\right)\right)}^2-12r\cos \left(q_2\left(t\right)\right)\frac{d}{dt}{q}_1\left(t\right)\frac{d}{dt}{q}_3\left(t\right)+10r\frac{d^2}{d{t}^2}{q}_2\left(t\right)\right)}{8}\\ \frac{3m{r}^2\left(\sin \left(q_2\left(t\right)\right)\frac{d^2}{d{t}^2}{q}_1\left(t\right)+\cos \left(q_2\left(t\right)\right)\frac{d}{dt}{q}_1\left(t\right)\frac{d}{dt}{q}_2\left(t\right)+\frac{d^2}{d{t}^2}{q}_3\left(t\right)\right)}{2}\end{array}\right]$$

【lrhs】为拉格朗日方程的右侧

$$\left[\begin{array}{c}\frac{d}{dt}{q}_1\left(t\right)\\ \frac{d}{dt}{q}_2\left(t\right)\\ \frac{d}{dt}{q}_3\left(t\right)\\ -2\left(2\tan \left(q_2\left(t\right)\right)\frac{d}{dt}{q}_1\left(t\right)+\frac{3\frac{d}{dt}{q}_3\left(t\right)}{\cos \left(q_2\left(t\right)\right)}\right)\frac{d}{dt}{q}_2\left(t\right)\\ -\frac{4g\sin \left(q_2\left(t\right)\right)}{5r}+\frac{\sin \left(2q_2\left(t\right)\right){\left(\frac{d}{dt}{q}_1\left(t\right)\right)}^2}{2}+\frac{6\cos \left(q_2\left(t\right)\right)\frac{d}{dt}{q}_1\left(t\right)\frac{d}{dt}{q}_3\left(t\right)}{5}\\ \left(-5\cos \left(q_2\left(t\right)\right)\frac{d}{dt}{q}_1\left(t\right)+6\tan \left(q_2\left(t\right)\right)\frac{d}{dt}{q}_3\left(t\right)+\frac{4\frac{d}{dt}{q}_1\left(t\right)}{\cos \left(q_2\left(t\right)\right)}\right)\frac{d}{dt}{q}_2\left(t\right)\end{array}\right]$$