112exam.mws

Ejercicio 2, 4º parcial 13/06/2000

José Mª Goicolea, junio 2000

> restart:

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> T:=(1/2)*2*m*diff(x(t),t)^2
+(1/2)*m*(diff(x(t),t)^2+(R/sqrt(3))^2*diff(theta(t),t)^2+2*(R/sqrt(3))*diff(x(t),t)*diff(theta(t),t)*cos(theta(t)))
+(1/2)*((m/12)*R^2)*diff(theta(t),t)^2;

T := m*diff(x(t),t)^2+1/2*m*(diff(x(t),t)^2+1/3*R^2...

> V:=2*(1/2)*k*x(t)^2-m*g*(R/sqrt(3))*cos(theta(t));

V := k*x(t)^2-1/3*m*g*R*sqrt(3)*cos(theta(t))

> L:=T-V;

L := m*diff(x(t),t)^2+1/2*m*(diff(x(t),t)^2+1/3*R^2...
L := m*diff(x(t),t)^2+1/2*m*(diff(x(t),t)^2+1/3*R^2...

> L||1:=subs({diff(x(t),t)=x||1,diff(theta(t),t)=theta||1,x(t)=x,theta(t)=theta},L);

L1 := m*x1^2+1/2*m*(x1^2+1/3*R^2*theta1^2+2/3*R*sqr...

> t1:=diff(L||1,x||1);

t1 := 2*m*x1+1/2*m*(2*x1+2/3*R*sqrt(3)*theta1*cos(t...

> t2:=subs({x||1=diff(x(t),t),theta||1=diff(theta(t),t),x=x(t),theta=theta(t)},t1);

t2 := 2*m*diff(x(t),t)+1/2*m*(2*diff(x(t),t)+2/3*R*...

> t3:=diff(t2,t)-diff(L||1,x);

t3 := 2*m*diff(x(t),`$`(t,2))+1/2*m*(2*diff(x(t),`$...

> ecu||1:=subs({diff(x(t),t,t)=x||2,diff(theta(t),t,t)=theta||2,diff(x(t),t)=x||1,diff(theta(t),t)=theta||1,x(t)=x,theta(t)=theta},t3);

ecu1 := 2*m*x2+1/2*m*(2*x2+2/3*R*sqrt(3)*theta2*cos...

> ecu||1:=simplify(ecu||1);

ecu1 := 3*m*x2+1/3*m*R*sqrt(3)*theta2*cos(theta)-1/...

> t4:=diff(L||1,theta||1);

t4 := 1/2*m*(2/3*R^2*theta1+2/3*R*sqrt(3)*x1*cos(th...

> t5:=subs({x||1=diff(x(t),t),theta||1=diff(theta(t),t),x=x(t),theta=theta(t)},t4);

t5 := 1/2*m*(2/3*R^2*diff(theta(t),t)+2/3*R*sqrt(3)...

> t6:=diff(t5,t)-diff(L||1,theta);

t6 := 1/2*m*(2/3*R^2*diff(theta(t),`$`(t,2))+2/3*R*...
t6 := 1/2*m*(2/3*R^2*diff(theta(t),`$`(t,2))+2/3*R*...

> ecu||2:=subs({diff(x(t),t,t)=x||2,diff(theta(t),t,t)=theta||2,diff(x(t),t)=x||1,diff(theta(t),t)=theta||1,x(t)=x,theta(t)=theta},t6);

ecu2 := 1/2*m*(2/3*R^2*theta2+2/3*R*sqrt(3)*x2*cos(...

> ecu||2:=simplify(ecu||2);

ecu2 := 5/12*m*R^2*theta2+1/3*m*R*sqrt(3)*x2*cos(th...

> ### WARNING: persistent store makes one-argument readlib obsolete
readlib(mtaylor):

> ecu1l:=mtaylor(ecu||1,[x||1,theta||1,x,theta],2);

ecu1l := 3*m*x2+1/3*m*R*sqrt(3)*theta2+2*k*x

> ecu2l:=mtaylor(ecu||2,[x||1,theta||1,x,theta],2);

ecu2l := 5/12*m*R^2*theta2+1/3*m*R*sqrt(3)*x2+1/3*m...

> M:=matrix(
[[coeff(ecu1l,x||2),coeff(ecu1l,theta||2)],
[coeff(ecu2l,x||2),coeff(ecu2l,theta||2)]]
);

M := matrix([[3*m, 1/3*m*R*sqrt(3)], [1/3*m*R*sqrt(...

> K:=matrix(
[[coeff(ecu1l,x),coeff(ecu1l,theta)],
[coeff(ecu2l,x),coeff(ecu2l,theta)]]
);

K := matrix([[2*k, 0], [0, 1/3*m*g*R*sqrt(3)]])

> K||1:=map2(subs,{k=m*g,R=1},K);
M||1:=map2(subs,{k=m*g,R=1},M);

K1 := matrix([[2*m*g, 0], [0, 1/3*m*g*sqrt(3)]])

M1 := matrix([[3*m, 1/3*m*sqrt(3)], [1/3*m*sqrt(3),...

> ecu_car:=det(K||1-lambda*M||1);

ecu_car := 2/3*m^2*g^2*sqrt(3)-5/6*m^2*g*lambda-lam...

> s_lambda:=solve(ecu_car,lambda);

s_lambda := (5/11+6/11*sqrt(3)+1/11*sqrt(133-28*sqr...

> evalf(s_lambda);

2.234984836*g, .5636160460*g

> s_omega:=map(sqrt,[s_lambda]);

s_omega := [sqrt((5/11+6/11*sqrt(3)+1/11*sqrt(133-2...

> map(evalf,s_omega);

[1.494986567*sqrt(g), .7507436620*sqrt(g)]

> evalm(M||1&*inverse(K||1));

matrix([[3/2*1/g, 1/g], [1/6*sqrt(3)/g, 5/12*sqrt(3...

> eigenvals(K||1 &* inverse(M||1));

(5/11+6/11*sqrt(3)+1/11*sqrt(133-28*sqrt(3)))*g, (5...

> modos:=eigenvectors(K||1 &* inverse(M||1));

modos := [(5/11+6/11*sqrt(3)+1/11*sqrt(133-28*sqrt(...
modos := [(5/11+6/11*sqrt(3)+1/11*sqrt(133-28*sqrt(...

> evalf(modos);

[2.234984836*g, 1., {vector([-.9500558750, 1.])}], ...

> lambda[1]:=modos[1,1];omega[1]:=sqrt(lambda[1]);

lambda[1] := (5/11+6/11*sqrt(3)+1/11*sqrt(133-28*sq...

omega[1] := sqrt((5/11+6/11*sqrt(3)+1/11*sqrt(133-2...

> evalf(lambda[1]);evalf(omega[1]);

2.234984836*g

1.494986567*sqrt(g)

> lambda[2]:=modos[2,1];omega[2]:=sqrt(lambda[2]);

lambda[2] := (5/11+6/11*sqrt(3)-1/11*sqrt(133-28*sq...

omega[2] := sqrt((5/11+6/11*sqrt(3)-1/11*sqrt(133-2...

> evalf(lambda[2]);evalf(omega[2]);

.5636160460*g

.7507436620*sqrt(g)

> modo[1]:=simplify(modos[1,3][1]);

modo[1] := vector([-5/4+3/2*sqrt(3)-1/4*sqrt(133-28...

> map(evalf,modo[1]);

vector([-.950055874, 1.])

> modo[2]:=simplify(modos[2,3][1]);

modo[2] := vector([-5/4+3/2*sqrt(3)+1/4*sqrt(133-28...

> map(evalf,modo[2]);

vector([3.646208298, 1.])

>