使用拉格朗日方法的滚动圆盘

在这里，滚动的圆盘从接触点向上形成，消除了引入广义速度的需要。只需要3个配置变量和3个速度变量来描述这个系统，以及圆盘的质量和半径，以及局部重力。:

>>> from sympy import symbols, cos, sin
>>> 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')

运动学由一系列简单的旋转构成。每个简单旋转创建一个新的坐标系，下一个旋转由新坐标系的基向量定义。此示例使用3-1-2系列的旋转，或Z、X、Y系列的旋转。角速度使用第二个坐标系的基（倾斜坐标系）来定义。:

>>> 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中创建一个没有速度的点；这是圆盘与地面接触的点。接下来，我们形成从接触点到圆盘质心的位置向量。最后，我们形成圆盘的速度和加速度。

>>> C = Point('C')
>>> C.set_vel(N, 0)
>>> Dmc = C.locatenew('Dmc', r * L.z)
>>> Dmc.v2pt_theory(C, N, R)
r*(sin(q2)*q1' + q3')*L.x - r*q2'*L.y

形成惯性偶极子。:

>>> I = inertia(L, m / 4 * r**2, m / 2 * r**2, m / 4 * r**2)
>>> mprint(I)
m*r**2/4*(L.x|L.x) + m*r**2/2*(L.y|L.y) + m*r**2/4*(L.z|L.z)
>>> BodyD = RigidBody('BodyD', Dmc, R, m, (I, Dmc))

然后我们设定势能并确定滚动圆盘的拉格朗日量。

>>> BodyD.potential_energy = - m * g * r * cos(q2)
>>> Lag = Lagrangian(N, BodyD)

然后通过初始化 LagrangesMethod 对象来生成运动方程。最后，我们使用 rhs 方法求解广义加速度（q 双点）。:

>>> q = [q1, q2, q3]
>>> l = LagrangesMethod(Lag, q)
>>> le = l.form_lagranges_equations()
>>> le.simplify(); le
Matrix([
[m*r**2*(6*sin(q2)*q3'' + 5*sin(2*q2)*q1'*q2' + 6*cos(q2)*q2'*q3' - 5*cos(2*q2)*q1''/2 + 7*q1''/2)/4],
[                      m*r*(4*g*sin(q2) - 5*r*sin(2*q2)*q1'**2/2 - 6*r*cos(q2)*q1'*q3' + 5*r*q2'')/4],
[                                                 3*m*r**2*(sin(q2)*q1'' + cos(q2)*q1'*q2' + q3'')/2]])
>>> lrhs = l.rhs(); lrhs.simplify(); lrhs
Matrix([
[                                                          q1'],
[                                                          q2'],
[                                                          q3'],
[                       -2*(2*tan(q2)*q1' + 3*q3'/cos(q2))*q2'],
[-4*g*sin(q2)/(5*r) + sin(2*q2)*q1'**2/2 + 6*cos(q2)*q1'*q3'/5],
[         (-5*cos(q2)*q1' + 6*tan(q2)*q3' + 4*q1'/cos(q2))*q2']])