Movimiento de un cuerpo rígido.sagews

Movimiento de un cuerpo rígido

---

En el presente trabajo consideramos un sistema de partículas que es sometido a un tipo especial de movimiento conocido como _movimiento de un cuerpo rígido_ (entendemos por cuerpo rígido a un objeto material que puede moverse sin deformación). Tal movimiento conserva la distancia entre cualquier par de partículas en el sistema. Este tipo de movimiento nos conduce a importantes conceptos físicos y geométricos, tales como velocidad angular y momento de inercia.

Suponemos que hay $N$ partículas y que $F_i: J \times O \rightarrow \mathbb{R}^3,\; i=1,...,N$, es un sistema general de posición, velocidad y fuerzas dependientes del tiempo.

Definición: _Movimiento de un cuerpo rígido._ Una curva en $\mathbb{R}^{3N}$ $$\textbf{r}=(\textbf{r}_1,...,\textbf{r}_N):I\rightarrow\mathbb{R}^{3N}$$ se dice que es un movimiento de cuerpo rígido si $$\begin{equation} \vert \textbf{r}_i(t)-\textbf{r}_j(t)\vert=\vert\textbf{r}_i(0)-\textbf{r}_j(0)\vert, \end{equation}$$ para todo $i,j=1,...,N$ y para todo $t\in I$.

En general, sólo nos ocuparemos de movimientos de cuerpos rígidos que satisfagan las ecuaciones movimiento $m_i\ddot{\textbf{r}}_i=F_i(t,\textbf{r},\dot{\textbf{r}})$. Por lo tanto, además de satisfacer las ecuaciones de movimiento, los $\textbf{r}_i$ deben mantener sus respectivas distancias (la misma que sus distancias iniciales). Es evidente que no todos los sistemas de fuerzas permiten movimientos de cuerpos rígidos, pero para aquellos que lo hacen, el movimiento $\textbf{r} = (\textbf{r}_1,..., \textbf{r}_N)$ tiene una forma muy particular, como muestra el siguiente teorema importante.

Teorema:Sea $\textbf{r}=(\textbf{r}_1,...,\textbf{r}_N):I\rightarrow\mathbb{R}^{3N}$ el movimiento de un cuerpo rígido. Entonces existe una función $Q:I\subset \mathbb{R}\rightarrow SO(3)$, con $SO(3)=\lbrace O\in \mathbb{R}^{3\times 3} : O^T O=I; \; det(O)=1\rbrace$, as'i $Q$ es una matriz ortogonal: $$\begin{equation} Q(t)Q(t)^T=I=Q(t)^T Q(t) \end{equation}$$ para cada $t\in I$ con $Q(0)=I$ y $N$-vectores $\textbf{u}_i\in \mathbb{R}^{3}$ $i=1,...,N$ que satisfacen $\sum_{i=1}^{N}m_i\textbf{u}_i=0$, tal que $$\begin{equation} \textbf{r}_i(t) = Q(t)\textbf{u}_i + \textbf{R}(t), \end{equation}$$ para cada $i=1,...,N$ y para todo $t\in I$. Donde $\textbf{R}(t)$ es la posición del centro de masa en el tiempo $t$.

El teorema dice que, en un movimiento de cuerpo rígido, cada partícula del cuerpo sigue el movimiento del centro de masa, mientras que gira alrededor de este centro de masa. Este es el contenido de las ecuaciones $$\begin{equation} \textbf{r}_i = Q\textbf{u}_i + \textbf{R}, \end{equation}$$ para $i = 1,. . . , N$, que expresan los vectores de posición de las partículas en términos del vector posición $\textbf{R} = \textbf{R} (t)$ del centro de masa, la matriz de rotación $Q = Q (t)$, y los vectores $\textbf{u}_i$ (constantes), que son las posiciones iniciales de los cuerpos en relación con el centro de la masa inicial ($\textbf{r}_i (0) = \textbf{u}_i + \textbf{R} (0)$).

---

Ecuaciones diferenciales de un cuerpo rígido

---

El teorema anterior nos muestra que cualquier movimiento de un cuerpo rígido $\textbf{r}: I \rightarrow\mathbb{R}^{3N}$ de $N$ partículas está completamente determinado por el movimiento de su centro de masa $\textbf{R}: I \rightarrow\mathbb{R}^3$ y una cierta familia $\lbrace Q (t)\rbrace_{t\in I}$ de matrices ortogonales. Si ahora necesitamos que el movimiento de un cuerpo rígido también satisfaga las ecuaciones de movimiento (segunda ley de Newton), entonces $\textbf{R}$ y $Q$ satisfacen un sistema de ecuaciones diferenciales, al que llamamos _ecuaciones diferenciales de cuerpos rígidos_. En esta sección obtendremos estas ecuaciones y, de paso, explicaremos el significado geométrico y dinámico de $Q$.

Para la obtención de las ecuaciones necesitaremos introducir algunos conceptos importantes, tales como la velocidad angular y el operador de momento de inercia.

A las fuerzas individuales $F_i=F_i(t,\textbf{r}_1,...,\textbf{r}_N,\dot{\textbf{r}}_1,...,\dot{\textbf{r}}_N)$ las podemos reescribir en términos de $\textbf{R}$, $Q$ y sus derivadas usando $$\textbf{r}_i=Q\textbf{u}_i+\textbf{R}$$ y en consecuencia $$\begin{equation} \dot{\textbf{r}}_i=\dot{Q}\textbf{u}_i+\dot{\textbf{R}}. \end{equation}$$

Habiendo hecho esto, la fuerza total $\textbf{F}$ y el torque $\textbf{T}$ se convierten en funciones de $\textbf{R}$, $Q$ y sus derivadas. Entonces, ya que la posición $\textbf{R}$ del centro de masa satisface $$M\ddot{\textbf{R}}=\textbf{F},$$ podemos interpretar esta ecuación como una parte del sistema de ecuaciones diferenciales para $\textbf{R}$ y $Q$ que estamos buscando.

Las ecuaciones diferenciales para $Q$ involucran su derivada $\dot{Q}$ y como $\lbrace Q (t)\rbrace_{t\in I}$ es una familia 1-paramétrica de rotaciones, está compuesta por una familia 1-paramétrica $\lbrace \Omega(t)\rbrace_{t\in I}$ de matrices antisimétricas.

Proposición: Sea $Q:J\rightarrow O(n)$ una familia 1-paremétrica (suave) de matrices ortogonales $n\times n$. Entonces, existe una familia 1-paramétrica (suave) $\Omega:J\rightarrow o(n)$ de matrices antisimétricas $n\times n$ tal que $Q$ satisface la ecuación diferencial $$\dot{Q}=Q\Omega.$$ Además, si $Q(0)=I$, entonces $det(Q(t))=1$, es decir, $Q(t)$ es una matriz de rotación, para cada $t\in J$.

Inversamente, si $\Omega: J\rightarrow o(n)$ es dada, una familia 1-paramétricas de matrices antisimétricas, entonces una solución $Q$ del problema de valores iniciales $$\dot{Q}=Q\Omega$$ $$Q(0)=I$$ es una familia 1-paramétricas de matrices de rotación.

Observamos que el resultado anterior nos permite considerar a $Q$ y $\Omega$ como intercambiables en términos de determinar el movimiento, ya que conociedo uno conocemos el otro. Por otro lado, debido a la condición de antisimetría, $\Omega$ está completamente determinada por tres funciones: $\omega_1, \omega_2, \omega_3: J \rightarrow \mathbb{R}$. Es decir, $\Omega$ tiene la forma $$\begin{equation} \Omega=\begin{bmatrix}{0}&{-\omega_3}&{\omega_2}\\{\omega_3}&{0}&{-\omega_1}\\ {-\omega_2}&{\omega_1}&{0}\end{bmatrix}, (01) \end{equation}$$

Hay una razón para expresar $\Omega$ en términos de $\omega_1, \omega_2, \omega_3$ de esta manera muy particular. Es decir, si $$\begin{equation} \omega\equiv(\omega_1, \omega_2, \omega_3), (02) \end{equation}$$ entonces para cualquier vector $v = (v_1, v_2, v_3)$, tenemos $\Omega v = \omega \times v$. Por lo tanto, la acción de la matriz antisimétrica $\Omega$ sobre los vectores $v$ es la misma que la acción del producto vectorial de $\omega$ en $v$. Esta relación entre $\Omega$ y $\omega$ vale en general: cualquier matriz antisimétrica $B$, $3\times 3$ tiene su acción $v \rightarrow Bv$ representado por la acción del producto vectorial $v \rightarrow b \times v$, con el vector fijo $b= (-B_{23}, B_{13}, -B_{12})$. Para la mecánica, la representación en producto vectorial de $\Omega$ es importante para la interpretación geométrica del movimiento.

Definición: Velocidad angular. Para el movimiento de un cuerpo rígido, (dependiente del tiempo) el vector $\omega$ en (02) se denomina _velocidad angular_ y a la matriz antisimétrica $\Omega$ (dependiente del tiempo) en (01) se denomina _operador de velocidad angular_.

El punto de vista actual es que $\omega$ es una de las incógnitas fundamentales en el movimiento de un cuerpo rígido y buscamos la ecuación diferencial que la determina. Con la determinación de $\omega$, obtenemos $\Omega$ y de esto tenemos las rotaciones $Q$.

Podemos escribir fácilmente la velocidad de la $i$-ésima partícula en términos de la velocidad angular de la siguiente manera: como $$\textbf{r}_i=Q \textbf{u}_i+\textbf{R}$$ tenemos $$\begin{darray}{rcl} \dot{\textbf{r}}_i &=& \dot{Q}\textbf{u}_i+\dot{\textbf{R}} \nonumber \\ &=& Q\Omega \textbf{u}_i +\dot{\textbf{R}}\nonumber \\ &=& Q(\omega\times \textbf{u}_i)+\dot{\textbf{R}} \nonumber \end{darray}$$

Esto expresa la velocidad de la $i$-ésima partícula del cuerpo como la suma de una velocidad de traslación $\textbf{R}$ y una velocidad de rotación $Q (\omega \times \textbf{u}_i) = Q\omega \times Q\textbf{u}_i$.

El primero es visto como una traslación instantánea, en la dirección $\dot{\textbf{R}}$, de todos las partículas del sistema, mientras que el último es visto como una rotación instantánea de todas las partículas del sistema sobre el eje a través de $Q\omega$.

Usando las expresiones anteriores, también podemos reescribir el momento angular total $\textbf{L}$ del sistema.

Nota: Una mejor simplificación en $\textbf{L}$ se debe a que $\sum_{i=1}^N m_i \textbf{u}_i=0$. Por lo tanto,

$$\begin{darray}{rcl} \textbf{L} &=& \sum_{i=1}^N m_i \textbf{r}_i \times \dot{\textbf{r}}_i \nonumber \\ &=& \sum_{i=1}^N m_i[Q \textbf{u}_i + \textbf{R}]\times [Q(\omega\times\textbf{u}_i)+\dot{\textbf{R}}] \nonumber \\ &=& Q \left(\sum_{i=1}^N m_i \textbf{u}_i\times (\omega\times \textbf{u}_i)\right)+ M\textbf{R}\times \dot{\textbf{R}} \nonumber \end{darray}$$

Aquí hemos utilizado el hecho de que las matrices de rotación se pueden distribuir en el producto vectorial: $Q (v \times w) = Qv \times Qw$. Si denotamos por $A\omega$ a la expresión en los paréntesis grandes en la última línea de arriba, obtenemos la siguiente expresión para el momento angular total del sistema:

**Momento angular de un cuerpo rígido:** $$\begin{equation} \textbf{L}= Q A\omega +M\textbf{R} \times \dot{\textbf{R}}, (03) \end{equation}$$

La notación $A\omega$ introduce un operador importante $A$ en la teoría. Este operador es definido por lo siguiente:

**Definición:** _Operador de inercia $A$._ Sean $m_1,...,m_N$ números positivos dados (las masas) y $\textbf{u}_1,...,\textbf{u}_N$ vectores dados (las posiciones iniciales de las partículas relativas al centro de masa), con $\sum_{i=1}^N m_i\textbf{u}_i=0$. Estos datos comprenden un "cuerpo rígido". Asociamos a este cuerpo rígido el operador lineal $A:\mathbb{R}^3\rightarrow \mathbb{R}^3$ definido por $$\begin{equation} Av\equiv \sum_{i=1}^N m_i\textbf{u}_i \times (v\textbf{u}_i), (04) \end{equation}$$ para $v\in \mathbb{R}^3$. Este operador se conoce como operador de inercia o tensor de inercia para el cuerpo rígido.

El operador $A$ captura las características geométricas del cuerpo rígido que son importantes para la naturaleza de su movimiento. Observemos que $A$ no depende del tiempo, pero si lo hace de las posiciones iniciales $\textbf{u}_1,...,\textbf{u}_N$ y de las masas de las partículas.

La ecuación diferencial central para $\omega$ proviene de reescribir la ecuación del momento angular $\dot{\textbf{L}}=\textbf{T}$ en término de la notación usada arriba. Primero notemos que diferenciando ambos lados de la ecuación (03) tenemos $$\begin{darray}{rcl} \dot{\textbf{L}} &=& QA\dot{\omega}+ \dot{Q}A\omega+ M\textbf{R}\times \ddot{\textbf{R}} \nonumber \\ &=& Q(A\dot{\omega}+\omega\times A\omega)+ M\textbf{R}\times \ddot{\textbf{R}}, (05) \end{darray}$$

Entonces calculando el torque total obtenemos $$\begin{darray}{rcl} \textbf{T} &=& \sum_{i=1}^N\textbf{r}_i\times F_i \nonumber \\ &=& \sum_{i=1}^N (Q\textbf{u}_i+\textbf{R})\times F_i \nonumber \\ &=& \left( \sum_{i=1}^N Q\textbf{u}_i \times F_i\right) +\textbf{R}\times \textbf{F}, (06) \end{darray}$$

Sustituyendo estas expresiones en $\dot{\textbf{L}}=\textbf{T}$ y reduciendo cada lado usando que $M \ddot{\textbf{R}}=\textbf{F}$ da $$Q(A\dot{\omega}+\omega\times A\omega)=\sum_{i=1}^N Q\textbf{u}_i\times F_i,$$ i.e., $$Q(A\dot{\omega}+\omega\times A\omega)=Q\left( \sum_{i=1}^N \textbf{u}_i\times Q^T F_i\right).$$

Aquí, la factorización de $Q$ en el lado izquierdo de la ecuación se deduce de la identidad general: $Qv\times Qw=Q(v\times w)$. Canculando $Q$ en la última ecuación de arriba nos da la ecuación diferencial deseada para $\omega$. Poniendo junto todo lo último llegamos a lo siguiente:

\textbf{Las ecuaciones de movimiento de un cuerpo rígido:} $$\begin{equation} M\ddot{\textbf{R}}=\textbf{F}, (07) \end{equation}$$

$$\begin{equation} A\dot{\omega}+\omega\times A \omega=\sum_{i=1}^N \textbf{u}_i\times Q^T F_i, (08) \end{equation}$$$$\begin{equation} \dot{Q}=Q\Omega, (09) \end{equation}$$

Las incógnitas aquí son $\textbf{R}$, $\omega$ y $Q$ (recordar que $\Omega$ en su expresión involucra a $\omega$ como en (01). Por lo tanto, lo anterior es un sistema de ecuaciones diferenciales con 15 ecuaciones y 15 incógnitas $R_i$, $\omega_i$, $Q_{ij}$, $i,j=1,2,3$. Las \emph{condiciones iniciales} son (en término de las incógnitas originales)

$$\textbf{R}(0)=\sum_{i=0}^N \frac{m_i}{M}\textbf{r}_i(0)$$ $$\dot{\textbf{R}}(0)=\sum_{i=0}^N \frac{m_i}{M}\dot{\textbf{r}}_i(0)$$ $$Q(0)=I$$ $$A\omega(0)=\textbf{L}(0)-M\textbf{R}(0)\times \dot{\textbf{R}}(0).$$

La última condición inicial sólo determina el valor inicial de $\omega$ cuando el operador de inercia $A$ es invertible. El caso en que $A$ no es invertible se produce sólo cuando todos los cuerpos se encuentran inicialmente en una línea recta. Esto se llama el caso degenerado.

La primera ecuación diferencial (07) determina el movimiento del centro de masa, mientras que la ecuación diferencial (08) determina el vector de velocidad angular $\omega$, que a su vez da $\Omega$ y por lo tanto, a través de la última ecuación diferencial (09), obtenemos la rotación $Q$ de las partículas sobre el centro de masa.

Para ser más específicos sobre el movimiento real, vamos a ver algunos casos especiales y ejemplos a continuación. Antes de hacerlo discutimos unos resultados adicionales para el movimiento de cuerpo rígido en general.

---

Energía cinética y momento de inercia

---

El movimiento general del cuerpo rígido que consideramos se caracteriza por la configuración inicial de las partículas que comprenden el sistema, es decir, por el de posiciones iniciales $\textbf{u}_1,...,\textbf{u}_N$ (con respecto al centro de masa inicial $\textbf{R}(0)$) y por las masas $m_1,. . ., m_N$, tal que $\sum_{i=1}^N m_i\textbf{u}_i=0$ y $\sum_{i=1}^N m_i=M$. Este sistema es el del "cuerpo rígido" y se mueve como una unidad tras el centro de masa, mientras que gira alrededor de él. Si bien el sistema es un sistema discreto, de un número finito de partículas (aunque posiblemente un gran número), la discusión y los conceptos se extienden al caso continuo, con un número infinito de partículas, que es lo que normalmente consideramos como un cuerpo rígido.

La influencia de la "forma" del cuerpo rígido en el movimiento se codifica en el operador de inercia $A$ y la siguiente proposición describe la naturaleza de $A$ y muestra cómo entra en la expresión para la energía cinética del sistema.

**Proposición:** Sea $A$ un operador de inercia definido por $$Av=\sum_{i=1}^N m_i\textbf{u}_i\times (v\times \textbf{u}_i),$$ para $v\in \mathbb{R}^3$. Entonces $$\begin{equation} Av\cdot w= \sum_{i=1}^N m_i(w\times \textbf{u}_i)\cdot (v\times \textbf{u}_i), (10) \end{equation}$$ para todo $v,w\in \mathbb{R}^3$. En particular $$\begin{equation} Av\cdot v=\sum_{i=1}^N m_i\vert v\times \textbf{u}_i\vert ^2\geq 0, (11) \end{equation}$$ para todo $v\in \mathbb{R}^3$. Por lo tanto $A$ es una matriz simétrica, semidefinida positiva. Es diagonalizable y tiene tres autovalores $I_1$, $I_2$, $I_3$ no negativos: $$0\leq I_1 \leq I_2\leq I_3$$ y los tres autovectores correspondientes $$Ae_i=I_ie_i, \;\;\;\;i=1,2,3$$ pueden ser elegidos de modo que sean ortogonales $e_i\cdot e_j=0$, $i\neq j$.

Si $K=\sum_{i=1}^N m_i\vert \dot{\textbf{r}}\vert ^2$ es la _energía cinética total del sistema_, entonces $$\begin{equation} K=\frac{1}{2}\omega\cdot \omega +\frac{1}{2}M\vert\dot{\textbf{R}}\vert ^2, (12) \end{equation}$$ Por lo tanto, la energía cinética se divide en dos partes. La primera parte, llamada _energía cinética rotacional_, es debida a la rotación instantánea, o giro, alrededor del eje determinado por el vector velocidad angular $\omega$. La segunda parte es llamada _energía cinética traslacional_.

Al igual que en la proposición, siempre vamos a etiquetar y ordenar los autovalores de $A$ como $I_1\leq I_2\leq I_3$ y seleccionaremos los tres autovectores correspondientes $e_1, e_2, e_3$, que son ortogonales, de longitud uno, $\vert e_i \vert = 1$, $i = 1, 2, 3$, y tal que $\lbrace e_1, e_2, e_3\rbrace$ está orientado positivamente. Cuando todos los autovalores son distintos esta selección es única. Cuando sólo hay dos autovalores distintos, digamos $I_1 = I_2 \neq I_3$, la selección de dos vectores ortonormales $e_1, e_2$ desde el espacio característico $E_{I_{1}}$ puede hacerse en un número infinito de formas, pero habiendo hecho una elección, entonces $e_3$ se determina de forma única. Cuando todos los autovalores son los mismos el espacio característico $E_{I_{1}}$ está en $\mathbb{R}^3$ y así cualquier orientación positiva puede ser seleccionada en la base ortonormal $\lbrace e_1, e_2, e_3\rbrace$ para $\mathbb{R}^3$.

Como es habitual, identificamos la base estándar de $\mathbb{R}^3$ por $\lbrace \xi_1, \xi_2, \xi_3\rbrace$ donde $$\xi_1 = (1, 0, 0),\;\; \xi_2 = (0, 1, 0),\;\; \xi_3 = (0, 0, 1).$$

Por lo tanto, el operador de inercia $A: \mathbb{R}^3\rightarrow\mathbb{R}^3$ se identifica con la matriz $3\times 3$ $A = \lbrace A_{ij}\rbrace_{i, j = 1,2,3}$, donde $$A_{ij} \equiv A\xi_j\cdot \xi_i.$$

Usando la definición de $A$, es fácil obtener la siguiente representación explícita de $A$.

**Proposición:** El operador inercia $A = \lbrace A_{ij}\rbrace_{i, j = 1,2,3}$ tiene entradas $$\begin{equation} A_{ij}=\sum_{k=1}^N m_k [\vert \textbf{u}_k \vert^2 \delta_{ij}-u_{ki}u_{kj}], (13) \end{equation}$$ donde $\delta_{ij}$ es la delta Kronecker (es decir, $\delta_{ij} = 0$, para $i\neq j$, $\delta_{ii} = 1$), y $$\textbf{u}_k = (u_{k1}, u_{k2}, u_{k3}),$$ para $k = 1,. . .,N$. Alternativamente, sin el uso la delta de Kronecker esto puede escribirse como $$\begin{equation} A_{ii}=\sum_{k=1}^N m_k[\vert \textbf{u}_k\vert^2 - u_{ki}^2], (14) \end{equation}$$ $$\begin{equation} A_{ij}=-\sum_{k=1}^N m_k u_{ki}u_{kj} \;\;\;\; i\neq j, (15) \end{equation}$$

Ejemplo:

Supongamos que hay tres partículas de igual masa, digamos $m_k=1$, $k=1,2,3$, que componen el cuerpo rígido, y las posiciones iniciales de esas partículas son $$\textbf{r}_1(0)=(0,0,0)$$ $$\textbf{r}_2(0)=(1,0,0)$$ $$\textbf{r}_3(0)=(0,1,0)$$

Para obtener el tensor de inercia, debemos primero calcular las posiciones relativas $\textbf{u}_k$, $k=1,2,3$. Para estos datos el centro de masa inicial es $$\textbf{R}(0)=\frac{1}{3}(0,0,0)+\frac{1}{3}(1,0,0)+\frac{1}{3}(0,1,0)=(\frac{1}{3},\frac{1}{3},0).$$ Entonces tenemos $$\textbf{u}_1=\textbf{r}_1(0)-\textbf{R}(0)=(-\frac{1}{3},-\frac{1}{3},0)$$ $$\textbf{u}_2=\textbf{r}_2(0)-\textbf{R}(0)=(\frac{2}{3},-\frac{1}{3},0)$$ $$\textbf{u}_3=\textbf{r}_3(0)-\textbf{R}(0)=(-\frac{1}{3},\frac{2}{3},0).$$

Para calcular el operador de inercia $A$ usamos las fórmulas (14) y (15). Tengamos en cuenta que la fórmula (14) dice que $A_{11} = \sum_{k=1}^N m_k (u_{k2}^2+u_{k3}^2)$ entonces en este ejemplo

$$A_{11}=\frac{1}{9}+\frac{1}{9}+\frac{4}{9}=\frac{2}{3}$$

Del mismo modo $A_{22}$ implica la suma de los cuadrados de la primera y tercera componentes de $u_k$, mientras $A_{33}$ implica la suma de los cuadrados de la primera y segunda componenete. Entonces tenemos,

$$A_{22}=\frac{1}{9}+\frac{4}{9}+\frac{1}{9}=\frac{2}{3}$$ $$A_{33}=\frac{1}{9}+\frac{1}{9}+\frac{4}{9}+\frac{1}{9}+\frac{1}{9}+\frac{4}{9}=\frac{4}{3}.$$

Para calcular $A_{ij}$ para $i\neq j$ notemos por ejemplo que $u_{k3}=0$ para todo $k$ y entonces $$A_{13}=0=A_{23}$$

Para $A_{12}$ tenemos $$A_{12}=-[u_{11}u_{12} + u_{21}u_{22} + u_{31}u_{32}] = -[\frac{1}{9}-\frac{2}{9}-\frac{2}{9}]=\frac{1}{3}.$$

Como la matriz $A$ es simétrica, estos son los único calculos que necesitamos. Por lo tanto $$A=\begin{bmatrix}{\frac{2}{3}}&{\frac{1}{3}}& {0}\\{\frac{1}{3}}&{\frac{2}{3}}& {0}\\ {0}&{0}&{\frac{4}{3}}\end{bmatrix}.$$

Una vez hecho esto, nos es fácil calcular los autovalores y autovectores de $A$: $$I_1=\frac{1}{3},\;\;\;\;\;e_1=\frac{1}{\sqrt{2}}(-1,1,0)$$ $$I_2=1,\;\;\;\;\;e_2=\frac{1}{\sqrt{2}}(1,1,0)$$ $$I_3=\frac{4}{3},\;\;\;\;\;e_3=\frac{1}{\sqrt{2}}(0,0,1)$$

#Podemos chequear este ejemplo con el siguiente algoritmo

M=3;
r=[vector([0,0,0]), vector([1,0,0]), vector([0,1,0]) ] #posiciones de las partículas respecto a un origen.
N=len(r) #cantidad de partículas
CentroMasa=sum([vector(h)/N for h in r ]) #calculamos el centro de masa.
U=[vector(h)-CentroMasa for h in r] #posiciones relativas al centro de masa.
R0=[1/3,1/3,0] #condición inicial
show(U)

[$\displaystyle \left(-\frac{1}{3},\,-\frac{1}{3},\,0\right)$, $\displaystyle \left(\frac{2}{3},\,-\frac{1}{3},\,0\right)$, $\displaystyle \left(-\frac{1}{3},\,\frac{2}{3},\,0\right)$]

#matriz inercia
A11 = sum([u.norm()^2-u[0]^2 for u in U])*M/N
A12 = -sum([u[0]*u[1] for u in U])*M/N
A13 = -sum([u[0]*u[2] for u in U])*M/N
A21 = -sum([u[1]*u[0] for u in U])*M/N
A22 = sum([u.norm()^2-u[1]^2 for u in U])*M/N
A23 = -sum([u[1]*u[2] for u in U])*M/N
A31 = -sum([u[2]*u[0] for u in U])*M/N
A32 = -sum([u[2]*u[1] for u in U])*M/N
A33 = sum([u.norm()^2-u[2]^2 for u in U])*M/N
A = Matrix([[A11,A12,A13],[A21,A22,A23],[A31,A32,A33]])
autovalores = A.eigenvalues()
show(A)
show(autovalores)

$\displaystyle \left(\begin{array}{rrr} \frac{2}{3} & \frac{1}{3} & 0 \\ \frac{1}{3} & \frac{2}{3} & 0 \\ 0 & 0 & \frac{4}{3} \end{array}\right)$

[$\displaystyle 1$, $\displaystyle \frac{1}{3}$, $\displaystyle \frac{4}{3}$]

Para interpretar aún más la importancia geométrica de los autovalores y autovectores del operador inercia, nos fijamos en el concepto de los momentos de inercia de un cuerpo rígido.

Cada matriz simétrica B da lugar a una función $f (v) \equiv \frac{(Bv \cdot v)}{(v \cdot v)}$ cuyos valores son, respectivamente, autovalores máximos y mínimos de $B$ y estos valores extremos son asumidos en los respectivos autovectores $v$ de $B$. En estos casos, la función $f$ se conoce como el _momento de inercia de la función_.

**Definición:**_Momento de Inercia._ La función $I: \mathbb{R}^3-\lbrace 0\rbrace \rightarrow \mathbb{R}$ definida por $$I(v)=\dfrac{v\cdot Av}{v\cdot v},$$ se llama _momento de inercia de la función_. El número $I (v)$ se llama _momento de inercia del sistema_ sobre el eje a través del centro de masa en la dirección de $v$. Debemos tener en cuenta que $I(cv) = I (v)$ para todos los escalares no nulos $c$ y por eso nos da el mismo valor para todos los vectores no nulos en el eje determinado por $v$. Los autovalores $I_1\leq I_2 \leq I_3$ de $A$ son llamados los _momentos principales de inercia_.

Se desprende de la teoría general, que los momentos principales de inercia son los valores de $I$ en los autovectores de $A$: $$I_i = I (ei),$$ $i = 1, 2, 3$, y que $I_1 \leq I (v) \leq I_3$, para cada $v \neq 0$ en $\mathbb{R}^3$. Por lo tanto, $I_1$ es el momento de inercia _mínimo_ del sistema. Así mismo, $I_3$ es el momento de inercia _máximo_ del sistema.

Esto se ilustra por la configuración elemental en el ejemplo anterior. Rotar alrededor del eje a través de $e_1$ da el menor momento de inercia $I_1 = \frac{1}{3}$, mientras que un giro alrededor del eje a través de $e_2$ da como resultado un mayor momento de inercia $I_2 = 1$. Esto corresponde al hecho de que los cuerpos están más cerca del eje de revolución en el primer caso que en el segundo (y todas las masas son iguales). El mayor momento de inercia $I_3 = \frac{4}{3}$ se produce para una revolución alrededor del eje a través de $e_3$, donde los cuerpos se separan más alejado del eje de revolución.

Una expresión alternativa para el momento de inercia $I (v)$ ayuda a clarificar la idea que se alude en el último párrafo. Esta expresión surge de observar que $\vert v \times \textbf{u}_i \vert$ es el área del paralelogramo determinado por $v$ y $\textbf{u}_i$. Esta zona es $\vert v \vert d_i (v)$, donde $d_i (v)$ denota la distancia desde la punta de $\textbf{u}_i$ por $v$.

Así, encontramos $$Av \cdot v =\sum_{i=1}^N m_i\vert v\times \textbf{u}_i\vert^2=\vert v \vert ^2 \sum_{i=1}^N m_i d_i(v)^2.$$

Esto nos dice que el momento de inercia de la función también se puede calcular como $$\begin{equation} I(v)= \sum_{i=1}^N m_i d_i(v)^2 \end{equation}$$

Esta es la definición tradicional del momento de inercia del sistema alrededor del eje $v$.

**Nota:** Aquí, $u_i$, $i = 1,. . . , N$, son las posiciones de las partículas con respecto al centro de la masa inicial y por lo tanto el momento de inercia de arriba se da alrededor de un eje a través del centro de masa.

El nombre para el momento de inercia surge de la forma en que esta cantidad entra en la expresión para la energía cinética del sistema. Reescribiendo la energía cinética $K$ de la ecuación (12), en términos de $I$ da la siguiente

**Energía cinética del cuerpo rígido:** $$\begin{equation} K=\frac{1}{2}I(\omega)\vert \omega\vert ^2 +\frac{1}{2}M \vert \dot{\textbf{R}}\vert^2. \end{equation}$$

Por lo tanto, con respecto a la velocidad angular $\omega$ del cuerpo, $I (\omega)$ desempeña el mismo papel que la masa con respecto a la velocidad lineal. Cuanto mayor es el valor de $I (\omega)$, más difícil es detener el giro del cuerpo alrededor del eje por $\omega$.

Hasta aquí hemos tratado de desarrollar una comprensión intuitiva de la velocidad angular $\omega$ y su papel en el movimiento de un cuerpo rígido. Una comprensión más precisa vendrá de la solución de las ecuaciones de movimiento en un número de casos especiales.

#El algoritmo que sigue permite aproximar soluciones de las ecuaciones diferenciales de un cuerpo rígido.
a=1
b=2
c=3
g=9.8
M=10;
L=srange(0,1.5,.5)
r=[ [a*i,b*j,c*k] for i in L for j in L for k in L] #posiciones de las partículas respecto a un origen.
N=len(r) #cantidad de partículas
F=vector([0,0,-M/N*g]) # Vamos suponer que la fuerza es la gravedad
CentroMasa=sum([vector(h)/N for h in r ]) #calculamos el centro de masa.
U=[vector(h)-CentroMasa for h in r] #posiciones relativas al centro de masa.
R0=[0,0,0] #condición inicial
Rpunto0=[0,0,0] #condición inicial
Q0=[1,0,0,0,1,0,0,0,1] #condición inicial
omega0=[0,1,1] #condición inicial

#matriz inercia
A11 = sum([u.norm()^2-u[0]^2 for u in U])*M/N
A12 = -sum([u[0]*u[1] for u in U])*M/N
A13 = -sum([u[0]*u[2] for u in U])*M/N
A21 = -sum([u[1]*u[0] for u in U])*M/N
A22 = sum([u.norm()^2-u[1]^2 for u in U])*M/N
A23 = -sum([u[1]*u[2] for u in U])*M/N
A31 = -sum([u[2]*u[0] for u in U])*M/N
A32 = -sum([u[2]*u[1] for u in U])*M/N
A33 = sum([u.norm()^2-u[2]^2 for u in U])*M/N
A = matrix([[A11,A12,A13],[A21,A22,A23],[A31,A32,A33]])
Ainversa = A.inverse()
autovalores = A.eigenvalues()
autovalores

<string>:1: UserWarning: Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues.

[21.6666666666667, 16.6666666666667, 8.33333333333333]

#momentos principales de inercia

e1=vector([1,0,0]) #autovectores de A
e2=vector([0,1,0])
e3=vector([0,0,1])
I1=e1.dot_product(A*e1)/e1.dot_product(e1)
I2=e2.dot_product(A*e2)/e2.dot_product(e2)
I3=e3.dot_product(A*e3)/e3.dot_product(e3)

T =  ode_solver()
def CuerpoRigido(t,y):
    omega=vector([y[6],y[7],y[8]]) #velocidad angular
    Omega=Matrix([[0,-y[8],y[7]],[y[8],0,-y[6]],[-y[7],y[6],0]]) #operador velocidad angular
    Q=Matrix([[y[9],y[10],y[11]],[y[12],y[13],y[14]],[y[15],y[16],y[17]]]) #matriz de rotación 
    Qpunto=Q*Omega
    Suma=sum([u.cross_product(Q.transpose()*F)    for u in U] )
    wpunto=Ainversa*(- omega.cross_product(A*omega)+Suma)
    return [y[3],y[4],y[5],F[0]/M,F[1]/M,F[2]/M,wpunto[0],wpunto[1],wpunto[2],Qpunto[0,0],Qpunto[0,1],Qpunto[0,2],Qpunto[1,0],Qpunto[1,1],Qpunto[1,2],Qpunto[2,0],Qpunto[2,1],Qpunto[2,2]]
T.function=CuerpoRigido
T.algorithm="rk4"

#Condiciones iniciales
T.ode_solve(y_0=R0+Rpunto0+omega0+Q0,t_span=[0,20], num_points=100)

#gráfico velocidad angular omega
a=T.solution
p=[[soln[1][6],soln[1][7],soln[1][8]] for soln in a]
B=list_plot(p)
show(B)

3D rendering not yet implemented

---

La ecuación de Euler

---

Consideremos el caso en que la fuerza total y el torque desaparecen, esto es, $\textbf{F}=0$ y $\textbf{T}=0$. De $\textbf{T}=0$ se tiene que $\sum_{i=1}^N Q\textbf{u}_i\times F_i=0$ (ver ecuación (06)). Entonces las ecuaciones de movimiento de un cuerpo rígido (07)-(09) se reducen a $$\begin{equation} M\ddot{\textbf{R}}=0, (19) \end{equation}$$ $$\begin{equation} A\dot{\omega}+\omega\times A\omega=0, (20) \end{equation}$$ $$\begin{equation} \dot{Q}=Q\Omega, (21) \end{equation}$$

Como era de esperar la primera ecuación da $\textbf{R} (t) = \dot{\textbf{R}} (0) t + \textbf{R} (0)$, por lo que el centro de masa se mueve en una línea recta con velocidad constante. La ecuación (20) se conoce como _ecuación de Euler_, su solución esta completa y la interpretación se da como sigue.

La ecuación de Euler es de primer orden, es un sistema no lineal de tres ecuaciones con tres incógnitas, a saber, las componentes de $\omega$. Como de costumbre, la primera parte del análisis de un sistema de ecuaciones diferenciales consiste en encontrar y clasificar puntos fijos.

**Puntos Fijos:** Los puntos fijos de la ecuación de Euler son vectores constantes de velocidad angular $\omega$ y deben satisfacer la ecuación algebraica:

$$\omega\times A\omega = 0.$$

Esto es equivalente a decir que $\omega$ y $A\omega$ se encuentran en la misma línea, es decir, que $\omega$ es un autovector del operador inercia $A$:

$$A\omega = \lambda\omega,$$

donde $\lambda = I_1, I_2$ o $I_3$. Por lo tanto, si los momentos principales de inercia $I_i$, $i = 1, 2, 3$ son todos distintos, entonces los puntos fijos de la ecuación de Euler son precisamente los puntos en los tres ejes principales: $\mathrm{gen}\lbrace e_i\rbrace, i = 1, 2, 3$. Si dos de los momentos de inercia son los mismos, por ejemplo $I_1 = I_2$, entonces el espacio característico correspondiente $E_{I_1} = \mathrm{gen} \lbrace e_1, e_2\rbrace$ constituye un plano de puntos fijos y el eje principal restante $E_{I_3} = gen\lbrace e_3\rbrace$ es una línea de puntos fijos perpendicular a este plano. Del mismo modo, si $I_2 = I_3$. Si todos los momentos de inercia son iguales, $I_1 = I_2 = I_3$, entonces el cuerpo se llama _perfectamente simétrico_. En este caso todos los puntos de $\mathbb{R}^3$ son puntos fijos de la ecuación de Euler y de este modo los puntos fijos son los únicos tipos de soluciones.

Vamos a determinar la estabilidad de los puntos fijos más adelante, pero primero debemos examinar qué tipo de movimiento de un cuerpo rígido corresponde a un punto fijo de la ecuación de Euler. Cuando $\omega = (\omega_1, \omega_2, \omega_3)$ es un punto fijo de las ecuaciones de Euler, el correspondiente movimiento del cuerpo rígido es particularmente simple.

Para ver esto, observemos que al ser $\omega$ constante, la matriz de velocidad angular $$\begin{equation} \Omega=\begin{bmatrix} 0&-\omega_3&\omega_2\\ \omega_3&0&-\omega_1\\ -\omega_2&\omega_1&0 \end{bmatrix}, (22) \end{equation}$$

es constante también. Por lo tanto la solución de $\dot{Q} = Q\Omega$ es simplemente $$Q(t)=e^{\Omega t}$$

Utilizando las técnicas sobre los sistemas lineales, podemos demostrar que, para cada $t$, $Q (t)$ es una rotación alrededor del eje a través de $\omega$ y que la rapidez de rotación (rapidez angular) es $\vert\omega\vert$. Los detalles de esto son los siguientes y, la discusión es válida para cualquier matriz antisimétrica $\Omega$, $3\times 3$.

En primer lugar, determinamos la forma canónica de Jordan para $\Omega$ y la utilizamos para calcular $e^{\Omega t}$.

Con $\Omega$ expresada en la forma (22), es fácil calcular $$\mathrm{det}(\Omega - \lambda I)=\lambda(\lambda^2 \vert \omega\vert^2),$$ y por lo tanto los autovalores de $\Omega$ son $\lambda= 0$, $\pm\vert\omega\vert i$. Además, puesto que $\Omega\omega=\omega\times\omega=0$, se deduce que $\omega$ es un autovector correspondiente al autovalor $0$. Uno puede verificar que un autovector complejo $v$ correspondiente al autovalor $\vert\omega\vert i$ es $$v=\begin{bmatrix} \omega_1\omega_3+\vert\omega\vert\omega_2 i \\ \omega_2\omega_3-\vert\omega\vert\omega_11 i \\ -(\omega_1+\omega_2)^2 \end{bmatrix}$$

Escribimos a $v$ en término de su parte real e imaginaria $v=u+qi$ y llamando $P = (q, u, \omega)$ a la matriz formada a partir de colocar a $q$, $u$ y $\omega$ como columnas, obtenemos de la teoría general (o cálculo directo) que $$P^{-1}\Omega P=J=\begin{bmatrix} 0&-\vert \omega\vert &0 \\ \vert \omega\vert & 0& 0\\ 0&0&0 \end{bmatrix}$$ que es la forma canónica de Jordan para $\Omega$. Entonces, para punto $c\in \mathbb{R}^3$, sea $b=P^{-1}c$. Tenemos $$c=Pb=b_1q+b_2u+b_3\omega,$$ donde $b=(b_1,b_2,b_3)$. Por lo tanto, $$\begin{darray}{rcl} Q(t)c &=& e^{\Omega t}c= e^{\Omega t} Pb=Pe^{J t}P^{-1} Pb\nonumber \\ &=& (b_1q+b_2u)\cos (\vert\omega\vert t)+(-b_2q+b_1u) \mathrm{sen} (\vert\omega\vert t)+b_3\omega, (23) \end{darray}$$

El movimiento de $c$ descrito por la ecuación (23) es un movimiento circular sobre el eje: $\mathrm{gen}\lbrace \omega\rbrace$, con rapidez angular constante.

---

La solución general de la ecuación de Euler

---

La solución general de la ecuación de Euler es más fácil de describir si escribimos esta ecuación en términos de los ejes de los momentos principales de inercia.

**Nota:** Utilizaremos $\omega_1,\omega_2,\omega_3$ para denotar las componentes de $\omega$ relativas a la base propia para $A$ $$\begin{equation} \omega(t)=\omega_1(t)e_1+\omega_2(t)e_2+\omega_3(t)e_3. (24) \end{equation}$$ Esperemos que esto no cause confusión con lo que hicimos anteriormente, donde los $\omega_i'\mathrm{s}$ eran las componentes de $\omega$ con respecto a la base estándar para $\mathbb{R}^3$ (y eran constantes). Ahora de $Ae_i = I_ie_i$, $i = 1, 2, 3$, se deduce que $$\begin{equation} A\omega=I_1\omega_1e_1+I_2\omega_2e_2+I_3\omega_3e_3, (25) \end{equation}$$ y $$\begin{equation} A\dot{\omega}=I_1\dot{\omega}_1e_1+I_2\dot{\omega}_2e_2+I_3\dot{\omega}_3e_3, (26) \end{equation}$$

Haciendo uso de esto y de un cálculo directo de $\omega \times A\omega$ en términos de los $e_i'\mathrm{s}$, se puede calcular fácilmente la ecuación de Euler: $$A\dot{\omega}+\omega\times A\omega=0,$$ en término de las expresiones anteriores y obtener la siguiente:

Ecuación de Euler en los ejes principales Coordenados: $$\begin{equation} \dot{\omega}_1=\left( \dfrac{I_2-I_3}{I_1}\right) \omega_2\omega_3, (27) \end{equation}$$ $$\begin{equation} \dot{\omega}_2=\left( \dfrac{I_3-I_1}{I_2}\right) \omega_3\omega_1, (28) \end{equation}$$ $$\begin{equation} \dot{\omega}_3=\left( \dfrac{I_1-I_2}{I_3}\right) \omega_1\omega_2, (29) \end{equation}$$

Tengamos en cuenta que, en cada ecuación, la expresión en los paréntesis grandes, es constante ya que los momentos principales de inercia $I_1, I_2, I_3$ son constantes.

Este sistema tiene una forma particularmente agradable y, de hecho, uno puede deducir fácilmente dos leyes de conservación que se siguen directamente de ella. Estas dos leyes determinan completamente las soluciones del sistema, como veremos más adelante.

Supongamos, entonces, que $\omega :I \rightarrow\mathbb{R}^3$ es una solución (curva integral) de las ecuaciones de Euler$'$s. Una ley de conservación se obtiene simplemente multiplicando la $i$-ésima ecuación en el sistema por $I_i\omega_i$ (para $i = 1, 2, 3$) y luego sumamos las tres ecuaciones juntas y tenemos, $$I_1\omega_1\dot{\omega}_1+I_2\omega_2\dot{\omega}_2+I_3\omega_3\dot{\omega}_3=(I_2-I_3)\omega_1\omega_2\omega_3+(I_3-I_1)\omega_1\omega_2\omega_3+(I_1-I_2)\omega_1\omega_2\omega_3 = 0.$$

Por lo tanto, existe una constante $k$, que depende de la curva integral $\omega$, tal que:

Ley de Conservación I (energía cinética rotacional): $$\begin{equation} I_1\omega_1^2+I_2\omega_2^2+I_3\omega_3^2=2k, (30) \end{equation}$$ para todo $t \in I$.

Una segunda ley de conservación resulta, algo menos transparente, de multiplicar la $i$-ésima ecuación en el sistema por $I_2^2\omega_i$ (para $i = 1, 2, 3$) y luego sumar las tres ecuaciones juntas. Ahora tenemos la expresión $$I_1^2\omega_1\dot{\omega}_1+I_2^2\omega_2\dot{\omega}_2+I_3^2\omega_3\dot{\omega}_3=[(I_2-I_3)I_1+(I_3-I_1)I_2+(I_1-I_2)I_3+\omega_1\omega_2\omega_3 = 0.$$ Por lo tanto, existe una constante $a> 0$, que depende de la curva integral $\omega$, tal que:

Ley de Conservación II (momento angular): $$\begin{equation} I_1^2\omega_1^2+I_2^2\omega_2^2+I_3^2\omega_3^2=a^2, (31) \end{equation}$$ para todo $t \in I$.

Es fácil ver a partir de las ecuaciones (24) - (25) que si $\dot{\textbf{R}}(0) = 0$, entonces la constante $k = A\omega\cdot\omega$ es la energía cinética de rotación del sistema. Si además asumimos sin pérdida de generalidad que $\textbf{R} (0) = 0$, entonces la constante $a = \vert A\omega \vert = \vert\textbf{L}\vert$ es la _magnitud del momento angular_ del sistema. Así, las dos ecuaciones (30) - (31) son las leyes para la _conservación de la energía cinética rotacional y el momento angular_ del cuerpo rígido.

Las ecuaciones (30) - (31) pueden ser consideradas como las ecuaciones para dos elipsoides en $\mathbb{R}^3$ con respecto a los ejes principales $e_1, e_2, e_3$. A pesar de que en las ecuaciones $\omega_1, \omega_2, \omega_3$ son funciones dependientes del tiempo, cuando no hay peligro de confusión, usamos $\omega_1, \omega_2, \omega_3$ para denotar las coordenadas de puntos en $\mathbb{R}^3$. Debemos tener en cuenta que los tres ejes de cada elipsoide coinciden con estos ejes principales de inercia. Estos dos elipsoides son llamados _elipsoide de energía cinética_ y _elipsoide de momento angular_, respectivamente. Históricamente, el elipsoide de energía cinética es también llamado _elipsoide Poinsot_.

El contenido de las dos leyes de conservación es que cualquier curva integral $\omega: I\rightarrow\mathbb{R}^3$ de la ecuación de Euler debe yacer sobre la superficie de cada elipsoide y por lo tanto en la curva de intersección de estas dos superficies. Por lo tanto, las curvas de intersección son, con excepción de la dependencia explícita del tiempo, las mismas que las curvas integrales.

La visualización de las curvas de intersección de estos elipsoides, y por lo tanto las curvas integrales de la ecuación de Euler, es difícil de hacer a mano. Sin embargo, hay varias maneras de reescribir las leyes de conservación que conducen a una visualización más fácil. Una manera es multiplicar la ecuación (30) por $\frac{a^2}{2k}$ y luego restar la ecuación (31) de él para conseguir$$\begin{equation} I_1\left( \dfrac{a^2}{2k}-I_1\right)\omega_1^2+I_2\left( \dfrac{a^2}{2k}-I_2\right)\omega_2^2+I_3\left( \dfrac{a^2}{2k}-I_3\right)\omega_3^2=0. (32) \end{equation}$$

Esta ecuación describe un cono doble en $\mathbb{R}^3$ en función de las constantes $I_1, I_2, I_3, k, a$. El razonamiento que conduce a esto nos da que cada curva integral se encuentra en este cono, así como en los dos elipsoides. La intersección de dos de estas superficies es suficiente para determinar las curvas integrales. Los ejes de los elipsoides coinciden con los ejes principales de inercia y el eje más largo, de longitud $\frac{2k}{I_1}$ corresponde al momento de inercia más pequeño. Fijamos este elipsoide y consideramos cómo el cono (32), para varios valores de $a$, se cruza con este elipsoide.

Si denotamos por $b_i = I_i (\frac{a^2}{2k} - I_i)$, para $i = 1, 2, 3$, entonces cada uno de ellos es negativo cuando $\frac{a^2}{2k} [removed] I_3$. Por lo tanto la ecuación para el cono $$b_1\omega_1^2+b_2\omega_2^2+b_3\omega_3^2=0$$ no tiene solución en cualquiera de estos casos, y así no hay cono. Otra forma de decir esto es la siguiente: las magnitudes del momento angular y la energía cinética rotacional, $a$, $k$ de cualquier movimiento correspondiente a una solución de la ecuación de Euler debe satisfacer $$\begin{equation} I_1 \leq \frac{a^2}{2k} \leq I_3. (33) \end{equation}$$

#Vamos a graficar la intersección del elipsoide de energía cinética y el cono.

#Datos
a=1
b=2
c=3
g=9.8
M=10;
L=srange(0,1.5,.5)
r=[ [a*i,b*j,c*k] for i in L for j in L for k in L] #posiciones de las partículas respecto a un origen.
N=len(r) #cantidad de partículas
F=vector([0,0,0]) # La fuerza es cero.
CentroMasa=sum([vector(h)/N for h in r ]) #calculamos el centro de masa.
U=[vector(h)-CentroMasa for h in r] #posiciones relativas al centro de masa.

#matriz inercia
A11 = sum([u.norm()^2-u[0]^2 for u in U])*M/N
A12 = -sum([u[0]*u[1] for u in U])*M/N
A13 = -sum([u[0]*u[2] for u in U])*M/N
A21 = -sum([u[1]*u[0] for u in U])*M/N
A22 = sum([u.norm()^2-u[1]^2 for u in U])*M/N
A23 = -sum([u[1]*u[2] for u in U])*M/N
A31 = -sum([u[2]*u[0] for u in U])*M/N
A32 = -sum([u[2]*u[1] for u in U])*M/N
A33 = sum([u.norm()^2-u[2]^2 for u in U])*M/N
A = matrix([[A11,A12,A13],[A21,A22,A23],[A31,A32,A33]])
Ainversa = A.inverse()
autovalores = A.eigenvalues()
I1 = vector(autovalores)
I1

(21.6666666666667, 16.6666666666667, 8.33333333333333)

#elipsoide de energía cinética
x, y, z = var('x,y,z')
k=2
elipsoide = I1[2]*x^2 + I1[1]*y^2 + I1[0]*z^2 - 2*k
E=implicit_plot3d(elipsoide, (x, -(2*k/I1[2])^(1/2), (2*k/I1[2])^(1/2)), (y, -(2*k/I1[1])^(1/2), (2*k/I1[1])^(1/2)), (z, -(2*k/I1[0])^(1/2), (2*k/I1[0])^(1/2)),opacity=.2,color="blue")

#cono doble
x, y, z = var('x,y,z')
a=4
b1 = I1[2]*(a^2/2*k - I1[2])
b2 = I1[1]*(a^2/2*k - I1[1])
b3 = I1[0]*(a^2/2*k - I1[0])
cono = b1*x^2 + b2*y^2 + b3*z^2
C=implicit_plot3d(cono, (x, -(2*k/I1[2])^(1/2), (2*k/I1[2])^(1/2)), (y, -(2*k/I1[1])^(1/2), (2*k/I1[1])^(1/2)), (z, -(2*k/I1[0])^(1/2), (2*k/I1[0])^(1/2)),opacity=.3,color="red")

show(E+C)

3D rendering not yet implemented

Así, con $k$ fijo, restringimos de modo que $\frac{a^2}{2k}$ se encuentre en el intervalo $[I_1, I_3]$ y miramos como varía la familia de conos resultante $C = C_a$ durante ese intervalo. Dividimos en casos, ya que la descripción es diferente dependiendo de si los momentos de inercia son todos distintos o no.

(1) **Tres momentos de inercia distintos:** Supongamos que $0 <I_1 <I_2 [removed] 0$, $b_3> 0$ y, por supuesto, esto significa que $\omega_2 = 0$, $\omega_3 = 0$ y el "cono" $C$ se reduce a una línea a lo largo de $e_1$. Por lo tanto, $C$ es el primer eje principal. A continuación $C \cap E=\lbrace(\pm \sqrt{\frac{2k}{I_1}}, 0, 0)\rbrace$ es un par de puntos, que son puntos fijos de las ecuaciones de Euler.

A continuación, para $a$ tal que $I_1 <\frac{a^2}{2k} [removed] 0$, mientras que $b_2 <0$, $b_3 <0$, y por lo tanto la ecuación para $C$ puede ser escrita como $$\omega_1^2=\frac{-b_2}{b_1}\omega_2^2+\frac{-b_3}{b_1}\omega_3^2.$$ Por lo tanto, $C$ es un cono elíptico con su eje coincidiendo con el primer eje principal. La intersección $C \cap E$ es bastante fácil de visualizar ya que consiste en un par de curvas circulares centradas sobre el primer eje principal.

Estos pares de curvas incrementan su tamaño a medida que $a$ aumenta, pero luego, cuando $a$ alcanza el valor para el cual $\frac{a^2}{2k} = I_2$, $C \cap E$ cambia naturalmente. Para este valor de $a$, los coeficientes $b_1>0$, $b_2 = 0$, $b_3 <0$ y por lo que la ecuación para $C$ puede escribirse como $$\omega_1=\pm \sqrt{\frac{-b_3}{b_1}}\omega_3.$$

Por lo tanto, $C$ es un par de planos, cada uno de los cuales es perpendicular al plano $\lbrace e_1, e_3\rbrace$. La intersección de este par de planos con el elipsoide $E$ es fácil de visualizar y da un par de curvas cruzadas. El par realmente se compone de cuatro curvas integrales, cada uno de los cuales va desde uno de los puntos fijos en el eje $e_2$ al otro punto fijo en este eje.

Para un valor ligeramente mayor de $a$, de modo que $I_2 <\frac{a^2}{2k}[removed] 0$, $b_2> 0$, mientras que $b_3 <0$. Así la ecuación para el cono se puede escribir como $$\omega_3^2=\frac{b_1}{b_3}\omega_1^2+\frac{b_2}{b_3}\omega_2^2.$$

Por lo tanto, el cono $C$ es un cono elíptico con su eje a lo largo del tercer eje director. Como antes, la intersección de $C \cap E$ es un par de curvas cerradas circulares, pero ahora centrado alrededor del tercer eje principal.

Finalmente, cuando $a$ es lo más grande posible, $\frac{a^2}{2k} = I_3$ y el cono $C$ degenera en una línea recta correspondiente al tercer eje principal. Entonces $C \cap E$ es un par de puntos fijos para la ecuación de Euler.

(2)**Dos momentos distintos de inercia:** Supongamos que $I_1 = I_2 <I_3$. La otra posibilidad: $I_1 <I_2 = I_3$, tiene un análisis similar. Primero notemos que el elipsoide es una superficie de revolución obtenida, por ejemplo, girando una elipse en el plano $e_2-e_3$ alrededor del eje $e_3$. El análisis de cómo el cono $C_a$ intersecta este elipsoide es como sigue.

Para $a$ mínimo, tenemos $\frac{a^2}{2k}= I_1 = I_2$, por lo que el "cono" $C$ se degenera en plano $w_3 = 0$. De hecho, para este valor de $a$, $$C = \mathrm{gen}\lbrace e_1, e_2\rbrace = E_{I_1}$$ es el espacio característico bidimensional correspondiente al autovalor $I_1 = I_2$. Por lo tanto la intersección $C \cap E$ es un círculo y cada punto de este círculo es un punto fijo de la ecuación de Euler.

#momentos de inercia
k=1/2
I1=3
I2=3
I3=9
a=(2*k*I1)^(1/2)

#elipsoide de energía cinética
x, y, z = var('x,y,z')
elipsoide = I1*x^2 + I2*y^2 + I3*z^2 - 2*k
E=implicit_plot3d(elipsoide, (x, -(2*k/I1)^(1/2), (2*k/I1)^(1/2)), (y, -(2*k/I2)^(1/2), (2*k/I2)^(1/2)), (z, -(2*k/I3)^(1/2), (2*k/I3)^(1/2)),opacity=.2,color="blue")

#cono doble
x, y, z = var('x,y,z')
b1 = I1*(a^2/(2*k) - I1)
b2 = I2*(a^2/(2*k) - I2)
b3 = I3*(a^2/(2*k) - I3)
#cono = b1*x^2 + b2*y^2 + b3*z^2 reducimos la ecuación a z=0, ya que sage no grafica el plano cuando z^2=0
cono= z==0
C=implicit_plot3d(cono, (x, -(2*k/I1)^(1/2), (2*k/I1)^(1/2)), (y, -(2*k/I2)^(1/2), (2*k/I2)^(1/2)), (z, -(2*k/I3)^(1/2), (2*k/I3)^(1/2)),opacity=.3,color="red")

show(E+C)

3D rendering not yet implemented

Cuando $a$ tiene un valor un poco más grande de modo que $I_2 <\frac{a^2}{2k} <I_3$, el cono $C$ es un cono circular cuyo eje coincide con el tercer eje principal. Por lo tanto, la intersección $C \cap E$ es un par de curvas circulares cerradas centradas en el tercer eje principal.

#momentos de inercia
k=1/2
a=2
I1=3
I2=3
I3=9
a^2/(2*k)

4

#elipsoide de energía cinética
x, y, z = var('x,y,z')
elipsoide = I1*x^2 + I2*y^2 + I3*z^2 - 2*k
E=implicit_plot3d(elipsoide, (x, -(2*k/I1)^(1/2), (2*k/I1)^(1/2)), (y, -(2*k/I2)^(1/2), (2*k/I2)^(1/2)), (z, -(2*k/I3)^(1/2), (2*k/I3)^(1/2)),opacity=.2,color="blue")

#cono doble
x, y, z = var('x,y,z')
b1 = I1*(a^2/(2*k) - I1)
b2 = I2*(a^2/(2*k) - I2)
b3 = I3*(a^2/(2*k) - I3)
cono = (b1*x^2 + b2*y^2 + b3*z^2==0).simplify()

C=implicit_plot3d(cono, (x, -(2*k/I1)^(1/2), (2*k/I1)^(1/2)), (y, -(2*k/I2)^(1/2), (2*k/I2)^(1/2)), (z, -(2*k/I3)^(1/2), (2*k/I3)^(1/2)),opacity=.3,color="red")

show(E+C)

3D rendering not yet implemented

simplify(z^2==0)

z^2 == 0

-54*z == 0

Finalmente, cuando $a$ es tal que $\frac{a^2}{2k} = I_3$, el cono $C_a$ se degenera en una línea recta que coincide con el tercer eje principal. De ahí $C\cap E$ es un par de puntos y estos puntos son puntos fijos de la ecuación de Euler.

Cada punto fijo desde el plano $E_{I_1} = \mathrm{gen} \lbrace e_1, e_2\rbrace$ de puntos fijos es inestable, mientras que cada punto fijo de la línea $E_{I_1} = \mathrm{gen}\lbrace e_3\rbrace$ de puntos fijos es estable.

(3) **Sólo un momento de inercia:** En el caso en el que los tres momentos de inercia coinciden: $I_1 = I_2 = I_3$, está claro que sólo hay un valor de $a$ para el cual $C$ no está vacío, es decir, para $a =\sqrt{2kI_3}$, y que para este valor, $C = \mathbb{R}^3$. En este caso todo el espacio $\mathbb{R}^3$ consta de puntos fijos (estables) de la ecuación de Euler. Entonces también, $C\cap E = E$ es una esfera de puntos fijos. Esto es de esperar, ya que, como hemos mencionado, el cuerpo rígido se llama perfectamente simétrico cuando todos sus momentos principales de inercia son los mismos. Los momentos de inercia alrededor de cualquier eje son los mismos y los únicos movimientos posibles del cuerpo (con fuerza total y torque igual a cero) son giros alrededor de dichos ejes con velocidad angular uniforme.

#momentos de inercia
I1=1
I2=1
I3=1

#elipsoide de energía cinética
x, y, z = var('x,y,z')
k=1/2
elipsoide = I1*x^2 + I2*y^2 + I3*z^2 - 2*k
E=implicit_plot3d(elipsoide, (x, -(2*k/I1)^(1/2), (2*k/I1)^(1/2)), (y, -(2*k/I2)^(1/2), (2*k/I2)^(1/2)), (z, -(2*k/I3)^(1/2), (2*k/I3)^(1/2)),opacity=.2,color="blue")

#cono doble
x, y, z = var('x,y,z')
a=(2*k*I3)^(1/2)
b1 = I1*(a^2/(2*k) - I1)
b2 = I2*(a^2/(2*k) - I2)
b3 = I3*(a^2/(2*k) - I3)
cono = b1*x^2 + b2*y^2 + b3*z^2
C=implicit_plot3d(cono, (x, -(2*k/I1)^(1/2), (2*k/I1)^(1/2)), (y, -(2*k/I2)^(1/2), (2*k/I2)^(1/2)), (z, -(2*k/I3)^(1/2), (2*k/I3)^(1/2)),opacity=.3,color="red")

show(E+C)

3D rendering not yet implemented

El siguiente algoritmo aproxima las soluciones de la ecuación de Euler, con tres momentos de inercia distintos.

Se comprueba el contenido de las dos leyes de conservación. El mismo es, como ya mencionamos, que cualquier curva integral $\omega: I\rightarrow\mathbb{R}^3$ de la ecuación de Euler debe yacer sobre la superficie de cada elipsoide y por lo tanto en la curva de intersección de estas dos superficies. Por lo tanto, las curvas de intersección son, con excepción de la dependencia explícita del tiempo, las mismas que las curvas integrales.

#Datos conocidos
a=1
b=2
c=3
g=9.8
M=10;
L=srange(0,1.5,.5)
r=[ [a*i,b*j,c*k] for i in L for j in L for k in L] #posiciones de las partículas respecto a un origen.
N=len(r) #cantidad de partículas
F=vector([0,0,0]) # La fuerza es cero.
CentroMasa=sum([vector(h)/N for h in r ]) #calculamos el centro de masa.
U=[vector(h)-CentroMasa for h in r] #posiciones relativas al centro de masa.

#matriz inercia
A11 = sum([u.norm()^2-u[0]^2 for u in U])*M/N
A12 = -sum([u[0]*u[1] for u in U])*M/N
A13 = -sum([u[0]*u[2] for u in U])*M/N
A21 = -sum([u[1]*u[0] for u in U])*M/N
A22 = sum([u.norm()^2-u[1]^2 for u in U])*M/N
A23 = -sum([u[1]*u[2] for u in U])*M/N
A31 = -sum([u[2]*u[0] for u in U])*M/N
A32 = -sum([u[2]*u[1] for u in U])*M/N
A33 = sum([u.norm()^2-u[2]^2 for u in U])*M/N
A = matrix([[A11,A12,A13],[A21,A22,A23],[A31,A32,A33]])
Ainversa = A.inverse()
autovalores = A.eigenvalues()
I1 = vector(autovalores)
I1

(21.6666666666667, 16.6666666666667, 8.33333333333333)

T =  ode_solver()
def CuerpoRigido(t,y):
    omega=vector([y[6],y[7],y[8]]) #velocidad angular
    Omega=Matrix([[0,-y[8],y[7]],[y[8],0,-y[6]],[-y[7],y[6],0]]) #operador velocidad angular
    Q=Matrix([[y[9],y[10],y[11]],[y[12],y[13],y[14]],[y[15],y[16],y[17]]]) #matriz de rotación 
    Qpunto=Q*Omega
    Suma=vector([0,0,0]) #condición de la ecuación de Euler.
    wpunto=Ainversa*(- omega.cross_product(A*omega)+Suma)
    return [y[3],y[4],y[5],0,0,0,wpunto[0],wpunto[1],wpunto[2],Qpunto[0,0],Qpunto[0,1],Qpunto[0,2],Qpunto[1,0],Qpunto[1,1],Qpunto[1,2],Qpunto[2,0],Qpunto[2,1],Qpunto[2,2]]
T.function=CuerpoRigido
T.algorithm="rk4"

#elipsoide de energía cinética rotacional
x, y, z = var('x,y,z')
k=1/2
elipsoide = I1[2]*x^2 + I1[1]*y^2 + I1[0]*z^2 - 2*k
E=implicit_plot3d(elipsoide, (x, -(2*k/I1[2])^(1/2), (2*k/I1[2])^(1/2)), (y, -(2*k/I1[1])^(1/2), (2*k/I1[1])^(1/2)), (z, -(2*k/I1[0])^(1/2), (2*k/I1[0])^(1/2)),opacity=.2,color="red")
show(E)

3D rendering not yet implemented

R0=[0,0,0] #condición inicial
Rpunto0=[0,0,0] #condición inicial
Q0=[1,0,0,0,1,0,0,0,1] #condición inicial
x1=(2*k/I1[2])^(1/2)
Dx=(x1+x1)/10
for om1 in srange(-x1+Dx,x1,Dx):
    x2=((2*k-I1[2]*om1^2)/I1[1])^(1/2)
    Dx2=(x2+x2)/10
    for om2 in srange(-x2+Dx2,x2,Dx2):
        om3=((2*k-I1[2]*om1^2-I1[1]*om2^2)/I1[0])^(1/2)
        omega0=[om1,om2,om3] #condición inicial
        T.ode_solve(y_0=R0+Rpunto0+omega0+Q0,t_span=[0,40], num_points=100)
        a=T.solution
        p=[[soln[1][6],soln[1][7],soln[1][8]] for soln in a]
        E+=list_plot(p)

show(E)

3D rendering not yet implemented