SymPy 1.13 中文文档（三十三）

原文：docs.sympy.org/latest/index.html

物理/力学中的线性化

原文链接：docs.sympy.org/latest/modules/physics/mechanics/linearize.html

sympy.physics.mechanics 包括了关于操作点（也称为平衡条件）附近生成的运动方程（EOM）的线性化方法。注意，这个操作点不必是一个平衡位置，它只需满足运动方程。

通过对绕操作点的 EOM 进行一阶泰勒展开来完成线性化。当没有依赖坐标或速度时，这只是关于 (q) 和 (u) 的右手边雅可比矩阵。然而，在存在约束条件时，需要更加小心。此处提供的线性化方法可以正确处理这些约束。

背景

在 sympy.physics.mechanics 中，我们假设所有系统可以用以下一般形式表示：

[\begin{split}f_{c}(q, t) &= 0_{l \times 1}\ f_{v}(q, u, t) &= 0_{m \times 1}\ f_{a}(q, \dot{q}, u, \dot{u}, t) &= 0_{m \times 1}\ f_{0}(q, \dot{q}, t) + f_{1}(q, u, t) &= 0_{n \times 1}\ f_{2}(q, u, \dot{u}, t) + f_{3}(q, \dot{q}, u, r, t) + f_{4}(q, \lambda, t) &= 0_{(o-m+k) \times 1}\end{split}]

其中

[\begin{split}q, \dot{q} & \in \mathbb{R}^n\ u, \dot{u} & \in \mathbb{R}^o\ r & \in \mathbb{R}^s\ \lambda & \in \mathbb{R}^k\end{split}]

在这种形式中，

• (f_{c}) 表示配置约束方程。

• (f_{v}) 表示速度约束方程。

• (f_{a}) 表示加速度约束方程。

• (f_{0}) 和 (f_{1}) 形成了运动学微分方程。

• (f_{2}), (f_{3}) 和 (f_{4}) 形成了动态微分方程。

• (q) 和 (\dot{q}) 是广义坐标及其导数。

• (u) 和 (\dot{u}) 是广义速度及其导数。

• (r) 是系统输入。

• (\lambda) 是拉格朗日乘子。

这种广义形式包含在 Linearizer 类中，它执行实际的线性化。KanesMethod 和 LagrangesMethod 对象都有使用 to_linearizer 类方法形成线性化器的方法。

一旦系统被强制转换为广义形式，可以解出线性化的 EOM。sympy.physics.mechanics 中提供的方法允许两种不同形式的线性化 EOM：

(M), (A), 和 (B)

在这种形式中，强制矩阵被线性化为两个独立的矩阵 (A) 和 (B)。这是线性化 EOM 的默认形式。得到的方程为：

[\begin{split}M \begin{bmatrix} \delta \dot{q} \ \delta \dot{u} \ \delta \lambda \end{bmatrix} = A \begin{bmatrix} \delta q_i \ \delta u_i \end{bmatrix} + B \begin{bmatrix} \delta r \end{bmatrix}\end{split}]

其中

[\begin{split}M &\in \mathbb{R}^{(n+o+k) \times (n+o+k)}\ A &\in \mathbb{R}^{(n+o+k) \times (n-l+o-m)}\ B &\in \mathbb{R}^{(n+o+k) \times s}\end{split}]

注意(q_i)和(u_i)只是独立坐标和速度，而(q)和(u)包含独立和依赖坐标和速度。

(A)和(B)

在此形式中，线性化的运动方程被带入明确的一阶形式，仅以独立坐标和速度来表示。此形式经常用于稳定性分析或控制理论。得到的方程如下：

[\begin{split}\begin{bmatrix} \delta \dot{q_i} \ \delta \dot{u_i} \end{bmatrix} = A \begin{bmatrix} \delta q_i \ \delta u_i \end{bmatrix} + B \begin{bmatrix} \delta r \end{bmatrix}\end{split}]

其中

[\begin{split}A &\in \mathbb{R}^{(n-l+o-m) \times (n-l+o-m)}\ B &\in \mathbb{R}^{(n-l+o-m) \times s}\end{split}]

要使用此形式，请在linearize类方法中设置A_and_B=True。

线性化 Kane 方程

初始化KanesMethod对象并使用kanes_equations类方法形成(F_r)和(F_r^*)后，可以通过几种方法完成线性化。不同的方法将在简单摆系统中进行演示：

>>> from sympy import symbols, Matrix
>>> from sympy.physics.mechanics import *
>>> q1 = dynamicsymbols('q1')                     # Angle of pendulum
>>> u1 = dynamicsymbols('u1')                     # Angular velocity
>>> q1d = dynamicsymbols('q1', 1)
>>> L, m, t, g = symbols('L, m, t, g')

>>> # Compose world frame
>>> N = ReferenceFrame('N')
>>> pN = Point('N*')
>>> pN.set_vel(N, 0)

>>> # A.x is along the pendulum
>>> A = N.orientnew('A', 'axis', [q1, N.z])
>>> A.set_ang_vel(N, u1*N.z)

>>> # Locate point P relative to the origin N*
>>> P = pN.locatenew('P', L*A.x)
>>> vel_P = P.v2pt_theory(pN, N, A)
>>> pP = Particle('pP', P, m)

>>> # Create Kinematic Differential Equations
>>> kde = Matrix([q1d - u1])

>>> # Input the force resultant at P
>>> R = m*g*N.x

>>> # Solve for eom with kanes method
>>> KM = KanesMethod(N, q_ind=[q1], u_ind=[u1], kd_eqs=kde)
>>> fr, frstar = KM.kanes_equations([pP], [(P, R)]) 

1. 直接使用Linearizer类：

可以使用to_linearizer类方法创建一个linearizer对象。这将强制KanesMethod对象中找到的表示形式转换为上述的广义形式。由于在创建KanesMethod对象时已经指定了独立和依赖坐标和速度，因此无需在此处指定它们。

>>> linearizer = KM.to_linearizer() 

然后可以使用Linearizer对象的linearize方法形成线性化的运动方程：

>>> M, A, B = linearizer.linearize()
>>> M
Matrix([
[1,       0],
[0, -L**2*m]])
>>> A
Matrix([
[                 0, 1],
[L*g*m*cos(q1(t)), 0]])
>>> B
Matrix(0, 0, []) 

或者，可以通过指定A_and_B=True来生成(A)和(B)形式：

>>> A, B = linearizer.linearize(A_and_B=True)
>>> A
Matrix([
[                0, 1],
[-g*cos(q1(t))/L, 0]])
>>> B
Matrix(0, 0, []) 

还可以将操作点指定为字典或字典集合。这将在返回矩阵之前在指定点评估线性化形式：

>>> op_point = {q1: 0, u1: 0}
>>> A_op, B_op = linearizer.linearize(A_and_B=True, op_point=op_point)
>>> A_op
Matrix([
[     0, 1],
[-g/L, 0]]) 

注意，通过在生成的矩阵上应用msubs而不使用op_point关键字也可以达到相同的效果：

>>> assert msubs(A, op_point) == A_op 

有时返回的矩阵可能不是最简化的形式。可以在事后执行简化，或者通过将simplify关键字设置为True使Linearizer对象在内部执行简化。

2. 使用linearize类方法：

KanesMethod类的linearize方法作为一个便捷的包装器提供，内部调用to_linearizer，执行线性化并返回结果。请注意，linearize方法中描述的所有关键字参数在此处也是可用的：

>>> A, B, inp_vec = KM.linearize(A_and_B=True, op_point=op_point, new_method=True)
>>> A
Matrix([
[     0, 1],
[-g/L, 0]]) 

附加输出inp_vec是一个包含所有未包含在广义坐标或速度向量中的dynamicsymbols的向量。这些被假定为系统的输入，形成了上述背景中描述的向量(r)。在这个例子中没有输入，因此向量为空：

>>> inp_vec
Matrix(0, 0, []) 

线性化 Lagrange 方程

Lagrange 方程的线性化过程与 Kane 方程的线性化过程基本相同。与之前一样，该过程将通过一个简单的摆系统进行演示：

>>> # Redefine A and P in terms of q1d, not u1
>>> A = N.orientnew('A', 'axis', [q1, N.z])
>>> A.set_ang_vel(N, q1d*N.z)
>>> P = pN.locatenew('P', L*A.x)
>>> vel_P = P.v2pt_theory(pN, N, A)
>>> pP = Particle('pP', P, m)

>>> # Solve for eom with Lagrange's method
>>> Lag = Lagrangian(N, pP)
>>> LM = LagrangesMethod(Lag, [q1], forcelist=[(P, R)], frame=N)
>>> lag_eqs = LM.form_lagranges_equations() 

1. 直接使用Linearizer类：

可以使用to_linearizer类方法从LagrangesMethod对象形成一个Linearizer对象。这个过程与KanesMethod类的过程唯一的区别在于，LagrangesMethod对象在内部没有指定其独立和依赖坐标和速度。这些必须在调用to_linearizer时指定。在这个例子中没有依赖的坐标和速度，但如果有的话，它们将包含在q_dep和qd_dep关键字参数中：

>>> linearizer = LM.to_linearizer(q_ind=[q1], qd_ind=[q1d]) 

一旦进入这种形式，一切就像之前的KanesMethod示例一样：

>>> A, B = linearizer.linearize(A_and_B=True, op_point=op_point)
>>> A
Matrix([
[     0, 1],
[-g/L, 0]]) 

2. 使用linearize类方法：

与KanesMethod类似，LagrangesMethod类还提供了一个linearize方法作为一个方便的包装器，内部调用to_linearizer，执行线性化并返回结果。与之前一样，唯一的区别是必须在调用中指定独立和依赖的坐标和速度：

>>> A, B, inp_vec = LM.linearize(q_ind=[q1], qd_ind=[q1d], A_and_B=True, op_point=op_point)
>>> A
Matrix([
[     0, 1],
[-g/L, 0]]) 

潜在问题

虽然Linearizer类应该能够线性化所有系统，但可能会出现一些潜在问题。下面讨论这些问题，并提供一些解决方法的故障排除提示。

1. 使用A_and_B=True的符号线性化速度较慢

这可能由多种原因引起，但最有可能的原因是，符号求解一个大型线性系统是一项昂贵的操作。指定一个操作点将减小表达式的大小并加快速度。然而，如果需要纯符号解（例如稍后应用多个操作点的应用），可以通过在A_and_B=False下评估，然后在应用操作点后手动求解来解决这个问题：

>>> M, A, B = linearizer.linearize()
>>> M_op = msubs(M, op_point)
>>> A_op = msubs(A, op_point)
>>> perm_mat = linearizer.perm_mat
>>> A_lin = perm_mat.T * M_op.LUsolve(A_op)
>>> A_lin
Matrix([
[     0, 1],
[-g/L, 0]]) 

在求解前，A和M中的符号越少，解决方案的速度就越快。因此，对于大型表达式，延迟到大多数符号被替换为它们的数值之后再转换为(A)和(B)形式可能会对你有利。

2. 线性化形式的矩阵元素为nan、zoo或oo

这有两个潜在原因。首先（也是你应该首先检查的原因）是，在某些操作点上，某些依赖坐标的选择会导致奇点。系统性地进行坐标分区以避免这种情况超出了本指南的范围；详细信息请参见[Blajer1994]。

另一个潜在的原因是在操作点替换之前，矩阵可能没有处于最简形式。这种行为的简单示例是：

>>> from sympy import sin, tan
>>> expr = sin(q1)/tan(q1)
>>> op_point = {q1: 0}
>>> expr.subs(op_point)
nan 

注意，如果在替换之前简化了这个表达式，将得到正确的值：

>>> expr.simplify().subs(op_point)
1 

目前还没有找到避免这种情况的好方法。对于大小合理的表达式，使用msubs和smart=True将应用一种算法，试图避免这些条件。但对于大型表达式，这会非常耗时。

>>> msubs(expr, op_point, smart=True)
1 

更多示例

上述摆例子很简单，但没有包含任何依赖坐标或速度。为了更详细的例子，同一个摆通过凯恩和拉格朗日方法使用依赖坐标进行了线性化：

• 非最小坐标摆

  • 凯恩方法

  • 拉格朗日方法

非极小坐标摆

原文链接：docs.sympy.org/latest/modules/physics/mechanics/examples/lin_pend_nonmin_example.html

在这个示例中，我们演示了使用 sympy.physics.mechanics 提供的功能来推导具有非极小坐标系的摆的运动方程（EOM）。由于摆是一个自由度系统，可以用一个坐标和一个速度（即摆角和角速度）描述。然而，选择使用质量的 (x) 和 (y) 坐标来描述系统则需要施加约束。系统如下所示：

image/svg+xml L m g q1 q2 u1 u2 Nx Ny Ax Ay

该系统将使用凯恩和拉格朗日方法建模，并线性化得到结果的运动方程。虽然这是一个简单的问题，但它应该能够展示在线性化方法中处理约束条件的应用。

凯恩方法

首先，我们需要创建用于描述上图所示系统的dynamicsymbols。在这种情况下，广义坐标(q_1)和(q_2)表示惯性坐标系(N)中的质量(x)和(y)坐标。同样，广义速度(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 坐标的形式指定其位置来轻松定位摆动质量。然后创建一个Particle对象来表示该位置处的质量。

>>> # Locate the pendulum mass
>>> P = pN.locatenew('P1', q1*N.x + q2*N.y)
>>> pP = Particle('pP', P, m) 

运动学微分方程（KDEs）将广义坐标的导数与广义速度相关联。在这种情况下，速度是导数，因此这些很简单。还创建了一个字典来映射(\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关键字参数指示求解反转(M)矩阵并仅解出显式线性化的(A)和(B)矩阵。simplify=True关键字参数指示在线性化调用内简化，并返回预简化的矩阵。对于简单系统而言，这样做的成本很小，但对于更大的系统来说，这可能是一项昂贵的操作，应该避免使用。

>>> # 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)矩阵的维度为 2 x 2，而总状态数为len(q) + len(u) = 2 + 2 = 4。这是因为对于约束系统，得到的A_and_B形式具有仅包含独立坐标和速度的分区状态向量。在数学上写出，围绕这一点线性化的系统将被写为：

[\begin{split}\begin{bmatrix} \dot{q_2} \ \dot{u_2} \end{bmatrix} = \begin{bmatrix} 0 & 1 \ \frac{-g}{L} & 0 \end{bmatrix} \begin{bmatrix} q_2 \ u_2 \end{bmatrix}\end{split}]

拉格朗日方法

使用拉格朗日方法进行的推导与上述使用 Kane 方法的方法非常相似。与之前一样，首先创建描述系统所需的dynamicsymbols。在本例中，广义坐标(q_1)和(q_2)表示惯性(N)框架中的质量(x)和(y)坐标。这导致时间导数(\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 坐标指定其位置，很容易找到摆质量的位置。然后创建一个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)矩阵的维度为 2 x 2，总状态数为2*len(q) = 4。这是因为对于约束系统，得到的A_and_B形式具有一个分隔的状态向量，仅包含独立坐标及其导数。在数学上写出来，关于该点线性化的系统如下所示：

[\begin{split}\begin{bmatrix} \dot{q_2} \ \ddot{q_2} \end{bmatrix} = \begin{bmatrix} 0 & 1 \ \frac{-g}{L} & 0 \end{bmatrix} \begin{bmatrix} q_2 \ \dot{q_2} \end{bmatrix}\end{split}]

物理/力学示例

原文：docs.sympy.org/latest/modules/physics/mechanics/examples.html

这里有一些例子，说明了如何典型地使用这个模块。我们按照难度递增的顺序排列了这些示例。如果您使用此模块做了其他人可能会觉得有用或有趣的事情，请考虑在这里添加！

• 滚动圆盘

• 自行车

• 非最小坐标摆

• 多自由度全纯系统

• 四连杆

一个滚动的圆盘

原文：docs.sympy.org/latest/modules/physics/mechanics/examples/rollingdisc_example.html

假设圆盘是无限薄的，在地面上仅接触 1 点，并且在地面上滚动而不打滑。请看下面的图片。

image/svg+xml R N nz nx ny rx ry rz

我们以三种不同的方式对滚动圆盘进行建模，以展示此模块的更多功能。

• 使用凯恩方法的滚动圆盘

• 使用凯恩方法和约束力的滚动圆盘

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

采用凯恩方法的滚动圆盘

原文链接：docs.sympy.org/latest/modules/physics/mechanics/examples/rollingdisc_example_kane.html

这里定义了滚动圆盘的运动学，从接触点向上形成，无需引入广义速度。仅需 3 个配置和三个速度变量来描述此系统，以及圆盘的质量和半径，以及局部重力（注意质量将被抵消）。

>>> from sympy import symbols, sin, cos, tan
>>> from sympy.physics.mechanics import *
>>> q1, q2, q3, u1, u2, u3  = dynamicsymbols('q1 q2 q3 u1 u2 u3')
>>> q1d, q2d, q3d, u1d, u2d, u3d = dynamicsymbols('q1 q2 q3 u1 u2 u3', 1)
>>> r, m, g = symbols('r m g')
>>> mechanics_printing(pretty_print=False) 

运动学是由一系列简单的旋转构成的。每次简单旋转都会创建一个新的参考系，下一次旋转的定义由新参考系的基向量决定。本例中使用的是 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])
>>> w_R_N_qd = R.ang_vel_in(N)
>>> R.set_ang_vel(N, u1 * L.x + u2 * L.y + u3 * L.z) 

这是平动运动学。我们在 N 中创建一个无速度点；这是圆盘和地面之间的接触点。接下来我们形成从接触点到圆盘质心的位置矢量。最后我们形成圆盘的速度和加速度。

>>> C = Point('C')
>>> C.set_vel(N, 0)
>>> Dmc = C.locatenew('Dmc', r * L.z)
>>> Dmc.v2pt_theory(C, N, R)
r*u2*L.x - r*u1*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) 

运动学微分方程；广义坐标的时间导数与广义速度的关系。

>>> kd = [dot(R.ang_vel_in(N) - w_R_N_qd, uv) for uv in L] 

创建力列表；这是圆盘质心处的重力。然后我们通过将一个点分配给质心属性、一个参考系分配给参考系属性、以及质量和惯性来创建圆盘。然后我们形成物体列表。

>>> ForceList = [(Dmc, - m * g * Y.z)]
>>> BodyD = RigidBody('BodyD', Dmc, R, m, (I, Dmc))
>>> BodyList = [BodyD] 

最后，我们形成运动方程，使用与之前相同的步骤。指定惯性参考系，提供广义坐标和速度，提供运动微分方程字典，计算来自力列表的 Fr 和来自物体列表的 Fr*，计算质量矩阵和强制项，然后求解广义速度的时间导数 u 点。

>>> KM = KanesMethod(N, q_ind=[q1, q2, q3], u_ind=[u1, u2, u3], kd_eqs=kd)
>>> (fr, frstar) = KM.kanes_equations(BodyList, ForceList)
>>> MM = KM.mass_matrix
>>> forcing = KM.forcing
>>> rhs = MM.inv() * forcing
>>> kdd = KM.kindiffdict()
>>> rhs = rhs.subs(kdd)
>>> rhs.simplify()
>>> mprint(rhs)
Matrix([
[(4*g*sin(q2) + 6*r*u2*u3 - r*u3**2*tan(q2))/(5*r)],
[                                       -2*u1*u3/3],
[                          (-2*u2 + u3*tan(q2))*u1]]) 

使用 Kane 方法和约束力的滚动圆盘

原文链接：docs.sympy.org/latest/modules/physics/mechanics/examples/rollingdisc_example_kane_constraints.html

现在我们将重新讨论滚动圆盘的示例，但这次我们将把非贡献力（约束力）也带入视野。详细说明请参见[Kane1985]。在这里，我们将打开进行向量操作时的自动简化。这使得小问题的输出更美观，但可能导致更大的向量操作挂起。

>>> from sympy import symbols, sin, cos, tan
>>> from sympy.physics.mechanics import *
>>> mechanics_printing(pretty_print=False)
>>> q1, q2, q3, u1, u2, u3  = dynamicsymbols('q1 q2 q3 u1 u2 u3')
>>> q1d, q2d, q3d, u1d, u2d, u3d = dynamicsymbols('q1 q2 q3 u1 u2 u3', 1)
>>> r, m, g = symbols('r m g') 

这两行引入了寻找约束力所需的额外量。

>>> u4, u5, u6, f1, f2, f3 = dynamicsymbols('u4 u5 u6 f1 f2 f3') 

大部分主要代码与之前相同。

>>> 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])
>>> w_R_N_qd = R.ang_vel_in(N)
>>> R.set_ang_vel(N, u1 * L.x + u2 * L.y + u3 * L.z) 

不滑动滚动的定义要求接触点的速度为零；在引入约束力的过程中，我们必须引入此点的速度，其定义始终为零。它们与地面垂直，沿着圆盘滚动的路径，在地面上沿着垂直方向。

>>> C = Point('C')
>>> C.set_vel(N, u4 * L.x + u5 * cross(Y.z, L.x) + u6 * Y.z)
>>> Dmc = C.locatenew('Dmc', r * L.z)
>>> vel = Dmc.v2pt_theory(C, N, R)
>>> I = inertia(L, m / 4 * r**2, m / 2 * r**2, m / 4 * r**2)
>>> kd = [dot(R.ang_vel_in(N) - w_R_N_qd, uv) for uv in L] 

就像我们之前作为这个过程的一部分引入了三个速度一样，我们也引入了三个力；它们与速度的方向相同，并代表这些方向上的约束力。

>>> ForceList = [(Dmc, - m * g * Y.z), (C, f1 * L.x + f2 * cross(Y.z, L.x) + f3 * Y.z)]
>>> BodyD = RigidBody('BodyD', Dmc, R, m, (I, Dmc))
>>> BodyList = [BodyD]

>>> KM = KanesMethod(N, q_ind=[q1, q2, q3], u_ind=[u1, u2, u3], kd_eqs=kd,
...           u_auxiliary=[u4, u5, u6])
>>> (fr, frstar) = KM.kanes_equations(BodyList, ForceList)
>>> MM = KM.mass_matrix
>>> forcing = KM.forcing
>>> rhs = MM.inv() * forcing
>>> kdd = KM.kindiffdict()
>>> rhs = rhs.subs(kdd)
>>> rhs.simplify()
>>> mprint(rhs)
Matrix([
[(4*g*sin(q2) + 6*r*u2*u3 - r*u3**2*tan(q2))/(5*r)],
[                                       -2*u1*u3/3],
[                          (-2*u2 + u3*tan(q2))*u1]])
>>> from sympy import trigsimp, signsimp, collect, factor_terms
>>> def simplify_auxiliary_eqs(w):
...     return signsimp(trigsimp(collect(collect(factor_terms(w), f2), m*r)))
>>> mprint(KM.auxiliary_eqs.applyfunc(simplify_auxiliary_eqs))
Matrix([
[                                      -m*r*(u1*u3 + u2') + f1],
[-m*r*u1**2*sin(q2) - m*r*u2*u3/cos(q2) + m*r*cos(q2)*u1' + f2],
[                -g*m + m*r*(u1**2*cos(q2) + sin(q2)*u1') + f3]]) 

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

原文链接：docs.sympy.org/latest/modules/physics/mechanics/examples/rollingdisc_example_lagrange.html

该滚动盘是从接触点向上形成的，无需引入广义速度。只需要 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 double dots）。

>>> 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']]) 

一辆自行车

原文链接：docs.sympy.org/latest/modules/physics/mechanics/examples/bicycle_example.html

自行车是一个有趣的系统，因为它有多个刚体、非完整约束和完整约束。运动的线性化动力学方程在[Meijaard2007]中有介绍。本例将详细介绍在sympy.physics.mechanics中构建运动方程。

>>> from sympy import *
>>> from sympy.physics.mechanics import *
>>> print('Calculation of Linearized Bicycle \"A\" Matrix, '
...       'with States: Roll, Steer, Roll Rate, Steer Rate')
Calculation of Linearized Bicycle "A" Matrix, with States: Roll, Steer, Roll Rate, Steer Rate 

请注意，此代码已从 Autolev 粗略移植，这也是一些不寻常命名约定的原因。目的是尽可能相似，以帮助初始移植和调试。

>>> mechanics_printing(pretty_print=False) 

坐标和速度的声明：在这段代码中，q 点的简单定义 q = u 被使用。速度包括：偏航框架角速率、滚动框架角速率、后轮框架角速率（旋转运动）、框架角速率（俯仰运动）、转向框架角速率以及前轮角速率（旋转运动）。车轮位置是可以忽略的坐标，因此没有被介绍。

>>> q1, q2, q3, q4, q5 = dynamicsymbols('q1 q2 q3 q4 q5')
>>> q1d, q2d, q4d, q5d = dynamicsymbols('q1 q2 q4 q5', 1)
>>> u1, u2, u3, u4, u5, u6 = dynamicsymbols('u1 u2 u3 u4 u5 u6')
>>> u1d, u2d, u3d, u4d, u5d, u6d = dynamicsymbols('u1 u2 u3 u4 u5 u6', 1) 

系统参数的声明：下面的符号应该相对清晰易懂。

>>> WFrad, WRrad, htangle, forkoffset = symbols('WFrad WRrad htangle forkoffset')
>>> forklength, framelength, forkcg1 = symbols('forklength framelength forkcg1')
>>> forkcg3, framecg1, framecg3, Iwr11 = symbols('forkcg3 framecg1 framecg3 Iwr11')
>>> Iwr22, Iwf11, Iwf22, Iframe11 = symbols('Iwr22 Iwf11 Iwf22 Iframe11')
>>> Iframe22, Iframe33, Iframe31, Ifork11 = \
...     symbols('Iframe22 Iframe33 Iframe31 Ifork11')
>>> Ifork22, Ifork33, Ifork31, g = symbols('Ifork22 Ifork33 Ifork31 g')
>>> mframe, mfork, mwf, mwr = symbols('mframe mfork mwf mwr') 

为系统设置参考框架：N - 惯性 Y - 偏航 R - 滚动 WR - 后轮，旋转角度是可以忽略的坐标，因此未定向。Frame - 自行车框架 TempFrame - 静态旋转框架，用于更容易地参考惯性定义 Fork - 自行车叉 TempFork - 静态旋转框架，用于更容易地参考惯性定义 WF - 前轮，同样具有可以忽略的坐标。

>>> N = ReferenceFrame('N')
>>> Y = N.orientnew('Y', 'Axis', [q1, N.z])
>>> R = Y.orientnew('R', 'Axis', [q2, Y.x])
>>> Frame = R.orientnew('Frame', 'Axis', [q4 + htangle, R.y])
>>> WR = ReferenceFrame('WR')
>>> TempFrame = Frame.orientnew('TempFrame', 'Axis', [-htangle, Frame.y])
>>> Fork = Frame.orientnew('Fork', 'Axis', [q5, Frame.x])
>>> TempFork = Fork.orientnew('TempFork', 'Axis', [-htangle, Fork.y])
>>> WF = ReferenceFrame('WF') 

自行车的运动学：第一段代码形成相关点的位置，后轮接触点 -> 后轮质心 -> 框架质心 + 框架/叉连接 -> 叉子质心 + 前轮质心 -> 前轮接触点。

>>> WR_cont = Point('WR_cont')
>>> WR_mc = WR_cont.locatenew('WR_mc', WRrad * R.z)
>>> Steer = WR_mc.locatenew('Steer', framelength * Frame.z)
>>> Frame_mc = WR_mc.locatenew('Frame_mc', -framecg1 * Frame.x + framecg3 * Frame.z)
>>> Fork_mc = Steer.locatenew('Fork_mc', -forkcg1 * Fork.x + forkcg3 * Fork.z)
>>> WF_mc = Steer.locatenew('WF_mc', forklength * Fork.x + forkoffset * Fork.z)
>>> WF_cont = WF_mc.locatenew('WF_cont', WFrad*(dot(Fork.y, Y.z)*Fork.y - \
...                                             Y.z).normalize()) 

设置每个框架的角速度：当首次需要时，角加速度将自动通过不同角速度的微分计算。:: u1 是偏航速率 u2 是滚动速率 u3 是后轮速率 u4 是框架俯仰速率 u5 是叉子转向速率 u6 是前轮速率

>>> Y.set_ang_vel(N, u1 * Y.z)
>>> R.set_ang_vel(Y, u2 * R.x)
>>> WR.set_ang_vel(Frame, u3 * Frame.y)
>>> Frame.set_ang_vel(R, u4 * Frame.y)
>>> Fork.set_ang_vel(Frame, u5 * Fork.x)
>>> WF.set_ang_vel(Fork, u6 * Fork.y) 

形成点的速度，使用两点定理。当首次需要时，加速度将自动计算。

>>> WR_cont.set_vel(N, 0)
>>> WR_mc.v2pt_theory(WR_cont, N, WR)
WRrad*(u1*sin(q2) + u3 + u4)*R.x - WRrad*u2*R.y
>>> Steer.v2pt_theory(WR_mc, N, Frame)
WRrad*(u1*sin(q2) + u3 + u4)*R.x - WRrad*u2*R.y + framelength*(u1*sin(q2) + u4)*Frame.x - framelength*(-u1*sin(htangle + q4)*cos(q2) + u2*cos(htangle + q4))*Frame.y
>>> Frame_mc.v2pt_theory(WR_mc, N, Frame)
WRrad*(u1*sin(q2) + u3 + u4)*R.x - WRrad*u2*R.y + framecg3*(u1*sin(q2) + u4)*Frame.x + (-framecg1*(u1*cos(htangle + q4)*cos(q2) + u2*sin(htangle + q4)) - framecg3*(-u1*sin(htangle + q4)*cos(q2) + u2*cos(htangle + q4)))*Frame.y + framecg1*(u1*sin(q2) + u4)*Frame.z
>>> Fork_mc.v2pt_theory(Steer, N, Fork)
WRrad*(u1*sin(q2) + u3 + u4)*R.x - WRrad*u2*R.y + framelength*(u1*sin(q2) + u4)*Frame.x - framelength*(-u1*sin(htangle + q4)*cos(q2) + u2*cos(htangle + q4))*Frame.y + forkcg3*((sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2))*u1 + u2*sin(htangle + q4)*sin(q5) + u4*cos(q5))*Fork.x + (-forkcg1*((-sin(q2)*sin(q5) + cos(htangle + q4)*cos(q2)*cos(q5))*u1 + u2*sin(htangle + q4)*cos(q5) - u4*sin(q5)) - forkcg3*(-u1*sin(htangle + q4)*cos(q2) + u2*cos(htangle + q4) + u5))*Fork.y + forkcg1*((sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2))*u1 + u2*sin(htangle + q4)*sin(q5) + u4*cos(q5))*Fork.z
>>> WF_mc.v2pt_theory(Steer, N, Fork)
WRrad*(u1*sin(q2) + u3 + u4)*R.x - WRrad*u2*R.y + framelength*(u1*sin(q2) + u4)*Frame.x - framelength*(-u1*sin(htangle + q4)*cos(q2) + u2*cos(htangle + q4))*Frame.y + forkoffset*((sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2))*u1 + u2*sin(htangle + q4)*sin(q5) + u4*cos(q5))*Fork.x + (forklength*((-sin(q2)*sin(q5) + cos(htangle + q4)*cos(q2)*cos(q5))*u1 + u2*sin(htangle + q4)*cos(q5) - u4*sin(q5)) - forkoffset*(-u1*sin(htangle + q4)*cos(q2) + u2*cos(htangle + q4) + u5))*Fork.y - forklength*((sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2))*u1 + u2*sin(htangle + q4)*sin(q5) + u4*cos(q5))*Fork.z
>>> WF_cont.v2pt_theory(WF_mc, N, WF)
- WFrad*((-sin(q2)*sin(q5)*cos(htangle + q4) + cos(q2)*cos(q5))*u6 + u4*cos(q2) + u5*sin(htangle + q4)*sin(q2))/sqrt((-sin(q2)*cos(q5) - sin(q5)*cos(htangle + q4)*cos(q2))*(sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2)) + 1)*Y.x + WFrad*(u2 + u5*cos(htangle + q4) + u6*sin(htangle + q4)*sin(q5))/sqrt((-sin(q2)*cos(q5) - sin(q5)*cos(htangle + q4)*cos(q2))*(sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2)) + 1)*Y.y + WRrad*(u1*sin(q2) + u3 + u4)*R.x - WRrad*u2*R.y + framelength*(u1*sin(q2) + u4)*Frame.x - framelength*(-u1*sin(htangle + q4)*cos(q2) + u2*cos(htangle + q4))*Frame.y + (-WFrad*(sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2))*((-sin(q2)*sin(q5) + cos(htangle + q4)*cos(q2)*cos(q5))*u1 + u2*sin(htangle + q4)*cos(q5) - u4*sin(q5))/sqrt((-sin(q2)*cos(q5) - sin(q5)*cos(htangle + q4)*cos(q2))*(sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2)) + 1) + forkoffset*((sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2))*u1 + u2*sin(htangle + q4)*sin(q5) + u4*cos(q5)))*Fork.x + (forklength*((-sin(q2)*sin(q5) + cos(htangle + q4)*cos(q2)*cos(q5))*u1 + u2*sin(htangle + q4)*cos(q5) - u4*sin(q5)) - forkoffset*(-u1*sin(htangle + q4)*cos(q2) + u2*cos(htangle + q4) + u5))*Fork.y + (WFrad*(sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2))*(-u1*sin(htangle + q4)*cos(q2) + u2*cos(htangle + q4) + u5)/sqrt((-sin(q2)*cos(q5) - sin(q5)*cos(htangle + q4)*cos(q2))*(sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2)) + 1) - forklength*((sin(q2)*cos(q5) + sin(q5)*cos(htangle + q4)*cos(q2))*u1 + u2*sin(htangle + q4)*sin(q5) + u4*cos(q5)))*Fork.z 

设置每个体的惯性。使用惯性框架构建惯性偶极子。轮毂惯性仅由主惯性矩定义，并且实际上是在框架和叉齿轮参考框架中恒定的；因此，轮毂方向的定义是不需要的。框架和叉齿轮的惯性定义在固定到适当的体框架的“Temp”框架中；这是为了更轻松地输入基准论文的参考值。请注意，由于稍微不同的方向，惯性矩的积需要翻转它们的符号；这在输入数值时稍后完成。

>>> Frame_I = (inertia(TempFrame, Iframe11, Iframe22, Iframe33, 0, 0,
...                                                   Iframe31), Frame_mc)
>>> Fork_I = (inertia(TempFork, Ifork11, Ifork22, Ifork33, 0, 0, Ifork31), Fork_mc)
>>> WR_I = (inertia(Frame, Iwr11, Iwr22, Iwr11), WR_mc)
>>> WF_I = (inertia(Fork, Iwf11, Iwf22, Iwf11), WF_mc) 

声明 RigidBody 容器的开始。

>>> BodyFrame = RigidBody('BodyFrame', Frame_mc, Frame, mframe, Frame_I)
>>> BodyFork = RigidBody('BodyFork', Fork_mc, Fork, mfork, Fork_I)
>>> BodyWR = RigidBody('BodyWR', WR_mc, WR, mwr, WR_I)
>>> BodyWF = RigidBody('BodyWF', WF_mc, WF, mwf, WF_I)

>>> print('Before Forming the List of Nonholonomic Constraints.')
Before Forming the List of Nonholonomic Constraints. 

运动学微分方程；它们定义得非常简单。列表中的每个条目都等于零。

>>> kd = [q1d - u1, q2d - u2, q4d - u4, q5d - u5] 

非完整约束是前轮接触点的速度与 X、Y 和 Z 方向的点乘；由于偏航框架更接近前轮（少 1 个连接它们的 DCM），所以使用偏航框架。这些约束强制前轮接触点在惯性框架中的速度为零；X 和 Y 方向的约束强制“无滑移”条件，而 Z 方向的约束强制前轮接触点不离开地面框架，实质上复制了不允许框架俯仰以无效方式变化的完整约束。

>>> conlist_speed = [dot(WF_cont.vel(N), Y.x),
...                  dot(WF_cont.vel(N), Y.y),
...                  dot(WF_cont.vel(N), Y.z)] 

完整约束是指，从后轮接触点到前轮接触点的位置，当点乘到垂直于地面平面方向时必须为零；实际上是前后轮接触点始终接触地面平面。这实际上不是动态方程的一部分，而是线性化过程中必要的。

>>> conlist_coord = [dot(WF_cont.pos_from(WR_cont), Y.z)] 

力列表；每个体在其质心处施加适当的重力力。

>>> FL = [(Frame_mc, -mframe * g * Y.z), (Fork_mc, -mfork * g * Y.z),
...       (WF_mc, -mwf * g * Y.z), (WR_mc, -mwr * g * Y.z)]
>>> BL = [BodyFrame, BodyFork, BodyWR, BodyWF] 

惯性框架是 N 框架，坐标按独立、依赖坐标的顺序提供。在这里也输入了运动学微分方程。这里指定了独立速度，然后是依赖速度，以及非完整约束。还提供了依赖坐标，具有完整约束。这仅在线性化过程中起作用，但对于正确工作的线性化来说是必要的。

>>> KM = KanesMethod(N, q_ind=[q1, q2, q5],
...           q_dependent=[q4], configuration_constraints=conlist_coord,
...           u_ind=[u2, u3, u5],
...           u_dependent=[u1, u4, u6], velocity_constraints=conlist_speed,
...           kd_eqs=kd)
>>> print('Before Forming Generalized Active and Inertia Forces, Fr and Fr*')
Before Forming Generalized Active and Inertia Forces, Fr and Fr*
>>> (fr, frstar) = KM.kanes_equations(BL, FL)
>>> print('Base Equations of Motion Computed')
Base Equations of Motion Computed 

这是从基准论文中输入数值的开始，以验证从该模型的线性化方程到参考特征值的特征值。查看前述论文以获取更多信息。其中一些是中间值，用于将论文中的值转换为此模型中使用的坐标系。

>>> PaperRadRear  =  0.3
>>> PaperRadFront =  0.35
>>> HTA           =  evalf.N(pi/2-pi/10)
>>> TrailPaper    =  0.08
>>> rake          =  evalf.N(-(TrailPaper*sin(HTA)-(PaperRadFront*cos(HTA))))
>>> PaperWb       =  1.02
>>> PaperFrameCgX =  0.3
>>> PaperFrameCgZ =  0.9
>>> PaperForkCgX  =  0.9
>>> PaperForkCgZ  =  0.7
>>> FrameLength   =  evalf.N(PaperWb*sin(HTA) - (rake - \
...                         (PaperRadFront - PaperRadRear)*cos(HTA)))
>>> FrameCGNorm   =  evalf.N((PaperFrameCgZ - PaperRadRear - \
...                          (PaperFrameCgX/sin(HTA))*cos(HTA))*sin(HTA))
>>> FrameCGPar    =  evalf.N((PaperFrameCgX / sin(HTA) + \
...                          (PaperFrameCgZ - PaperRadRear - \
...                           PaperFrameCgX / sin(HTA) * cos(HTA)) * cos(HTA)))
>>> tempa         =  evalf.N((PaperForkCgZ - PaperRadFront))
>>> tempb         =  evalf.N((PaperWb-PaperForkCgX))
>>> tempc         =  evalf.N(sqrt(tempa**2 + tempb**2))
>>> PaperForkL    =  evalf.N((PaperWb*cos(HTA) - \
...                          (PaperRadFront - PaperRadRear)*sin(HTA)))
>>> ForkCGNorm    =  evalf.N(rake + (tempc * sin(pi/2 - \
...                          HTA - acos(tempa/tempc))))
>>> ForkCGPar     =  evalf.N(tempc * cos((pi/2 - HTA) - \
...                          acos(tempa/tempc)) - PaperForkL) 

这是数值值的最终汇编。符号‘v’是自行车的前向速度（这个概念只在静态平衡的直立情况下有意义？）。这些值将以后会被替换进字典中。再次提醒，在这里乘积惯性值的符号被翻转了，因为坐标系的不同方向。

>>> v = Symbol('v')
>>> val_dict = {
...       WFrad: PaperRadFront,
...       WRrad: PaperRadRear,
...       htangle: HTA,
...       forkoffset: rake,
...       forklength: PaperForkL,
...       framelength: FrameLength,
...       forkcg1: ForkCGPar,
...       forkcg3: ForkCGNorm,
...       framecg1: FrameCGNorm,
...       framecg3: FrameCGPar,
...       Iwr11: 0.0603,
...       Iwr22: 0.12,
...       Iwf11: 0.1405,
...       Iwf22: 0.28,
...       Ifork11: 0.05892,
...       Ifork22: 0.06,
...       Ifork33: 0.00708,
...       Ifork31: 0.00756,
...       Iframe11: 9.2,
...       Iframe22: 11,
...       Iframe33: 2.8,
...       Iframe31: -2.4,
...       mfork: 4,
...       mframe: 85,
...       mwf: 3,
...       mwr: 2,
...       g: 9.81,
...       q1: 0,
...       q2: 0,
...       q4: 0,
...       q5: 0,
...       u1: 0,
...       u2: 0,
...       u3: v/PaperRadRear,
...       u4: 0,
...       u5: 0,
...       u6: v/PaperRadFront}
>>> kdd = KM.kindiffdict()
>>> print('Before Linearization of the \"Forcing\" Term')
Before Linearization of the "Forcing" Term 

线性化强迫向量；方程式设置为 MM udot = forcing，其中 MM 是质量矩阵，udot 是表示广义速度时间导数的向量，forcing 是一个向量，其中包含外部强迫项和内部强迫项，如离心力或科里奥利力。实际上，这返回一个行数等于总坐标和速度的矩阵，但列数仅等于独立坐标和速度的数量（请注意，在此下面的内容已被注释掉，因为执行 doctests 时运行时间较长，这不利于执行）。

>>> # forcing_lin = KM.linearize()[0].subs(sub_dict) 

如前所述，线性化强迫项的大小被扩展，包括 q 和 u，因此质量矩阵也必须这样做。这可能会被更改为线性化过程的一部分，以备将来参考。

>>> MM_full = (KM._k_kqdot).row_join(zeros(4, 6)).col_join(
...           (zeros(6, 4)).row_join(KM.mass_matrix))
>>> print('Before Substitution of Numerical Values')
Before Substitution of Numerical Values 

我认为这非常容易理解。不过，花费的时间确实很长。我尝试过使用 evalf 进行替换，但由于递归深度超过了最大限制而失败；我还尝试过将其转换为 lambda 表达式，但也没有成功（同样由于速度原因注释掉了）。

>>> # MM_full = MM_full.subs(val_dict)
>>> # forcing_lin = forcing_lin.subs(val_dict)
>>> # print('Before .evalf() call')

>>> # MM_full = MM_full.evalf()
>>> # forcing_lin = forcing_lin.evalf() 

最后，我们构建了一个“A”矩阵，用于形式为 xdot = A x 的系统（其中 x 为状态向量，尽管在这种情况下，大小有点不同）。以下行仅提取了用于特征值分析所需的最小条目，这些条目对应于倾斜、转向、倾斜率和转向率的行和列（由于依赖于上述代码，这些都被注释掉了）。

>>> # Amat = MM_full.inv() * forcing_lin
>>> # A = Amat.extract([1,2,4,6],[1,2,3,5])
>>> # print(A)
>>> # print('v = 1')
>>> # print(A.subs(v, 1).eigenvals())
>>> # print('v = 2')
>>> # print(A.subs(v, 2).eigenvals())
>>> # print('v = 3')
>>> # print(A.subs(v, 3).eigenvals())
>>> # print('v = 4')
>>> # print(A.subs(v, 4).eigenvals())
>>> # print('v = 5')
>>> # print(A.subs(v, 5).eigenvals()) 

运行以上代码后，取消注释的行，比较计算得到的特征值与参考文献中的特征值。这就完成了自行车示例。

多自由度完整约束系统

原文：docs.sympy.org/latest/modules/physics/mechanics/examples/multi_degree_freedom_holonomic_system.html

在这个例子中，我们演示了sympy.physics.mechanics提供的功能，用于推导包括粒子和刚体的完整约束系统的运动方程（EOM），其中包括贡献的力和扭矩，其中一些是指定的力和扭矩。系统如下图所示：

系统将使用System进行建模。首先，我们需要创建dynamicsymbols()来描述上述图中显示的系统。在这种情况下，广义坐标(q_1)表示块离墙的横向距离，(q_2)表示复合摆杆相对垂直的角度，(q_3)表示简单摆杆相对复合摆杆的角度。广义速度(u_1)表示块的横向速度，(u_2)表示复合摆杆的横向速度，(u_3)表示 C 相对于 B 的角速度。

我们还创建一些symbols()来表示摆杆的长度和质量，以及重力和其他因素。

>>> from sympy import zeros, symbols
>>> from sympy.physics.mechanics import *
>>> q1, q2, q3, u1, u2, u3 = dynamicsymbols('q1, q2, q3, u1, u2, u3')
>>> F, T = dynamicsymbols('F, T')
>>> l, k, c, g, kT = symbols('l, k, c, g, kT')
>>> ma, mb, mc, IBzz= symbols('ma, mb, mc, IBzz') 

定义了所有符号后，我们现在可以定义各个物体，并初始化我们的System实例。

>>> wall = RigidBody('N')
>>> block = Particle('A', mass=ma)
>>> compound_pend = RigidBody('B', mass=mb)
>>> compound_pend.central_inertia = inertia(compound_pend.frame, 0, 0, IBzz)
>>> simple_pend = Particle('C', mass=mc)
>>> system = System.from_newtonian(wall)
>>> system.add_bodies(block, compound_pend, simple_pend) 

接下来，我们使用关节将物体连接起来以建立运动学。请注意，我们为两个粒子指定了中间框架，因为粒子本身没有相关联的框架。

>>> block_frame = ReferenceFrame('A')
>>> block.masscenter.set_vel(block_frame, 0)
>>> slider = PrismaticJoint('J1', wall, block, coordinates=q1, speeds=u1,
...                         child_interframe=block_frame)
>>> rev1 = PinJoint('J2', block, compound_pend, coordinates=q2, speeds=u2,
...                 joint_axis=wall.z, child_point=l*2/3*compound_pend.y,
...                 parent_interframe=block_frame)
>>> simple_pend_frame = ReferenceFrame('C')
>>> simple_pend.masscenter.set_vel(simple_pend_frame, 0)
>>> rev2 = PinJoint('J3', compound_pend, simple_pend, coordinates=q3,
...                 speeds=u3, joint_axis=compound_pend.z,
...                 parent_point=-l/3*compound_pend.y,
...                 child_point=l*simple_pend_frame.y,
...                 child_interframe=simple_pend_frame)

>>> system.add_joints(slider, rev1, rev2) 

现在我们可以对物体施加载荷（力和扭矩），重力作用于所有物体，线性弹簧和阻尼器作用于块和墙，旋转线性弹簧作用于 C 相对于 B，指定的扭矩 T 作用于复合摆杆和块，指定的力 F 作用于块。

>>> system.apply_uniform_gravity(-g * wall.y)
>>> system.add_loads(Force(block, F * wall.x))
>>> spring_damper_path = LinearPathway(wall.masscenter, block.masscenter)
>>> system.add_actuators(
...     LinearSpring(k, spring_damper_path),
...     LinearDamper(c, spring_damper_path),
...     TorqueActuator(T, wall.z, compound_pend, wall),
...     TorqueActuator(kT * q3, wall.z, compound_pend, simple_pend_frame),
... ) 

系统设置完成后，我们现在可以在后台使用KanesMethod来形成运动方程。

>>> system.form_eoms(explicit_kinematics=True)
Matrix([
[                                -c*u1(t) - k*q1(t) + 2*l*mb*u2(t)**2*sin(q2(t))/3 - l*mc*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t)))*Derivative(u3(t), t) - l*mc*(-sin(q2(t))*cos(q3(t)) - sin(q3(t))*cos(q2(t)))*(u2(t) + u3(t))**2 + l*mc*u2(t)**2*sin(q2(t)) - (2*l*mb*cos(q2(t))/3 + mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))))*Derivative(u2(t), t) - (ma + mb + mc)*Derivative(u1(t), t) + F(t)],
[-2*g*l*mb*sin(q2(t))/3 - g*l*mc*(sin(q2(t))*cos(q3(t)) + sin(q3(t))*cos(q2(t))) - g*l*mc*sin(q2(t)) + l**2*mc*(u2(t) + u3(t))**2*sin(q3(t)) - l**2*mc*u2(t)**2*sin(q3(t)) - mc*(l**2*cos(q3(t)) + l**2)*Derivative(u3(t), t) - (2*l*mb*cos(q2(t))/3 + mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))))*Derivative(u1(t), t) - (IBzz + 4*l**2*mb/9 + mc*(2*l**2*cos(q3(t)) + 2*l**2))*Derivative(u2(t), t) + T(t)],
[                                                                                                                                                                        -g*l*mc*(sin(q2(t))*cos(q3(t)) + sin(q3(t))*cos(q2(t))) - kT*q3(t) - l**2*mc*u2(t)**2*sin(q3(t)) - l**2*mc*Derivative(u3(t), t) - l*mc*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t)))*Derivative(u1(t), t) - mc*(l**2*cos(q3(t)) + l**2)*Derivative(u2(t), t)]])

>>> system.mass_matrix_full
Matrix([
[1, 0, 0,                                                                                            0,                                                                                            0,                                                     0],
[0, 1, 0,                                                                                            0,                                                                                            0,                                                     0],
[0, 0, 1,                                                                                            0,                                                                                            0,                                                     0],
[0, 0, 0,                                                                                 ma + mb + mc, 2*l*mb*cos(q2(t))/3 + mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))), l*mc*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t)))],
[0, 0, 0, 2*l*mb*cos(q2(t))/3 + mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))),                                         IBzz + 4*l**2*mb/9 + mc*(2*l**2*cos(q3(t)) + 2*l**2),                           mc*(l**2*cos(q3(t)) + l**2)],
[0, 0, 0,                                        l*mc*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))),                                                                  mc*(l**2*cos(q3(t)) + l**2),                                               l**2*mc]])

>>> system.forcing_full
Matrix([
[                                                                                                                                                                           u1(t)],
[                                                                                                                                                                           u2(t)],
[                                                                                                                                                                           u3(t)],
[                  -c*u1(t) - k*q1(t) + 2*l*mb*u2(t)**2*sin(q2(t))/3 - l*mc*(-sin(q2(t))*cos(q3(t)) - sin(q3(t))*cos(q2(t)))*(u2(t) + u3(t))**2 + l*mc*u2(t)**2*sin(q2(t)) + F(t)],
[-2*g*l*mb*sin(q2(t))/3 - g*l*mc*(sin(q2(t))*cos(q3(t)) + sin(q3(t))*cos(q2(t))) - g*l*mc*sin(q2(t)) + l**2*mc*(u2(t) + u3(t))**2*sin(q3(t)) - l**2*mc*u2(t)**2*sin(q3(t)) + T(t)],
[                                                                                -g*l*mc*(sin(q2(t))*cos(q3(t)) + sin(q3(t))*cos(q2(t))) - kT*q3(t) - l**2*mc*u2(t)**2*sin(q3(t))]]) 

一个四杆机构

原文：docs.sympy.org/latest/modules/physics/mechanics/examples/four_bar_linkage_example.html

四杆机构是力学中常见的示例，可以仅使用两个完整约束来描述。本示例将利用sympy.physics.mechanics提供的关节功能。总结来说，我们将使用刚体和关节定义开环系统。接下来，我们定义配置约束以闭环。System 将用于“整体系统的账务”，KanesMethod 作为后端处理。

链接 3 链接 2 链接 4 自由度（DoF） 活动约束 约束自由度 运动体 固定体 链接 1 q1, u1 q3, u3 q2, u2

首先，我们需要创建dynamicsymbols()来描述如上图所示的系统。在本例中，广义坐标 (q_1), (q_2) 和 (q_3) 表示链接之间的角度。同样，广义速度 (u_1), (u_2) 和 (u_3) 表示链接之间的角速度。我们还创建一些symbols()来表示链接的长度和密度。

>>> from sympy import Matrix, linear_eq_to_matrix, pi, simplify, symbols
>>> from sympy.physics.mechanics import *
>>> mechanics_printing(pretty_print=False)
>>> q1, q2, q3, u1, u2, u3 = dynamicsymbols('q1:4, u1:4')
>>> l1, l2, l3, l4, rho = symbols('l1:5, rho') 

随着所有符号的定义，我们现在可以定义物体，并初始化我们的System实例。

>>> N = ReferenceFrame('N')
>>> mass_centers = [Point(f'mc{i}') for i in range(1, 5)]
>>> inertias = [Inertia.from_inertia_scalars(P, N, 0, 0, rho*l**3/12)
...             for P, l in zip(mass_centers, (l1, l2, l3, l4))]
>>> link1 = RigidBody('Link1', frame=N, mass=rho*l1,
...                   masscenter=mass_centers[0], inertia=inertias[0])
>>> link2 = RigidBody('Link2', mass=rho*l2, masscenter=mass_centers[1],
...                   inertia=inertias[1])
>>> link3 = RigidBody('Link3', mass=rho*l3, masscenter=mass_centers[2],
...                   inertia=inertias[2])
>>> link4 = RigidBody('Link4', mass=rho*l4, masscenter=mass_centers[3],
...                   inertia=inertias[3])
>>> system = System.from_newtonian(link1) 

接下来，我们还定义了前三个关节，这些关节创建了开环摆，然后将它们添加到系统中。

>>> joint1 = PinJoint('J1', link1, link2, coordinates=q1, speeds=u1,
...                   parent_point=l1/2*link1.x,
...                   child_point=-l2/2*link2.x, joint_axis=link1.z)
>>> joint2 = PinJoint('J2', link2, link3, coordinates=q2, speeds=u2,
...                   parent_point=l2/2*link2.x,
...                   child_point=-l3/2*link3.x, joint_axis=link2.z)
>>> joint3 = PinJoint('J3', link3, link4, coordinates=q3, speeds=u3,
...                   parent_point=l3/2*link3.x,
...                   child_point=-l4/2*link4.x, joint_axis=link3.z)
>>> system.add_joints(joint1, joint2, joint3) 

现在我们可以制定会闭合运动学环路的完整约束。

>>> start_point = link1.masscenter.locatenew('start_point', -l1/2*link1.x)
>>> end_point = link4.masscenter.locatenew('end_point', l4/2*link4.x)
>>> loop = end_point.pos_from(start_point)
>>> system.add_holonomic_constraints(loop.dot(link1.x), loop.dot(link1.y)) 

在生成运动方程之前，我们需要指定哪些广义坐标和速度是独立的，哪些是依赖的。之后，我们可以运行validate_system()进行一些基本的一致性检查。

>>> system.q_ind = [q1]
>>> system.u_ind = [u1]
>>> system.q_dep = [q2, q3]
>>> system.u_dep = [u2, u3]
>>> system.validate_system() 

由于整个系统已经准备好，我们现在可以使用KanesMethod作为后端来形成运动方程。

>>> simplify(system.form_eoms())
 Matrix([[l2*rho*(-2*l2**2*sin(q3)*u1' + 3*l2*l3*u1**2*sin(q2 + q3)*sin(q2) + 3*l2*l3*sin(q2)*cos(q2 + q3)*u1' - 3*l2*l3*sin(q3)*u1' + 3*l2*l4*u1**2*sin(q2 + q3)*sin(q2) + 3*l2*l4*sin(q2)*cos(q2 + q3)*u1' + 3*l3**2*u1**2*sin(q2)*sin(q3) + 6*l3**2*u1*u2*sin(q2)*sin(q3) + 3*l3**2*u2**2*sin(q2)*sin(q3) + 2*l3**2*sin(q2)*cos(q3)*u1' + 2*l3**2*sin(q2)*cos(q3)*u2' - l3**2*sin(q3)*cos(q2)*u1' - l3**2*sin(q3)*cos(q2)*u2' + 3*l3*l4*u1**2*sin(q2)*sin(q3) + 6*l3*l4*u1*u2*sin(q2)*sin(q3) + 3*l3*l4*u2**2*sin(q2)*sin(q3) + 3*l3*l4*sin(q2)*cos(q3)*u1' + 3*l3*l4*sin(q2)*cos(q3)*u2' + l4**2*sin(q2)*u1' + l4**2*sin(q2)*u2' + l4**2*sin(q2)*u3')/(6*sin(q3))]]) 

揭示非贡献力

要揭示闭合关节处的非贡献力，我们必须在端点引入 x 和 y 方向的辅助速度。

>>> uaux1, uaux2 = dynamicsymbols('uaux1:3')
>>> end_point_aux = end_point.locatenew('end_point_aux', 0)
>>> end_point_aux.set_vel(N, end_point.vel(N) + uaux1*N.x + uaux2*N.y) 

要确保速度包含在速度约束中，我们必须手动覆盖速度约束，因为这些默认情况下指定为完整约束的时间导数。

>>> system.velocity_constraints = [
...    end_point_aux.vel(N).dot(N.x), end_point_aux.vel(N).dot(N.y)] 

在添加非贡献力时，我们需要它们仅依赖于辅助速度，而不是被约束消除的速度。这可以通过对非辅助端点施加等大而相反的力来实现。

>>> faux1, faux2 = dynamicsymbols('faux1:3')
>>> noncontributing_forces = [
...   Force(end_point_aux, faux1*N.x + faux2*N.y),
...   Force(end_point, -(faux1*N.x + faux2*N.y)),
... ] 

或者，我们可以指定一个新点，该点已经减去了被约束消除的速度。

>>> end_point_forces = end_point.locatenew('end_point_forces', 0)
>>> end_point_forces.set_vel(N, uaux1*N.x + uaux2*N.y)
>>> noncontributing_forces = [Force(end_point_forces, faux1*N.x + faux2*N.y)] 

接下来，我们可以将辅助速度和非贡献力添加到系统中。

>>> system.add_loads(*noncontributing_forces)
>>> system.u_aux = [uaux1, uaux2] 

要在验证系统并形成运动方程之前包含重力，我们可以使用apply_uniform_gravity()。

>>> g = symbols('g')
>>> system.apply_uniform_gravity(-g*N.y)
>>> system.validate_system()
>>> eoms = system.form_eoms() 

当形成运动方程后，我们可以解辅助方程以计算一个简单配置的非贡献力的值。

>>> auxiliary_eqs = system.eom_method.auxiliary_eqs
>>> forces_eqs = Matrix.LUsolve(
...   *linear_eq_to_matrix(auxiliary_eqs, [faux1, faux2]))
>>> subs = {
...   l1: 2, l2: 1, l3: 2, l4: 1,
...   rho: 5, g: 9.81,
...   q1: pi/2, q2: pi/2, q3: pi/2,
...   u1: 0, u2: 0, u3: 0, u1.diff(): 0, u2.diff(): 0, u3.diff(): 0,
... }
>>> forces_eqs.xreplace(subs)
Matrix([
[    0],
[-98.1]]) 

物理/力学中的潜在问题/高级主题/未来特性

原文链接：docs.sympy.org/latest/modules/physics/mechanics/advanced.html

本文将描述此模块提供的一些更高级的功能，但这些功能不是“官方”接口的一部分。在这里，还将涵盖一些将来将实施的特性，以及关于正确功能的未解答问题。同时，还将讨论常见问题，并提供一些解决方案。

常见问题

在这里，数值积分代码、选择dynamicsymbols作为坐标和速度表示、打印、差分和替换的问题将会出现。

数值积分代码

参见未来特性：代码输出

差分

在 SymPy 中，差分非常大的表达式可能需要一些时间；大表达式可能需要几分钟才能评估出导数。这通常会出现在线性化中。

坐标和速度的选择

Kane 对象被设置为假设广义速度与广义坐标的时间导数不是相同的符号。这并不意味着它们不能是相同的，只是它们必须有不同的符号。如果你这样做了：

>> KM.coords([q1, q2, q3])
>> KM.speeds([q1d, q2d, q3d]) 

你的代码将无法工作。目前，运动微分方程需要提供。希望用户在这一点上能发现他们不应尝试上述代码中显示的行为。

这种行为对于形成运动方程的其他方法可能不适用。

打印

默认的打印选项是对Vector和Dyad测量数使用排序，并且从mprint、mpprint和mlatex函数中输出未排序的内容。如果要打印大量内容，请使用这些函数之一，因为排序可能会将打印时间从几秒增加到几分钟。

替换

在力学中有两个常见的替代问题：

• 当为dynamicsymbols进行表达式替换时，sympy 的正常subs将会为动态符号的导数进行替换：

  >>> from sympy.physics.mechanics import dynamicsymbols
  >>> x = dynamicsymbols('x')
  >>> expr = x.diff() + x
  >>> sub_dict = {x: 1}
  >>> expr.subs(sub_dict)
  Derivative(1, t) + 1 
  

  在这种情况下，x也被替换为Derivative内部的 1，这是不希望的。

• 对大表达式进行替换可能会很慢。

如果你的替换很简单（直接用其他表达式替换表达式，例如在操作点评估时），建议使用提供的msubs函数，因为它的速度显著更快，并适当处理了导数问题：

>>> from sympy.physics.mechanics import msubs
>>> msubs(expr, sub_dict)
Derivative(x(t), t) + 1 

线性化

目前，线性化方法不支持存在非坐标、非速度动态符号的情况外的“动态方程”情况。它也不支持这些类型动态符号的时间导数出现的情况。这意味着如果您有具有非坐标、非速度动态符号的运动微分方程，它将无法工作。这也意味着如果您已将系统参数（如长度、距离或质量）定义为动态符号，则其时间导数可能会出现在动态方程中，并且这将阻止线性化。

点的加速度

至少需要定义点的速度，因为在同一框架内，可以通过速度的时间导数计算加速度。如果使用了一点或两点定理来计算速度，则速度表达式的时间导数很可能比使用加速度一级和两级定理更复杂。在这一点上使用加速度级别的方法可以导致更短的表达式，这将在后续形成 Kane 方程时导致更短的表达式。

高级接口

高级功能

记住，Kane对象支持具有时间变化质量和惯性的物体，尽管这种功能与线性化方法并不完全兼容。

之前讨论过运算符作为在Vector和Dyad对象上进行数学运算的潜在方式。该模块中的大多数代码实际上是用它们编写的，因为这可以（主观地）导致更干净、更短、更可读的代码。如果在您的代码中使用此接口，请记住小心使用括号；Python 中的默认运算顺序导致在一些向量积之前发生加法，因此应大量使用括号。

未来特性

这将涵盖计划添加到此子模块的功能。

代码输出

实现数值积分代码输出函数是下一个要实现的最重要功能。这里有许多考虑因素。

C 代码输出（使用 GSL 库）、Fortran 90（使用 LSODA）、MATLAB 和 SciPy 是目标。需要考虑的事项包括：在 MATLAB 和 SciPy 上对大表达式使用cse，这是解释性的。目前不清楚编译语言是否会从常见子表达式消除中受益，特别是考虑到它是编译器优化的常见部分，并且在调用cse时可能会有显著的时间惩罚。

在构建这些表达式的字符串时需要小心，以及处理输入参数和其他动态符号。在进行积分时如何处理输出量也需要决定，考虑到可能有多个选项。

物理/力学参考文献

原文链接：docs.sympy.org/latest/modules/physics/mechanics/reference.html

[Blajer1994]

Blajer, Wojciech, Werner Schiehlen, and Walter Schirm. “A projective criterion to the coordinate partitioning method for multibody dynamics.” Archive of Applied Mechanics 64: 86-98. Print.

[Kane1983]

Kane, Thomas R., Peter W. Likins, and David A. Levinson. Spacecraft Dynamics. New York: McGraw-Hill Book, 1983. Print.

[Kane1985]

Kane, Thomas R., and David A. Levinson. Dynamics, Theory and Applications. New York: McGraw-Hill, 1985. Print.

[Meijaard2007]

Meijaard, J.P., Jim M. Papadopoulos, Andy Ruina, and A.L. Schwab. “Linearized Dynamics Equations for the Balance and Steer of a Bicycle: a Benchmark and Review.” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 463.2084 (2007): 1955-982. Print.

Autolev 解析器

原文链接：docs.sympy.org/latest/modules/physics/mechanics/autolev_parser.html

介绍

Autolev（现在被 MotionGenesis 取代）是用于符号多体动力学的特定领域语言。SymPy 力学模块现在具有足够的能力和功能，可以作为一个完整的符号动力学模块。该解析器通过利用 SymPy 的数学库和力学模块，将 Autolev（版本 4.1）代码解析为 SymPy 代码。

该解析器是使用 ANTLR 框架构建的，其主要目的是帮助 Autolev 的旧用户熟悉 SymPy 中的多体动力学。

下面的部分将讨论解析器的细节，如使用方法、常见问题、问题和未来的改进。如果你想详细比较 Autolev 和 SymPy Mechanics，可以参考 SymPy Mechanics for Autolev Users guide。

使用方法

我们首先从一个 Autolev 代码文件开始。

让我们以这个例子为例（注释 % 用于展示 Autolev 的响应）：

% double_pendulum.al
%-------------------
MOTIONVARIABLES' Q{2}', U{2}'
CONSTANTS L,M,G
NEWTONIAN N
FRAMES A,B
SIMPROT(N, A, 3, Q1)
% -> N_A = [COS(Q1), -SIN(Q1), 0; SIN(Q1), COS(Q1), 0; 0, 0, 1]
SIMPROT(N, B, 3, Q2)
% -> N_B = [COS(Q2), -SIN(Q2), 0; SIN(Q2), COS(Q2), 0; 0, 0, 1]
W_A_N>=U1*N3>
% -> W_A_N> = U1*N3>
W_B_N>=U2*N3>
% -> W_B_N> = U2*N3>
POINT O
PARTICLES P,R
P_O_P> = L*A1>
% -> P_O_P> = L*A1>
P_P_R> = L*B1>
% -> P_P_R> = L*B1>
V_O_N> = 0>
% -> V_O_N> = 0>
V2PTS(N, A, O, P)
% -> V_P_N> = L*U1*A2>
V2PTS(N, B, P, R)
% -> V_R_N> = L*U1*A2> + L*U2*B2>
MASS P=M, R=M
Q1' = U1
Q2' = U2
GRAVITY(G*N1>)
% -> FORCE_P> = G*M*N1>
% -> FORCE_R> = G*M*N1>
ZERO = FR() + FRSTAR()
% -> ZERO[1] = -L*M*(2*G*SIN(Q1)+L*(U2²*SIN(Q1-Q2)+2*U1'+COS(Q1-Q2)*U2'))
% -> ZERO[2] = -L*M*(G*SIN(Q2)-L*(U1²*SIN(Q1-Q2)-U2'-COS(Q1-Q2)*U1'))
KANE()
INPUT M=1,G=9.81,L=1
INPUT Q1=.1,Q2=.2,U1=0,U2=0
INPUT TFINAL=10, INTEGSTP=.01
CODE DYNAMICS() some_filename.c 

解析器可以如下使用：

>>> from sympy.parsing.autolev import parse_autolev
>>> sympy_code = parse_autolev(open('double_pendulum.al'), include_numeric=True)

# The include_pydy flag is False by default. Setting it to True will
# enable PyDy simulation code to be outputted if applicable.

>>> print(sympy_code)
import sympy.physics.mechanics as me
import sympy as sm
import math as m
import numpy as np

q1, q2, u1, u2 = me.dynamicsymbols('q1 q2 u1 u2')
q1d, q2d, u1d, u2d = me.dynamicsymbols('q1 q2 u1 u2', 1)
l, m, g=sm.symbols('l m g', real=True)
frame_n=me.ReferenceFrame('n')
frame_a=me.ReferenceFrame('a')
frame_b=me.ReferenceFrame('b')
frame_a.orient(frame_n, 'Axis', [q1, frame_n.z])
# print(frame_n.dcm(frame_a))
frame_b.orient(frame_n, 'Axis', [q2, frame_n.z])
# print(frame_n.dcm(frame_b))
frame_a.set_ang_vel(frame_n, u1*frame_n.z)
# print(frame_a.ang_vel_in(frame_n))
frame_b.set_ang_vel(frame_n, u2*frame_n.z)
# print(frame_b.ang_vel_in(frame_n))
point_o=me.Point('o')
particle_p=me.Particle('p', me.Point('p_pt'), sm.Symbol('m'))
particle_r=me.Particle('r', me.Point('r_pt'), sm.Symbol('m'))
particle_p.point.set_pos(point_o, l*frame_a.x)
# print(particle_p.point.pos_from(point_o))
particle_r.point.set_pos(particle_p.point, l*frame_b.x)
# print(particle_p.point.pos_from(particle_r.point))
point_o.set_vel(frame_n, 0)
# print(point_o.vel(frame_n))
particle_p.point.v2pt_theory(point_o,frame_n,frame_a)
# print(particle_p.point.vel(frame_n))
particle_r.point.v2pt_theory(particle_p.point,frame_n,frame_b)
# print(particle_r.point.vel(frame_n))
particle_p.mass = m
particle_r.mass = m
force_p = particle_p.mass*(g*frame_n.x)
# print(force_p)
force_r = particle_r.mass*(g*frame_n.x)
# print(force_r)
kd_eqs = [q1d - u1, q2d - u2]
forceList = [(particle_p.point,particle_p.mass*(g*frame_n.x)), (particle_r.point,particle_r.mass*(g*frame_n.x))]
kane = me.KanesMethod(frame_n, q_ind=[q1,q2], u_ind=[u1, u2], kd_eqs = kd_eqs)
fr, frstar = kane.kanes_equations([particle_p, particle_r], forceList)
zero = fr+frstar
# print(zero)
#---------PyDy code for integration----------
from pydy.system import System
sys = System(kane, constants = {l:1, m:1, g:9.81},
specifieds={},
initial_conditions={q1:.1, q2:.2, u1:0, u2:0},
times = np.linspace(0.0, 10, 10/.01))

y=sys.integrate() 

注释代码不属于输出代码的一部分。打印语句展示了如何获取类似 Autolev 文件中的响应。请注意，我们需要在许多情况下使用 SymPy 函数，如 .ang_vel_in()、.dcm() 等，而不是直接打印变量如 zero。如果你完全是 SymPy 力学的新手，SymPy Mechanics for Autolev Users guide 应该会有所帮助。你可能还需要使用基本的 SymPy 简化和操作，如 trigsimp()、expand()、evalf() 等，以获取类似 Autolev 的输出。更多相关信息请参考 SymPy Tutorial。 ## 注意事项

• 不要使用与 Python 保留字冲突的变量名。这是一个违反规则的例子：

  %Autolev Code
  %------------
  LAMBDA = EIG(M) 
  

  #SymPy Code
  #----------
  lambda = sm.Matrix([i.evalf() for i in (m).eigenvals().keys()]) 
  

---

• 确保矢量和标量的名称不同。Autolev 将这些视为不同的内容，但在 Python 中会被覆盖。目前解析器允许体和标量/矢量的名称相同，但不允许标量和矢量之间相同。这在将来可能需要改变。

  %Autolev Code
  %------------
  VARIABLES X,Y
  FRAMES A
  A> = X*A1> + Y*A2>
  A = X+Y 
  

  #SymPy Code
  #----------
  x, y = me.dynamicsymbols('x y')
  frame_a = me.ReferenceFrame('a')
  a = x*frame_a.x + y*frame_a.y
  a = x + y
  # Note how frame_a is named differently so it doesn't cause a problem.
  # On the other hand, 'a' gets rewritten from a scalar to a vector.
  # This should be changed in the future. 
  

---

• 当处理函数返回的矩阵时，必须检查值的顺序，因为它们可能与 Autolev 中的不同。特别是对于特征值和特征向量而言。

  %Autolev Code
  %------------
  EIG(M, E1, E2)
  % -> [5; 14; 13]
  E2ROW = ROWS(E2, 1)
  EIGVEC> = VECTOR(A, E2ROW) 
  

  #SymPy Code
  #----------
  e1 = sm.Matrix([i.evalf() for i in m.eigenvals().keys()])
  # sm.Matrix([5;13;14]) different order
  e2 = sm.Matrix([i[2][0].evalf() for i in m.eigenvects()]).reshape(m.shape[0], m.shape[1])
  e2row = e2.row(0)
  # This result depends on the order of the vectors in the eigenvecs.
  eigenvec = e2row[0]*a.x + e2row[1]*a.y + e2row[2]*a.y 
  

---

• 当使用 EVALUATE 时，如 90*UNITS(deg,rad) 用于角替代，因为在 SymPy 中弧度是默认的。你也可以直接在 SymPy 代码中添加 np.deg2rad()。

  在解析 CODE 命令生成的输出代码时，不需要这样做，因为解析器在 INPUT 声明中给出 deg 单位时会处理这个问题。

  另一方面，DEGREES 设置仅在某些情况下有效，例如在需要角度的 SIMPROT 中。

  %Autolev Code
  %------------
  A> = Q1*A1> + Q2*A2>
  B> = EVALUATE(A>, Q1:30*UNITS(DEG,RAD)) 
  

  #SymPy Code
  #----------
  a = q1*a.frame_a.x + q2*frame_a.y
  b = a.subs({q1:30*0.0174533})
  # b = a.subs({q1:np.deg2rad(30)} 
  

---

• 大多数 Autolev 设置尚未解析，对解析器没有影响。唯一部分起作用的是 COMPLEX 和 DEGREES。建议查找 SymPy 和 Python 中的替代方案。

---

• REPRESENT 命令不受支持。请改用 MATRIX、VECTOR 或 DYADIC 命令。Autolev 4.1 也建议使用这些命令而不是 REPRESENT，虽然仍允许使用，但解析器不解析它。

---

• 不要使用类似 WO{3}RD{2,4} 的变量声明类型。解析器只能处理一个变量名后跟一对花括号和任意数量的 '。如果要实现类似 WO{3}RD{2,4} 的效果，必须手动声明所有情况。

---

• 解析器可以处理大多数命令的普通版本，但在大多数情况下可能无法正确解析带有矩阵参数的函数。例如：

  M=COEF([E1;E2],[U1,U2,U3])

  这将计算 E1 和 E2 中 U1、U2 和 U3 的系数。最好使用这些命令的常规版本手动构造一个矩阵。

  %Autolev Code
  %------------
  % COEF([E1;E2],[U1,U2,U3])
  M = [COEF(E1,U1),COEF(E1,U2),COEF(E1,U3) &
      ;COEF(E2,U1),COEF(E2,U2),COEF(E2,U3)] 
  

---

• 必须使用 MOTIONVARIABLE 声明广义坐标和速度，并且必须在常规的 VARIABLE 声明中声明所有其他变量。解析器需要这样做以区分它们，并将正确的参数传递给 Kane 方法对象。

  还建议始终声明与坐标对应的速度，并传递运动微分方程。解析器能够处理某些情况，即使没有这些也引入自己的一些虚拟变量，但 SymPy 自身确实需要它们。

  还要注意，像 VARIABLES U{3}' 这样的旧 Autolev 声明也不受支持。

  %Autolev Code
  %------------
  MOTIONVARIABLES' Q{2}', U{2}'
  % ----- OTHER LINES ----
  Q1' = U1
  Q2' = U2
  ----- OTHER LINES ----
  ZERO = FR() + FRSTAR() 
  

  #SymPy Code
  #----------
  q1, q2, u1, u2 = me.dynamicsymbols('q1 q2 u1 u2')
  q1d, q2d, u1d, u2d = me.dynamicsymbols('q1 q2 u1 u2', 1)
  
  # ------- other lines -------
  
  kd_eqs = [q1d - u1, q2d - u2]
  kane = me.KanesMethod(frame_n, q_ind=[q1,q2], u_ind=[u1, u2], kd_eqs = kd_eqs)
  fr, frstar = kane.kanes_equations([particle_p, particle_r], forceList)
  zero = fr+frstar 
  

---

• 需要在所有出现的 Kane 方程中将 me.dynamicsymbols._t 更改为 me.dynamicsymbols('t')。例如，查看此 弹簧阻尼器示例 的第 10 行。此方程用于形成 Kane 方程，因此在这种情况下需要将 me.dynamicsymbols._t 更改为 me.dynamicsymbols('t')。

  这样做的主要原因是因为 PyDy 要求明确列出时间相关的指定变量，而 Autolev 则通过方程自动处理方程中的时间变量。

  问题在于 PyDy 的 System 类不接受 dynamicsymbols._t 作为指定。参见问题 #396。这种改变实际上并不理想，因此未来应该找到更好的解决方案。

---

• 解析器通过解析 Autolev 代码中的变量声明来创建 SymPy 的 symbols 和 dynamicsymbols。

  对于直接初始化的中间表达式，解析器不会创建 SymPy 符号，只是将它们分配给表达式。

  另一方面，当声明变量被赋予表达式时，解析器将表达式存储在字典中，以避免将其重新分配给完全不同的实体。这个约束是由 Python 的固有特性及其与 Autolev 等语言的不同之处所决定的。

  另外，Autolev 似乎能够假定在某些情况下使用变量或 rhs 表达式，即使在表达式中没有明确的RHS()调用。然而，为了解析器能够正常工作，最好在变量的 rhs 表达式应使用RHS()。

  %Autolev Code
  %------------
  VARIABLES X, Y
  E = X + Y
  X = 2*Y
  
  RHS_X = RHS(X)
  
  I1 = X
  I2 = Y
  I3 = X + Y
  
  INERTIA B,I1,I2,I3
  % -> I_B_BO>> = I1*B1>*B1> + I2*B2>*B2> + I3*B3>*B3> 
  

  #SymPy Code
  #----------
  x,y = me.dynamicsymbols('x y')
  e = x + y  # No symbol is made out of 'e'
  
  # an entry like {x:2*y} is stored in an rhs dictionary
  
  rhs_x = 2*y
  
  i1 = x  # again these are not made into SymPy symbols
  i2 = y
  i3 = x + y
  
  body_b.inertia = (me.inertia(body_b_f, i1, i2, i3), b_cm)
  # This prints as:
  # x*b_f.x*b_f.x + y*b_f.y*b_f.y + (x+y)*b_f.z*b_f.z
  # while Autolev's output has I1,I2 and I3 in it.
  # Autolev however seems to know when to use the RHS of I1,I2 and I3
  # based on the context. 
  

---

• 解析SOLVE命令的方法如下：

  %Autolev Code
  %------------
  SOLVE(ZERO,X,Y)
  A = RHS(X)*2 + RHS(Y) 
  

  #SymPy Code
  #----------
  print(sm.solve(zero,x,y))
  # Behind the scenes the rhs of x
  # is set to sm.solve(zero,x,y)[x].
  a = sm.solve(zero,x,y)[x]*2 + sm.solve(zero,x,y)[y] 
  

  像[x]和[y]这样的索引并不总是有效，因此您可能需要查看 solve 返回的底层字典并正确索引它。

---

• 在解析器的上下文中，惯性声明和惯性函数的工作方式略有不同。这一点起初可能很难理解，但由于 SymPy 和 Autolev 之间的差异，必须这样做以弥合差距。以下是它们的一些要点：

  1. 惯性声明(INERTIA B,I1,I2,I3)设置刚体的惯性。

  2. 形式为I_C_D>> = expr的惯性设置器只有在 C 是物体时才设置惯性。如果 C 是粒子，则I_C_D>> = expr仅简单地解析为i_c_d = expr，并且i_c_d表现得像一个普通变量。

  3. 当涉及惯性获取器（在表达式中使用的I_C_D>>或INERTIA命令中使用）时，这些必须与EXPRESS命令一起使用以指定框架，因为 SymPy 需要此信息来计算惯性二阶张量。

  %Autolev Code
  %------------
  INERTIA B,I1,I2,I3
  I_B_BO>> = X*A1>*A1> + Y*A2>*A2>  % Parser will set the inertia of B
  I_P_Q>> = X*A1>*A1> + Y²*A2>*A2> % Parser just parses it as i_p_q = expr
  
  E1 = 2*EXPRESS(I_B_O>>,A)
  E2 =  I_P_Q>>
  E3 = EXPRESS(I_P_O>>,A)
  E4 = EXPRESS(INERTIA(O),A)
  
  % In E1 we are using the EXPRESS command with I_B_O>> which makes
  % the parser and SymPy compute the inertia of Body B about point O.
  
  % In E2 we are just using the dyadic object I_P_Q>> (as I_P_Q>> = expr
  % doesn't act as a setter) defined above and not asking the parser
  % or SymPy to compute anything.
  
  % E3 asks the parser to compute the inertia of P about point O.
  
  % E4 asks the parser to compute the inertias of all bodies wrt about O. 
  

---

• 在物体的惯性声明中，如果惯性是围绕除了质心以外的点设置的，则需要确保该点的位置向量设置器和质心在惯性声明之前出现，否则 SymPy 会抛出错误。

  %Autolev Code
  %------------
  P_SO_O> = X*A1>
  INERTIA S_(O) I1,I2,I3 
  

---

• 注意，并非所有的 Autolev 命令都已实现。解析器现在覆盖了它们的基本形式中的重要部分。如果您对是否包含某个命令感到怀疑，请查看此文件中的源代码。搜索“”以验证此事。查看特定命令的代码也将有助于了解预期的工作形式。 ## 限制和问题

• 很多问题已经在“陷阱”部分讨论过了。其中一些是：

  • 在 Python 中，向量名称与标量名称重合时会被覆盖。

  • 某些方便的变量声明没有被解析。

  • 某些返回矩阵的方便形式的函数没有被解析。

  • 设置没有被解析。

  • Python 中的符号和 rhs 表达式工作方式非常不同，这可能导致不良结果。

  • 解析SOLVE命令的代码的字典索引在许多情况下并不正确。

  • 需要将dynamicsymbols._t更改为dynamicsymbols('t')以使 PyDy 仿真代码正常工作。

这里还有其他一些：

• 特征向量似乎没有按预期工作。在许多情况下，Autolev 和 SymPy 中的值并不相同。

• 解析器不能解析块矩阵。实际上，通过对 SymPy 进行更改以允许矩阵接受其他矩阵作为参数将更容易实现这一点。

• SymPy 中TAYLOR命令的等效.series()与dynamicsymbols()不兼容。

• 只有当前的DEPENDENT约束条件被解析。需要解析AUXILIARY约束条件。这应该很快完成，因为这并不是很困难。

• 目前没有正确解析任何能量和动量函数。将这些功能也搞定会很好。可能需要对 SymPy 进行一些更改。例如，SymPy 没有等效于NICHECK()的函数。

• 数值积分部分目前只在没有参数的KANE命令的情况下正常工作。像KANE(F1,F2)这样的事情目前不起作用。

• 此外，PyDy 数值仿真代码仅适用于像ZERO = FR() + FRSTAR()这样的情况下。当矩阵插入其他方程时效果不佳。在实现这一点时遇到的一个障碍是 PyDy 的 System 类自动接受forcing_full和mass_matrix_full并解决它们，而不给用户指定方程的灵活性。希望能向 System 类添加此功能。## 未来改进

1. 完成在线动力学

解析器是通过参考和解析来自Autolev 教程和书籍Dynamics Online: Theory and Implementation Using Autolev中的代码构建的。基本上，这个过程涉及通过每个代码，验证解析器的结果，并改进规则（如果需要）以确保代码解析良好。

这些解析代码可在 GitLab 这里找到。仓库是私有的，因此需要请求访问权限。截至目前，大部分Dynamics Online第四章的代码已经解析完成。

完成书中所有剩余的代码（即2-10、2-11、第四章剩余部分、第五章和第六章（较不重要））将使解析器更加完整。

2. 修复问题

第二件事是按照优先级和易用性顺序解决上述“Gotchas”和“Limitations and Issues”部分中描述的问题。其中许多问题需要对解析器代码进行更改，而有些问题最好通过向 SymPy 添加一些功能来解决。

3. 切换至 AST

当前解析器使用一种具体语法树（CST）来构建，使用ANTLR框架。理想情况下，从 CST 切换到抽象语法树（AST）会更好。这样，解析器代码将不再依赖于 ANTLR 语法，使其更加灵活。同时，更改语法和解析器规则也会更加容易。

适用于 Autolev 用户的 SymPy Mechanics

原文链接：docs.sympy.org/latest/modules/physics/mechanics/sympy_mechanics_for_autolev_users.html

介绍

Autolev（现在被 MotionGenesis 取代）是一种特定领域的编程语言，用于符号多体动力学。SymPy mechanics 模块现在具有足够的功能和功能，可以成为一个完全功能的符号动力学模块。PyDy 包将 SymPy 输出扩展到数值领域，用于仿真、分析和可视化。Autolev 和 SymPy Mechanics 有很多共同点，但它们之间也有许多不同之处。本页将扩展这些差异，旨在成为 Autolev 用户转向 SymPy Mechanics 的参考。

在浏览本页之前，了解 SymPy 和 SymPy Mechanics 的基本理解将非常有帮助。如果您完全是 Python 的新手，可以查看官方Python 教程。尤其是查看 SymPy 文档中的教程，以便对 SymPy 有所了解。此外，对于想要了解 Python 中多体动力学的介绍，这个讲座非常有帮助。

您可能还会发现 Autolev 解析器，它是 SymPy 的一部分，非常有帮助。

一些关键差异

Autolev	SymPy Mechanics
Autolev 是一种特定领域的编程语言，旨在执行多体动力学。由于它是一种独立的语言，它有非常严格的语言规范。它根据输入代码预定义、假设和计算许多内容。因此，它的代码更加清晰简洁。	SymPy 是用通用目的语言 Python 编写的库。虽然 Autolev 的代码更紧凑，但 SymPy（作为 Python 的补充）更加灵活。用户可以更多地控制他们可以做的事情。例如，可以在他们的代码中创建一个类，用于定义具有共同属性的刚体类型。还有广泛的科学 Python 库可供使用，这也是一个重要的优势。
Autolev 从一小部分符号数学生成 Matlab、C 或 Fortran 代码。	SymPy 从使用 SymPy 创建的大量符号数学生成数值 Python、C 或 Octave/Matlab 代码。它还建立在流行的科学 Python 堆栈上，如 NumPy、SciPy、IPython、matplotlib、Cython 和 Theano。
Autolev 使用基于 1 的索引。序列的初始元素使用 a[1]找到。	Python 使用基于 0 的索引。序列的初始元素使用 a[0]找到。
Autolev 不区分大小写。	SymPy 代码作为 Python 代码区分大小写。
在 Autolev 中，可以通过创建 .R 和 .A 文件来定义自己的命令，这些文件可以在程序中使用。	SymPy 代码是 Python 代码，因此可以在代码中定义函数。这要方便得多。
Autolev 是专有的。	SymPy 是开源的。

粗略的 Autolev-SymPy 等价物

下表提供了一些常见 Autolev 表达式的粗略等价物。这些并非完全等价，而是应视为指导以便朝正确方向前进。更多细节请参阅内置文档 SymPy vectors、SymPy mechanics 和 PyDy 。

在下表中，假设您已在 Python 中执行了以下命令：

import sympy.physics.mechanics as me
import sympy as sm 

数学等价物

Autolev	SymPy	备注
Constants A, B	a, b = sm.symbols('a b', real=True)	注意，符号的名称可以与其分配的变量的名称不同。我们可以定义 a, b = symbols('b a')，但遵循惯例是个好习惯。
Constants C+	c = sm.symbols('c', real=True, nonnegative=True)	更多信息请参阅 SymPy assumptions。
Constants D-	d = sm.symbols('d', real=True, nonpositive=True)
Constants K{4}	k1, k2, k3, k4 = sm.symbols('k1 k2 k3 k4', real=True)
Constants a{2:4}	a2, a3, a4 = sm.symbols('a2 a3 a4', real=True)
Constants b{1:2, 1:2}	b11, b12, b21, b22 = sm.symbols('b11 b12 b21 b22', real=True)
Specified Phi	phi = me.dynamicsymbols('phi')
Variables q, s	q, s = me.dynamicsymbols('q', 's')
Variables x''	x = me.dynamicsymbols('x')``xd = me.dynamicsymbols('x', 1)``xd2 = me.dynamicsymbols('x', 2)
Variables y{2}'	y1 = me.dynamicsymbols('y1')``y2 = me.dynamicsymbols('y2')``y1d = me.dynamicsymbols('y1', 1)``y2d = me.dynamicsymbols('y2', 1)
MotionVariables u{2}	u1 = me.dynamicsymbols('u1')``u2 = me.dynamicsymbols('u2')	SymPy 在声明过程中不区分变量、运动变量和指定项。相反，它将这些不同类型的变量作为参数传递给像 KanesMethod 这样的对象。
Imaginary j	j = sm.I	I 是一个 sympy 对象，代表虚数单位。可以使用它定义复数，例如 z = x + I*y，其中 x、y 和 z 是符号。
Tina = 2*pi``s = u*t + a*t²/2	tina = 2*sm.pi``tina = tina.evalf()``t = me.dynamicsymbols._t``s = u*t + a*t**2/2	使用 .evalf() 将得到数值结果。
abs(x)³ + sin(x)² + acos(x)	sm.abs(x)**3 + sm.sin(x)**2 + sm.acos(x)
E = (x+2*y)² + 3*(7+x)*(x+y) Expand(E) Factor(E, x) Coef(y, x) Replace(y, sin(x)=3) Exclude(E,x) Include(E,x) Arrange(E,2,y)	E = (x+2*y)**2 + 3*(7+x)*(x+y) sm.expand(E) sm.horner(E, wrt=x) y.coeff(x) y.subs({sm.sin(x): 3}) e.collect(x).coeff(x, 0) e.collect(x).coeff(x, 1) e.collect(y)要了解更多信息，请参阅 simplification.这些 SymPy 函数不会原地操作。它们只返回表达式。如果要覆盖原始表达式，可以像这样操作：y = y.subs({sm.sin(x): 3})
Dy = D(E, y) Dt = Dt(E) Dt2 = Dt(V, A)，其中 V 是向量，A 是框架。Dy2 = D(V, y, A)	E.diff(y) E.diff(me.dynamicsymbols._t)如果表达式由动力符号组成，则有效。dt2 = v.dt(A) dy2 = v.diff(y, A)要了解更多信息，请参阅 calculus.
E = COS(X*Y) TY = Taylor(E, 0:2, x=0, y=0)	e = sm.cos(x*y) b = e.series(x, 0, 2).removeO().series(y, 0, 2).removeO()要了解更多信息，请参阅 series.
F = Evaluate(E, x=a, y=2)	E.subs([(x, a), (y, 2)])要从数值表达式中获得浮点数，请使用.evalf()。E.evalf((a + sm.pi).subs({a: 3}))
P = Polynomial([a, b, c], x)	p = sm.Poly(sm.Matrix([a, b, c]).reshape(1, 3), x)要了解更多信息，请参阅 polys.
Roots(Polynomial(a*x² + b*x + c, x, 2) Roots([1;2;3])	sm.solve(sm.Poly(a*x**2 + b*x + c)) sm.solve(sm.Poly(sm.Matrix([1,2,3]).reshape(3, 1), x), x)要了解更多信息，请参阅 Solvers。有关与多项式和根相关的数值计算，请参阅mpmath/calculus.
解(A, x1, x2)其中 A 是表示线性方程的增广矩阵，x1, x2 是要解的变量。	sm.linsolve(A, (x1, x2))其中 A 是增广矩阵。要了解更多信息，请参阅：:ref: solvers/solveset. <solveset>要查看非线性求解器，请参阅 solvers.
RowMatrix = [1, 2, 3, 4] ColMatrix = [1; 2; 3; 4] MO = [a, b; c, 0] MO[2, 2] := d A + B*C Cols(A) Cols(A, 1) Rows(A) Rows(A, 1) Det(A) Element(A, 2, 3) Inv(A) Trace(A) Transpose(A) Diagmat(4, 1) Eig(A) Eig(A, EigVal, EigVec)	row_matrix = sm.Matrix([[1],[2], [3],[4]]) col_matrix = sm.Matrix([1, 2, 3, 4]) MO = sm.Matrix([[a, b], [c, 0]]) MO[1, 1] = d A + B*C A.cols A.col(0) A.rows A.row(0) M.det() M[2, 3] M**-1 sm.trace(A) A.T sm.diag(1,1,1,1) A.eigenvals() eigval = A.eigenvals() eigvec = A.eigenvects()要了解更多信息，请参阅 matrices.

物理等效

Autolev	SymPy	Notes
Bodies A声明了 A、其质心 Ao 以及在 A 中固定的正交向量 A1>、A2>和 A3>。	m =sm.symbols(‘m’)``Ao = sm.symbols(‘Ao’)``Af = me.ReferenceFrame(‘Af’ )``I = me.outer(Af.x,Af.x)``P = me.Point(‘P’)``A =me.RigidBody(‘A’, Ao, Af, m, (I, P))Af.x、Af.y 和 Af.z 等价于 A1>、A2>和 A3>。	Autolev 在声明后指定质量和惯性，第四和第五个参数用于此。可以传递参数占位符，并使用设置器A.mass = \_和A.inertia = \_来稍后设置它们。更多信息请参考力学/质量。
Frames A``V1> = X1*A1> + X2*A2>	A = me.ReferenceFrame(‘A’ )``v1 = x1*A.x + x2*A.y	更多信息请参考物理/矢量。
Newtonian N	N = me.ReferenceFrame(‘N’ )	SymPy 在声明时没有明确指定参考系为惯性参考系。许多函数如set_ang_vel()需要惯性参考系作为参数。
Particles C	m = sm.symbols(‘m’)``Po = me.Point(‘Po’)``C = me.Particle(‘C’, Po, m)	第二个和第三个参数是用于点和质量。在 Autolev 中，这些在声明后指定。可以传递占位符，并使用设置器（A.point = \_和A.mass = \_）稍后设置它们。
Points P, Q	P = me.Point(‘P’)``Q = me.Point(‘Q’)
Mass B=mB	mB = symbols(‘mB’)``B.mass = mB
Inertia B, I1, I2, I3, I12, I23, I31	I = me.inertia(Bf, i1, i2, i3, i12, i23, i31)``B.inertia = (I, P) B 是一个刚体，Bf 是相关的参考框架，P 是 B 的质心。也可以使用矢量外积形成惯性二元。I = me.outer(N.x, N.x)	更多信息请参考力学 API。
vec> = P_O_Q>/L``vec> = u1*N1> + u2*N2>``Cross(a>, b>)``Dot(a>, b>)``Mag(v>)``Unitvec(v>)``DYAD>> = 3*A1>*A1> + A2>*A2> + 2*A3>*A3>	vec  = (Qo.pos_from(O))/L``vec = u1*N.x + u2*N.y``cross(a, b)``dot(a, b)``v.magnitude()``v.normalize()``dyad = 3*me.outer(a.x ,a.x) + me.outer(a.y, a.y) + 2*me.outer(a.z ,a.z)	更多信息请参考物理/矢量。
P_O_Q> = LA*A1>``P_P_Q> = LA*A1>	Q.point = O.locatenew(‘Qo’, LA*A.x)where A is a reference frame.Q.point = P.point.locatenew(‘Qo ’, LA*A.x)	更多信息请参考运动学 API。这些矢量和运动学函数都应用在Point对象上而不是Particle对象上，因此必须使用.point`来表示粒子。
V_O_N> = u3*N.1> + u4*N.2>``Partials(V_O_N>, u3)	O.set_vel(N, u1*N.x + u2*N.y)``O.partial_velocity(N , u3)	获取器为O.vel(N)。
A_O_N> = 0> 点 O 在参考框架 N 中的加速度。	O.set_acc(N, 0)	获取器为O.acc(N)。
W_B_N> = qB’*B3>	B.set_ang_vel(N, qBd*Bf.z)，这里 Bf 是与体 B 相关的框架。	获取器应为 B.ang_vel_in(N)。
ALF_B_N> =Dt(W_B_N>, N)	B.set_ang_acc(N, diff(B.ang_vel_in(N) )	获取器应为 B.ang_acc_in(N)。
Force_O> = F1*N1> + F2*N2>	Torque_A> = -c*qA’*A3>	在 SymPy 中，需要使用包含所有力和力矩的列表。fL.append((O, f1*N.x + f2*N.y))，这里 fL 是力列表。fl.append((A, -c*qAd*A.z))
A_B = M where M is a matrix and A, B are frames.D = A_B*2 + 1	B.orient(A, 'DCM', M)，其中 M 是 SymPy 矩阵。D = A.dcm(B)*2 + 1
CM(B)	B.masscenter
Mass(A,B,C)	A.mass + B.mass + C.mass
V1pt(A,B,P,Q)	Q.v1pt_theory(P, A, B)	这里假设 P 和 Q 是 Point 对象。记得使用 .point 表示粒子。
V2pts(A,B,P,Q)	Q.v2pt_theory(P, A, B)
A1pt(A,B,P,Q)	Q.a1pt_theory(P, A, B)
A2pts(A,B,P,Q)	Q.a2pt_theory(P, A, B)
Angvel(A,B)	B.ang_vel_in(A)
Simprot(A, B, 1, qA)	B.orient(A, ‘Axis’, qA, A.x)
Gravity(G*N1>)	fL.extend(gravity( g*N.x, P1, P2, ...))	在 SymPy 中，我们必须使用包含形式为 (point, force_vector) 的 forceList（这里是 fL）。这将传递给 KanesMethod 对象的 kanes_equations() 方法。
CM(O,P1,R)	me.functions. center_of_mass(o, p1, r)
Force(P/Q, v>)	fL.append((P, -1*v), (Q, v))
Torque(A/B, v>)	fL.append((A, -1*v), (B, v))
Kindiffs(A, B ...)	KM.kindiffdict()
Momentum(option)	linear_momentum(N, B1, B2 ...)参考框架后跟一个或多个物体 angular_momentum(O, N, B1, B2 ...)点、参考框架后跟一个或多个物体
KE()	kinetic_energy(N, B1, B2 ...)参考框架后跟一个或多个物体
Constrain(...)	velocity_constraints = [...]``u_dependent = [...]``u_auxiliary = [...]这些列表会传递给 KanesMethod 对象。	更多详细信息，请参阅 mechanics/kane 和 kane api.
Fr() FrStar()	KM = KanesMethod(f, q_ind, u_ind, kd_eqs, q_dependent, configura tion_constraints, u_de pendent, velocity_cons traints, acceleration_ constraints, u_auxilia ry)KanesMethod 对象接受一个参考框架，后面跟着多个列表作为参数。(fr, frstar) = KM.kanes_equations(fL, bL)，其中 fL 和 bL 是力和刚体的列表。	更多详细信息，请参阅 mechanics/kane 和 kane api.

数值评估与可视化

Autolev 的 CODE Option() 命令允许生成 Matlab、C 或 Fortran 代码用于数值评估和可视化。Option 可以是 Dynamics、ODE、Nonlinear 或 Algebraic。

可以使用 PyDy 进行动力学的数值评估。可以将 KanesMethod 对象传递给 System 类，以及常数、指定值、初始条件和时间步长。然后可以积分运动方程。使用 matlplotlib 进行绘图。以下是来自PyDy 文档的示例：

from numpy import array, linspace, sin
from pydy.system import System

sys = System(kane,
             constants = {mass: 1.0, stiffness: 1.0,
                          damping: 0.2, gravity: 9.8},
             specifieds = {force: lambda x, t: sin(t)},
             initial_conditions = {position: 0.1, speed:-1.0},
             times = linspace(0.0, 10.0, 1000))

y = sys.integrate()

import matplotlib.pyplot as plt
plt.plot(sys.times, y)
plt.legend((str(position), str(speed)))
plt.show() 

有关 PyDy 可以完成的所有信息，请参阅PyDy 文档。

PyDy 工作流中的工具包括：

• SymPy：SymPy 是 Python 的一个库，用于

  符号计算。它提供计算机代数能力，可以作为独立应用程序、其他应用程序的库，或作为 SymPy Live 或 SymPy Gamma 上的 Web 实时存在。

• NumPy：NumPy 是一个库，用于

  Python 编程语言，添加对大型多维数组和矩阵的支持，以及大量高级数学函数集合，用于操作这些数组。

• SciPy：SciPy 是一个开源的

  用于科学计算和技术计算的 Python 库。SciPy 包含优化、线性代数、积分、插值、特殊函数、FFT、信号和图像处理、ODE 求解器以及其他在科学与工程中常见的任务模块。

• IPython：IPython 是一个命令行外壳

  用于多种编程语言的交互计算，最初为 Python 编程语言开发，提供内省、丰富的媒体、shell 语法、选项完成和历史记录。

• Aesara：Aesara 是

  Python 的数值计算库。在 Aesara 中，使用类似 NumPy 的语法表达计算，并编译以在 CPU 或 GPU 架构上高效运行。

• Cython：Cython 是一个超集

  Python 编程语言，旨在通过大部分 Python 编写的代码提供类似 C 的性能。Cython 是一种编译语言，用于生成 CPython 扩展模块。

• matplotlib：matplotlib 是一个

  用于 Python 编程语言及其数学扩展 NumPy 的绘图库。

通过使用这些科学计算工具，可以编写与 Autolev 生成的 Matlab、C 或 Fortran 代码等效的代码。建议详细了解这些模块，以便理解 Python 科学计算。

链接

SymPy 入门教程

SymPy 文档

SymPy 物理向量文档

SymPy 力学文档

PyDy 文档

使用 Python 进行多体动力学