float在python中的书写形式错误的是_python – 使用(numpy)浮点数时Sympy的结果不正确...

最新推荐文章于 2023-10-20 12:41:35 发布

贵妃奶黄包

最新推荐文章于 2023-10-20 12:41:35 发布

阅读量339

点赞数

文章标签： float在python中的书写形式错误的是

本文链接：https://blog.csdn.net/weixin_34159899/article/details/113565559

版权

我试图从时间相关的旋转矩阵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, [])

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

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

所以这是一个错误.让开发人员修复它.来源：http://www.icode9.com/content-1-210701.html

贵妃奶黄包

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
float在python中的书写形式错误的是_python – 使用(numpy)浮点数时Sympy的结果不正确...

我试图从时间相关的旋转矩阵RE(t)(即纬度48.3°的地球自转)计算velocity tensor.这是通过确定偏斜对称矩阵SE(t)= dRE(t)/ dt * RE.T来实现的.使用float而不是Sympy表达式时,我得到的结果不正确,如下例所示：from IPython.display import displayimport sympy as sysy.init_printing() ...
复制链接

扫一扫