Diferencia entre revisiones de «Movimiento de partículas. Grupo 22C»

Revisión del 16:47 5 dic 2014

Consideramos un conjunto de 20 partículas que inicialmente se encuentran en los puntos de coordenadas (x_i,y_i,z_i) = (cos(2πi/10),sin(2πi/10), i/10), i=1, 2, ..., 20 respecto a la base ortonormal {e₁,e₂,e₃}. Supondremos que las partículas están unidas por alambres de masa despreciable de manera que su posición relativa no cambia.

1 . VISUALIZACIÓN DE UN SISTEMA DE PARTÍCULAS

Creamos un programa para obtener representación de la posición de las partículas Los ejes de la siguiente figura están fijados en la región [-2,2]x[-2,2]x[0,2]. construimos los vectores con las coordenadas X, Y, Z de los 20 puntos y los representamos.

Visualizacion de un sistema de particulas

%la primera parte del programa calcula las coordenadas de los puntos
%materiales.
%crea una matriz de 3 filas y 20 columnas, donde cada columna son las
%tres coordenadas de los puntos.
coordenadas=zeros(3,20);
for i=1:20
coordenadas(1,i)=(cos((2*pi*i)/10));
coordenadas(2,i)=(sin((2*pi*i)/10));
coordenadas(3,i)=(i/10);
end
clear('i');
%dibujamos la helice en 3d con los ejes pedidos, 
%creamos los vectores de coordeandas X, Y, Z
coorX=zeros(1,20);
coorY=zeros(1,20);
coorZ=zeros(1,20);
for i=1:20
coorX(1,i)=coordenadas(1,i);
coorY(1,i)=coordenadas(2,i);
coorZ(1,i)=coordenadas(3,i);
end
clear('i');
%pintamos esos vectores
figure(1);
plot3(coorX,coorY,coorZ, '--. b');
axis([-2,2,-2,2,0,2]);
hold on

2 . CENTRO DE MASAS DE UN SISTEMA DE PARTÍCULAS

Sabiendo que las partículas tienen masa m_i=10+i/10, para calcular el centro de masas hacemos uso de la siguiente fórmula: (∑_iⁿr_im_i)/M donde M es la masa total y ri= la distancia del punto al eje correspondiente X, Y ó Z. Calculamos los valores ponderados que cada punto material aporta al centro de masas, los sumamos y representamos las coordenadas (XG,YG,ZG) en color verde.

Visualización del centro de masas

%calculamos una matriz con las masas
masas=zeros(1,20);
for i=1:20
masas(1,i)=(10+i/10);
end
clear('i');
%calculamos cuanto aporta cada particula al centro de masas
%y el centro de masas
auxXG=zeros(1,20);
auxYG=zeros(1,20);
auxZG=zeros(1,20);
for i=1:20
auxXG(1,i)=(masas(1,i)*coorX(1,i));
auxYG(1,i)=(masas(1,i)*coorY(1,i));
auxZG(1,i)=(masas(1,i)*coorZ(1,i));
end
clear('i');
masatotal=sum(masas);
auxXG=sum(auxXG);
auxYG=sum(auxYG);
auxZG=sum(auxZG);
%calculamos las coordenads del centro de masas y lo pintamos
XG=((1/masatotal)*auxXG);
YG=((1/masatotal)*auxYG);
ZG=((1/masatotal)*auxZG);
plot3(XG,YG,ZG,'. g ','linewidth',10)
hold off
figure(2);

3 ROTACIÓN DE UN SISTEMA DE PARTÍCULAS

3.1 Matriz rotación con eje ω=e₃ y angulo θ= π/16

Con la definición tensorial de la rotación, hacemos un programa que calcule las componentes de dicho tensor para un determinado ángulo θ=π/16 alrededor de un eje generado por el tensor unitario w=e₃. La matriz de rotación obtenida por el programa es la siguiente: [math]\begin{pmatrix} 0,9808 & -0,1951 &0 \\ 0,1951&0,9808 &0 \\ 0&0 & 1 \end{pmatrix}[/math] El código que nos permite obtener dicha matriz de componentes de la rotación es este:

%esta parte del programa calcula la primera matriz de rotación
%calculamos la matriz de la rotación a partir de un eje de vector w
%y un ángulo theta
w=[0, 0 , 1];
theta=pi/16;
matriz1=eye(3);
matriz2=zeros(3,3);
matriz3=zeros(3,3);
for i=1:3
matriz2(1,i)=(w(1)*w(i));
matriz2(2,i)=(w(2)*w(i));
matriz2(3,i)=(w(3)*w(i));
end
clear('i');
matriz3(1,1)=0;
matriz3(1,2)=(-1*w(3));
matriz3(1,3)=(w(2));
matriz3(2,1)=(w(3));
matriz3(2,2)=0;
matriz3(2,3)=(-1*w(1));
matriz3(3,1)=(-1*w(2));
matriz3(3,2)=(w(1));
matriz3(3,3)=0;
%construimos la matriz de rotación
rotacion1=zeros(3,3);
matrizaux1=(cos(theta)*matriz1);
matrizaux2=((1-cos(theta))*matriz2);
matrizaux3=(sin(theta)*matriz3);
rotacion1=matrizaux1+matrizaux2+matrizaux3;

Una vez tenemos la matriz del tensor rotación, realizamos el producto contraído de dicho tensor con los radiovectores de los puntos materiales para realizarles el giro pedido. Representamos las nuevas coordenadas de los puntos junto con las originales para que el lector pueda apreciar la transformación

Visualización de un sistema de puntos rotados con eje ω=e₃ y angulo θ= π/16

%separamos los vectores directores de cada particula
%y les hacemos la rotación
coorxR1=zeros(1,20);
cooryR1=zeros(1,20);
coorzR1=zeros(1,20);
vaux=zeros(3,1);
for i=1:20
vaux(1,1)=coorX(1,i);
vaux(2,1)=coorY(1,i);
vaux(3,1)=coorZ(1,i);
rotpunto=rotacion1*vaux;
coorxR1(1,i)=rotpunto(1);
cooryR1(1,i)=rotpunto(2);
coorzR1(1,i)=rotpunto(3);
clear('rotpunto');
end
clear('i');
plot3(coorxR1,cooryR1, coorzR1, '--. g')
axis([-2,2,-2,2,0,2])
hold on
plot3(coorX, coorY, coorZ, '--p b')
hold off

3.2 Matriz rotación con eje ω=e₁ y angulo θ= π/16

Código:

Visualización de un sistema de puntos rotados con eje ω=e₁ y angulo θ= π/16

%segunda matriz de rotacion, eje e1
w=[1, 0 , 0];
theta=pi/16;
matriz1=eye(3);
matriz2=zeros(3,3);
matriz3=zeros(3,3);
for i=1:3
matriz2(1,i)=(w(1)*w(i));
matriz2(2,i)=(w(2)*w(i));
matriz2(3,i)=(w(3)*w(i));
end
clear('i');
matriz3(1,1)=0;
matriz3(1,2)=(-1*w(3));
matriz3(1,3)=(w(2));
matriz3(2,1)=(w(3));
matriz3(2,2)=0;
matriz3(2,3)=(-1*w(1));
matriz3(3,1)=(-1*w(2));
matriz3(3,2)=(w(1));
matriz3(3,3)=0;
%construimos la matriz de rotación
rotacion2=zeros(3,3);
matrizaux1=(cos(theta)*matriz1);
matrizaux2=((1-cos(theta))*matriz2);
matrizaux3=(sin(theta)*matriz3);
rotacion2=matrizaux1+matrizaux2+matrizaux3;
clear ('vaux');
%separamos los vectores directores de cada particula y les hacemos la rotación
coorxR2=zeros(1,20);
cooryR2=zeros(1,20);
coorzR2=zeros(1,20);
vaux=zeros(3,1);
for i=1:20
vaux(1,1)=coorX(1,i);
vaux(2,1)=coorY(1,i);
vaux(3,1)=coorZ(1,i);
rotpunto=rotacion2*vaux;
coorxR2(1,i)=rotpunto(1);
cooryR2(1,i)=rotpunto(2);
coorzR2(1,i)=rotpunto(3);
clear('rotpunto');
end
clear('i');
plot3(coorxR2,cooryR2, coorzR2, '--p r')
axis([-2,2,-2,2,0,2])
hold on
plot3(coorX, coorY, coorZ, '--. b')
hold off

La matriz de componentes del tensor de la segunda rotación es:

   1.0000         0         0
        0    0.9808   -0.1951
        0    0.1951    0.9808

3.3 Matriz rotación con eje ω=e₂ y angulo θ= π/16

Calculamos la matriz de rotación a partir de un eje ω=e₂ y un ángulo θ= π/16. El código que nos permite obtener dicha matriz de rotación es esta:

Visualización de un sistema de puntos rotados con eje ω=e₂ y ángulo θ= π/16

%calculamos la matriz de la tercera rotación
%a partir de un vector w y un eje angulo theta
w=[0, 1, 0];
theta=pi/16;
matriz1=eye(3);
matriz2=zeros(3,3);
matriz3=zeros(3,3);
for i=1:3
matriz2(1,i)=(w(1)*w(i));
matriz2(2,i)=(w(2)*w(i));
matriz2(3,i)=(w(3)*w(i));
end
clear('i');
matriz3(1,1)=0;
matriz3(1,2)=(-1*w(3));
matriz3(1,3)=(w(2));
matriz3(2,1)=(w(3));
matriz3(2,2)=0;
matriz3(2,3)=(-1*w(1));
matriz3(3,1)=(-1*w(2));
matriz3(3,2)=(w(1));
matriz3(3,3)=0;
%construimos la matriz de rotación
rotacion3=zeros(3,3);
matrizaux1=(cos(theta)*matriz1);
matrizaux2=((1-cos(theta))*matriz2);
matrizaux3=(sin(theta)*matriz3);
rotacion3=matrizaux1+matrizaux2+matrizaux3;
clear('vaux');
%separamos los vectores directores de cada particula y les hacemos la rotación
coorxR3=zeros(1,20);
cooryR3=zeros(1,20);
coorzR3=zeros(1,20);
vaux=zeros(3,1);
for i=1:20
vaux(1,1)=coorX(1,i);
vaux(2,1)=coorY(1,i);
vaux(3,1)=coorZ(1,i);
rotpunto=rotacion3*vaux;
coorxR3(1,i)=rotpunto(1);
cooryR3(1,i)=rotpunto(2);
coorzR3(1,i)=rotpunto(3);
clear('rotpunto');
end
clear('i');
plot3(coorxR3,cooryR3, coorzR3, '--p k')
axis([-2,2,-2,2,0,2])
hold on
plot3(coorX, coorY, coorZ, '--. b')
hold off

La matriz de componentes del tensor de la tercera rotación es:

   0.9808         0    0.1951
        0    1.0000         0
  -0.1951         0    0.9808

3.4 Matriz rotación con eje ω=e₁+e₂+e₃ y ángulo θ= π/16

Visualizacion de un sistema de puntos rotados con eje ω=e₁+e₂+e₃ y ángulo θ= π/16

%calculamos la matriz de la cuarta rotación
%a partir de un vector w y un eje angulo theta
a=1/sqrt(3);
w=[a, a, a];clc
theta=pi/16;
matriz1=eye(3);
matriz2=zeros(3,3);
matriz3=zeros(3,3);
for i=1:3
matriz2(1,i)=(w(1)*w(i));
matriz2(2,i)=(w(2)*w(i));
matriz2(3,i)=(w(3)*w(i));
end
clear('i');
matriz3(1,1)=0;
matriz3(1,2)=(-1*w(3));
matriz3(1,3)=(w(2));
matriz3(2,1)=(w(3));
matriz3(2,2)=0;
matriz3(2,3)=(-1*w(1));
matriz3(3,1)=(-1*w(2));
matriz3(3,2)=(w(1));
matriz3(3,3)=0;
%construimos la matriz de rotación
rotacion4=zeros(3,3);
matrizaux1=(cos(theta)*matriz1);
matrizaux2=((1-cos(theta))*matriz2);
matrizaux3=(sin(theta)*matriz3);
rotacion4=matrizaux1+matrizaux2+matrizaux3;
clear('vaux');
%separamos los vectores directores de cada partícula y les hacemos la rotación
coorxR4=zeros(1,20);
cooryR4=zeros(1,20);
coorzR4=zeros(1,20);
vaux=zeros(3,1);
for i=1:20
vaux(1,1)=coorX(1,i);
vaux(2,1)=coorY(1,i);
vaux(3,1)=coorZ(1,i);
rotpunto=rotacion4*vaux;
coorxR4(1,i)=rotpunto(1);
cooryR4(1,i)=rotpunto(2);
coorzR4(1,i)=rotpunto(3);
clear('rotpunto');
end
clear('i');
plot3(coorxR4,cooryR4, coorzR4, '--p m')
axis([-2,2,-2,2,0,2])
hold on
plot3(coorX, coorY, coorZ, '--. b')
clear ('matriz1','matriz2','matriz3','w','matrizaux1','matrizaux2','matrizaux3');
hold off
%fin del apartado 3

La matriz de componentes del tensor de la tercera rotación es:

   0.9872   -0.1062    0.1190
   0.1190    0.9872   -0.1062
  -0.1062    0.1190    0.9872

4 DEFINICIÓN DE LA VELOCIDAD COMO UN TENSOR ANTISIMÉTRICO

Dado un sistema que gira alrededor de un eje genérico con vector de dirección unitario [math]\overrightarrow{\omega}[/math] cuya variación angular viene dada por la función θ(t) donde t ∈ [0, π] es el tiempo, comprobamos analíticamente que la velocidad de los puntos del sistema está relacionada con su posición mediante un tensor antisimétrico A de la siguiente forma [math]v_{i}(t)=\frac{dr_{i}(t)}{dt}=A\cdot r_{i}(t)[/math] para todo i = 1, 2, ..., 20.

Para poder expresar la velocidad angular [math]\Omega [/math] en forma vectorial, primero multiplicamos su módulo [math]\dot{\Theta}(t)[/math] por su vector director unitario [math]\overrightarrow{\omega}[/math] [math]\rightarrow \Omega =\left | \dot{\Theta}(t) \right |\cdot \overrightarrow{\omega}[/math].

Por el campo de velocidades sabemos que [math]\overrightarrow{V_{i}}=\overrightarrow{\Omega}\times \overrightarrow{OP}=\overrightarrow{\Omega}\times \overrightarrow{r_{i}}[/math].

En forma tensorial [math]\overrightarrow{V_{i}}=(\overrightarrow{\Omega }\times )\cdot \overrightarrow{r_{i}}[/math].

A continuación se demuestra que [math]\overrightarrow{\Omega }\times [/math] es un tensor antisimétrico donde [math]\Omega =\overrightarrow{\omega}\times [/math] :

[math]\Omega _{ij}=\overrightarrow{g_{i}}\cdot [\Omega \cdot \overrightarrow{g_{j}}]=\overrightarrow{g_{i}}\cdot [(\overrightarrow{\omega}\times )\cdot \overrightarrow{g_{j}}]=\overrightarrow{g_{i}}\cdot (\overrightarrow{\omega}\times \overrightarrow{g_{j}})=[\overrightarrow{g_{i}},\overrightarrow{\omega},\overrightarrow{g_{j}}]=-[\overrightarrow{g_{i}},\overrightarrow{g_{j}},\overrightarrow{\omega}]=-\omega ^{k}[\overrightarrow{g_{i}},\overrightarrow{g_{j}},\overrightarrow{g_{k}}]=-\omega ^{k}\pm \sqrt{g}\varepsilon _{ijk}=-\omega ^{1}\cdot \varepsilon _{ij1}-\omega ^{2}\cdot \varepsilon _{ij2}-\omega ^{3}\cdot \varepsilon _{ik3}[/math] [math]-\omega^{1}\cdot \bigl(\begin{smallmatrix} 0 & 0 & 0\\ 0 & 0 & 1\\ 0 & -1 & 0 \end{smallmatrix}\bigr)-\omega ^{2}\cdot \bigl(\begin{smallmatrix} 0 & 0 & -1\\ 0 & 0 & 0\\ 1 & 0 & 0 \end{smallmatrix}\bigr)-\omega ^{3}\cdot \bigl(\begin{smallmatrix} 0 & 1 & 0\\ -1 & 0 & 0\\ 0 & 0 & 0 \end{smallmatrix}\bigr)=\bigl(\begin{smallmatrix} 0 & -\omega ^{3} & \omega ^{2}\\ \omega ^{3} & 0 & -\omega ^{1}\\ -\omega ^{2} & \omega ^{1} & 0 \end{smallmatrix}\bigr)[/math]

Por consiguiente se puede demostrar que el vector axial asociado a [math]\overrightarrow{\Omega}\times [/math] es [math]\dot{\Theta }(t)\cdot \overrightarrow{\omega } [/math] multiplicando la matriz de coordenadas de [math]\Omega [/math] por [math]\dot{\Theta }(t)[/math].

5 VISUALIZACIÓN DE LOS VECTORES VELOCIDAD

Visualización de los vectores velocidad

véase que los vectores son tangentes a la trayectoria

%pintamos los las partículas y calculamos
%los vectores de velocidad
plot3(coorX,coorY,coorZ, '-. b')
axis([-2,2,-2,2,0,2]);
hold on
%calculamos las coordenadas de los vectores velocidad
%de los puntos materiales
coorVx=zeros(1,20);
coorVy=zeros(1,20);
coorVz=zeros(1,20);
omega=[0,0,1];
radioaux=zeros(1,3);
veloaux=zeros(1,3);
for i=1:20
radioaux(1,1)=coorX(1,i);
radioaux(1,2)=coorY(1,i);
radioaux(1,3)=coorZ(1,i);
veloaux=cross(omega,radioaux);
coorVx(1,i)=veloaux(1);
coorVy(1,i)=veloaux(2);
coorVz(1,i)=veloaux(3);
clear('veloaux','radioaux');
radioaux=zeros(1,3);
veloaux=zeros(1,3);
end
for i=1:20
xaux=coorX(1,i);
yaux=coorY(1,i);
zaux=coorZ(1,i);
vxaux=coorVx(1,i);
vyaux=coorVy(1,i);
vzaux=coorVz(1,i);
%pintamos el vector velocidad con el comando quiver
quiver3(xaux,yaux,zaux,vxaux,vyaux,vzaux)
clear('xaux','yaux','zaux','vxaux','vyaux','vzaux');
end
hold off

6 TENSOR DE INERCIA

El momento angular del sistema se define por [math]L=\sum_{i=1}^{20}\overrightarrow {r_{i}}\times m_{i}\cdot\overrightarrow {v_{i}}[/math]. Sustituyendo el valor de la velocidad [math]\overrightarrow{v_{i}}=\overrightarrow{\omega}\times \overrightarrow{r_{i}} [/math], donde [math]\overrightarrow{\omega}[/math] es la velocidad angular, comprobamos analíticamente que podemos escribir [math]L=I\cdot \overrightarrow{\omega}[/math] donde I es un tensor de orden 2 conocido como tensor de inercia y sus componentes momentos de inercia.

[math]L=\sum_{i=1}^{20}\overrightarrow {r_{i}}\times m_{i}\cdot\overrightarrow {v_{i}}[/math] [math]=\sum_{i=1}^{20}\overrightarrow{r_{i}}\cdot (m_{i}\cdot \overrightarrow{\omega}\times \overrightarrow{r_{i}})[/math] [math]=\sum_{i=1}^{20}\begin{vmatrix} \overrightarrow{\omega}\cdot m_{i} & \overrightarrow{r_{i}}\\ \overrightarrow{r_{i}}\cdot \overrightarrow{\omega}\cdot m_{i} & \overrightarrow{r_{i}}\cdot \overrightarrow{r_{i}} \end{vmatrix}[/math] [math]=\sum_{i=1}^{20}\left [ \overrightarrow{\omega}\cdot m_{i}\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{r_{i}})\cdot \overrightarrow{r_{i}}\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{\omega}\cdot m_{i}) \right ][/math] [math]=\sum_{i=1}^{20}\left [m_{i} \cdot (\overrightarrow{\omega}\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{r_{i}}))-\overrightarrow{r_{i}}(\overrightarrow{r_{i}}\cdot \overrightarrow{\omega})\right ][/math] [math]=\sum_{i=1}^{20}[m_{i}[\overrightarrow{\omega}\cdot \left | \overrightarrow{r_{i}} \right |^{2}-\overrightarrow{r_{i}}(\overrightarrow{r_{i}}\cdot \overrightarrow{\omega})]][/math] [math]=\sum_{i=1}^{20}m_{i}\cdot [\left | \overrightarrow{r_{i}} \right |^{2}\cdot 1-\overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}} ]\cdot \overrightarrow{\omega}[/math] [math]=\sum_{i=1}^{20}[m_{i}\cdot \left | \overrightarrow{r_{i}} \right |^{2}\cdot 1-m_{i}\cdot \overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}}]\overrightarrow{\omega}=L[/math] [math]=\frac{1}{2}m_{i}[(\overrightarrow{\omega}\cdot \overrightarrow{\omega})\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{r_{i}})-(\overrightarrow{\omega}\cdot \overrightarrow{r_{i}})\cdot (\overrightarrow{\omega}\cdot \overrightarrow{r_{i}})]=\frac{1}{2}m_{i}\cdot \overrightarrow{\omega}[\overrightarrow{\omega}\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{r_{i}})-\overrightarrow{r_{i}}\cdot (\overrightarrow{\omega}\cdot \overrightarrow{r_{i}})]=\frac{1}{2}m_{i}\cdot \overrightarrow{\omega}[\left | \overrightarrow{r_{i}} \right |^{2}\cdot \overrightarrow{\omega}-(\overrightarrow{r_{i}}\cdot \overrightarrow{\omega})\cdot \overrightarrow{r_{i}}][/math] [math]=\frac{1}{2}\cdot \overrightarrow{\omega}[m_{i}\cdot (\left | \overrightarrow{r_{i}} \right |^{2}\cdot \overrightarrow{\omega})-m_{i}\cdot [(\overrightarrow{r_{i}}\cdot \overrightarrow{\omega})\cdot \overrightarrow{r_{i}}]]=\frac{1}{2}\overrightarrow{\omega}\cdot [m_{i}\cdot \left | \overrightarrow{r_{i}} \right |^{2}\cdot 1-m_{i}\cdot \overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}}]\cdot \overrightarrow{\omega}=\frac{1}{2}\cdot \overrightarrow{\omega}\cdot I\cdot \overrightarrow{\omega}[/math]

6.1 Cálculo de las componentes del tensor I de inercia en la base {e₁,e₂,e₃}

%calculamos el tensor de inercia
%primera matriz del tensor% %base i,j,k
%calculamos la primera matriz
clear ('vaux','i');
primeramatriz=zeros(3,3);
identidad=eye(3);
for i=1:20
primeramatrizaux=zeros(3,3);
vaux=zeros(1,3);
vaux(1,1)=coorX(1,i);
vaux(1,2)=coorY(1,i);
vaux(1,3)=coorZ(1,i);
modulocuadrado=vaux*vaux';
masa=masas(i);
momento=masa*modulocuadrado;
primeramatrizaux=momento*identidad;
primeramatriz=primeramatriz+primeramatrizaux;
clear('masa','primeramatrizaux','vaux','modulocuadrado','momento');
end
clear('i');
%calculamos la segunda matriz
segundamatriz=zeros(3,3);
clear ('vaux');
clear ('masa');
for i=1:20
segundamatrizaux=zeros(3,3);
segundamatrizauxaux=zeros(3,3);
vaux=zeros(1,3);
vaux(1,1)=coorX(1,i);
vaux(1,2)=coorY(1,i);
vaux(1,3)=coorZ(1,i);
masa=masas(i);
segundamatrizaux(1,1)=vaux(1,1)*vaux(1,1);
segundamatrizaux(1,2)=vaux(1,1)*vaux(1,2);
segundamatrizaux(1,3)=vaux(1,1)*vaux(1,3);
segundamatrizaux(2,1)=vaux(1,2)*vaux(1,1);
segundamatrizaux(2,2)=vaux(1,2)*vaux(1,2);
segundamatrizaux(2,3)=vaux(1,2)*vaux(1,3);
segundamatrizaux(3,1)=vaux(1,3)*vaux(1,1);
segundamatrizaux(3,2)=vaux(1,3)*vaux(1,2);
segundamatrizaux(3,3)=vaux(1,3)*vaux(1,3);
segundamatrizauxaux=masa*segundamatrizaux;
segundamatriz=segundamatriz+segundamatrizauxaux;
clear ('vaux','segundamatrizaux','segundamatrizauxaux');
end
clear('i');
tensordeinercia=zeros(3,3);
tensordeinercia=primeramatriz+segundamatriz

Ahora lo calculamos suponiendo que ω=e₃:

w=[0;0;1];
momentoejez=tensordeinercia*w;
momentoejez

7 RELACIÓN DEL TENSOR DE INERCIA Y LA ENERGÍA CINÉTICA

La energía cinética total del sistema se define como [math]E_{c}=\sum_{i=1}^{20}\frac{1}{2}\cdot m_{i}\cdot \left | \overrightarrow{v_{i}} \right |^{2}[/math]. Sustituyendo el valor de la velocidad [math]\overrightarrow{v_{i}}=\overrightarrow{\omega}\times \overrightarrow{r_{i}} [/math], comprobamos que podemos escribir [math]E_{c}=\frac{1}{2}\cdot \overrightarrow{\omega}\cdot I\cdot \overrightarrow{\omega}[/math], donde I es el tensor de inercia.

[math]L=\sum_{i=1}^{20}\overrightarrow {r_{i}}\times m_{i}\cdot\overrightarrow {v_{i}}[/math] [math]=\sum_{i=1}^{20}\overrightarrow{r_{i}}\cdot (m_{i}\cdot \overrightarrow{\omega}\times \overrightarrow{r_{i}})[/math] [math]=\sum_{i=1}^{20}\begin{vmatrix} \overrightarrow{\omega}\cdot m_{i} & \overrightarrow{r_{i}}\\ \overrightarrow{r_{i}}\cdot \overrightarrow{\omega}\cdot m_{i} & \overrightarrow{r_{i}}\cdot \overrightarrow{r_{i}} \end{vmatrix}[/math] [math]=\sum_{i=1}^{20}\left [ \overrightarrow{\omega}\cdot m_{i}\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{r_{i}})\cdot \overrightarrow{r_{i}}\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{\omega}\cdot m_{i}) \right ][/math] [math]=\sum_{i=1}^{20}\left [m_{i} \cdot (\overrightarrow{\omega}\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{r_{i}}))-\overrightarrow{r_{i}}(\overrightarrow{r_{i}}\cdot \overrightarrow{\omega})\right ][/math] [math]=\sum_{i=1}^{20}[m_{i}[\overrightarrow{\omega}\cdot \left | \overrightarrow{r_{i}} \right |^{2}-\overrightarrow{r_{i}}(\overrightarrow{r_{i}}\cdot \overrightarrow{\omega})]][/math] [math]=\sum_{i=1}^{20}m_{i}\cdot [\left | \overrightarrow{r_{i}} \right |^{2}\cdot 1-\overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}} ]\cdot \overrightarrow{\omega}[/math] [math]=\sum_{i=1}^{20}[m_{i}\cdot \left | \overrightarrow{r_{i}} \right |^{2}\cdot 1-m_{i}\cdot \overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}}]\overrightarrow{\omega}=L[/math] [math]=\frac{1}{2}m_{i}[(\overrightarrow{\omega}\cdot \overrightarrow{\omega})\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{r_{i}})-(\overrightarrow{\omega}\cdot \overrightarrow{r_{i}})\cdot (\overrightarrow{\omega}\cdot \overrightarrow{r_{i}})]=\frac{1}{2}m_{i}\cdot \overrightarrow{\omega}[\overrightarrow{\omega}\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{r_{i}})-\overrightarrow{r_{i}}\cdot (\overrightarrow{\omega}\cdot \overrightarrow{r_{i}})]=\frac{1}{2}m_{i}\cdot \overrightarrow{\omega}[\left | \overrightarrow{r_{i}} \right |^{2}\cdot \overrightarrow{\omega}-(\overrightarrow{r_{i}}\cdot \overrightarrow{\omega})\cdot \overrightarrow{r_{i}}][/math] [math]=\frac{1}{2}\cdot \overrightarrow{\omega}[m_{i}\cdot (\left | \overrightarrow{r_{i}} \right |^{2}\cdot \overrightarrow{\omega})-m_{i}\cdot [(\overrightarrow{r_{i}}\cdot \overrightarrow{\omega})\cdot \overrightarrow{r_{i}}]]=\frac{1}{2}\overrightarrow{\omega}\cdot [m_{i}\cdot \left | \overrightarrow{r_{i}} \right |^{2}\cdot 1-m_{i}\cdot \overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}}]\cdot \overrightarrow{\omega}=\frac{1}{2}\cdot \overrightarrow{\omega}\cdot I\cdot \overrightarrow{\omega}[/math]

%calculamos la energía cinética del sistema de partículas
energiaux=tensordeinercia*w;
W=(1/2)*w;
energiacinetica=(W')*energiaux;
energiacinetica

8 TEOREMA DE STEINER

El Teorema de Steiner permite relacionar el tensor de inercia respecto a un eje de giro que pasa por el centro de masas I_G con el tensor de inercia correspondiente a un eje de giro que pasa por cualquier otro punto I_P [math]\Rightarrow I_{P}=I_{G}+m\left \| \overrightarrow{a} \right \|^{2}\cdot 1-\overrightarrow{a}\otimes \overrightarrow{a}[/math] siendo [math]\overrightarrow{a}[/math] el vector que une el centro de masas G con el punto P.

Para empezar la demostración definimos I_P=I_G+T, donde T es un tensor de orden 2. El teorema quedará demostrado si [math]T=m(\left \| \overrightarrow{a} \right \|^{2}\cdot 1-\overrightarrow{a}\otimes \overrightarrow{a})[/math].

Definimos r'=(r-a) para demostrar lo siguiente: [math]I_{P}=\sum_{i=1}^{20}m_{i}(\left | \overrightarrow{{r_{i}}'} \right |^{2}\cdot 1-\overrightarrow{{r_{i}}'}\otimes \overrightarrow{{r_{i}}'})=\sum_{i=1}^{20}m_{i}(\left | \overrightarrow{r_{i}}- \overrightarrow{a}\right |^{2}\cdot 1-(\overrightarrow{r_{i}}-\overrightarrow{a})\otimes (\overrightarrow{r_{i}}-\overrightarrow{a}))=\sum_{i=1}^{20}m_{i}([\left | \overrightarrow{r_{i}} \right |^{2}+\left | \overrightarrow{a} \right |^{2}-2\cdot (\overrightarrow{r_{i}}-\overrightarrow{a})]\cdot 1-[\overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}}-\overrightarrow{r_{i}}\otimes \overrightarrow{a}-\overrightarrow{a}\otimes \overrightarrow{r_{i}}+\overrightarrow{a}\otimes \overrightarrow{a}])[/math] [math]I_{P}=\sum_{i=1}^{20}(m_{i}\cdot \left | \overrightarrow{r_{i}} \right |^{2}\cdot 1+m_{i}\left | \overrightarrow{a} \right |^{2}\cdot 1-m_{i}\cdot 2\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{a})\cdot 1-m_{i}\cdot \overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}}+m_{i}\overrightarrow{r_{i}}\otimes \overrightarrow{a}+m_{i}\cdot \overrightarrow{a}\otimes \overrightarrow{r_{i}}-m_{i}\cdot \overrightarrow{a}\otimes \overrightarrow{a})[/math] [math]I_{P}=\sum_{i=1}^{20}(m_{i}\cdot [\left | \overrightarrow{r_{i}} \right |^{2}\cdot 1-\overrightarrow{r_{i}}\otimes \overrightarrow{r_{i}}])+\sum_{i=1}^{20}(m_{i}\cdot [\left | \overrightarrow{a} \right |^{2}\cdot 1-2\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{a})\cdot 1+\overrightarrow{r_{i}}\otimes \overrightarrow{a}+\overrightarrow{a}\otimes \overrightarrow{r_{i}}-\overrightarrow{a}\otimes \overrightarrow{a}])[/math]

Steiner quedará demostrado si [math]\sum_{i=1}^{20}m_{i}\cdot (\left | \overrightarrow{a} \right |^{2}\cdot 1+\overrightarrow{r_{i}}\otimes \overrightarrow{a}+\overrightarrow{a}\otimes \overrightarrow{r_{i}}-\overrightarrow{a}\otimes \overrightarrow{a})=\sum_{i=1}^{20}m_{i}\cdot (\left | \overrightarrow{a} \right |^{2}\cdot 1+\overrightarrow{a}\otimes \overrightarrow{a})[/math] [math]\sum_{i=1}^{20}m_{i}\cdot (\left | \overrightarrow{a} \right |^{2}\cdot 1-2\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{a} )\cdot 1+\overrightarrow{r_{i}}\otimes\overrightarrow{a}+\overrightarrow{a}\otimes \overrightarrow{r_{i}}-\overrightarrow{a}\otimes \overrightarrow{a} - \left | \overrightarrow{a} \right |^{2}\cdot 1+\overrightarrow{a}\otimes \overrightarrow{a})=0[/math] [math]\sum_{i=1}^{20}m_{i}\cdot(-2\cdot (\overrightarrow{r_{i}}\cdot \overrightarrow{a} )\cdot 1+\overrightarrow{r_{i}}\otimes \overrightarrow{a}+\overrightarrow{a}\otimes \overrightarrow{r_{i}})=0[/math] [math](-2\cdot [(\sum_{i=1}^{20}m_{i}\cdot \overrightarrow{r_{i}})\cdot \overrightarrow{a})]\cdot 1+(\sum_{i=1}^{20}m_{i}\cdot \overrightarrow{r_{i}})\otimes \overrightarrow{a}+\overrightarrow{a}\otimes (\sum_{i=1}^{20}m_{i}\cdot\overrightarrow{r_{i}} )=0[/math]. Como por definición del centro de masas [math]\sum_{i=1}^{20}m_{i}\cdot\overrightarrow{r_{i}}=0[/math]. [math](-2\cdot (0\cdot \overrightarrow{a})\cdot 1+0\otimes\overrightarrow{a}+\overrightarrow{a}\otimes 0=0 [/math] tal como se quería demostrar.