MateWiki - Contribuciones del usuario [es]

Sistema resorte-masa

2013-03-08T15:59:45Z

Angelcf:

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot u=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h*k3x)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h*k3x)-((k2+k3)/m2)*(yy(n)+k3y*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

t=t0:h:tN;
axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;

for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2y;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3y;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
end

% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
}}

[[Archivo:Run4.png|marco|centro]]

===Energía del sistema===
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::

<math>E=\frac{1}{2}m_{1}\dot x ^2+ \frac{1}{2}k_{1}x^2</math>:

[[Archivo:Enerlog.png|marco|centro]]

===Influencia del coeficiente <math>\mu</math>===

A continuación se muestran una serie de casos con valores de <math>\mu</math> comprendidos entre 0 y 5. Las gráficas muestran las trayectorias y el decrecimiento de la energía para cada uno de los casos:

'''<big>Caso 1 <math>\mu=0</math></big>'''
Caso limite sin amortiguamiento en el que la energía es constante
[[Archivo:4c0.png|marco|centro]]
'''<big>Caso 2 <math>\mu=0.5</math></big>'''
[[Archivo:4c05.png|marco|centro]]
'''<big>Caso 3 <math>\mu=1</math></big>'''
[[Archivo:4c1.png|marco|centro]]
'''<big>Caso 4 <math>\mu=2</math></big>'''
[[Archivo:4c2.png|marco|centro]]
'''<big>Caso 5 <math>\mu=3</math></big>'''
[[Archivo:4c3.png|marco|centro]]
'''<big>Caso 6 <math>\mu=5</math></big>'''
[[Archivo:4c5.png|marco|centro]]

La conclusión que se extrae en vista de las gráficas es que la influencia del coeficiente de amortiguamiento <math>\mu</math> es muy clara, cuanto mayor es este coeficiente más rapidamente decrece la energía pues la masa se estabiliza antes en torno a la posición de equilibrio.

[[Usuario:Angelcf|Angelcf]] ([[Usuario discusión:Angelcf|discusión]]) 18:39 1 mar 2013 (CET)

Sistema resorte-masa

2013-03-03T18:35:28Z

Angelcf: /* Método de Runge-Kutta */

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot u=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3y*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

t=t0:h:tN;
axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;

for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2y;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3y;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
end

% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
}}

[[Archivo:Run4.png|marco|centro]]

===Energía del sistema===
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::

<math>E=\frac{1}{2}m_{1}\dot x ^2+ \frac{1}{2}k_{1}x^2</math>:

[[Archivo:Enerlog.png|marco|centro]]

===Influencia del coeficiente <math>\mu</math>===

A continuación se muestran una serie de casos con valores de <math>\mu</math> comprendidos entre 0 y 5. Las gráficas muestran las trayectorias y el decrecimiento de la energía para cada uno de los casos:

'''<big>Caso 1 <math>\mu=0</math></big>'''
Caso limite sin amortiguamiento en el que la energía es constante
[[Archivo:4c0.png|marco|centro]]
'''<big>Caso 2 <math>\mu=0.5</math></big>'''
[[Archivo:4c05.png|marco|centro]]
'''<big>Caso 3 <math>\mu=1</math></big>'''
[[Archivo:4c1.png|marco|centro]]
'''<big>Caso 4 <math>\mu=2</math></big>'''
[[Archivo:4c2.png|marco|centro]]
'''<big>Caso 5 <math>\mu=3</math></big>'''
[[Archivo:4c3.png|marco|centro]]
'''<big>Caso 6 <math>\mu=5</math></big>'''
[[Archivo:4c5.png|marco|centro]]

La conclusión que se extrae en vista de las gráficas es que la influencia del coeficiente de amortiguamiento <math>\mu</math> es muy clara, cuanto mayor es este coeficiente más rapidamente decrece la energía pues la masa se estabiliza antes en torno a la posición de equilibrio.

[[Usuario:Angelcf|Angelcf]] ([[Usuario discusión:Angelcf|discusión]]) 18:39 1 mar 2013 (CET)

Sistema resorte-masa

2013-03-03T18:04:17Z

Angelcf: /* Método de Runge-Kutta */

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot u=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3y*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;

for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2y;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3y;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
end

% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
}}

[[Archivo:Run4.png|marco|centro]]

===Energía del sistema===
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::

<math>E=\frac{1}{2}m_{1}\dot x ^2+ \frac{1}{2}k_{1}x^2</math>:

[[Archivo:Enerlog.png|marco|centro]]

===Influencia del coeficiente <math>\mu</math>===

A continuación se muestran una serie de casos con valores de <math>\mu</math> comprendidos entre 0 y 5. Las gráficas muestran las trayectorias y el decrecimiento de la energía para cada uno de los casos:

'''<big>Caso 1 <math>\mu=0</math></big>'''
Caso limite sin amortiguamiento en el que la energía es constante
[[Archivo:4c0.png|marco|centro]]
'''<big>Caso 2 <math>\mu=0.5</math></big>'''
[[Archivo:4c05.png|marco|centro]]
'''<big>Caso 3 <math>\mu=1</math></big>'''
[[Archivo:4c1.png|marco|centro]]
'''<big>Caso 4 <math>\mu=2</math></big>'''
[[Archivo:4c2.png|marco|centro]]
'''<big>Caso 5 <math>\mu=3</math></big>'''
[[Archivo:4c3.png|marco|centro]]
'''<big>Caso 6 <math>\mu=5</math></big>'''
[[Archivo:4c5.png|marco|centro]]

La conclusión que se extrae en vista de las gráficas es que la influencia del coeficiente de amortiguamiento <math>\mu</math> es muy clara, cuanto mayor es este coeficiente más rapidamente decrece la energía pues la masa se estabiliza antes en torno a la posición de equilibrio.

[[Usuario:Angelcf|Angelcf]] ([[Usuario discusión:Angelcf|discusión]]) 18:39 1 mar 2013 (CET)

Sistema resorte-masa

2013-03-03T18:02:52Z

Angelcf:

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot u=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3y*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;

for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2x;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3x;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
end

% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
}}

[[Archivo:Run4.png|marco|centro]]

===Energía del sistema===
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::

<math>E=\frac{1}{2}m_{1}\dot x ^2+ \frac{1}{2}k_{1}x^2</math>:

[[Archivo:Enerlog.png|marco|centro]]

===Influencia del coeficiente <math>\mu</math>===

A continuación se muestran una serie de casos con valores de <math>\mu</math> comprendidos entre 0 y 5. Las gráficas muestran las trayectorias y el decrecimiento de la energía para cada uno de los casos:

'''<big>Caso 1 <math>\mu=0</math></big>'''
Caso limite sin amortiguamiento en el que la energía es constante
[[Archivo:4c0.png|marco|centro]]
'''<big>Caso 2 <math>\mu=0.5</math></big>'''
[[Archivo:4c05.png|marco|centro]]
'''<big>Caso 3 <math>\mu=1</math></big>'''
[[Archivo:4c1.png|marco|centro]]
'''<big>Caso 4 <math>\mu=2</math></big>'''
[[Archivo:4c2.png|marco|centro]]
'''<big>Caso 5 <math>\mu=3</math></big>'''
[[Archivo:4c3.png|marco|centro]]
'''<big>Caso 6 <math>\mu=5</math></big>'''
[[Archivo:4c5.png|marco|centro]]

La conclusión que se extrae en vista de las gráficas es que la influencia del coeficiente de amortiguamiento <math>\mu</math> es muy clara, cuanto mayor es este coeficiente más rapidamente decrece la energía pues la masa se estabiliza antes en torno a la posición de equilibrio.

[[Usuario:Angelcf|Angelcf]] ([[Usuario discusión:Angelcf|discusión]]) 18:39 1 mar 2013 (CET)

Sistema resorte-masa

2013-03-03T11:15:58Z

Angelcf:

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot u=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3y*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>:
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>:

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;

for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2x;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3x;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
end

% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
}}

[[Archivo:Run4.png|marco|centro]]

===Energía del sistema===
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::

<math>E=\frac{1}{2}m_{1}\dot x ^2+ \frac{1}{2}k_{1}x^2</math>:

[[Archivo:Enerlog.png|marco|centro]]

===Influencia del coeficiente <math>\mu</math>===

A continuación se muestran una serie de casos con valores de <math>\mu</math> comprendidos entre 0 y 5. Las gráficas muestran las trayectorias y el decrecimiento de la energía para cada uno de los casos:

'''<big>Caso 1 <math>\mu=0</math></big>'''
Caso limite sin amortiguamiento en el que la energía es constante
[[Archivo:4c0.png|marco|centro]]
'''<big>Caso 2 <math>\mu=0.5</math></big>'''
[[Archivo:4c05.png|marco|centro]]
'''<big>Caso 3 <math>\mu=1</math></big>'''
[[Archivo:4c1.png|marco|centro]]
'''<big>Caso 4 <math>\mu=2</math></big>'''
[[Archivo:4c2.png|marco|centro]]
'''<big>Caso 5 <math>\mu=3</math></big>'''
[[Archivo:4c3.png|marco|centro]]
'''<big>Caso 6 <math>\mu=5</math></big>'''
[[Archivo:4c5.png|marco|centro]]

La conclusión que se extrae en vista de las gráficas es que la influencia del coeficiente de amortiguamiento <math>\mu</math> es muy clara, cuanto mayor es este coeficiente más rapidamente decrece la energía pues la masa se estabiliza antes en torno a la posición de equilibrio.

[[Usuario:Angelcf|Angelcf]] ([[Usuario discusión:Angelcf|discusión]]) 18:39 1 mar 2013 (CET)

Sistema resorte-masa

2013-03-02T18:32:07Z

Angelcf:

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3y*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>:
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>:

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;

for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2x;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3x;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
end

% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
}}

[[Archivo:Run4.png|marco|centro]]

===Energía del sistema===
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::

<math>E=\frac{1}{2}m_{1}\dot x ^2+ \frac{1}{2}k_{1}x^2</math>:

[[Archivo:Enerlog.png|marco|centro]]

===Influencia del coeficiente <math>\mu</math>===

A continuación se muestran una serie de casos con valores de \mu comprendidos entre 0 y 5. Las gráficas muestran las trayectorias y el decrecimiento de la energía para cada uno de los casos:

'''<big>Caso 1 <math>\mu=0</math></big>'''
Caso limite sin amortiguamiento en el que la energía es constante
[[Archivo:4c0.png|marco|centro]]
'''<big>Caso 2 <math>\mu=0.5</math></big>'''
[[Archivo:4c05.png|marco|centro]]
'''<big>Caso 3 <math>\mu=1</math></big>'''
[[Archivo:4c1.png|marco|centro]]
'''<big>Caso 4 <math>\mu=2</math></big>'''
[[Archivo:4c2.png|marco|centro]]
'''<big>Caso 5 <math>\mu=3</math></big>'''
[[Archivo:4c3.png|marco|centro]]
'''<big>Caso 6 <math>\mu=5</math></big>'''
[[Archivo:4c5.png|marco|centro]]

La conclusión que se extrae en vista de las gráficas es que la influencia del coeficiente de amortiguamiento <math>\mu</math> es muy clara, cuanto mayor es este coeficiente más rapidamente decrece la energía pues la masa se estabiliza antes en torno a la posición de equilibrio.

[[Usuario:Angelcf|Angelcf]] ([[Usuario discusión:Angelcf|discusión]]) 18:39 1 mar 2013 (CET)

Sistema resorte-masa

2013-03-01T17:39:43Z

Angelcf:

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3y*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>:
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>:

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;

for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2x;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3x;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
end

% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
}}

[[Archivo:Run4.png|marco|centro]]

===Energía del sistema===
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::

<math>E=\frac{1}{2}m_{1}\dot x ^2+ \frac{1}{2}k_{1}x^2</math>:

[[Archivo:Enerlog.png|marco|centro]]

===Influencia del coeficiente <math>\mu</math>===

A continuación se muestran una serie de casos con valores de \mu comprendidos entre 0 y 5. Las gráficas muestran las trayectorias y el decrecimiento de la energía para cada uno de los casos:

'''<big>Caso 1 <math>\mu=0</math></big>'''
Caso limite sin amortiguamiento en el que la energía es constante
[[Archivo:4c0.png|marco|centro]]
'''<big>Caso 2 <math>\mu=0.5</math></big>'''
[[Archivo:4c05.png|marco|centro]]
'''<big>Caso 3 <math>\mu=1</math></big>'''
[[Archivo:4c1.png|marco|centro]]
'''<big>Caso 4 <math>\mu=2</math></big>'''
[[Archivo:4c2.png|marco|centro]]
'''<big>Caso 5 <math>\mu=3</math></big>'''
[[Archivo:4c3.png|marco|centro]]
'''<big>Caso 6 <math>\mu=5</math></big>'''
[[Archivo:4c5.png|marco|centro]]

La conclusión que se extrae en vista de las gráficas es que la influencia del coeficiente de amortiguamiento \mu es muy clara, cuanto mayor es este coeficiente más rapidamente decrece la energía pues la masa se estabiliza antes en torno a la posición de equilibrio.

[[Usuario:Angelcf|Angelcf]] ([[Usuario discusión:Angelcf|discusión]]) 18:39 1 mar 2013 (CET)

Sistema resorte-masa

2013-03-01T17:30:42Z

Angelcf:

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3y*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>:
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>:

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;

for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2x;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3x;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
end

% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
}}

[[Archivo:Run4.png|marco|centro]]

===Energía del sistema===
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::

<math>E=\frac{1}{2}m_{1}\dot x ^2+ \frac{1}{2}k_{1}x^2</math>:

[[Archivo:Enerlog.png|marco|centro]]

===Influencia del coeficiente <math>\mu</math>===

A continuación se muestran una serie de casos con valores de \mu comprendidos entre 0 y 5. Las gráficas muestran las trayectorias y el decrecimiento de la energía para cada uno de los casos:

'''<big>Caso 1 <math>\mu=0</math></big>'''
Caso limite sin amortiguamiento en el que la energía es constante
[[Archivo:4c0.png|marco|centro]]
'''<big>Caso 2 <math>\mu=0.5</math></big>'''
[[Archivo:4c05.png|marco|centro]]
'''<big>Caso 3 <math>\mu=1</math></big>'''
[[Archivo:4c1.png|marco|centro]]
'''<big>Caso 4 <math>\mu=2</math></big>'''
[[Archivo:4c2.png|marco|centro]]
'''<big>Caso 5 <math>\mu=3</math></big>'''
[[Archivo:4c3.png|marco|centro]]
'''<big>Caso 6 <math>\mu=5</math></big>'''
[[Archivo:4c5.png|marco|centro]]

La conclusión que se extrae en vista de las gráficas es que la influencia del coeficiente de amortiguamiento \mu es muy clara, cuanto mayor es este coeficiente más rapidamente decrece la energía pues la masa se estabiliza antes en torno a la posición de equilibrio.

Archivo:4c5.png

2013-03-01T17:25:42Z

Angelcf:

Archivo:4c3.png

2013-03-01T17:24:23Z

Angelcf:

Archivo:4c2.png

2013-03-01T17:23:47Z

Angelcf:

Archivo:4c1.png

2013-03-01T17:23:17Z

Angelcf:

Archivo:4c05.png

2013-03-01T17:22:48Z

Angelcf:

Archivo:4c0.png

2013-03-01T17:21:57Z

Angelcf:

Sistema resorte-masa

2013-03-01T17:02:47Z

Angelcf:

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3y*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>:
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>:

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;

for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2x;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3x;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
end

% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
}}

[[Archivo:Run4.png|marco|centro]]

===Energía del sistema===
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::

<math>E=\frac{1}{2}m_{1}\dot x ^2+ \frac{1}{2}k_{1}x^2</math>:

[[Archivo:Enerlog.png|marco|centro]]

El

Sistema resorte-masa

2013-03-01T17:02:09Z

Angelcf:

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3y*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>:
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>:

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;

for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2x;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3x;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
end

% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
}}

[[Archivo:Run4.png|marco|centro]]

===Energía del sistema===
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::

<math>E=\frac{1}{2}m_{1}(\dot x)^2+ \frac{1}{2}k_{1}x^2</math>:

[[Archivo:Enerlog.png|marco|centro]]

El

Archivo:Enerlog.png

2013-03-01T17:01:44Z

Angelcf:

Archivo:Ener4.png

2013-03-01T16:45:06Z

Angelcf:

Archivo:Run4.png

2013-03-01T16:35:36Z

Angelcf:

Sistema resorte-masa

2013-03-01T16:14:47Z

Angelcf:

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
==Introducción==
El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.
[[Archivo:Muellemasa.png|marco|centro]]

==Ecuaciones diferenciales==
La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante '''m1''' y '''m2'''). Los grados de libertad usados son '''x''' e '''y''', siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

==Resolucion==
Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:
*Método de Newmark para sistemas de orden 2
*Método de Runge-Kutta de cuarto orden
Los valores constantes usados para la resolución son los siguientes::
<math>m_{1}=2~[kg]</math>:
<math>m_{2}=1~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>k_{2}=2~[N/m]</math>:
<math>k_{3}=1~[N/m]</math>:

===Método de Newmark===
Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable::
<math>\dot x=u</math>:
<math>\dot y=v</math>
Una vez realizado el cambio de variable el sistema resultante es::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
Con las condiciones iniciales <math>x_{0}=1~[m]</math> e <math>y_{0}=1,5~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all

% Intervalo de tiempo
t0=0;
tN=10;

% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;

%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

% Parametros Newmark
beta=1/4;
gamma=1/2;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;

% Vector solucion S
S=[x(1);y(1);u(1);v(1)];

% Metodo de Newmark

for n=1:N
% Matrices de coeficientes
A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
-beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
-h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];

B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
(1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
-h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

% Sistema
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
u(n+1)=S(3);
v(n+1)=S(4);
end

% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:New2.png|marco|centro]]

===Método de Runge-Kutta===
Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones::
<math>m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y</math>:
<math>m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y</math>:
<math>\dot x=u</math>:
<math>\dot y=v</math>
El código en MATLAB resulta:
{{matlab|codigo=

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);
% Segunda K
k2x=uu(n)+1/2*h;
k2y=vv(n)+1/2*h;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h)+k2/m1*(yy(n)+1/2*k1u*h);
k2v=k2/m2*(xx(n)+1/2*h)-((k2+k3)/m2)*(yy(n)+1/2*k1v*h);

% Tercera K
k3x=uu(n)+1/2*h;
k3y=vv(n)+1/2*h;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h)+k2/m1*(yy(n)+1/2*k2u*h);
k3v=k2/m2*(xx(n)+1/2*h)-((k2+k3)/m2)*(yy(n)+1/2*k2v*h);

% Cuarta K
k4x=uu(n)+h;
k4y=vv(n)+h;
k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3u*h);
k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3v*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off
}}

[[Archivo:Run2.png|marco|centro]]

===Efecto de la velocidad inicial en la trayectoria===
En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos
* Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
* Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.
Estableciendo las condiciones iniciales:
<math> Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} </math>:<math> Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} </math>
Las masas describen estos movimientos: 
'''<big>Caso 1</big>'''
[[Archivo:New3a.png|marco|centro]]
'''<big>Caso 2</big>'''
[[Archivo:New3b.png|marco|centro]]

==Muelle y amortiguador==
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.

[[Archivo:Mra.png|marco|centro]]

Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:

:<math>m_{1}\ddot x+c \dot x +k_{1}x=0</math>

Tomando los parámetros anteriores::

<math>m_{1}=2~[kg]</math>:
<math>k_{1}=4~[N/m]</math>:
<math>\mu=c=1</math>:
La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta

===Método de Newmark===
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable <math>\dot x=y</math> resulta un sistema de ecuaciones diferenciales de orden uno::
<math>\dot x=y</math>:

<math>m_{1}\dot y=-cy-k_{1}x</math>:

Con las condiciones iniciales <math>x_{0}=1~[m]</math> el código en MATLAB del método es:

{{matlab|codigo=

% sistema masa-muelle con amortiguador viscoso
clear all

%tiempo
t0=0;
tN=10;

%datos
k1=4;
m1=2;
c=1;

%condiciones iniciales
x0=1;
y0=0;

%parametros newmark
beta=1/4;
gamma=1/2;

N=1000;
h=(tN-t0)/N;

%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;

%vector solucion
S=[x(1) ; y(1)];

for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end

t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,10])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off

}}

[[Archivo:New4.png|marco|centro]]

===Método de Runge-Kutta===

Archivo:New4.png

2013-03-01T16:13:37Z

Angelcf:

Archivo:Mra.png

2013-03-01T15:52:13Z

Angelcf:

Archivo:Muellemasa.png

2013-03-01T15:40:46Z

Angelcf:

Sistema resorte-masa

2013-03-01T15:14:14Z

Angelcf: /* Efecto de la velocidad inicial */

Sistema resorte-masa

2013-03-01T15:13:33Z

Angelcf: Página creada con «Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas. ==Introducción== El caso considerado está formado por tres muelles y dos masas, unidos a dos ...»

Archivo:New3b.png

2013-03-01T15:08:08Z

Angelcf:

Archivo:New3a.png

2013-03-01T15:07:57Z

Angelcf:

Archivo:Run2.png

2013-03-01T14:47:30Z

Angelcf:

Archivo:New2.png

2013-03-01T14:47:01Z

Angelcf: