Modelo de vibraciones con muelles y masas (Grupo 3B)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelo de vibraciones con muelles y masas (Grupo 3B) |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2016-17 |
| Autores | Alejandro García García - 2006; Cindy Devia Preciado - 2134; Juan Carlos López Segovia - 2183; Lays Oscarina de Souza Soares - 2333; Rafael Ventura Márquez de Prado Arrarás - 2019; Sofía de Miguel González - 2132 |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
En este trabajo vamos a estudiar el comportamiento de tres masas unidas entre sí y con dos paredes en los extremos mediante muelles, que deslizan libremente sobre un plano horizontal.
Contenido
1 Modelización del sistema
Llamaremos [math] p_i0 [/math] a la posición de equilibrio para cada masa y [math] x_i [/math] al desplazamiento de cada muelle, el cual consideraremos positivo hacia la derecha. Entonces [math] P_i [/math] que es la posición de cada masa se define así::
[math] P_1(t)=p_10+x_1(t) [/math]: [math] P_2(t)=p_20+x_2(t) [/math]: [math] P_3(t)=p_30+x_3(t) [/math]:
Llamando [math] k_i [/math] a las constantes de recuperación de cada muelle, considerando nulo el rozamiento y ningún amortiguamiento se puede ver que las fuerzas que actúan sobre cada masa son: Masa 1: [math] F_1=-k_1x_1 [/math]:
[math] F_2=k_2(x_2-x_1) [/math]
Masa 2: [math] F_3=-k_2(x_2-x_1) [/math]:
[math] F_4=k_3(x_3-x_2) [/math]
Masa 3: [math] F_5=-k_3(x_3-x_2) [/math]:
[math] F_6=-k_4x_3 [/math]
Aplicando la Ley de Hooke y la Segunda Ley de Newton, se obtiene el siguiente sistema de ecuaciones:: [math]\begin{cases}m_1x_1''=-k_1x_1+k_2(x_2-x_1)\\m_2x_2''=-k_2(x_2-x_1)+k_3(x_3-x_2)\\m_3x_3''=-k_3(x_3-x_2)-k_4x_3\end{cases}[/math]
Ahora reducimos para un sistema de orden 1:: [math]\begin{cases}x_1'=v_1\\x_2'=v_2\\x_3'=v_3\\v_1'=\frac{1}{m_1}(-k_1x_1+k_2(x_2-x_1))\\v_2'=\frac{1}{m_2}(-k_2(x_2-x_1)+k_3(x_3-x_2))\\v_3'=\frac{1}{m_3}(-k_3(x_3-x_2)-k_4x_3)\end{cases}[/math]
Para resolverlo necesitaremos seis condiciones iniciales que serán las posiciones y velocidades iniciales de cada masa.
2 Resolución del sistema
En primer lugar resolveremos el sistema tomando los desplazamientos iniciales de las masas como 0.4, 0.9 y 0.6 metros hacia la derecha de la posición de equilibrio y con velocidad inicial nula en las tres masas. Para eso vamos a considerar que [math] k_1 = 3.5 N/m [/math], [math] k_2 = 2.1 N/m [/math], [math] k_3 = 2.5 N/m [/math], [math] k_4 = 3.1 N/m [/math], [math] m_1 = 2.3 kg [/math], [math] m_2 = 1.1 kg [/math], [math] m_3 = 3.3 kg [/math], considerando también que la distancia entre las paredes es de 13 metros y las posiciones de equilibrio de las tres masas son 3, 6 y 9 metros.
De esta forma el problema queda así:: [math]\begin{cases}x_1'=v_1\\x_2'=v_2\\x_3'=v_3\\v_1'=\frac{1}{2.3}(-3.5x_1+2.1(x_2-x_1))\\v_2'=\frac{1}{1.1}(-2.1(x_2-x_1)+2.5(x_3-x_2))\\v_3'=\frac{1}{3.3}(-2.5(x_3-x_2)-3.1x_3)\end{cases}[/math]
Resolveremos utilizando los siguientes métodos:
- Método de Runge-Kutta
- Método de Euler modificado
Para eso utilizaremos [math] h = 0.025 m [/math] y [math] h = 0.1 m [/math]
2.1 Runge-Kutta
El primer código de Matlab (Con paso [math] h = 0.025 [/math]) es este:
% 1º - Datos del programa:
t0=0; %Tiempo Inicial
tf=20; %Tiempo Final
y0=[0.4; 0.9;0.6;0;0;0]; %El orden con el que vamos a obtener las soluciones es x1, x2, x3, x1', x2', x3'
% 2º - Discretizacion:
h=0.025; %Tamaño del paso
N=(tf-t0)/h; %Número de intervalos
t=t0:h:tf; %Vector tiempo
% 3º - Metodo numerico (Metodo de Runge-Kutta):
%Creamos el vector de soluciones (vacío)
y=zeros(6,N+1);
%La solución en el inicio es dato
y(:,1)=y0;
%Convertimos el sistema en una matriz por el vector 'y', obteniéndose la función F
F=@(t,y) [y(4);y(5);y(6);(1/2.3)*(-3.5*y(1)+2.1*(y(2)-y(1)));(1/1.1)*(-2.1*(y(2)-y(1))+2.5*(y(3)-y(2)));(1/3.3)*(-2.5*(y(3)-y(2))-3.1*y(3))];
%Aplicamos el método Runge-Kutta
for n=1:N
k1=F(t(n),y(:,n));
k2=F(t(n)+0.5*h,y(:,n)+0.5*k1*h);
k3=F(t(n)+0.5*h,y(:,n)+0.5*k2*h);
k4=F(t(n)+h,y(n)+k3*h);
y(:,n+1)=y(:,n)+(h/6)*(k1+2*k2+2*k3+k4);
end
%Ahora para calcular los desplazamientos entorno a las posiciones de
%equilibrio de cada masa, sumamos la propia posición de equilibrio a cada
%una
y(1,:)=y(1,:)+ones(1,N+1)*3;
y(2,:)=y(2,:)+ones(1,N+1)*6;
y(3,:)=y(3,:)+ones(1,N+1)*9;
% 4º - Post Proceso:
%Ahora dibujamos el gráfico y añadimos los nombres de los ejes y un título
plot(t,y(1:3,:))
ylabel('Posicion de la masa')
xlabel('Tiempo')
title('Runge-Kutta')
legend('x1(t)','x2(t)','x3(t)');
Y se obtiene la siguiente gráfica:
En esta gráfica podemos observar que las tres masas no van a oscilar de la misma manera; la uno y la dos lo harán de manera irregular mientras que la tres se asemeja más a un movimiento armónico. Además se ve que las masas no tienden a detenerse para un tiempo infinito sino que oscilan en torno a la posición de equilibrio.
El segundo código de Matlab (Con paso [math] h = 0.1 [/math]) es este:
% 1º - Datos del programa:
t0=0; %Tiempo Inicial
tf=20; %Tiempo Final
y0=[0.4; 0.9;0.6;0;0;0]; %El orden con el que vamos a obtener las soluciones es x1, x2, x3, x1', x2', x3'
% 2º - Discretizacion:
h=0.1; %Tamaño del paso
N=(tf-t0)/h; %Número de intervalos
t=t0:h:tf; %Vector tiempo
% 3º - Metodo numerico (Metodo de Runge-Kutta):
%Creamos el vector de soluciones (vacío)
y=zeros(6,N+1);
%La solución en el inicio es dato
y(:,1)=y0;
%Convertimos el sistema en una matriz por el vector 'y', obteniéndose la función F
F=@(t,y) [y(4);y(5);y(6);(1/2.3)*(-3.5*y(1)+2.1*(y(2)-y(1)));(1/1.1)*(-2.1*(y(2)-y(1))+2.5*(y(3)-y(2)));(1/3.3)*(-2.5*(y(3)-y(2))-3.1*y(3))];
%Aplicamos el método Runge-Kutta
for n=1:N
k1=F(t(n),y(:,n));
k2=F(t(n)+0.5*h,y(:,n)+0.5*k1*h);
k3=F(t(n)+0.5*h,y(:,n)+0.5*k2*h);
k4=F(t(n)+h,y(n)+k3*h);
y(:,n+1)=y(:,n)+(h/6)*(k1+2*k2+2*k3+k4);
end
%Ahora para calcular los desplazamientos entorno a las posiciones de
%equilibrio de cada masa, sumamos la propia posición de equilibrio a cada
%una
y(1,:)=y(1,:)+ones(1,N+1)*3;
y(2,:)=y(2,:)+ones(1,N+1)*6;
y(3,:)=y(3,:)+ones(1,N+1)*9;
% 4º - Post Proceso:
%Ahora dibujamos el gráfico y añadimos los nombres de los ejes y un título
plot(t,y(1:3,:))
ylabel('Posicion de la masa')
xlabel('Tiempo')
title('Runge-Kutta 2')
legend('x1(t)','x2(t)','x3(t)');
Utilizando ese código, hemos obtenido la siguiente gráfica:
Podemos observar que las tres masas siguen sin oscilar de la misma manera, y a pesar de que las masas uno y dos siguieren oscilando de una manera irregular, las tres oscilaciones se han "suavizado" muy minimamente, en otras palabras, están más parecidas con un movimiento armónico que en el caso anterior, pero de una manera muy leve, casi imperceptible a simple vista. Además, seguimos viendo que las masas no tienden a detenerse para un tiempo infinito sino que oscilan en torno a la posición de equilibrio.
2.2 Euler Modificado
El primer código de Matlab (Con paso [math] h = 0.025 [/math]) es este:
% 1º - Datos del programa:
t0=0; %Tiempo Inicial
tf=20; %Tiempo Final
y0=[0.4; 0.9;0.6;0;0;0]; %El orden con el que vamos a obtener las soluciones es x1, x2, x3, x1', x2', x3'
% 2º - Discretizacion:
h=0.025; %Tamaño del paso
N=(tf-t0)/h; %Número de intervalos
t=t0:h:tf; %Vector tiempo
% 3º - Metodo numerico (Metodo de Euler Modificado):
%Creamos el vector de soluciones (vacío)
y=zeros(6,N+1);
%La solución en el inicio es dato
y(:,1)=y0;
%Convertimos el sistema en una matriz por el vector 'y', obteniéndose la función F
F=@(t,y) [y(4);y(5);y(6);(1/2.3)*(-3.5*y(1)+2.1*(y(2)-y(1)));(1/1.1)*(-2.1*(y(2)-y(1))+2.5*(y(3)-y(2)));(1/3.3)*(-2.5*(y(3)-y(2))-3.1*y(3))];
%Aplicamos el método de Euler Modificado
for n=1:N
k1=F(t(n),y(:,n));
y(:,n+1)=y(:,n)+h*F(t(n)+0.5*h,y(:,n)+0.5*k1*h);
end
%Ahora para calcular los desplazamientos entorno a las posiciones de
%equilibrio de cada masa, sumamos la propia posición de equilibrio a cada
%una
y(1,:)=y(1,:)+ones(1,N+1)*3;
y(2,:)=y(2,:)+ones(1,N+1)*6;
y(3,:)=y(3,:)+ones(1,N+1)*9;
% 4º - Post Proceso:
%Ahora dibujamos el gráfico y añadimos los nombres de los ejes y un título
plot(t,y(1:3,:))
ylabel('Posicion de la masa')
xlabel('Tiempo')
title('Euler Modificado')
legend('x1(t)','x2(t)','x3(t)');
La gráfica generada es la siguiente:
El segundo código de Matlab (Con paso [math] h = 0.1 [/math]) es este:
% 1º - Datos del programa:
t0=0; %Tiempo Inicial
tf=20; %Tiempo Final
y0=[0.4; 0.9;0.6;0;0;0]; %El orden con el que vamos a obtener las soluciones es x1, x2, x3, x1', x2', x3'
% 2º - Discretizacion:
h=0.1; %Tamaño del paso
N=(tf-t0)/h; %Número de intervalos
t=t0:h:tf; %Vector tiempo
% 3º - Metodo numerico (Metodo de Euler Modificado):
%Creamos el vector de soluciones (vacío)
y=zeros(6,N+1);
%La solución en el inicio es dato
y(:,1)=y0;
%Convertimos el sistema en una matriz por el vector 'y', obteniéndose la función F
F=@(t,y) [y(4);y(5);y(6);(1/2.3)*(-3.5*y(1)+2.1*(y(2)-y(1)));(1/1.1)*(-2.1*(y(2)-y(1))+2.5*(y(3)-y(2)));(1/3.3)*(-2.5*(y(3)-y(2))-3.1*y(3))];
%Aplicamos el método de Euler Modificado
for n=1:N
k1=F(t(n),y(:,n));
y(:,n+1)=y(:,n)+h*F(t(n)+0.5*h,y(:,n)+0.5*k1*h);
end
%Ahora para calcular los desplazamientos entorno a las posiciones de
%equilibrio de cada masa, sumamos la propia posición de equilibrio a cada
%una
y(1,:)=y(1,:)+ones(1,N+1)*3;
y(2,:)=y(2,:)+ones(1,N+1)*6;
y(3,:)=y(3,:)+ones(1,N+1)*9;
% 4º - Post Proceso:
%Ahora dibujamos el gráfico y añadimos los nombres de los ejes y un título
plot(t,y(1:3,:))
ylabel('Posicion de la masa')
xlabel('Tiempo')
title('Euler Modificado 2')
legend('x1(t)','x2(t)','x3(t)');
La gráfica generada es la siguiente:
Al cambiar el tamaño del paso para ese método seguimos notando una ligera modificación, y sabemos que eso ocurre por que el error global es proporcional al tamaño del paso elegido. Pero el cambio mayor que se puede notar es en las gráficas. Las irregularidades han aumentado de una manero notable con el cambio de método. Eso se debe a que el método de Runge-Kutta es un método de 4º orden, y por eso bastante más preciso que el método de Euler modificado, que es de 1º orden.
3 Efecto de la velocidad
3.1 Supuesto a) - Misma velocidad para las 3 masas
Ahora tenemos resolver el sistema tomando los desplazamientos iniciales de las masas como 0.5 metros hacia la izquierda para la primera masa, 0.8 metros hacia la derecha para la masa 2 y 0.7 metros hacia la izquierda para la masa 3. Las tres masas están ahora con velocidad inicial de 1 m/s en dirección hacia la derecha. Las masas, las constantes elásticas de los muelles y las posiciones de equilibrio siguen las mismas.
De esta forma el problema queda así:: [math]\begin{cases}x_1'=v_1=1 m/s\\x_2'=v_2=1 m/s\\x_3'=v_3=1 m/s\\v_1'=\frac{1}{2.3}(-3.5x_1+2.1(x_2-x_1))\\v_2'=\frac{1}{1.1}(-2.1(x_2-x_1)+2.5(x_3-x_2))\\v_3'=\frac{1}{3.3}(-2.5(x_3-x_2)-3.1x_3)\end{cases}[/math]
Resolveremos utilizando el siguiente método:
- Método del Trapezio
Para eso utilizaremos [math] h = 0.20 m [/math]
Método del Trapecio
El código de Matlab obtenido (Con paso [math] h = 0.20 [/math]) es este:
% 1º - Datos del programa:
t0=0; %Tiempo Inicial
tf=20; %Tiempo Final
y0=[-0.5;0.8;-0.7;1;1;1]; %El orden con el que vamos a obtener las soluciones es x1, x2, x3, x1', x2', x3'
%Los desplazamientos de las masas 1 y 3 son negativos porque son hacia la izquierda,
%y hemos admitido como positiva la dirección hacia la derecha
% 2º - Discretizacion:
h=0.2; %Tamaño del paso
N=round(tf-t0)/h; %Número de intervalos
t=t0:h:tf; %Vector tiempo
% 3º - Metodo numerico (Metodo del Trapezio):
%Creamos el vector de soluciones (vacío)
y=zeros(6,N+1);
%La solución en el inicio es dato
y(:,1)=y0;
%la matriz F, que contiene las ecuaciones del sistema es formada por:
%F(t,y)=A(t)*y+B(t), A y B son matrizes de coeficientes.
%En nuestro caso, la matriz B(t) no existe.
%Por ser un metodo implicito, necesitamos despejar y(n+1).
%Ecuaciones auxiliares para crear la matriz A(t)
a=(3.5+2.1)/2.3;
b=2.1/2.3;
c=2.1/1.1;
d=(2.1+2.5)/1.1;
e=2.5/1.1;
f=2.5/3.2;
g=(2.5+3.1)/3.3;
A=@(t) [0,0,0,1,0,0;0,0,0,0,1,0;0,0,0,0,0,1;-a,b,0,0,0,0;c,-d,e,0,0,0;0,f,-g,0,0,0];
for n=1:N
y(:,n+1)=(eye(6)-h/2*A(t(n+1)))\((eye(6)+h/2*A(t(n)))*y(:,n)+h/2);
end
%Ahora para calcular los desplazamientos entorno a las posiciones de
%equilibrio de cada masa, sumamos la propia posición de equilibrio a cada una
y(1,:)=y(1,:)+ones(1,N+1)*3;
y(2,:)=y(2,:)+ones(1,N+1)*6;
y(3,:)=y(3,:)+ones(1,N+1)*9;
%Ahora añadimos el gráfico y los nombres de los ejes y el título
plot(t,y(1:3,:))
ylabel('Posición de la masa')
xlabel('Tiempo')
title('Método del Trapecio')
legend('x1(t)','x2(t)','x3(t)')
Al cual obtenemos la siguiente gráfica:
Podemos ver en la gráfica anterior que las masas experimentales oscilan de manera irregular. Las tres comienzan el movimiento con la misma velocidad, pero parten de posiciones distintas. De ello derivan dos consecuencias: en primer lugar, las oscilaciones son muy diferentes, si se comparan entre sí, sin embargo, la masa tres se puede asemejar a un movimiento armónico. Por otro lado, se observa que la amplitud de onda es significativamente elevada.
- Para obtener las posiciones relativas de las masas en los primeros 10 segundos, hemos utilizado el Método de Euler modificado con [math] h = 0.05 m [/math]:
Método de Euler modificado
El código de Matlab obtenido con ese método y con el paso elegido es el siguiente:
% 1º - Datos del programa:
t0=0; %Tiempo Inicial
tf=10; %Tiempo Final
y0=[0.4;0.8;-0.7;1;1;1]; %El orden con el que vamos a obtener las soluciones es x1, x2, x3, x1', x2', x3'
%Los desplazamientos de las masas 1 y 3 son negativos porque son hacia la izquierda,
%y hemos admitido como positiva la dirección hacia la derecha
% 2º - Discretizacion:
h=0.05; %Tamaño del paso
N=round(tf-t0)/h; %Número de intervalos
t=t0:h:tf; %Vector tiempo
% 3º - Metodo numerico (Metodo del Trapezio):
%Creamos el vector de soluciones (vacío)
y=zeros(6,N+1);
%La solución en el inicio es dato
y(:,1)=y0;
%Convertimos el sistema en una matriz por el vector 'y', obteniéndose la función F
F=@(t,y) [y(4);y(5);y(6);(1/2.3)*(-3.5*y(1)+2.1*(y(2)-y(1)));(1/1.1)*(-2.1*(y(2)-y(1))+2.5*(y(3)-y(2)));(1/3.3)*(-2.5*(y(3)-y(2))-3.1*y(3))];
%Aplicamos el método de Euler Modificado
for n=1:N
k1=F(t(n),y(:,n));
y(:,n+1)=y(:,n)+h*F(t(n)+0.5*h,y(:,n)+0.5*k1*h);
end
%Ahora para calcular los desplazamientos entorno a las posiciones de
%equilibrio de cada masa, sumamos la propia posición de equilibrio a cada una
y(1,:)=y(1,:)+ones(1,N+1)*3;
y(2,:)=y(2,:)+ones(1,N+1)*6;
y(3,:)=y(3,:)+ones(1,N+1)*9;
%Ahora añadimos el gráfico y los nombres de los ejes y el título
figure(1)
plot(t,y(1:1,:),t,y(3:3,:))
ylabel('Posición relativa de la masa')
xlabel('Tiempo')
title('Método de Euler modificado')
legend('x1(t)','x3(t)')
figure(2)
plot(t,y(1:3,:))
ylabel('Posición relativa de la masa')
xlabel('Tiempo')
title('Método de Euler modificado')
legend('x1(t)','x2(t)','x3(t)')
Y hemos obtenido como gráficas para las masas 1 y 3 la siguiente gráfica:
Y para las tres masas, la siguiente:
En estos gráficos se representa el movimiento de las tres masas durante los primeros 10 segundos de movimiento, sendo el primero especifico de las masas 1 y 3. Esta vez se ha aplicado el método del Euler modificado que tiene un orden de error menor que el método del trapecio por lo que estos datos son más precisos que los obtenidos en el apartado anterior. Aparentemente se observan dos masas que realizan movimientos parejos a los que describe una onda armónica. En el primero de ellos las dos masas presentan similitudes, por el contrario en el segundo gráfico se pone de manifiesto la discordancia entre las dos anteriores y la tercera en cuestión.
3.2 Supuesto b) - Velocidades distintas para las 3 masas
Ahora tomando los desplazamientos iniciales de las masas como 0.4 metros hacia la derecha, y velocidad de 1 m/s hacia la izquierda para la primera masa, 1.1 metros hacia la derecha, sin velocidad para la masa 2 y 0.5 metros hacia la izquierda, con velocidad de 0,5 m/s hacia la derecha para la masa 3. Las masas, las constantes elásticas de los muelles y las posiciones de equilibrio siguen las mismas.
De esta forma el problema queda así:: [math]\begin{cases}x_1'=v_1=-1 m/s\\x_2'=v_2=0 m/s\\x_3'=v_3=0.5 m/s\\v_1'=\frac{1}{2.3}(-3.5x_1+2.1(x_2-x_1))\\v_2'=\frac{1}{1.1}(-2.1(x_2-x_1)+2.5(x_3-x_2))\\v_3'=\frac{1}{3.3}(-2.5(x_3-x_2)-3.1x_3)\end{cases}[/math]
Como anteriormente hemos admitido como dirección positiva hacia la derecha, la velocidad de la masa 1 es negativa por está hacia la izquierda.
Resolveremos, como en el supuesto a), utilizando el siguiente método:
- Método del Trapezio
Para eso utilizaremos [math] h = 0.20 m [/math]
Método del Trapezio
El código de Matlab obtenido (Con paso [math] h = 0.20 [/math]) es este:
% 1º - Datos del programa:
t0=0; %Tiempo Inicial
tf=20; %Tiempo Final
y0=[0.4;1.1;-0.5;-1;0;0.5]; %El orden con el que vamos a obtener las soluciones es x1, x2, x3, x1', x2', x3'
%La velocidad de la masa 1 y el desplazamiento de la masa 3
%son negativos porque son hacia la izquierda,
%y hemos admitido como positiva la dirección hacia la derecha
% 2º - Discretizacion:
h=0.2; %Tamaño del paso
N=round(tf-t0)/h; %Número de intervalos
t=t0:h:tf; %Vector tiempo
% 3º - Metodo numerico (Metodo del Trapezio):
%Creamos el vector de soluciones (vacío)
y=zeros(6,N+1);
%La solución en el inicio es dato
y(:,1)=y0;
%la matriz F, que contiene las ecuaciones del sistema es formada por:
%F(t,y)=A(t)*y+B(t), A y B son matrizes de coeficientes.
%En nuestro caso, la matriz B(t) no existe.
%Por ser un metodo implicito, necesitamos despejar y(n+1).
%Ecuaciones auxiliares para crear la matriz A(t)
a=(3.5+2.1)/2.3;
b=2.1/2.3;
c=2.1/1.1;
d=(2.1+2.5)/1.1;
e=2.5/1.1;
f=2.5/3.2;
g=(2.5+3.1)/3.3;
A=@(t) [0,0,0,1,0,0;0,0,0,0,1,0;0,0,0,0,0,1;-a,b,0,0,0,0;c,-d,e,0,0,0;0,f,-g,0,0,0];
for n=1:N
y(:,n+1)=(eye(6)-h/2*A(t(n+1)))\((eye(6)+h/2*A(t(n)))*y(:,n)+h/2);
end
%Ahora para calcular los desplazamientos entorno a las posiciones de
%equilibrio de cada masa, sumamos la propia posición de equilibrio a cada una
y(1,:)=y(1,:)+ones(1,N+1)*3;
y(2,:)=y(2,:)+ones(1,N+1)*6;
y(3,:)=y(3,:)+ones(1,N+1)*9;
%Ahora añadimos el gráfico y los nombres de los ejes y el título
plot(t,y(1:3,:))
ylabel('Posición de la masa')
xlabel('Tiempo')
title('Método del Trapecio')
legend('x1(t)','x2(t)','x3(t)')
Al cual obtenemos la siguiente gráfica:
En la gráfica anterior se refleja que las masas experimentales oscilan de manera irregular, pero de manera mas homogênea que en el supuesto a). Es posible notar que los picos y los valles de la gráfica no son tan grandes como en el supuesto anterior, indicando que el cambio de posición no se ha dado de manera tan brusca. Importante recordad que las tres comienzan el movimiento con velocidades distintas y posiciones distintas.
- Para obtener las posiciones relativas de las masas en los primeros 10 segundos, hemos utilizado el Método de Euler modificado con [math] h = 0.05 m [/math]:
Método de Euler modificado
El código de Matlab obtenido con ese método y con el paso elegido es el siguiente:
% 1º - Datos del programa:
t0=0; %Tiempo Inicial
tf=20; %Tiempo Final
y0=[0.4;1.1;-0.5;-1;0;0.5]; %El orden con el que vamos a obtener las soluciones es x1, x2, x3, x1', x2', x3'
%La velocidad de la masa 1 y el desplazamiento de la masa 3
%son negativos porque son hacia la izquierda,
%y hemos admitido como positiva la dirección hacia la derecha
% 2º - Discretizacion:
h=0.2; %Tamaño del paso
N=round(tf-t0)/h; %Número de intervalos
t=t0:h:tf; %Vector tiempo
% 3º - Metodo numerico (Metodo del Trapezio):
%Creamos el vector de soluciones (vacío)
y=zeros(6,N+1);
%La solución en el inicio es dato
y(:,1)=y0;
%Convertimos el sistema en una matriz por el vector 'y', obteniéndose la función F
F=@(t,y) [y(4);y(5);y(6);(1/2.3)*(-3.5*y(1)+2.1*(y(2)-y(1)));(1/1.1)*(-2.1*(y(2)-y(1))+2.5*(y(3)-y(2)));(1/3.3)*(-2.5*(y(3)-y(2))-3.1*y(3))];
%Aplicamos el método de Euler Modificado
for n=1:N
k1=F(t(n),y(:,n));
y(:,n+1)=y(:,n)+h*F(t(n)+0.5*h,y(:,n)+0.5*k1*h);
end
%Ahora para calcular los desplazamientos entorno a las posiciones de
%equilibrio de cada masa, sumamos la propia posición de equilibrio a cada una
y(1,:)=y(1,:)+ones(1,N+1)*3;
y(2,:)=y(2,:)+ones(1,N+1)*6;
y(3,:)=y(3,:)+ones(1,N+1)*9;
%Ahora añadimos el gráfico y los nombres de los ejes y el título
figure(1)
plot(t,y(1:1,:),t,y(3:3,:))
ylabel('Posición relativa de la masa')
xlabel('Tiempo')
title('Método de Euler modificado')
legend('x1(t)','x3(t)')
figure(2)
plot(t,y(1:3,:))
ylabel('Posición relativa de la masa')
xlabel('Tiempo')
title('Método de Euler modificado')
legend('x1(t)','x2(t)','x3(t)')
Y hemos obtenido como gráficas para las masas 1 y 3 la siguiente gráfica:
Y para las tres masas, la siguiente:
En los gráficos anteriores se representa el movimiento realizado por las tres masas experimentales durante los primeros 10 segundos. Al aplicar el método de Euler modificado el error que cabe esperar es muy inferior al obtenido en la primera parte del supuesto que nos ocupa. Por otro lado, hay que destacar que, aunque las condiciones de partida son diferentes respecto al apartado anterior, la similitud de las trayectorias descritas por la masa1 y la masa 3 se mantiene. Ambas siguen diferenciándose claramente de la masa restante. Además, se observa que el periodo de oscilación y la amplitud de onda ha aumentado considerablemente.
4 Sistema sumergido en un medio viscoso
Al sumergir el sistema en un medio viscoso, aparece una fuerza de amortiguamiento de resistencia al avance. Esta fuerza tendrá como módulo el producto de una constante de proporcionalidad y la derivada del desplazamiento [math]x_i[/math]. La fuerza de amortiguamiento lleva siempre el signo sentido contrario a la velocidad. La fuerza de resistencia [math] f_0 [/math] vendrá dada por la expresión::
[math] R_i=-μxi [/math]: Donde [math] μ [/math] es la constante de amortiguamiento [N(s/m)], [math] x_i'[/math] es la velocidad de cada masa [m/s] y [math] R_i[/math] es igual a la fuerza de amortiguamiento [N].
Por tanto, el sistema a resolver será similar a los realizados en apartados anteriores, únicamente habrá que añadir la nueva fuerza a las ya presentes en las ecuaciones de la dinámica del sistema.
[math]\begin{cases}m_1x_1''=-k_1x_1+k_2(x_2-x_1)-μx_1'\\m_2x_2''=-k_2(x_2-x_1)+k_3(x_3-x_2)-μx_2'\\m_3x_3''=-k_3(x_3-x_2)-k_4x_3-μx_3'\end{cases}[/math]
Apreciamos que sigue tratándose de un sistema de ecuaciones diferenciales con tres funciones incógnita [math]x_1(t), x_2(t), x_3(t)[/math] que debemos resolver. En este caso aparecen los términos del amortiguamiento que engloban a las primeras derivadas de estas funciones. Para su resolución, hemos reducido a un sistema de primer orden como se hizo en el apartado anterior.
Aplicando el cambio de variable: [math] v_1=x_1' , v_2=x_2' , v_3=x_3' [/math] con lo que se obtiene: [math] v_1'=x_1'' , v_2'=x_2'' , v_3'=x_3'' [/math] De este modo, se ha llegado a un sistema de dimensión 6, de carácter lineal.
[math]\begin{cases}x_1'=v_1\\x_2'=v_2\\x_3'=v_3\\v_1'=\frac{1}{m_1}[-(k_1+k_2)x_1+k_2x_2-μv_1]\\v_2'=\frac{1}{m_2}[μ_2x_1-(k_2+k_3)x_2+k_3x_3-μv_2]\\v_3'=\frac{1}{m_3}[k_3x_2-(k_2+k_3)*x_3-μv_3]\end{cases}[/math]
Para su resolución debemos plantear el método del trapecio con [math] h=0,2 [/math] Primeramente tratamos de escribir el sistema de la forma [math] y'=Ay+T[/math] se aplica el sistema numérico adecuado al método: [math] y_n+1=(I-\frac{h}{2A})\(y_n+\frac{h}{2}(f(t_n,y_n)+T_n+1)[/math] donde la función [math]f[/math] es la matriz de funciones.
Método del Trapecio
Se han tomado las condiciones iniciales del apartado 2.
El código de Matlab utilizado para obtener las posiciones de las masas es el siguiente:
%Primeramente, definiremos todas los datos de masas y constantes recuperadoras
%de los muelles del problema.
m1=2.3; m2=1.1; m3=3.3;
k1=3.5; k2=2.1; k3=2.5; k4=3.1;
%Definimos la longitud del espacio
L=13;
%Las condiciones iniciales son:
x1=0.4; x2=0.9; x3=0.6; v1=0; v2=0; v3=0;
xiniciales=[x1;x2;x3;v1;v2;v3];
%Definimos el coeficiente de amortiguamiento del medio viscoso
c=1;
%El número de ecuaciones con el que trabajaremos será
ne=6;
%Definimos el tamaño de paso a emplear
h=0.2;
%Discretizamos el intervalo de tiempo
t0=0; tN=60;
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%El vector de la solución, que incluye las soluciones de las seis funciones, será de la forma:
x=zeros(6,length(t));
%Al cual, imponemos las condiciones iniciales
x(:,1)=xiniciales;
%Definimos la matriz de funciones F
F=@(t,x)[x(4);x(5);x(6);(-(k1+k2)*x(1)+k2*x(2)-c*x(4))/m1; (k2*x(1)-(k2+k3)*x(2)+k3*x(3)-c*x(5))/m2;(k3*x(2)-(k3+k4)*x(3)-c*x(6))/m3];
%Ahora debemos aplicar el método del trapecio:
%Primeramente tomamos la matriz identidad, la cual nos será útil
I=eye(ne);
%Definimos la matriz A de los coeficientes del sistema
A1=[0 0 0 1 0 0];
A2=[0 0 0 0 1 0];
A3=[0 0 0 0 0 1];
A4=[-(k1+k2)/m1 k2/m1 0 -c/m1 0 0];
A5=[k2/m2 -(k2+k3)/m2 k3/m2 0 -c/m2 0];
A6=[0 k3/m3 -(k3+k4)/m3 0 0 -c/m3];
A=[A1;A2;A3;A4;A5;A6];
for i=1:N
T=zeros(6,1);
x(:,i+1)=(I-(h/2)*A)\(x(:,i)+(h/2)*(F(t(i),x(:,i))+T));
end
%Ahora en la matriz x ya tenemos las seis funciones solución.
%Obtenemos las posiciones de cada masa sumando a cada desplazamiento la
%posición de equilibrio, como hicimos en apartado anteriores
p1=x(1,:)+3;
p2=x(2,:)+6;
p3=x(3,:)+9;
%Gráfica de la solución
hold on
plot(t,p1,'b')
plot(t,p2,'r')
plot(t,p3,'g')
xlabel('Paso del tiempo')
ylabel('Posición de las masas')
legend('Posición de la masa 1','Posición de la masa 2', 'Posición de la masa 3')
hold off
Con el código, hemos obtenido la siguiente gráfica:
En ella se observa que al no haber energía cinética inicial (1) comienza su movimiento acercándose a la posición de equilibrio (2). Para este caso de amortiguamiento, se producen oscilaciones al rededor de la posición de equilibrio. La solución analítica nos da una función armónica multiplicada por una función exponencial decreciente. Esto justifica como se producen las oscilaciones a medida que disminuye la amplitud gracias a la acción del amortiguador.
(1) y estar todas ellas desplazadas a la derecha (2) desplazándose hacia la izquierda Rápidamente se amortigua el movimiento para estas condiciones iniciales, que se hace despreciable a partir de 15 segundos pese a que la oscilación continua infinitamente de manera teórica.
En el caso de dar al sistema mecánico las condiciones iniciales descritas en el apartado 3.a) obtenemos un comienzo de las oscilaciones distintas. En este caso, para una mayor interpretación visual tomamos el eje x como el eje de la posición de las masas y el eje y como el del paso del tiempo.
%Primeramente, definiremos todas los datos de masas y constantes recuperadoras
%de los muelles del problema.
m1=2.3; m2=1.1; m3=3.3;
k1=3.5; k2=2.1; k3=2.5; k4=3.1;
%Definimos la longitud del espacio
L=13;
%Las condiciones iniciales son:
x1=-0.5; x2=0.8; x3=-0.7; v1=1; v2=1; v3=1;
xiniciales=[x1;x2;x3;v1;v2;v3];
%Definimos el coeficiente de amortiguamiento del medio viscoso
c=1;
%El número de ecuaciones con el que trabajaremos será
ne=6;
%Definimos el tamaño de paso a emplear
h=0.2;
%Discretizamos el intervalo de tiempo
t0=0; tN=60;
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%La matriz de la solución, que incluye las soluciones de las seis funciones, será de la forma:
x=zeros(6,length(t));
%Al cual, imponemos las condiciones iniciales
x(:,1)=xiniciales;
%Definimos la matriz de funciones F
F=@(t,x)[x(4);x(5);x(6);(-(k1+k2)*x(1)+k2*x(2)-c*x(4))/m1; (k2*x(1)-(k2+k3)*x(2)+k3*x(3)-c*x(5))/m2;(k3*x(2)-(k3+k4)*x(3)-c*x(6))/m3];
%Ahora debemos aplicar el método del trapecio:
%Primeramente tomamos la matriz identidad, la cual nos será útil
I=eye(ne);
%Definimos la matriz A de los coeficientes del sistema
A1=[0 0 0 1 0 0];
A2=[0 0 0 0 1 0];
A3=[0 0 0 0 0 1];
A4=[-(k1+k2)/m1 k2/m1 0 -c/m1 0 0];
A5=[k2/m2 -(k2+k3)/m2 k3/m2 0 -c/m2 0];
A6=[0 k3/m3 -(k3+k4)/m3 0 0 -c/m3];
A=[A1;A2;A3;A4;A5;A6];
for i=1:N
T=zeros(6,1);
x(:,i+1)=(I-(h/2)*A)\(x(:,i)+(h/2)*(F(t(i),x(:,i))+T));
end
%Ahora en la matriz x ya tenemos las seis funciones solución.
%Obtenemos las posiciones de cada masa sumando a cada desplazamiento la
%posición de equilibrio, como hicimos en apartado anteriores
p1=x(1,:)+3;
p2=x(2,:)+6;
p3=x(3,:)+9;
%Gráfica de la solución
hold on
plot(t,p1,'b')
plot(t,p2,'r')
plot(t,p3,'g')
xlabel('Paso del tiempo')
ylabel('Posición de las masas')
legend('Posición de la masa 1','Posición de la masa 2', 'Posición de la masa 3')
hold off
%Creamos una nueva gráfica para una mejor interpretación visual
figure
hold on
plot(p1,t,'b')
plot(p2,t,'r')
plot(p3,t,'g')
ylabel('Paso del tiempo')
xlabel('Posición de las masas')
legend('Posición de la masa 1','Posición de la masa 2', 'Posición de la masa 3')
hold off
Vemos ahora de mejor forma, como la masa 1 comienza su oscilación 0.5m a la izquierda de su posición de equilibrio, lo cual añadido a su velocidad inicial de 1m/s hacia la derecha, hace que comience su movimiento desplazándose rápidamente a derechas acercándose a la posición de equilibrio [math] p_1=3 [/math].
Para la masa 2, tenemos un desplazamiento de 0.8m hacia la derecha. Sin embargo, su otra condición inicial le proporciona una velocidad de 1m/s también hacia la derecha. Esto hace que, junto con las fuerzas debidas a los muelles y amortiguamiento, esta masa comience su movimiento alejándose de la posición de equilibrio, al contrario de como en un principio el movimiento de las masas a estima pudiera sugerir. Rápidamente ese movimiento se ve frenado, y comienza a acercarse a dicha posición hasta que, tras una serie de oscilaciones de distinta amplitud y 'eje', el movimiento se hace despreciable en la posición de equilibrio.
Finalmente, para la masa 3, por estas desplazada 0.7m hacia la izquierda y tener una velocidad inicial de 1m/s hacia la derecha, comienza su movimiento análogamente hacia la derecha acercándose a la posición de equilibrio hasta que, nuevamente, su oscilación se hace despreciable alrededor de ella con el paso del tiempo.
4.1 Energía Asociada
La energía viene dada por la fórmula: [math] E=E_c+E_p=\frac{1}{2}mx'^{2}+\frac{1}{2}kx^{2} [/math], que particularizada en este caso queda: [math] E=\frac{1}{2}(x'^{2})+x^{2} [/math].
clear all, clc
t0=0;
tf=60;
h=0.2; %Tamaño de paso
t=t0:h:tf;
N=(tf-t0)/h;
y=zeros(2,N+1);
y0=[0.3,0.3]';
y(:,1)=y0;
A=[0 1;-2 -1];
C=inv(eye(size(A))-h/2*A);
for n=1:N
sv=y(2,n)/abs(y(2,n)); % Signo de la velocidad
y(:,n+1)=C*(y(:,n)+h/2*(A*y(:,n)+[0 sv*abs(2*exp(-0.01*t(n))*sin(2*t(n)))]'+[0 sv*abs(2*exp(-0.01*t(n+1))*sin(2*t(n+1)))]'));
end
x=y(1,:); %Posición
v=y(2,:); %Velocidad
e= 0.5*((v.*v)+2*(x.*x)); %Energía
%Gráfica de la solución
hold on
plot(t,e)
semilogy(t,e); %Escala logarítmica
xlabel('Energía')
ylabel('Tiempo')
title('Energía asociada al sistema')
hold off
Obteniendo la siguiente gráfica:
Sí, al aumentar la viscosidad del sistema disminuye la energía. Esto se debe a que obtenemos una solución armónica que depende de la constante de viscosidad, la cual, al aumentar, hace que x disminuya. Como x disminuye, también lo hace la energía total. Además, desde el punto de vista físico esta situación tiene sentido, ya que a mayor viscosidad más oposición al movimiento.
