ENUNCIADO DEL EJEMPLO 9
    

   Una particula de masa concentrada m se haya insertada en un semiaro homogeneo de masa M y radio R.

Los extremos del semiaro se encuentran ligados al eje OX 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:

> with(mecapac3d):

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

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

cg := [x, theta, phi]

> R:=1;

R := 1

   Paso 2. Definición mediante variables de los elementos que forman el sistema mecánico.

El sistema consta de un semiaro pesado que se desplaza articulado por sus extremos en el eje X, sobre el que puede girar. Ademas se encuentra ensartada una particula pesada.

Defino primero el centro de masa del semiaro, su rotacion y un subsistema de ejes ligados a el.

> xgsemia:=[x,2*R/Pi*cos(theta),-2*R/Pi*sin(theta)];

xgsemia := [x, 2*cos(theta)/Pi, -2*sin(theta)/Pi]

> rotsemia:=rota(-theta,1);

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

> submovil:=[subsistema2,[x,0,0],rotsemia,[part]];

submovil := [subsistema2, [x, 0, 0], rotsemia, [part]]

Ahora defino el semiaro y la particula.

> semia:=[semiaro,xgsemia,rotsemia,M,R];

semia := [semiaro, [x, 2*cos(theta)/Pi, -2*sin(theta)/Pi], rotsemia, M, 1]

> part:=[punto,R*cos(phi),R*sin(phi),0,m];

part := [punto, cos(phi), sin(phi), 0, m]

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

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

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

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

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

> TX := [texto,[20,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:=[semia,ejex,ejey,ejez,TO,TX,TY,TZ,submovil];

sistema := [[semiaro, [x, 2*cos(theta)/Pi, -2*sin(theta)/Pi], rotsemia, M, 1], [vector, [0, 0, 0], [20, 0, 0], red], [vector, [0, 0, 0], [0, 5, 0], green], [vector, [0, 0, 0], [0, 0, 5], blue], [texto...sistema := [[semiaro, [x, 2*cos(theta)/Pi, -2*sin(theta)/Pi], rotsemia, M, 1], [vector, [0, 0, 0], [20, 0, 0], red], [vector, [0, 0, 0], [0, 5, 0], green], [vector, [0, 0, 0], [0, 0, 5], blue], [texto...

   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*(x1^2+4*sin(theta)^2*theta1^2/Pi^2+4*cos(theta)^2*theta1^2/Pi^2)+1/4*theta1^2*M*(Pi^2-8)/Pi^2+1/2*m*((x1-sin(phi)*phi1)^2+(-sin(theta)*theta1*sin(phi)+cos(theta)*cos(phi)*phi1)^2+(-cos(thet...

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

> V:=fV(sistema);

V := -2*M*g*sin(theta)/Pi-m*g*sin(theta)*sin(phi)

   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/4*(2*M*x1^2*Pi+theta1^2*M*Pi+2*m*Pi*x1^2-4*m*Pi*x1*sin(phi)*phi1+2*m*Pi*phi1^2+2*m*Pi*theta1^2-2*m*Pi*theta1^2*cos(phi)^2+8*M*g*sin(theta)+4*m*g*sin(theta)*sin(phi)*Pi)/Pi

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

> ecua:=ec_lag();

ecua := [1/4*(4*M*diff(x(t), `$`(t, 2))*Pi+4*m*Pi*diff(x(t), `$`(t, 2))-4*m*Pi*cos(phi(t))*diff(phi(t), t)^2-4*m*Pi*sin(phi(t))*diff(phi(t), `$`(t, 2)))/Pi, 1/4*(2*diff(theta(t), `$`(t, 2))*M*Pi+4*m*P...ecua := [1/4*(4*M*diff(x(t), `$`(t, 2))*Pi+4*m*Pi*diff(x(t), `$`(t, 2))-4*m*Pi*cos(phi(t))*diff(phi(t), t)^2-4*m*Pi*sin(phi(t))*diff(phi(t), `$`(t, 2)))/Pi, 1/4*(2*diff(theta(t), `$`(t, 2))*M*Pi+4*m*P...ecua := [1/4*(4*M*diff(x(t), `$`(t, 2))*Pi+4*m*Pi*diff(x(t), `$`(t, 2))-4*m*Pi*cos(phi(t))*diff(phi(t), t)^2-4*m*Pi*sin(phi(t))*diff(phi(t), `$`(t, 2)))/Pi, 1/4*(2*diff(theta(t), `$`(t, 2))*M*Pi+4*m*P...

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

> g:=9.8;M:=2;m:=1;

g := 9.8

M := 2

m := 1

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

> res:=fint([0,3,0,0.1,Pi/4,0.1]):

   Paso 11. Representación gráfica de las evoluciones temporales x, theta y phi mediante odeplot.

> p1:=odeplot(res,[t,x(t)],0..4,color=red):

> p2:=odeplot(res,[t,theta(t)],0..4,color=green):

> p3:=odeplot(res,[t,phi(t)],0..4,color=blue):

> display({p1,p2,p3});

[Plot]

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

> dibu3(4,70);

[Plot]

>