python - Incorrect results with Sympy when utilizing (numpy's) floats -

i trying calculate velocity tensor time dependent rotationmatrix re(t) (namely earth rotation @ latitude 48.3°). achieved determining skew symmetric matrix se(t) = dre(t)/dt * re.t. obtaining incorrect results when utilizing float instead of sympy expression, shown in following example:

from ipython.display import display import sympy sy  sy.init_printing()  # latex pretty printing 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 +                          "' not in ['x', 'y', 'z']!")   t, lat = sy.symbols("t, lat", real=true)  # time , 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 float print("\nlat_sy - lat_fl = {}".format((lat_sy - lat_fl).evalf()))  # earth rotation matrix @ latitiude 48.232°: re = (mk_rotmatrix(ome*t, "z") * mk_rotmatrix(lat - sy.pi/2, "y")) # substitute latitude sympy , float value: re_sy, re_fl = re.subs(lat, lat_sy), re.subs(lat, lat_fl)  # angular velocity in world coordinates 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 sympy latitude ({}):".format(lat_sy)) display(se_sy)  # correct result print("\nangular velocity float latitude ({}):".format(lat_fl)) display(se_fl)  # incorrect result 

the result is:

for float latitude result totally wrong in spite of difference of -3e-17 sympy value. not clear me, why happens. numerically, calculation not seem problematic.

my question is, how work around such deficits. should avoid mixing sympy , float/numpy data types? quite difficult detect more complex settings.

ps: sympy version 0.7.6.

tl; dr

it bug. if don't believe it, try this:

in [1]: sympy import factor, symbol  in [2]: factor(1e-20*symbol('t')-7.292115e-5) out[2]: -2785579325.00000 

---

two years ago, default value parameter tol in realfield.__init__ changed none false in commit polys: disabled automatic reduction 0 in rr , cc.
later, tol reverted none fix simplification issue, in commit changed tol on complex , real field none.
seems developers didn't expect reversion bring other issue.

if modify tol=none @ realfield.__init__ in realfield.py, tol=false, correct result 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]]) 

---

the change of tol can explain why you've got wrong result, don't thint root of issue.
imho, there deficiency in polynomial factorization in sympy. i'll illustrate deficiency.
convenience, let preparation work.
add followings example.

from sympy import simplify, expand, s sympy.polys import factor sympy.polys.domains import qq, rr, realfield sympy.polys.factortools import dup_convert sympy.polys.polytools import poly sympy.polys.polytools import _symbolic_factor_list, _poly_from_expr sympy.polys.polyerrors import polificationfailed sympy.polys import polyoptions options 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 exc:         print(exc)         print(exc.expr)     else:         _coeff, _factors = poly.factor_list()         print(poly)         print(_coeff, _factors)         return poly 

we don't need study whole matrices. let focus on 1 element, element @ row 1 , column 2. has shown result incorrect.

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 

when call simplify, in case, it's equivalent combination of 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 

now, know wrong results produced during invocation of factor().
actually, factor wrapper, major work done _symbolic_factor_list.

in [20]: _symbolic_factor_list(expr_fl, opt, 'factor') out[20]: (-2785579325.00000, []) 

let take @ _symbolic_factor_list. key part is:

        try:             poly, _ = _poly_from_expr(base, opt)         except polificationfailed exc:             factors.append((exc.expr, exp))         else:             func = getattr(poly, method + '_list')              _coeff, _factors = func() 

we use above my_symbolic_factor_list simulate procedure.

in [22]: expand(expr_sy) out[22]: -7.29211500000000e-5  in [23]: my_symbolic_factor_list(expr_sy) can't construct polynomial -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 polynomial 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, []) 

by design, constant polynomial should execute except polificationfailed exc: suite, while other polynomials should execute else: suite.
expr_sy, number after expand(), , 1 both constant polynomials, polificationfaileds thrown.
poly_fl -7.292115e-5 * sin(7.292115e-5*t) ** 0, namely, -7.292115e-5, constant polynomial, whereas expr_fl not. supposed same polynomial, different representation. not.
deficiency mentioned.

where missing 1.35525271560688e-20*sin(7.292115e-5*t)**2?
let recall: tol reverted none, means automatic reduction 0 in rr enabled again.
1.35525271560688e-20 reduced zero. thus, poly_fl became constant polynomial.
if tol false, won't happen.

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 

now, can explain why you've got -2785579325.0. in else: suite, poly.factor_list called.
according docs:

factor_list(f)[source]

returns list of irreducible factors of f.

poly_fl supposed non constant polynomial, number. thus, sympy tring use rational number approximate poly_fl. numerator kept, while denominator discarded.

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, []) 

i don't think should blame mixing sympy , float/numpy data types. problem not caused pitfalls sympy mentioned.
simple factorization can produce counterintuitive result.

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 

so bug. let developers fix it.