El péndulo doble

El péndulo doble como se muestra en la figura, está formado por dos péndulos simples de longitudes l1 y l2, de los que cuelgan partículas de masas m1 y m2. En un instante determinado t, los hilos inextensibles forman ángulos θ1 y θ2 con la vertical.

Energía

La energía del sistema es la suma de la energía potencial y de la energía cinética de las dos partículas.

Energía potencial

Situamos el nivel cero de energía potencial en el punto de suspensión del primer péndulo. La energía potencial es

V=-m1·gl1cos θ1-m2g(l1cos θ1+l2cos θ2)=-(m1+m2)·gl1cos θ1-m2gl2cos θ2

Energía cinética

El módulo de la velocidad del primer péndulo es v1=l1·dθ1/dt (velocidad angular por el radio de la circunferencia que describe), sus componentes son

v 1 =( l 1 d θ 1 dt cos θ 1 i+ l 1 d θ 1 dt sin θ 1 j )

El módulo de la velocidad de la segunda partícula en un sistema de referencia que se mueve con la velocidad v1 de la primera partícula es v2=l2·dθ2/dt. Sus componentes son

v 2 =( l 2 d θ 2 dt cos θ 2 i+ l 2 d θ 2 dt sin θ 2 j )

La velocidad de la segunda partícula respecto al sistema de referencia inercial situado en O es la suma vectorial de ambas velocidades (véase la figura de la derecha)

v 2 + v 1 =( l 2 d θ 2 dt cos θ 2 + l 1 d θ 1 dt cos θ 1 )i+( l 2 d θ 2 dt sin θ 2 + l 1 d θ 1 dt sin θ 1 )j

Calculamos los módulos de las velocidades de las dos partículas. La energía cinética del sistema es

T = 1 2 m 1 l 1 2 ( d θ 1 dt ) 2 + 1 2 m 2 l 1 2 ( d θ 1 dt ) 2 + 1 2 m 2 l 2 2 ( d θ 2 dt ) 2 + m 2 l 1 l 2 ( d θ 1 dt )( d θ 2 dt )cos( θ 1 θ 2 )

Ecuaciones del movimiento

Las ecuaciones del movimiento de Lagrange nos llevan al sistema de dos ecuaciones diferenciales de segundo orden.La lagrangiana L=T-V

L= 1 2 ( m 1 + m 2 ) l 1 2 ( d θ 1 dt ) 2 + 1 2 m 2 l 2 2 ( d θ 2 dt ) 2 + m 2 l 1 l 2 ( d θ 1 dt )( d θ 2 dt )cos( θ 1 θ 2 ) +( m 1 + m 2 )g l 1 cos θ 1 + m 2 g l 2 cos θ 2

Las primera ecuación del movimiento se obtiene

d dt L θ ˙ 1 L θ 1 =0 θ ˙ 1 = d θ 1 dt L θ ˙ 1 =( m 1 + m 2 ) l 1 2 ( d θ 1 dt )+ m 2 l 1 l 2 ( d θ 2 dt )cos( θ 1 θ 2 ) L θ 1 = m 2 l 1 l 2 ( d θ 1 dt )( d θ 2 dt )sin( θ 1 θ 2 )( m 1 + m 2 )g l 1 sin θ 1 d dt ( L θ ˙ 1 )=( m 1 + m 2 ) l 1 2 ( d 2 θ 1 d t 2 )+ m 2 l 1 l 2 ( d 2 θ 2 d t 2 )cos( θ 1 θ 2 ) m 2 l 1 l 2 ( d θ 2 dt )sin( θ 1 θ 2 )( d θ 1 dt + d θ 2 dt ) ( m 1 + m 2 ) l 1 2 ( d 2 θ 1 d t 2 )+ m 2 l 1 l 2 ( d 2 θ 2 d t 2 )cos( θ 1 θ 2 )+ m 2 l 1 l 2 ( d θ 2 dt ) 2 sin( θ 1 θ 2 )+( m 1 + m 2 )g l 1 sin θ 1 =0

La segunda ecuación del movimiento, se obtiene

d dt L θ ˙ 2 L θ 2 =0 θ ˙ 2 = d θ 2 dt L θ ˙ 2 = m 2 l 2 2 ( d θ 2 dt )+ m 2 l 1 l 2 ( d θ 1 dt )cos( θ 1 θ 2 ) L θ 1 = m 2 l 1 l 2 ( d θ 1 dt )( d θ 2 dt )sin( θ 1 θ 2 ) m 2 g l 2 sin θ 2 d dt ( L θ ˙ 1 )= m 2 l 2 2 ( d 2 θ 2 d t 2 )+ m 2 l 1 l 2 ( d 2 θ 1 d t 2 )cos( θ 1 θ 2 ) m 2 l 1 l 2 ( d θ 1 dt )sin( θ 1 θ 2 )( d θ 1 dt + d θ 2 dt ) m 2 l 2 2 ( d 2 θ 2 d t 2 )+ m 2 l 1 l 2 ( d 2 θ 1 d t 2 )cos( θ 1 θ 2 ) m 2 l 1 l 2 ( d θ 1 dt ) 2 sin( θ 1 θ 2 )+ m 2 g l 2 sin θ 2 =0

Dividimos ambas ecuaciones por m2l1l2 y llamamos m=m1/m2

(m+1) l 1 l 2 ( d 2 θ 1 d t 2 )+( d 2 θ 2 d t 2 )cos( θ 1 θ 2 )+ ( d θ 2 dt ) 2 sin( θ 1 θ 2 )+(m+1) g l 2 sin θ 1 =0 l 2 l 1 ( d 2 θ 2 d t 2 )+( d 2 θ 1 d t 2 )cos( θ 1 θ 2 ) ( d θ 1 dt ) 2 sin( θ 1 θ 2 )+ g l 1 sin θ 2 =0

Tenemos un sistema de dos ecuaciones con dos incógnitas en las que tenemos que despejar la derivada segunda de la posición angular de θ1 y θ2 de cada uno de los péndulos para proceder a su integración numérica.

{ a 1 d 2 θ 1 d t 2 +b d 2 θ 2 d t 2 + c 1 =0 a 2 d 2 θ 2 d t 2 +b d 2 θ 1 d t 2 + c 2 =0 d 2 θ 1 d t 2 = a 2 c 1 b c 2 b 2 a 1 a 2 d 2 θ 2 d t 2 = a 1 c 2 b c 1 b 2 a 1 a 2

El resultado es

d 2 θ 1 d t 2 = cos( θ 1 θ 2 )( g l 1 sin θ 2 ( d θ 1 dt ) 2 sin( θ 1 θ 2 ) ) l 2 l 1 ( (m+1) g l 2 sin θ 1 + ( d θ 2 dt ) 2 sin( θ 1 θ 2 ) ) m+ sin 2 ( θ 1 θ 2 ) d 2 θ 2 d t 2 = cos( θ 1 θ 2 )( (m+1) g l 2 sin θ 1 + ( d θ 2 dt ) 2 sin( θ 1 θ 2 ) )(m+1) l 1 l 2 ( g l 1 sin θ 2 ( d θ 1 dt ) 2 sin( θ 1 θ 2 ) ) m+ sin 2 ( θ 1 θ 2 )

Se resuelve el sistema de ecuaciones diferenciales por procedimientos numéricos con las siguientes condiciones iniciales: en el instante t=0, θ1=θ10 y 1/dt=0, θ2=θ20 y 2/dt=0

La energía total del sistema de dos partículas es constante, lo que nos permitirá verificar si el procedimiento numérico aplicado da buenos resultados

E= 1 2 ( m 1 + m 2 ) l 1 2 ( d θ 1 dt ) 2 + 1 2 m 2 l 2 2 ( d θ 2 dt ) 2 + m 2 l 1 l 2 ( d θ 1 dt )( d θ 2 dt )cos( θ 1 θ 2 ) ( m 1 + m 2 )g l 1 cos θ 1 m 2 g l 2 cos θ 2

m=1;  %relación m1/m2
L1=1;
L2=2;
x0=[pi/3, 0, -pi/6, 0];
tspan=[0,10];

fg=@(t,x)[x(2);(cos(x(1)-x(3))*(9.8*sin(x(3))/L1-x(2)^2*sin(x(1)-x(3)))-
(L2*((m+1)*9.8*sin(x(1))/L2+x(4)^2*sin(x(1)-x(3)))/L1))
/(m+sin(x(1)-x(3))*sin(x(1)-x(3))); x(4);
(cos(x(1)-x(3))*((m+1)*9.8*sin(x(1))/L2+x(4)^2*sin(x(1)-x(3)))
-((m+1)*L1*(9.8*sin(x(3))
/L1-x(2)^2*sin(x(1)-x(3)))/L2))/(m+sin(x(1)-x(3))*sin(x(1)-x(3)))];
[t,x]=ode45(fg,tspan,x0);
hold on
plot(x(:,1),x(:,2),'r') 
plot(x(:,3),x(:,4),'b')
hold off
grid on
legend('1','2');
xlabel('\theta_1, \theta_2')
ylabel('\omega_1, \omega_2')
title('Péndulo doble')

Einicial=-(m+1)*9.8*cos(x0(1))/L2-9.8*cos(x0(3))/L1;
E=(m+1)*L1*x(:,2).^2/(2*L2)+L2*x(:,4).^2/(2*L1)+(x(:,2).*x(:,4)).
*cos(x(:,1)-x(:,3))-(m+1)*9.8*cos(x(:,1))/L2-9.8*cos(x(:,3))/L1;
max(abs(100*(Einicial-E)/Einicial)) 
%error máximo en el intervalo de tiempo especificado
ans =    0.7792

Representación de la trayectoria de las dos partículas en el espacio de las fases: θ1, 1/dt en color rojo y θ2, 2/dt en color azul

Actividades

Se introduce

Se pulsa el botón titulado Nuevo

Se observa el movimiento de los dos partículas en el espacio de las fases: θ1, 1/dt en color rojo y θ2, 2/dt en color azul.

Se calcula el tanto por ciento de error,

| E E 0 E 0 |·100 .

Donde E0=-(m1+m2)·gl1cos θ10-m2gl2cos θ20, es la energía inicial y E es la energía en el instante t. Ambas cantidades deberían ser iguales. Cuando el tanto por ciento de error supera el 1%, el proceso de integración se detiene, se considera el procedimiento numérico no resuelve adecuadamente el sistema de ecuaciones diferenciales y se invita al lector a disminuir el Paso de integración, que por defecto, se ha establecido en h=0.01

Se puede controlar la trayectoria en el control titulado Trayectoria, n. puntos. Se dibuja un punto por cada paso de integración de la ecuación diferencial. Por defecto, se muestran los últimos 400 puntos. Se puede cambiar este número para ver con más detalle la trayectoria en el espacio de las fases



Aproximaciones

Si nos restringimos a pequeños valores de los ángulos θ1 y θ2, las ecuaciones del movimiento se hacen mucho más simples. Desarrollamos en serie de Taylor cosθ y tomamos los dos primeros términos no nulos del desarrollo en serie

cosθ1 1 2 θ 2

La energía potencial con esta aproximación, se expresa

E p ( m 1 + m 2 )g l 1 ( 1 1 2 θ 1 2 ) m 2 g l 2 ( 1 1 2 θ 2 2 )= 1 2 ( m 1 + m 2 )g l 1 θ 1 2 + 1 2 m 2 g l 2 θ 2 2 +cte

La energía cinética con esta aproximación, se expresa

E k 1 2 ( m 1 + m 2 ) l 1 2 ( d θ 1 dt ) 2 + 1 2 m 2 l 2 2 ( d θ 2 dt ) 2 + m 2 l 1 l 2 ( d θ 1 dt )( d θ 2 dt )

Las ecuaciones del movimiento de Lagrange nos llevan al sistema de dos ecuaciones diferenciales de segundo orden. La lagrangiana L=Ek-Ep con el símbolo θ ˙ = dθ dt

L= 1 2 ( m 1 + m 2 ) l 1 2 θ ˙ 1 2 + 1 2 m 2 l 2 2 θ ˙ 2 2 + m 2 l 1 l 2 θ ˙ 1 θ ˙ 2 1 2 ( m 1 + m 2 )g l 1 θ 1 2 1 2 m 2 g l 2 θ 2 2 d dt L θ ˙ 1 L θ 1 =0 d dt L θ ˙ 2 L θ 2 =0

Las ecuaciones del movimiento son

( m 1 + m 2 ) l 1 d 2 θ 1 d t 2 +( m 1 + m 2 )g θ 1 + m 2 l 2 d 2 θ 2 d t 2 =0 l 2 d 2 θ 2 d t 2 +g θ 2 + l 1 d 2 θ 1 d t 2 =0

Calculamos las derivadas segundas respecto del tiempo de estas dos ecuaciones diferenciales

( m 1 + m 2 ) l 1 d 4 θ 1 d t 4 +( m 1 + m 2 )g d 2 θ 1 d t 2 + m 2 l 2 d 4 θ 2 d t 4 =0 l 2 d 4 θ 2 d t 4 +g d 2 θ 2 d t 2 + l 1 d 4 θ 1 d t 4 =0

Entre estas cuatro ecuaciones, eliminamos θ2, obteniendo la ecuación diferencial de cuarto orden en θ1.

m 1 l 1 d 4 θ 1 d t 4 +( m 1 + m 2 )g( 1+ l 1 l 2 ) d 2 θ 1 d t 2 + ( m 1 + m 2 ) g 2 l 2 θ 1 =0

Suponiendo una solución de la forma

θ1=Asin(ωt)+Bcos(ωt) o bien, θ1=Asin(ωt+φ)

Calculamos la derivada segunda y la derivada cuarta de θ1 y las introducimos en la ecuación diferencial de cuarto orden, obteniendo la siguiente ecuación bicuadrada en ω.

m 1 l 1 ω 4 ( m 1 + m 2 )g( 1+ l 1 l 2 ) ω 2 + ( m 1 + m 2 ) g 2 l 2 =0

Las dos raíces reales de esta ecuación son las frecuencias ω1 y ω2 de los modos normales de vibración

ω 1 2 = ( m 1 + m 2 )g 2 m 1 l 1 { ( 1+ l 1 l 2 ) ( 1+ l 1 l 2 ) 2 4 m 1 ( m 1 + m 2 ) l 1 l 2 } ω 2 2 = ( m 1 + m 2 )g 2 m 1 l 1 { ( 1+ l 1 l 2 )+ ( 1+ l 1 l 2 ) 2 4 m 1 ( m 1 + m 2 ) l 1 l 2 }

La forma general del ángulo θ1 en función del tiempo t es una combinación de los dos modos normales de vibración

θ1(t)=A1sin(ω1t)+ B1cos(ω1t)+ C1sin(ω2t)+ D1cos(ω2t)

Lo mismo para θ2

θ2(t)=A2sin(ω1t)+ B2cos(ω1t)+ C2sin(ω2t)+ D2cos(ω2t)

Las velocidades angulares son

d θ 1 dt = A 1 ω 1 cos( ω 1 t) B 1 ω 1 sin( ω 1 t)+ C 1 ω 2 cos( ω 2 t) D 1 ω 2 sin( ω 2 t) d θ 2 dt = A 2 ω 1 cos( ω 1 t) B 2 ω 1 sin( ω 1 t)+ C 2 ω 2 cos( ω 2 t) D 2 ω 2 sin( ω 2 t)

Relacionamos los coeficientes A1 y A2, B1 y B2, C1 y C2, D1 y D2, calculando la derivada segunda de θ2, la derivada segunda de θ1 e introduciéndolos en una de las dos ecuaciones diferenciales del movimiento

l 2 d 2 θ 2 d t 2 +g θ 2 + l 1 d 2 θ 1 d t 2 =0

Tenemos cuatro relaciones

( ω 1 2 l 2 +g) A 2 ω 1 2 l 1 A 1 =0 ( ω 1 2 l 2 +g) B 2 ω 1 2 l 1 B 1 =0 ( ω 2 2 l 2 +g) C 2 ω 2 2 l 1 C 1 =0 ( ω 2 2 l 2 +g) D 2 ω 2 2 l 1 D 1 =0

Las condiciones iniciales determinan los valores de los cuatro coeficientes A1, B1, C1 y D1 y por tanto, de los coeficientes A2, B2, C2 y D2, a través de las relaciones establecidas anteriormente.

Desplazamos el primer péndulo un ángulo θ10 con respecto a la vertical y el segundo péndulo un ángulo θ20 a continuación, los soltamos. La velocidad inicial es  1/dt=0, 2/dt=0 en el instante t=0.

Resolvemos el sistema de cuatro ecuaciones junto a las relaciones establecidas anteriormente.

θ10=B1+D1
0=ω1·A1+ ω2·C1
θ20
=B2+D2
0=ω1·A2+ ω2·C2

La solución de este sistema es

A1=C1=A2=C2=0

B 2 = ω 1 2 ( ω 2 2 l 1 θ 10 + ω 2 2 l 2 θ 20 g θ 20 ) g( ω 2 2 ω 1 2 ) B 1 =( ω 1 2 l 2 +g ω 1 2 l 1 ) B 2 D 2 = ω 2 2 ( ω 1 2 l 1 θ 10 ω 1 2 l 2 θ 20 +g θ 20 ) g( ω 2 2 ω 1 2 ) D 1 =( ω 2 2 l 2 +g ω 2 2 l 1 ) D 2

La posición angular en función del tiempo t de cada una de las partículas es

θ1(t)=B1cos(ω1t)+ D1cos(ω2t)
θ2
(t)=B2cos(ω1t)+D2cos(ω2t)

Modos normales de vibración

Cuando los péndulos son iguales

Cuando m1=m2 y l1=l2 las frecuencias ω1 y ω2 de los modos normales de vibración tienen una expresión mucho más simple

ω 1 2 = g l (2 2 ) ω 2 2 = g l (2+ 2 )

Los coeficientes B2, B1, D2 y D1 valen

B 2 = 2 1 2 ( (2+ 2 )( θ 10 + θ 20 ) θ 20 ) B 1 = 2 2 B 2 D 2 = 2 +1 2 ( (2 2 )( θ 10 + θ 20 )+ θ 20 ) D 1 = 2 2 D 2

Es importante analizar el caso de que θ20=0. Se desplaza el primer péndulo un ángulo θ10 respecto de la vertical y se suelta

B 2 = 2 2 θ 10 B 1 = 1 2 θ 10   D 2 = 2 2 θ 10 D 1 = 1 2 θ 10

Las ecuaciones del movimiento de cada unos de las partículas se escriben

θ 1 (t)= 1 2 θ 10 ( cos( ω 1 t)+cos( ω 2 t) )= θ 10 cos( ω 2 ω 1 2 t )cos( ω 2 + ω 1 2 t ) θ 2 (t)= 2 2 θ 10 ( cos( ω 1 t)cos( ω 2 t) )= 2 θ 10 sin( ω 2 ω 1 2 t )sin( ω 2 + ω 1 2 t )

En la figura, se representa en color rojo θ1(t) y en color azul la amplitud modulada θ 10 cos( ω 2 ω 1 2 t )

m1=1; %masa del primer péndulo
m2=1;  %masa del segundo péndulo
lon1=1;  %longitud del primer péndulo
lon2=1; %longitud del segundo péndulo
ang1=pi/18; %angulo inicial del primer péndulo
ang2=0; %angulo inicial del segundo péndulo

w1=((m1+m2)*9.8/(2*m1*lon1))*((1+lon1/lon2)
-sqrt((1+lon1/lon2)*(1+lon1/lon2)-4*m1*lon1/((m1+m2)*lon2)));
w2=((m1+m2)*9.8/(2*m1*lon1))*((1+lon1/lon2)
+sqrt((1+lon1/lon2)*(1+lon1/lon2)-4*m1*lon1/((m1+m2)*lon2)));
%coeficientes
B2=w1*(w2*lon1*ang1+w2*lon2*ang2-9.8*ang2)/(9.8*(w2-w1));
D2=w2*(-w1*lon1*ang1-w1*lon2*ang2+9.8*ang2)/(9.8*(w2-w1));
B1=(-w1*lon2+9.8)*B2/(w1*lon1);
D1=(-w2*lon2+9.8)*D2/(w2*lon1);

t=linspace(0,10,200);
y1=B1*cos(sqrt(w1)*t)+D1*cos(sqrt(w2)*t);
y2=B2*cos(sqrt(w1)*t)+D2*cos(sqrt(w2)*t);
  
figure
amplitud1=ang1*cos((sqrt(w2)-sqrt(w1))*t/2);
plot(t,y1,'r',t,amplitud1,'b')
grid on
xlabel('t')
ylabel('\theta_1')
title('Posición angular del primer péndulo')

figure
amplitud2=sqrt(2)*ang1*sin((sqrt(w2)-sqrt(w1))*t/2);
plot(t,y2,'r',t,amplitud2,'b')
grid on
xlabel('t')
ylabel('\theta_2')
title('Posición angular del segundo péndulo')

En la figura, se representa en color rojo θ2(t) y en color azul la amplitud modulada 2 θ 10 sin( ω 2 ω 1 2 t )

Actividades

Se introduce

Se pulsa el botón titulado Nuevo

Se observa el movimiento de los dos péndulos y la representación gráfica del ángulo de desviación de cada péndulo θ1 y θ2 en función del tiempo t.

En la parte superior derecha, se muestra el valor de la energía total que permanece constante durante el movimiento de los péndulos.

Para observar los modos normales de vibración, se introducen los siguientes datos. Para el primer modo: Angulo 1, 7.07, Angulo 2, 10. Para el segundo modo: Angulo 1, -7.07, Angulo 2, 10



Referencias

Lee S. M. The double-simple pendulum problem. Am. J. Phys. 38 (1970) pp. 536-537