SymPy Vector 模块：向量与 ReferenceFrame 参考系的数学原理及代码实现详解

SymPy Vector 模块向量与 ReferenceFrame 参考系的数学原理及代码实现详解【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy在刚体动力学与多体系统建模中向量 参考系是最底层的描述语言。本文基于 SymPy 的sympy.physics.vector模块官方说明文档系统讲解向量代数加法、点积、叉积、基分解、向量微积分参考系内导数、方向余弦矩阵 DCM、多参考系求导的数学定义并结合 vector.py 与 frame.py 源码说明 SymPy 是如何把这些数学关系编码为可符号计算的Vector与ReferenceFrame对象的。读完后你将能够独立完成参考系定向Axis/Body/Space/Quaternion、多参考系下的时间与变量导数运算以及向量矩阵化等符号动力学基础操作。一、向量的定义与向量代数1.1 相等性与基本运算向量是具有大小模和方向的几何对象。文档给出的基本约定是两个向量相等当且仅当它们具有相同的大小和方向。向量代数支持三类运算向量加法基于平行四边形法则满足交换律与结合律a b b a(a b) c a (b c)。标量乘法向量与标量相乘结果方向不变、模按标量缩放乘以 -1 等价于绕垂直平面内任意轴旋转 180 度。向量乘法sympy.physics.vector实现了三种——点积、叉积与外积外积留待惯性章节讨论。点积将两个向量映射为标量a·b |a||b|cos(θ)其中 θ 为两向量夹角。两个互相垂直的向量点积为零点积满足交换律a·b b·a。叉积将两个向量映射为向量a × b c其中c同时垂直于a与b模为|c| |a||b|sin(θ)方向由右手定则确定。叉积的关键性质不满足交换律a × b ≠ b × a且a × b -b × a不满足结合律(a × b) × c ≠ a × (b × c)两平行向量叉积为零。其他常用恒等式文档原文公式的完整继承α(a b) αa αb a·(b c) a·b a·c a × (b c) a × b a × c (a × b)·c 为标量三重积a × (b·c) 无意义不能叉积标量 (a × b)·c a·(b × c) (b × c)·a (c × a)·b (a × b) × c b(a·c) - a(b·c) a × (b × c) b(a·c) - c(a·b)1.2 基Basis表示正交基与度量数取三个不共面的单位向量n_x, n_y, n_z任何向量都可表示为a a_x n_x a_y n_y a_z n_z其中n_x, n_y, n_z称为基basisa_x, a_y, a_z称为度量数measure numbers。若三个单位向量两两正交则称为正交基orthonormal basis通常取右手系。在基表示下两向量a与b相等当且仅当a_xb_x, a_yb_y, a_zb_z加法a b (a_xb_x)n_x (a_yb_y)n_y (a_zb_z)n_z标量乘αb αb_x n_x αb_y n_y αb_z n_z点积a·b a_x b_x a_y b_y a_z b_z叉积a × b等于以n_x, n_y, n_z为第一行、a与b分量为后两行的 3×3 行列式标量三重积(a × b)·c等于a, b, c分量构成的行列式任意向量可写为a (a·n_x)n_x (a·n_y)n_y (a·n_z)n_z用点积投影提取度量数。文档中的数值示例α为标量a n_x 5 n_y b n_y α n_z a b n_x 6 n_y α n_z a·b 5 a·n_y 5 a·n_z 0 a × b 5α n_x - α n_y n_z b × a -5α n_x α n_y - n_z二、向量微积分参考系、导数与方向余弦矩阵2.1 为什么必须引入参考系对运动物体做向量微积分时同一个向量的导数取决于观察者的参考系。文档用经典例子说明车厢内两人静坐相对速度为零但车外观察者看到两人都有速度。严格地说参考系是一个虚拟的观察平台。若向量a在参考系N中的任何性质模、方向观测上都不变则称a固定在N中。每个参考系通常配有一组固定的正交单位基向量如N系的n_x, n_y, n_z。由此导数必须把参考系写进符号^A d e / dθ ≠ 0在 A 系中对 e 关于 θ 求导而^B d e / dθ 0e 固定在 B 系中。参考系特定导数的性质包括^A d/dt (a b) ^A da/dt ^A db/dt ^A d/dt (γa) dγ/dt · a γ · ^A da/dt ^A d/dt (a × b) ^A da/dt × b a × ^A db/dt2.2 方向余弦矩阵DCM两个参考系的基向量之间的关系用方向余弦矩阵表达[â_x] [^A C^B] [b̂_x] [â_y] [ ] [b̂_y] [â_z] [ ] [b̂_z]当两参考系初始对齐、再绕与某基向量重合的轴旋转后称两系之间为简单旋转simple rotation。绕 Z 轴旋转 θ 的 DCM 为^A C^B [[ cos θ, -sin θ, 0], [ sin θ, cos θ, 0], [ 0, 0, 1]]绕 X 轴与 Y 轴的简单旋转 DCM 分别为绕 X 轴: [[1, 0, 0 ], [0, cos θ, -sin θ ], [0, sin θ, cos θ ]] 绕 Y 轴: [[cos θ, 0, sin θ ], [ 0, 1, 0 ], [-sin θ, 0, cos θ ]]正方向旋转按右手定则定义。DCM 还可以直接用两组基向量的点积定义矩阵第 i 行第 j 列元素是â_i · b̂_j。DCM 是正交矩阵满足^A C^B (^B C^A)^{-1} (^B C^A)^T。换基示例设 A、B 两系之间发生绕 Z 轴的简单旋转 θ定义a â_x â_y â_z、b b̂_x b̂_y b̂_z。要把b用 A 系表示就分别取b与â_x, â_y, â_z的点积作为度量数得到b (cosθ - sinθ)â_x (sinθ cosθ)â_y â_z反向将a用 B 系表示则得a (cosθ sinθ)b̂_x (-sinθ cosθ)b̂_y b̂_z。2.3 多参考系下的导数与文档算例要把向量b b_x b̂_x b_y b̂_y b_z b̂_z在参考系 A 中求导必须先把它用 A 系的基表示再对各度量数求导^A db/dx d(b·â_x)/dx â_x d(b·â_y)/dx â_y d(b·â_z)/dx â_z文档的完整算例两个刚体各附一个参考系θ 与 x 都是时间的函数。定义c x b̂_x l b̂_y在 B 系中的定义θ 为绕 Y 轴的旋转角。在 B 系中l 为常量^B dc/dt dx/dt b̂_x ẋ b̂_x在 A 系中先写出绕 Y 轴的 DCM^A C^B将 c 用 A 系表示后再求导^A dc/dt (-θ̇ sinθ · x cosθ · ẋ) â_x (θ̇ cosθ · x sinθ · ẋ) â_z注意这是A 系中的时间导数且以 A 系基表达。同一个导数也可以改用 B 系基表达^A dc/dt ẋ b̂_x - θ̇ x b̂_z文档原文此处 θx 实为 θ̇x 的排版简写。两种表达等价但后者形式简洁得多。文档特别强调用更复杂的基形式定义向量会显著拖慢运动方程的建立过程、并使表达式膨胀到无法在屏幕上显示——这正是符号动力学工具必须内置换系表达能力的原因。三、代码实现从ReferenceFrame到向量运算3.1 创建参考系与访问基向量文档指出在sympy.physics.vector中开始任何问题第一步就是定义参考系并且必须先创建参考系才能访问基向量 from sympy.physics.vector import * N ReferenceFrame(N) N.x N.x N.y N.y N.z N.z从源码看ReferenceFrame.__init__见 frame.py除了存储名称还支持indices用方括号下标如O[1]访问基向量与latexs自定义 LaTeX 输出两个可选参数此外N.x, N.y, N.z本质上是不可变Vector其在本系中的度量数分别为[1,0,0], [0,1,0], [0,0,1]与文档 How Vectors are Coded 一节的描述一致。3.2 向量代数操作与接口约定基本代数 N.x N.x True N.x N.y False N.x N.y N.x N.y 2 * N.x N.y 2*N.x N.y注意标量不能加到向量上N.x 5会报错。向量分量中可以使用 SymPy 的Symbol from sympy import Symbol, symbols x Symbol(x) x * N.x x*N.x x*(N.x N.y) x*N.x x*N.y向量乘法提供三个层级的接口运算符、^、方法dot/cross与函数dot(...)/cross(...)。官方推荐函数接口因为运算符的优先级不符合数学直觉混用时必须小心括号函数实现位于 functions.py N.x.dot(N.x) 1 N.x.dot(N.y) 0 dot(N.x, N.x) 1 dot(N.x, N.y) 0 N.x.cross(N.x) 0 N.x.cross(N.z) - N.y cross(N.x, N.y) N.z cross(N.x, (N.y N.z)) - N.y N.z归一化与取模 (N.x N.y).normalize() sqrt(2)/2*N.x sqrt(2)/2*N.y (N.x N.y).magnitude() sqrt(2)矩阵化输出——注意矩阵形式不含参考系信息必须显式提供参考系才能提取度量数 (x * N.x 2 * x * N.y 3 * x * N.z).to_matrix(N) Matrix([ [ x], [2*x], [3*x]])源码印证Vector.to_matrix见 vector.py的实现正是对目标参考系的三个基向量逐一做点积后 reshape 成 3×1 矩阵与上文 1.2 节的投影公式a (a·n_x)n_x ...完全对应。而Vector.dot的实现vector.py在跨参考系点积时会取出两参考系间的 DCM 做矩阵乘法这正是DCM 由基向量点积定义的代码形态Vector.crossvector.py则因 SymPy 的Matrix不能容纳Vector元素专门写了一个 3×3 行列式辅助函数_det来计算叉积的行列式定义。3.3 单参考系求导diff方法与dt简写 (x * N.x N.y).diff(x, N) N.x文档特别警告SymPy 通用的diff函数目前不适用于sympy.physics.vector的Vector因为向量求导除了对谁求导还必须指定在哪个参考系中求导通用diff的签名无法容纳这一参数。请一律使用Vector.diff(var, frame)。diff的第三个参数var_in_dcm默认True控制算法若请求已知求导变量不出现在方向余弦矩阵中可传False跳过换系重构以提升性能见 vector.py 中的分支逻辑——若frame.dcm(component_frame).diff(var)为零矩阵就直接对度量数求导否则先把分量换表到目标参考系再求导。时间导数可用简写方法dt内部调用全局函数time_derivative见 vector.py (B.y*q2 B.z).dt(N) (-q1 q2)*B.y q2*q1*B.z重要行为导数的输出向量保持在输入向量所在的参考系中表达即使输入向量混合了多个参考系的基向量 (B.y*q2 B.z).diff(q2, N) B.y (B.y*q2 B.z q2*N.x).diff(q2, N) N.x B.y3.4 定向orientAxis / Body / Space / Quaternion两个未定向的参考系之间可以相加但不能做向量乘法——必须先定义朝向关系 A ReferenceFrame(A) A.x N.x N.x A.xReferenceFrame.orient方法提供朝向定义Axis 类型表示绕指定轴做简单旋转 A.orient(N, Axis, [x, N.y])上式表示 A 相对 N 绕 Y 轴旋转角度 x。任意时刻可用dcm方法查看两系间的方向余弦矩阵A.dcm(N)返回^A C^N。更复杂的旋转类型包括 Body 旋转、Space 旋转、四元数Quaternion / Euler 参数以及任意轴旋转。Body 旋转等价于连续三次简单旋转每次都绕新旋转后参考系的基向量。文档给出了三次 Axis 依次定向与Body 三步定向结果完全一致的验证示例 N ReferenceFrame(N) Bp ReferenceFrame(Bp) Bpp ReferenceFrame(Bpp) B ReferenceFrame(B) q1,q2,q3 symbols(q1 q2 q3) Bpp.orient(N,Axis, [q1, N.x]) Bp.orient(Bpp,Axis, [q2, Bpp.y]) B.orient(Bp,Axis, [q3, Bp.z]) N.dcm(B) Matrix([ [ cos(q2)*cos(q3), -sin(q3)*cos(q2), sin(q2)], [sin(q1)*sin(q2)*cos(q3) sin(q3)*cos(q1), -sin(q1)*sin(q2)*sin(q3) cos(q1)*cos(q3), -sin(q1)*cos(q2)], [sin(q1)*sin(q3) - sin(q2)*cos(q1)*cos(q3), sin(q1)*cos(q3) sin(q2)*sin(q3)*cos(q1), cos(q1)*cos(q2)]]) B.orient(N,Body,[q1,q2,q3],XYZ) N.dcm(B) # 输出与上一次 N.dcm(B) 完全相同Space 定向与 Body 定向类似但旋转施加的次序相反从母系出发旋转序列可以是 XYZ、YZX、ZXZ、YXY 等两轴或三轴序列。关键约束相邻的简单旋转必须绕不同的轴——ZZX 这样的序列不能完整定向一个三维基。四元数与任意轴旋转详见orient/orientnew方法帮助文档或 Kane, 1983 参考书。源码中这些接口都有对应的专用方法orient_axis、orient_explicit、orient_dcm、orient_body_fixed、orient_space_fixed、orient_quaternion见 frame.pyorient是它们的统一分发入口。orientnew则把创建 定向一步完成封装了orient的全部能力 C N.orientnew(C, Axis, [q1, N.x])3.5dynamicsymbols与向量打印多参考系微分运算前文档引入了dynamicsymbols——创建时间的未定义函数的快捷方式。这类dynamicsymbol对时间的导数会自动显示为带撇号的记号 from sympy import diff q1, q2, q3 dynamicsymbols(q1 q2 q3) diff(q1, Symbol(t)) Derivative(q1(t), t) q1 q1(t) q1d diff(q1, Symbol(t)) vprint(q1) q1 vprint(q1d) q1非交互会话用vprint交互会话用init_vprinting此外 SymPy 的vprint、pprint、latex也各有对应版本vprint、vpprint、vlatex from sympy.physics.vector import init_vprinting init_vprinting(pretty_printFalse) q1 q1 q1d q1文档的约定在sympy.physics.vector中任何随时间变化的量——坐标、速度、力——都应当用 dynamicsymbol 表示其主要用途是速度speeds与广义坐标。3.6 结合 dynamicsymbol 的完整求导示例 N ReferenceFrame(N) B N.orientnew(B, Axis, [q1, N.x]) (B.y*q2 B.z).diff(q2, N) B.y (B.y*q2 B.z).dt(N) (-q1 q2)*B.y q2*q1*B.z (B.y*q2 B.z q2*N.x).diff(q2, N) N.x B.y对比可见diff(q2, N)只求关于 q2 的偏导而dt(N)是完整时间导数——B 系的基向量本身随 q1(t) 旋转因此dt结果中出现了q1项这与 2.3 节手写算例中^A dc/dt的推导结构一致。四、Vector 的内部编码结构进阶阅读文档最后专门用一节说明内部实现供想理解模块机制的读者参考。核心要点args属性是向量的全部主要信息它是一个列表元素个数等于该向量分量中涉及的不同ReferenceFrame的数量——若向量含 A、B 两系基向量则args长度为 2再加 C 系则为 3。每个元素是一个二元组第一项是 SymPyMatrix存该参考系下各基向量的度量数第二项是对应的ReferenceFrame把度量数与该参考系关联起来。构造时同名参考系的度量数矩阵会先相加合并见 vector.py 中d[inp[1]] inp[0]的合并逻辑全零分量会被丢弃。ReferenceFrame自身存储创建时的name、以及orient/orientnew时建立的 DCM。DCM 用 SymPyMatrix表示存放在以ReferenceFrame为键、Matrix为值的字典中且双向设置把 A 定向到 N 时A 的定向字典加入 N 及其矩阵同时 N 的定向字典也加入 A 及其矩阵该 DCM 是前者的转置。从源码看见 frame.pydcm()方法还维护了双向的_dcm_cacheoutdcm.T会被缓存到对方参考系中从而保证^A C^B (^B C^A)^T这一正交性以 O(1) 查询成立若 A 已定向到 N 之后再次调用orient源码还会清空相关缓存并检查引用循环。在ReferenceFrame创建之前代码中不存在任何向量参考系创建时x/y/z属性即成为本系度量数为[1,0,0]、[0,1,0]、[0,0,1]的不可变向量此后一切新向量都由基向量的代数运算派生而来且天然允许跨参考系混合分量。相关测试位于 test_vector.py、test_frame.py 与 test_functions.py可用于对照本文示例验证各接口行为。五、小结与使用建议先建参考系再做向量ReferenceFrame(N)是一切操作的起点基向量只从参考系上访问。函数接口优先dot(vec1, vec2)、cross(vec1, vec2)比运算符与方法更安全能避免运算符优先级陷阱。求导必须指定参考系用Vector.diff(var, frame)时间导数用dt不要用 SymPy 通用diff已知求导变量不在 DCM 中时可传var_in_dcmFalse提速。换系表达是表达简洁性的来源同一导数在不同参考系基下表达复杂度差异巨大符号建模时应始终选择最简基形式。时变量一律用 dynamicsymbol坐标、速度、力都用dynamicsymbols生成并配合vprint/vpprint/vlatex获得紧凑的可读输出。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考