SharedMovimiento de un cuerpo rígido.sagewsOpen in CoCalc
Download, Raw, Embed
Author: Fernando Mazzone
Views : 8
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 NN partículas y que Fi:J×O→R3,  i=1,...,NF_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 R3N\mathbb{R}^{3N} r=(r1,...,rN):I→R3N\textbf{r}=(\textbf{r}_1,...,\textbf{r}_N):I\rightarrow\mathbb{R}^{3N} se dice que es un movimiento de cuerpo rígido si para todo i,j=1,...,Ni,j=1,...,N y para todo t∈It\in I.
En general, sólo nos ocuparemos de movimientos de cuerpos rígidos que satisfagan las ecuaciones movimiento mir¨i=Fi(t,r,r˙)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 ri\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 r=(r1,...,rN)\textbf{r} = (\textbf{r}_1,..., \textbf{r}_N) tiene una forma muy particular, como muestra el siguiente teorema importante.
Teorema:Sea r=(r1,...,rN):I→R3N\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⊂R→SO(3)Q:I\subset \mathbb{R}\rightarrow SO(3), con SO(3)={O∈R3×3:OTO=I;  det(O)=1}SO(3)=\lbrace O\in \mathbb{R}^{3\times 3} : O^T O=I; \; det(O)=1\rbrace, as'i QQ es una matriz ortogonal: para cada t∈It\in I con Q(0)=IQ(0)=I y NN-vectores ui∈R3\textbf{u}_i\in \mathbb{R}^{3} i=1,...,Ni=1,...,N que satisfacen ∑i=1Nmiui=0\sum_{i=1}^{N}m_i\textbf{u}_i=0, tal que para cada i=1,...,Ni=1,...,N y para todo t∈It\in I. Donde R(t)\textbf{R}(t) es la posición del centro de masa en el tiempo tt.
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 para i=1,...,Ni = 1,. . . , N, que expresan los vectores de posición de las partículas en términos del vector posición R=R(t)\textbf{R} = \textbf{R} (t) del centro de masa, la matriz de rotación Q=Q(t)Q = Q (t), y los vectores ui\textbf{u}_i (constantes), que son las posiciones iniciales de los cuerpos en relación con el centro de la masa inicial (ri(0)=ui+R(0)\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 r:I→R3N\textbf{r}: I \rightarrow\mathbb{R}^{3N} de NN partículas está completamente determinado por el movimiento de su centro de masa R:I→R3\textbf{R}: I \rightarrow\mathbb{R}^3 y una cierta familia {Q(t)}t∈I\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 R\textbf{R} y QQ 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 QQ.
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 Fi=Fi(t,r1,...,rN,r˙1,...,r˙N)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 R\textbf{R}, QQ y sus derivadas usando ri=Qui+R\textbf{r}_i=Q\textbf{u}_i+\textbf{R} y en consecuencia
Habiendo hecho esto, la fuerza total F\textbf{F} y el torque T\textbf{T} se convierten en funciones de R\textbf{R}, QQ y sus derivadas. Entonces, ya que la posición R\textbf{R} del centro de masa satisface MR¨=F,M\ddot{\textbf{R}}=\textbf{F}, podemos interpretar esta ecuación como una parte del sistema de ecuaciones diferenciales para R\textbf{R} y QQ que estamos buscando.
Las ecuaciones diferenciales para QQ involucran su derivada Q˙\dot{Q} y como {Q(t)}t∈I\lbrace Q (t)\rbrace_{t\in I} es una familia 1-paramétrica de rotaciones, está compuesta por una familia 1-paramétrica {Ω(t)}t∈I\lbrace \Omega(t)\rbrace_{t\in I} de matrices antisimétricas.
Proposición: Sea Q:J→O(n)Q:J\rightarrow O(n) una familia 1-paremétrica (suave) de matrices ortogonales n×nn\times n. Entonces, existe una familia 1-paramétrica (suave) Ω:J→o(n)\Omega:J\rightarrow o(n) de matrices antisimétricas n×nn\times n tal que QQ satisface la ecuación diferencial Q˙=QΩ.\dot{Q}=Q\Omega. Además, si Q(0)=IQ(0)=I, entonces det(Q(t))=1det(Q(t))=1, es decir, Q(t)Q(t) es una matriz de rotación, para cada t∈Jt\in J.

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

Observamos que el resultado anterior nos permite considerar a QQ 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: ω1,ω2,ω3:J→R\omega_1, \omega_2, \omega_3: J \rightarrow \mathbb{R}. Es decir, Ω\Omega tiene la forma
Hay una razón para expresar Ω\Omega en términos de ω1,ω2,ω3\omega_1, \omega_2, \omega_3 de esta manera muy particular. Es decir, si entonces para cualquier vector v=(v1,v2,v3)v = (v_1, v_2, v_3), tenemos Ωv=ω×v\Omega v = \omega \times v. Por lo tanto, la acción de la matriz antisimétrica Ω\Omega sobre los vectores vv es la misma que la acción del producto vectorial de ω\omega en vv. Esta relación entre Ω\Omega y ω\omega vale en general: cualquier matriz antisimétrica BB, 3×33\times 3 tiene su acción v→Bvv \rightarrow Bv representado por la acción del producto vectorial v→b×vv \rightarrow b \times v, con el vector fijo b=(−B23,B13,−B12)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 QQ.
Podemos escribir fácilmente la velocidad de la ii-ésima partícula en términos de la velocidad angular de la siguiente manera: como ri=Qui+R\textbf{r}_i=Q \textbf{u}_i+\textbf{R} tenemos
Esto expresa la velocidad de la ii-ésima partícula del cuerpo como la suma de una velocidad de traslación R\textbf{R} y una velocidad de rotación Q(ω×ui)=Qω×QuiQ (\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 R˙\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ωQ\omega.
Usando las expresiones anteriores, también podemos reescribir el momento angular total L\textbf{L} del sistema.

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

Aquí hemos utilizado el hecho de que las matrices de rotación se pueden distribuir en el producto vectorial: Q(v×w)=Qv×QwQ (v \times w) = Qv \times Qw. Si denotamos por Aω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:**
La notación AωA\omega introduce un operador importante AA en la teoría. Este operador es definido por lo siguiente:
**Definición:** _Operador de inercia AA._ Sean m1,...,mNm_1,...,m_N números positivos dados (las masas) y u1,...,uN\textbf{u}_1,...,\textbf{u}_N vectores dados (las posiciones iniciales de las partículas relativas al centro de masa), con ∑i=1Nmiui=0\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:R3→R3A:\mathbb{R}^3\rightarrow \mathbb{R}^3 definido por para v∈R3v\in \mathbb{R}^3. Este operador se conoce como operador de inercia o tensor de inercia para el cuerpo rígido.
El operador AA captura las características geométricas del cuerpo rígido que son importantes para la naturaleza de su movimiento. Observemos que AA no depende del tiempo, pero si lo hace de las posiciones iniciales u1,...,uN\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 L˙=T\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
Entonces calculando el torque total obtenemos
Sustituyendo estas expresiones en L˙=T\dot{\textbf{L}}=\textbf{T} y reduciendo cada lado usando que MR¨=FM \ddot{\textbf{R}}=\textbf{F} da Q(Aω˙+ω×Aω)=∑i=1NQui×Fi,Q(A\dot{\omega}+\omega\times A\omega)=\sum_{i=1}^N Q\textbf{u}_i\times F_i, i.e., Q(Aω˙+ω×Aω)=Q(∑i=1Nui×QTFi).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 QQ en el lado izquierdo de la ecuación se deduce de la identidad general: Qv×Qw=Q(v×w)Qv\times Qw=Q(v\times w). Canculando QQ 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:}

Las incógnitas aquí son R\textbf{R}, ω\omega y QQ (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 RiR_i, ωi\omega_i, QijQ_{ij}, i,j=1,2,3i,j=1,2,3. Las \emph{condiciones iniciales} son (en término de las incógnitas originales)
R(0)=∑i=0NmiMri(0)\textbf{R}(0)=\sum_{i=0}^N \frac{m_i}{M}\textbf{r}_i(0) R˙(0)=∑i=0NmiMr˙i(0)\dot{\textbf{R}}(0)=\sum_{i=0}^N \frac{m_i}{M}\dot{\textbf{r}}_i(0) Q(0)=IQ(0)=I Aω(0)=L(0)−MR(0)×R˙(0).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 AA es invertible. El caso en que AA 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 QQ 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 u1,...,uN\textbf{u}_1,...,\textbf{u}_N (con respecto al centro de masa inicial R(0)\textbf{R}(0)) y por las masas m1,...,mNm_1,. . ., m_N, tal que ∑i=1Nmiui=0\sum_{i=1}^N m_i\textbf{u}_i=0 y ∑i=1Nmi=M\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 AA y la siguiente proposición describe la naturaleza de AA y muestra cómo entra en la expresión para la energía cinética del sistema.
**Proposición:** Sea AA un operador de inercia definido por Av=∑i=1Nmiui×(v×ui),Av=\sum_{i=1}^N m_i\textbf{u}_i\times (v\times \textbf{u}_i), para v∈R3v\in \mathbb{R}^3. Entonces para todo v,w∈R3v,w\in \mathbb{R}^3. En particular para todo v∈R3v\in \mathbb{R}^3. Por lo tanto AA es una matriz simétrica, semidefinida positiva. Es diagonalizable y tiene tres autovalores I1I_1, I2I_2, I3I_3 no negativos: 0≤I1≤I2≤I30\leq I_1 \leq I_2\leq I_3 y los tres autovectores correspondientes Aei=Iiei,        i=1,2,3Ae_i=I_ie_i, \;\;\;\;i=1,2,3 pueden ser elegidos de modo que sean ortogonales ei⋅ej=0e_i\cdot e_j=0, i≠ji\neq j.
Si K=∑i=1Nmi∣r˙∣2K=\sum_{i=1}^N m_i\vert \dot{\textbf{r}}\vert ^2 es la _energía cinética total del sistema_, entonces 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 AA como I1≤I2≤I3I_1\leq I_2\leq I_3 y seleccionaremos los tres autovectores correspondientes e1,e2,e3e_1, e_2, e_3, que son ortogonales, de longitud uno, ∣ei∣=1\vert e_i \vert = 1, i=1,2,3i = 1, 2, 3, y tal que {e1,e2,e3}\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 I1=I2≠I3I_1 = I_2 \neq I_3, la selección de dos vectores ortonormales e1,e2e_1, e_2 desde el espacio característico EI1E_{I_{1}} puede hacerse en un número infinito de formas, pero habiendo hecho una elección, entonces e3e_3 se determina de forma única. Cuando todos los autovalores son los mismos el espacio característico EI1E_{I_{1}} está en R3\mathbb{R}^3 y así cualquier orientación positiva puede ser seleccionada en la base ortonormal {e1,e2,e3}\lbrace e_1, e_2, e_3\rbrace para R3\mathbb{R}^3.
Como es habitual, identificamos la base estándar de R3\mathbb{R}^3 por {ξ1,ξ2,ξ3}\lbrace \xi_1, \xi_2, \xi_3\rbrace donde ξ1=(1,0,0),    ξ2=(0,1,0),    ξ3=(0,0,1).\xi_1 = (1, 0, 0),\;\; \xi_2 = (0, 1, 0),\;\; \xi_3 = (0, 0, 1).
Por lo tanto, el operador de inercia A:R3→R3A: \mathbb{R}^3\rightarrow\mathbb{R}^3 se identifica con la matriz 3×33\times 3 A={Aij}i,j=1,2,3A = \lbrace A_{ij}\rbrace_{i, j = 1,2,3}, donde Aij≡Aξj⋅ξi.A_{ij} \equiv A\xi_j\cdot \xi_i.
Usando la definición de AA, es fácil obtener la siguiente representación explícita de AA.
**Proposición:** El operador inercia A={Aij}i,j=1,2,3A = \lbrace A_{ij}\rbrace_{i, j = 1,2,3} tiene entradas donde δij\delta_{ij} es la delta Kronecker (es decir, δij=0\delta_{ij} = 0, para i≠ji\neq j, δii=1\delta_{ii} = 1), y uk=(uk1,uk2,uk3),\textbf{u}_k = (u_{k1}, u_{k2}, u_{k3}), para k=1,...,Nk = 1,. . .,N. Alternativamente, sin el uso la delta de Kronecker esto puede escribirse como

Ejemplo:

Supongamos que hay tres partículas de igual masa, digamos mk=1m_k=1, k=1,2,3k=1,2,3, que componen el cuerpo rígido, y las posiciones iniciales de esas partículas son r1(0)=(0,0,0)\textbf{r}_1(0)=(0,0,0) r2(0)=(1,0,0)\textbf{r}_2(0)=(1,0,0) r3(0)=(0,1,0)\textbf{r}_3(0)=(0,1,0)
Para obtener el tensor de inercia, debemos primero calcular las posiciones relativas uk\textbf{u}_k, k=1,2,3k=1,2,3. Para estos datos el centro de masa inicial es R(0)=13(0,0,0)+13(1,0,0)+13(0,1,0)=(13,13,0).\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 u1=r1(0)−R(0)=(−13,−13,0)\textbf{u}_1=\textbf{r}_1(0)-\textbf{R}(0)=(-\frac{1}{3},-\frac{1}{3},0) u2=r2(0)−R(0)=(23,−13,0)\textbf{u}_2=\textbf{r}_2(0)-\textbf{R}(0)=(\frac{2}{3},-\frac{1}{3},0) u3=r3(0)−R(0)=(−13,23,0).\textbf{u}_3=\textbf{r}_3(0)-\textbf{R}(0)=(-\frac{1}{3},\frac{2}{3},0).
Para calcular el operador de inercia AA usamos las fórmulas (14) y (15). Tengamos en cuenta que la fórmula (14) dice que A11=∑k=1Nmk(uk22+uk32)A_{11} = \sum_{k=1}^N m_k (u_{k2}^2+u_{k3}^2) entonces en este ejemplo
A11=19+19+49=23A_{11}=\frac{1}{9}+\frac{1}{9}+\frac{4}{9}=\frac{2}{3}
Del mismo modo A22A_{22} implica la suma de los cuadrados de la primera y tercera componentes de uku_k, mientras A33A_{33} implica la suma de los cuadrados de la primera y segunda componenete. Entonces tenemos,
A22=19+49+19=23A_{22}=\frac{1}{9}+\frac{4}{9}+\frac{1}{9}=\frac{2}{3} A33=19+19+49+19+19+49=43.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 AijA_{ij} para i≠ji\neq j notemos por ejemplo que uk3=0u_{k3}=0 para todo kk y entonces A13=0=A23A_{13}=0=A_{23}
Para A12A_{12} tenemos A12=−[u11u12+u21u22+u31u32]=−[19−29−29]=13.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 AA es simétrica, estos son los único calculos que necesitamos. Por lo tanto A=[23130132300043].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 AA: I1=13,          e1=12(−1,1,0)I_1=\frac{1}{3},\;\;\;\;\;e_1=\frac{1}{\sqrt{2}}(-1,1,0) I2=1,          e2=12(1,1,0)I_2=1,\;\;\;\;\;e_2=\frac{1}{\sqrt{2}}(1,1,0) I3=43,          e3=12(0,0,1)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)
[(−13, −13, 0)\displaystyle \left(-\frac{1}{3},\,-\frac{1}{3},\,0\right), (23, −13, 0)\displaystyle \left(\frac{2}{3},\,-\frac{1}{3},\,0\right), (−13, 23, 0)\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)
(23130132300043)\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)
[1\displaystyle 1, 13\displaystyle \frac{1}{3}, 43\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)≡(Bv⋅v)(v⋅v)f (v) \equiv \frac{(Bv \cdot v)}{(v \cdot v)} cuyos valores son, respectivamente, autovalores máximos y mínimos de BB y estos valores extremos son asumidos en los respectivos autovectores vv de BB. En estos casos, la función ff se conoce como el _momento de inercia de la función_.
**Definición:**_Momento de Inercia._ La función I:R3−{0}→RI: \mathbb{R}^3-\lbrace 0\rbrace \rightarrow \mathbb{R} definida por I(v)=v⋅Avv⋅v,I(v)=\dfrac{v\cdot Av}{v\cdot v}, se llama _momento de inercia de la función_. El número I(v)I (v) se llama _momento de inercia del sistema_ sobre el eje a través del centro de masa en la dirección de vv. Debemos tener en cuenta que I(cv)=I(v)I(cv) = I (v) para todos los escalares no nulos cc y por eso nos da el mismo valor para todos los vectores no nulos en el eje determinado por vv. Los autovalores I1≤I2≤I3I_1\leq I_2 \leq I_3 de AA son llamados los _momentos principales de inercia_.
Se desprende de la teoría general, que los momentos principales de inercia son los valores de II en los autovectores de AA: Ii=I(ei),I_i = I (ei), i=1,2,3i = 1, 2, 3, y que I1≤I(v)≤I3I_1 \leq I (v) \leq I_3, para cada v≠0v \neq 0 en R3\mathbb{R}^3. Por lo tanto, I1I_1 es el momento de inercia _mínimo_ del sistema. Así mismo, I3I_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 e1e_1 da el menor momento de inercia I1=13I_1 = \frac{1}{3}, mientras que un giro alrededor del eje a través de e2e_2 da como resultado un mayor momento de inercia I2=1I_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 I3=43I_3 = \frac{4}{3} se produce para una revolución alrededor del eje a través de e3e_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)I (v) ayuda a clarificar la idea que se alude en el último párrafo. Esta expresión surge de observar que ∣v×ui∣\vert v \times \textbf{u}_i \vert es el área del paralelogramo determinado por vv y ui\textbf{u}_i. Esta zona es ∣v∣di(v)\vert v \vert d_i (v), donde di(v)d_i (v) denota la distancia desde la punta de ui\textbf{u}_i por vv.
Así, encontramos Av⋅v=∑i=1Nmi∣v×ui∣2=∣v∣2∑i=1Nmidi(v)2.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
Esta es la definición tradicional del momento de inercia del sistema alrededor del eje vv.
**Nota:** Aquí, uiu_i, i=1,...,Ni = 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 KK de la ecuación (12), en términos de II da la siguiente
**Energía cinética del cuerpo rígido:**
Por lo tanto, con respecto a la velocidad angular ω\omega del cuerpo, I(ω)I (\omega) desempeña el mismo papel que la masa con respecto a la velocidad lineal. Cuanto mayor es el valor de I(ω)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, F=0\textbf{F}=0 y T=0\textbf{T}=0. De T=0\textbf{T}=0 se tiene que ∑i=1NQui×Fi=0\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
Como era de esperar la primera ecuación da R(t)=R˙(0)t+R(0)\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:
ω×Aω=0.\omega\times A\omega = 0.
Esto es equivalente a decir que ω\omega y AωA\omega se encuentran en la misma línea, es decir, que ω\omega es un autovector del operador inercia AA:
Aω=λω,A\omega = \lambda\omega,
donde λ=I1,I2\lambda = I_1, I_2 o I3I_3. Por lo tanto, si los momentos principales de inercia IiI_i, i=1,2,3i = 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: gen{ei},i=1,2,3\mathrm{gen}\lbrace e_i\rbrace, i = 1, 2, 3. Si dos de los momentos de inercia son los mismos, por ejemplo I1=I2I_1 = I_2, entonces el espacio característico correspondiente EI1=gen{e1,e2}E_{I_1} = \mathrm{gen} \lbrace e_1, e_2\rbrace constituye un plano de puntos fijos y el eje principal restante EI3=gen{e3}E_{I_3} = gen\lbrace e_3\rbrace es una línea de puntos fijos perpendicular a este plano. Del mismo modo, si I2=I3I_2 = I_3. Si todos los momentos de inercia son iguales, I1=I2=I3I_1 = I_2 = I_3, entonces el cuerpo se llama _perfectamente simétrico_. En este caso todos los puntos de R3\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 ω=(ω1,ω2,ω3)\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
es constante también. Por lo tanto la solución de Q˙=QΩ\dot{Q} = Q\Omega es simplemente Q(t)=eΩtQ(t)=e^{\Omega t}
Utilizando las técnicas sobre los sistemas lineales, podemos demostrar que, para cada tt, Q(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×33\times 3.
En primer lugar, determinamos la forma canónica de Jordan para Ω\Omega y la utilizamos para calcular eΩte^{\Omega t}.
Con Ω\Omega expresada en la forma (22), es fácil calcular det(Ω−λI)=λ(λ2∣ω∣2),\mathrm{det}(\Omega - \lambda I)=\lambda(\lambda^2 \vert \omega\vert^2), y por lo tanto los autovalores de Ω\Omega son λ=0\lambda= 0, ±∣ω∣i\pm\vert\omega\vert i. Además, puesto que Ωω=ω×ω=0\Omega\omega=\omega\times\omega=0, se deduce que ω\omega es un autovector correspondiente al autovalor 00. Uno puede verificar que un autovector complejo vv correspondiente al autovalor ∣ω∣i\vert\omega\vert i es v=[ω1ω3+∣ω∣ω2iω2ω3−∣ω∣ω11i−(ω1+ω2)2]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 vv en término de su parte real e imaginaria v=u+qiv=u+qi y llamando P=(q,u,ω)P = (q, u, \omega) a la matriz formada a partir de colocar a qq, uu y ω\omega como columnas, obtenemos de la teoría general (o cálculo directo) que P−1ΩP=J=[0−∣ω∣0∣ω∣00000]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∈R3c\in \mathbb{R}^3, sea b=P−1cb=P^{-1}c. Tenemos c=Pb=b1q+b2u+b3ω,c=Pb=b_1q+b_2u+b_3\omega, donde b=(b1,b2,b3)b=(b_1,b_2,b_3). Por lo tanto,
El movimiento de cc descrito por la ecuación (23) es un movimiento circular sobre el eje: gen{ω}\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 ω1,ω2,ω3\omega_1,\omega_2,\omega_3 para denotar las componentes de ω\omega relativas a la base propia para AA Esperemos que esto no cause confusión con lo que hicimos anteriormente, donde los ωi′s\omega_i'\mathrm{s} eran las componentes de ω\omega con respecto a la base estándar para R3\mathbb{R}^3 (y eran constantes). Ahora de Aei=IieiAe_i = I_ie_i, i=1,2,3i = 1, 2, 3, se deduce que y
Haciendo uso de esto y de un cálculo directo de ω×Aω\omega \times A\omega en términos de los ei′se_i'\mathrm{s}, se puede calcular fácilmente la ecuación de Euler: Aω˙+ω×Aω=0,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:

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 I1,I2,I3I_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 ω:I→R3\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 ii-ésima ecuación en el sistema por IiωiI_i\omega_i (para i=1,2,3i = 1, 2, 3) y luego sumamos las tres ecuaciones juntas y tenemos, I1ω1ω˙1+I2ω2ω˙2+I3ω3ω˙3=(I2−I3)ω1ω2ω3+(I3−I1)ω1ω2ω3+(I1−I2)ω1ω2ω3=0.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 kk, que depende de la curva integral ω\omega, tal que:

Ley de Conservación I (energía cinética rotacional): para todo t∈It \in I.

Una segunda ley de conservación resulta, algo menos transparente, de multiplicar la ii-ésima ecuación en el sistema por I22ωiI_2^2\omega_i (para i=1,2,3i = 1, 2, 3) y luego sumar las tres ecuaciones juntas. Ahora tenemos la expresión I12ω1ω˙1+I22ω2ω˙2+I32ω3ω˙3=[(I2−I3)I1+(I3−I1)I2+(I1−I2)I3+ω1ω2ω3=0.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>0a> 0, que depende de la curva integral ω\omega, tal que:

Ley de Conservación II (momento angular): para todo t∈It \in I.

Es fácil ver a partir de las ecuaciones (24) - (25) que si R˙(0)=0\dot{\textbf{R}}(0) = 0, entonces la constante k=Aω⋅ω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 R(0)=0\textbf{R} (0) = 0, entonces la constante a=∣Aω∣=∣L∣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 R3\mathbb{R}^3 con respecto a los ejes principales e1,e2,e3e_1, e_2, e_3. A pesar de que en las ecuaciones ω1,ω2,ω3\omega_1, \omega_2, \omega_3 son funciones dependientes del tiempo, cuando no hay peligro de confusión, usamos ω1,ω2,ω3\omega_1, \omega_2, \omega_3 para denotar las coordenadas de puntos en R3\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 ω:I→R3\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 a22k\frac{a^2}{2k} y luego restar la ecuación (31) de él para conseguir
Esta ecuación describe un cono doble en R3\mathbb{R}^3 en función de las constantes I1,I2,I3,k,aI_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 2kI1\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 aa, se cruza con este elipsoide.
Si denotamos por bi=Ii(a22k−Ii)b_i = I_i (\frac{a^2}{2k} - I_i), para i=1,2,3i = 1, 2, 3, entonces cada uno de ellos es negativo cuando $\frac{a^2}{2k} I_3.Porlotantolaecuacioˊnparaelcono. Por lo tanto la ecuación para el cono b1ω12+b2ω22+b3ω32=0b_1\omega_1^2+b_2\omega_2^2+b_3\omega_3^2=0notienesolucioˊnencualquieradeestoscasos,yasıˊnohaycono.Otraformadedecirestoeslasiguiente:lasmagnitudesdelmomentoangularylaenergıˊacineˊticarotacional, 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 
#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 kk fijo, restringimos de modo que a22k\frac{a^2}{2k} se encuentre en el intervalo [I1,I3][I_1, I_3] y miramos como varía la familia de conos resultante C=CaC = 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 0,, b_3> 0y,porsupuesto,estosignificaque y, por supuesto, esto significa que \omega_2 = 0,, \omega_3 = 0yel"cono" y el "cono" Csereduceaunalıˊneaalolargode se reduce a una línea a lo largo de e_1.Porlotanto,. Por lo tanto, Ceselprimerejeprincipal.Acontinuacioˊn 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 aa tal que $I_1 <\frac{a^2}{2k} 0,mientrasque, mientras que b_2 <0,, b_3 <0,yporlotantolaecuacioˊnpara, y por lo tanto la ecuación para Cpuedeserescritacomo puede ser escrita como ω12=−b2b1ω22+−b3b1ω32.\omega_1^2=\frac{-b_2}{b_1}\omega_2^2+\frac{-b_3}{b_1}\omega_3^2.Porlotanto, Por lo tanto, Cesunconoelıˊpticoconsuejecoincidiendoconelprimerejeprincipal.Lainterseccioˊn 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 aa aumenta, pero luego, cuando aa alcanza el valor para el cual a22k=I2\frac{a^2}{2k} = I_2, C∩EC \cap E cambia naturalmente. Para este valor de aa, los coeficientes b1>0b_1>0, b2=0b_2 = 0, b3<0b_3 <0 y por lo que la ecuación para CC puede escribirse como ω1=±−b3b1ω3.\omega_1=\pm \sqrt{\frac{-b_3}{b_1}}\omega_3.
Por lo tanto, CC es un par de planos, cada uno de los cuales es perpendicular al plano {e1,e3}\lbrace e_1, e_3\rbrace. La intersección de este par de planos con el elipsoide EE 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 e2e_2 al otro punto fijo en este eje.
Para un valor ligeramente mayor de aa, de modo que $I_2 <\frac{a^2}{2k} 0,, b_2> 0,mientrasque, mientras que b_3 <0.Asıˊlaecuacioˊnparaelconosepuedeescribircomo. Así la ecuación para el cono se puede escribir como ω32=b1b3ω12+b2b3ω22.\omega_3^2=\frac{b_1}{b_3}\omega_1^2+\frac{b_2}{b_3}\omega_2^2.$
Por lo tanto, el cono CC es un cono elíptico con su eje a lo largo del tercer eje director. Como antes, la intersección de C∩EC \cap E es un par de curvas cerradas circulares, pero ahora centrado alrededor del tercer eje principal.
Finalmente, cuando aa es lo más grande posible, a22k=I3\frac{a^2}{2k} = I_3 y el cono CC degenera en una línea recta correspondiente al tercer eje principal. Entonces C∩EC \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
Para aa mínimo, tenemos a22k=I1=I2\frac{a^2}{2k}= I_1 = I_2, por lo que el "cono" CC se degenera en plano w3=0w_3 = 0. De hecho, para este valor de aa, C=gen{e1,e2}=EI1C = \mathrm{gen}\lbrace e_1, e_2\rbrace = E_{I_1} es el espacio característico bidimensional correspondiente al autovalor I1=I2I_1 = I_2. Por lo tanto la intersección C∩EC \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 aa tiene un valor un poco más grande de modo que $I_2 <\frac{a^2}{2k} 
#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 aa es tal que a22k=I3\frac{a^2}{2k} = I_3, el cono CaC_a se degenera en una línea recta que coincide con el tercer eje principal. De ahí C∩EC\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 EI1=gen{e1,e2}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 EI1=gen{e3}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: I1=I2=I3I_1 = I_2 = I_3, está claro que sólo hay un valor de aa para el cual CC no está vacío, es decir, para a=2kI3a =\sqrt{2kI_3}, y que para este valor, C=R3C = \mathbb{R}^3. En este caso todo el espacio R3\mathbb{R}^3 consta de puntos fijos (estables) de la ecuación de Euler. Entonces también, C∩E=EC\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 ω:I→R3\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