ENUNCIADO DEL EJEMPLO 2
    Un aro de masa uniforme M y radio R rueda en torno al eje vertical Z al que se encuentra unido por medio de una cuerda de masa despreciable.. Sobre dicho aro se encuentra insertada una particula de masa puntual m.

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

> restart:

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

Warning, the name changecoords has been redefined

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

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:=[phi,omega,theta]:

   Paso 2. Definición mediante variables de los elementos que forman el sistema mecánico. Estos son la masa puntual p1 y el aro con masa uniformemente repartida a1.

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

> r2:=rota(phi,2):

> r3:=rota(omega,3):

> r4:=rota(Pi/2,3):

> rot1:=evalm(r1 &* r2):rot2:=evalm(r4 &* rot1):rottot:=evalm(rot2 &* r3):

> cdg_a1:=[l*cos(phi),l*sin(phi),R]:

> a1:=[aro,cdg_a1,rottot,M,R]:

> cdg_p11:=sqrt(l^2+R^2*(cos(omega+theta))^2)*cos(phi+arctan(R*cos(omega+theta)/l)):

> cdg_p12:=sqrt(l^2+R^2*(cos(omega+theta))^2)*sin(phi+arctan(R*cos(omega+theta)/l)):

> cdg_p13:=R+R*sin(omega+theta):

> p1:=[punto,cdg_p11,cdg_p12,cdg_p13,m]:

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

> ejex:=[vector,[0,0,0],[5,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,[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:=[a1,p1,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*(l^2*sin(phi)^2*phi1^2+l^2*cos(phi)^2*phi1^2)+1/4*sin(omega)^2*phi1^2*M*R^2+1/4*phi1^2*cos(omega)^2*M*R^2+1/2*omega1^2*M*R^2+1/2*m*((-cos(phi+arctan(R*cos(omega+theta)/l))*R^2*cos(omega+the...T := 1/2*M*(l^2*sin(phi)^2*phi1^2+l^2*cos(phi)^2*phi1^2)+1/4*sin(omega)^2*phi1^2*M*R^2+1/4*phi1^2*cos(omega)^2*M*R^2+1/2*omega1^2*M*R^2+1/2*m*((-cos(phi+arctan(R*cos(omega+theta)/l))*R^2*cos(omega+the...T := 1/2*M*(l^2*sin(phi)^2*phi1^2+l^2*cos(phi)^2*phi1^2)+1/4*sin(omega)^2*phi1^2*M*R^2+1/4*phi1^2*cos(omega)^2*M*R^2+1/2*omega1^2*M*R^2+1/2*m*((-cos(phi+arctan(R*cos(omega+theta)/l))*R^2*cos(omega+the...

   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+m*g*(R+R*sin(omega+theta))

   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/2*omega1^2*M*R^2-M*g*R-m*g*R-m*g*R*sin(omega+theta)+1/2*l^2*M*phi1^2+1/4*phi1^2*M*R^2+1/2*l^2*m*phi1^2+1/2*m*phi1^2*R^2*cos(omega+theta)^2+m*R^2*omega1*theta1+1/2*m*R^2*omega1^2+1/2*m*R^2*theta...

   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.

> M:=12; l:=2.5; R:=1.5; m:=0.4; g:=9.8;

M := 12

l := 2.5

R := 1.5

m := .4

g := 9.8

A su vez comprobamos que las condiciones iniciales son las deseadas:

> fG([evalf(Pi/4),evalf(Pi/4),evalf(Pi/4)]);

[Plot]

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

> res:=fint([0,1,0,0,evalf(Pi/4),0]):

   Paso 11. Representación gráfica de las evoluciones temporales de la evolución temporal de phi y theta mediante odeplot.

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

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

> 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(6,70);

[Plot]

>