非最小坐标摆#
在本例中，我们演示了中提供的功能的用法 sympy.physics.mechanics 用非最小坐标系推导摆的运动方程。由于摆锤是一个单自由度系统，可以用一个坐标和一个速度（分别是摆角和角速度）来描述。而是选择使用 \(x\) 和 \(y\) 质量坐标导致需要约束。系统显示如下：
系统将使用凯恩和拉格朗日方法建模，并将得到的EOM线性化。虽然这是一个简单的问题，但它应该说明在存在约束的情况下线性化方法的使用。
凯恩方法#
首先我们需要创建 dynamicsymbols 需要如上图所示描述系统。在这种情况下，广义坐标 \(q_1\) 和 \(q_2\) 代表质量 \(x\) 和 \(y\) 惯性坐标系 \(N\) 框架。同样，广义速度 \(u_1\) 和 \(u_2\) 表示这些方向上的速度。我们也创造了一些 symbols 表示摆锤的长度和质量，以及重力和时间。:
>>> from sympy.physics.mechanics import *
>>> from sympy import symbols, atan, Matrix, solve
>>> # Create generalized coordinates and speeds for this non-minimal realization
>>> # q1, q2 = N.x and N.y coordinates of pendulum
>>> # u1, u2 = N.x and N.y velocities of pendulum
>>> q1, q2 = dynamicsymbols('q1:3')
>>> q1d, q2d = dynamicsymbols('q1:3', level=1)
>>> u1, u2 = dynamicsymbols('u1:3')
>>> u1d, u2d = dynamicsymbols('u1:3', level=1)
>>> L, m, g, t = symbols('L, m, g, t')
接下来，我们创建一个世界坐标系 \(N\) ，及其原点 \(N^*\) . 原点的速度设定为0。第二个坐标系 \(A\) 其方向应使其x轴沿着摆锤（如上图所示）。:
>>> # Compose world frame
>>> N = ReferenceFrame('N')
>>> pN = Point('N*')
>>> pN.set_vel(N, 0)
>>> # A.x is along the pendulum
>>> theta1 = atan(q2/q1)
>>> A = N.orientnew('A', 'axis', [theta1, N.z])
确定摆锤质量的位置很容易，只要用它在世界坐标系中的x和y坐标来指定它的位置。A Particle 然后创建对象来表示该位置的质量。:
>>> # Locate the pendulum mass
>>> P = pN.locatenew('P1', q1*N.x + q2*N.y)
>>> pP = Particle('pP', P, m)
运动微分方程将广义坐标的导数与广义速度联系起来。在这个例子中，速度是导数，所以很简单。还创建了一个字典来映射 \(\dot{{q}}\) 到 \(u\) ：：
>>> # Calculate the kinematic differential equations
>>> kde = Matrix([q1d - u1,
... q2d - u2])
>>> dq_dict = solve(kde, [q1d, q2d])
质量的速度就是这个位置相对于原点的时间导数 \(N^*\) ：：
>>> # Set velocity of point P
>>> P.set_vel(N, P.pos_from(pN).dt(N).subs(dq_dict))
由于该系统的坐标比自由度多，因此需要约束。配置约束将坐标相互关联。在这种情况下，约束条件是从原点到质量的距离始终是长度 \(L\) （钟摆不会变长）。同样，速度约束是 A.x 方向始终为0（无径向速度）。:
>>> f_c = Matrix([P.pos_from(pN).magnitude() - L])
>>> f_v = Matrix([P.vel(N).express(A).dot(A.x)])
>>> f_v.simplify()
系统所受的力就是重力 P . ：：
>>> # Input the force resultant at P
>>> R = m*g*N.x
在问题设置中，可以使用 KanesMethod 班级。由于存在约束，因此需要向类提供从属坐标和独立坐标。在这种情况下，我们将使用 \(q_2\) 和 \(u_2\) 作为独立的坐标：
>>> # Derive the equations of motion using the KanesMethod class.
>>> KM = KanesMethod(N, q_ind=[q2], u_ind=[u2], q_dependent=[q1],
... u_dependent=[u1], configuration_constraints=f_c,
... velocity_constraints=f_v, kd_eqs=kde)
>>> (fr, frstar) = KM.kanes_equations([pP],[(P, R)])
为了线性化，可以在调用时指定操作点，或者在以后替换。在本例中，我们将在调用中提供它们，以列表形式提供。这个 A_and_B=True kwarg表示要求解 \(M\) 矩阵与显式线性化的求解 \(A\) 和 \(B\) 矩阵。这个 simplify=True kwarg表示在线性化调用中进行简化，并返回预简化的矩阵。对于简单的系统来说，这样做的成本很小，但是对于较大的系统来说，这可能是一个昂贵的操作，应该避免。:
>>> # Set the operating point to be straight down, and non-moving
>>> q_op = {q1: L, q2: 0}
>>> u_op = {u1: 0, u2: 0}
>>> ud_op = {u1d: 0, u2d: 0}
>>> # Perform the linearization
>>> A, B, inp_vec = KM.linearize(op_point=[q_op, u_op, ud_op], A_and_B=True,
... new_method=True, simplify=True)
>>> A
Matrix([
[ 0, 1],
[-g/L, 0]])
>>> B
Matrix(0, 0, [])
结果 \(A\) 矩阵的维数为2x2，而总状态数为 len(q) + len(u) = 2 + 2 = 4 . 这是因为对于受约束的系统 A_and_B 表单有一个只包含独立坐标和速度的分区状态向量。在数学上，关于这一点线性化的系统可以写成：
拉格朗日方法#
拉格朗日方法的推导与上述凯恩方法的推导非常相似。和以前一样，我们首先创建 dynamicsymbols 需要描述系统。在这种情况下，广义坐标 \(q_1\) 和 \(q_2\) 代表质量 \(x\) 和 \(y\) 惯性坐标系 \(N\) 框架。这就产生了时间导数 \(\dot{{q_1}}\) 和 \(\dot{{q_2}}\) 表示这些方向上的速度。我们也创造了一些 symbols 表示摆锤的长度和质量，以及重力和时间。:
>>> from sympy.physics.mechanics import *
>>> from sympy import symbols, atan, Matrix
>>> q1, q2 = dynamicsymbols('q1:3')
>>> q1d, q2d = dynamicsymbols('q1:3', level=1)
>>> L, m, g, t = symbols('L, m, g, t')
接下来，我们创建一个世界坐标系 \(N\) ，及其原点 \(N^*\) . 原点的速度设定为0。第二个坐标系 \(A\) 其方向应使其x轴沿着摆锤（如上图所示）。:
>>> # Compose World Frame
>>> N = ReferenceFrame('N')
>>> pN = Point('N*')
>>> pN.set_vel(N, 0)
>>> # A.x is along the pendulum
>>> theta1 = atan(q2/q1)
>>> A = N.orientnew('A', 'axis', [theta1, N.z])
确定摆锤质量的位置很容易，只要用它在世界坐标系中的x和y坐标来指定它的位置。A Particle 然后创建对象来表示该位置的质量。:
>>> # Create point P, the pendulum mass
>>> P = pN.locatenew('P1', q1*N.x + q2*N.y)
>>> P.set_vel(N, P.pos_from(pN).dt(N))
>>> pP = Particle('pP', P, m)
由于该系统的坐标比自由度多，因此需要约束。在这种情况下，只需要一个完整约束：从原点到质量的距离总是长度 \(L\) （钟摆不会变长）。:
>>> # Holonomic Constraint Equations
>>> f_c = Matrix([q1**2 + q2**2 - L**2])
系统所受的力就是重力 P . ：：
>>> # Input the force resultant at P
>>> R = m*g*N.x
通过问题设置，可以计算拉格朗日，并形成运动方程。请注意 LagrangesMethod 包括拉格朗日、广义坐标、约束（由 hol_coneqs 或 nonhol_coneqs )，列表（物体，力）和惯性系。与 KanesMethod 初始值设定项、独立坐标和从属坐标在 LagrangesMethod 对象。稍后将提供这样的分区。:
>>> # Calculate the lagrangian, and form the equations of motion
>>> Lag = Lagrangian(N, pP)
>>> LM = LagrangesMethod(Lag, [q1, q2], hol_coneqs=f_c, forcelist=[(P, R)], frame=N)
>>> lag_eqs = LM.form_lagranges_equations()
接下来，我们编写了操作点字典，设置在悬挂静止位置：：：
>>> # Compose operating point
>>> op_point = {q1: L, q2: 0, q1d: 0, q2d: 0, q1d.diff(t): 0, q2d.diff(t): 0}
由于公式中有约束条件，会有相应的拉格朗日乘子。这些也可能出现在线性化形式中，因此也应包括在操作点字典中。幸运的是 LagrangesMethod 类提供了一种简单的方法来求解给定操作点处的乘法器 solve_multipliers 方法。：：
>>> # Solve for multiplier operating point
>>> lam_op = LM.solve_multipliers(op_point=op_point)
利用这个解，可以完成线性化。请注意，与 KanesMethod 接近，那个 LagrangesMethod.linearize 方法还要求将广义坐标及其时间导数分解为独立和相依向量。这和传入 KanesMethod 以上施工单位：
>>> op_point.update(lam_op)
>>> # Perform the Linearization
>>> A, B, inp_vec = LM.linearize([q2], [q2d], [q1], [q1d],
... op_point=op_point, A_and_B=True)
>>> A
Matrix([
[ 0, 1],
[-g/L, 0]])
>>> B
Matrix(0, 0, [])
结果 \(A\) 矩阵的维数为2x2，而总状态数为 2*len(q) = 4 . 这是因为对于受约束的系统 A_and_B form有一个只包含独立坐标及其导数的分块状态向量。在数学上，关于这一点线性化的系统可以写成：