El péndulo esférico

Consideramos el movimiento de un péndulo de longitud l que se puede mover en todas las direcciones del espacio bajo la acción de su propio peso. En un instante dado t, la partícula dista l del centro de oscilación y hace un ángulo φ con el eje X y un ángulo θ con el eje Z, tal como se muestra en la figura. Describimos el movimiento de la partícula en coordenadas esféricas (r, φ, θ) con el radio r=l constante.

La posición de la partícula es

{ x=lsinθcosφ y=lsinθsinφ z=lcosθ

Derivando con respecto del tiempo obtenemos las componentes de la velocidad

{ dx dt =lcosθ dθ dt cosφlsinθsinφ dφ dt dy dt =lcosθ dθ dt sinφ+lsinθcosφ dφ dt dz dt =lsinθ dθ dt

La energía cinética es

T= 1 2 m( ( dx dt ) 2 + ( dy dt ) 2 + ( dz dt ) 2 )= 1 2 m l 2 ( ( dθ dt ) 2 + sin 2 θ ( dφ dt ) 2 )

La energía potencial es V=-mgz=-mglcosθ

La lagrangiana L=T-V es

L= 1 2 m l 2 ( ( dθ dt ) 2 + sin 2 θ ( dφ dt ) 2 )+mglcosθ L= 1 2 m l 2 ( θ ˙ 2 + sin 2 θ· φ ˙ 2 )+mglcosθ

Ecuaciones del movimiento

La primera ecuación del movimiento es

d dt ( L θ ˙ ) L θ =0 l d 2 θ d t 2 lsinθcosθ ( dφ dt ) 2 +gsinθ=0

La segunda ecuación del movimiento es

d dt ( L φ ˙ ) L φ =0 d dt ( m l 2 sin 2 θ dφ dt )=0

Principios de conservación

Como la lagrangiana L es independiente de φ, obtenemos una cantidad h que se conserva

sin 2 θ dφ dt =h

Comprobamos que h está relacionado con la componente Z del momento angular L = r ×m v

L z =m( x dy dt y dx dt )=m l 2 sin 2 θ dφ dt

Llamamos ω 0 2 =g/l a la frecuencia natural de la oscilación de un péndulo simple. La primera ecuación se convierte en

d 2 θ d t 2 h 2 sin 3 θ cosθ+ ω 0 2 sinθ=0

que a su vez es la derivada con respecto del tiempo de una cantidad que se conserva

d dt ( 1 2 ( dθ dt ) 2 + h 2 2 sin 2 θ ω 0 2 cosθ )=0

La conservación de la energía ε se expresa

1 2 ( dθ dt ) 2 + h 2 2 sin 2 θ ω 0 2 cosθ=ε

Donde V(θ) es la energía potencial efectiva.

V(θ)= h 2 2 sin 2 θ ω 0 2 cosθ

Para un valor dado de ε el movimiento tiene lugar entre los ángulos θ1 y θ2 dados por las soluciones de la ecuación

ε= h 2 2 sin 2 θ ω 0 2 cosθ

R=1; %longitud del péndulo
w0=sqrt(9.8/R); %frecuencia natural de oscilación
h2=21.1; % cuadrado del momento angular, Lz
energia=37.3; %energía
f=@(theta) h2/(2*sin(theta)^2)-w0^2*cos(theta);
fplot(f,[pi/12,pi-pi/12],'color','r')
line([0,pi],[energia, energia],'color','k')
g=@(theta) f(theta)-energia;

%raíces
theta1=fzero(g,1);
theta2=fzero(g,2);
line([theta1,theta1],[0,energia], 'lineStyle','--')
line([theta2,theta2],[0,energia],'lineStyle','--')
ylim([0,40])
grid on
xlabel('\theta')
ylabel('V(\theta)')
title('Energía potencial')

Trayectoria de la partícula

Para representar la trayectoria seguida por la partícula, partimos de la posición angular θ=θ1 en el que se cumple que dθ/dt=0, establecemos φ=0. Determinamos la constante h a partir de los datos de θ1 y θ2. Igualando V(θ1)=V(θ2) obtenemos

ε= h 2 2 sin 2 θ 1 ω 0 2 cos θ 1 = h 2 2 sin 2 θ 2 ω 0 2 cos θ 2 h 2 = ω 0 2 sin 2 θ 1 sin 2 θ 2 cos( θ 1 + θ 2 2 )cos( θ 1 θ 2 2 )

Integramos las dos ecuaciones diferenciales con las condiciones iniciales señaladas para obtener la trayectoria θ=θ(φ) en la superficie esférica de radio l

d 2 θ d t 2 h 2 sin 3 θ cosθ+ ω 0 2 sinθ=0 sin 2 θ dφ dt =h

R=1; %longitud del péndulo
w0=sqrt(9.8/R); %frecuencia natural de oscilación

theta1=pi/6;
theta2=pi/3; %angulos mínimo y máximo
h2=w0^2*sin(theta1)^2*sin(theta2)^2/(cos((theta1+theta2)/2)
*cos((theta1-theta2)/2));
x0=[theta1,0,0]; %condiciones iniciales
tspan=[0,5];
% x(1)=theta, x(2)=dtheta/dt,x(3)=phi
fg=@(t,x)[x(2); h2*cos(x(1))/(sin(x(1))^3)-w0^2*sin(x(1));
 sqrt(h2)/sin(x(1))^2];
[t,x]=ode45(fg,tspan,x0);
xp=R*sin(x(:,1)).*cos(x(:,3));
yp=R*sin(x(:,1)).*sin(x(:,3));
zp=R*cos(x(:,1));

%esfera
phi=linspace(0,pi,30);
theta=linspace(0,2*pi,40);
[phi,theta]=meshgrid(phi,theta);
x=R*sin(phi).*cos(theta);
y=R*sin(phi).*sin(theta);
z=R*cos(phi);
h1=mesh(x,y,z);
set(h1,'EdgeColor',[0.6,0.6,0.6],'EdgeAlpha',0.5,'FaceAlpha',0.5)
axis equal

%trayectoria
h1=line(xp,yp,zp);
set(h1,'Color',[.7,0,0],'LineWidth',1.5)
set(gca,'zDir','reverse')
view (-40,20)
hold off
axis equal
grid on
xlabel('x')
ylabel('y')
zlabel('z')
title('Movimiento en la superficie de una esfera')

Si φ=φ0=constante, dφ/dt=0, h=0

d 2 θ d t 2 + ω 0 2 sinθ=0

Obtenemos la ecuación diferencial del péndulo simple que oscila en el plano φ0

Trayectorias casi circulares

Un péndulo cónico gira alrededor del eje Z con velocidad angular dφ/dt constante, haciendo un ángulo θ0 con la vertical. La partícula describirá una trayectoria circular de radio lsinθ0 perpendicular al eje Z.

Si θ=θ0=constante, dφ/dt=constante.

h 2 = ω 0 2 sin 4 θ 0 cos θ 0 dφ dt = h sin 2 θ 0 = ω 0 cos θ 0 = g lcos θ 0

Supongamos ahora que θ=θ0+δ(t). Donde δ es una cantidad pequeña frente a θ0

Para determinar δ(t), en la ecuación diferencial

d 2 θ d t 2 h 2 sin 3 θ cosθ+ ω 0 2 sinθ=0 d 2 θ d t 2 + ω 0 2 ( sin 4 θ 0 cos θ 0 cosθ sin 3 θ +sinθ )=0

aproximamos el término entre paréntesis utilizando Math Symbolic.

>> syms x x0;
>> y=sin(x)-sin(x0)^4*cos(x)/(cos(x0)*sin(x)^3);
>> z=taylor(y,x,x0,'Order',2);
>> simplify(z)
((x - x0)*(3*cos(x0)^2 + 1))/cos(x0)

La ecuación diferencial se transforma en otra mucho más simple

d 2 δ d t 2 + ω 0 2 1+3 cos 2 θ 0 cos θ 0 δ=0

que corresponde a una oscilación con frecuencias angular

Ω= ω 0 1+3 cos 2 θ 0 cos θ 0

En el script anterior cambiamos las líneas

theta1=pi/6;
theta2=pi/4; %ángulos mínimo y máximo
....

para obtener mediante integración numérica la trayectoria de la partícula, y compararla con la trayectoria que se obtiene al suponer una trayectoria casi circular

R=1; %longitud del péndulo
w0=sqrt(9.8/R); %frecuencia natural de oscilación

theta1=pi/4;
theta2=pi/3; %angulos mínimo y máximo
W=w0*sqrt((1+3*cos(theta1)^2)/cos(theta1));
w=w0/sqrt(cos(theta1));
% x(1)=theta, x(2)=dtheta/dt,x(3)=phi
t=0:0.01:5;
theta=(theta1+theta2)/2+(theta2-theta1)*sin(W*t)/2;
phi=w*t;
xp=R*sin(theta).*cos(phi);
yp=R*sin(theta).*sin(phi);
zp=R*cos(theta);

%esfera
phi=linspace(0,pi,30);
theta=linspace(0,2*pi,40);
[phi,theta]=meshgrid(phi,theta);
x=R*sin(phi).*cos(theta);
y=R*sin(phi).*sin(theta);
z=R*cos(phi);
h1=mesh(x,y,z);
set(h1,'EdgeColor',[0.6,0.6,0.6],'EdgeAlpha',0.5,'FaceAlpha',0.5)
axis equal

%trayectoria
h1=line(xp,yp,zp);
set(h1,'Color',[.7,0,0],'LineWidth',1.5)
set(gca,'zDir','reverse')
view (-40,20)
hold off
axis equal
grid on
xlabel('x')
ylabel('y')
zlabel('z')
title('Movimiento en la superficie de una esfera')