Sympy中linsolve与solve求解线性系统得到不同结果的问题排查

核心原因分析 #

• 浮点数系数的有理化转换误差：SymPy的linsolve默认会将输入的浮点数系数自动转换为有理数进行精确符号运算，若你的浮点数是无理数的近似值（比如π、sqrt(2)的截断值），转换后的有理数会和实际系数存在微小偏差，52阶线性系统的求解过程会放大这种偏差，最终得到错误的解。而solve在你设置simplify=False时不会做强制有理化转换，直接基于输入的浮点系数求解，因此结果正确。
• 高条件数带来的数值稳定性问题：即使矩阵满秩，若52阶系数矩阵的条件数过高，linsolve默认的求解路径（比如LU分解）会放大浮点误差，而solve的通用求解逻辑对线性系统会自动选择更稳定的分解策略，因此结果更准确。
• 方程解析逻辑差异：linsolve直接解析方程列表时的标准化逻辑和solve存在差异，可能会错误处理某些带特殊写法的浮点系数方程，导致转换后的矩阵和原方程不一致。
• 旧版本SymPy的已知缺陷：部分旧版本的SymPy（<1.10）的linsolve在处理大于40阶的浮点系数线性系统时存在已知的求解bug，无报错但输出错误结果。

修复方案 #

• 直接传入矩阵而非方程列表求解，避免方程解析误差：

from sympy import linear_eq_to_matrix, linsolve
A, b = linear_eq_to_matrix(concrete_focs, lagran)
linresult = list(linsolve((A, b), lagran))

• 关闭linsolve的默认有理化转换，直接用浮点运算求解：

linresult = list(linsolve(concrete_focs, lagran, rational=False))

• 若矩阵条件数过高，使用更高精度的浮点运算求解：

A, b = linear_eq_to_matrix(concrete_focs, lagran)
# 用30位有效数字运算，可根据需求调整精度
A_high = A.evalf(30)
b_high = b.evalf(30)
linresult = list(linsolve((A_high, b_high), lagran))

• 升级SymPy到最新稳定版本，排除已知bug影响。

内容的提问来源于stack exchange，提问作者econbernardo