MateWiki - Contribuciones del usuario [es]

Sistemas resorte-masa

2013-03-05T10:57:58Z

Alvarocr: /* Método Runge-Kutta */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M.

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]
===Comparación===
A pesar de la diferente disposición de las gráficas en el método de Newmark y en el de Runge-Kutta, es apreciable con la propia vista que el resultado obtenido por ambos métodos es muy parecido, aun así hemos planteado un código de Matlab de tal forma que se superpongan ambas gráficas. Al realizar esto se puede apreciar que hasta que no ampliamos la imagen hasta obtener escalas centesimales no se parecia la variación, por lo que ambos métodos dan resultados muy similares. Además, con el fin de asegurarnos de realizar una comparación precisa, se ha planteado en el código un cociente entre las posiciones obtenidas según ambos métodos, obteniéndose como resultado siempre la unidad.
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

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

%RUNGE-KUTTA
y=[1;0;1.5;0];

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];

for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

%NEWMARK
beta=1/4;
gamma=1/2;
S=[1;1.5;0;0];

A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]

for n=1:N
S(:,n+1)= A\(B*S(:,n));
end

%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off

figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.

%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.}}
Al activar el código aquí se planteado se obtienen las siguientes gráficas:
[[Archivo:Comparacion.jpg|marco|centro|Posiciones de ambas masas obtenidas por los dos métodos superpuestas]]
[[Archivo:Comparacionzoom.jpg|marco|centro|Al hacer zoom en la gráfica se acaba apreciando una separación entre ellas]]
[[Archivo:Cociente.jpg|marco|centro|Cociente entre las posiciones de la masa 1 según ambos métodos]]

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off
%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

A partir de ambos métodos se obtienen semejantes gráficas de la velocidad y posición que es la siguiente:
[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
%lo que vamos a hacer es crear un bucle que me va a coger 10 valores de mu,
%desde 1 hasta 5, aumentando de 0.5 en 0.5
g=5;
for i=1:0.5:g
M=[0,1;-2,-i/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
subplot(g*2,3,1+(((i*2)-1)*3)) %el desplazamiento lo colocamos en la primera columna
plot(t,y(1,:),'red');
subplot(g*2,3,2+(((i*2)-1)*3))%la velocidad en la segunda columna
plot(t,y(2,:),'green');
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v));
subplot(g*2,3,3+(((i*2)-1)*3))%la energia en la tercera
plot(t,E)
end

}}

Se muestran a continuación los cambios que produce la μ en la posición, la velocidad y la energía:
[[Archivo:Comparacionmu.jpg|marco|centro|Desplazamiento, velocidad y energía aumentando μ de 1 a 5 en intervalos de 0,5]]

Discusión:Página principal

2013-03-05T10:39:35Z

Alvarocr: /* Comentarios */

= Comentarios =
Somos el grupo 18 y nuestro trabajo Sistemas resorte-masa no aparece en la página principal. Gracias.
Añade un comentario debajo si tu artículo no aparece en la portada. Indica el nombre de tu artículo, y no te olvides de firmar (mira el icono con la pluma/lápiz al editar esta página). --[[Usuario:Herraiz|Herraiz]] ([[Usuario discusión:Herraiz|discusión]]) 23:09 4 mar 2013 (CET)

Sistemas resorte-masa

2013-03-05T10:31:08Z

Alvarocr: /* Influencia de μ */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]
===Comparación===
A pesar de la diferente disposición de las gráficas en el método de Newmark y en el de Runge-Kutta, es apreciable con la propia vista que el resultado obtenido por ambos métodos es muy parecido, aun así hemos planteado un código de Matlab de tal forma que se superpongan ambas gráficas. Al realizar esto se puede apreciar que hasta que no ampliamos la imagen hasta obtener escalas centesimales no se parecia la variación, por lo que ambos métodos dan resultados muy similares. Además, con el fin de asegurarnos de realizar una comparación precisa, se ha planteado en el código un cociente entre las posiciones obtenidas según ambos métodos, obteniéndose como resultado siempre la unidad.
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

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

%RUNGE-KUTTA
y=[1;0;1.5;0];

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];

for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

%NEWMARK
beta=1/4;
gamma=1/2;
S=[1;1.5;0;0];

A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]

for n=1:N
S(:,n+1)= A\(B*S(:,n));
end

%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off

figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.

%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.}}
Al activar el código aquí se planteado se obtienen las siguientes gráficas:
[[Archivo:Comparacion.jpg|marco|centro|Posiciones de ambas masas obtenidas por los dos métodos superpuestas]]
[[Archivo:Comparacionzoom.jpg|marco|centro|Al hacer zoom en la gráfica se acaba apreciando una separación entre ellas]]
[[Archivo:Cociente.jpg|marco|centro|Cociente entre las posiciones de la masa 1 según ambos métodos]]

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off
%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

A partir de ambos métodos se obtienen semejantes gráficas de la velocidad y posición que es la siguiente:
[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
%lo que vamos a hacer es crear un bucle que me va a coger 10 valores de mu,
%desde 1 hasta 5, aumentando de 0.5 en 0.5
g=5;
for i=1:0.5:g
M=[0,1;-2,-i/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
subplot(g*2,3,1+(((i*2)-1)*3)) %el desplazamiento lo colocamos en la primera columna
plot(t,y(1,:),'red');
subplot(g*2,3,2+(((i*2)-1)*3))%la velocidad en la segunda columna
plot(t,y(2,:),'green');
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v));
subplot(g*2,3,3+(((i*2)-1)*3))%la energia en la tercera
plot(t,E)
end

}}

Se muestran a continuación los cambios que produce la μ en la posición, la velocidad y la energía:
[[Archivo:Comparacionmu.jpg|marco|centro|Desplazamiento, velocidad y energía aumentando μ de 1 a 5 en intervalos de 0,5]]

Archivo:Comparacionmu.jpg

2013-03-05T10:28:52Z

Alvarocr: Comparación de los desplazamientos, velocidades y energías

Comparación de los desplazamientos, velocidades y energías

Sistemas resorte-masa

2013-03-05T10:27:12Z

Alvarocr: /* Método de Newmark */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]
===Comparación===
A pesar de la diferente disposición de las gráficas en el método de Newmark y en el de Runge-Kutta, es apreciable con la propia vista que el resultado obtenido por ambos métodos es muy parecido, aun así hemos planteado un código de Matlab de tal forma que se superpongan ambas gráficas. Al realizar esto se puede apreciar que hasta que no ampliamos la imagen hasta obtener escalas centesimales no se parecia la variación, por lo que ambos métodos dan resultados muy similares. Además, con el fin de asegurarnos de realizar una comparación precisa, se ha planteado en el código un cociente entre las posiciones obtenidas según ambos métodos, obteniéndose como resultado siempre la unidad.
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

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

%RUNGE-KUTTA
y=[1;0;1.5;0];

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];

for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

%NEWMARK
beta=1/4;
gamma=1/2;
S=[1;1.5;0;0];

A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]

for n=1:N
S(:,n+1)= A\(B*S(:,n));
end

%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off

figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.

%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.}}
Al activar el código aquí se planteado se obtienen las siguientes gráficas:
[[Archivo:Comparacion.jpg|marco|centro|Posiciones de ambas masas obtenidas por los dos métodos superpuestas]]
[[Archivo:Comparacionzoom.jpg|marco|centro|Al hacer zoom en la gráfica se acaba apreciando una separación entre ellas]]
[[Archivo:Cociente.jpg|marco|centro|Cociente entre las posiciones de la masa 1 según ambos métodos]]

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off
%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

A partir de ambos métodos se obtienen semejantes gráficas de la velocidad y posición que es la siguiente:
[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
%lo que vamos a hacer es crear un bucle que me va a coger 10 valores de mu,
%desde 1 hasta 5, aumentando de 0.5 en 0.5
g=5;
for i=1:0.5:g
M=[0,1;-2,-i/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
subplot(g*2,3,1+(((i*2)-1)*3)) %el desplazamiento lo colocamos en la primera columna
plot(t,y(1,:),'red');
subplot(g*2,3,2+(((i*2)-1)*3))%la velocidad en la segunda columna
plot(t,y(2,:),'green');
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v));
subplot(g*2,3,3+(((i*2)-1)*3))%la energia en la tercera
plot(t,E)
end

}}

Sistemas resorte-masa

2013-03-05T10:26:02Z

Alvarocr: /* Comparación */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]
===Comparación===
A pesar de la diferente disposición de las gráficas en el método de Newmark y en el de Runge-Kutta, es apreciable con la propia vista que el resultado obtenido por ambos métodos es muy parecido, aun así hemos planteado un código de Matlab de tal forma que se superpongan ambas gráficas. Al realizar esto se puede apreciar que hasta que no ampliamos la imagen hasta obtener escalas centesimales no se parecia la variación, por lo que ambos métodos dan resultados muy similares. Además, con el fin de asegurarnos de realizar una comparación precisa, se ha planteado en el código un cociente entre las posiciones obtenidas según ambos métodos, obteniéndose como resultado siempre la unidad.
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

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

%RUNGE-KUTTA
y=[1;0;1.5;0];

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];

for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

%NEWMARK
beta=1/4;
gamma=1/2;
S=[1;1.5;0;0];

A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]

for n=1:N
S(:,n+1)= A\(B*S(:,n));
end

%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off

figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.

%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.}}
Al activar el código aquí se planteado se obtienen las siguientes gráficas:
[[Archivo:Comparacion.jpg|marco|centro|Posiciones de ambas masas obtenidas por los dos métodos superpuestas]]
[[Archivo:Comparacionzoom.jpg|marco|centro|Al hacer zoom en la gráfica se acaba apreciando una separación entre ellas]]
[[Archivo:Cociente.jpg|marco|centro|Cociente entre las posiciones de la masa 1 según ambos métodos]]

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off
%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
%lo que vamos a hacer es crear un bucle que me va a coger 10 valores de mu,
%desde 1 hasta 5, aumentando de 0.5 en 0.5
g=5;
for i=1:0.5:g
M=[0,1;-2,-i/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
subplot(g*2,3,1+(((i*2)-1)*3)) %el desplazamiento lo colocamos en la primera columna
plot(t,y(1,:),'red');
subplot(g*2,3,2+(((i*2)-1)*3))%la velocidad en la segunda columna
plot(t,y(2,:),'green');
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v));
subplot(g*2,3,3+(((i*2)-1)*3))%la energia en la tercera
plot(t,E)
end

}}

Sistemas resorte-masa

2013-03-05T10:25:33Z

Alvarocr: /* Influencia de μ */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]
===Comparación===
A pesar de la diferente disposición de las gráficas en el método de Newmark y en el de Runge-Kutta, es apreciable con la propia vista que el resultado obtenido por ambos métodos es muy parecido, aun así hemos planteado un código de Matlab de tal forma que se superpongan ambas gráficas. Al realizar esto se puede apreciar que hasta que no ampliamos la imagen hasta obtener escalas centesimales no se parecia la variación, por lo que ambos métodos dan resultados muy similares. Además, con el fin de asegurarnos de realizar una comparación precisa, se ha planteado en el código un cociente entre las posiciones obtenidas según ambos métodos, obteniéndose como resultado siempre la unidad.
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

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

%RUNGE-KUTTA
y=[1;0;1.5;0];

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];

for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

%NEWMARK
beta=1/4;
gamma=1/2;
S=[1;1.5;0;0];

A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]

for n=1:N
S(:,n+1)= A\(B*S(:,n));
end

%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off

figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.

%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.}}

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off
%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
%lo que vamos a hacer es crear un bucle que me va a coger 10 valores de mu,
%desde 1 hasta 5, aumentando de 0.5 en 0.5
g=5;
for i=1:0.5:g
M=[0,1;-2,-i/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
subplot(g*2,3,1+(((i*2)-1)*3)) %el desplazamiento lo colocamos en la primera columna
plot(t,y(1,:),'red');
subplot(g*2,3,2+(((i*2)-1)*3))%la velocidad en la segunda columna
plot(t,y(2,:),'green');
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v));
subplot(g*2,3,3+(((i*2)-1)*3))%la energia en la tercera
plot(t,E)
end

}}

Archivo:Cociente.jpg

2013-03-05T10:25:27Z

Alvarocr: Cociente entre las posiciones de x

Cociente entre las posiciones de x

Sistemas resorte-masa

2013-03-05T10:21:05Z

Alvarocr: /* Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]
===Comparación===
A pesar de la diferente disposición de las gráficas en el método de Newmark y en el de Runge-Kutta, es apreciable con la propia vista que el resultado obtenido por ambos métodos es muy parecido, aun así hemos planteado un código de Matlab de tal forma que se superpongan ambas gráficas. Al realizar esto se puede apreciar que hasta que no ampliamos la imagen hasta obtener escalas centesimales no se parecia la variación, por lo que ambos métodos dan resultados muy similares. Además, con el fin de asegurarnos de realizar una comparación precisa, se ha planteado en el código un cociente entre las posiciones obtenidas según ambos métodos, obteniéndose como resultado siempre la unidad.
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

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

%RUNGE-KUTTA
y=[1;0;1.5;0];

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];

for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

%NEWMARK
beta=1/4;
gamma=1/2;
S=[1;1.5;0;0];

A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]

for n=1:N
S(:,n+1)= A\(B*S(:,n));
end

%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off

figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.

%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.}}

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off
%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;
t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
%RUNGE-KUTTA
y=[1;0;1.5;0];
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

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

S=[1;1.5;0;0];
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
S(:,n+1)= A\(B*S(:,n));
end

%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off
figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.
%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.

}}

Archivo:Comparacionzoom.jpg

2013-03-05T10:18:50Z

Alvarocr: Zoom de las comparaciones

Zoom de las comparaciones

Sistemas resorte-masa

2013-03-05T10:18:10Z

Alvarocr: /* Influencia de μ */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]
===Comparación===
A pesar de la diferente disposición de las gráficas en el método de Newmark y en el de Runge-Kutta, es apreciable con la propia vista que el resultado obtenido por ambos métodos es muy parecido, aun así hemos planteado un código de Matlab de tal forma que se superpongan ambas gráficas. Al realizar esto se puede apreciar que hasta que no ampliamos la imagen hasta obtener escalas centesimales no se parecia la variación, por lo que ambos métodos dan resultados muy similares. Además, con el fin de asegurarnos de realizar una comparación precisa, se ha planteado en el código un cociente entre las posiciones obtenidas según ambos métodos, obteniéndose como resultado siempre la unidad.
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

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

%RUNGE-KUTTA
y=[1;0;1.5;0];

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];

for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

%NEWMARK
beta=1/4;
gamma=1/2;
S=[1;1.5;0;0];

A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]

for n=1:N
S(:,n+1)= A\(B*S(:,n));
end

%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off

figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.

%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.}}

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;
t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
%RUNGE-KUTTA
y=[1;0;1.5;0];
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

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

S=[1;1.5;0;0];
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
S(:,n+1)= A\(B*S(:,n));
end

%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off
figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.
%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.

}}

Archivo:Comparacion.jpg

2013-03-05T10:18:08Z

Alvarocr: Superposición de las gráficas siguiendo Newmark y Runge-Kutta

Superposición de las gráficas siguiendo Newmark y Runge-Kutta

Sistemas resorte-masa

2013-03-05T10:17:35Z

Alvarocr: /* Aproximación de la posición mediante métodos numéricos */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]
===Comparación===
A pesar de la diferente disposición de las gráficas en el método de Newmark y en el de Runge-Kutta, es apreciable con la propia vista que el resultado obtenido por ambos métodos es muy parecido, aun así hemos planteado un código de Matlab de tal forma que se superpongan ambas gráficas. Al realizar esto se puede apreciar que hasta que no ampliamos la imagen hasta obtener escalas centesimales no se parecia la variación, por lo que ambos métodos dan resultados muy similares. Además, con el fin de asegurarnos de realizar una comparación precisa, se ha planteado en el código un cociente entre las posiciones obtenidas según ambos métodos, obteniéndose como resultado siempre la unidad.
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

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

%RUNGE-KUTTA
y=[1;0;1.5;0];

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];

for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

%NEWMARK
beta=1/4;
gamma=1/2;
S=[1;1.5;0;0];

A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]

for n=1:N
S(:,n+1)= A\(B*S(:,n));
end

%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off

figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.

%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.}}

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.
{{matlab|codigo=%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;t0=0; tN=10; N=400; h=(tN-t0)/N;t=t0:h:tN;%RUNGE-KUTTAy=[1;0;1.5;0];M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; for n=1:N k1=M*y(:,n); k2=M*(y(:,n)+1/2*k1*h); k3=M*(y(:,n)+1/2*k2*h); k4=M*(y(:,n)+k3*h); y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);end%NEWMARKbeta=1/4;gamma=1/2;S=[1;1.5;0;0]; A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1] for n=1:N S(:,n+1)= A\(B*S(:,n));end%COMPARACION DE DESPLAZAMIENTOSfigure(1)hold ony(3,:)=y(3,:)+4S(2,:)=S(2,:)+4plot(y(1,:),t,'black');plot(y(3,:),t,'black');plot(S(1,:),t,'red');plot(S(2,:),t,'red');hold offfigure(2)z=y(1,:)/S(1,:);plot(t,z);%como podemos ver el la gráfica que nos divide las x de rk con las de%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere%decir que la diferencia entre ambos métodos en prácticamente nula.%NOTA:se podría hacer lo mismo comparando las velocidades y saldría lo mismo.}}

Sistemas resorte-masa

2013-03-05T10:16:09Z

Alvarocr: /* Influencia de μ */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.
{{matlab|codigo=%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;t0=0; tN=10; N=400; h=(tN-t0)/N;t=t0:h:tN;%RUNGE-KUTTAy=[1;0;1.5;0];M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; for n=1:N k1=M*y(:,n); k2=M*(y(:,n)+1/2*k1*h); k3=M*(y(:,n)+1/2*k2*h); k4=M*(y(:,n)+k3*h); y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);end%NEWMARKbeta=1/4;gamma=1/2;S=[1;1.5;0;0]; A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1] for n=1:N S(:,n+1)= A\(B*S(:,n));end%COMPARACION DE DESPLAZAMIENTOSfigure(1)hold ony(3,:)=y(3,:)+4S(2,:)=S(2,:)+4plot(y(1,:),t,'black');plot(y(3,:),t,'black');plot(S(1,:),t,'red');plot(S(2,:),t,'red');hold offfigure(2)z=y(1,:)/S(1,:);plot(t,z);%como podemos ver el la gráfica que nos divide las x de rk con las de%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere%decir que la diferencia entre ambos métodos en prácticamente nula.%NOTA:se podría hacer lo mismo comparando las velocidades y saldría lo mismo.}}

Sistemas resorte-masa

2013-03-05T10:12:23Z

Alvarocr: /* Influencia de μ */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.

Sistemas resorte-masa

2013-03-05T10:09:22Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al decrecer mayor velocidad en el mismo tiempo pues pierde energía cinética.

Sistemas resorte-masa

2013-03-05T10:08:52Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-C/2\dot x</math>

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
<math>\dot x=u</math>
Con lo que el sistema nos queda:
<math>\dot x=u\\\dot u=-2x-(C/2)u</math>

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al decrecer mayor velocidad en el mismo tiempo pues pierde energía cinética.

Sistemas resorte-masa

2013-03-05T10:04:56Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.
== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:2m3mu.jpg|marco|centro|Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas
[[Archivo:Nemark2posicion.jpg|marco|centro|Masa 1 en color azul, masa 2 en negro]]

[[Archivo:Newmark2velo.jpg|marco|centro|La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde]]

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off}}
Las gráfica que define el movimiento en este caso es la siguiente:
[[Archivo:Rk3a.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:
[[Archivo:Rk3b.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

[[Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png|marco|centro|Sistema con amortiguamiento]]

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
<math>-C\dot x-Kx=m\ddot x\\\ddot x=-2x-C/2\dot x</math>

====Método Runge-Kutta====
{{matlab|codigo=

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)
}}

====Método de Newmark====
{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
B=[1-(h^2)/4, h-(h^2)/8; ...
-h/2, 1-h/4];
Y=A\(B*Y);
y(n+1)=Y(1);
z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)
}}

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
(cv)
Con lo que el sistema nos queda:
(sistema)
Aplicando los métodos, obtenemos:
(métodos)
[[Archivo:Rk4pos.jpg|marco|centro|Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo]]

===Energía===
Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.
[[Archivo:Rk4energia.jpg|marco|centro|Energía mecánica de la partícula a lo largo del tiempo]]

===Influencia de μ===
En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al decrecer mayor velocidad en el mismo tiempo pues pierde energía cinética.

Sistemas resorte-masa

2013-03-05T10:03:20Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T10:00:34Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T09:57:17Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T09:56:33Z

Alvarocr:

Sistemas resorte-masa

2013-03-05T09:56:29Z

Alvarocr: /* Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas */

Sistemas resorte-masa

2013-03-05T09:48:31Z

Alvarocr: /* Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas */

Sistemas resorte-masa

2013-03-05T09:46:46Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T09:44:16Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T09:43:40Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T09:43:10Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T09:42:35Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T09:42:16Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T09:41:47Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T09:40:35Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Sistemas resorte-masa

2013-03-05T09:39:51Z

Alvarocr:

Sistemas resorte-masa

2013-03-05T09:37:32Z

Alvarocr: /* Sistema de una masa en medio viscoso */

Archivo:Captura de pantalla 2013-03-05 a la(s) 10.31.43.png

2013-03-05T09:36:32Z

Alvarocr: Sistema con amortiguamiento

Sistema con amortiguamiento

Sistemas resorte-masa

2013-03-05T09:33:50Z

Alvarocr: /* Sistemas resorte-masa */

Archivo:2m3mu.jpg

2013-03-05T09:33:05Z

Alvarocr: Sistema de dos masas y tres muelles

Sistema de dos masas y tres muelles

Sistemas resorte-masa

2013-03-05T09:31:36Z

Alvarocr: /* Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas */

Sistemas resorte-masa

2013-03-05T09:30:45Z

Alvarocr: /* Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas */

Sistemas resorte-masa

2013-03-05T09:29:21Z

Alvarocr: /* Método de Newmark */

Sistemas resorte-masa

2013-03-05T09:28:26Z

Alvarocr: /* Método de Newmark */

Sistemas resorte-masa

2013-03-05T09:27:57Z

Alvarocr: /* Método de Newmark */

Sistemas resorte-masa

2013-03-05T09:24:42Z

Alvarocr: /* Método de Newmark */

Sistemas resorte-masa

2013-03-05T09:24:00Z

Alvarocr: /* Método de Newmark */

Sistemas resorte-masa

2013-03-05T09:21:34Z

Alvarocr:

Sistemas resorte-masa

2013-03-05T09:20:04Z

Alvarocr: /* Sistemas resorte-masa (grupo 18) */

== Sistemas resorte-masa==

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

[[Archivo:Sistema doble|miniaturadeimagen|centro]]

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.</math>

== Aproximación de la posición mediante métodos numéricos ==
Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

<math>\dot x=u\\\dot y=v</math>

Y así obtenemos el sistema:
<math>\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\\\dot x=u\\\dot y=v\end{matrix}\right.</math>

En estas ecuaciones sustituimos los datos que plantean en el problema:

<math>m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m</math>

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación:
<math>x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0</math>

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

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M:

Sabiendo que el método de Runge-Kutta consiste en:
<math>y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)</math>
En nuestro, el programa queda de la siguiente manera:
{{matlab|codigo=

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0

M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema

%Aplicación del método Runge-Kutta
for n=1:N
k1=M*y(:,n);
k2=M*(y(:,n)+1/2*k1*h);
k3=M*(y(:,n)+1/2*k2*h);
k4=M*(y(:,n)+k3*h);
y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off

%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off }}

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:
[[Archivo:Rk2.jpg|marco|centro|Las posiciones se muestran en negro y las velocidades en verde]]

===Método de Newmark===

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

<math>x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)</math>

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

{{matlab|codigo=

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red'); %w=y'
hold off }}

==Posiciones de las masas en función del tiempo para los siguientes casos==
===Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas===
Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

<math>u(0)= 1m/s\\v(0)= 1m/s</math>

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo.
Las gráfica que define el movimiento en este caso es la siguiente:

===Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas===
Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

<math>u(0)= 1m/s\\v(0)= -1m/s</math>

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:

==Sistema de una masa en medio viscoso==

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1.
Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:
(suma de fuerzas)
Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable:
(cv)
Con lo que el sistema nos queda:
(sistema)
Aplicando los métodos, obtenemos:
(métodos)
(graficas)

Archivo:Newmark2velo.jpg

2013-03-05T09:19:33Z

Alvarocr: Velocidades de las masas 1 y 2 siguiendo el método de Newmark

Velocidades de las masas 1 y 2 siguiendo el método de Newmark

Archivo:Nemark2posicion.jpg

2013-03-05T09:19:07Z

Alvarocr: Posiciones de las masas 1 y 2 según el método Newmark

Posiciones de las masas 1 y 2 según el método Newmark

Archivo:Rk4energia.jpg

2013-03-05T09:18:33Z

Alvarocr: Energía mecánica de la masa en el apartado 4

Energía mecánica de la masa en el apartado 4

Archivo:Rk4pos.jpg

2013-03-05T09:18:09Z

Alvarocr: Posición y velocidad de la masa en el apartado 4

Posición y velocidad de la masa en el apartado 4