Se dispone de un semiaro de masa M y radio R . En su punto más bajo hay unido un muelle de constante k y longitud natural l, estando el  otro extremo del muelle en el suelo(en el origen de coordenadas).El centro de gravedad del semiaro se mueve en todo momento segun el eje z. El semiaro  puede girar alrededor del eje z. Sobre el semiaro, y en su mismo plano, se mueve un disco de masa m y radio r,que rueda sin deslizar. Obtener las ecuaciones del movimiento del sistema.

>

restart:

Cargamos los paquetes de Maple que vamos a necesitar.

>	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

Indicamos el directorio en el que se encuentra la librería mecapac3d y lo cargamos.

>	libname:="C:\",libname:

>	with(mecapac3d):

El sistema formado por el semiaro y disco tiene tres grados de libertad. Los identificamos con las coordenadas generalizadas theta, phi y z, que representan, respectivamente, el ángulo que forma con la vertical el radiovector que une el centro del semiaro con el centro de gravedad del disco, el giro del semiaro alrededor de la vertical, y la coordenada z del punto más bajo del semiaro. Con estas tres coordenadas generalizadas queda definida la posición de todos los elementos del sistema en cualquier instante de tiempo.

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

Definimos ahora el semiaro, obteniendo primero la posición del centro de gravedad y la matriz de rotación entre los ejes fijos y los ejes propios. Para obtener la matriz de rotación realizamos el producto por la izquierda, ya que los sucesivos giros que se realizan al semiaro para llevarlo a su posición final son giros absolutos, referidos al sistema de referencia de ejes fijos.

>	xgvar1:=[0,0,z+R-2*R/Pi]:

>	rotp1:=rota(-Pi/2,3):

>	rotp2:=rota(Pi/2,2):

>	rotp4:=rota(phi,3):

>	rotp:=evalm(rotp4&*rotp2&* rotp1):

>	m1:=[semiaro,xgvar1,rotp,M,R];

De la misma forma, definimos el disco, indicando las coordenadas de su centro de gravedad en función de las coordenadas generalizadas, así como la matriz de rotación, teniendo en cuenta que el disco rueda sin deslizar.

>	xgvar2:=[-(R-r)*sin(theta)*sin(phi),(R-r)*sin(theta)*cos(phi),z+R-(R-r)*cos(theta)]:

>	rotp3:=rota(Pi/2,2):

>	rotp6:=rota(phi,3):

>	rotp7:=rota(-R/r*theta,1):

>	rotdisco:=evalm(rotp6&*rotp7&*rotp3):

>	m2:=[disco,xgvar2,rotdisco,m,r];

A continuación definimos el muelle, indicando sus dos extremos, su constante y su longitud natural.

>	x1:=[0.,0.,0.]:

>	x2:=[0.,0.,z]:

>	c1:=[muelle,x1,x2,k,l];

También incorporamos al sistema unos elementos gráficos para visualizarlo mejor. Los elementos son segmentos para representar los ejes, dos ángulos que simbolicen las coordenadas generalizadas phi y theta, y textos que nombren los ejes.

>	ejex:=[segmento,[0,0,0],[15,0,0],red]:

>	ejey:=[segmento,[0,0,0],[0,15,0],blue]:

>	ejez:=[segmento,[0,0,0],[0,0,15],green]:

>	antheta:=[angulo,[0,0,z],[0,0,z+R],[-(R-r)*sin(theta)*sin(phi),(R-r)*sin(theta)*cos(phi),z+R-(R-r)*cos(theta)],R/2]:

>	anphi:=[angulo,[0,R,z+R],[0,0,z+R],[-R*sin(phi),R*cos(phi),z+R],R/2]:

>	textox:=[texto,[15,2,0],"X"]:

>	textoy:=[texto,[2,15,0],"Y"]:

>	textoz:=[texto,[2,2,15],"Z"]:

A continuación, definimos el sistema, indicando todos los elementos que lo conforman.

>	sistema:=[m1,m2,c1,ejex,ejey,ejez,antheta,anphi,textox,textoy,textoz];

Definimos las componentes de la Lagrangiana para después obtener las ecuaciones de Lagrange o del movimiento.

>	T:=fT(sistema):

>	V:=fV(sistema):

>	L:=simplify(T-V);

Ecuaciones del movimiento

>	ecua:=simplify(ec_lag());

Damos valores a los parámetros que anteriormente habíamos definido por nombre, para obtener una representación del sistema, comprobar que el sistema se ha definido correctamente y poder realizar la integración numérica.

>	M:=2:R:=10:m:=1:l:=10:g:=9.8:k:=10:r:=2.5:

>	fG([evalf(Pi/3),evalf(Pi/8),5.]);

Integración numérica, indicando valores iniciales y velocidades iniciales.

>	res:=fint([evalf(Pi/3),0.3,evalf(Pi/8),0.6,5.,0.1]);

Animación del movimiento durante 5 segundos con 40 imágenes.

>	dibu3(4.8,40);

A continuación representamos las gráficas de las coordenadas generalizadas en el tiempo, así como la relación entre las distintas coordenadas generalizadas a lo largo del tiempo.

>	odeplot(res,[t,theta(t)],t=0..5,color=red);

>	odeplot(res,[t,phi(t)],t=0..5,color=green);

>	odeplot(res,[t,z(t)],t=0..5,color=blue);

>	odeplot(res,[z(t),phi(t)],0..4);

>	odeplot(res,[theta(t),phi(t)],0..4);

>	odeplot(res,[z(t),theta(t)],0..4);

>