使用（numpy 的）浮点数时 Sympy 的结果不正确
Incorrect results with Sympy when utilizing (numpy's) floats
我正在尝试从依赖于时间的旋转矩阵 RE(t)（即纬度 48.3° 处的地球自转）计算 velocity tensor。这是通过确定偏斜对称矩阵 SE(t) = dRE(t)/dt * RE.T 来实现的。当使用浮点数而不是 Sympy 表达式时，我得到了不正确的结果，如下例所示：
from IPython.display import display
import sympy as sy
sy.init_printing() # LaTeX like pretty printing for IPython
def mk_rotmatrix(alpha, coord_ax="x"):
""" Rotation matrix around coordinate axis """
ca, sa = sy.cos(alpha), sy.sin(alpha)
if coord_ax == "x":
return sy.Matrix([[1, 0, 0],
[0, ca, -sa],
[0, sa, +ca]])
elif coord_ax == 'y':
return sy.Matrix([[+ca, 0, sa],
[0, 1, 0],
[-sa, 0, ca]])
elif coord_ax == 'z':
return sy.Matrix([[ca, -sa, 0],
[sa, +ca, 0],
[0, 0, 1]])
else:
raise ValueError("Parameter coord_ax='" + coord_ax +
"' is not in ['x', 'y', 'z']!")
t, lat = sy.symbols("t, lat", real=True) # time and latitude
omE = 7.292115e-5 # rad/s -- earth rotation rate (15.04107 °/h)
lat_sy = 48.232*sy.pi/180 # latitude in rad
lat_fl = float(lat_sy) # latitude as float
print("\nlat_sy - lat_fl = {}".format((lat_sy - lat_fl).evalf()))
# earth rotation matrix at latitiude 48.232°:
RE = (mk_rotmatrix(omE*t, "z") * mk_rotmatrix(lat - sy.pi/2, "y"))
# substitute latitude with sympy and float value:
RE_sy, RE_fl = RE.subs(lat, lat_sy), RE.subs(lat, lat_fl)
# Angular velocity in world coordinates as skew symmetric matrix:
SE_sy = sy.simplify(RE_sy.diff(t) * RE_sy.T)
SE_fl = sy.simplify(RE_fl.diff(t) * RE_fl.T)
print("\nAngular velocity with Sympy latitude ({}):".format(lat_sy))
display(SE_sy) # correct result
print("\nAngular velocity with float latitude ({}):".format(lat_fl))
display(SE_fl) # incorrect result
结果是：
对于浮动纬度，结果是完全错误的，尽管与 Sympy 值仅相差 -3e-17。我不清楚为什么会这样。从数值上看，这个计算似乎没有问题。
我的问题是，如何解决此类缺陷。我应该避免混合使用 Sympy 和 float/Numpy 数据类型吗？对于更复杂的设置，它们很难检测到。
PS: Sympy 版本为 0.7.6.
我认为这可能是 Sympy 中的一个错误；当我 运行 你的脚本在我的系统上时（Ubuntu 14.04 64 位，Python 2.7，Sympy 0.7.4.1），我得到
lat_sy - lat_fl = -2.61291277482447e-17
Angular velocity with Sympy latitude (0.267955555555556*pi):
Matrix([
[ 0, -7.292115e-5, 0],
[7.292115e-5, 0, 0],
[ 0, 0, 0]])
Angular velocity with float latitude (0.841807204822):
Matrix([
[3.3881317890172e-21*sin(0.0001458423*t), -7.29211495242194e-5, 0],
[ 7.29211495242194e-5, -3.3881317890172e-21*sin(0.0001458423*t), 0],
[ 0, 0, 0]])
看起来不错。
我不确定该建议什么：您可以尝试比 0.7.6 更旧的 Sympy 版本，或者 Github 的最新版本。
[回复评论]至于为什么对角线不为零，我的第一个评论是 3e-21/7e-5 大约是 4e-17； IEEE754 64 位 ("float")
数值精度约为 2e-16。在 3e-21 rad/s 一次公转将需要 60 万亿年（约 2e21 秒）。不用担心。
我不完全确定这里发生了什么，但是在将它添加到您的脚本之后
def matrix_product_element(a, b, i, j):
v = a[3*i:3*i+3]
w = b[j::3]
summand_list = [v[k]*w[k]
for k in range(3)]
print('element ({},{})'.format(i, j))
print(' summand_list: {}'.format(summand_list))
print(' sum(summand_list): {}'.format(sum(summand_list)))
print(' sum(summand_list).simplify(): {}'.format(sum(summand_list)))
matrix_product_element(RE_fl.diff(t), RE_fl.T, 0, 0)
matrix_product_element(RE_fl.diff(t), RE_fl.T, 1, 0)
matrix_product_element(RE_fl.diff(t), RE_fl.T, 2, 0)
sumlist=[sy.Float(-4.05652668591092e-5,15), sy.Float(7.292115e-5,15), sy.Float(-3.23558831408908e-5,14)]
display(sumlist)
display(sum(sumlist))
我明白了
element (0,0)
summand_list: [-4.05652668591092e-5*sin(7.292115e-5*t)*cos(7.292115e-5*t), 7.292115e-5*sin(7.292115e-5*t)*cos(7.292115e-5*t), -3.23558831408908e-5*sin(7.292115e-5*t)*cos(7.292115e-5*t)]
sum(summand_list): 6.7762635780344e-21*sin(7.292115e-5*t)*cos(7.292115e-5*t)
sum(summand_list).simplify(): 6.7762635780344e-21*sin(7.292115e-5*t)*cos(7.292115e-5*t)
element (1,0)
summand_list: [4.05652668591092e-5*cos(7.292115e-5*t)**2, 7.292115e-5*sin(7.292115e-5*t)**2, 3.23558831408908e-5*cos(7.292115e-5*t)**2]
sum(summand_list): 7.292115e-5*sin(7.292115e-5*t)**2 + 7.292115e-5*cos(7.292115e-5*t)**2
sum(summand_list).simplify(): 7.292115e-5*sin(7.292115e-5*t)**2 + 7.292115e-5*cos(7.292115e-5*t)**2
element (2,0)
summand_list: [0, 0, 0]
sum(summand_list): 0
sum(summand_list).simplify(): 0
[-4.05652668591092e-5, 7.29211500000000e-5, -3.2355883140891e-5]
6.77626357803440e-21
第一次求和的系数总和应该为零，但事实并非如此。在最后几行中，我设法通过以较低的精度重新创建系数来伪造这种效果（这只是运气，可能并不那么重要）。它是 "sort-of"，因为列表中的第三个值 (-3.2355883140891e-5) 与被加数列表 (-3.23558831408908e-5) 中的系数不匹配，后者被赋予 15 个位置。
Sympy 文档在这里 http://docs.sympy.org/dev/gotchas.html#evaluating-expressions-with-floats-and-rationals 讨论了这类问题，并提供了一些关于如何缓解问题的建议。这是您的代码的一个直接变体，将浮点数的替换推迟到最后：
# encoding:utf-8
from IPython.display import display
import sympy as sy
sy.init_printing() # LaTeX like pretty printing for IPython
def mk_rotmatrix(alpha, coord_ax="x"):
""" Rotation matrix around coordinate axis """
ca, sa = sy.cos(alpha), sy.sin(alpha)
if coord_ax == "x":
return sy.Matrix([[1, 0, 0],
[0, ca, -sa],
[0, sa, +ca]])
elif coord_ax == 'y':
return sy.Matrix([[+ca, 0, sa],
[0, 1, 0],
[-sa, 0, ca]])
elif coord_ax == 'z':
return sy.Matrix([[ca, -sa, 0],
[sa, +ca, 0],
[0, 0, 1]])
else:
raise ValueError("Parameter coord_ax='" + coord_ax +
"' is not in ['x', 'y', 'z']!")
# time [s], latitude [rad], earth rate [rad/s]
t, lat, omE = sy.symbols("t, lat, omE", real=True)
RE = (mk_rotmatrix(omE*t, "z") * mk_rotmatrix(lat - sy.pi/2, "y"))
SE = sy.simplify(RE.diff(t) * RE.T)
display(SE)
display(SE.subs({lat: 48.232*sy.pi/180, omE: 7.292115e-5}))
这给出：
Matrix([
[ 0, -omE, 0],
[omE, 0, 0],
[ 0, 0, 0]])
Matrix([
[ 0, -7.292115e-5, 0],
[7.292115e-5, 0, 0],
[ 0, 0, 0]])
无论数值优势如何，我都更喜欢这个，因为人们可以从符号解的形式中学到一些东西。
TL;博士
这是一个错误。不信你试试这个：
In [1]: from sympy import factor, Symbol
In [2]: factor(1e-20*Symbol('t')-7.292115e-5)
Out[2]: -2785579325.00000
两年前，RealField.__init__ 中参数 tol 的默认值在提交 polys: Disabled automatic reduction to zero in RR and CC 中从 None 更改为 False。
后来，tol 恢复为 None 以修复简化问题，提交 Changed tol on Complex and Real field to None.
开发者似乎没想到这次回归会带来一些其他问题。
如果您将 realfield.py 中的 RealField.__init__ 处的 tol=None 修改为 tol=False，您将得到 SE_fl 的正确结果。
Matrix([
[3.3881317890172e-21*sin(0.0001458423*t), -7.29211495242194e-5, 0],
[ 7.29211495242194e-5, -3.3881317890172e-21*sin(0.0001458423*t), 0],
[ 0, 0, 0]])
tol的变化可以解释为什么你得到了错误的结果，但我不认为这是问题的根源。
恕我直言，SymPy 中的多项式分解存在缺陷。我会说明这个不足。
为了方便起见，让我们做一些准备工作。
将以下内容添加到您的示例中。
from sympy import simplify, expand, S
from sympy.polys import factor
from sympy.polys.domains import QQ, RR, RealField
from sympy.polys.factortools import dup_convert
from sympy.polys.polytools import Poly
from sympy.polys.polytools import _symbolic_factor_list, _poly_from_expr
from sympy.polys.polyerrors import PolificationFailed
from sympy.polys import polyoptions as options
from sympy.simplify.fu import TR6
def new_opt():
args = dict()
options.allowed_flags(args, [])
opt = options.build_options((), args)
return opt
def my_symbolic_factor_list(base):
opt = new_opt()
try:
poly, _ = _poly_from_expr(base, opt)
except PolificationFailed as exc:
print(exc)
print(exc.expr)
else:
_coeff, _factors = poly.factor_list()
print(poly)
print(_coeff, _factors)
return poly
我们不需要研究整个矩阵。让我们关注一个元素，第1行第2列的元素。它已经显示结果不正确。
In [8]: elm_sy = (RE_sy.diff(t) * RE_sy.T)[1]
In [9]: elm_fl = (RE_fl.diff(t) * RE_fl.T)[1]
In [10]: elm_sy
Out[10]: -7.292115e-5*sin(0.267955555555556*pi)**2*sin(7.292115e-5*t)**2 - 7.292115e-5*sin(7.292115e
-5*t)**2*cos(0.267955555555556*pi)**2 - 7.292115e-5*cos(7.292115e-5*t)**2
In [11]: elm_fl
Out[11]: -7.292115e-5*sin(7.292115e-5*t)**2 - 7.292115e-5*cos(7.292115e-5*t)**2
In [12]: simplify(elm_sy)
Out[12]: -7.29211500000000e-5
In [13]: simplify(elm_fl)
Out[13]: -2785579325.00000
当我们调用simplify时，在这种情况下，它几乎相当于TR6和factor的组合。
In [15]: expr_sy = TR6(elm_sy)
In [16]: expr_fl = TR6(elm_fl)
In [17]: expr_fl
Out[17]: 1.35525271560688e-20*sin(7.292115e-5*t)**2 - 7.292115e-5
In [18]: factor(expr_fl)
Out[18]: -2785579325.00000
现在，我们知道调用 factor() 时会产生错误的结果。
实际上，factor只是一个包装器，主要工作由_symbolic_factor_list完成。
In [20]: _symbolic_factor_list(expr_fl, opt, 'factor')
Out[20]: (-2785579325.00000, [])
让我们来看看_symbolic_factor_list。关键部分是：
try:
poly, _ = _poly_from_expr(base, opt)
except PolificationFailed as exc:
factors.append((exc.expr, exp))
else:
func = getattr(poly, method + '_list')
_coeff, _factors = func()
我们用上面的my_symbolic_factor_list来模拟这个程序。
In [22]: expand(expr_sy)
Out[22]: -7.29211500000000e-5
In [23]: my_symbolic_factor_list(expr_sy)
can't construct a polynomial from -7.292115e-5*sin(0.267955555555556*pi)**2*sin(7.292115e-5*t)**2 -
7.292115e-5*(-sin(0.267955555555556*pi)**2 + 1)*sin(7.292115e-5*t)**2 + 7.292115e-5*sin(7.292115e-5*
t)**2 - 7.292115e-5
-7.29211500000000e-5
In [24]: my_symbolic_factor_list(S(1))
can't construct a polynomial from 1
1
In [25]: expr_fl
Out[25]: 1.35525271560688e-20*sin(7.292115e-5*t)**2 - 7.292115e-5
In [26]: poly_fl = my_symbolic_factor_list(expr_fl)
Poly(-7.292115e-5, sin(7.292115e-5*t), domain='RR')
(-2785579325.00000, [])
根据设计，常数多项式应执行 except PolificationFailed as exc: 组，而其他多项式应执行 else: 组。
expr_sy是expand()之后的一个数，1都是常数多项式，因此抛出了PolificationFailed。
poly_fl是-7.292115e-5 * sin(7.292115e-5*t) ** 0，即-7.292115e-5，一个常数多项式，而expr_fl不是。它们应该是相同的多项式，只是表示形式不同。现在他们不是了。
这就是我说的不足
失踪的1.35525271560688e-20*sin(7.292115e-5*t)**2在哪里？
让我们回忆一下：tol 被还原为 None，这意味着 RR 中的自动归零再次启用。
1.35525271560688e-20 减少为零。这样，poly_fl就变成了常数多项式。
如果 tol 是 False，则不会发生这种情况。
In [31]: arg2 = expr_fl.args[1].args[0]
In [32]: arg2
Out[32]: 1.35525271560688e-20
In [33]: RR.from_sympy(arg2)
Out[33]: 0.0
In [34]: R = RealField(tol=False)
In [35]: R.from_sympy(arg2)
Out[35]: 1.35525271560688e-20
现在，我们可以解释为什么您有 -2785579325.0。在 else: 套件中，调用 Poly.factor_list。
根据 docs:
factor_list(f)[source]
Returns a list of irreducible factors of f.
poly_fl 应该是一个非常数多项式，但它只是一个数字。
因此，SymPy 试图使用有理数来近似 poly_fl。保留分子，舍去分母。
In [42]: poly_fl.factor_list()
Out[42]: (-2785579325.00000, [])
In [43]: dup_convert(poly_fl.coeffs(), RR, QQ)
Out[43]: [-2785579325/38199881995827]
In [44]: Poly([S(1.25)], t, domain='RR').factor_list()
Out[44]: (5.00000000000000, [])
In [45]: dup_convert(Poly([S(1.25)], t, domain='RR').coeffs(), RR, QQ)
Out[45]: [5/4]
In [46]: Poly((RE_fl.diff(t) * RE_fl.T)[3].args[0].args[0], t).factor_list()
Out[46]: (1767051195.00000, [])
我认为我们不应该责怪混合使用 Sympy 和 float/Numpy 数据类型。此问题不是由提到的 pitfalls SymPy 引起的。
即使是非常简单的因式分解也会产生违反直觉的结果。
In [47]: factor(1e-20*t-1.2345e-5)
Out[47]: -539023891.000000
In [48]: factor(S(1e-20)*t-S(1.2345e-5))
Out[48]: -539023891.000000
所以这是一个错误。让开发人员修复它。
我正在尝试从依赖于时间的旋转矩阵 RE(t)（即纬度 48.3° 处的地球自转）计算 velocity tensor。这是通过确定偏斜对称矩阵 SE(t) = dRE(t)/dt * RE.T 来实现的。当使用浮点数而不是 Sympy 表达式时，我得到了不正确的结果，如下例所示：
from IPython.display import display
import sympy as sy
sy.init_printing() # LaTeX like pretty printing for IPython
def mk_rotmatrix(alpha, coord_ax="x"):
""" Rotation matrix around coordinate axis """
ca, sa = sy.cos(alpha), sy.sin(alpha)
if coord_ax == "x":
return sy.Matrix([[1, 0, 0],
[0, ca, -sa],
[0, sa, +ca]])
elif coord_ax == 'y':
return sy.Matrix([[+ca, 0, sa],
[0, 1, 0],
[-sa, 0, ca]])
elif coord_ax == 'z':
return sy.Matrix([[ca, -sa, 0],
[sa, +ca, 0],
[0, 0, 1]])
else:
raise ValueError("Parameter coord_ax='" + coord_ax +
"' is not in ['x', 'y', 'z']!")
t, lat = sy.symbols("t, lat", real=True) # time and latitude
omE = 7.292115e-5 # rad/s -- earth rotation rate (15.04107 °/h)
lat_sy = 48.232*sy.pi/180 # latitude in rad
lat_fl = float(lat_sy) # latitude as float
print("\nlat_sy - lat_fl = {}".format((lat_sy - lat_fl).evalf()))
# earth rotation matrix at latitiude 48.232°:
RE = (mk_rotmatrix(omE*t, "z") * mk_rotmatrix(lat - sy.pi/2, "y"))
# substitute latitude with sympy and float value:
RE_sy, RE_fl = RE.subs(lat, lat_sy), RE.subs(lat, lat_fl)
# Angular velocity in world coordinates as skew symmetric matrix:
SE_sy = sy.simplify(RE_sy.diff(t) * RE_sy.T)
SE_fl = sy.simplify(RE_fl.diff(t) * RE_fl.T)
print("\nAngular velocity with Sympy latitude ({}):".format(lat_sy))
display(SE_sy) # correct result
print("\nAngular velocity with float latitude ({}):".format(lat_fl))
display(SE_fl) # incorrect result
结果是：
对于浮动纬度，结果是完全错误的，尽管与 Sympy 值仅相差 -3e-17。我不清楚为什么会这样。从数值上看，这个计算似乎没有问题。
我的问题是，如何解决此类缺陷。我应该避免混合使用 Sympy 和 float/Numpy 数据类型吗？对于更复杂的设置，它们很难检测到。
PS: Sympy 版本为 0.7.6.
我认为这可能是 Sympy 中的一个错误；当我 运行 你的脚本在我的系统上时（Ubuntu 14.04 64 位，Python 2.7，Sympy 0.7.4.1），我得到
lat_sy - lat_fl = -2.61291277482447e-17
Angular velocity with Sympy latitude (0.267955555555556*pi):
Matrix([
[ 0, -7.292115e-5, 0],
[7.292115e-5, 0, 0],
[ 0, 0, 0]])
Angular velocity with float latitude (0.841807204822):
Matrix([
[3.3881317890172e-21*sin(0.0001458423*t), -7.29211495242194e-5, 0],
[ 7.29211495242194e-5, -3.3881317890172e-21*sin(0.0001458423*t), 0],
[ 0, 0, 0]])
看起来不错。
我不确定该建议什么：您可以尝试比 0.7.6 更旧的 Sympy 版本，或者 Github 的最新版本。
[回复评论]至于为什么对角线不为零，我的第一个评论是 3e-21/7e-5 大约是 4e-17； IEEE754 64 位 ("float") 数值精度约为 2e-16。在 3e-21 rad/s 一次公转将需要 60 万亿年（约 2e21 秒）。不用担心。
我不完全确定这里发生了什么，但是在将它添加到您的脚本之后
def matrix_product_element(a, b, i, j):
v = a[3*i:3*i+3]
w = b[j::3]
summand_list = [v[k]*w[k]
for k in range(3)]
print('element ({},{})'.format(i, j))
print(' summand_list: {}'.format(summand_list))
print(' sum(summand_list): {}'.format(sum(summand_list)))
print(' sum(summand_list).simplify(): {}'.format(sum(summand_list)))
matrix_product_element(RE_fl.diff(t), RE_fl.T, 0, 0)
matrix_product_element(RE_fl.diff(t), RE_fl.T, 1, 0)
matrix_product_element(RE_fl.diff(t), RE_fl.T, 2, 0)
sumlist=[sy.Float(-4.05652668591092e-5,15), sy.Float(7.292115e-5,15), sy.Float(-3.23558831408908e-5,14)]
display(sumlist)
display(sum(sumlist))
我明白了
element (0,0)
summand_list: [-4.05652668591092e-5*sin(7.292115e-5*t)*cos(7.292115e-5*t), 7.292115e-5*sin(7.292115e-5*t)*cos(7.292115e-5*t), -3.23558831408908e-5*sin(7.292115e-5*t)*cos(7.292115e-5*t)]
sum(summand_list): 6.7762635780344e-21*sin(7.292115e-5*t)*cos(7.292115e-5*t)
sum(summand_list).simplify(): 6.7762635780344e-21*sin(7.292115e-5*t)*cos(7.292115e-5*t)
element (1,0)
summand_list: [4.05652668591092e-5*cos(7.292115e-5*t)**2, 7.292115e-5*sin(7.292115e-5*t)**2, 3.23558831408908e-5*cos(7.292115e-5*t)**2]
sum(summand_list): 7.292115e-5*sin(7.292115e-5*t)**2 + 7.292115e-5*cos(7.292115e-5*t)**2
sum(summand_list).simplify(): 7.292115e-5*sin(7.292115e-5*t)**2 + 7.292115e-5*cos(7.292115e-5*t)**2
element (2,0)
summand_list: [0, 0, 0]
sum(summand_list): 0
sum(summand_list).simplify(): 0
[-4.05652668591092e-5, 7.29211500000000e-5, -3.2355883140891e-5]
6.77626357803440e-21
第一次求和的系数总和应该为零，但事实并非如此。在最后几行中，我设法通过以较低的精度重新创建系数来伪造这种效果（这只是运气，可能并不那么重要）。它是 "sort-of"，因为列表中的第三个值 (-3.2355883140891e-5) 与被加数列表 (-3.23558831408908e-5) 中的系数不匹配，后者被赋予 15 个位置。
Sympy 文档在这里 http://docs.sympy.org/dev/gotchas.html#evaluating-expressions-with-floats-and-rationals 讨论了这类问题，并提供了一些关于如何缓解问题的建议。这是您的代码的一个直接变体，将浮点数的替换推迟到最后：
# encoding:utf-8
from IPython.display import display
import sympy as sy
sy.init_printing() # LaTeX like pretty printing for IPython
def mk_rotmatrix(alpha, coord_ax="x"):
""" Rotation matrix around coordinate axis """
ca, sa = sy.cos(alpha), sy.sin(alpha)
if coord_ax == "x":
return sy.Matrix([[1, 0, 0],
[0, ca, -sa],
[0, sa, +ca]])
elif coord_ax == 'y':
return sy.Matrix([[+ca, 0, sa],
[0, 1, 0],
[-sa, 0, ca]])
elif coord_ax == 'z':
return sy.Matrix([[ca, -sa, 0],
[sa, +ca, 0],
[0, 0, 1]])
else:
raise ValueError("Parameter coord_ax='" + coord_ax +
"' is not in ['x', 'y', 'z']!")
# time [s], latitude [rad], earth rate [rad/s]
t, lat, omE = sy.symbols("t, lat, omE", real=True)
RE = (mk_rotmatrix(omE*t, "z") * mk_rotmatrix(lat - sy.pi/2, "y"))
SE = sy.simplify(RE.diff(t) * RE.T)
display(SE)
display(SE.subs({lat: 48.232*sy.pi/180, omE: 7.292115e-5}))
这给出：
Matrix([
[ 0, -omE, 0],
[omE, 0, 0],
[ 0, 0, 0]])
Matrix([
[ 0, -7.292115e-5, 0],
[7.292115e-5, 0, 0],
[ 0, 0, 0]])
无论数值优势如何，我都更喜欢这个，因为人们可以从符号解的形式中学到一些东西。
TL;博士
这是一个错误。不信你试试这个：
In [1]: from sympy import factor, Symbol
In [2]: factor(1e-20*Symbol('t')-7.292115e-5)
Out[2]: -2785579325.00000
两年前，RealField.__init__ 中参数 tol 的默认值在提交 polys: Disabled automatic reduction to zero in RR and CC 中从 None 更改为 False。
后来，tol 恢复为 None 以修复简化问题，提交 Changed tol on Complex and Real field to None.
开发者似乎没想到这次回归会带来一些其他问题。
如果您将 realfield.py 中的 RealField.__init__ 处的 tol=None 修改为 tol=False，您将得到 SE_fl 的正确结果。
Matrix([
[3.3881317890172e-21*sin(0.0001458423*t), -7.29211495242194e-5, 0],
[ 7.29211495242194e-5, -3.3881317890172e-21*sin(0.0001458423*t), 0],
[ 0, 0, 0]])
tol的变化可以解释为什么你得到了错误的结果，但我不认为这是问题的根源。
恕我直言，SymPy 中的多项式分解存在缺陷。我会说明这个不足。
为了方便起见，让我们做一些准备工作。
将以下内容添加到您的示例中。
from sympy import simplify, expand, S
from sympy.polys import factor
from sympy.polys.domains import QQ, RR, RealField
from sympy.polys.factortools import dup_convert
from sympy.polys.polytools import Poly
from sympy.polys.polytools import _symbolic_factor_list, _poly_from_expr
from sympy.polys.polyerrors import PolificationFailed
from sympy.polys import polyoptions as options
from sympy.simplify.fu import TR6
def new_opt():
args = dict()
options.allowed_flags(args, [])
opt = options.build_options((), args)
return opt
def my_symbolic_factor_list(base):
opt = new_opt()
try:
poly, _ = _poly_from_expr(base, opt)
except PolificationFailed as exc:
print(exc)
print(exc.expr)
else:
_coeff, _factors = poly.factor_list()
print(poly)
print(_coeff, _factors)
return poly
我们不需要研究整个矩阵。让我们关注一个元素，第1行第2列的元素。它已经显示结果不正确。
In [8]: elm_sy = (RE_sy.diff(t) * RE_sy.T)[1]
In [9]: elm_fl = (RE_fl.diff(t) * RE_fl.T)[1]
In [10]: elm_sy
Out[10]: -7.292115e-5*sin(0.267955555555556*pi)**2*sin(7.292115e-5*t)**2 - 7.292115e-5*sin(7.292115e
-5*t)**2*cos(0.267955555555556*pi)**2 - 7.292115e-5*cos(7.292115e-5*t)**2
In [11]: elm_fl
Out[11]: -7.292115e-5*sin(7.292115e-5*t)**2 - 7.292115e-5*cos(7.292115e-5*t)**2
In [12]: simplify(elm_sy)
Out[12]: -7.29211500000000e-5
In [13]: simplify(elm_fl)
Out[13]: -2785579325.00000
当我们调用simplify时，在这种情况下，它几乎相当于TR6和factor的组合。
In [15]: expr_sy = TR6(elm_sy)
In [16]: expr_fl = TR6(elm_fl)
In [17]: expr_fl
Out[17]: 1.35525271560688e-20*sin(7.292115e-5*t)**2 - 7.292115e-5
In [18]: factor(expr_fl)
Out[18]: -2785579325.00000
现在，我们知道调用 factor() 时会产生错误的结果。
实际上，factor只是一个包装器，主要工作由_symbolic_factor_list完成。
In [20]: _symbolic_factor_list(expr_fl, opt, 'factor')
Out[20]: (-2785579325.00000, [])
让我们来看看_symbolic_factor_list。关键部分是：
try:
poly, _ = _poly_from_expr(base, opt)
except PolificationFailed as exc:
factors.append((exc.expr, exp))
else:
func = getattr(poly, method + '_list')
_coeff, _factors = func()
我们用上面的my_symbolic_factor_list来模拟这个程序。
In [22]: expand(expr_sy)
Out[22]: -7.29211500000000e-5
In [23]: my_symbolic_factor_list(expr_sy)
can't construct a polynomial from -7.292115e-5*sin(0.267955555555556*pi)**2*sin(7.292115e-5*t)**2 -
7.292115e-5*(-sin(0.267955555555556*pi)**2 + 1)*sin(7.292115e-5*t)**2 + 7.292115e-5*sin(7.292115e-5*
t)**2 - 7.292115e-5
-7.29211500000000e-5
In [24]: my_symbolic_factor_list(S(1))
can't construct a polynomial from 1
1
In [25]: expr_fl
Out[25]: 1.35525271560688e-20*sin(7.292115e-5*t)**2 - 7.292115e-5
In [26]: poly_fl = my_symbolic_factor_list(expr_fl)
Poly(-7.292115e-5, sin(7.292115e-5*t), domain='RR')
(-2785579325.00000, [])
根据设计，常数多项式应执行 except PolificationFailed as exc: 组，而其他多项式应执行 else: 组。
expr_sy是expand()之后的一个数，1都是常数多项式，因此抛出了PolificationFailed。
poly_fl是-7.292115e-5 * sin(7.292115e-5*t) ** 0，即-7.292115e-5，一个常数多项式，而expr_fl不是。它们应该是相同的多项式，只是表示形式不同。现在他们不是了。
这就是我说的不足
失踪的1.35525271560688e-20*sin(7.292115e-5*t)**2在哪里？
让我们回忆一下：tol 被还原为 None，这意味着 RR 中的自动归零再次启用。
1.35525271560688e-20 减少为零。这样，poly_fl就变成了常数多项式。
如果 tol 是 False，则不会发生这种情况。
In [31]: arg2 = expr_fl.args[1].args[0]
In [32]: arg2
Out[32]: 1.35525271560688e-20
In [33]: RR.from_sympy(arg2)
Out[33]: 0.0
In [34]: R = RealField(tol=False)
In [35]: R.from_sympy(arg2)
Out[35]: 1.35525271560688e-20
现在，我们可以解释为什么您有 -2785579325.0。在 else: 套件中，调用 Poly.factor_list。
根据 docs:
factor_list(f)[source]
Returns a list of irreducible factors of f.
poly_fl 应该是一个非常数多项式，但它只是一个数字。
因此，SymPy 试图使用有理数来近似 poly_fl。保留分子，舍去分母。
In [42]: poly_fl.factor_list()
Out[42]: (-2785579325.00000, [])
In [43]: dup_convert(poly_fl.coeffs(), RR, QQ)
Out[43]: [-2785579325/38199881995827]
In [44]: Poly([S(1.25)], t, domain='RR').factor_list()
Out[44]: (5.00000000000000, [])
In [45]: dup_convert(Poly([S(1.25)], t, domain='RR').coeffs(), RR, QQ)
Out[45]: [5/4]
In [46]: Poly((RE_fl.diff(t) * RE_fl.T)[3].args[0].args[0], t).factor_list()
Out[46]: (1767051195.00000, [])
我认为我们不应该责怪混合使用 Sympy 和 float/Numpy 数据类型。此问题不是由提到的 pitfalls SymPy 引起的。
即使是非常简单的因式分解也会产生违反直觉的结果。
In [47]: factor(1e-20*t-1.2345e-5)
Out[47]: -539023891.000000
In [48]: factor(S(1e-20)*t-S(1.2345e-5))
Out[48]: -539023891.000000
所以这是一个错误。让开发人员修复它。