Un péndulo

El péndulo compuesto está formado por una varilla de masa m y longitud l, que está suspendido de su extremo. El momento de inercia de una varilla respecto de un eje perpendicular a la misma y que pasa por su centro de masa es I_c=ml²/12

Las coordenadas del centro de masa (punto de color rojo) son x₁=(l/2)sinθ₁, y₁=(l/2)cosθ₁. Derivando con respecto del tiempo obtenemos la velocidad del centro de masas.

${\displaystyle \overrightarrow{v_1}=\frac{d{x}_1}{dt}\widehat{i}+\frac{d{y}_1}{dt}\widehat{j}=\frac{l}{2}\left(\cos {\theta }_1\widehat{i}-\sin {\theta }_1\widehat{j}\right)\frac{d{\theta }_1}{dt}}$

La energía cinética, T₁ es la suma de la energía cinética de traslación del c.m. con velocidad v₁ y la energía cinética de rotación alrededor de un eje perpendicular a la varilla y que pasa por el c.m. con velocidad angular dθ₁/dt

${\displaystyle T_1=\frac{1}{2}m{v}_1^2+\frac{1}{2}\left(\frac{1}{12}m{l}^2\right){\left(\frac{d{\theta }_1}{dt}\right)}^2=\frac{1}{2}\left\{\frac{1}{3}m{l}^2{\left(\frac{d{\theta }_1}{dt}\right)}^2\right\}}$

que coincide con la energía cinética de rotación alrededor de un eje perpendicular a la varilla y que pasa por su extremo. El momento de inercia de la varilla respecto de dicho eje, aplicando Steiner es, I_o=ml²/12+m(l/2)²=ml²/3

La variación de energía potencial del c.m. es

${\displaystyle V_1=mg\frac{l}{2}-mg{y}_1=mg\frac{l}{2}\left(1-\cos {\theta }_1\right)\approx mg\frac{l}{4}{\theta }_1^2}$

Para ángulos θ pequeños, hacemos la aproximación 1-cosθ≈θ²/2

La lagrangiana L₁=T₁-V₁ es

${\displaystyle L_1=\frac{1}{2}\left\{\frac{1}{3}m{l}^2{\left(\frac{d{\theta }_1}{dt}\right)}^2\right\}-\frac{1}{4}mgl{\theta }_1^2}$

La ecuación del movimiento

${\displaystyle \begin{array}{l}\frac{d}{dt}\left(\frac{\partial L_1}{\partial {\dot{\theta }}_1}\right)-\frac{\partial L_1}{\partial {\theta }_1}=0\\ \frac{1}{3}m{l}^2\frac{d^2{\theta }_1}{d{t}^2}+\frac{1}{2}mlg{\theta }_1=0\\ \frac{d^2{\theta }_1}{d{t}^2}+\frac{3g}{2l}{\theta }_1=0\end{array}}$

La solución de esta ecuación diferencial es θ₁=Asin(ωt)+Bcos(ωt), con ω²=3g/(2l). Los coeficientes A y B se determinan a partir de las condiciones iniciales (posición y velocidad)

Dos péndulos acoplados

Las coordenadas del c.m. de la segunda varilla son x₂=lsinθ₁+(l/2)sinθ₂, y₂=lcosθ₁+(l/2)cosθ₂. Derivando con respecto del tiempo obtenemos la velocidad del c.m.

${\displaystyle \begin{array}{l}\overrightarrow{v_2}=\frac{d{x}_2}{dt}\widehat{i}+\frac{d{y}_2}{dt}\widehat{j}=\\ \left(l\cos {\theta }_1\frac{d{\theta }_1}{dt}+\frac{l}{2}\cos {\theta }_2\frac{d{\theta }_2}{dt}\right)\widehat{i}-\left(l\sin {\theta }_1\frac{d{\theta }_1}{dt}+\frac{l}{2}\sin {\theta }_2\frac{d{\theta }_2}{dt}\right)\widehat{j}\end{array}}$

La energía cinética, T₂ es la suma de la energía cinética de traslación del c.m. con velocidad v₂ y la energía cinética de rotación alrededor de un eje perpendicular a la varilla y que pasa por el c.m. con velocidad angular dθ₂/dt

${\displaystyle \begin{array}{l}{T}_2=\frac{1}{2}m{v}_2^2+\frac{1}{2}\left(\frac{1}{12}m{l}^2\right){\left(\frac{d{\theta }_2}{dt}\right)}^2=\\ \frac{1}{2}m{l}^2\left\{{\left(\frac{d{\theta }_1}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_2}{dt}\right)}^2+\cos \left({\theta }_1-{\theta }_2\right)\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}\right\}\approx \frac{1}{2}m{l}^2\left\{{\left(\frac{d{\theta }_1}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_2}{dt}\right)}^2+\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}\right\}\end{array}}$

Para ángulos θ₁ y θ₂ pequeños, hacemos la aproximación cos(θ₁-θ₂)≈1

La variación de energía potencial del c.m. es

${\displaystyle V_2=mg\left(l+\frac{l}{2}\right)-mg{y}_2=mg\left\{l\left(1-\cos {\theta }_1\right)+\frac{l}{2}\left(1-\cos {\theta }_2\right)\right\}\approx mg\frac{l}{2}\left({\theta }_1^2+\frac{1}{2}{\theta }_2^2\right)}$

Para ángulos θ pequeños, hacemos la aproximación 1-cosθ≈θ²/2

La lagrangiana L₂=L₁+T₂-V₂ es

${\displaystyle \begin{array}{l}{L}_2=L_1+T_2-V_2=\\ \frac{1}{2}m{l}^2\left\{\frac{1}{3}{\left(\frac{d{\theta }_1}{dt}\right)}^2\right\}-\frac{1}{4}mgl{\theta }_1^2+\frac{1}{2}m{l}^2\left\{{\left(\frac{d{\theta }_1}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_2}{dt}\right)}^2+\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}\right\}-mg\frac{l}{2}\left({\theta }_1^2+\frac{1}{2}{\theta }_2^2\right)=\\ \frac{1}{2}m{l}^2\left\{\frac{4}{3}{\left(\frac{d{\theta }_1}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_2}{dt}\right)}^2+\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}\right\}-mg\frac{l}{2}\left(\frac{3}{2}{\theta }_1^2+\frac{1}{2}{\theta }_2^2\right)\end{array}}$

Las ecuaciones del movimiento

${\displaystyle \{\begin{array}{l}\frac{4}{3}\frac{d^2{\theta }_1}{d{t}^2}+\frac{1}{2}\frac{d^2{\theta }_2}{d{t}^2}+\frac{3}{2}\frac{g}{l}{\theta }_1=0\\ \frac{1}{2}\frac{d^2{\theta }_1}{d{t}^2}+\frac{1}{3}\frac{d^2{\theta }_2}{d{t}^2}+\frac{1}{2}\frac{g}{l}{\theta }_2=0\end{array}}$

La forma general del ángulo θ₁ en función del tiempo t es una combinación de los dos modos normales de vibración de frecuencias ω₁ y ω₂ que calcularemos más adelante.

θ₁(t)=A₁sin(ω₁t)+ B₁cos(ω₁t)+ C₁sin(ω₂t)+ D₁cos(ω₂t)

Lo mismo para θ₂

θ₂(t)=A₂sin(ω₁t)+ B₂cos(ω₁t)+ C₂sin(ω₂t)+ D₂cos(ω₂t)

Las dos ecuaciones diferenciales proporcionan cuatro relaciones entre los coeficientes y las condiciones iniciales (posición y velocidad) otras cuatro de las que se despejan, A₁, A₂, B₁, B₂, C₁, C₂. Como hemos visto en el péndulo doble, siempre es posible encontrar las condiciones iniciales que hacen que algunos coeficientes sean cero. De modo que, la solución más simple sea θ₁=A₁sin(ω₁t), θ₂=A₂sin(ω₁t), para ω₁, y θ₁=C₁sin(ω₂t), θ₂=C₂sin(ω₂t), para ω₂.

Vamos a calcular y representar los modos normales de vibración. Suponiendo una solución de la forma θ ₁=A₁sin(ωt), θ₂=A₂sin(ωt), obtenemos el sistema de ecuaciones

${\displaystyle \{\begin{array}{l}\left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)A_1-\frac{1}{2}{\omega }^2{A}_2=0\\ -\frac{1}{2}{\omega }^2{A}_1+\left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)A_2=0\end{array}}$

En este sistema homogéneo de dos ecuaciones con dos incógnitas, el determinante de los coeficientes ha de ser cero

${\displaystyle \begin{array}{l}\left|\begin{array}{cc}\left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& -\frac{1}{2}{\omega }^2\\ -\frac{1}{2}{\omega }^2& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|=0\\ \frac{7}{36}{\omega }^4-\frac{7}{6}\frac{g}{l}{\omega }^2+\frac{3}{4}\frac{g^2}{l^2}=0\{\begin{array}{l}{\omega }_1^2=\left(3-\frac{6\sqrt{7}}{7}\right)\frac{g}{l}\\ {\omega }_2^2=\left(3+\frac{6\sqrt{7}}{7}\right)\frac{g}{l}\end{array}\end{array}}$

En el código x es ω² y g es g/l

>> syms g x;
>> A=[(3*g/2-4*x/3),-x/2;-x/2,(g/2-x/3)];
>> eq=det(A)
eq = (3*g^2)/4 - (7*g*x)/6 + (7*x^2)/36
>> solve(eq,x)
ans = 
 3*g - (6*7^(1/2)*g)/7
 3*g + (6*7^(1/2)*g)/7

Tomamos la segunda ecuación del sistema homogéneo de dos ecuaciones, para calcular la amplitud A₂ en términos de la amplitud A₁

${\displaystyle \begin{array}{l}-\frac{1}{2}{\omega }^2{A}_1+\left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)A_2=0\\ A_2=\frac{{\omega }^2}{\frac{g}{l}-\frac{2}{3}{\omega }^2}{A}_1\end{array}}$

Calculamos A₂ para el primer modo de vibración de frecuencia angular ω=ω₁. El ángulo θ₁ que forma el primer péndulo con la vertical será proporcional a A₁ y el ángulo θ₂ que forma el segundo péndulo con la vertical será proporcional a A₂. Representamos el primer modo de vibración.

L=1; %longitud total
n=2; %número de péndulos
w2=[3-6*sqrt(7)/7,3+6*sqrt(7)/7]*9.8/(L/n); %cuadrado frecuencias
A1=0.2; %escala
iModo=1; %modo de vibración
A2=w2(iModo)*A1/(9.8/(L/n)-2*w2(iModo)/3);
angulos=[A1,A2];

hold on
x1=0; y1=0; %origen
plot(x1,y1,'ko','markersize',4,'markerfacecolor','k')
for j=1:n
    x2=x1+(L/n)*sin(angulos(j));
    y2=y1-(L/n)*cos(angulos(j));
    line([x1,x2],[y1,y2], 'lineWidth',1.5)
    if j==n
        break;
    end
    plot(x2,y2,'ro','markersize',3,'markerfacecolor','r')
    x1=x2;
    y1=y2;
end
text(0.1,-0.1, sprintf('frecuencia angular, %1.3f',sqrt(w2(iModo))))
hold off
axis equal
xlabel('x')
ylabel('y')
grid on
title(sprintf('Modo de vibración, %i',iModo))

Repetimos para el segundo modo de vibración de frecuencia angular ω=ω₂, cambiando el valor de la variable iModo=2

Tres péndulos

Las coordenadas del c.m. de la tercera varilla son x₃=lsinθ₁+lsinθ₂+(l/2)sinθ₃, y₃=lcosθ₁+lcosθ₂+(l/2)cosθ₃. Derivando con respecto del tiempo obtenemos la velocidad del c.m.

${\displaystyle \begin{array}{l}\overrightarrow{v_3}=\frac{d{x}_3}{dt}\widehat{i}+\frac{d{y}_3}{dt}\widehat{j}=\\ \left(l\cos {\theta }_1\frac{d{\theta }_1}{dt}+l\cos {\theta }_2\frac{d{\theta }_2}{dt}+\frac{l}{2}\cos {\theta }_3\frac{d{\theta }_3}{dt}\right)\widehat{i}-\left(l\sin {\theta }_1\frac{d{\theta }_1}{dt}+l\sin {\theta }_2\frac{d{\theta }_2}{dt}+\frac{l}{2}\sin {\theta }_3\frac{d{\theta }_3}{dt}\right)\widehat{j}\end{array}}$

La energía cinética, T₃ es la suma de la energía cinética de traslación del c.m. con velocidad v₃ y la energía cinética de rotación alrededor de un eje perpendicular a la varilla y que pasa por el c.m. con velocidad angular dθ₃/dt

${\displaystyle \begin{array}{l}{T}_3=\frac{1}{2}m{v}_3^2+\frac{1}{2}\left(\frac{1}{12}m{l}^2\right){\left(\frac{d{\theta }_3}{dt}\right)}^2=\\ \frac{1}{2}m\left\{l^2{\left(\frac{d{\theta }_1}{dt}\right)}^2+l^2{\left(\frac{d{\theta }_2}{dt}\right)}^2+\frac{l^2}{4}{\left(\frac{d{\theta }_3}{dt}\right)}^2+2l^2\cos \left({\theta }_1-{\theta }_2\right)\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}+l^2\cos \left({\theta }_1-{\theta }_3\right)\frac{d{\theta }_1}{dt}\frac{d{\theta }_3}{dt}+l^2\cos \left({\theta }_2-{\theta }_3\right)\frac{d{\theta }_2}{dt}\frac{d{\theta }_3}{dt}\right\}+\frac{1}{2}\left(\frac{1}{12}m{l}^2\right){\left(\frac{d{\theta }_3}{dt}\right)}^2\approx \\ \approx \frac{1}{2}m{l}^2\left\{{\left(\frac{d{\theta }_1}{dt}\right)}^2+{\left(\frac{d{\theta }_2}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_3}{dt}\right)}^2+2\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}+\frac{d{\theta }_1}{dt}\frac{d{\theta }_3}{dt}+\frac{d{\theta }_2}{dt}\frac{d{\theta }_3}{dt}\right\}\end{array}}$

Para ángulos θ₁, θ₂ y θ₃ pequeños, hacemos la aproximación cos(θ₁-θ₂)≈1, cos(θ₁-θ₃)≈1 y cos(θ₂-θ₃)≈1

La variación de energía potencial del c.m. es

${\displaystyle \begin{array}{l}{V}_3=mg\left(2l+\frac{l}{2}\right)-mg{y}_3=mg\left\{l\left(1-\cos {\theta }_1\right)+l\left(1-\cos {\theta }_2\right)+\frac{l}{2}\left(1-\cos {\theta }_3\right)\right\}\\ \approx mg\frac{l}{2}\left({\theta }_1^2+{\theta }_2^2+\frac{1}{2}{\theta }_3^2\right)\end{array}}$

La lagrangiana es

${\displaystyle \begin{array}{l}{L}_3=L_2+T_3-V_3=\frac{1}{2}m{l}^2\left\{\frac{4}{3}{\left(\frac{d{\theta }_1}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_2}{dt}\right)}^2+\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}\right\}-mg\frac{l}{2}\left(\frac{3}{2}{\theta }_1^2+\frac{1}{2}{\theta }_2^2\right)+\\ \frac{1}{2}m{l}^2\left\{{\left(\frac{d{\theta }_1}{dt}\right)}^2+{\left(\frac{d{\theta }_2}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_3}{dt}\right)}^2+2\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}+\frac{d{\theta }_1}{dt}\frac{d{\theta }_3}{dt}+\frac{d{\theta }_2}{dt}\frac{d{\theta }_3}{dt}\right\}-mg\frac{l}{2}\left({\theta }_1^2+{\theta }_2^2+\frac{1}{2}{\theta }_3^2\right)=\\ \frac{1}{2}m{l}^2\left\{\frac{7}{3}{\left(\frac{d{\theta }_1}{dt}\right)}^2+\frac{4}{3}{\left(\frac{d{\theta }_2}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_3}{dt}\right)}^2+3\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}+\frac{d{\theta }_1}{dt}\frac{d{\theta }_3}{dt}+\frac{d{\theta }_2}{dt}\frac{d{\theta }_3}{dt}\right\}-mg\frac{l}{2}\left(\frac{5}{2}{\theta }_1^2+\frac{3}{2}{\theta }_2^2+\frac{1}{2}{\theta }_3^2\right)\end{array}}$

Las ecuaciones del movimiento

${\displaystyle \{\begin{array}{l}\frac{7}{3}\frac{d^2{\theta }_1}{d{t}^2}+\frac{3}{2}\frac{d^2{\theta }_2}{d{t}^2}+\frac{1}{2}\frac{d^2{\theta }_3}{d{t}^2}+\frac{5}{2}\frac{g}{l}{\theta }_1=0\\ \frac{3}{2}\frac{d^2{\theta }_1}{d{t}^2}+\frac{4}{3}\frac{d^2{\theta }_2}{d{t}^2}+\frac{1}{2}\frac{d^2{\theta }_3}{d{t}^2}+\frac{3}{2}\frac{g}{l}{\theta }_2=0\\ \frac{1}{2}\frac{d^2{\theta }_1}{d{t}^2}+\frac{1}{2}\frac{d^2{\theta }_2}{d{t}^2}+\frac{1}{3}\frac{d^2{\theta }_3}{d{t}^2}+\frac{1}{2}\frac{g}{l}{\theta }_3=0\end{array}}$

Vamos a calcular y representar los modos normales de vibración. Suponiendo una solución θ₁=A₁sin(ωt), θ₂=A₂sin(ωt), θ₃=A₃sin(ωt), obtenemos el sistema de ecuaciones

${\displaystyle \{\begin{array}{l}\left(\frac{5}{2}\frac{g}{l}-\frac{7}{3}{\omega }^2\right)A_1-\frac{3}{2}{\omega }^2{A}_2-\frac{1}{2}{\omega }^2{A}_3=0\\ -\frac{3}{2}{\omega }^2{A}_1+\left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)A_2-\frac{1}{2}{\omega }^2{A}_3=0\\ -\frac{1}{2}{\omega }^2{A}_1-\frac{1}{2}{\omega }^2{A}_2+\left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)A_3=0\end{array}}$

En este sistema homogéneo de tres ecuaciones con tres incógnitas, el determinante de los coeficientes ha de ser cero

${\displaystyle \left|\begin{array}{ccc}\left(\frac{5}{2}\frac{g}{l}-\frac{7}{3}{\omega }^2\right)& -\frac{3}{2}{\omega }^2& -\frac{1}{2}{\omega }^2\\ -\frac{3}{2}{\omega }^2& \left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& -\frac{1}{2}{\omega }^2\\ -\frac{1}{2}{\omega }^2& -\frac{1}{2}{\omega }^2& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|=0}$

Obtenemos la ecuación cúbica en ω². En el código x es ω² y g es g/l

>> syms g x;
>> A=[(5*g/2-7*x/3),-3*x/2,-x/2;-3*x/2,(3*g/2-4*x/3),-x/2;-x/2,-x/2,(g/2-x/3)];
>> eq=det(A)
eq = (15*g^3)/8 - (14*g^2*x)/3 + (41*g*x^2)/24 - (13*x^3)/108

${\displaystyle -\frac{13}{108}{\omega }^6+\frac{41}{24}\frac{g}{l}{\omega }^4-\frac{14}{3}\frac{g^2}{l^2}{\omega }^2+\frac{15}{8}\frac{g^3}{l^3}=0}$

Los cuadrados de las frecuencias de los modos normales de vibración son las raíces de esta ecuación

>> n=3; %número de péndulos
>> L=1; %longitud total
>> roots([-13/108,41*9.8/(24*(L/n)), -14*9.8^2/(3*(L/n)^2),15*9.8^3/(8*(L/n)^3)])
ans =
  314.7970
   88.2000
   14.2569

Tomamos la segunda y tercera ecuación del sistema homogéneo de tres ecuaciones, para calcular las amplitudes A₂ y A₃ en términos de la amplitud A₁

${\displaystyle \begin{array}{l}\{\begin{array}{l}\left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)A_2-\frac{1}{2}{\omega }^2{A}_3=\frac{3}{2}{\omega }^2{A}_1\\ -\frac{1}{2}{\omega }^2{A}_2+\left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)A_3=\frac{1}{2}{\omega }^2{A}_1\end{array}\\ A_2=\frac{\left|\begin{array}{cc}\frac{3}{2}& -\frac{1}{2}{\omega }^2\\ \frac{1}{2}& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|}{\left|\begin{array}{cc}\left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& -\frac{1}{2}{\omega }^2\\ -\frac{1}{2}{\omega }^2& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|}{\omega }^2{A}_1=\frac{-\frac{1}{4}{\omega }^2+\frac{3}{4}\frac{g}{l}}{\frac{7}{36}{\omega }^4-\frac{7}{6}\frac{g}{l}{\omega }^2+\frac{3}{4}\frac{g^2}{l^2}}{\omega }^2{A}_1\\ A_3=\frac{\left|\begin{array}{cc}\left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& \frac{3}{2}\\ -\frac{1}{2}{\omega }^2& \frac{1}{2}\end{array}\right|}{\left|\begin{array}{cc}\left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& -\frac{1}{2}{\omega }^2\\ -\frac{1}{2}{\omega }^2& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|}{\omega }^2{A}_1=\frac{\frac{1}{12}{\omega }^2+\frac{3}{4}\frac{g}{l}}{\frac{7}{36}{\omega }^4-\frac{7}{6}\frac{g}{l}{\omega }^2+\frac{3}{4}\frac{g^2}{l^2}}{\omega }^2{A}_1\end{array}}$

>> syms g x;
>> A=[(3*g/2-4*x/3),-x/2;-x/2,(g/2-x/3)];
>> det(A)
ans = (3*g^2)/4 - (7*g*x)/6 + (7*x^2)/36 %denominador
>> A=[3/2,-x/2;1/2,(g/2-x/3)];
>> det(A)
ans = (3*g)/4 - x/4 %numerador A2
>> A=[(3*g/2-4*x/3),3/2;-x/2,1/2];
>> det(A)
ans = (3*g)/4 + x/12 %numerador A3

Calculamos A₂ y A₃ para el primer modo de vibración de frecuencia angular ω=ω₁. El ángulo θ₁ que forma el primer péndulo con la vertical será proporcional a A₁, el ángulo θ₂ que forma el segundo péndulo con la vertical será proporcional a A₂ y el ángulo θ₃ que forma el tercer péndulo con la vertical será proporcional a A₃. Representamos el primer modo de vibración.

L=1; %longitud total
n=3; %número de péndulos
w2=sort(roots([-13/108,41*9.8/(24*(L/n)), -14*9.8^2/(3*(L/n)^2),
15*9.8^3/(8*(L/n)^3)])); 
A1=0.2; %escala
iModo=3; %modo de vibración
den=7*w2(iModo)^2/36-7*w2(iModo)*9.8/(6*L/n)+3*9.8^2/(4*(L/n)^2);
A2=(-w2(iModo)/4+3*9.8/(4*L/n))*w2(iModo)*A1/den;
A3=(w2(iModo)/12+3*9.8/(4*L/n))*w2(iModo)*A1/den;
angulos=[A1,A2, A3];

hold on
x1=0; y1=0; %origen
plot(x1,y1,'ko','markersize',4,'markerfacecolor','k')
for j=1:n
    x2=x1+(L/n)*sin(angulos(j));
    y2=y1-(L/n)*cos(angulos(j));
    line([x1,x2],[y1,y2], 'lineWidth',1.5)
    if j==n
        break;
    end
    plot(x2,y2,'ro','markersize',3,'markerfacecolor','r')
    x1=x2;
    y1=y2;
end
text(0.1,-0.1, sprintf('frecuencia angular, %1.3f',sqrt(w2(iModo))))
hold off
axis equal
xlabel('x')
ylabel('y')
grid on
title(sprintf('Modo de vibración, %i',iModo))

Repetimos para el segundo modo de vibración de frecuencia angular ω=ω₂, cambiando el valor de la variable iModo=2

Repetimos para el tercer modo de vibración de frecuencia angular ω=ω₃, cambiando el valor de la variable iModo=3

Cuatro péndulos

Las coordenadas del c.m. de la tercera varilla son x₄=lsinθ₁+lsinθ₂+lsinθ₃+(l/2)sinθ₄, y₄=lcosθ₁+lcosθ₂+lcosθ₃+(l/2)cosθ₄. Derivando con respecto del tiempo obtenemos la velocidad del c.m.

${\displaystyle \begin{array}{l}\overrightarrow{v_4}=\frac{d{x}_4}{dt}\widehat{i}+\frac{d{y}_4}{dt}\widehat{j}=\\ \left(l\cos {\theta }_1\frac{d{\theta }_1}{dt}+l\cos {\theta }_2\frac{d{\theta }_2}{dt}+l\cos {\theta }_3\frac{d{\theta }_3}{dt}+\frac{l}{2}\cos {\theta }_4\frac{d{\theta }_4}{dt}\right)\widehat{i}-\left(l\sin {\theta }_1\frac{d{\theta }_1}{dt}+l\sin {\theta }_2\frac{d{\theta }_2}{dt}+l\sin {\theta }_3\frac{d{\theta }_3}{dt}+\frac{l}{2}\sin {\theta }_4\frac{d{\theta }_4}{dt}\right)\widehat{j}\end{array}}$

La energía cinética, T₄ es la suma de la energía cinética de traslación del c.m. con velocidad v₄ y la energía cinética de rotación alrededor de un eje perpendicular a la varilla y que pasa por el c.m. con velocidad angular dθ₄/dt

${\displaystyle \begin{array}{l}{T}_4=\frac{1}{2}m{l}^2\left\{\begin{array}{l}{\left(\frac{d{\theta }_1}{dt}\right)}^2+{\left(\frac{d{\theta }_2}{dt}\right)}^2+{\left(\frac{d{\theta }_3}{dt}\right)}^2+\frac{1}{4}{\left(\frac{d{\theta }_4}{dt}\right)}^2+2\cos \left({\theta }_1-{\theta }_2\right)\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}+2\cos \left({\theta }_1-{\theta }_3\right)\frac{d{\theta }_1}{dt}\frac{d{\theta }_3}{dt}+\\ \cos \left({\theta }_1-{\theta }_4\right)\frac{d{\theta }_1}{dt}\frac{d{\theta }_4}{dt}+2\cos \left({\theta }_2-{\theta }_3\right)\frac{d{\theta }_2}{dt}\frac{d{\theta }_3}{dt}+\cos \left({\theta }_2-{\theta }_4\right)\frac{d{\theta }_2}{dt}\frac{d{\theta }_4}{dt}+\cos \left({\theta }_3-{\theta }_4\right)\frac{d{\theta }_3}{dt}\frac{d{\theta }_4}{dt}\end{array}\right\}+\frac{1}{2}\left(\frac{1}{12}m{l}^2\right){\left(\frac{d{\theta }_4}{dt}\right)}^2\approx \\ \frac{1}{2}m{l}^2\left\{{\left(\frac{d{\theta }_1}{dt}\right)}^2+{\left(\frac{d{\theta }_2}{dt}\right)}^2+{\left(\frac{d{\theta }_3}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_4}{dt}\right)}^2+2\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}+2\frac{d{\theta }_1}{dt}\frac{d{\theta }_3}{dt}+\frac{d{\theta }_1}{dt}\frac{d{\theta }_4}{dt}+2\frac{d{\theta }_2}{dt}\frac{d{\theta }_3}{dt}+\frac{d{\theta }_2}{dt}\frac{d{\theta }_4}{dt}+\frac{d{\theta }_3}{dt}\frac{d{\theta }_4}{dt}\right\}\end{array}}$

Para ángulos θ₁, θ₂, θ₃ y θ₄ pequeños, hacemos la aproximación cos(θ₁-θ₂)≈1, cos(θ₁-θ₃)≈1, cos(θ₁-θ₄)≈1, cos(θ₂-θ₃)≈1, cos(θ₂-θ₄)≈1 y cos(θ₃-θ₄)≈1

La variación de energía potencial del c.m. es

${\displaystyle \begin{array}{l}{V}_4=mg\left(3l+\frac{l}{2}\right)-mg{y}_4=mg\left\{l\left(1-\cos {\theta }_1\right)+l\left(1-\cos {\theta }_2\right)+l\left(1-\cos {\theta }_3\right)+\frac{l}{2}\left(1-\cos {\theta }_4\right)\right\}\\ \approx mg\frac{l}{2}\left({\theta }_1^2+{\theta }_2^2+{\theta }_3^2+\frac{1}{2}{\theta }_4^2\right)\end{array}}$

La lagrangiana es

${\displaystyle \begin{array}{l}{L}_4=L_3+T_4-V_4=\frac{1}{2}m{l}^2\left\{\frac{7}{3}{\left(\frac{d{\theta }_1}{dt}\right)}^2+\frac{4}{3}{\left(\frac{d{\theta }_2}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_3}{dt}\right)}^2+3\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}+\frac{d{\theta }_1}{dt}\frac{d{\theta }_3}{dt}+\frac{d{\theta }_2}{dt}\frac{d{\theta }_3}{dt}\right\}-mg\frac{l}{2}\left(\frac{5}{2}{\theta }_1^2+\frac{3}{2}{\theta }_2^2+\frac{1}{2}{\theta }_3^2\right)\\ \frac{1}{2}m{l}^2\left\{{\left(\frac{d{\theta }_1}{dt}\right)}^2+{\left(\frac{d{\theta }_2}{dt}\right)}^2+{\left(\frac{d{\theta }_3}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_4}{dt}\right)}^2+2\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}+2\frac{d{\theta }_1}{dt}\frac{d{\theta }_3}{dt}+\frac{d{\theta }_1}{dt}\frac{d{\theta }_4}{dt}+2\frac{d{\theta }_2}{dt}\frac{d{\theta }_3}{dt}+\frac{d{\theta }_2}{dt}\frac{d{\theta }_4}{dt}+\frac{d{\theta }_3}{dt}\frac{d{\theta }_4}{dt}\right\}-mg\frac{l}{2}\left({\theta }_1^2+{\theta }_2^2+{\theta }_3^2+\frac{1}{2}{\theta }_4^2\right)\\ \frac{1}{2}m{l}^2\left\{\frac{10}{3}{\left(\frac{d{\theta }_1}{dt}\right)}^2+\frac{7}{3}{\left(\frac{d{\theta }_2}{dt}\right)}^2+\frac{4}{3}{\left(\frac{d{\theta }_3}{dt}\right)}^2+\frac{1}{3}{\left(\frac{d{\theta }_4}{dt}\right)}^2+5\frac{d{\theta }_1}{dt}\frac{d{\theta }_2}{dt}+3\frac{d{\theta }_1}{dt}\frac{d{\theta }_3}{dt}+3\frac{d{\theta }_2}{dt}\frac{d{\theta }_3}{dt}+\frac{d{\theta }_1}{dt}\frac{d{\theta }_4}{dt}+\frac{d{\theta }_2}{dt}\frac{d{\theta }_4}{dt}+\frac{d{\theta }_3}{dt}\frac{d{\theta }_4}{dt}\right\}-mg\frac{l}{2}\left(\frac{7}{2}{\theta }_1^2+\frac{5}{2}{\theta }_2^2+\frac{3}{2}{\theta }_3^2+\frac{1}{2}{\theta }_4^2\right)\end{array}}$

Las ecuaciones del movimiento

${\displaystyle \{\begin{array}{l}\frac{10}{3}\frac{d^2{\theta }_1}{d{t}^2}+\frac{5}{2}\frac{d^2{\theta }_2}{d{t}^2}+\frac{3}{2}\frac{d^2{\theta }_3}{d{t}^2}+\frac{1}{2}\frac{d^2{\theta }_4}{d{t}^2}+\frac{7}{2}\frac{g}{l}{\theta }_1=0\\ \frac{5}{2}\frac{d^2{\theta }_1}{d{t}^2}+\frac{7}{3}\frac{d^2{\theta }_2}{d{t}^2}+\frac{3}{2}\frac{d^2{\theta }_3}{d{t}^2}+\frac{1}{2}\frac{d^2{\theta }_4}{d{t}^2}+\frac{5}{2}\frac{g}{l}{\theta }_2=0\\ \frac{3}{2}\frac{d^2{\theta }_1}{d{t}^2}+\frac{3}{2}\frac{d^2{\theta }_2}{d{t}^2}+\frac{4}{3}\frac{d^2{\theta }_3}{d{t}^2}+\frac{1}{2}\frac{d^2{\theta }_4}{d{t}^2}+\frac{3}{2}\frac{g}{l}{\theta }_3=0\\ \frac{1}{2}\frac{d^2{\theta }_1}{d{t}^2}+\frac{1}{2}\frac{d^2{\theta }_2}{d{t}^2}+\frac{1}{2}\frac{d^2{\theta }_3}{d{t}^2}+\frac{1}{3}\frac{d^2{\theta }_4}{d{t}^2}+\frac{1}{2}\frac{g}{l}{\theta }_4=0\end{array}}$

Vamos a calcular y representar los modos normales de vibración. Suponiendo una solución θ₁=A₁sin(ωt), θ₂=A₂sin(ωt), θ₃=A₃sin(ωt), θ₄=A₄sin(ωt) obtenemos el sistema de ecuaciones

${\displaystyle \{\begin{array}{l}\left(\frac{7}{2}\frac{g}{l}-\frac{10}{3}{\omega }^2\right)A_1-\frac{5}{2}{\omega }^2{A}_2-\frac{3}{2}{\omega }^2{A}_3-\frac{1}{2}{\omega }^2{A}_4=0\\ -\frac{5}{2}{\omega }^2{A}_1+\left(\frac{5}{2}\frac{g}{l}-\frac{7}{3}{\omega }^2\right)A_2-\frac{3}{2}{\omega }^2{A}_3-\frac{1}{2}{\omega }^2{A}_4=0\\ -\frac{3}{2}{\omega }^2{A}_1-\frac{3}{2}{\omega }^2{A}_2+\left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)A_3-\frac{1}{2}{\omega }^2{A}_4=0\\ -\frac{1}{2}{\omega }^2{A}_1-\frac{1}{2}{\omega }^2{A}_2-\frac{1}{2}{\omega }^2{A}_3+\left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)A_4=0\end{array}}$

En este sistema homogéneo de cuatro ecuaciones con cuatro incógnitas, el determinante de los coeficientes ha de ser cero

${\displaystyle \left|\begin{array}{cccc}\left(\frac{7}{2}\frac{g}{l}-\frac{10}{3}{\omega }^2\right)& -\frac{5}{2}{\omega }^2& -\frac{3}{2}{\omega }^2& -\frac{1}{2}{\omega }^2\\ -\frac{5}{2}{\omega }^2& \left(\frac{5}{2}\frac{g}{l}-\frac{7}{3}{\omega }^2\right)& -\frac{3}{2}{\omega }^2& -\frac{1}{2}{\omega }^2\\ -\frac{3}{2}{\omega }^2& -\frac{3}{2}{\omega }^2& \left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& -\frac{1}{2}{\omega }^2\\ -\frac{1}{2}{\omega }^2& -\frac{1}{2}{\omega }^2& -\frac{1}{2}{\omega }^2& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|=0}$

Obtenemos la ecuación cuarto orden en ω². En el código x es ω² y g es g/l

>> syms g x;
>> A=[(7*g/2-10*x/3),-5*x/2, -3*x/2,-x/2;-5*x/2,(5*g/2-7*x/3),-3*x/2,-x/2;
-3*x/2,-3*x/2,(3*g/2-4*x/3),-x/2; -x/2,-x/2,-x/2,(g/2-x/3)];
>> det(A)
ans =(105*g^4)/16 - (271*g^3*x)/12 + (943*g^2*x^2)/72 - (215*g*x^3)/108 + 
(97*x^4)/1296

${\displaystyle \frac{97}{1296}{\omega }^8-\frac{215}{108}\frac{g}{l}{\omega }^6+\frac{943}{72}\frac{g^2}{l^2}{\omega }^4-\frac{271}{12}\frac{g^3}{l^3}{\omega }^2+\frac{105}{16}\frac{g^4}{l^4}=0}$

Los cuadrados de las frecuencias de los modos normales de vibración son las raíces de esta ecuación

>> n=4; %número de péndulos
>> L=1; %longitud total
>> roots([97/1296,-215*9.8/(108*L/n),943*9.8^2/(72*(L/n)^2),
-271*9.8^3/(12*(L/n)^3),105*9.8^4/(16*(L/n)^4)])
ans =
  690.8821
  254.8446
   82.6924
   14.2201

Tomamos la segunda, tercera y cuarta ecuación del sistema homogéneo de cuatro ecuaciones, para calcular las amplitudes A₂, A₃ y A₄ en términos de la amplitud A₁

${\displaystyle \begin{array}{l}\{\begin{array}{l}\left(\frac{5}{2}\frac{g}{l}-\frac{7}{3}{\omega }^2\right)A_2-\frac{3}{2}{\omega }^2{A}_3-\frac{1}{2}{\omega }^2{A}_4=\frac{5}{2}{\omega }^2{A}_1\\ -\frac{3}{2}{\omega }^2{A}_2+\left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)A_3-\frac{1}{2}{\omega }^2{A}_4=\frac{3}{2}{\omega }^2{A}_1\\ -\frac{1}{2}{\omega }^2{A}_2-\frac{1}{2}{\omega }^2{A}_3+\left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)A_4=\frac{1}{2}{\omega }^2{A}_1\end{array}\\ A_2=\frac{\left|\begin{array}{ccc}\frac{5}{2}& -\frac{3}{2}{\omega }^2& -\frac{1}{2}{\omega }^2\\ \frac{3}{2}& \left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& -\frac{1}{2}{\omega }^2\\ \frac{1}{2}& -\frac{1}{2}{\omega }^2& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|}{\left|\begin{array}{ccc}\left(\frac{5}{2}\frac{g}{l}-\frac{7}{3}{\omega }^2\right)& -\frac{3}{2}{\omega }^2& -\frac{1}{2}{\omega }^2\\ -\frac{3}{2}{\omega }^2& \left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& -\frac{1}{2}{\omega }^2\\ -\frac{1}{2}{\omega }^2& -\frac{1}{2}{\omega }^2& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|}{\omega }^2{A}_1=\frac{\frac{11}{72}{\omega }^4-\frac{17}{12}\frac{g}{l}{\omega }^2+\frac{15}{8}\frac{g^2}{l^2}}{-\frac{13}{108}{\omega }^6+\frac{41}{24}\frac{g}{l}{\omega }^4-\frac{14}{3}\frac{g^2}{l^2}{\omega }^2+\frac{15}{8}\frac{g^3}{l^3}}{\omega }^2{A}_1\\ A_3=\frac{\left|\begin{array}{ccc}\left(\frac{5}{2}\frac{g}{l}-\frac{7}{3}{\omega }^2\right)& \frac{5}{2}& -\frac{1}{2}{\omega }^2\\ -\frac{3}{2}{\omega }^2& \frac{3}{2}& -\frac{1}{2}{\omega }^2\\ -\frac{1}{2}{\omega }^2& \frac{1}{2}& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|}{\left|\begin{array}{ccc}\left(\frac{5}{2}\frac{g}{l}-\frac{7}{3}{\omega }^2\right)& -\frac{3}{2}{\omega }^2& -\frac{1}{2}{\omega }^2\\ -\frac{3}{2}{\omega }^2& \left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& -\frac{1}{2}{\omega }^2\\ -\frac{1}{2}{\omega }^2& -\frac{1}{2}{\omega }^2& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|}{\omega }^2{A}_1=\frac{-\frac{1}{24}{\omega }^4-\frac{1}{2}\frac{g}{l}{\omega }^2+\frac{15}{8}\frac{g^2}{l^2}}{-\frac{13}{108}{\omega }^6+\frac{41}{24}\frac{g}{l}{\omega }^4-\frac{14}{3}\frac{g^2}{l^2}{\omega }^2+\frac{15}{8}\frac{g^3}{l^3}}{\omega }^2{A}_1\\ A_4=\frac{\left|\begin{array}{ccc}\left(\frac{5}{2}\frac{g}{l}-\frac{7}{3}{\omega }^2\right)& -\frac{3}{2}{\omega }^2& \frac{5}{2}\\ -\frac{3}{2}{\omega }^2& \left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& \frac{3}{2}\\ -\frac{1}{2}{\omega }^2& -\frac{1}{2}{\omega }^2& \frac{1}{2}\end{array}\right|}{\left|\begin{array}{ccc}\left(\frac{5}{2}\frac{g}{l}-\frac{7}{3}{\omega }^2\right)& -\frac{3}{2}{\omega }^2& -\frac{1}{2}{\omega }^2\\ -\frac{3}{2}{\omega }^2& \left(\frac{3}{2}\frac{g}{l}-\frac{4}{3}{\omega }^2\right)& -\frac{1}{2}{\omega }^2\\ -\frac{1}{2}{\omega }^2& -\frac{1}{2}{\omega }^2& \left(\frac{1}{2}\frac{g}{l}-\frac{1}{3}{\omega }^2\right)\end{array}\right|}{\omega }^2{A}_1=\frac{\frac{1}{72}{\omega }^4+\frac{1}{3}\frac{g}{l}{\omega }^2+\frac{15}{8}\frac{g^2}{l^2}}{-\frac{13}{108}{\omega }^6+\frac{41}{24}\frac{g}{l}{\omega }^4-\frac{14}{3}\frac{g^2}{l^2}{\omega }^2+\frac{15}{8}\frac{g^3}{l^3}}{\omega }^2{A}_1\end{array}}$

Calculamos los determinantes. En el código x es ω² y g es g/l

>> syms g x;
%denominador
>> A=[(5*g/2-7*x/3),-3*x/2,-x/2;-3*x/2,(3*g/2-4*x/3),-x/2; -x/2,-x/2,(g/2-x/3)]; 
>> det(A)
ans = (15*g^3)/8 - (14*g^2*x)/3 + (41*g*x^2)/24 - (13*x^3)/108
%numerdador A2
>> A=[5/2,-3*x/2,-x/2;3/2,(3*g/2-4*x/3),-x/2; 1/2,-x/2,(g/2-x/3)]; 
>> det(A)
ans = (15*g^2)/8 - (17*g*x)/12 + (11*x^2)/72 
%numerdador A3
>> A=[(5*g/2-7*x/3),5/2,-x/2;-3*x/2,3/2,-x/2; -x/2,1/2,(g/2-x/3)];
>> det(A) 
ans = (15*g^2)/8 - (g*x)/2 - x^2/24
%numerdador A4
>> A=[(5*g/2-7*x/3),-3*x/2,5/2;-3*x/2,(3*g/2-4*x/3),3/2; -x/2,-x/2,1/2];
>> det(A)
ans = (15*g^2)/8 + (g*x)/3 + x^2/72

Calculamos A₂, A₃ y A₄ para el primer modo de vibración de frecuencia angular ω=ω₁. El ángulo θ₁ que forma el primer péndulo con la vertical será proporcional a A₁, el ángulo θ₂ que forma el segundo péndulo con la vertical será proporcional a A₂, el ángulo θ₃ que forma el tercer péndulo con la vertical será proporcional a A₃ y el ángulo θ₄ que forma el cuarto péndulo con la vertical será proporcional a A₄. Representamos el primer modo de vibración.

L=1; %longitud total
n=4; %número de péndulos
w2=sort(roots([97/1296,-215*9.8/(108*L/n),943*9.8^2/(72*(L/n)^2),
-271*9.8^3/(12*(L/n)^3),105*9.8^4/(16*(L/n)^4)])); 
A1=0.2; %escala
iModo=1; %modo de vibración
den=-13*w2(iModo)^3/108+41*w2(iModo)^2*9.8/(24*L/n)-14*w2(iModo)*9.8^2
/(3*(L/n)^2)+15*9.8^3/(8*(L/n)^3);
A2=(11*w2(iModo)^2/72-17*w2(iModo)*9.8/(12*L/n)+15*9.8^2/(8*(L/n)^2))
*w2(iModo)*A1/den;
A3=(-w2(iModo)^2/24-w2(iModo)*9.8/(2*L/n)+15*9.8^2/(8*(L/n)^2))*w2(iModo)*A1/den;
A4=(w2(iModo)^2/72+w2(iModo)*9.8/(3*L/n)+15*9.8^2/(8*(L/n)^2))*w2(iModo)*A1/den;
angulos=[A1,A2, A3, A4];

hold on
x1=0; y1=0; %origen
plot(x1,y1,'ko','markersize',4,'markerfacecolor','k')
for j=1:n
    x2=x1+(L/n)*sin(angulos(j));
    y2=y1-(L/n)*cos(angulos(j));
    line([x1,x2],[y1,y2], 'lineWidth',1.5)
    if j==n
        break;
    end
    plot(x2,y2,'ro','markersize',3,'markerfacecolor','r')
    x1=x2;
    y1=y2;
end
text(0.1,-0.1, sprintf('frecuencia angular, %1.3f',sqrt(w2(iModo))))
hold off
axis equal
xlabel('x')
ylabel('y')
grid on
title(sprintf('Modo de vibración, %i',iModo))

Repetimos para el segundo modo de vibración de frecuencia angular ω=ω₂, cambiando el valor de la variable iModo=2

n péndulos

Para n= 4 péndulos, escribimos el sistema de ecuaciones diferenciales en forma matricial. A partir de este este ejemplo, escribiremos el código para calcular las frecuencias y representar los modos normales de vibración para un sistema de n péndulos

${\displaystyle \left(\begin{array}{cccc}10/3& 5/2& 3/2& 1/2\\ 5/2& 7/3& 3/2& 1/2\\ 3/2& 3/2& 4/3& 1/2\\ 1/2& 1/2& 1/2& 1/3\end{array}\right)\left(\begin{array}{c}\frac{d^2{\theta }_1}{d{t}^2}\\ \frac{d^2{\theta }_2}{d{t}^2}\\ \frac{d^2{\theta }_3}{d{t}^2}\\ \frac{d^2{\theta }_4}{d{t}^2}\end{array}\right)+\frac{g}{l}\left(\begin{array}{cccc}7/2& 0& 0& 0\\ 0& 5/2& 0& 0\\ 0& 0& 3/2& 0\\ 0& 0& 0& 1/2\end{array}\right)\left(\begin{array}{c}{\theta }_1\\ {\theta }_2\\ {\theta }_3\\ {\theta }_4\end{array}\right)=0}$

Valores propios

Se calculan las frecuencias de los modos normales de vibración resolviendo el sistema homogéneo

${\displaystyle \begin{array}{l}\begin{array}{l}-{\omega }^{2}M·x+K·x=0\hfill \\ \left(K-{\omega }^{2}M\right)·x=0\hfill \end{array}\\ \left(M^{-1}K-{\omega }^{2}I\right)x=0\end{array}}$

Los valores propios de la matriz M^-1·K son los cuadrados de las frecuencias ω² de los modos normales de vibración

Para n=4 péndulos de longitud total L=1.12 m

L=1; %longitud total
n=4; %número de péndulos

A=[4/3,1/2;1/2,1/3];
M=A;
for j=3:n
    M=zeros(j);
    M(2:j,2:j)=A;
    M(1,1)=((j-1)*3+1)/3;
    v=(2*j-3:-2:1)/2;
    M(2:j,1)=v;
    M(1,2:j)=v;
    A=M;
end

K=diag(2*n-1:-2:1)*9.8/(2*L/n);
[V,D]=eig(K/M);
w2=sort(diag(D)); %vector de frecuencias propias al cuadrado, ordenados

El cuadrado de las frecuencias de los modos normales de vibración, coinciden con las calculadas en el apartado anterior

>> w2
w2 =
   14.2201
   82.6924
  254.8446
  690.8821

Tomamos la segunda, tercera ... n ecuación del sistema homogéneo de n ecuaciones, para calcular las amplitudes A₂, A₃ ... A_n en términos de la amplitud A₁. En forma matricial se escribe para n=4, que es fácil generalizar para cualquier n

${\displaystyle \left\{-{\omega }^2\left(\begin{array}{ccc}\frac{7}{3}& \frac{3}{2}& \frac{1}{2}\\ \frac{3}{2}& \frac{4}{3}& \frac{1}{2}\\ \frac{1}{2}& \frac{1}{2}& \frac{1}{3}\end{array}\right)+\frac{g}{l}\left(\begin{array}{ccc}\frac{5}{2}& 0& 0\\ 0& \frac{3}{2}& 0\\ 0& 0& \frac{1}{2}\end{array}\right)\right\}\left(\begin{array}{c}{A}_2\\ A_3\\ A_4\end{array}\right)={\omega }^2\left(\begin{array}{c}\frac{5}{2}\\ \frac{3}{2}\\ \frac{1}{2}\end{array}\right)A_1}$

La primera matriz es la submatriz de M de 3×3 inferior derecha. Añadimos estas líneas de código para calcular los coeficienets A₂, A₃ y A₄

....
A=zeros(n-1);
for i=1:length(w2)
    C=-w2(i)*M(2:n,2:n)+(9.8/(L/n))*diag(n-1:-1:1); %incógnitas
    D=w2(i)*(n-1:-1:1)'; %términos independientes
    for j=1:n-1
        num=C;
        num(1:end,j)=D;
        A(i,j)=det(num)/det(C);
    end
end
....

>> A
A =
    1.2489    1.5264    1.8064  %A2, A3 y A4 para el primer modo w1
    0.6109   -0.8944   -3.7196  %A2, A3 y A4 para el segundo modo w2
   -0.5351   -1.9121    2.8219  %A2, A3 y A4 para el tercer modo w3
   -2.0574    1.6788   -1.0188  %A2, A3 y A4 para el cuarto modo w4

El ángulo θ₁ que forma el primer péndulo con la vertical será proporcional a A₁, el ángulo θ₂ que forma el segundo péndulo con la vertical será proporcional a A₂,... el ángulo θ_n que forma el péndulo n con la vertical será proporcional a A_n . Representamos el primer modo de vibración para n=8 péndulos. El código completo es

L=1; %longitud total
n=10; %número de péndulos

A=[4/3,1/2;1/2,1/3];
M=A;
for j=3:n
    M=zeros(j);
    M(2:j,2:j)=A;
    M(1,1)=((j-1)*3+1)/3;
    v=(2*j-3:-2:1)/2;
    M(2:j,1)=v;
    M(1,2:j)=v;
    A=M;
end

K=diag(2*n-1:-2:1)*9.8/(2*L/n);
[V,D]=eig(K/M);
w2=sort(diag(D)); %vector de frecuencias propias al cuadrado, ordenados

A=zeros(n-1);
for i=1:length(w2)
    C=-w2(i)*M(2:n,2:n)+(9.8/(L/n))*diag(2*n-3:-2:1)/2; %incógnitas
    D=(w2(i)*(2*n-3:-2:1)/2)'; %términos independientes
    for j=1:n-1
        num=C;
        num(1:end,j)=D;
        A(i,j)=det(num)/det(C);
    end
end
A1=0.2; %escala
iModo=1; %modo normal de vibración
angulos=[A1,A(iModo,1:end)*A1];

hold on
x1=0; y1=0; %origen
plot(x1,y1,'ko','markersize',4,'markerfacecolor','k')
for j=1:n
    x2=x1+(L/n)*sin(angulos(j));
    y2=y1-(L/n)*cos(angulos(j));
    line([x1,x2],[y1,y2], 'lineWidth',1.5)
    if j==n
        break;
    end
    plot(x2,y2,'ro','markersize',3,'markerfacecolor','r')
    x1=x2;
    y1=y2;
end
text(0.1,-0.1, sprintf('frecuencia angular, %1.3f',sqrt(w2(iModo))))
hold off
axis equal
xlabel('x')
ylabel('y')
grid on
title(sprintf('Modo de vibración, %i',iModo))

Repetimos para el segundo modo de vibración de frecuencia angular ω=ω₂, cambiando el valor de la variable iModo=2

Las frecuencias de los modos de vibración de una cuerda que cuelga son

${\displaystyle {\omega }_n=\frac{1}{2}{\xi }_n\sqrt{\frac{g}{L}}}$

Donde ξ_n son los ceros de de la función J₀(x) de Bessel

>> L=1;
>> x=linspace(0,20,20);
>> w2=raices(@(x) besselj(0,x),x)*sqrt(9.8/L)/2 %frecuencias 
w2 =    3.7641    8.6403   13.5452   18.4567   23.3706   28.2857

Para n=5,10,20, 40,... péndulos, las primeras frecuencias se van aproximado a las de los modos de vibración de una cuerda que cuelga (n=∞) de la misma longitud L