Modos normales de vibración (I)

Supongamos que la primera partícula se desplaza x₁ de la posición de equilibrio, que la segunda se desplaza x₂ y la tercera, se desplaza x_3.

En las figuras, se muestra las fuerzas sobre cada una de las partículas

Ecuación del movimiento de la primera partícula

$m_{1} \frac{d^{2} x_{1}}{d t^{2}} = - k_{1} x_{1} + k_{2} (x_{2} - x_{1})$

Ecuación del movimiento de la segunda partícula

$m_{2} \frac{d^{2} x_{2}}{d t^{2}} = - k_{2} (x_{2} - x_{1}) + k_{3} (x_{3} - x_{2})$

Ecuación del movimiento de la tercera partícula

$m_{3} \frac{d^{2} x_{3}}{d t^{2}} = - k_{3} (x_{3} - x_{2}) - k_{4} x_{3}$

En forma matricial

$\begin{array}{l} (\begin{matrix} m_{1} & 0 & 0 \\ 0 & m_{2} & 0 \\ 0 & 0 & m_{3} \end{matrix}) \frac{d^{2}}{d t^{2}} (\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \end{matrix}) + (\begin{matrix} (k_{1} + k_{2}) & - k_{2} & 0 \\ - k_{2} & (k_{2} + k_{3}) & - k_{3} \\ 0 & - k_{3} & (k_{3} + k_{4}) \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \end{matrix}) = 0 \\ M \frac{d^{2} x}{d t^{2}} + K x = 0 \end{array}$

En general, tendremos para n partículas unidas a n+1 muelles

$(\begin{matrix} m_{1} & 0 & 0 & 0 & ... & 0 \\ 0 & m_{2} & 0 & 0 & ... & 0 \\ 0 & 0 & m_{3} & 0 & ... & 0 \\ 0 & 0 & 0 & m_{4} & .. & 0 \\ ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & 0 & ... & m_{n} \end{matrix}) \frac{d^{2}}{d t^{2}} (\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ ... \\ x_{n} \end{matrix}) + (\begin{matrix} (k_{1} + k_{2}) & - k_{2} & 0 & 0 & ... & 0 \\ - k_{2} & (k_{2} + k_{3}) & - k_{3} & 0 & ... & 0 \\ 0 & - k_{3} & (k_{3} + k_{4}) & - k_{4} & ... & 0 \\ 0 & 0 & - k_{4} & (k_{4} + k_{5}) & .. & 0 \\ ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & 0 & ... & (k_{n} + k_{n + 1}) \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ ... \\ x_{n} \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ ... \\ 0 \end{matrix})$

Los modos normales de vibración y sus propiedades

Supongamos un sistema formado por tres partículas de la misma masa m unidas por muelles elásticos iguales de constante k, tal como se muestra en la figura. Este ejemplo, nos va a servir para elaborar un script que nos permita estudiar cualquier configuración de muelles y partículas.

En este caso, las tres ecuaciones diferenciales del movimiento se escriben

$\begin{array}{l} m \frac{d^{2} x_{1}}{d t^{2}} + 2 k x_{1} - k x_{2} = 0 \\ m \frac{d^{2} x_{2}}{d t^{2}} - k x_{1} + 2 k x_{2} - k x_{3} = 0 \\ m \frac{d^{2} x_{3}}{d t^{2}} - k x_{2} + 2 k x_{3} = 0 \end{array}$

En forma matricial

$\begin{array}{l} (\begin{matrix} m & 0 & 0 \\ 0 & m & 0 \\ 0 & 0 & m \end{matrix}) (\begin{matrix} \frac{d^{2} x_{1}}{d t^{2}} \\ \frac{d^{2} x_{2}}{d t^{2}} \\ \frac{d^{2} x_{3}}{d t^{2}} \end{matrix}) + (\begin{matrix} 2 k & - k & 0 \\ - k & 2 k & - k \\ 0 & - k & 2 k \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \end{matrix}) = 0 \\ M \ddot{x} + K x = 0 \end{array}$

La primera es la matriz M de las masas y la segunda es la matriz K de las constantes de los muelles. Fijamos el valor de la masa m=1 y de la constante k=1 de los muelles.

$M = (\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) K = (\begin{matrix} 2 & - 1 & 0 \\ - 1 & 2 & - 1 \\ 0 & - 1 & 2 \end{matrix})$

Valores y vectores propios

Buscamos una solución de la forma

x₁=X₁sin(ωt+φ₁), x₂=X₂sin(ωt+φ₂), x₃=X₃sin(ωt+φ₃)

que representan MAS de amplitud X₁, X₂, X₃ , frecuencia angular ω.

Las frecuencias de los modos normales de vibración se calculan resolviendo el sistema homogéneo

$\begin{array}{l} - ω^{2} M \cdot X + K \cdot X = 0 \\ (K - ω^{2} M) \cdot X = 0 \end{array}$ ;

Como la matriz M es diagonal su inversa M^-1 es tambien diagonal

$\begin{array}{l} M = (\begin{matrix} m_{1} & 0 & 0 \\ 0 & m_{2} & 0 \\ 0 & 0 & m_{3} \end{matrix}) M^{- 1} = (\begin{matrix} \frac{1}{m_{1}} & 0 & 0 \\ 0 & \frac{1}{m_{2}} & 0 \\ 0 & 0 & \frac{1}{m_{3}} \end{matrix}) \\ (M^{- 1} K - ω^{2} I) X = 0 \end{array}$

>> K=sym('[2,-1,0;-1,2,-1;0,-1,2]');
>> [V,D]=eig(K)
V =
[       1,        1, -1]
[ 2^(1/2), -2^(1/2),  0]
[       1,        1,  1]
D =
[ 2 - 2^(1/2),           0, 0]
[           0, 2^(1/2) + 2, 0]
[           0,           0, 2]

Los cuadrados de las frecuencias de los modos normales de vibración están en la diagonal de la matriz D. El vector propio correspondiente a cada uno de los valores propios son los vectores columna X⁽¹⁾, X⁽²⁾ y X⁽³⁾ de la matriz V, tal como se indica a continuación:

$\begin{array}{l} D = (\begin{matrix} ω_{1}^{2} & 0 & 0 \\ 0 & ω_{2}^{2} & 0 \\ 0 & 0 & ω_{3}^{2} \end{matrix}) {\begin{cases} ω_{1}^{2} = 2 - \sqrt{2} \\ ω_{2}^{2} = 2 + \sqrt{2} \\ ω_{3}^{2} = 2 \end{cases} \\ V = (\begin{matrix} X_{1}^{(1)} & X_{1}^{(2)} & X_{1}^{(3)} \\ X_{2}^{(1)} & X_{2}^{(2)} & X_{2}^{(3)} \\ X_{3}^{(1)} & X_{3}^{(2)} & X_{3}^{(3)} \end{matrix}) \\ X^{(1)} = (\begin{matrix} 1 \\ \sqrt{2} \\ 1 \end{matrix}) X^{(2)} = (\begin{matrix} 1 \\ - \sqrt{2} \\ 1 \end{matrix}) X^{(3)} = (\begin{matrix} - 1 \\ 0 \\ 1 \end{matrix}) \end{array}$

Propiedades de los vectores propios

Comprobamos las siguientes propiedades de los vectores propios X⁽¹⁾, X⁽²⁾ y X⁽³⁾

$\begin{array}{l} {(X^{(i)})}^{T} M X^{(j)} = 0 i \neq j \\ {(X^{(i)})}^{T} M X^{(i)} \neq 0 \\ {(X^{(i)})}^{T} K X^{(j)} = 0 i \neq j \\ {(X^{(i)})}^{T} K X^{(i)} \neq 0 \end{array}$

donde el símbolo T indica traspuesta (cambiar filas por columnas)

>> X1=V(:,1);
>> X2=V(:,2);
>> X3=V(:,3);

>> r=X1'*K*X2
r =4 - 2^(1/2)*(2*2^(1/2) - 2) - 2*2^(1/2)
>> simplify(r)
ans =0
>> r=X3'*K*X1 
r =0
>> r=X2'*K*X3
r =0
>> r=X2'*M*X1
r =0
>> r=X3'*M*X2
r =0

>> r=X1'*M*X1
r =4
>> r=X2'*M*X2
r =4
>> r=X3'*M*X3
r =2
>> r=X1'*K*X1
r =8 - 4*2^(1/2)
>> r=X2'*K*X2
r =4*2^(1/2) + 8
>> r=X3'*K*X3
r =4

El cociente

$\frac{{(X^{(i)})}^{T} K X^{(i)}}{{(X^{(i)})}^{T} M X^{(i)}} = ω_{i}^{2} i = 1,2,3$

es el cuadrado de las frecuencias de cada uno de los modos normales de vibración

Vamos a dividir los vectores X⁽¹⁾, X⁽²⁾ y X⁽³⁾ por un factor de escala de modo que

${(X^{(i)})}^{T} M X^{(i)} = 1$

De este modo, los cuadrados de las frecuencias de cada cada uno de los modos normales de vibración se calcularán a partir de la matriz K de las constantes de los muelles.

$ω_{i}^{2} = {(X^{(i)})}^{T} K X^{(i)}$

Como

$\begin{array}{l} {(X^{(1)})}^{T} M X^{(1)} = 4 \\ {(X^{(2)})}^{T} M X^{(2)} = 4 \\ {(X^{(3)})}^{T} M X^{(3)} = 2 \end{array}$

Al vector X⁽¹⁾ tendremos que dividirlo entre dos, X⁽²⁾ habrá que dividirlo entre dos y a X⁽³⁾entre raíz de dos

>> X1=X1/2;
>> X2=X2/2;
>> X3=X3/sym('sqrt(2)');

>> X1'*M*X1
ans =1
>> X2'*M*X2
ans =1 
>> X3'*M*X3
ans =1
 
>> X1'*K*X1 
ans =2 - 2^(1/2)
>> X2'*K*X2
ans =2^(1/2) + 2 
>> X3'*K*X3
ans =2

Los nuevos vectores X⁽¹⁾, X⁽²⁾ y X⁽³⁾son

$X^{(1)} = (\begin{matrix} \frac{1}{2} \\ \frac{\sqrt{2}}{2} \\ \frac{1}{2} \end{matrix}) X^{(2)} = (\begin{matrix} \frac{1}{2} \\ - \frac{\sqrt{2}}{2} \\ \frac{1}{2} \end{matrix}) X^{(3)} = (\begin{matrix} - \frac{\sqrt{2}}{2} \\ 0 \\ \frac{\sqrt{2}}{2} \end{matrix})$

y la nueva matriz V, formada por los vectores columna X⁽¹⁾, X⁽²⁾ y X⁽³⁾ es

$V = (\begin{matrix} \frac{1}{2} & \frac{1}{2} & - \frac{\sqrt{2}}{2} \\ \frac{\sqrt{2}}{2} & - \frac{\sqrt{2}}{2} & 0 \\ \frac{1}{2} & \frac{1}{2} & \frac{\sqrt{2}}{2} \end{matrix})$

Comprobamos que

$\begin{array}{l} M_{g} = {(V)}^{T} \cdot M \cdot V = (\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) \\ K_{g} = {(V)}^{T} \cdot K \cdot V = (\begin{matrix} ω_{1}^{2} & 0 & 0 \\ 0 & ω_{2}^{2} & 0 \\ 0 & 0 & ω_{3}^{2} \end{matrix}) \end{array}$

donde M_g y K_g la nueva matriz de masas y la nueva matriz diagonal de las constantes de los muelles. En la diagonal de esta matriz aparecen los cuadrados de las frecuencias de los modos normales de vibración

>> V=[X1,X2,X3] 
V =
[       1/2,        1/2, -2^(1/2)/2]
[ 2^(1/2)/2, -2^(1/2)/2,          0]
[       1/2,        1/2,  2^(1/2)/2]
 
>> V'*M*V
ans =
[ 1, 0, 0]
[ 0, 1, 0]
[ 0, 0, 1]
>> V'*K*V
ans = 
[ 2 - 2^(1/2),           0, 0]
[           0, 2^(1/2) + 2, 0]
[           0,           0, 2]

Superposición

Como vimos en la página anterior, el desplazamiento x_i(t) de cada una de las partículas es combinación la lineal

$\begin{array}{l} x_{1} (t) = X_{1}^{(1)} u_{1} (t) + X_{1}^{(2)} u_{2} (t) + X_{1}^{(3)} u_{3} (t) \\ x_{2} (t) = X_{2}^{(1)} u_{1} (t) + X_{2}^{(2)} u_{2} (t) + X_{2}^{(3)} u_{3} (t) \\ x_{3} (t) = X_{3}^{(1)} u_{1} (t) + X_{3}^{(2)} u_{2} (t) + X_{3}^{(3)} u_{3} (t) \end{array}$

En forma matricial escribimos

$x (t) = V \cdot u (t)$

Ecuaciones del movimiento desacopladas

En la ecuación diferencial del movimiento $M \ddot{x} + K x = 0$ reemplazamos

$\begin{array}{l} x (t) = V \cdot u (t) \\ \ddot{x} (t) = V \cdot \ddot{u} (t) \end{array}$

y multiplicamos por la matriz traspuesta de V.

$\begin{array}{l} {(V)}^{T} \cdot M \cdot V \cdot \ddot{u} + {(V)}^{T} \cdot K \cdot u = 0 \\ \ddot{u} + K_{g} \cdot u = 0 \\ (\begin{matrix} {\ddot{u}}_{1} \\ {\ddot{u}}_{2} \\ {\ddot{u}}_{3} \end{matrix}) + (\begin{matrix} ω_{1}^{2} & 0 & 0 \\ 0 & ω_{1}^{2} & 0 \\ 0 & 0 & ω_{1}^{2} \end{matrix}) (\begin{matrix} u_{1} \\ u_{2} \\ u_{3} \end{matrix}) = 0 {\begin{cases} \frac{d^{2} u_{1}}{d t^{2}} + ω_{1}^{2} u_{1} = 0 \\ \frac{d^{2} u_{2}}{d t^{2}} + ω_{2}^{2} u_{2} = 0 \\ \frac{d^{2} u_{3}}{d t^{2}} + ω_{3}^{2} u_{3} = 0 \end{cases} \end{array}$

Las nuevas ecuaciones del movimiento expresadas en términos de las coordenadas u(t) están desacopladas y su solución es bien conocida, un Movimiento Armónico Simple. Las soluciones de estas ecuaciones diferenciales son:

${\begin{cases} u_{1} (t) = A_{1} \cos (ω_{1} t) + B_{1} \sin (ω_{1} t) \\ u_{2} (t) = A_{2} \cos (ω_{2} t) + B_{2} \sin (ω_{2} t) \\ u_{3} (t) = A_{3} \cos (ω_{3} t) + B_{3} \sin (ω_{3} t) \end{cases}$

donde las constantes A_i y B_i se determinan a partir de las condiciones iniciales

${\begin{cases} u_{i} (0) = A_{i} \\ {\dot{u}}_{i} (0) = ω_{i} B_{i} \end{cases}$

Condiciones iniciales

Las condiciones inicales vienes determinadas por el desplazamiento de la posición de equilibrio y velocidad de cada una de las partículas en el instante t=0. Por ejemplo, desplazamos las partículas de sus posición de equilibrio y las soltamos

$t = 0 {\begin{cases} x_{1} = x_{01} \\ x_{2} = x_{02} \\ x_{3} = x_{03} \end{cases} {\begin{cases} \frac{d x_{1}}{d t} = 0 \\ \frac{d x_{2}}{d t} = 0 \\ \frac{d x_{3}}{d t} = 0 \end{cases}$

Teniendo en cuenta la relación entre x_i(t) y u_i(t)

$\begin{array}{l} x (0) = V \cdot u (0) \\ \dot{x} (0) = V \cdot \dot{u} (0) \\ u (0) = V^{- 1} x (0) \\ \dot{u} (0) = V^{- 1} \dot{x} (0) \end{array}$

En vez de calcular la matriz inversa de V, la sustituimos por un producto de matrices, que es más conveniente

$\begin{array}{l} {(V)}^{T} \cdot M \cdot V = I \\ V^{- 1} = {(V)}^{T} \cdot M \end{array}$

Las nuevas condiciones iniciales se escriben en forma matricial

$\begin{array}{l} u (0) = {(V)}^{T} \cdot M \cdot x (0) \\ \dot{u} (0) = {(V)}^{T} \cdot M \cdot \dot{x} (0) \end{array}$

Como ejemplo tomemos las condiciones iniciales siguientes

$x (0) = (\begin{matrix} 1 \\ 0 \\ 0 \end{matrix}) \dot{x} (0) = (\begin{matrix} 0 \\ 0 \\ 0 \end{matrix})$

>> x0=[1;0;0]; %posición
>> xp0=[0;0;0]; %velocidad (x prima)
>> u0=V'*M*x0   %posición transformada
u0 =
        1/2
        1/2
 -2^(1/2)/2
 
>> up0=V'*M*xp0  %velocidad (u prima) transformada 
up0 =
 0
 0
 0

A partir de las condiciones iniciales en el espacio u

${\begin{cases} u_{1} (0) = \frac{1}{2} \\ u_{2} (0) = \frac{1}{2} \\ u_{3} (0) = - \frac{\sqrt{2}}{2} \end{cases} {\begin{cases} {\dot{u}}_{1} (0) = 0 \\ {\dot{u}}_{2} (0) = 0 \\ {\dot{u}}_{3} (0) = 0 \end{cases}$

determinamos las constantes A_i y B_i

$\begin{array}{l} {\begin{cases} A_{1} = \frac{1}{2} \\ A_{2} = \frac{1}{2} \\ A_{3} = - \frac{\sqrt{2}}{2} \end{cases} {\begin{cases} B_{1} = 0 \\ B_{2} = 0 \\ B_{3} = 0 \end{cases} \\ u_{1} (t) = \frac{1}{2} \cos (ω_{1} t) \\ u_{2} (t) = \frac{1}{2} \cos (ω_{2} t) \\ u_{3} (t) = - \frac{\sqrt{2}}{2} \cos (ω_{3} t) \end{array}$

Finalmente, la posición de cada una de las partículas respecto del tiempo, el vector x se obtiene, multiplicando la matriz cuadrada V por el vector u

$\begin{array}{l} x (t) = V \cdot u (t) = (\begin{matrix} \frac{1}{2} & \frac{1}{2} & - \frac{\sqrt{2}}{2} \\ \frac{\sqrt{2}}{2} & - \frac{\sqrt{2}}{2} & 0 \\ \frac{1}{2} & \frac{1}{2} & \frac{\sqrt{2}}{2} \end{matrix}) (\begin{matrix} \frac{1}{2} \cos (ω_{1} t) \\ \frac{1}{2} \cos (ω_{2} t) \\ - \frac{\sqrt{2}}{2} \cos (ω_{3} t) \end{matrix}) = \\ (\begin{matrix} \frac{1}{4} \cos (ω_{1} t) + \frac{1}{4} \cos (ω_{2} t) + \frac{1}{2} \cos (ω_{3} t) \\ \frac{\sqrt{2}}{4} \cos (ω_{1} t) - \frac{\sqrt{2}}{4} \cos (ω_{2} t) \\ \frac{1}{4} \cos (ω_{1} t) + \frac{1}{4} \cos (ω_{2} t) - \frac{1}{2} \cos (ω_{3} t) \end{matrix}) \end{array}$

>> syms t;
>> u1=cos(sqrt(D(1,1))*t)/2;
>> u2=cos(sqrt(D(2,2))*t)/2;
>> u3=-cos(sqrt(D(3,3))*t)*sqrt(2)/2;
>> u=[u1;u2;u3];
>> x=V*u
 cos(t*(2^(1/2) + 2)^(1/2))/4 + cos(2^(1/2)*t)/2 + 
    cos(t*(2 - 2^(1/2))^(1/2))/4
    -(2^(1/2)*(cos(t*(2^(1/2) + 2)^(1/2)) -
    cos(t*(2 - 2^(1/2))^(1/2))))/4
    cos(t*(2^(1/2) + 2)^(1/2))/4 - cos(2^(1/2)*t)/2 + 
    cos(t*(2 - 2^(1/2))^(1/2))/4
>> hold on
>> ezplot(x(1),[0 20])
>> g=ezplot(x(2),[0 20]);
>> set(g,'color','red')
>> g=ezplot(x(3),[0 20]);
>> set(g,'color','green')
>> grid on

Representación gráfica del movimiento de cada partícula

Solamente nos queda unir las distintas porciones de código MATLAB en un script que nos permita resolver el sistema formado por tres partículas y optimizar el código para que pueda ser generalizable a cualquier configuración

syms t;
%matriz de las constantes de los muelles
K=sym('[2,-1,0;-1,2,-1;0,-1,2]'); 
M=sym('[1,0,0;0,1,0;0,0,1]'); %matriz de las masas

[V,D]=eig(inv(M)*K); %valores y vectores propios
w=diag(sqrt(D)); %vector de frecuencias propias
for i=1:length(w) %para hacer diagonal K 
    r=V(:,i)'*M*V(:,i);
    V(:,i)=V(:,i)/sqrt(r);
end

%condiciones iniciales
x0=sym('[1;0;0]'); %posición inicial
xp0=sym('[0;0;0]'); %velocidad (x prima) inicial
u0=V'*M*x0; 
up0=V'*M*xp0;

%movimiento de cada una de las partículas 
A=u0;
B=up0./diag(D);
u=diag(A)*cos(w*t)+diag(B)*sin(w*t);
x=V*u;
x=simplify(x)

%representación gráfica
color=['b','r','g'];
hold on
for i=1:length(w)
    h=ezplot(x(i),[0,20]);
    set(h,'color',color(i))
end
title('Tres osciladores acoplados')
ylabel('x_1,x_2,x_3')
xlabel('t')
grid on
hold off

En la ventana de comandos obtenemos las frecuencias de los modos normales de vibración y el vector x(t), posición de cada partícula en función del tiempo, cuya representación gráfica se muestra en la figura: la primera en color azul, la segunda en color rojo y la tercera en color verde..

w =
 (2 - 2^(1/2))^(1/2)
 (2^(1/2) + 2)^(1/2)
             2^(1/2)
x = cos(t*(2^(1/2) + 2)^(1/2))/4 + cos(2^(1/2)*t)/2 
    + cos(t*(2 - 2^(1/2))^(1/2))/4
         -(2^(1/2)*(cos(t*(2^(1/2) + 2)^(1/2)) - 
         cos(t*(2 - 2^(1/2))^(1/2))))/4
 cos(t*(2^(1/2) + 2)^(1/2))/4 - cos(2^(1/2)*t)/2 
 + cos(t*(2 - 2^(1/2))^(1/2))/4

Cambiamos la matriz K de las constantes de los muelles elásticos, la matriz M de las masas y los vectores posición inicial x0 y velocidad inicial xp0. Si se incrementa la dimensión de la matrices hay que añadir más colores al vector color.

El siguiente script es similar al anterior pero que funciona con la versión básica de MATLAB mientras que en el anterior es necesario disponer del complemento Math Symbolic por lo que se aconseja su utilización, salvo en aquellos casos en los que sea preciso obtener una respuesta analítica exacta.

K=[2,-1,0;-1,2,-1;0,-1,2]; %matriz constante de los muelles
M=[1,0,0;0,1,0;0,0,1]; %matriz de las masas

[V,D]=eig(inv(M)*K); %valores y vectores propios
w=diag(sqrt(D)) %vector de frecuencias propias
for i=1:length(w)
    r=V(:,i)'*M*V(:,i);
    V(:,i)=V(:,i)/sqrt(r);
end

%condiciones iniciales
x0=[1;0;0]; %posición inicial
xp0=[0;0;0]; %velocidad (x prima) inicial
u0=V'*M*x0;
up0=V'*M*xp0;

%movimiento de cada una de las partículas 
A=u0;
B=up0./diag(D);
t=linspace(0,20,200);
u=diag(A)*cos(w*t)+diag(B)*sin(w*t);
x=V*u;

%representación gráfica
plot(t,x);
title('Tres osciladores acoplados')
ylabel('x_1,x_2,x_3')
xlabel('t')
grid on

Obtenemos una gráfica similar, pero los colores de las representación gráfica de la posición x_i(t) de cada partícula lo establece MATLAB.

En la ventana de comandos obtenemos las frecuencias de los modos normales de vibración

Frecuencias de los modos normales de vibración

Sistema formado por dos partículas

Las ecuaciones del movimiento son

$\begin{array}{l} m \frac{d^{2} x_{1}}{d t^{2}} = - k x_{1} + k (x_{2} - x_{1}) \\ m \frac{d^{2} x_{2}}{d t^{2}} = - k (x_{2} - x_{1}) - k x_{2} \end{array}$

expresadas en forma matricial

$(\begin{matrix} m & 0 \\ 0 & m \end{matrix}) (\begin{matrix} \frac{d^{2} x_{1}}{d t^{2}} \\ \frac{d^{2} x_{2}}{d t^{2}} \end{matrix}) + (\begin{matrix} 2 k & - k \\ - k & 2 k \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \end{matrix}) = 0$

La primera es la matriz M de las masas y la segunda es la matriz K de las constantes de los muelles. Las frecuencias de los modos normales de vibración ω² son los valores propios de la matriz M^-1K. El vector propio correspondiente a cada uno de los valores propios son los vectores columna de la matriz V

>> syms k m;
>> M=[m,0;0,m];
>> K=[2*k,-k;-k,2*k];
>> [V,D]=eig(M^-1*K);
>> w2 =diag(D)
w2 =
     k/m
 (3*k)/m
 >> V
V =
[ 1, -1]
[ 1,  1]

Valores propios ω²	Vectores propios
$\frac{k}{m}$	(1, 1)
$3 \frac{k}{m}$	(-1,1)

Sistema formado por tres partículas

Las ecuaciones del movimiento son

$\begin{array}{l} m \frac{d^{2} x_{1}}{d t^{2}} = - k x_{1} + k (x_{2} - x_{1}) \\ m \frac{d^{2} x_{2}}{d t^{2}} = - k (x_{2} - x_{1}) + k (x_{3} - x_{2}) \\ m \frac{d^{2} x_{3}}{d t^{2}} = - k (x_{3} - x_{2}) - k x_{3} \end{array}$

expresada en forma matricial

$(\begin{matrix} m & 0 & 0 \\ 0 & m & 0 \\ 0 & 0 & m \end{matrix}) (\begin{matrix} \frac{d^{2} x_{1}}{d t^{2}} \\ \frac{d^{2} x_{2}}{d t^{2}} \\ \frac{d^{2} x_{3}}{d t^{2}} \end{matrix}) + (\begin{matrix} 2 k & - k & 0 \\ - k & 2 k & - k \\ 0 & - k & 2 k \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \end{matrix}) = 0$

>> syms k m;
>> M=[m,0,0;0,m,0; 0,0,m];
>> K=[2*k,-k,0;-k,2*k,-k;0,-k,2*k];
>> [V,D]=eig(M^-1*K);
>> w2 =diag(D)
w2 =
  (k*(2^(1/2) + 2))/m
 -(k*(2^(1/2) - 2))/m
              (2*k)/m
>> V
V =
[        1,       1, -1]
[ -2^(1/2), 2^(1/2),  0]
[        1,       1,  1]

Valores propios ω²	Vectores propios
$(2 + \sqrt{2}) \frac{k}{m}$	$(1, - \sqrt{2}, 1)$
$(2 - \sqrt{2}) \frac{k}{m}$	$(1, \sqrt{2},1)$
$2 \frac{k}{m}$	(-1,0,1)

Sistema formado por cuatro partículas

Las ecuaciones del movimiento son

expresadas en forma matricial

$(\begin{matrix} m & 0 & 0 & 0 \\ 0 & m & 0 & 0 \\ 0 & 0 & m & 0 \\ 0 & 0 & 0 & m \end{matrix}) (\begin{matrix} \frac{d^{2} x_{1}}{d t^{2}} \\ \frac{d^{2} x_{2}}{d t^{2}} \\ \frac{d^{2} x_{3}}{d t^{2}} \\ \frac{d^{2} x_{4}}{d t^{2}} \end{matrix}) + (\begin{matrix} 2 k & - k & 0 & 0 \\ - k & 2 k & - k & 0 \\ 0 & - k & 2 k & - k \\ 0 & 0 & - k & 2 k \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \end{matrix}) = 0$

>> syms k m;
>> M=[m,0,0,0;0,m,0,0; 0,0,m,0; 0,0,0,m];
>> K=[2*k,-k,0,0;-k,2*k,-k,0;0,-k,2*k,-k;0,0,-k,2*k];
>> [V,D]=eig(M^-1*K);
>> w2 =diag(D)
w2 =
   (k*(5^(1/2) + 3))/(2*m)
 -(k*(5^(1/2) - 3))/(2*m)
  (k*(5^(1/2) + 5))/(2*m)
 -(k*(5^(1/2) - 5))/(2*m)
>> V
V =
[               1,               1,                -1,              -1]
[ 1/2 - 5^(1/2)/2, 5^(1/2)/2 + 1/2,   5^(1/2)/2 + 1/2, 1/2 - 5^(1/2)/2]
[ 1/2 - 5^(1/2)/2, 5^(1/2)/2 + 1/2, - 5^(1/2)/2 - 1/2, 5^(1/2)/2 - 1/2]
[               1,               1,                 1,               1]

Valores propios ω²	Vectores propios
$(\frac{3 + \sqrt{5}}{2}) \frac{k}{m}$	$(1, \frac{1 - \sqrt{5}}{2}, \frac{1 - \sqrt{5}}{2},1)$
$(\frac{5 - \sqrt{5}}{2}) \frac{k}{m}$	$(- 1, \frac{1 - \sqrt{5}}{2}, - \frac{1 - \sqrt{5}}{2},1)$
$(\frac{3 - \sqrt{5}}{2}) \frac{k}{m}$	$(1, \frac{1 + \sqrt{5}}{2}, \frac{1 + \sqrt{5}}{2},1)$
$(\frac{5 + \sqrt{5}}{2}) \frac{k}{m}$	$(- 1, \frac{1 + \sqrt{5}}{2}, \frac{1 - \sqrt{5}}{2},1)$

Actividades

Se ha fijado

La masa de las partículas, m=1 kg
La constante elástica de los muelles, k=1 N/m

Se introduce

El número N de partículas, en el control titulado Partículas

Se pulsa el botón titulado Nuevo

Observamos el primer modo normal de vibración, en la parte superior izquierda se proporciona la frecuencia angular.

Se pulsa el botón titulado >>, observamos los siguientes modos normales

Se pulsa el botón titulado << para observar el modo normal de vibración anterior.

En la parte inferior, se muestra el desplazamiento de las partículas

Partículas:

Ejemplos

Probamos con los ejemplos estudiados en la página "Modos normales de vibración":

Ejemplo 1

El sistema está formado por dos partículas de masa m=1. La constante de los muelles es k_1,3=10 y la constante del acoplamiento k₂=0.5. Se desvía la primera partícula x₀₁=1 de la posición de eqilibrio y se suelta. La posición inicial de la segunda partícula x₀₂=0, velocidad inicial de ambas v₀₁=0, v₀₂=0.

k1=10; k2=0.5; k3=10;
m1=1; m2=1;
K=[(k1+k2),-k2;-k2,(k2+k3)]; %sistema de dos partículas
M=[m1,0;0,m2];
....
%condiciones iniciales
x0=[1;0];
xp0=[0;0];
....

En la ventana de comandos obtenemos las frecuencias de los modos normales de vibración

w =
    3.1623
    3.3166

Utilizamos alternativamente, el script de cálculo simbólico

syms k1 k2 k3 m1 m2;
K=[(k1+k2),-k2;-k2,(k2+k3)];
M=[m1,0;0,m2];
K=subs(K,{k1,k2,k3},{10,0.5,10});
M=subs(M,{m1,m2},{1,1});
....
%condiciones iniciales
x0=sym('[1;0]');
xp0=sym('[0;0]');
.....

En la ventana de comandos obtenemos la expresión de las posiciones de las partículas en función del tiempo, x₁(t) y x₂(t). Los argumentos de las funciones coseno son las frecuencias de los modos normales de vibración

$ω_{1} = \sqrt{10} ω_{2} = \sqrt{11}$

x =
 cos(10^(1/2)*t)/2 + cos(11^(1/2)*t)/2
 cos(10^(1/2)*t)/2 - cos(11^(1/2)*t)/2

Ejemplo 2

El sistema formado por dos partículas de masas m₁=10 y m₂=1 unidas por muelles de constantes k₁=30, k₂=5 y k₃=0. Las condiciones iniciales en el instante t=0, son las siguientes: posición inicial x₀₁=1, x₀₂=0, velocidad inicial v₀₁=0, v₀₂=0.

k1=30; k2=5; k3=0;m1=10; m2=1;
K=[(k1+k2),-k2;-k2,(k2+k3)];
M=[m1,0;0,m2];
....
%condiciones iniciales
x0=[1;0];
xp0=[0;0];
...

En la ventana de comandos obtenemos las frecuencias de los modos normales de vibración

w =
    1.5811
    2.4495

Utilizamos alternativamente, el script de cálculo simbólico

syms k1 k2 k3 m1 m2;
K=[(k1+k2),-k2;-k2,(k2+k3)]; %para un sistema de dos partículas
M=[m1,0;0,m2];
K=subs(K,{k1,k2,k3},{30,5,0});
M=subs(M,{m1,m2},{10,1});
....
%condiciones iniciales
x0=sym('[1;0]');
xp0=sym('[0;0]');
.....

$ω_{1} = \sqrt{6} ω_{2} = \frac{\sqrt{10}}{2}$

X =
   (2*cos(6^(1/2)*t))/7 + (5*cos((10^(1/2)*t)/2))/7
 (10*cos((10^(1/2)*t)/2))/7 - (10*cos(6^(1/2)*t))/7

Ejemplo 3

El sistema formado por dos partículas de masas m₁=2 y m₂=1 unidas por muelles de constantes k₁=6, k₂=3 y k₃=0. Las condiciones iniciales en el instante t=0, son las siguientes: posición inicial x₀₁=0, x₀₂=0, velocidad inicial v₀₁=1, v₀₂=0.

En la ventana de comandos obtenemos las frecuencias de los modos normales de vibración

w =
    2.4495
    1.2247