MateWiki - Contribuciones del usuario [es]

Sistema resorte-masa (Grupo 8)

2013-03-06T14:16:14Z

Belen.jn: /* Energía mecánica */

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación. 
Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son: :
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2: :
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>:
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=0; z01=1;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta482.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo en escala logarítmica:

[[archivo:Energia6585.jpg]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posición de equilibrio.

 

--[[Usuario:Belen.jn|Belen.jn]] ([[Usuario discusión:Belen.jn|discusión]]) 14:45 4 mar 2013 (CET)

Archivo:Energia6585.jpg

2013-03-06T14:15:38Z

Belen.jn:

Sistema resorte-masa (Grupo 8)

2013-03-06T14:13:07Z

Belen.jn: /* Método de Runge-Kutta */

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación. 
Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son: :
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2: :
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>:
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=0; z01=1;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta482.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo en escala logarítmica:

[[archivo:Energiaaaa.jpg]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posición de equilibrio.

 

--[[Usuario:Belen.jn|Belen.jn]] ([[Usuario discusión:Belen.jn|discusión]]) 14:45 4 mar 2013 (CET)

Archivo:Rungekutta482.jpg

2013-03-06T14:12:15Z

Belen.jn:

Sistema resorte-masa (Grupo 8)

2013-03-06T14:07:34Z

Belen.jn: /* Energía mecánica */

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación. 
Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son: :
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2: :
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>:
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=1; z01=0;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta48.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo en escala logarítmica:

[[archivo:Energiaaaa.jpg]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posición de equilibrio.

 

--[[Usuario:Belen.jn|Belen.jn]] ([[Usuario discusión:Belen.jn|discusión]]) 14:45 4 mar 2013 (CET)

Archivo:Energiaaaa.jpg

2013-03-06T14:07:00Z

Belen.jn:

Sistema resorte-masa (Grupo 8)

2013-03-05T10:32:13Z

Belen.jn: /* Energía mecánica */

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación. 
Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son: :
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2: :
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>:
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=1; z01=0;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta48.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo en escala logarítmica:

[[archivo:Energia8.jpg]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posición de equilibrio.

 

--[[Usuario:Belen.jn|Belen.jn]] ([[Usuario discusión:Belen.jn|discusión]]) 14:45 4 mar 2013 (CET)

Archivo:Energia8.jpg

2013-03-05T10:31:30Z

Belen.jn: Belen.jn subió una nueva versión de «Archivo:Energia8.jpg»

Sistema resorte-masa (Grupo 8)

2013-03-04T13:46:14Z

Belen.jn: /* Energía mecánica */

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación. 
Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son: :
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2: :
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>:
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=1; z01=0;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta48.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo:

[[archivo:Energia8.jpg ]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posición de equilibrio.

 

--[[Usuario:Belen.jn|Belen.jn]] ([[Usuario discusión:Belen.jn|discusión]]) 14:45 4 mar 2013 (CET)

Sistema resorte-masa (Grupo 8)

2013-03-04T13:45:53Z

Belen.jn:

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación. 
Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son: :
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2: :
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>:
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=1; z01=0;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta48.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo:

[[archivo:Energia8.jpg ]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posición de equilibrio.
--[[Usuario:Belen.jn|Belen.jn]] ([[Usuario discusión:Belen.jn|discusión]]) 14:45 4 mar 2013 (CET)

Sistema resorte-masa (Grupo 8)

2013-03-04T13:42:31Z

Belen.jn: /* Energía mecánica */

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación. 
Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son: :
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2: :
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>:
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=1; z01=0;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta48.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo:

[[archivo:Energia8.jpg ]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posición de equilibrio.

Sistema resorte-masa (Grupo 8)

2013-03-04T13:41:16Z

Belen.jn: /* Método de Runge-Kutta */

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación. 
Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son: :
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2: :
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>:
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=1; z01=0;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta48.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo:

[[archivo:Energia8.jpg ]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posción de equilibrio.

Sistema resorte-masa (Grupo 8)

2013-03-04T13:39:35Z

Belen.jn: /* Modelización del sistema de ecuaciones */

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación. 
Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son: :
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2: :
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=1; z01=0;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta48.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo:

[[archivo:Energia8.jpg ]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posción de equilibrio.

Sistema resorte-masa (Grupo 8)

2013-03-04T13:39:17Z

Belen.jn: /* Modelización del sistema de ecuaciones */

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación.
Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son: :
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2: :
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=1; z01=0;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta48.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo:

[[archivo:Energia8.jpg ]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posción de equilibrio.

Sistema resorte-masa (Grupo 8)

2013-03-04T13:38:12Z

Belen.jn: Página creada con «Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos a dos paredes q...»

Muchos modelos de vibraciones se modelizan con sistemas de
muelles y masas. Consideramos aquí un ejemplo simple formado por tres muelles y dos masas unidos
a dos paredes que se deslizan libremente sobre un plano horizontal, como se muestra en la figura:

[[Archivo: Resortemasa.jpg]]

=='''Modelización del sistema de ecuaciones'''==
Modelizamos el sistema de ecuaciones diferenciales para el desplazamiento de ambas masas desde la posición de equilibrio. El sistema quedará en función de las posiciones de las masas, x1 y x2.
Suponemos que deslizan sobre una superficie horizontal y lisa, por lo que despreciamos la fuerza de rozamiento. Por tanto, las únicas fuerzas a tener en cuenta son las fuerzas que ejercen los muelles.
Según la Ley de Hooke, la fuerza ocasionada por un muelle es proporcional a su elongación. Así para nuestro sistema las fuerzas que actúan sobre la masa 1 son:
<math> F1=-k1x1 </math> :
<math> F2=k2(x2-x1) </math>
Sobre la masa 2:
<math> F3=-k2(x2-x1) </math> :
<math> F4=-k3x2 </math>
Donde k1, k2 y k3 son las constantes elásticas de cada muelle.

Aplicando la segunda ley de Newton, <math> F=mx'' </math>, nos queda un sistema de ecuaciones diferenciales de segundo orden:
<math> m1x1''=-k1x1+k2(x2-x1)</math>:
<math> m2x2''=-k2(x2-x1)-k3x2</math>

=='''Resolución del sistema'''==
Suponemos: 
* Las masas m1=2kg, m2=1kg. 
* Las constantes de los muelles k1=4N/m, k2=2N/m, k3=1N/m. 
* Las posiciones en el equilibrio x1=2, X2=4.
El sistema que se obtiene :
<math> x1''=-3x1+x2</math>:
<math> x2''=2x1-3x2</math>
===Supuesto 1===
Para t=0, se desplazan las masas 1 y 1,5 respectivamente y se sueltan repentinamente. Lo resolvemos mediante dos métodos.
====''Método de Newmark''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=1; z10=0; x20=1.5; z20=0;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Newmark2.jpg]]
====''Método Runge-Kutta''====
Para aplicar el Método de Runge-Kutta es necesario descomponer el sistema en cuatro ecuaciones de primer orden:
<math> x1'=z1</math>:
<math>z1'=x1''=-3x1+x2</math>:
<math>x2'=z2</math>:
<math>z2'=x2''=2x1-3x2</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=1; z01=0; x02=1.5; z02=0;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((3*sqrt(2)+4)/4)*cos(sqrt(3-sqrt(2))*t)-((3*sqrt(2)-4)/4)*cos(sqrt(3+sqrt(2))*t)];
Y2=((3*sqrt(2)+4)/(4*sqrt(2)))*cos(sqrt(3-sqrt(2))*t)+((3*sqrt(2)-4)/(4*sqrt(2)))*cos(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}
Obtenemos la siguiente gráfica tiempo-posición:

[[Archivo:Rungekutta281.jpg ‎]]

====''Observaciones''====
Se observa que la gráfica es igual para ambos métodos. La diferencia a simple vista no se aprecia, pero Runge-Kutta es un método de cuarto orden y por tanto más preciso que el de Newmark, que es de segundo orden.
===Supuesto 2===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en el mismo sentido. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;

A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)&Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3a.jpg ]]
====''Método de Runge-Kutta''====
{{matlab| codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3a8.jpg ]]
===Supuesto 3===
En este caso, partimos de las masas en la posición de equilibrio y con velocidad inicial v=1m/s en sentidos opuestos. Aplicamos ambos métodos de la misma forma, sólo cambian las condiciones iniciales.
====''Método de Newmark''====
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

a1=-3; b1=0; c1=1; d1=0;
a2=2; b2=0; c2=-3; d2=0;
A=[0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0 ];
B=[a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2; 0 0 0 0];
C=[0 0 0 0; a1 b1 c1 d1; 0 0 0 0; a2 b2 c2 d2];
D=((h^2)*beta)*B+(h*gamma)*C;
F=eye(4)-D;

x10=0; z10=1; x20=0; z20=-1;
X(:,1)=[x10; z10; x20; z20];
Y=X(:,1);

for n=1:N
Y=(eye(4)+(h*A)+((h^2)*((1/2)-beta))*B+(h*(1-gamma))*C)*Y;
X(:,n+1)=F\Y;
Y=X(:,n+1);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+X(1,:);
e=4+X(3,:);
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark3b.jpg ‎]]
====''Método de Runge-Kutta''====
{{matlab|codigo=
clear all
t0=0;
tN=10;
N=100;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1 0 0;-3 0 1 0;0 0 0 1;2 0 -3 0];

x01=0; z01=1; x02=0; z02=-1;
x1(1)=x01; z1(1)=z01; x2(1)=x02; z2(1)=z02;
X=[x01;z01;x02;z02];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
x2(n+1)=X(3);
z2(n+1)=X(4);
end

%Solución real
Y1=1/2*[((-sqrt(2)+2)/(2*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)-((-sqrt(2)-2)/(2*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t)];
Y2=((-sqrt(2)+2)/(2*sqrt(2)*sqrt(3-sqrt(2))))*sin(sqrt(3-sqrt(2))*t)+((-sqrt(2)-2)/(2*sqrt(2)*sqrt(3+sqrt(2))))*sin(sqrt(3+sqrt(2))*t);

%Representación gráfica
figure(1)
hold on
d=2+x1;
e=4+x2;
plot(d,t,'x')%Solución real
plot(e,t,'x')%Solución real
f=2+Y1;
g=4+Y2;
plot(f,t)%Solución aproximada
plot(g,t)%Solución aproximada
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Rungekutta3b8.jpg ]]
=='''Sistema de una masa'''==
Suponemos ahora que el sistema tiene una masa, es decir desenganchamos la masa 2 de manera que nos queda un muelle y una masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento proporcional a la velocidad de la masa con coeficiente c=1.
Modelizamos este sistema con los mismos parámetros para la masa y el muelle. Nos queda la ecuación:
<math>x''=-(1/2)x'-2x</math>
===''Método de Newmark''===
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

beta=1/4;
gamma=1/2;

A=[1+2*h^2*beta h^2*beta/2; 2*gamma*h 1+gamma*h/2];
B=[1+h^2+2*h^2*beta h-(h^2)/4+(1/2)*h^2*beta; -2*h+2*gamma*h 1-h/2+gamma*h/2];

y0=0; z0=1;
y(1)=y0; z(1)=z0;
Y=[y(1);z(1)];

for n=1:N
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

figure(1)
d=2+y;
plot(d,t);
hold on
}}

Obteniendo la gráfica tiempo-posición:

[[archivo: Newmark4.jpg ]]
===''Método de Runge-Kutta''===
Descomponemos la ecuación de segundo orden en dos ecuaciones de primer orden, quedando el sistema:
<math>x'=z</math>
<math>z'=x''=-2x-(1/2)z</math>
{{matlab|codigo=
clear all

t0=0;
tN=10;
N=1000;
h=(tN-t0)/N;
t=t0:h:tN;

A=[0 1; -2 -1/2];

x01=1; z01=0;
x1(1)=x01; z1(1)=z01;
X=[x01;z01];

for n=1:N
k1=A*X;
k2=A*(X+1/2*h*k1);
k3=A*(X+1/2*h*k2);
k4=A*(X+h*k3);
X=X+h/6*(k1+2*k2+2*k3+k4);
x1(n+1)=X(1);
z1(n+1)=X(2);
end

figure(1)
hold on
d=2+x1;
plot(d,t)
hold off
}}

Obteniendo la gráfica tiempo-posición:

[[archivo:Rungekutta48.jpg ‎]]

===''Energía mecánica''===
Estudiamos ahora la evolución de la energía mecánica en el tiempo.
La fórmula general es:
<math> E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2</math>
Particularizando en este caso:
<math>E=x'^2+2x^2=z^2+2x^2</math>

Utilizamos la solución del Runge-Kutta anterior y obtenemos la gráfica energía-tiempo:

[[archivo:Energia8.jpg ]]

Al incrementar el coeficiente de amortiguamiento "c", la energía decae más rápidamente. Esto se debe a que cuanto mayor es la viscosidad del medio, más rápido se disipa la energía, volviendo antes la masa a la posción de equilibrio.

Archivo:Energia8.jpg

2013-03-04T13:34:04Z

Belen.jn:

Archivo:Rungekutta48.jpg

2013-03-04T13:19:31Z

Belen.jn:

Archivo:Newmark4.jpg

2013-03-04T13:15:23Z

Belen.jn:

Archivo:Rungekutta3b8.jpg

2013-03-04T13:02:40Z

Belen.jn:

Archivo:Newmark3b.jpg

2013-03-04T13:00:05Z

Belen.jn:

Archivo:Rungekutta3a8.jpg

2013-03-04T12:54:44Z

Belen.jn:

Archivo:Newmark3a.jpg

2013-03-04T12:49:33Z

Belen.jn:

Archivo:Rungekutta281.jpg

2013-03-04T12:37:37Z

Belen.jn: grupo 8

grupo 8

Archivo:Rungekutta2.jpg

2013-03-04T12:36:07Z

Belen.jn: Belen.jn subió una nueva versión de «Archivo:Rungekutta2.jpg»: Rungekutta2 del grupo 8

Archivo:Rungekutta2.jpg

2013-03-04T12:34:07Z

Belen.jn: Belen.jn subió una nueva versión de «Archivo:Rungekutta2.jpg»

Archivo:Newmark2.jpg

2013-03-04T12:22:45Z

Belen.jn:

Archivo:Resortemasa.jpg

2013-03-04T11:06:15Z

Belen.jn: