Movimiento de un sistema de partículas (Grupo G9)

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Movimiento de un sistema de partículas. Grupo G9
Asignatura Teoría de Campos
Curso 2013-14
Autores Ramón Fernández-Reyes Doallo, Antón Fernández Arriola, Pablo García González, Adrián Piñeiro Novas
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Visualización de un sistema de partículas

Para empezar a dibujar un sistema de partículas primero tenemos que generarlos. En nuestro caso, tenemos 21 partículas que inicialmente se encuentran en las coordenadas (xi ,yi ,zi) = (sin(π(i − 1)/10,cos(π(i − 1)/4), cos(π(i −1)/10) para i = 1, 2, ..., 21 respecto a la base ortonormal {e1,e2,e3}. Supondremos que nuestras partículas están unidas por un alambre de masa despreciable de manera que su posición relativa no cambia.

Visualización de la curva dada
















t = (1:1:21);              % Generamos un vector t con los 21 puntos 
x = sin(pi*(t-1)/10);      % Vectores de las componenetes de los puntos de la curva
y = cos(pi*(t-1)/4);
z = cos(pi*(t-1)/10);
plot3(x,y,z,'-x')          % Representamos la curva
axis equal                 % Igualamos los ejes 
axis ([-2,2,-2,2,-2,2])    % Definimos los ejes entre -2 y 2


2 Centro de masas de un sistema de particulas

2.1 Formula del centro de masas

En nuestro caso tenemos 21 particulas y cada particula tiene la misma masa m=10. Para calcular el centro de masas hacemos uso de la formula de centro de masas : (∑inrimasai)/M donde M es la masa total (21×10=210),"masa" la masa de cada particula y ri= la distancia del punto al eje correspondiente x, y ó z. El codigo que utilizaremos es: 

masa = 10;              %masa de cada particula
Mtotal = 21*masa;       %La masa total la obtenemos al sumar la masa de las 21 particulas
X = zeros(1,21);              
Y = zeros(1,21);        %Creamos 3 matrices de ceros que posteriormente completaremos
Z = zeros(1,21);

for i = 1:21
   X(1,i) = masa*x(i);  %Creamaos 3 matrices fila con las coordenadas de cada punto multiplicadas por la masa del mismo
   Y(1,i) = masa*y(i);                 
   Z(1,i) = masa*z(i);
end
X = sum(X);                            
Y = sum(Y);             % X Y Z serán el sumatorio de los productos de las coordenadas de los puntos por su masa      
Z = sum(Z);

XG = (X/Mtotal);        %Dividimos XYZ entre la masa total y obtenemos las coordenadas del centro de masas.
YG = (Y/Mtotal);
ZG = (Z/Mtotal);


2.2 Visualicacion del centro de masas

A continuación mostraremos un grafico con la curva y su centro de masas: 

%Representacion de la curva
t = (1:1:21);                      
x = (sin(pi*(t-1)/10));
y = (cos (pi*(t-1)/4));
z = cos(pi*(t-1)/10);

plot3(x,y,z,'-x')
axis equal 
axis ([-2,2,-2,2,-2,2])
hold on

%Representación del centro de masas
masa = 10;
Mtotal = 21*masa;
X = zeros(1,21);
Y = zeros(1,21);
Z = zeros(1,21);

for i = 1:21
   X(1,i) = masa*x(i);
   Y(1,i) = masa*y(i);
   Z(1,i) = masa*z(i);
end
XG = sum(X)/Mtotal;
YG = sum(Y)/Mtotal;
ZG = sum(Z)/Mtotal;

plot3(XG,YG,ZG,'.g','linewidth',5)
hold off
Visualización de la curva dada en azul y el centro de masas en verde
























3 Rotacion de un sistema de particulas

3.1 Matriz de componentes de una Rotación

En primer lugar, vamos a calcular analíticamente la matriz de componentes de una rotación. La fórmula de una rotación es: izquierda



donde w es el vector unitario: izquierda





Refiriéndonos a la base ortonormal {e1,e2,e3}, podemos expresar la rotación como: izquierda















Por lo que la matriz de componentes es: izquierda





3.2 Matriz Rotacion con eje ω=e3 y angulo θ= π/16

Para generar la matriz de rotacion necesitamos saber el eje y el angulo con el que queremos rotarlo. En nuestro caso estos valores seran eje ω=e3 y angulo θ= π/16. Para generar esta matriz utilizaremos: 

We3= [0,0,1];            %Eje de rotacion 
mod = sqrt(sum(We3.^2)); %Hallamos el modulo del vector rotación, aunque en este caso sea innecesario, lo haremos para un eje generico
We3 = We3./mod;          %Por ultimo hacemos que el eje sea un vector unitario
tt = pi/16;              %Definimos el angulo de rotación 

% La Matriz de rotacion,segun hemos explicado será: R=cos(θ)I + (1-cos(θ)u×u +sin(θ)ux

R1=eye(3)*cos(tt);                                % Genera la matriz corespondiente a cos(θ)I
R2=kron(We3,We3')*(1-cos(tt));                    % Genera la matriz corespondiete a (1-cos(θ)we3×we3
R3=[0 -We3(3) We3(2); We3(3) 0 -We3(1);-We3(2) We3(1) 0]*sin(tt);
                                                  % Genera la matriz corespondiete a sin(θ)We3x

R=R1+R2+R3;                                       % Suma las matrices anteriores para calcula la matriz de rotacion
disp(R);

R3=R= [0.98079 -0.19509 0.00000 ;0.19509  0.98709 0.00000 ;0.00000  0.00000 1.00000] %Esta será el tensor de rotación segun el eje e3


3.3 Visualizacion de un sistema de puntos rotados con eje ω=e3 y angulo θ= π/16

Para poder visualizar el sistema de puntos rotado tenemos que multiplicar cada coordenada de nuestras particulas por la matriz de rotacion. De esta manera conseguimos sistema de puntos rotado. 

R3=R=[0.98079 -0.19509 0.00000 ;0.19509 0.98709 0.00000 ;0.00000 0.00000 1.00000]; %Componentes del tensor rotación de eje w3
i=(1:21); 
x=sin (pi*(i-1)/10);                                     %Coordenadas x de las particulas segun i
y=cos(pi*(i-1)/4);                                       %Coordenadas y de las particulas segun i
z=cos(pi*(i-1)/10);                                      %Coordenadas z de las particulas segun i
f=[x;y;z];                                               %Creamos una matriz cuyas filas seran las coordenadas xyz de los puntos
g=([x(i); y(i) ;z(i)]'*R3)';                             %La matriz g nos da las particulas rotadas segun R3
plot3(g(1,:),g(2,:),g(3,:),'o-','MarkerFaceColor','r')   %Dibuja las partículas rotadas
hold on
plot3(x,y,z,'x-','MarkerFaceColor','y')            %Dibujamos tambien los puntos de la curva
axis([-2,2,-2,2,-2,2])                             % Ejes fijados en la region [-2,2]*[-2,2]*[-2,2]
hold off


Visualización de los puntos de la curva y sus rotados segun el tensor rotación de eje e3 y angulo π/16















Los puntos azules son los iniciales, los puntos rojos son los rotados

3.4 Calculo del tensor de rotacion de un sistema de puntos con ejes ω = e1, ω = e2, ω = e1+e2+e3 y angulo θ= π/16

Para rotar el conjunto de puntos por estos ejes de giros usaremos el mismo codigo que en el apartado 3.1. Para esto calculamos de nuevo las matrices de rotacion para los ejes: ω = e1, ω = e2 y ω = e1+e2+e3. El procedimiento es el mismo, solo hay que cambiar el vector introducido por v=[1 0 0], v=[0 1 0] y v=[1 1 1] respectivamente.

Asi pues, para no repetir el codigo utilizado en apartados anteriores, escribiremos directamente el tensor rotación obtenido para cada eje, siendo R1 para [1 0 0], R2 para [0 1 0] y R4 para [1 1 1]. 

R1=[1 0 0; 0 0.980785 -0.19509; 0 0.19509 0.980785]                                  %Componentes del tensor de rotación de eje e1
R2=[0.980785 0 0.19509; 0 1 0; -0.19509 0 0.980785]                                  %Componentes del tensor de rotación de eje e2
R4=[0.98719 -0.106231 0.11904; 0.11904 0.98719 -0.106231; -0.106231 0.11904 0.98719] %Componentes del tensor de rotacion de eje e1+e2+e3


3.5 Visualizacion de un sistema de puntos rotados con ejes ω = e1, ω = e2, ω = e1+e2+e3 y angulo θ= π/16

Para visualizar los sistemas de puntos rotados tenemos que multiplicar cada coordenada de nuestras particulas por cada matriz de rotacion como hizimos en el apartado 3.2 . De esta manera conseguimos el sistema de puntos rotados para cada uno. Para evitar volver a escribir el codigo tres veces, mostraremos las graficas directamente.


(Los puntos azules son los iniciales, los puntos rojos son los rotados)

Visualización de los puntos de la curva y sus rotados segun el tensor rotación de eje e1 y angulo π/16
Visualización de los puntos de la curva y sus rotados segun el tensor rotación de eje e2 y angulo π/16
Visualización de los puntos de la curva y sus rotados segun el tensor rotación de eje e1+e2+e3 y angulo π/16























4 Velocidad de los puntos del sistema

4.1 Relación entre la velocidad de los puntos del sistema con su posición mediante un tensor antisimétrico A

Supongamos ahora que el sistema gira alrededor de un eje genérico con vector de dirección unitario ϖ de manera que la variación angular viene dada por la función θ(t) donde t ∈ [0, π] es el tiempo.

En primer lugar, vamos a comprobar analíticamente que la velocidad de los puntos del sistema está relacionada con su posición mediante un tensor antisimétrico A, es decir:

izquierda




Un tensor es antisimétrico si

izquierda



En (1) A es la velocidad angular. Para demostrar que la velocidad angular de una rotación es un tensor antisimético, vamos a derivar la siguiente relación:

izquierda



de donde se obtiene la matriz antisimétrica definida como:

izquierda



que coincide con la definición de tensor velocidad angular.


4.2 Demostracción que el vector axial asociado a A es θ(t)ϖ

La velocidad angular se define como el ángulo girado por unidad de tiempo.

izquierda




Para un objeto que gira alrededor de un eje, cada punto del objeto tiene la misma velocidad angular.

La forma matricial para representar la velocidad angular, puede ser deducida a partir de matrices de rotación.

Cualquier vector que gira alrededor de un eje con velocidad angular ϖ satisface:

izquierda




donde ω× se conoce como vector axial.


Podemos definir el tensor velocidad angular asociado con la velocidad angular ϖ como:

izquierda






Este tensor antisimétrico A(t) actúa como si ω× fuera un operador:

izquierda



Dada una matriz de rotación R(t) se puede obtener en cada instante el tensor velocidad angular A(t) como se muestra a continuación.


Se cumple que:

izquierda




Como la velocidad angular debe ser la misma para los tres vectores de un mismo sistema de referencia, si la matriz R(t) cuyas tres columnas son tres vectores unitarios simultáneamente perpendiculares, podemos escribir la relación:

izquierda




Y por tanto la velocidad angular se puede definir como:

izquierda




De donde se deduce que

izquierda



5 Visualización de vectores velocidad

Hallamos las componentes de la velocidad de las partículas al girar con velocidad angular ω=e3 con la formula que hemos demostrado en el apartado 4. Una vez encontramos estas componentes representaremos sus vectores aplicados en las coordenadas de cada punto. Para ello utilizamos: 

p=zeros(21,3);                                                  %Vector de posición de las partículas
for i=1:21
    p(i,:)=[cos(pi*(i-1)/10),cos(pi*(i-1)/4),cos(pi*(i-1)/10)]; %Definimos un vector con las componentes de los puntos de la curva
    w(i,:)=[0,0,1];                                             %Definimos el vector velocidad angular
end
v=cross(w,p);                                                   % Producto vectorial de los vectores posición y ω nos da el vector velocidad

hold on
figure(1)
quiver3(p(:,1),p(:,2),p(:,3),v(:,1),v(:,2),v(:,3))              %Representación tanto de la curva como de la velocidad
axis([-2,2,-2,2,-2,2])                                          %Ejes en los que definimos el grafico
plot3(p(:,1),p(:,2),p(:,3),'o','MarkerFaceColor','b') 
hold off


Visualización de los puntos de la curva y sus respectivos vectores velocidad, con el fin de no dificultar la visión de dichos vectores, hemos decidido no unir los puntos

















Los valores de las componentes de la velocidad van a ser distintos de 0 para el caso de e1 y e2, ya que en el caso de e3 sabemos que el producto vectorial de e3 por e3 es nulo. Por otro lado, el primer punto del sistema de partículas tiene como coordenadas (0,0,0), por tanto su velocidad es nula.

La gráfica obtenida sería la siguiente:

6 Momento angular de un sistema de partículas

El momento angular de un sistema de 21 partículas, que tienen un vector de desplazamiento ri, masa igual a mi, las cuales son todas idénticas y un vector de velocidad vi; este momento lineal se define como

izquierda





El momento lineal se puede expresar como:

izquierda



Para ello se supone que los puntos se mueven con una velocidad angular w y además se sustituye el valor de la velocidad vi por w x r.

izquierda











6.1 Momento angular

El momento angular es una magnitud física que bajo ciertas condiciones de simetría rotacional de los sistemas se mantiene constante con el tiempo a medida que el sistema evoluciona, que da lugar la ley de conservación conocida como Ley de conservación del momento angular. El momento angular de un conjunto de partículas es la suma de los momentos angulares de cada una:

Para calcular las componentes del tensor de inercia I lo haremos de la siguiente forma: 

L= zeros(3,3);
AUX= zeros(3,3);
m=10;                % masa de cada particula
for i=1:21;
    vpos=[cos(pi*(i-1)/10),cos(pi*(i-1)/4),cos(pi*(i-1)/10)];   %vector posición de las particulas
    for t=1:3;
        for h=1:3;
            if t==h;   % Para modificar la diagonal principal que ya calculamos directamente sus valores
                AUX(1,1)=(vpos(2)^2)+(vpos(3)^2);    % Elementos de la diagonal principal
                AUX(2,2)=(vpos(3)^2)+(vpos(1)^2);    
                AUX(3,3)=(vpos(2)^2)+(vpos(1)^2);   
            else
                AUX(t,h)=-(vpos(t)*vpos(h));         % Elementos t,h son el producto de la componentes t,h del vector de posicion
            end
        end
    end
     L=L+AUX;
end
L=m*L          % Multiplicamos la matriz por la masa de cada particula obteniendo el tensor momento de inercia.
w=[0 0 1];
Vangular= L*w'


7 Energía cinética del sistema

La energía cinética de un cuerpo es aquella energía que posee debido a su movimiento. Se define como el trabajo necesario para acelerar una masa determinada desde el reposo hasta alcanzar una velocidad . Una vez conseguida esta energía durante la aceleración, el cuerpo mantiene su energía cinética salvo que cambie su velocidad. Para un sistema de partículas, la energía cinética total se escribe como: 

centro

Sustituyendo en la fórmula anterior el valor de la velocidad que es: centro

Nos da:

centro

llegando así,a la expresión pedida.Hemos usado las propiedades del producto mixto y la regla de la expulsión. Para obtener la energía cinética del sistema utilizamos: 

L=[220 0 -110; 0 220 0; -110 0 220];
%velocidad angular
w=[0 0 1];                              
Ec=0.5*(w*L*w')
Ec= 110


8 Teorema de Steiner

El teorema de Steiner es un teorema usado en la determinación del momento de inercia de un sólido rígido sobre cualquier eje, dado el momento de inercia del objeto sobre el eje paralelo que pasa a través del centro de masa y de la distancia perpendicular (r) entre ejes. Se nos pide demostrar la siguiente fórmula, definiendo el vector 'a' como el vector que une el centro de masas G con el punto P:

centro

A continuación escribimos el codigo para usar el teorema de Steiner de modo que: 

m=10;                                    % masa de cada partícula
p=[-1.16633e-17 6.76707e-17 0.047619];   % Punto en el que queremos hallar el momento de inercia
G=[-1.16633e-17 6.76707e-17 0.047619];   % Coordenadas del centro de masas
t=G-p;                                   % t es el vector que uno G y p
B=m*(((t*t')*eye(3))-t'*t);              % Tensor el momento de Inercia para el vector t
x=zeros(1,21);
y=zeros(1,21);
z=zeros(1,21);
for i=1:21                               % Bucle para obtener las coordenadas de los puntos ya visto en apartados anteriores
 x(i)=sin(pi*(i-1)/10);                                                                     
 y(i)=cos(pi*(i-1)/4);                                       
 z(i)=cos(pi*(i-1)/10); 
end
vpos=[x' y' z'];                         %vpos sera el la matriz cuyas columnas son las componentes de los puntos
Ig=zeros(3,3);
for i=1:21
 Ig=m*(vpos(i,:))*(vpos(i,:)')*eye(3)-m*(vpos(i,:)')*(vpos(i,:))+Ig; % Tensor el momento de Inercia para el vector de posición
end
Ip=Ig+B                                  %Por ultimo obtenemos el momento de inercia en P utilizando el Teorema de Steiner
w=[0 0 1];
Ec=0.5*(w*Ig*w');                        %Demostración de que el la energia cinetica también se puede calcular de esta manera


9 Ejes principales de inercia

Llamamos ejes principales de inercia a los vectores propios de Ig. Todo sólido que gira libremente alrededor de uno de estos ejes no varía su orientación en el espacio. Cuando un sólido gira alrededor de uno de sus ejes principales, el momento angular L y la velocidad angular son vectores paralelos por estar ambos alineados con una dirección principal:

L=Ig·omega=λ·omega

Donde λ es una magnitud escalar que coincide con el momento de inercia correspondiente a dicho eje. En general, un cuerpo rígido tiene tres momentos principales de inercia diferentes. Y si dos ejes principales se corresponden a momentos principales de inercia diferentes, dichos ejes son perpendiculares.

Lo primero que hacemos es realizar el cálculo de autovalores y autovectores del Tensor Ig 

x=zeros(1,21);% Matrices de ceros para las coordenadas
y=zeros(1,21);
z=zeros(1,21);
G=[-1.16633e-17 6.76707e-17 0.047619];%Coordenadas del centro de masas
m=10; %Masa de cada partícula
for i=1:21 %Creamos de nuevo el bucle para obtener los puntos de la curva
x(i)=sin(pi*(i-1)/10);                                                                     
y(i)=cos(pi*(i-1)/4);                                       
z(i)=cos(pi*(i-1)/10); 
end
pos=[x' y' z']; %Generamos una matriz cuyas columnas seran las coordenadas x y z de los puntos
Ig=zeros(3,3); %Matriz de ceros para introducir Ig
for i=1:21 %Bucle del tensor diada que genera el termino de Ig cuya formula ya se ha deducido
Ig=m*(pos(i,:)')*(pos(i,:))+Ig;
end
[a,b]=eig(Ig);% Obtenemos los autovectores de la matriz Ig obteniendo las direcciones principales de Inercia
vector1=a(:,1)';
vector2=a(:,2)';  %Como se nos indica en el enunciado, estos autovectores seran las direcciones de los ejes
vector3=a(:,3)';
l=linspace(-1,0.5,2); %creamos un parametro landa que introducimos en las paramétricas.
E1x=G(1)+l*vector1(1); %A cada recta viene dada por el punto G centro de masas y la dirección del autovector de Ig.
E1y=G(2)+l*vector1(2);
E1z=G(3)+l*vector1(3);
E2x=G(1)+l*vector2(1);
E2y=G(2)+l*vector2(2);
E2z=G(3)+l*vector2(3);
E3x=G(1)+l*vector3(1);
E3y=G(2)+l*vector3(2);
E3z=G(3)+l*vector3(3);%De esta manera obtenemos las 3 rectas cada una con sus coordenadas.
hold on
plot3(x,y,z,"-*b") %Representación de la curva en azul. 
plot3(E1x,E1y,E1z,"r") %Representación del Eje 1 en rojo.
plot3(E2x,E2y,E2z,"g") %Representación del Eje 2 en verde.
plot3(E3x,E3y,E3z,"k") %Representación del Eje 3 en negro.
hold off

Quedando representado en el sistema de forma:

Visualización de los puntos de la curva y los tres ejes principales de inercia en verde rojo y negro, se puede observar que son perpendiculares por definición de eje principal de inercia

























10 Calculo del tensor de inercia de un solido

Nuestro solido esta definido de la siguiente manera: Su base es la placa plana con forma de anillo circular centrado en el origen y comprendido entre los radios 100 y 200 (centimetros). Luego tiene 20 centimetros de grosor y su densidad d depende del radio, es decir que d(ρ, θ, z) =|cosθ| + 1/ρ Kg/cm3.

10.1 Calculo de la masa total

Para calcular la masa total del solido tenemos que integrar su volumen por su densidad. Esto nos da una integral triple que resolveremos en cilindricas. Como hacemos el cambio a cilíndricas tenemos que multiplicar la integral por el jacobiano que es ρ por lo que la integral a calcular es ρcosθ+1. Al plantear dicha integral vemos que los extremos de integración para el termino cosθ van a ser [0,2*pi] por lo que ese factor quedaria anulado, resultando ser la integral 1 la que debemos calcular. Para Calcular esta integral en matlab hemos descompuesto la integral triple en una simple y otra doble. 

N=1000;                                       %Numero de puntos intermedios
ro1=100; ro2=200; tt1=0; tt2=2*pi; z1=0; z2=20;   %Extremos de intervalos expresados en cm
h=(z2-z1)/N; 
h1=(ro2-ro1)/N;  %Longitud de cada particion
h2=(tt2-tt1)/N;      
r=z1:h:z2; s=ro1:h1:ro2; t=tt1:h2:tt2;            %Vector de la partición para cada variable con intervalos cada h,h1,h2               
fun1=ones(N+1,1);                             %La funcion de la integral simple que es 1
fun2=ones(N+1,N+1);                        %Funcion de la integral doble que tambien es 1                                         
w=ones(N+1,1);                 
w(1)=1/2; w(N+1)=1/2;                      %Siguiendo el metodo del trapecio el primer y ultimo termino de w es 1/2
simple=h*w'*fun1;                                   %Resultado de la integral simple
doble=h1*h2*w'*fun2*w;                       %Resultado de la integral doble 
triple=simple*doble                                %Resultado de la integral triple (masa total)


La masa en kg es igual a 

  1.2566e+04

10.2 Calculo del centro de masas

Para calcular el centro de masas hemos tenido que calcular las coordenadas de Xg,Yg,Zg. Para esto hemos tenido que calcular tres integrales triples que a su vez hemos descompuesto en tres integrales simples. Las coordenadas del centro de masas vienen dadas por la integral centro


donde dm=densidad*dxdydz respectivamente asi pues nuestros integrandos en cilidriacas y habiendolos multiplicado por el jacobiano serán:

para Xg= ρ^2*cos^2(θ) teta+ro*cosθ

para Yg=ρ^2*senθ*cosθ+ρ*senθ

para Zg=ρ*cosθ+1 


N=1000;                                  %Numero de puntos intermedios
M=12566;                                %Masa total
ro1=100; ro2=200; teta1=0; teta2=2*pi; z1=0; z2=20;   %Extremos de intervalos en cm
h=(z2-z1)/N; h1=(ro2-ro1)/N; h2=(teta2-teta1)/N;         %Longitud de cada particion
r=z1:h:z2; s=ro1:h1:ro2; t=teta1:h2:teta2;                      %Coordenadas de la particion
[ro,teta]=meshgrid(s,t);
w=ones(N+1,1);                 
w(1)=1/2; w(N+1)=1/2;
%Funciones para la coordenada de Xg
fx1=ones(N+1,1);                     %Funcion de la integral simple de z (1)
fx2=ro;                                      %Funcion de la integral simple de ro (ro)
fx3=cos(teta);                           %Funcion de la integral simple de teta (cos(teta))     
fx4=ro^2;
fx5=cos(teta)*cos(teta);                                    
%Funciones para la coordenada de Yg
fy1=ones(N+1,1);                    %Funcion de la integral simple de z (1)
fy2=ro;                                     %Funcion de la integral simple de ro (ro)
fy3=sin(teta);                           %Funcion de la integral simple de teta (sin(teta))
fy4=ro^2;
fy5=sin(teta)*cos(teta);
%Funciones para la coordenada de Zg
fz1=r';                                       %Funcion de la integral simple de z (z)
fz2=ones(N+1,1);                     %Funcion de la integral simple de ro (1)
fz3=ones(N+1,1);                     %Funcion de la integral simple de teta (1)
fz4=ro;
fz5=cos(teta);
%Calculo de la coordenada de Xg
intx1=h*w'*fx1;                        %Resultado de la integral simple de z
intx2=h1*w'*fx2*w;                 %Resultado de la integral simple de ro
intx3=h2*w'*fx3*w;                 %Resultado de la integral simple de teta
intx4=h1*w'*fx4*w;
intx5=h2*w'*fx5*w;
xg1=intx1*intx2*intx3/M;           
xg2=intx4*intx5/M;
xG=xg1+xg2                       %Coordeanda de Xg   

%Calculo de la coordenada de Yg
inty1=h*w'*fy1;                        %Resultado de la integral simple de z
inty2=h1*w'*fy2*w;                 %Resultado de la integral simple de ro
inty3=h2*w'*fy3*w;                 %Resultado de la integral simple de teta
inty4=h1*w'*fy4*w;
inty5=h2*w'*fy5*w;
yg1=inty1*inty2*inty3/M;
yg2=inty4*inty5/M;
yG=yg1+yg2                         %Coordeanda de Yg
%Calculo de la coordenada de Zg
intz1=h*w'*fz1;                         %Resultado de la integral simple de z
intz2=h1*w'*fz2;                       %Resultado de la integral simple de ro
intz3=h2*w'*fz3;                       %Resultado de la integral simple de teta
intz4=h1*w'*fz4*w;
intz5=h2*w'*fz5*w;
zg1=intz1*intz2*intz3/M;
zg2=intz4*intz5/M;
zG=zg1+zg2                          %Coordeanda de Zg


El resultado nos da 

  Xg =-1.3719e-004
  Yg =-6.7179e-004
  Zg =10.0000

Como hemos usado la Formual del Trapecio en nuestras integrales el resultado es aproximado pero podriamos concluir que es 

  Xg=0
  Yg=0
  Zg=10

10.3 Calculo del tensor de inercia respecto al centro de masas

Una vez calculado la masa y el centro de masas es facil calcular el tensor de inercia. Este ultimo esta definido como I=M[(r·r)*1-r×r]. Siendo M la masa total; r el vector que une el centro con el centro de masas; r·r el producto escalar de r; 1 el tensor metrico y r×r el producto tensorial entre r y r . Para calcular este tensor hemos usado el siguiente codigo: 

Xg=0;                   %coordenada de centro de masas en X
Yg=0;                   %coordenada de centro de masas en Y
Zg=10;                  %coordenada de centro de masas en Z
M=12566;             %masa total
v=[Xg,Yg,Zg];     %vector que une el centro con el centro de masas
I1=(v*v')*eye(3); %Parte del tensor que coresponde con (r·r)*1
I2=kron(v,v');       %Parte del tensor que coresponde con r×r
I=(I1-I2)*M          %Tensor de Inercia


El tensor obtenido es el 

    1256600           0           0
          0     1256600           0
          0           0           0
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Movimiento_de_un_sistema_de_partículas_(Grupo_G9)&oldid=34032»