ENUNCIADO DEL EJEMPLO 5

    Un disco homogeneo de radio R y masa M esta contenido en un plano vertical que corresponde con el plano YOZ estandole permitida la rodadura sin deslizamiento sobre el eje OY.

Al centro de este disco se unen otros dos elementos. El primero de ellos es una varilla homogenea de masa m y longitud l que actua como un pendulo tridimensional (se supone que no existe interaccion posible entre el disco y la varilla excepto en su punto de union). El segundo es un muelle lineal que une el centro del disco con el origen de coordenadas.

   Paso 0. Reiniciación de las variables del sistema y llamada a los paquetes linalg, plots y plottools.

> restart:

> with(linalg):with(plots):with(plottools):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

> libname:="C:/",libname;

libname :=

> with(mecapac3d):

   Paso 1. Definimos las coordenadas generalizadas del sistema en una lista que se denominará cg.

> cg := [theta,phi,s];  

cg := [theta, phi, s]

   Paso 2. Definición mediante variables de los elementos que forman el sistema mecánico. Es decir, la varilla, el disco y el muelle.

Comenzamos por la varilla de masa m y longitud l.

> xg := [l*cos(phi)*cos(theta)/2,l*cos(phi)*sin(theta)/2,-l*sin(phi)/2];

xg := [1/2*l*cos(phi)*cos(theta), 1/2*l*cos(phi)*sin(theta), -1/2*l*sin(phi)]

> r1 := rota(-Pi/2-phi,1);

r1 := matrix([[1, 0, 0], [0, -sin(phi), cos(phi)], [0, -cos(phi), -sin(phi)]])

> r2 := rota(-Pi/2+theta,3);

r2 := matrix([[sin(theta), cos(theta), 0], [-cos(theta), sin(theta), 0], [0, 0, 1]])

> rtot := evalm(r2 &*r1):

> va := [varilla,xg,rtot,m,l];

va := [varilla, [1/2*l*cos(phi)*cos(theta), 1/2*l*cos(phi)*sin(theta), -1/2*l*sin(phi)], rtot, m, l]

Definimos el subsistema formado por esta varilla:

> sist:=[subsistema2,[0.,s,R],rota(0,1),[va]];

sist := [subsistema2, [0., s, R], matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]]), [[varilla, [1/2*l*cos(phi)*cos(theta), 1/2*l*cos(phi)*sin(theta), -1/2*l*sin(phi)], rtot, m, l]]]

Y a continuacion definimos el disco con masa M y radio R.

> dis:=[disco,[0.,s,R],evalm(rota(-s/R,1)&*rota(Pi/2,2)),M,R];

dis := [disco, [0., s, R], matrix([[0, 0, 1], [-sin(s/R), cos(s/R), 0], [-cos(s/R), -sin(s/R), 0]]), M, R]

Lo mismo para el muelle indicando su rigidez k y una masa despreciable 0.

> m1:=[muelle,[0.,0.,R],[0.,s,R],k,0];

m1 := [muelle, [0., 0., R], [0., s, R], k, 0]

   Paso 3. Definición de los elementos gráficos que definiran nuestro sistema de ejes.

> ejex:=[vector,[0,0,0],[2.5,0,0],red]:

> ejey:=[vector,[0,-5,0],[0,5,0],green]:

> ejez:=[vector,[0,0,0],[0,0,5],blue]:

> TO := [texto,[0,0,-1],"O"]:

> TX := [texto,[2.5,0,-1],"X"]:

> TY := [texto,[0,5,-1],"Y"]:

> TZ := [texto,[0,0,6],"Z"]:

   Paso 4. Definición de la variable sistema que agrupa en una lista todos los elementos anteriores.

> sistema:=[sist,m1,dis,ejex,ejey,ejez,TO,TX,TY,TZ]:

   Paso 5. Obtención de la energía cinética del sistema mediante fT asignándola a la variable T.

> T:=fT(sistema);

T := 1/2*m*((-1/2*l*sin(phi)*phi1*cos(theta)-1/2*l*cos(phi)*sin(theta)*theta1)^2+(s1-1/2*l*sin(phi)*phi1*sin(theta)+1/2*l*cos(phi)*cos(theta)*theta1)^2+1/4*l^2*cos(phi)^2*phi1^2)+3/4*s1^2*M

   Paso 6. Obtención de la energía potencial del sistema mediante fV asignándola a la variable V.

> V:=fV(sistema);

V := m*g*(R-1/2*l*sin(phi))+1/2*k*s^2+M*g*R

   Paso 7. Obtención de la lagrangiana como diferencia de energías entre la energía cinética y la potencial.

> L:=simplify(T-V);

L := 1/8*m*l^2*cos(phi)^2*theta1^2+1/2*m*s1^2-1/2*m*s1*l*sin(phi)*phi1*sin(theta)+1/2*m*s1*l*cos(phi)*cos(theta)*theta1+1/8*m*l^2*phi1^2+3/4*s1^2*M-m*g*R+1/2*m*g*l*sin(phi)-1/2*k*s^2-M*g*R

   Paso 8. Obtención de las ecuaciones de lagrange para las dos coordenadas generalizadas mediante el operador Ec_lag

> ecua:=ec_lag():

   Paso 9. Asignación de valores numéricos a los parámetros que queden sun asignar para poder proceder a la integración numérica.

> g:=9.8;k:=100;l:=2;m:=10;M:=50;R:=1;

g := 9.8

k := 100

l := 2

m := 10

M := 50

R := 1

   Paso 10. Integración numérica del problema mediante la función fint asignando el resultado a la variable res.

> res := fint([0.,0.,0.,0.,5,0.]):

   Paso 11. Procedemos a realizar una animación del movimiento del conjunto por medio de la función dibu3.

> dibu3(5,75);

[Plot]

>