Siéntete en la libertad de compartir nuestra página y códigos en tus redes, necesitamos de tu ayuda para aumentar esta comunidad.
No hay necesidad de usar Eigensystem o Eigenvectors para encontrar el eje de una matriz de rotación. En su lugar, puede leer los componentes del vector del eje directamente de la matriz simétrica sesgada
$$ a equiv R ^ TR $$
En tres dimensiones (que se supone en la pregunta), aplicar esta matriz a un vector equivale a aplicar un producto cruzado con un vector formado por las tres componentes independientes de $ a $:
1, -1, 1Extract[a, 3, 2, 3, 1, 2, 1]
Este método de una línea para encontrar el eje se aplica en la siguiente función. Para obtener el ángulo de rotación, construyo dos vectores ovec, nvec perpendiculares al eje y entre sí, para encontrar el coseno y el seno del ángulo usando el Dot producto (igualmente podría haber usado Projection). Para obtener un primer vector ovec que no es paralelo al eje, permuto las componentes del vector del eje usando el hecho de que
Solve[x, -y, z == y, z, x, x, y, z]
(* ==> x -> 0, y -> 0, z -> 0 *)
lo que significa que la permutación anterior con cambio de signo de un vector de eje distinto de cero es siempre diferente del eje. Esto es suficiente para usar Orthogonalize y Cross para obtener los vectores ortogonales deseados.
axisAngle[m_] := Module[
axis, ovec, nvec
,
axis, ovec =
Orthogonalize[1, -1, 1 #, Permute[#, Cycles[1, 3, 2]]] &@
Extract[m - Transpose[m], 3, 2, 3, 1, 2, 1];
(* nvec is orthogonal to axis and ovec: *)
nvec = Cross[axis, ovec];
axis,
Arg[Complex @@ (((m.ovec).# &) /@ ovec, nvec)]
]
El ángulo se calcula con Arg en lugar de ArcTan[x, y] aquí porque este último arroja un error para x = y = 0.
Aquí pruebo los resultados de la función para 100 matrices de rotación aleatorias:
testRotation[] := Module[
m, a, axis, ovec, nvec,
v = Normalize[RandomReal[0, 1, 3]],
α = RandomReal[-Pi, Pi], angle
,
m = RotationMatrix[α, v];
axis, angle = axisAngle[m];
Chop[
angle Dot[v, axis] - α
] === 0
]
And @@ Table[testRotation[], 100]
(* ==> True *)
En la prueba, tengo que tener en cuenta el hecho de que si la función axisAngle definido el vector del eje con el signo opuesto como el vector de prueba aleatorio, tengo que invertir el signo del ángulo de rotación. Esto es lo que el factor Dot[v, axis] lo hace.
Explicación de cómo el eje resulta de una matriz simétrica sesgada
Si $ vec v $ es el eje de la matriz de rotación $ R $, entonces tenemos $ R vec v = vec v $ y $ R ^ T vec v = vec v $ porque $ R ^ T $ es solo la rotación inversa. Por lo tanto, con $ a equiv R ^ TR $ como arriba, obtenemos
$$ a vec v = vec 0 $$
Ahora, la propiedad de simetría sesgada $ a ^ T = -a $, que se puede ver en su definición, significa que hay exactamente tres elementos de matriz independientes en $ a $. Se pueden organizar en forma de un vector 3D $ vec w $ que debe tener la propiedad $ a vec w = 0 $. Este vector se obtiene en el Extract línea de arriba. De hecho, $ a vec x = vec w times vec x $ para todos $ vec x $, y por lo tanto si $ a vec x = 0 $ entonces $ vec x paralelo vec w $. Por lo tanto, el vector $ vec v $ también es paralelo a $ vec w $, y este último es una representación válida del eje de rotación.
Edición 2: consideraciones de velocidad
Dado que el algoritmo anterior implica solo operaciones elementales que se pueden compilar, tiene sentido que una aplicación práctica de este enfoque utilice Compile. Entonces, la función podría definirse de la siguiente manera (manteniendo los valores de retorno organizados como se muestra arriba):
Clear[axisAngle1, axisAngle]
axisAngle1 =
Compile[m, _Real, 2,
Module[axis, ovec, nvec, tvec, w, w1,
tvec = m[[3, 2]] - m[[2, 3]], m[[3, 1]] - m[[1, 3]],
m[[2, 1]] - m[[1, 2]];
If[tvec == 0., 0., 0.,
#/Sqrt[#.#] &[#[[[email protected][N[Abs[#]]]]] &[
1/2 (m + 1, 0, 0, 0, 1, 0, 0, 0, 1)]],
If[Sum[m[[i, i]], i, 3] == 3, 0, Pi] 1, 1, 1,
axis = 1, -1, 1 tvec;
axis = axis/Sqrt[axis.axis];
w = tvec[[2]], tvec[[3]], tvec[[1]];
ovec = w - axis Dot[w, axis];
nvec = Cross[axis, ovec];
w1 = m.ovec;
axis, 1, 1, 1 ArcTan[w1.ovec, w1.nvec]
]
]
];
axisAngle[m_] := #1, Last[#2] & @@ axisAngle1[m]
Los resultados son los mismos que para la definición anterior de axisAngle, pero ahora obtengo una ejecución mucho más rápida como se puede ver en esta prueba:
tab = RotationMatrix @@ # & /@
Table[RandomReal[-Pi, Pi],
Normalize[RandomReal[0, 1, 3]], 100];
timeAvg =
Function[func,
Do[If[# > 0.3, Return[#/5^i]] & @@ [email protected][func, 5^i], i, 0,
15], HoldFirst];
timeAvg[axisAngle /@ tab]
(* ==> 0.000801259 *)
Esto es más de un orden de magnitud más rápido que la versión no compilada. quite Orthogonalize del código porque no lo encontré en la lista de funciones compilables.
Tenga en cuenta que Eigensystem tampoco está en esa lista.
Editar 3
La primera versión de axisAngle demostró las matemáticas básicas, pero la versión compilada axisAngle1 (junto con el redefinido axisAngle como envoltorio) es más rápido. Una cosa que faltaba era el tratamiento correcto del caso de borde donde la rotación es exactamente $ pi $ en ángulo.
Agregué esa solución solo a la versión compilada (axisAngle1) porque creo que esa es la versión más práctica de todos modos. El caso trivial del ángulo de rotación cero ya estaba incluido en la versión anterior.
Para explicar el código agregado, primero tenga en cuenta que para el ángulo $ pi $ no puede leer el eje de $ R ^ T – R $ porque la matriz resultante desaparece. Para sortear este caso singular, podemos usar el hecho geométrico de que una rotación de $ pi $ es equivalente a una inversión en el plano perpendicular al eje de rotación dado por el vector unitario $ vec n $. Por lo tanto, si formamos la suma de un vector $ vec v $ y su contraparte rotada $ pi $, las componentes transversales al eje de rotación se cancelan y el resultado es siempre paralelo al eje. En forma de matriz,
$$ (R + 1) vec v = 2 vec n ( vec n cdot vec v) = 2 left ( vec n vec n ^ T right) vec v $$
Dado que esto es válido para todos los vectores, es una identidad matricial. El lado derecho contiene una matriz $ vec n vec n ^ T $ que debe tener al menos una fila que no sea cero. Esta fila es proporcional a $ vec n ^ T $, por lo que puede leer el vector del eje directamente desde $ (R + 1) $, nuevamente sin ningún cálculo de valor propio.
Aquí hay una reimplementación un poco más corta del método de Jens, que usa el “método Pixar” (en sí mismo una modificación del método de Frisvad) para generar una base ortogonal:
axisAngleC = Compile[m, _Real, 2,
Module[ang, axis, ovec, nn, nvec, rm, s, w, w1, wm, xx, yy, zz,
axis = m[[3, 2]] - m[[2, 3]], m[[1, 3]] - m[[3, 1]], m[[2, 1]] - m[[1, 2]];
nn = Norm[axis];
If[nn == 0., (* singular case *)
ang = If[Total[Transpose[m, 1, 1]] == 3., 0., π];
rm = 0.5 (m + 1., 0., 0., 0., 1., 0., 0., 0., 1.);
axis = rm[[Ordering[Max /@ Abs[rm], -1][[1]]]];
axis /= Norm[axis],
(* non-singular case *)
axis /= nn; xx, yy, zz = axis;
s = 2 UnitStep[zz] - 1; w = -1/(s + zz); w1 = xx yy w;
ovec = 1 + s w xx xx, s w1, -s xx;
nvec = w1, s + w yy yy, -yy;
wm = m.ovec; ang = Arg[wm.ovec + I wm.nvec]];
Append[axis, ang]]];
Verifique con algunos pares de eje-ángulo generados aleatoriamente:
(* https://mathematica.stackexchange.com/a/111693 *)
rodrigues[th_, axis_?VectorQ] :=
First[LinearAlgebra`MatrixPolynomial[1, Sin[th], 2 Sin[th/2]^2,
-LeviCivitaTensor[3, List].Normalize[axis]]]
Max[Table[With[v = Normalize[RandomVariate[NormalDistribution[], 3]],
th = RandomReal[π],
Norm[axisAngleC[rodrigues[th, v]] - Append[v, th], ∞]], 10^3]] // Chop
0
Marque los casos singulares:
axisAngleC[IdentityMatrix[3]]
0., 0., 1., 0.
Max[Table[With[v = Normalize[RandomVariate[NormalDistribution[], 3]],
Norm[Cross[Most[axisAngleC[RotationMatrix[π, v]]], v], ∞]], 10^3]] // Chop
0
Aquí hay una versión para la entrada exacta que devuelve un inactivo RotationMatrix[]:
axisAngle[m_?OrthogonalMatrixQ] /; Det[m] == 1 :=
Module[ang, axis, ovec, nn, nvec, rm, s, w, w1, wm, xx, yy, zz,
axis = m[[3, 2]] - m[[2, 3]], m[[1, 3]] - m[[3, 1]], m[[2, 1]] - m[[1, 2]];
nn = Simplify[Norm[axis]];
If[nn == 0, (* singular case *)
ang = π Boole[Total[Diagonal[m]] < 3];
rm = (m + IdentityMatrix[3])/2;
axis = Normalize[Extract[rm, Ordering[Max /@ Abs[rm], -1]]],
(* non-singular case *)
xx, yy, zz = Simplify[axis/nn];
s = 2 UnitStep[zz] - 1; w = -1/(s + zz); w1 = xx yy w;
ovec = 1 + s w xx xx, s w1, -s xx;
nvec = w1, s + w yy yy, -yy;
wm = m.ovec; ang = Arg[Simplify[wm.ovec + I wm.nvec]]];
Inactive[RotationMatrix][ang, axis]]
Dado que esta pregunta todavía parece estar viva después de un tiempo, estoy dando una solución de Presentations, que vendo. Tiene una rutina RotationAngleAndAxis y aquí se usa en el ejemplo.
<< Presentations`
r = 0.966496, -0.214612, 0.14081, 0.241415,
0.946393, -0.214612, -0.0872034, 0.241415, 0.966496;
RotationAngleAndAxis[r]
(* 0.349066, 0.666667, 0.333333, 0.666667 *)
Aquí hay un breve ejemplo de su uso que verifica el resultado. El vector del eje está normalizado. Observe que al invertir el ángulo de rotación y el eje se obtiene una rotación equivalente.
rotations =
Table[RandomReal[-π, π],
Normalize[Array[RandomReal[-1, 1] &, 3]], 2]
(* 2.90598, -0.596373, -0.74938,
0.287697, -2.44158, 0.331763, -0.943343, 0.00605196 *)
rotationMatrices = RotationMatrix @@@ rotations;
anglesAndAxes = RotationAngleAndAxis /@ rotationMatrices
(* 2.90598, -0.596373, -0.74938, 0.287697, 2.44158, -0.331763,
0.943343, -0.00605196 *)
resultingMatrices = RotationMatrix @@@ anglesAndAxes;
rotationMatrices - resultingMatrices // Chop
Total[Flatten[%]]
(* 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 *)
(* 0 *)
La siguiente es una prueba sobre 1000 especificaciones de rotación inicial. Observe también que estoy comparando las matrices de rotación resultantes mediante un procedimiento un poco arriesgado.
rotations =
Table[RandomReal[-π, π],
Normalize[Array[RandomReal[-1, 1] &, 3]], 1000];
rotationMatrices = RotationMatrix @@@ rotations;
Timing[anglesAndAxes = RotationAngleAndAxis /@ rotationMatrices;]
resultingMatrices = RotationMatrix @@@ anglesAndAxes;
Total[Flatten[rotationMatrices - resultingMatrices // Chop]]
(* 0.873606, Null *)
(* 0 *)
Presentaciones también tiene un EulerAngles rutina que devolverá los dos conjuntos de ángulos de Euler para cualquiera de las posibles secuencias de rotaciones especificadas por cadenas. Por ejemplo, aquí están los dos conjuntos de ángulos de Euler del ejemplo para dos secuencias de rotación diferentes.
EulerAngles[r, "ZXZ"]
EulerAngles[r, "XYZ"]
(* -0.346633, 0.259587, 0.580661, 2.79496, -0.259587, -2.56093 *)
(* 0.244775, 0.0873143, 0.244775, -2.89682, 3.05428, -2.89682 *)
Comentarios y calificaciones
Al final de la web puedes encontrar las observaciones de otros gestores de proyectos, tú incluso eres capaz dejar el tuyo si te apetece.