Y a-t-il quelque chose d'inexact dans cette implémentation d'une rotation 2-D en numpy?

Question

Je me suis gratté la tête à propos de la précision d'une rotation de 2 jours en numpy. La mise en œuvre est manuel, et mon application appelle un système gaucher, donc:

from numpy import sin, cos 

def rotate(pathx, pathy, r): 
    """ 
    pathx and pathy are lists of np.float64's recording x and y 
    coordinates of points to be rotated 
    """ 
    c = cos(r) 
    s = sin(r) 
    pathx = c*pathx + s*pathy 
    pathy = -s*pathx + c*pathy 

Pour tester cela, je nourris dans pathx=[1] et pathy=[1] et r=arctan2(1,1)~=pi/4~=0.78539816339744828 avec le résultat attendu que le vecteur (1,1) serait être aligné avec l'axe des x après avoir été tourné par pi/4

Je reçois [1.4142135623730949, -0.29289321881345221]. La raison pour laquelle cela semble absurde est parce que je m'attendais à beaucoup plus plus proche de 0.0 sur la coordonnée y. J'ai également essayé de booster les choses en faisant à la fois les entrées et les sorties sin, cos et arctan2 toutes dtype=float64, mais cela ne faisait aucune différence.

Est-ce que je fais une erreur stupide? Ou y a-t-il une instabilité numérique que j'aurais dû anticiper? Je ne peux pas croire que l'ampleur de l'erreur ...

Answer 1

Lorsque vous attribuez

pathy = -s*pathx + c*pathy 

ce utilise maintenant la variable pathx que vous juste affecté, ce qui a déjà été tourné. Vous souhaitez utiliser le tableau qui a été initialement transmis à la fonction en tant que pathx, avant la rotation.

Corrélativement, une solution facile de maintenir votre approche actuelle serait

def rotate(pathx, pathy, r): 
    c = np.cos(r) 
    s = np.sin(r) 
    pathx_new = c * pathx + s * pathy 
    pathy_new =-s * pathx + c * pathy 
    return pathx_new, pathy_new 

qui, pour votre exemple le cas retourne

(array([ 1.41421356]), array([ 1.11022302e-16])) 

Answer 2

Peut-être que cela va donner un sens de l'erreur:

pathx_new = c*pathx  + s*pathy 
pathy  = -s*pathx_new + c*pathy 

essentiellement, vous mettez à jour pathx et la définition pathy en fonction de sa nouvelle valeur

Vous voulez faire cela, car en place réattribution (plus difficile que prévu en fait):

def rotate(pathxy, r): 
    """ 
    pathx and pathy are lists of np.float64's recording x and y 
    coordinates of points to be rotated 
    """ 
    c = np.cos(r) 
    s = np.sin(r) 
    rot = np.array([[c, s], [-s, c]]).T 
    pathxy.dot(rot, out = pathxy) 

pathxy = np.array([pathx, pathy]).T.astype(float) 

puis pour pathx = pathy = 1,

rotate(pathxy, 0.78539816339744828) 
array([ 1.41421356, 0.  ]) 

et pathxy = np.ones((3,2))

rotate(pathxy, 0.78539816339744828) 

pathxy 

array([[ 1.41421356, 0.  ], 
     [ 1.41421356, 0.  ], 
     [ 1.41421356, 0.  ]])