我试图从时间相关的旋转矩阵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
所以这是一个错误.让开发人员修复它.