Solución numérica de las ecuaciones de Euler y verificación analítica mediante funciones elípticas de Jacobi — Landau & Lifshitz §36–37
Se estudia la rotación libre —sin torques externos— de un cuerpo rígido completamente asimétrico cuyos momentos de inercia principales satisfacen $I_1 < I_2 < I_3$. El problema consiste en determinar analíticamente la estabilidad de la rotación alrededor de cada uno de los tres ejes principales, resolver numéricamente las ecuaciones de Euler y verificar que la rotación alrededor del eje de momento de inercia intermedio $I_2$ es inestable: una pequeña perturbación provoca un volteo periódico de $180°$ en la orientación del cuerpo, fenómeno conocido como el Teorema de la Raqueta de Tenis (o efecto Dzhanibekov), cuya solución analítica exacta involucra funciones elípticas de Jacobi (Landau & Lifshitz §37).
Un primer borrador de la integración de las ecuaciones de movimiento, del cálculo
de la polhodia y de la comparación entre distintos métodos numéricos, escrito en
Python, puede consultarse en
este cuaderno de Google Colab.
La simulación interactiva presentada en la sección 3 corresponde a una migración de ese
código a JavaScript, con la rotación 3D del cuerpo rígido y su animación construidas
sobre la API Canvas del navegador mediante integración de cuaterniones;
tanto la migración como la animación se desarrollaron con apoyo del modelo Claude.
Se considera un cuerpo rígido asimétrico con tres momentos de inercia principales distintos $I_1 < I_2 < I_3$, referidos a los ejes principales del cuerpo $\hat{\mathbf{e}}_1$, $\hat{\mathbf{e}}_2$, $\hat{\mathbf{e}}_3$. En ausencia de torques, el vector de velocidad angular en el frame del cuerpo es $\boldsymbol{\Omega} = (\Omega_1,\Omega_2,\Omega_3)$.
La ecuación de movimiento para el momento angular en el frame inercial es $d\mathbf{L}/dt = \mathbf{N}$. En el frame del cuerpo, para un cuerpo libre ($\mathbf{N}=\mathbf{0}$), con $L_i = I_i\Omega_i$ en los ejes principales:
El sistema (1) es no lineal y acoplado, pero admite dos integrales primeras exactas.
Las dos cantidades conservadas se obtienen de las leyes de conservación:
$$2E \;=\; I_1\Omega_1^2 + I_2\Omega_2^2 + I_3\Omega_3^2 \qquad (\text{energía cinética}), \tag{2}$$ $$M^2 \;=\; I_1^2\Omega_1^2 + I_2^2\Omega_2^2 + I_3^2\Omega_3^2 \qquad (\text{momento angular}^2). \tag{3}$$En términos de $M_i = I_i\Omega_i$, las ecuaciones (2) y (3) describen respectivamente un elipsoide con semiejes $\sqrt{2EI_i}$ y una esfera de radio $M$. El extremo del vector $\mathbf{M}$ queda confinado a la intersección de ambas superficies: la polhodia (Landau §37, Fig. 51).
| Régimen | Condición | Forma de la polhodia | Estabilidad |
|---|---|---|---|
| Cerca de $\hat{\mathbf{e}}_1$ | $M^2 \approx 2EI_1$ | Curva cerrada alrededor del polo en $\hat{\mathbf{e}}_1$ | ESTABLE |
| Separatriz (eje 2) | $M^2 = 2EI_2$ | Dos semielipses que cruzan los polos de $\hat{\mathbf{e}}_2$ | INESTABLE |
| Cerca de $\hat{\mathbf{e}}_3$ | $M^2 \approx 2EI_3$ | Curva cerrada alrededor del polo en $\hat{\mathbf{e}}_3$ | ESTABLE |
Para el eje intermedio $\hat{\mathbf{e}}_2$, sea $\boldsymbol{\Omega} = (\delta\Omega_1,\,\Omega_0,\,\delta\Omega_3)$ con $\delta\Omega \ll \Omega_0$. Linearizando (1):
De (1a) y (1c), con $\Omega_2 \approx \Omega_0 = \text{const}$ al primer orden:
$$I_1\,\delta\dot\Omega_1 = (I_2-I_3)\,\Omega_0\,\delta\Omega_3, \qquad I_3\,\delta\dot\Omega_3 = (I_1-I_2)\,\Omega_0\,\delta\Omega_1.$$Derivando la primera respecto al tiempo y sustituyendo:
$$\delta\ddot\Omega_1 = \frac{(I_2-I_3)(I_1-I_2)}{I_1 I_3}\,\Omega_0^2\;\delta\Omega_1. \tag{4}$$Dado que $I_1 < I_2 < I_3$, el producto $(I_2-I_3)(I_1-I_2) > 0$. La ecuación (4) tiene la forma $\delta\ddot\Omega_1 = +\lambda^2\delta\Omega_1$ con:
$$\lambda = \Omega_0\sqrt{\frac{(I_3-I_2)(I_2-I_1)}{I_1 I_3}} \;>\; 0. \tag{5}$$La solución $\delta\Omega_1 \sim e^{\lambda t}$ crece exponencialmente: el eje intermedio es inestable. El mismo análisis para $\hat{\mathbf{e}}_1$ y $\hat{\mathbf{e}}_3$ produce oscilaciones acotadas: ambos son estables.
La reducción de las ecuaciones de Euler mediante las integrales (2) y (3) conduce a la solución exacta (para $M^2 > 2EI_2$):
donde $\mathrm{sn}$, $\mathrm{cn}$, $\mathrm{dn}$ son las funciones elípticas de Jacobi con módulo:
El período temporal de la solución es:
$$T = 4K(k)\sqrt{\frac{I_1 I_2 I_3}{(I_3-I_2)(M^2-2EI_1)}}, \tag{9}$$donde $K(k) = \int_0^1 ds/\sqrt{(1-s^2)(1-k^2s^2)}$ es la integral elíptica completa de primera especie.
Cuando la rotación inicial se dirige cerca del eje intermedio $\hat{\mathbf{e}}_2$, se tiene $M^2 \approx 2EI_2$ y $k^2 \to 1^-$. En este límite:
Comportamiento de $K(k)$: cuando $k\to 1$, $K(k)\to\infty$ logarítmicamente. Por tanto, el período $T\to\infty$: la separatriz tiene período infinito.
Comportamiento de las funciones elípticas: en el límite $k=1$, $\mathrm{cn}(\tau,1) = \mathrm{sech}\,\tau$. La función $\mathrm{cn}(\tau,k)$ cambia de signo durante cada período, lo que implica que $\Omega_1(t)$ pasa de $+A$ a $-A$: el cuerpo experimenta un volteo de $180°$.
Un cuerpo rígido asimétrico que rota con $\boldsymbol{\Omega}_0 \approx \Omega_0\hat{\mathbf{e}}_2$ experimenta volteos periódicos de $180°$ alrededor de los ejes $\hat{\mathbf{e}}_1$ y $\hat{\mathbf{e}}_3$, con período $T \to \infty$ cuando $\boldsymbol{\Omega}_0 \to \Omega_0\hat{\mathbf{e}}_2$ exacto. Este comportamiento corresponde al régimen $k\to 1$ en la solución (6).
Las ecuaciones de Euler (1) se integran numéricamente mediante RK4. Para seguir la orientación del cuerpo en el espacio se integra simultáneamente el cuaternión unitario $\mathbf{q} = (q_w, q_x, q_y, q_z)$, evitando la singularidad del bloqueo del cardán presente en los ángulos de Euler.
El vector de estado completo es $\mathbf{y} = (\Omega_1,\Omega_2,\Omega_3,q_w,q_x,q_y,q_z)^\top$, gobernado por el sistema de 7 EDOs:
El cuaternión se renormaliza $|\mathbf{q}|=1$ al finalizar cada paso RK4 para controlar el drift numérico.
| Parámetro | Símbolo | Rango / valor | Rol físico |
|---|---|---|---|
| $I_1$ | Momento mínimo | $0.2$–$3.0$ kg·m² | Eje estable $\hat{\mathbf{e}}_1$ |
| $I_2$ | Momento intermedio | $0.2$–$3.0$ kg·m² | Eje inestable: $k\to 1$ cuando $\boldsymbol{\Omega}_0 \parallel \hat{\mathbf{e}}_2$ |
| $I_3$ | Momento máximo | $0.2$–$5.0$ kg·m² | Eje estable $\hat{\mathbf{e}}_3$ |
| $\boldsymbol{\Omega}_0$ | Condición inicial | Presets: tennis, sep, estable 1/3 | Determina $M^2$, $E$ y por tanto $k^2$ según (8) |
| $k^2$ | Módulo elíptico | $0 \leq k^2 \leq 1$ (calculado) | $k\to 0$: rotación casi pura; $k\to 1$: volteo de 180° |
| $h$ | Paso RK4 | $10^{-3}$ s (fijo) | Precisión $\mathcal{O}(h^4)$ para resolver la escala $1/|\boldsymbol{\Omega}|$ |
La simulación integra numéricamente el sistema (10) mediante RK4 con paso $h = 10^{-3}$ s y hasta 80 pasos por cuadro. Se visualizan simultáneamente: (i) la rotación 3D del cuerpo rígido con sus tres ejes principales coloreados; (ii) las componentes $\Omega_i(t)$ en función del tiempo, cuya forma corresponde a las funciones $\mathrm{cn}$, $\mathrm{sn}$, $\mathrm{dn}$ de la solución analítica (6); y (iii) la polhodia, trayectoria del extremo de $\boldsymbol{\Omega}$ en el elipsoide de inercia. Los controles permiten variar $I_1$, $I_2$, $I_3$ y seleccionar presets que ubican al sistema en distintos regímenes de $k^2$.