SymPy 力学模块实战：用拉格朗日方法建模纯滚动圆盘并推导运动方程

SymPy 力学模块实战用拉格朗日方法建模纯滚动圆盘并推导运动方程【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy在 SymPy 的sympy.physics.mechanics模块中拉格朗日方法Lagranges Method提供了一条从系统能量出发、自动推导二阶运动方程的高效路径。本文以「无限薄圆盘在水平地面上无滑动纯滚动」这一经典非完整约束系统为例完整演示如何用Lagrangian与LagrangesMethod从接触点向上构建运动学、组装惯性张量、计算拉格朗日量并最终求出广义加速度的显式右端项right-hand side。读完本文你将掌握 SymPy 力学模块中拉格朗日方法的标准工作流并理解其底层质量矩阵、强迫向量与rhs()的求解机制。本文对应的原始教程为 doc/src/tutorials/physics/mechanics/rollingdisc_example_lagrange.rst其完整示例位于同一目录下的三篇姊妹教程Kane 方法版、带约束力的 Kane 方法版 与总览页它们从不同形式主义formalism出发建模同一个系统便于横向对比。sympy 滚动圆盘示意图一、问题设定从接触点向上建立模型滚动圆盘模型假设圆盘无限薄与地面仅有一个接触点且在地面上无滑动滚动。总览文档 rollingdisc_example.rst 明确指出官方会在该教程中用三种不同的方式建模同一个系统以展示sympy.physics.mechanics模块的多种功能Kane 方法版引入广义速度generalized speeds带约束力的 Kane 方法版额外引入辅助速度与约束力把非约束力constraint forces显式表达出来拉格朗日方法版本文主角从接触点向上构建无需引入广义速度仅用 3 个构型变量configuration variables和 3 个速度变量、加上圆盘质量、半径与局部重力即可完整描述该系统。在拉格朗日方法中圆盘的完整状态由三个广义坐标描述q1绕地面法线N.z的旋转角q2圆盘相对竖直方向的倾斜角lean angleq3圆盘绕自身对称轴的滚动角。这三个坐标通过一个3-1-2即 Z-X-Y欧拉角序列连接三个中间参考系从而避免引入额外的广义速度这正是原文档强调的“从接触点向上formed from the contact point up”建模思路的核心优势。二、符号定义与参考系frame搭建首先导入所需符号并开启 SymPy 力学模块的简洁打印模式 from sympy import symbols, cos, sin from sympy.physics.mechanics import * mechanics_printing(pretty_printFalse) q1, q2, q3 dynamicsymbols(q1 q2 q3) q1d, q2d, q3d dynamicsymbols(q1 q2 q3, 1) r, m, g symbols(r m g)这里的关键点dynamicsymbols(q1 q2 q3)创建的是时间的函数q1(t)、q2(t)、q3(t)第二个参数1表示对时间求一阶导即q1d对应Derivative(q1(t), t)后续代码会隐式利用这一点完成速度表达式的自动求导。r, m, g是常量符号分别代表圆盘半径、质量与局部重力加速度。mechanics_printing(pretty_printFalse)让后续输出以更适合复制粘贴的线性文本显示。运动学由一系列简单旋转simple rotation构成每次简单旋转都创建一个新参考系下一次旋转由新参考系的基向量定义。本示例采用 3-1-2 系列旋转即 Z、X、Y 序列并且角速度定义在第二个参考系倾斜系 L的基上这正是需要定义中间参考系而不是直接使用 body-three orientation 的原因 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——————惯性参考系地面/竖直YNN.zq1绕竖直轴的水平方位角LYY.xq2倾斜参考系lean frameRLL.yq3圆盘本体参考系滚动角关于 3-1-2 序列的选择Kane 方法版教程给出了同样的解释由于最终角速度需要用第二个参考系倾斜系L的基向量来表达因此必须显式构建中间参考系。三、平动运动学从接触点到质心接下来是平动运动学。首先创建一个在N中速度为零的点C它就是圆盘与地面的接触点然后从接触点指向圆盘质心形成位置向量最后利用v2pt_theory自动求取质心的速度 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这段代码的要点C.set_vel(N, 0)声明接触点在惯性系中静止。由于圆盘无滑动滚动接触点速度恒为零这是整个运动学推导的起点。C.locatenew(Dmc, r * L.z)沿倾斜系L的z轴偏移半径r得到质心Dmc。注意这里沿L.z而不是R.z因为圆盘关于倾斜系对称这样会得到更简洁的惯性表达式Kane 版教程明确指出“圆盘的惯量在滚动过程中于倾斜系内保持不变”。v2pt_theory(C, N, R)是两点速度理论two-point theoremv_Dmc_N v_C_N w_R_N × r_C-Dmc它需要提供固定点C、参考系N以及包含两点的刚体参考系R。输出的速度向量中同时出现了q1、q2、q3说明圆盘的纯滚动约束已经通过运动学结构自动编码进来——这正是“从接触点向上建模”的收益不需要显式写出滚动约束方程。四、组装惯性张量inertia dyadic与刚体形成惯性张量。圆盘视为无限薄绕直径轴的转动惯量为m*r**2/4绕对称轴L.z的转动惯量为m*r**2/2 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))关键点inertia(frame, ixx, iyy, izz)返回一个惯性并矢dyadic分量定义在倾斜系L中。由于圆盘的对称性交叉惯性积为零因此只需提供三个主轴分量。RigidBody(BodyD, Dmc, R, m, (I, Dmc))将名称、质心点Dmc、本体参考系R、质量m与惯性并矢 质心打包成一个刚体对象。注意质心位置取L.z方向、而本体参考系用R两者可以不同——这正是前面刻意让惯性在L中表达的原因。五、势能与拉格朗日量接着设置势能并计算滚动圆盘的拉格朗日量 BodyD.potential_energy - m * g * r * cos(q2) Lag Lagrangian(N, BodyD)质心相对接触点的高度为r*cos(q2)因此势能V m*g*r*cos(q2)代码中的负号来自 SymPy 约定取零势能参考点为圆盘竖直时质心所处高度向上为负方向。Lagrangian(frame, *body)定义在 sympy/physics/mechanics/functions.py返回系统中所有Particle与RigidBody的动能减去势能即L T - V是一个标量。它以N作为计算动能的参考惯性系。该函数会遍历传入的刚体利用其质量、质心速度与本体系角速度自动累加平动动能1/2*m*v**2与转动动能1/2*w·(I·w)。六、生成拉格朗日方程并求广义加速度运动方程通过初始化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]])输出结果中前三行是平凡的运动学关系q_i q_i速度定义后三行是广义加速度的显式表达式即把二阶方程整理成x f(x)的一阶显式形式后的右端项。例如q2 -4*g*sin(q2)/(5*r) ...展示了倾斜角加速度中重力项与陀螺/科氏耦合项的分离。6.1LagrangesMethod的构造参数从源码 sympy/physics/mechanics/lagrange.py 可以看到LagrangesMethod构造函数的完整签名是LagrangesMethod(Lagrangian, qs, forcelistNone, bodiesNone, frameNone, hol_coneqsNone, nonhol_coneqsNone)各参数含义摘自类文档字符串参数类型说明LagrangianSympifyable标量表达式为系统动能与势能之和的函数q 与 q 的函数qs时间函数组成的可迭代对象多体系统的广义坐标 qhol_coneqsExpr 可迭代对象可选完整约束holonomic constraint残差nonhol_coneqsExpr 可迭代对象可选非完整约束nonholonomic constraint残差forcelist可迭代对象可选施加在点上的力 / 施加在参考系上的力矩的 (Point, Vector) 或 (ReferenceFrame, Vector) 元组只应包含非保守力保守力已在拉格朗日量中体现bodies可迭代对象可选组成多体系统的Particle、RigidBody或Body对象frameReferenceFrame可选与构造拉格朗日量所用一致的惯性参考系仅在提供forcelist时必须如果提供约束方程拉格朗日乘子Lagrange multipliers会被自动生成其个数与约束方程数相等见源码中self.lam_vec Matrix(dynamicsymbols(lam1: str(m 1)))。6.2form_lagranges_equations的底层实现form_lagranges_equationslagrange.py内部把运动方程表示为四个项的组合EOM term1 - term2 - term3 - term4 0其中term1 d/dt(∂L/∂q)即对拉格朗日量求速度的雅可比后再对时间求导源码self._L.jacobian(qds).diff(t).Tterm2 ∂L/∂q源码self._L.jacobian(self.q).Tterm3若存在约束为约束雅可比的转置乘拉格朗日乘子Cᵀ·λ源码self._term3 self.lam_coeffs.T * self.lam_vec否则为零矩阵term4非保守广义力的贡献源码中通过v.diff(qd, N).dot(f)对每个广义速度计算广义力。最终eom term1 - term2 - term4 - term3并通过linear_eq_to_matrix从「不含乘子」的部分提取动态质量矩阵Md与强迫向量Fd满足Md*q gd 0。有约束时运动方程形如Md*q Cᵀ*λ gd 0。6.3rhs()如何求解广义加速度rhs()定义在基类 sympy/physics/mechanics/method.py其返回值是完整一阶形式的右端项rhs(u, q, t) Inv(M) F具体地当inv_method未指定时使用self.mass_matrix_full.LUsolve(self.forcing_full)LU 分解求解否则使用指定的矩阵求逆方法。mass_matrix_full与forcing_full属性定义在LagrangesMethod中lagrange.py无约束时mass_matrix_full [[I, 0], [0, Md]]forcing_full [q, Fd]有约束时扩展为含乘子的分块形式[[I, 0, 0], [0, Md, -Cᵀ], [0, C, 0]]与[q, Fd, Fc]。在滚动圆盘示例中不存在显式约束方程因此结果矩阵即q拼上由Md⁻¹·Fd得到的三个广义加速度与教程输出完全吻合。另外注意LagrangesMethod还提供mass_matrix、forcing、to_linearizer()、linearize()与solve_multipliers()等属性与方法可用于后续的线性化与拉格朗日乘子求解。七、与其他建模方法的对比同一个系统在官方教程中被建模了三次理解差异有助于在实际项目中选型方法状态变量约束处理输出形式参考教程拉格朗日方法3 个 q 3 个 q纯滚动约束由“从接触点向上”的运动学自动满足二阶 ODE 一阶显式 rhsrollingdisc_example_lagrange.rstKane 方法3 个 q 3 个广义速度 u通过运动微分方程kd关联 q 与 u一阶显式 rhsrollingdisc_example_kane.rstKane 方法 约束力3 个 q 3 个 u 3 个辅助速度显式引入 u4/u5/u6 与 f1/f2/f3将约束力暴露出来一阶 rhs 辅助方程rollingdisc_example_kane_constraints.rst对比 Kane 版的输出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]])可以看到u2的表达式-2*u1*u3/3与拉格朗日版中q2 -2*(2*tan(q2)*q1 3*q3/cos(q2))*q2在代入运动微分方程后是等价的前者以广义速度 u 表达后者以 q 表达。两种形式主义给出数学上等价、形式上不同的运动方程这正是官方教程刻意安排的对比目的。若希望把约束力也显式求解出来可参考 Kane 约束力版它引入垂直于地面、沿滚动路径、以及地面内垂直方向的三个辅助速度u4,u5,u6其速度恒为零并对应施加三个约束力f1,f2,f3通过KanesMethod(..., u_auxiliary[u4, u5, u6])后从KM.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]])八、后续步骤与延伸阅读得到rhs()的一阶显式方程后你可以用l.rhs()的输出直接做数值积分配合lambdify进行仿真或传递给Linearizer做局部线性化l.to_linearizer()/l.linearize()对有约束的系统调用l.solve_multipliers(op_point...)求出特定工作点下的拉格朗日乘子对比阅读本仓库中的源码实现 lagrange.py、functions.py 与 method.py以及sympy/physics/mechanics/tests/目录下针对LagrangesMethod的单元测试进一步验证上述推导。本篇教程是 Mechanics Tutorials 系列的一部分该系列还包含 duffing 振子、多自由度完整约束系统、四连杆机构、自行车模型等案例读者可在掌握滚动圆盘流程后进一步拓展到更复杂的多体系统。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考