企业官网建设流程全解析

Sympy SymbolicSystem：统一表达多体系统运动方程的三种标准形式

【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy

本文基于 Sympy 官方解释文档《Symbolic Systems in Physics/Mechanics》，完整讲解sympy.physics.mechanics中的SymbolicSystem类：它以统一的数据格式承载多体动力系统的运动方程，支持显式与隐式三种方程形态，并可选携带系统刚体、载荷、代数约束行号等元信息。读完本文，你将能够把任意多体系统（含约束的 DAE 系统）的运动方程手工输入到SymbolicSystem中，正确区分三种方程形式与参数约定，并通过属性与compute_explicit_form()在形式之间安全转换，为后续对接数值求解器（ODE/DAE 求解代码）做好准备。

一、SymbolicSystem 的定位与设计目标

SymbolicSystem是physics/mechanics模块中用于存放多体动力系统全部相关信息的容器类。按其文档说明（源文件），它的核心定位是：

• 最基本形态：包含系统的运动方程（equations of motion, EOM）；
• 可选扩展信息：系统所受载荷（loads）、组成系统的刚体/质点（bodies），以及用户认为重要的任何附加方程；
• 设计目标：为运动方程提供统一的输出格式，使数值分析代码可以围绕这一格式进行设计。

也就是说，SymbolicSystem并不负责“从模型推导方程”，而是负责把已经得到的方程（无论是解析推导还是其他工具生成）规范地组织起来，供数值积分、线性化等下游环节消费。

需要与同模块中另一个类System（定义）区分开：System是面向建模的高层类，通过add_bodies、add_joints、add_loads、add_kdes等接口搭建模型，再由form_eoms(eom_method=KanesMethod, **kwargs)（实现）调用KanesMethod或LagrangesMethod后端生成运动方程；而SymbolicSystem是面向方程本身的数据容器，两者都通过 模块导出 对外暴露，可以from sympy.physics.mechanics import SymbolicSystem或import sympy.physics.mechanics.system as system两种方式引入。

二、三种运动方程标准形式

SymbolicSystem支持三种等价的输入形式，这是理解整个类的关键。设x为状态量（如[q, u]），t为时间，r为指定的外部输入（exogenous inputs），p为常数，q为广义坐标，u为广义速度：

其中各符号的含义为：

• F_1：合并方程显式形式的右端；
• F_2：合并方程隐式形式的右端；
• F_3：动力学方程隐式形式的右端；
• M_2：合并方程隐式形式的质量矩阵；
• M_3：动力学方程隐式形式的质量矩阵（即广义惯量阵）；
• G：运动学微分方程的右端。

从源码看，__init__依据传入参数自动判定形式（形式判定逻辑）：

1. 若提供了coordinate_derivatives（即G），判为形式 [3]，此时right_hand_side被解释为F_3，mass_matrix被解释为M_3；
2. 否则若提供了mass_matrix，判为形式 [2]，right_hand_side被解释为F_2，mass_matrix为M_2；
3. 否则判为形式 [1]，right_hand_side即F_1。

这一“由参数组合决定解释方式”的机制，是后续所有属性可用性与报错行为的根源。

三、完整实例：单摆的笛卡尔坐标模型手工输入

下面的示例与官方文档 symsystem.rst 保持一致：以笛卡尔坐标描述单摆质量点位置作为广义坐标（而非通常的最小坐标形式），把运动方程手工输入SymbolicSystem。该模型与 lin_pend_nonmin_example 教程中的非最小坐标单摆等价——教程使用q1、q2作为质量点的水平/竖直坐标，而本文示例将其替换为x、y，且参考系相对教程旋转了 90 度（因此重力载荷沿N.x方向）。

3.1 导入与符号初始化

from sympy import atan, symbols, Matrix from sympy.physics.mechanics import (dynamicsymbols, ReferenceFrame, Particle, Point) import sympy.physics.mechanics.system as system from sympy.physics.vector import init_vprinting init_vprinting(pretty_print=False) # 动态符号（时间的函数）：位置 x, y；广义速度 u, v；约束乘子 lam x, y, u, v, lam = dynamicsymbols('x y u v lambda') # 常数符号：质量、摆长、重力加速度 m, l, g = symbols('m l g')

状态向量为(x, y, u, v, lam)共 5 个分量，其中lam是用于约束x**2 + y**2 = l**2的乘子——这正是方程呈微分代数方程（DAE）形式的原因。

3.2 用三种形式写出同一套运动方程

先定义形式 [3] 的动力学方程（M_3 u' = F_3）与形式 [2] 的合并隐式方程（M_2 x' = F_2）：

# 形式 [3]：3x3 质量矩阵 + 右端，覆盖 [u', v', lam'] dyn_implicit_mat = Matrix([[1, 0, -x/m], [0, 1, -y/m], [0, 0, l**2/m]]) dyn_implicit_rhs = Matrix([0, 0, u**2 + v**2 - g*y]) # 形式 [2]：5x5 块矩阵 + 右端，前两行即运动学方程 x' = u, y' = v comb_implicit_mat = Matrix([[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, -x/m], [0, 0, 0, 1, -y/m], [0, 0, 0, 0, l**2/m]]) comb_implicit_rhs = Matrix([u, v, 0, 0, u**2 + v**2 - g*y]) # 形式 [3] 的运动学右端 G kin_explicit_rhs = Matrix([u, v]) # 形式 [1]：对隐式系统做 LU 分解求解，得到显式右端 comb_explicit_rhs = comb_implicit_mat.LUsolve(comb_implicit_rhs)

注意形式 [2] 的comb_implicit_mat具有块对角结构：左上块（前两行）是运动学方程的隐式矩阵（单位阵），右下块（后三行）就是M_3。源码在从形式 [3] 自动构造comb_implicit_mat属性时，正是按此结构做零块拼接（见下文 5.2 节）。

3.3 建立参考系、质点与载荷

虽然方程已手工给出，仍可建立刚体/质点对象，以便将bodies与loads一并传入容器：

theta = atan(x/y) omega = dynamicsymbols('omega') N = ReferenceFrame('N') A = N.orientnew('A', 'Axis', [theta, N.z]) A.set_ang_vel(N, omega * N.z) O = Point('O') O.set_vel(N, 0) P = O.locatenew('P', l * A.x) P.v2pt_theory(O, N, A) # 得到速度 l*omega*A.y Pa = Particle('Pa', P, m) bodies = [Pa] loads = [(P, g * m * N.x)]

载荷的书写约定为：力用(作用点, 力矢量)的元组，力矩用(受作用的参考系, 力矩矢量)的元组（见 类文档字符串）。

3.4 标记代数约束行（alg_con）

本示例的运动方程是DAE形式：DAE 求解器需要知道哪些行是代数方程而非微分方程。这一信息通过alg_con以行索引列表传入SymbolicSystem，有一个重要的索引约定：

• 传入时，行索引对应你输入的矩阵——形式 [3] 下对应dyn_implicit_mat（3 行），本例为alg_con = [2]；
• 访问时，alg_con属性始终对应合并了运动学与动力学的完整方程（5 行），本例为alg_con_full = [4]。

alg_con = [2] # 形式 [3] 下的输入索引 alg_con_full = [4] # 合并形式下的索引

源码在__init__中会自动完成这一偏移：若同时提供了coordinate_derivatives且给出alg_con，则对每个索引加上运动学方程行数（偏移逻辑），因此形式 [3] 传入[2]后，symsystem.alg_con属性返回的就是[4]。

3.5 状态向量与坐标/速度索引

states = (x, y, u, v, lam) coord_idxs = (0, 1) # 前两个状态分量是广义坐标 speed_idxs = (2, 3) # 第 3、4 个分量是广义速度

当第一个参数coord_states传入的是完整状态（而非单独的坐标）时，coord_idxs/speed_idxs告诉SymbolicSystem哪些分量是坐标、哪些是速度；若不传入，coordinates与speeds属性将无法访问并抛出AttributeError。

3.6 创建三个等价的 SymbolicSystem 实例

三种形式分别构造，参数与形式的对应关系是本节重点：

# 形式 [1]：只给显式右端 symsystem1 = system.SymbolicSystem(states, comb_explicit_rhs, alg_con=alg_con_full, bodies=bodies, loads=loads) # 形式 [2]：右端 + 合并质量矩阵 symsystem2 = system.SymbolicSystem(states, comb_implicit_rhs, mass_matrix=comb_implicit_mat, alg_con=alg_con_full, coord_idxs=coord_idxs) # 形式 [3]：动力学右端 + 动力学质量矩阵 + 运动学右端 symsystem3 = system.SymbolicSystem(states, dyn_implicit_rhs, mass_matrix=dyn_implicit_mat, coordinate_derivatives=kin_explicit_rhs, alg_con=alg_con, coord_idxs=coord_idxs, speed_idxs=speed_idxs)

各实例属性的实际输出：

>>> symsystem1.states Matrix([ [x], [y], [u], [v], [lambda]]) >>> symsystem2.coordinates Matrix([ [x], [y]]) >>> symsystem3.speeds Matrix([ [u], [v]]) >>> symsystem1.comb_explicit_rhs Matrix([ [u], [v], [(-g*y + u**2 + v**2)*x/l**2], [(-g*y + u**2 + v**2)*y/l**2], [m*(-g*y + u**2 + v**2)/l**2]]) >>> symsystem2.comb_implicit_rhs Matrix([ [u], [v], [0], [0], [-g*y + u**2 + v**2]]) >>> symsystem2.comb_implicit_mat Matrix([ [1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, -x/m], [0, 0, 0, 1, -y/m], [0, 0, 0, 0, l**2/m]]) >>> symsystem3.dyn_implicit_rhs Matrix([ [0], [0], [-g*y + u**2 + v**2]]) >>> symsystem3.dyn_implicit_mat Matrix([ [1, 0, -x/m], [0, 1, -y/m], [0, 0, l**2/m]]) >>> symsystem3.kin_explicit_rhs Matrix([ [u], [v]]) >>> symsystem1.alg_con [4] >>> symsystem1.bodies (Pa,) >>> symsystem1.loads ((P, g*m*N.x),)

几个值得注意的行为：

• symsystem3.alg_con同样返回[4]——尽管构造时传入的是[2]，源码已自动完成偏移；
• bodies与loads在构造时被统一转换为tuple存储（转换代码），所以symsystem1.bodies输出为(Pa,)；
• 与coord_idxs同理，bodies/loads只有在初始化时指定过，属性才可访问；未指定则访问抛AttributeError。

四、最小示例：单自由度形式 [3] 输入

类文档字符串（示例部分）给出了一个更简洁的用法：单摆用角度theta作广义坐标、omega作广义速度，直接按形式 [3] 输入，坐标与速度通过位置参数分开传入（此时第一个参数是坐标，speeds作为第三个位置参数）：

from sympy import Matrix, sin, symbols from sympy.physics.mechanics import dynamicsymbols, SymbolicSystem l, m, g = symbols('l m g') theta, omega = dynamicsymbols('theta omega') kin_explicit_rhs = Matrix([omega]) # G dyn_implicit_mat = Matrix([l**2 * m]) # M_3 dyn_implicit_rhs = Matrix([-g * l * m * sin(theta)]) # F_3 symsystem = SymbolicSystem([theta], dyn_implicit_rhs, [omega], dyn_implicit_mat)

对比 3.6 节可知两种传参风格的差异：

源码对两种风格的分支处理见 初始化逻辑：传入speeds时，states由coordinates.col_join(speeds)拼成；不传时，第一个参数整体作为states。

五、源码级实现解析

5.1 构造函数与形式判定

SymbolicSystem的完整构造签名为（源码）：

def __init__(self, coord_states, right_hand_side, speeds=None, mass_matrix=None, coordinate_derivatives=None, alg_con=None, output_eqns={}, coord_idxs=None, speed_idxs=None, bodies=None, loads=None):

参数要点：

• coord_states：有序可迭代，时间的函数集合；是否表示“坐标”取决于speeds是否提供；
• right_hand_side：Matrix，其语义（F_1/F_2/F_3）由mass_matrix、coordinate_derivatives是否传入决定；
• speeds：提供后，第一参数被解释为广义坐标；
• mass_matrix：形式 [2]/[3] 的质量矩阵；
• coordinate_derivatives：提供即宣告形式 [3]，内容是运动学方程q' = G的右端；
• alg_con：代数约束行索引；形式 [3] 下按mass_matrix/right_hand_side的行编号解释，内部自动偏移为合并形式的行号；
• output_eqns：字典，键为输出方程名、值为符号表达式，用于追踪额外要跟踪的输出量（单摆示例中可为{PE: m*g*(l+y)}，见 测试）；
• coord_idxs/speed_idxs：当第一参数是完整状态时，指定坐标/速度对应的索引；
• bodies/loads：刚体对象与载荷（力为(point, force)、力矩为(frame, torque)）的可迭代对象，内部转为 tuple。

5.2 属性的自动推导与形式互转

各矩阵属性并非全部需要用户提供，源码按以下规则推导：

• comb_implicit_mat/comb_implicit_rhs：若以形式 [3] 输入，访问这两个属性时源码会用eye(n_kin)、zeros块自动拼接出合并矩阵/右端（实现），与 3.2 节手工写的块矩阵结构完全一致；
• comb_explicit_rhs：不能直接访问，必须先调用compute_explicit_form()计算（实现）。该方法内部即对隐式系统做LUsolve：形式 [3] 下先解M_3 u' = F_3再与kin_explicit_rhs纵向拼接，形式 [2] 下直接解M_2 x' = F_2。源码注释特别提示该计算“potentially take awhile”，因此被设计为显式调用而非惰性自动计算；
• states：始终可用；coordinates/speeds：仅在初始化时能确定对应分量时可用；
• 只读约束：所有属性均为只读 property，尝试赋值（如symsystem.bodies = 42）会抛AttributeError，见 属性只读测试。

5.3 未指定信息的错误行为

测试test_not_specified_errors（源码）系统覆盖了越界访问的行为，值得在集成时记住：

• 形式 [1] 的实例访问comb_implicit_mat、dyn_implicit_rhs、kin_explicit_rhs等均未定义，抛AttributeError；且因为comb_explicit_rhs已存在，compute_explicit_form()也会报错（守卫逻辑）；
• 形式 [2] 的实例访问dyn_implicit_mat、kin_explicit_rhs报错——源码不会从合并矩阵中反拆出动力学块；
• 仅传状态而未传coord_idxs/speed_idxs时访问coordinates/speeds报错；
• 未传bodies/loads时访问对应属性报错。

5.4 动态符号与常数符号的自动提取

SymbolicSystem还提供两个辅助方法，便于在对接数值求解器时确定变量清单：

• dynamic_symbols()（实现）：扫描运动方程表达式（显式右端，或隐式矩阵与右端的全部元素）中的Dynamicsymbol，再并入状态量，返回时间依赖符号的 tuple；
• constant_symbols()（实现）：取同样的表达式集合的free_symbols，剔除时间符号t，返回常数符号的 tuple。

对本文示例，两者分别得到{x, y, u, v, lam}与{m, l, g}，与 测试断言 一致。

六、API 速查表

构造函数参数

主要属性与方法

七、小结

SymbolicSystem用一套紧凑的接口解决了“运动方程以什么格式交给数值代码”的问题：三种形式覆盖了显式/隐式、运动学与动力学合并/分离的典型场景，alg_con索引约定让 DAE 求解器能正确识别代数约束，bodies/loads与output_eqns则为后续线性化、输出追踪预留了空间。其属性设计遵循“传入什么、自动推导什么、越界访问即报错”的明确边界，配合 test_system.py 中三种形式等价性的完整断言（三种形式构造的实例在comb_implicit_*与compute_explicit_form()之后输出一致），可作为对接第三方数值求解器时行为契约的可靠依据。实际使用时，建议从 3.6 节的最小构造方式起步：确定方程属于哪种形式、核对alg_con的行号基准、区分“坐标+速度”与“完整状态”两种传参风格，即可获得一个可供数值分析代码消费的标准化系统描述。

【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy