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 reduzimos 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
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]
3.1 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.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:
