ENUNCIADO DEL EJEMPLO 26

Sólido formado por un disco homogéneo pesado de masa m=2 y radio a=10 unido rígidamente a una varilla de masa despreciable de longitud L=2a(2^1/2) perpendicularmente por el centro de ambos. Un extremo de la varilla está obligado a moverse según una recta vertical lisa, y el otro extremo se mueve sobre otra recta horizontal lisa que dista 2a de la primera. Se pide obtener las ecuaciones del movimiento.

> restart;

Empezaremos cargando los paquetes necesarios para realizar el ejercicio, entre ellos el mecapac3d, para el cual tendremos que indicar la ubicación de la librería en el disco duro.

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

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

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

> with(mecapac3d):

El sistema tiene dos grados de libertad, que asociaremos a las coordenadas generalizadas [theta,phi]. Theta representa el ángulo que forma la varilla con la vertical, y phi es el ángulo girado por el sistema alrededor de la dirección de la varilla.

> cg:=[theta,phi];

cg := [theta, phi]

Pasamos ahora a definir los elementos del sistema. La varilla y el disco los definiré dentro de un sistema de referencia móvil, cuyo eje z coincide con la varilla y los ejes x e y se encuentran en el plano del círculo, pero no acompañan al sólido en su movimiento de giro alrededor de la dirección de la varilla. Por ello, al definir la matriz de rotación del disco, habrá que indicar un giro alrededor del eje 3 (el eje z del sistema de referencia móvil) de valor phi. Definimos a continuación la varilla y el disco.

> v1:=[varilla,[0,0,0],rota(0,1),0,l]:

> d1:=[disco,[0,0,0],rota(phi,3),m,a]:

Para definir el subsistema emplearé un ángulo Psi asociado con el Theta a partir de una relación trigonométrica, de forma que la definición de la matriz de rotación del subsistema (sistema de referencia móvil) con respecto a los ejes fijos se exprese con mayor sencillez.

> psi:=arcsin(1/(2^(1/2)*sin(theta))):

> rot1:=rota(-theta,2):

> rot2:=rota(psi,3):

> rottot:=evalm(rot2&*rot1):

Definida la matriz de rotación del subsistema, definimos las coordenadas de su centro de coordenadas respecto del sistema fijo.

> xg:=[a*2^(1/2)*sin(theta)*cos(psi),a,a*cos(theta)*2^(1/2)]:

Como al definir la varilla habíamos indicado una longitud l que podría considerarse como aleatoria, hay que señalar que en realidad depende del valor de a.

> l:=2*a*2^(1/2):

Con todas las definiciones anteriores ya podemos definir el subsistema.

> s1:=[subsistema2,xg,rottot,[v1,d1]]:

Antes de definir el sistema general, incorporamos unos elementos gráficos para representar los ejes y la recta auxiliar por la que desliza la varilla en el plano coordenado (x,y)

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

> ejey:=[segmento,[0,0,0],[0,22,0],green]:

> ejez:=[segmento,[0,0,-20],[0,0,22],blue]:

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

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

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

> TX := [texto,[23,0,0],"X"]:

> rectaaux:=[segmento,[0,2*a,0],[20,2*a,0],red]:

> sistema:=[s1,ejex,ejey,ejez,rectaaux,TO,TX,TY,TZ];

sistema := [[subsistema2, [1/2*a*2^(1/2)*sin(theta)*(4-2/sin(theta)^2)^(1/2), a, a*cos(theta)*2^(1/2)], rottot, [[varilla, [0, 0, 0], matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]]), 0, 2*a*2^(1/2)], [disco...sistema := [[subsistema2, [1/2*a*2^(1/2)*sin(theta)*(4-2/sin(theta)^2)^(1/2), a, a*cos(theta)*2^(1/2)], rottot, [[varilla, [0, 0, 0], matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]]), 0, 2*a*2^(1/2)], [disco...sistema := [[subsistema2, [1/2*a*2^(1/2)*sin(theta)*(4-2/sin(theta)^2)^(1/2), a, a*cos(theta)*2^(1/2)], rottot, [[varilla, [0, 0, 0], matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]]), 0, 2*a*2^(1/2)], [disco...

Damos un valor al parámentro geométrico a para obtener una representación del sistema y así comprobar que se ha definido correctamente.

> a:= 10:

> fG([evalf(Pi/3),evalf(Pi/2)]);

[Plot]

Volvemos a dejar la constante a como un parámetro para obtener las energías y la Lagrangiana de una forma más genérica.

> a:='a':

> V:=fV(sistema):

> T:=fT(sistema):

> L:=simplify(T-V):

Ecuaciones del movimiento

> ecua := map(simplify,ec_lag());

ecua := [-1/4*m*a*(-22*a*diff(theta(t), `$`(t, 2))*cos(theta(t))^8*((-1+2*cos(theta(t))^2)/(-1+cos(theta(t))^2))^(1/2)-8*a*diff(phi(t), `$`(t, 2))*cos(theta(t))^8+16*g*sin(theta(t))*2^(1/2)*((-1+2*cos...ecua := [-1/4*m*a*(-22*a*diff(theta(t), `$`(t, 2))*cos(theta(t))^8*((-1+2*cos(theta(t))^2)/(-1+cos(theta(t))^2))^(1/2)-8*a*diff(phi(t), `$`(t, 2))*cos(theta(t))^8+16*g*sin(theta(t))*2^(1/2)*((-1+2*cos...ecua := [-1/4*m*a*(-22*a*diff(theta(t), `$`(t, 2))*cos(theta(t))^8*((-1+2*cos(theta(t))^2)/(-1+cos(theta(t))^2))^(1/2)-8*a*diff(phi(t), `$`(t, 2))*cos(theta(t))^8+16*g*sin(theta(t))*2^(1/2)*((-1+2*cos...ecua := [-1/4*m*a*(-22*a*diff(theta(t), `$`(t, 2))*cos(theta(t))^8*((-1+2*cos(theta(t))^2)/(-1+cos(theta(t))^2))^(1/2)-8*a*diff(phi(t), `$`(t, 2))*cos(theta(t))^8+16*g*sin(theta(t))*2^(1/2)*((-1+2*cos...ecua := [-1/4*m*a*(-22*a*diff(theta(t), `$`(t, 2))*cos(theta(t))^8*((-1+2*cos(theta(t))^2)/(-1+cos(theta(t))^2))^(1/2)-8*a*diff(phi(t), `$`(t, 2))*cos(theta(t))^8+16*g*sin(theta(t))*2^(1/2)*((-1+2*cos...ecua := [-1/4*m*a*(-22*a*diff(theta(t), `$`(t, 2))*cos(theta(t))^8*((-1+2*cos(theta(t))^2)/(-1+cos(theta(t))^2))^(1/2)-8*a*diff(phi(t), `$`(t, 2))*cos(theta(t))^8+16*g*sin(theta(t))*2^(1/2)*((-1+2*cos...ecua := [-1/4*m*a*(-22*a*diff(theta(t), `$`(t, 2))*cos(theta(t))^8*((-1+2*cos(theta(t))^2)/(-1+cos(theta(t))^2))^(1/2)-8*a*diff(phi(t), `$`(t, 2))*cos(theta(t))^8+16*g*sin(theta(t))*2^(1/2)*((-1+2*cos...

Damos valores a los parámetros geométricos y másicos, así como a la gravedad para realizar la integración numérica, para lo cuál se requiere indicar los valores iniciales de las coordenadas generalizadas y sus velocidades.

> m := 2 : g := 9.8 : a := 10 :

> res := fint([evalf(Pi/3),0.,evalf(Pi/2),0.]):

Hecha la integración, obtenemos las gráficas de las coordenadas generalizadas en función del movimiento.

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

[Plot]

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

[Plot]

Animación del movimiento.

> dibu3(2.5,20);

[Plot]

>