lundi 29 juin 2015

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


I am trying to calculate a velocity tensor from a time dependent rotationmatrix RE(t) (Namely the earth rotation at latitude 48.3°). This is achieved by determining the skew symmetric matrix SE(t) = dRE(t)/dt * RE.T. I am obtaining incorrect results when utilizing a float instead of a Sympy expression, as shown in the following example:

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

The result is:

calculation results

For the float latitude the result is totally wrong in spite of the difference of only -3e-17 to the Sympy value. It is not clear to me, why this happens. Numerically, this calculation does not seem to be problematic.

My question is, how to work around such deficits. Should I avoid mixing Sympy and float/Numpy data types? They are quite difficult to detect for more complex settings.

PS: The Sympy version is 0.7.6.


Aucun commentaire:

Enregistrer un commentaire