ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

Poly solvers:SymPy 底层线性方程组求解模块深度解析

Poly solvers:SymPy 底层线性方程组求解模块深度解析 Poly solversSymPy 底层线性方程组求解模块深度解析【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy本文以 SymPy 官方文档 doc/src/modules/polys/solvers.rst 为核心系统讲解sympy.polys.solvers模块的设计初衷、全部公开/内部 API、核心算法流程及其在linsolve、solve、Risch 积分等高层功能中的实际调用场景帮助读者理解 SymPy 是如何在多项式环与系数域上高效求解线性方程组的。模块定位面向内部的高效线性求解器SymPy 是一个用纯 Python 编写的计算机代数系统CAS。在它的代码库中sympy.polys.solvers是一个专门负责“求解线性方程组”的低层low-level模块。官方文档开宗明义地指出This module provides functions for solving systems of linear equations that are used internally in sympy.也就是说这个模块的定位非常清晰它不是一个面向终端用户的高层 API而是 SymPy 内部大量算法共用的基础设施。它的核心优势在于输入是多项式环元素PolyElement而非普通表达式Expr系数直接落在Domain如整数环ZZ、有理数域QQ、参数有理函数域等上核心运算全部在DomainElement级别完成相比以Expr为载体的符号运算对于最常见输入有着更高的效率上层函数如linsolve、solve会把普通表达式先转换到多项式环再委托给本模块完成实际求解。从源码结构看sympy/polys/solvers.py模块公开了 3 个入口函数和 2 个下划线内部函数官方文档通过autofunction指令逐一收录函数可见性职责solve_lin_sys(eqs, ring, _rawTrue)公开偏内部求解一个来自多项式环的线性方程组eqs_to_matrix(eqs_coeffs, eqs_rhs, gens, domain)公开内部把 dict 形式的低层方程转换为增广矩阵DomainMatrixsympy_eqs_to_ring(eqs, symbols)公开内部把Expr形式的方程组转换为PolyRing上的元素_solve_lin_sys(eqs_coeffs, eqs_rhs, ring)内部将系统拆分为连通分量并逐一求解_solve_lin_sys_component(eqs_coeffs, eqs_rhs, ring)内部对单个连通分量做高斯-若尔当消元并回代此外模块还定义了一个异常类PolyNonlinearError专门用于在solve_lin_sys遇到非线性项时抛出。一、solve_lin_sys模块的主入口solve_lin_sys(eqs, ring, _rawTrue)是模块对外的主函数它的功能是求解一个由PolyRing多项式环元素给出的线性方程组其中每个方程都被视为等于零。参数说明参数类型含义eqslist[PolyElement]待求解的线性方程作为某个PolynomialRing的元素隐式等于 0ringPolynomialRing方程所在的多项式环环的生成元generators就是要求解的未知量环的域domain就是方程系数的系数域_rawbool若为False返回字典的键和值均为Expr类型并且会去掉键上的“单位”系数如1.0x若为True默认返回低层多项式类型如PolyElement: PythonRational返回值系统无解时返回None_rawFalse时返回dict[Symbol, Expr]_rawTrue时返回dict[Symbol, DomainElement]。官方示例 from sympy import symbols from sympy.polys.solvers import solve_lin_sys, sympy_eqs_to_ring x, y symbols(x, y) eqs [x - y, x y - 2] eqs_ring, ring sympy_eqs_to_ring(eqs, [x, y]) solve_lin_sys(eqs_ring, ring) {y: 1, x: 1}传入_rawFalse时结果相同但键会从低层多项式类型变为Expr solve_lin_sys(eqs_ring, ring, _rawFalse) {x: 1, y: 1}底层处理流程从源码实现sympy/polys/solvers.py#L133-L247可以看到solve_lin_sys内部做了如下预处理断言环的域是域Fieldassert ring.domain.is_Field。这是求解能正常进行的先决条件——只有在域上才能做除法消元时需要除以主元。把每个PolyElement转成 dict键为单项式monomial值为系数并从 dict 中弹出常数项作为右端项eqs_rhs剩余部分作为左端系数eqs_coeffs。逐项检查非线性对每个非零单项式检查sum(monom) ! 1一旦出现次数和不为 1 的单项式即非线性项立即抛出PolyNonlinearError(Nonlinear term encountered in solve_lin_sys)。剔除恒等式并检测矛盾若某方程无任何未知量eq_coeffs为空当常数项也为 0 时该方程是恒等式直接跳过当常数项非 0 时说明系统矛盾直接返回None。将预处理后的(eqs_coeffs, eqs_rhs)交给内部函数_solve_lin_sys。若_rawFalse最后还会把结果中的每个Symbol/DomainElement通过as_expr()或ring.domain.to_sympy(x)转换回Expr并处理类似1.0x的单位系数k.as_coeff_Mul()后除以系数。注意虽然solve_lin_sys是一个公开函数但官方文档明确提示它的接口“不一定方便”建议普通用户改用sympy.solvers.solveset.linsolve——该函数内部正是调用本函数实现的。二、sympy_eqs_to_ring从Expr到PolyRing的桥接高层函数如solve期望输入是Expr但solve_lin_sys需要低层的多项式类型。sympy_eqs_to_ring(eqs, symbols)正是完成这个转换的桥梁参数类型含义eqslist[Expr]以Expr实例表示的方程组symbolslist[Symbol]方程组中的未知量返回值为(list[PolyElement], Ring)方程作为PolyElement实例的列表以及这些方程所在的多项式环。官方示例 from sympy import symbols from sympy.polys.solvers import sympy_eqs_to_ring a, x, y symbols(a, x, y) eqs [x-y, xa*y] eqs_ring, ring sympy_eqs_to_ring(eqs, [x, y]) eqs_ring [x - y, x a*y] type(eqs_ring[0]) class sympy.polys.rings.PolyElement ring ZZ(a)[x,y]转换完成后结果可以直接交给solve_lin_sys from sympy.polys.solvers import solve_lin_sys solve_lin_sys(eqs_ring, ring) {y: 0, x: 0}注意上例中ring的输出是ZZ(a)[x,y]——表示环的域是ZZ(a)以a为参数的有理函数域生成元是x, y。这体现了多项式环求解器的通用性系数本身可以含参数。实现细节源码实现sympy/polys/solvers.py#L78-L130核心只有两行K, eqs_K sring(eqs, symbols, fieldTrue, extensionTrue) return eqs_K, K.to_domain()sring见 sympy/polys/rings.py会根据方程构造一个多项式环fieldTrue表示把系数域构造成域保证可除性extensionTrue表示自动把表达式中出现的非符号量如sqrt(2)、log(3)等代数数扩张进系数域。这里还有一个关键的容错分支如果sring抛出NotInvertible对应 GitHub issue sympy#18874 中的场景例如系数域构造失败则降级用domainEXEX 域即把整个表达式当作系数重新构造从而保证在特殊输入下也能给出结果而不是崩溃。三、eqs_to_matrix把方程写成增广矩阵eqs_to_matrix(eqs_coeffs, eqs_rhs, gens, domain)是一个内部辅助函数由solve_lin_sys在处理完方程后调用负责把 dict 形式的低层方程组装成增广矩阵参数类型含义eqs_coeffslist[dict[Symbol, DomainElement]]方程的左端从符号到系数的映射系数是DomainElement实例eqs_rhslist[DomainElement]方程的右端同样是DomainElement实例genslist[Symbol]方程组中的未知量domainDomain左右两端系数的公共定义域返回值为系统的增广矩阵表示类型是DomainMatrix。官方示例 from sympy import symbols, ZZ from sympy.polys.solvers import eqs_to_matrix x, y symbols(x, y) eqs_coeff [{x:ZZ(1), y:ZZ(1)}, {x:ZZ(1), y:ZZ(-1)}] eqs_rhs [ZZ(0), ZZ(-1)] eqs_to_matrix(eqs_coeff, eqs_rhs, [x, y], ZZ) DomainMatrix([[1, 1, 0], [1, -1, 1]], (2, 3), ZZ)实现要点源码实现sympy/polys/solvers.py#L20-L75的逻辑非常直白建立sym2index {x: n for n, x in enumerate(gens)}的符号→列号映射矩阵行数等于方程数nrows列数为len(gens) 1多出的最后一列是增广列对每个方程把eq_coeff中每个符号的系数填入对应列并把右端项取负放入最后一列row[-1] -domain.convert(eq_rhs)。把右端项取负这一步是标准的“将A·x b写成增广矩阵[A | -b]”的处理方式——因为DomainMatrix的rref()行最简形算法是按齐次形式A·x 0工作的把b挪到左边并取负后整个系统就化为了齐次形式。在测试文件 sympy/polys/tests/test_solvers.py#L108-L113 中验证了它与普通Matrix的等价性def test_eqs_to_matrix(): domain, x1,x2 ring(x1,x2, QQ) eqs_coeff [{x1: QQ(1), x2: QQ(1)}, {x1: QQ(2), x2: QQ(-1)}] eqs_rhs [QQ(-5), QQ(0)] M eqs_to_matrix(eqs_coeff, eqs_rhs, [x1, x2], QQ) assert M.to_Matrix() Matrix([[1, 1, 5], [2, -1, 0]])四、_solve_lin_sys连通分量拆分_solve_lin_sys(eqs_coeffs, eqs_rhs, ring)是内部函数由solve_lin_sys在预处理完成后调用。它的核心职责是把整个大方程组按未知量的关联关系拆分成若干个相互独立的“连通分量”connected components再分别求解。拆分思路从源码sympy/polys/solvers.py#L250-L308可以看出其图论思想把每个未知量当作图的一个顶点构造顶点集V ring.gens对每个方程取其左端出现的所有符号并用pairwise(syms)在相邻符号之间连边得到边集E用sympy.utilities.iterables.connected_components求出连通分量建立sym2comp映射把每个方程归入其未知量所在的分量形成若干个子系统依次对每个子系统调用_solve_lin_sys_component求解只要有一个分量无解整个系统就无解返回None否则合并所有分量的解返回。这种拆分对性能意义重大线性方程组的求解复杂度与未知量个数大致呈立方关系把一个大系统拆成若干互不相干的小系统后总求解成本会大幅下降而且天然适合并行化/模块化处理。官方示例 from sympy import symbols, sring from sympy.polys.solvers import _solve_lin_sys x, y symbols(x, y) R, (xr, yr) sring([x, y], [x, y]) eqs [{xr:R.one, yr:-R.one}, {xr:R.one, yr:R.one}] eqs_rhs [R.zero, -2*R.one] _solve_lin_sys(eqs, eqs_rhs, R) {y: 1, x: 1}这里对应方程组x - y 0与x y 2解为x 1, y 1。五、_solve_lin_sys_component高斯-若尔当消元与回代_solve_lin_sys_component(eqs_coeffs, eqs_rhs, ring)同样是内部函数是求解的“最后一公里”。官方文档对其算法做了明确说明The system of equations is solved using Gauss-Jordan elimination with division followed by back-substitution.使用带除法的 Gauss-Jordan 消元随后进行回代求解。求解步骤源码实现sympy/polys/solvers.py#L311-L381的流程是矩阵化调用eqs_to_matrix把该分量转成增广矩阵扩域如果矩阵的定义域不是域Field则通过matrix.to_field()转换到域上如从整数环ZZ转到有理数域QQ保证消元时可做除法行化简调用matrix.rref()得到行最简阶梯形reduced row echelon form以及主元列索引pivots判定无解如果主元的最后一个恰好落在增广列pivots[-1] len(keys)说明出现了0 常数的矛盾行返回None唯一解情形当主元个数等于未知量个数len(pivots) len(keys)直接取每行增广列的值作为对应未知量的解无穷多解情形当主元少于未知量个数非主元未知量作为自由参数主元未知量用自由参数线性表出。实现上先把阶梯形中的地面域系数转换回环元素ring.ring.ground_new或ring.ground_new再对每个主元列p计算v echelon[i][-1] - sum(echelon[i][j]*g[j] for j in range(p1, len(g)) if echelon[i][j])即“右端值减去所有自由项”从而得到形如x1 3 - x3的解析解。官方示例 from sympy import symbols, sring from sympy.polys.solvers import _solve_lin_sys_component x, y symbols(x, y) R, (xr, yr) sring([x, y], [x, y]) eqs [{xr:R.one, yr:-R.one}, {xr:R.one, yr:R.one}] eqs_rhs [R.zero, -2*R.one] _solve_lin_sys_component(eqs, eqs_rhs, R) {y: 1, x: 1}六、三种求解结果形态唯一解 / 无解 / 无穷多解通过测试文件 sympy/polys/tests/test_solvers.py可以清晰看到solve_lin_sys对三种结果的完整覆盖唯一解2 方程 2 未知量def test_solve_lin_sys_2x2_one(): domain, x1,x2 ring(x1,x2, QQ) eqs [x1 x2 - 5, 2*x1 - x2] sol {x1: QQ(5, 3), x2: QQ(10, 3)} _sol solve_lin_sys(eqs, domain) assert _sol sol and all(s.ring domain for s in _sol)注意其中的断言all(s.ring domain for s in _sol)返回的每个PolyElement解仍属于原来的环说明求解过程不会“污染”结果所在的代数结构。无解矛盾方程组def test_solve_lin_sys_2x4_none(): domain, x1,x2 ring(x1,x2, QQ) eqs [x1 - 1, x1 - x2, x1 - 2*x2, x2 - 1] assert solve_lin_sys(eqs, domain) is None无穷多解欠定方程组自由变量保留在表达式中def test_solve_lin_sys_3x3_inf(): domain, x1,x2,x3 ring(x1,x2,x3, QQ) eqs [x1 - x2 2*x3 - 1, 2*x1 x2 x3 - 8, x1 x2 - 5] sol {x1: -x3 3, x2: x3 2} assert solve_lin_sys(eqs, domain) sol此外测试还覆盖了方程数多于未知量的超定系统3x4、大稀疏系统4x7、5x5以及一个系数来自参数域field(d,r,e,..., ZZ)的 6×6 电路方程组test_solve_lin_sys_6x6_1/test_solve_lin_sys_6x6_2——后者充分展示了该求解器处理系数带参数的线性方程组的能力解中会自然出现以参数表示的分式表达式。七、模块在 SymPy 中的实际调用链sympy.polys.solvers不是孤立存在的它是 SymPy 多个核心功能的公共底座。通过代码检索可以确认以下真实调用关系1.linsolve—— 直接包装linsolvesympy.solvers.solveset.linsolve是面向用户求解线性方程组的高层入口。其源码sympy/solvers/solveset.py#L3149-L3169中有一段关键注释# This is just a wrapper for solve_lin_sys eqs [] rows A.tolist() ... eqs, ring sympy_eqs_to_ring(eqs, symbols) sol solve_lin_sys(eqs, ring, _rawFalse) if sol is None: return S.EmptySet sol FiniteSet(Tuple(*(sol.get(sym, sym) for sym in symbols)))也就是说linsolve会把矩阵形式的输入展开成方程列表经sympy_eqs_to_ring转换到多项式环然后调用solve_lin_sys(..., _rawFalse)求解无解时返回空集S.EmptySet。因此理解本模块就理解了linsolve的底层行为。2.solve通用求解器solve在解决策中会提取多项式部分并调用sympy_eqs_to_ring与solve_lin_sys见 sympy/solvers/solveset.py#L54 的导入以及后续_handle_poly分支。当用户用solve求解的方程恰好是线性组时最终会落到本模块。3. Risch 积分算法在符号积分领域heurisch启发式积分Risch 算法族之一在求解待定系数时直接调用了solve_lin_syssympy/integrals/heurisch.py#L39 导入L743 调用from sympy.polys.solvers import solve_lin_sys ... solution solve_lin_sys(numer.coeffs(), coeff_ring, _rawFalse)它把被积函数分母展开后的分子多项式系数提取出来构成一个以“积分系数”为未知量的线性方程组再交给solve_lin_sys求解从而确定原函数中各个待定系数的值。4. 性能基准sympy/polys/benchmarks/bench_solvers.py 中内置了若干针对solve_lin_sys的基准用例如R_165、R_49、R_8等不同规模的多项式环用于持续监控求解器的性能防止回归。八、与新一代稀疏求解器_linsolve的对比值得注意的是SymPy 中还存在一个更新的求解路径sympy.polys.matrices.linsolve._linsolve。该模块的头部注释sympy/polys/matrices/linsolve.py#L1-L27明确说明了二者的关系This is a replacement for solve_lin_sys in sympy.polys.solvers which is inefficient for large sparse systems due to the use of a PolyRing with many generators.也就是说solve_lin_sys基于PolyRing生成元数量等于未知量个数内部使用稠密矩阵实现DDM对于小型稠密系统更高效_linsolve基于SDM稀疏矩阵实现专门针对大型稀疏线性系统优化通过SDM.rref、SDM.nullspace等方法在系数域内直接完成求解。二者共享PolyNonlinearError异常类。官方文档中的这一定位也解释了为何solve_lin_sys虽为公开函数却始终强调“主要供内部使用”——它的接口面向多项式环普通用户应优先选择语义更友好的linsolve。九、实践建议与小结什么时候该用solve_lin_sys编写 SymPy 内部扩展/插件需要在一个已构造好的PolynomialRing上求解线性方程组时可以直接使用它省去反复在Expr与环元素之间转换的开销追求符号系数的精确性当系数来自参数域如ZZ(a)[x,y]时它能给出带参数的精确分式解需要None语义的无解判定None返回值可以直接映射到上层逻辑的“系统矛盾”分支。什么时候不该用普通用户求解线性方程组请使用sympy.solvers.solveset.linsolve处理大型稀疏系统可以关注sympy.polys.matrices.linsolve._linsolve所在的新路径。核心要点回顾sympy.polys.solvers是 SymPy 内部的低层线性方程组求解器核心优势在于把运算下沉到DomainElement/DomainMatrix级别完整调用链为Expr方程组 →sympy_eqs_to_ring转成PolyRing元素 →solve_lin_sys预处理非线性检查、常数项分离→_solve_lin_sys按图连通性拆分 →_solve_lin_sys_component用eqs_to_matrix构造增广矩阵并做rref消元 → 唯一解/无解/无穷多解三种结果官方文档要求环的域必须是域Field否则assert ring.domain.is_Field会失败遇到NotInvertible时实现会自动降级到EX域该模块为linsolve、solve和 Risch 启发式积分提供了统一、精确、可复用的线性求解能力是 SymPy 代数内核中承上启下的关键一环。相关参考文件模块官方文档核心源码实现单元测试高层入口 linsolve包装 solve_lin_sys稀疏求解替代实现 _linsolveRisch 启发式积分中的调用示例性能基准【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考
返回列表