python – 使用(numpy)浮点数时Sympy的结果不正确

我试图从时间相关的旋转矩阵RE(t)(即纬度48.3°的地球自转)计算velocity tensor.这是通过确定偏斜对称矩阵SE(t)= dRE(t)/ dt * RE.T来实现的.使用float而不是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

结果是:

对于浮动纬度,尽管只有-3e-17与Sympy值的差异,但结果完全错误.我不清楚为什么会这样.从数字上看,这种计算似乎没有问题.

我的问题是,如何解决这些赤字问题.我应该避免混合Sympy和float / Numpy数据类型吗?对于更复杂的设置,它们很难检测到.

PS:Sympy版本是0.7.6.

解决方法:

TL; DR

这是一个错误.如果你不相信,试试这个:

In [1]: from sympy import factor, Symbol

In [2]: factor(1e-20*Symbol('t')-7.292115e-5)
Out[2]: -2785579325.00000

两年前,在提交polys: Disabled automatic reduction to zero in RR and CC中,RealField .__ init__中参数tol的默认值从None更改为False.
后来,在提交Changed tol on Complex and Real field to None中,tol被恢复为无以修复简化问题.
看来开发人员没想到这种回归会带来一些其他问题.

如果在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

当我们称之为简化时,在这种情况下,它几乎等同于TR6和因子的组合.

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()期间会产生错误的结果.
实际上,因素只是一个包装,主要工作是在_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, [])

按照设计,常量多项式应该执行除PolificationFailed外的exc:suite,而其他多项式应该执行else:suite.
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

所以这是一个错误.让开发人员修复它.

上一篇:python – sympy – 象征性地用浮动力求解方程


下一篇:python – 使用符号反转矩阵3×3