Diferencia entre revisiones de «Sistema de masas y muelles»
De MateWiki
| (No se muestran 47 ediciones intermedias de 2 usuarios) | |||
| Línea 14: | Línea 14: | ||
Suponemos que en el instante t=0 las tres masas están desplazadas 0.5, 1 y 0.8 metros hacia la derecha de la posici ́on de equilibrio y se sueltan repentinamente, sin velocidad. Vamos a suponer que k1=4N/m, k2=2N/m, k3=1N/m, k4=3N/m, m1=2kg, m2=1kg, m3=3kg, que la distancia entre las paredes es de 12 metros y que en equilibrio las tres masas están en las posiciones 2.5, 4 y 8.Resolveremos el problema numéricamente utilizando los métodos de Euler modificado y Runge Kutta con pasos de h=0.1m y h=0.025m. | Suponemos que en el instante t=0 las tres masas están desplazadas 0.5, 1 y 0.8 metros hacia la derecha de la posici ́on de equilibrio y se sueltan repentinamente, sin velocidad. Vamos a suponer que k1=4N/m, k2=2N/m, k3=1N/m, k4=3N/m, m1=2kg, m2=1kg, m3=3kg, que la distancia entre las paredes es de 12 metros y que en equilibrio las tres masas están en las posiciones 2.5, 4 y 8.Resolveremos el problema numéricamente utilizando los métodos de Euler modificado y Runge Kutta con pasos de h=0.1m y h=0.025m. | ||
| − | Sustituyendo ahora los valores dados arriba en las ecuaciones obtenemos las siguientes, con las cuales operaremos | + | Sustituyendo ahora los valores dados arriba en las ecuaciones obtenemos las siguientes, con las cuales operaremos. Para aplicar el método reduciremos el sistema a uno de primer orden ayudándonos de las siguientes ecuaciones. |
| − | + | ||
| − | + | ||
| − | + | ||
| − | Para aplicar el método reduciremos el sistema a uno de primer orden ayudándonos de las siguientes ecuaciones. | + | |
<math>\dot x=u\\\dot y=v\\\dot z=w</math> | <math>\dot x=u\\\dot y=v\\\dot z=w</math> | ||
| Línea 26: | Línea 22: | ||
<math>\left\{\begin{matrix}m_{1}\ddot x=-k_{1}x+k_{2}(y-x)\\m_{2}\ddot y=-k_{2}(y-x)+k_{3}(z-y)\\m_{3}\ddot z=-k_{3}(z-y)-k_{4}(z)\\\dot x=u\\\dot y=v\\\dot z=w\end{matrix}\right.</math> | <math>\left\{\begin{matrix}m_{1}\ddot x=-k_{1}x+k_{2}(y-x)\\m_{2}\ddot y=-k_{2}(y-x)+k_{3}(z-y)\\m_{3}\ddot z=-k_{3}(z-y)-k_{4}(z)\\\dot x=u\\\dot y=v\\\dot z=w\end{matrix}\right.</math> | ||
| − | + | ===Método de Euler Modificado=== | |
{{matlab|codigo= | {{matlab|codigo= | ||
| − | + | % Euler Modificado h=0.1 | |
% Sistema de 3 masas y 4 muelles | % Sistema de 3 masas y 4 muelles | ||
clear all | clear all | ||
% Intervalo de tiempo | % Intervalo de tiempo | ||
| − | t0=0; | + | t0=0;tN=10; |
| − | tN=10; | + | |
% Datos | % Datos | ||
| − | k1=4; | + | k1=4;k2=2;k3=1;k4=3; |
| − | k2=2; | + | m1=2;m2=1;m3=3; |
| − | k3=1; | + | |
| − | k4=3; | + | |
| − | m1=2; | + | |
| − | m2=1; | + | |
| − | m3=3; | + | |
%Condiciones iniciales | %Condiciones iniciales | ||
| − | x0=0.5; | + | x0=0.5;y0=1;z0=0.8; |
| − | y0=1; | + | u0=0;v0=0;w0=0; |
| − | z0=0.8; | + | |
| − | u0=0; | + | |
| − | v0=0; | + | |
| − | w0=0 | + | |
| − | N=100; % Pasos | + | N=100; % Pasos para h=0.1 |
h=(tN-t0)/N; % Intervalo | h=(tN-t0)/N; % Intervalo | ||
| − | |||
% Primeros valores | % Primeros valores | ||
| − | xx(1)=x0; | + | xx(1)=x0;yy(1)=y0;zz(1)=z0;uu(1)=u0;vv(1)=v0;ww(1)=w0; |
| − | yy(1)=y0; | + | |
| − | zz(1)=z0; | + | |
| − | uu(1)=u0; | + | |
| − | vv(1)=v0; | + | |
| − | ww(1)=w0; | + | |
| − | % | + | for n=1:N |
| + | % Primera K | ||
| + | k1x=uu(n); | ||
| + | k1y=vv(n); | ||
| + | k1z=ww(n); | ||
| + | k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n); | ||
| + | k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n); | ||
| + | k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n); | ||
| + | xa=xx(n)+h/2*k1x; | ||
| + | ya=yy(n)+h/2*k1y; | ||
| + | za=zz(n)+h/2*k1z; | ||
| + | ua=uu(n)+h/2*k1u; | ||
| + | va=vv(n)+h/2*k1v; | ||
| + | wa=ww(n)+h/2*k1w; | ||
| + | |||
| + | % Segunda K | ||
| + | k2x=ua; | ||
| + | k2y=va; | ||
| + | k2z=wa; | ||
| + | k2u=-(k2+k1)/m1*xa+k2/m1*ya; | ||
| + | k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za; | ||
| + | k2w=k3/m3*ya-((k4+k3)/m3)*za; | ||
| + | xx(n+1)=xx(n)+h*k2x; | ||
| + | yy(n+1)=yy(n)+h*k2y; | ||
| + | zz(n+1)=zz(n)+h*k2z; | ||
| + | uu(n+1)=uu(n)+h*k2u; | ||
| + | vv(n+1)=vv(n)+h*k2v; | ||
| + | ww(n+1)=ww(n)+h*k2w; | ||
| + | end | ||
| + | hold on | ||
| + | t=t0:h:tN; | ||
| + | e=t0:0.01:tN | ||
| + | plot(xx+2.5,t,'-') | ||
| + | plot(yy+4,t,'-r') | ||
| + | plot(zz+8,t,'-k') | ||
| + | plot(2.5,e,'-g') | ||
| + | plot(4,e,'-g') | ||
| + | plot(8,e,'-g') | ||
| + | hold off | ||
| + | axis([0,12,0,10]) | ||
| + | hold on | ||
| + | title('Euler Modificado') | ||
| + | legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio') | ||
| + | hold off | ||
| + | }} | ||
| + | |||
| + | Aqui podemos ver los gráficos obtenidos. El de la izquierda corresponde a un paso de 0.1 y el derecho a uno de 0.025. El código es el mismo para ambos a excepción de N=100 para el primero y N=40 para el segundo. | ||
| + | |||
| + | {| | ||
| + | |- | ||
| + | | [[Archivo:em20.1.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,10) con paso de 0.1]] || [[Archivo:em20.25.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,10) con paso de 0.025|600px]] | ||
| + | |} | ||
| + | |||
| + | ===Método de Runge-Kutta=== | ||
| + | |||
| + | El código de Matlab es: | ||
| + | {{matlab|codigo= | ||
| + | |||
| + | % Sistema de 3 masas y 4 muelles | ||
| + | t0=0;tN=10; | ||
| + | % Datos | ||
| + | k1=4;k2=2;k3=1;k4=3; | ||
| + | m1=2;m2=1;m3=3; | ||
| + | %Condiciones iniciales | ||
| + | x0=0.5;y0=1;z0=0.8;u0=0;v0=0;w0=0; | ||
| + | N=100; % Pasos | ||
| + | h=(tN-t0)/N; % Intervalo | ||
| + | |||
| + | % Primeros valores | ||
| + | xx(1)=x0;yy(1)=y0;zz(1)=z0;uu(1)=u0;vv(1)=v0;ww(1)=w0; | ||
| + | |||
| + | % Metodo de Runge-Kutta | ||
for n=1:N | for n=1:N | ||
% Primera K | % Primera K | ||
| Línea 83: | Línea 134: | ||
k2w=k3/m3*(yy(n)+1/2*k1y*h)-((k4+k3)/m3)*(zz(n)+1/2*k1z*h); | k2w=k3/m3*(yy(n)+1/2*k1y*h)-((k4+k3)/m3)*(zz(n)+1/2*k1z*h); | ||
| − | |||
% Tercera K | % Tercera K | ||
k3x=uu(n)+1/2*h*k2u; | k3x=uu(n)+1/2*h*k2u; | ||
| Línea 99: | Línea 149: | ||
k4v=k2/m2*(xx(n)+h*k3x)-((k2+k3)/m2)*(yy(n)+k3y*h)+(k3/m2)*(zz(n)+k3z*h); | k4v=k2/m2*(xx(n)+h*k3x)-((k2+k3)/m2)*(yy(n)+k3y*h)+(k3/m2)*(zz(n)+k3z*h); | ||
k4w=k3/m3*(yy(n)+h*k3y)-((k4+k3)/m3)*(zz(n)+k3z*h); | k4w=k3/m3*(yy(n)+h*k3y)-((k4+k3)/m3)*(zz(n)+k3z*h); | ||
| − | |||
xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x); | xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x); | ||
| Línea 109: | Línea 158: | ||
end | end | ||
| + | e=t0:0.01:tN | ||
t=t0:h:tN; | t=t0:h:tN; | ||
axis([0,12,0,10]) | axis([0,12,0,10]) | ||
| Línea 115: | Línea 165: | ||
plot(xx+2.5,t,'-') | plot(xx+2.5,t,'-') | ||
plot(yy+4,t,'-r') | plot(yy+4,t,'-r') | ||
| − | plot(zz+ | + | plot(zz+8,t,'-k') |
| − | plot(2.5, | + | plot(2.5,e,'-g') |
| − | plot(4, | + | plot(4,e,'-g') |
| − | plot( | + | plot(8,e,'-g') |
legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio') | legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio') | ||
hold off | hold off | ||
}} | }} | ||
| + | Aqui podemos ver los gráficos obtenidos. El superior corresponde a un paso de 0.1 y el inferior a uno de 0.025. El código es el mismo para ambos a excepción de N=100 para el primero y N=40 para el segundo. | ||
| − | [[Archivo: | + | {| |
| + | |- | ||
| + | | [[Archivo:rk20.1.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,10) con paso de 0.1]] || [[Archivo:rk20.25.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,10) con paso de 0.025|600px]] | ||
| + | |} | ||
==Otras condiciones iniciales== | ==Otras condiciones iniciales== | ||
| Línea 129: | Línea 183: | ||
Ahora vamos a imponer unas nuevas CI para representar dos nuevas situaciones. Dichas condiciones son: | Ahora vamos a imponer unas nuevas CI para representar dos nuevas situaciones. Dichas condiciones son: | ||
| − | <math>\left\{\begin{matrix}x(0)=0.5 \\\dot x(0)=1\\y(0)=1 \\\dot y(0)=1\\z(0)=-0.5 \\\dot z=-1\\\end{matrix}\right.</math> | + | {| |
| − | + | |- | |
| − | <math>\left\{\begin{matrix}x(0)=0.5 \\\dot x(0)=-1\\y(0)=1 \\\dot y(0)=0\\z(0)=-0.5 \\\dot z=0.5\\\end{matrix}\right.</math> | + | | A<math>\left\{\begin{matrix}x(0)=0.5 \\\dot x(0)=1\\y(0)=1 \\\dot y(0)=1\\z(0)=-0.5 \\\dot z=-1\\\end{matrix}\right.</math> || B<math>\left\{\begin{matrix}x(0)=0.5 \\\dot x(0)=-1\\y(0)=1 \\\dot y(0)=0\\z(0)=-0.5 \\\dot z=0.5\\\end{matrix}\right.</math> |
| + | |} | ||
===Método del trapecio=== | ===Método del trapecio=== | ||
| + | Utilizando las condiciones iniciales A y con un h=0.2 metros durante los primeros 20 segundos obtenemos: | ||
| + | |||
| + | {{matlab|codigo= | ||
| + | %Sistema EDOS lineal | ||
| + | %Metodo de Trapecio Apartado A | ||
| + | %%datos del problema | ||
| + | t0=0; | ||
| + | tf=20; | ||
| + | y0=[-0.5 1 1 1 -0.5 1]'; | ||
| + | k1=4; | ||
| + | k2=2; | ||
| + | k3=1; | ||
| + | k4=3; | ||
| + | m1=2; | ||
| + | m2=1; | ||
| + | m3=3; | ||
| + | u=-1; | ||
| + | v=1; | ||
| + | w=1; | ||
| + | |||
| + | A=[0 1 0 0 0 0; -(k2+k1)/m1 0 k2/m1 0 0 0; 0 0 0 1 0 0; k2/m2 0 -((k2+k3)/m2) 0 (k3/m2) 0; 0 0 0 0 0 1 ; 0 0 k3/m3 0 -((k4+k3)/m3) 0]; | ||
| + | %%discretizacion | ||
| + | N=100; | ||
| + | h=(tf-t0)/N; | ||
| + | %%Vectores de tiempo y solucion approx. | ||
| + | t=t0:h:tf; | ||
| + | y=zeros(6,N+1);%%matriz de ceros con 3 filas (por no. de ecus)y N+1 columnas | ||
| + | %%Inicializacion | ||
| + | y(:,1)=y0; | ||
| + | yy=y0; | ||
| + | %%necesito gn y gn+1 | ||
| + | for n=1:N | ||
| + | gn=[0 0 0 0 0 0]';%vector con terminos libres de y | ||
| + | gnp1=[0 0 0 0 0 0]'; | ||
| + | yy=(eye(6)-h/2*A)\(((eye(6)+h/2*A)*yy)+h/2*(gn+gnp1)); | ||
| + | y(:,n+1)=yy; | ||
| + | end | ||
| + | clf | ||
| + | |||
| + | e=t0:0.05:tf | ||
| + | hold on | ||
| + | title('Trapecio A') | ||
| + | plot(y(1,:)+2.5,t,'-') | ||
| + | plot(y(1,:)+4,t,'-r') | ||
| + | plot(y(1,:)+8,t,'-k') | ||
| + | plot(2.5,e,'-g') | ||
| + | plot(4,e,'-g') | ||
| + | plot(8,e,'-g') | ||
| + | legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio') | ||
| + | hold off | ||
| + | }} | ||
| + | |||
| + | Si bien usando los mismos parámetros de h y tiempo con las condiciones B tenemos: | ||
| + | |||
| + | {{matlab|codigo= | ||
| + | |||
| + | %Sistema EDOS lineal | ||
| + | %Metodo de Trapecio Apartado B | ||
| + | %%datos del problema | ||
| + | t0=0; | ||
| + | tf=20; | ||
| + | y0=[0.5 -1 1 0 -0.5 0.5]'; | ||
| + | k1=4; | ||
| + | k2=2; | ||
| + | k3=1; | ||
| + | k4=3; | ||
| + | m1=2; | ||
| + | m2=1; | ||
| + | m3=3; | ||
| + | u=-1; | ||
| + | v=0; | ||
| + | w=0.5; | ||
| + | |||
| + | A=[0 1 0 0 0 0; -(k2+k1)/m1 0 k2/m1 0 0 0; 0 0 0 1 0 0; k2/m2 0 -((k2+k3)/m2) 0 (k3/m2) 0; 0 0 0 0 0 1 ; 0 0 k3/m3 0 -((k4+k3)/m3) 0]; | ||
| + | %%discretizacion | ||
| + | N=100; | ||
| + | h=(tf-t0)/N; | ||
| + | %%Vectores de tiempo y solucion approx. | ||
| + | t=t0:h:tf; | ||
| + | y=zeros(6,N+1);%%matriz de ceros con 3 filas (por no. de ecus)y N+1 columnas | ||
| + | %%Inicializacion | ||
| + | y(:,1)=y0; | ||
| + | yy=y0; | ||
| + | %%necesito gn y gn+1 | ||
| + | for n=1:N | ||
| + | gn=[0 0 0 0 0 0]';%vector con terminos libres de y | ||
| + | gnp1=[0 0 0 0 0 0]'; | ||
| + | yy=(eye(6)-h/2*A)\(((eye(6)+h/2*A)*yy)+h/2*(gn+gnp1)); | ||
| + | y(:,n+1)=yy; | ||
| + | end | ||
| + | clf | ||
| + | |||
| + | e=t0:0.05:tf | ||
| + | hold on | ||
| + | title('Trapecio B') | ||
| + | plot(y(1,:)+2.5,t,'-') | ||
| + | plot(y(1,:)+4,t,'-r') | ||
| + | plot(y(1,:)+8,t,'-k') | ||
| + | plot(2.5,e,'-g') | ||
| + | plot(4,e,'-g') | ||
| + | plot(8,e,'-g') | ||
| + | legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio') | ||
| + | hold off | ||
| + | }} | ||
| + | |||
| + | {| | ||
| + | |- | ||
| + | | [[Archivo:tresatrap.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,20) condiciones A]] || [[Archivo:tresbtrap.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,20) condiciones B|600px]] | ||
| + | |} | ||
| + | |||
===Método de Euler Modificado=== | ===Método de Euler Modificado=== | ||
| + | Ahora utilizaremos el método de Euler Modificado junto con las condiciones A con un h=0.05 metros y un tiempo de 10 segundos: | ||
| + | |||
| + | -->Para la primera y tercera masas tenemos: | ||
| + | |||
| + | |||
| + | {{matlab|codigo= | ||
| + | |||
| + | % Euler Modificado Apartado A con Masas 1 y 3 | ||
| + | % Sistema de 3 masas y 4 muelles | ||
| + | clear all | ||
| + | % Intervalo de tiempo | ||
| + | t0=0; | ||
| + | tN=10; | ||
| + | % Datos | ||
| + | k1=4; | ||
| + | k2=2; | ||
| + | k3=1; | ||
| + | k4=3; | ||
| + | m1=2; | ||
| + | m2=1; | ||
| + | m3=3; | ||
| + | %Condiciones iniciales | ||
| + | x0=-0.5; | ||
| + | y0=1; | ||
| + | z0=-0.5; | ||
| + | u0=1; | ||
| + | v0=1; | ||
| + | w0=1; | ||
| + | |||
| + | N=50; % Pasos | ||
| + | h=(tN-t0)/N; % Intervalo | ||
| + | |||
| + | |||
| + | % Primeros valores | ||
| + | xx(1)=x0; | ||
| + | yy(1)=y0; | ||
| + | zz(1)=z0; | ||
| + | uu(1)=u0; | ||
| + | vv(1)=v0; | ||
| + | ww(1)=w0; | ||
| + | |||
| + | for n=1:N | ||
| + | % Primera K | ||
| + | k1x=uu(n); | ||
| + | k1y=vv(n); | ||
| + | k1z=ww(n); | ||
| + | k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n); | ||
| + | k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n); | ||
| + | k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n); | ||
| + | |||
| + | xa=xx(n)+h/2*k1x; | ||
| + | ya=yy(n)+h/2*k1y; | ||
| + | za=zz(n)+h/2*k1z; | ||
| + | ua=uu(n)+h/2*k1u; | ||
| + | va=vv(n)+h/2*k1v; | ||
| + | wa=ww(n)+h/2*k1w; | ||
| + | |||
| + | % Segunda K | ||
| + | k2x=ua; | ||
| + | k2y=va; | ||
| + | k2z=wa; | ||
| + | k2u=-(k2+k1)/m1*xa+k2/m1*ya; | ||
| + | k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za; | ||
| + | k2w=k3/m3*ya-((k4+k3)/m3)*za; | ||
| + | |||
| + | |||
| + | |||
| + | xx(n+1)=xx(n)+h*k2x; | ||
| + | yy(n+1)=yy(n)+h*k2y; | ||
| + | zz(n+1)=zz(n)+h*k2z; | ||
| + | uu(n+1)=uu(n)+h*k2u; | ||
| + | vv(n+1)=vv(n)+h*k2v; | ||
| + | ww(n+1)=ww(n)+h*k2w; | ||
| + | end | ||
| + | hold on | ||
| + | t=t0:h:tN; | ||
| + | plot(zz+8-(xx+2),t,'-k') | ||
| + | hold off | ||
| + | axis([0,12,0,10]) | ||
| + | hold on | ||
| + | title('Euler Modificado Apartado A con Masas 1 y 3') | ||
| + | |||
| + | legend('Dif entre m3,m1') | ||
| + | hold off | ||
| + | }} | ||
| + | |||
| + | -->Para todas las masas tenemos: | ||
| + | |||
| + | {{matlab|codigo= | ||
| + | |||
| + | % Euler Modificado Apartado A con Todas las Masas | ||
| + | % Sistema de 3 masas y 4 muelles | ||
| + | clear all | ||
| + | % Intervalo de tiempo | ||
| + | t0=0; | ||
| + | tN=10; | ||
| + | % Datos | ||
| + | k1=4; | ||
| + | k2=2; | ||
| + | k3=1; | ||
| + | k4=3; | ||
| + | m1=2; | ||
| + | m2=1; | ||
| + | m3=3; | ||
| + | %Condiciones iniciales | ||
| + | x0=-0.5; | ||
| + | y0=1; | ||
| + | z0=-0.5; | ||
| + | u0=1; | ||
| + | v0=1; | ||
| + | w0=1; | ||
| + | |||
| + | N=50; % Pasos | ||
| + | h=(tN-t0)/N; % Intervalo | ||
| + | |||
| + | |||
| + | % Primeros valores | ||
| + | xx(1)=x0; | ||
| + | yy(1)=y0; | ||
| + | zz(1)=z0; | ||
| + | uu(1)=u0; | ||
| + | vv(1)=v0; | ||
| + | ww(1)=w0; | ||
| + | |||
| + | |||
| + | for n=1:N | ||
| + | % Primera K | ||
| + | k1x=uu(n); | ||
| + | k1y=vv(n); | ||
| + | k1z=ww(n); | ||
| + | k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n); | ||
| + | k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n); | ||
| + | k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n); | ||
| + | |||
| + | xa=xx(n)+h/2*k1x; | ||
| + | ya=yy(n)+h/2*k1y; | ||
| + | za=zz(n)+h/2*k1z; | ||
| + | ua=uu(n)+h/2*k1u; | ||
| + | va=vv(n)+h/2*k1v; | ||
| + | wa=ww(n)+h/2*k1w; | ||
| + | |||
| + | % Segunda K | ||
| + | k2x=ua; | ||
| + | k2y=va; | ||
| + | k2z=wa; | ||
| + | k2u=-(k2+k1)/m1*xa+k2/m1*ya; | ||
| + | k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za; | ||
| + | k2w=k3/m3*ya-((k4+k3)/m3)*za; | ||
| + | |||
| + | |||
| + | |||
| + | xx(n+1)=xx(n)+h*k2x; | ||
| + | yy(n+1)=yy(n)+h*k2y; | ||
| + | zz(n+1)=zz(n)+h*k2z; | ||
| + | uu(n+1)=uu(n)+h*k2u; | ||
| + | vv(n+1)=vv(n)+h*k2v; | ||
| + | ww(n+1)=ww(n)+h*k2w; | ||
| + | end | ||
| + | hold on | ||
| + | t=t0:h:tN; | ||
| + | plot((yy+4)-(xx+2.5),t,'-k') | ||
| + | plot((zz+8)-(yy+4)+((yy+4)-(xx+2.5)),t,'-g') | ||
| + | hold off | ||
| + | axis([0,12,0,10]) | ||
| + | hold on | ||
| + | title('Euler Modificado') | ||
| + | |||
| + | legend('Dif entre m1,m2','Dif entre m2,m3') | ||
| + | hold off | ||
| + | }} | ||
| + | |||
| + | {| | ||
| + | |- | ||
| + | | [[Archivo:tresa13.png|thumb|600px|left|Diferencia relativa entre las posiciones m1,m3]] || [[Archivo:tresa123.png|thumb|600px|left|Diferencia relativa entre las posiciones m1,m2,m3|600px]] | ||
| + | |} | ||
| + | |||
| + | Variando a continuación solo las condiciones iniciales utilizando las condiciones de B: | ||
| + | |||
| + | -->Para la primera y tercera masas tenemos: | ||
| + | |||
| + | {{matlab|codigo= | ||
| + | % Euler Modificado Apartado B con Masas 1 y 3 | ||
| + | % Sistema de 3 masas y 4 muelles | ||
| + | clear all | ||
| + | % Intervalo de tiempo | ||
| + | t0=0; | ||
| + | tN=10; | ||
| + | % Datos | ||
| + | k1=4; | ||
| + | k2=2; | ||
| + | k3=1; | ||
| + | k4=3; | ||
| + | m1=2; | ||
| + | m2=1; | ||
| + | m3=3; | ||
| + | %Condiciones iniciales | ||
| + | x0=0.5; | ||
| + | y0=1; | ||
| + | z0=-0.5; | ||
| + | u0=-1; | ||
| + | v0=0; | ||
| + | w0=0.5; | ||
| + | |||
| + | N=50; % Pasos | ||
| + | h=(tN-t0)/N; % Intervalo | ||
| + | |||
| + | |||
| + | % Primeros valores | ||
| + | xx(1)=x0; | ||
| + | yy(1)=y0; | ||
| + | zz(1)=z0; | ||
| + | uu(1)=u0; | ||
| + | vv(1)=v0; | ||
| + | ww(1)=w0; | ||
| + | |||
| + | |||
| + | |||
| + | for n=1:N | ||
| + | % Primera K | ||
| + | k1x=uu(n); | ||
| + | k1y=vv(n); | ||
| + | k1z=ww(n); | ||
| + | k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n); | ||
| + | k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n); | ||
| + | k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n); | ||
| + | |||
| + | xa=xx(n)+h/2*k1x; | ||
| + | ya=yy(n)+h/2*k1y; | ||
| + | za=zz(n)+h/2*k1z; | ||
| + | ua=uu(n)+h/2*k1u; | ||
| + | va=vv(n)+h/2*k1v; | ||
| + | wa=ww(n)+h/2*k1w; | ||
| + | |||
| + | % Segunda K | ||
| + | k2x=ua; | ||
| + | k2y=va; | ||
| + | k2z=wa; | ||
| + | k2u=-(k2+k1)/m1*xa+k2/m1*ya; | ||
| + | k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za; | ||
| + | k2w=k3/m3*ya-((k4+k3)/m3)*za; | ||
| + | |||
| + | |||
| + | |||
| + | xx(n+1)=xx(n)+h*k2x; | ||
| + | yy(n+1)=yy(n)+h*k2y; | ||
| + | zz(n+1)=zz(n)+h*k2z; | ||
| + | uu(n+1)=uu(n)+h*k2u; | ||
| + | vv(n+1)=vv(n)+h*k2v; | ||
| + | ww(n+1)=ww(n)+h*k2w; | ||
| + | end | ||
| + | hold on | ||
| + | t=t0:h:tN; | ||
| + | plot(zz+8-(xx+2),t,'-k') | ||
| + | hold off | ||
| + | axis([0,12,0,10]) | ||
| + | hold on | ||
| + | title('Euler Modificado Apartado B con Masas 1 y 3') | ||
| + | |||
| + | legend('Dif entre m3,m1') | ||
| + | hold off | ||
| + | }} | ||
| + | -->Para todas las masas tenemos: | ||
| + | |||
| + | {{matlab|codigo= | ||
| + | |||
| + | % Euler Modificado Apartado B con Todas las Masas | ||
| + | % Sistema de 3 masas y 4 muelles | ||
| + | clear all | ||
| + | % Intervalo de tiempo | ||
| + | t0=0; | ||
| + | tN=10; | ||
| + | % Datos | ||
| + | k1=4; | ||
| + | k2=2; | ||
| + | k3=1; | ||
| + | k4=3; | ||
| + | m1=2; | ||
| + | m2=1; | ||
| + | m3=3; | ||
| + | %Condiciones iniciales | ||
| + | x0=0.5; | ||
| + | y0=1; | ||
| + | z0=-0.5; | ||
| + | u0=-1; | ||
| + | v0=0; | ||
| + | w0=0.5; | ||
| + | |||
| + | N=50; % Pasos | ||
| + | h=(tN-t0)/N; % Intervalo | ||
| + | |||
| + | |||
| + | % Primeros valores | ||
| + | xx(1)=x0; | ||
| + | yy(1)=y0; | ||
| + | zz(1)=z0; | ||
| + | uu(1)=u0; | ||
| + | vv(1)=v0; | ||
| + | ww(1)=w0; | ||
| + | |||
| + | |||
| + | for n=1:N | ||
| + | % Primera K | ||
| + | k1x=uu(n); | ||
| + | k1y=vv(n); | ||
| + | k1z=ww(n); | ||
| + | k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n); | ||
| + | k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n); | ||
| + | k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n); | ||
| + | |||
| + | xa=xx(n)+h/2*k1x; | ||
| + | ya=yy(n)+h/2*k1y; | ||
| + | za=zz(n)+h/2*k1z; | ||
| + | ua=uu(n)+h/2*k1u; | ||
| + | va=vv(n)+h/2*k1v; | ||
| + | wa=ww(n)+h/2*k1w; | ||
| + | |||
| + | % Segunda K | ||
| + | k2x=ua; | ||
| + | k2y=va; | ||
| + | k2z=wa; | ||
| + | k2u=-(k2+k1)/m1*xa+k2/m1*ya; | ||
| + | k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za; | ||
| + | k2w=k3/m3*ya-((k4+k3)/m3)*za; | ||
| + | |||
| + | |||
| + | |||
| + | xx(n+1)=xx(n)+h*k2x; | ||
| + | yy(n+1)=yy(n)+h*k2y; | ||
| + | zz(n+1)=zz(n)+h*k2z; | ||
| + | uu(n+1)=uu(n)+h*k2u; | ||
| + | vv(n+1)=vv(n)+h*k2v; | ||
| + | ww(n+1)=ww(n)+h*k2w; | ||
| + | end | ||
| + | hold on | ||
| + | t=t0:h:tN; | ||
| + | plot((yy+4)-(xx+2.5),t,'-k') | ||
| + | plot((zz+8)-(yy+4)+((yy+4)-(xx+2.5)),t,'-g') | ||
| + | hold off | ||
| + | axis([0,12,0,10]) | ||
| + | hold on | ||
| + | title('Euler Modificado') | ||
| + | |||
| + | legend('Dif entre m1,m2','Dif entre m2,m3') | ||
| + | hold off | ||
| + | }} | ||
| + | |||
| + | {| | ||
| + | |- | ||
| + | | [[Archivo:tresb13.png|thumb|600px|left|Diferencia relativa entre las posiciones m1,m3]] || [[Archivo:tresb123.png|thumb|600px|left|Diferencia relativa entre las posiciones m1,m2,m3|600px]] | ||
| + | |} | ||
==Estudio de un sistema en medio viscoso y con fuerza exterior== | ==Estudio de un sistema en medio viscoso y con fuerza exterior== | ||
| Línea 146: | Línea 661: | ||
<math>\left\{\begin{matrix}\dot v+v+4x-2e^{(-0.01t)·sin(2t)}=0\\\dot x=v\end{matrix}\right.</math> | <math>\left\{\begin{matrix}\dot v+v+4x-2e^{(-0.01t)·sin(2t)}=0\\\dot x=v\end{matrix}\right.</math> | ||
| − | Si operamos el sistema aplicando la fórmula del trapecio y despejamos dejando como incógnita los términos en (n+1) resolvemos con el siguiente código: | + | Si operamos el sistema aplicando la fórmula del trapecio y despejamos dejando como incógnita los términos en (n+1) resolvemos con el siguiente código de paso h=0.2: |
{{matlab|codigo= | {{matlab|codigo= | ||
| − | a=0 | + | a=0 |
| − | N= | + | b=60 |
| + | N=300 | ||
h=(b-a)/N | h=(b-a)/N | ||
t=a:h:b; | t=a:h:b; | ||
| − | t0=a | + | t0=a |
| − | x(1)=x0; v(1)=v0; E(1)=0. | + | x0=0.3 |
| + | v0=0.3 | ||
| + | k=4; | ||
| + | c=1; | ||
| + | x(1)=x0;v(1)=v0;E(1)=2*(x0^2)+0.5*(v0^2) | ||
for n=1:N | for n=1:N | ||
| − | A=[x(n)-h/2*v(n)+h/2*(2*exp(-0.01*t(n))*sin(2*t(n)))+2*exp(-0.01*t(n+1)*sin(2*t(n+1))); | + | A=[x(n)-c*h/2*v(n)+h/2*k*x(n)+h/2*(2*exp(-0.01*t(n))*sin(2*t(n)))+2*exp(-0.01*t(n+1)*sin(2*t(n+1))); |
| − | (1+h/2)*v(n)] | + | c*(1+h/2)*v(n)] |
| − | B=[1, | + | B=[1-h/2*k,c*h/2;0,c] |
rest=A\B; | rest=A\B; | ||
x(n+1)=rest(1); | x(n+1)=rest(1); | ||
v(n+1)=rest(2); | v(n+1)=rest(2); | ||
| − | |||
end | end | ||
| − | + | E=2*(x.^2)+0.5*(v.^2); | |
| + | semilogx(t,E,'g') | ||
}} | }} | ||
| + | Ahora vemos la ganancia de energía del sistema a medida que transcurre el tiempo, en escala logarítmica. Primero vemos para c=1 y a continuación para c=3. Es bastante razonable que se disipe más energía con c más alta porque la viscosidad del medio es mayor y sería necesaria una mayor energía para mantener el movimiento. Sin embargo en nuestro caso se aprecia que aumenta y se puede deber a que la fuerza externa aporta más energía al sistema de la que disipa el medio viscoso. | ||
| + | |||
| + | [[Archivo:energiac1.png|marco|centro|Energía del sistema en el primer minuto para c=1]] || [[Archivo:energiac3.png|marco|centro|Energía del sistema en el primer minuto para c=3]] | ||
| + | |||
[[Categoría:Ecuaciones Diferenciales]] | [[Categoría:Ecuaciones Diferenciales]] | ||
[[Categoría:ED13/14]] | [[Categoría:ED13/14]] | ||
[[Categoría:Trabajos 2013-14]] | [[Categoría:Trabajos 2013-14]] | ||
Revisión actual del 00:55 6 mar 2014
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Sistema de masas y muelles. Grupo 11-B |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2013-14 |
| Autores | Javier Mellado, Jacobo Campos, Javier Chamizo, Miguel García, Alfonso Ascanio |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Sistemas de muelles y masas. Estos sistemas son bastante útiles a la hora de estudiar vibraciones, consiguiendo traducirlo a un lenguaje matemático que sea fácil de estudiar e interpretar.
Contenido
1 Introducción
En este artículo vamos a estudiar un sistema de 3 masas y 4 muelles en disposición horizontal, entre dos paredes rígidas. Empezaremos estudiándolo primero desde la posición de equilibrio en el caso de no haber rozamiento ni amortiguamiento, por lo que solo actúan las fuerzas restauradoras de los muelles. Tomamos como variables las posiciones de cada masa (x,y,z) en función del tiempo. A partir de ellas escribimos las ecuaciones diferenciales correspondientes a una interpretación según la mecánica de Newton.
Sistema formado por 3 masas y 4 muelles,anclado a las paredes laterales
[math]\left\{\begin{matrix}m_{1}\ddot x=-k_{1}x+k_{2}(y-x)\\m_{2}\ddot y=-k_{2}(y-x)+k_{3}(z-y)\\m_{3}\ddot z=-k_{3}(z-y)-k_{4}(z)\end{matrix}\right.[/math]
2 Condiciones iniciales y caso particular
Suponemos que en el instante t=0 las tres masas están desplazadas 0.5, 1 y 0.8 metros hacia la derecha de la posici ́on de equilibrio y se sueltan repentinamente, sin velocidad. Vamos a suponer que k1=4N/m, k2=2N/m, k3=1N/m, k4=3N/m, m1=2kg, m2=1kg, m3=3kg, que la distancia entre las paredes es de 12 metros y que en equilibrio las tres masas están en las posiciones 2.5, 4 y 8.Resolveremos el problema numéricamente utilizando los métodos de Euler modificado y Runge Kutta con pasos de h=0.1m y h=0.025m.
Sustituyendo ahora los valores dados arriba en las ecuaciones obtenemos las siguientes, con las cuales operaremos. Para aplicar el método reduciremos el sistema a uno de primer orden ayudándonos de las siguientes ecuaciones.
[math]\dot x=u\\\dot y=v\\\dot z=w[/math]
Y lo juntamos con el ya existente para obtener:
[math]\left\{\begin{matrix}m_{1}\ddot x=-k_{1}x+k_{2}(y-x)\\m_{2}\ddot y=-k_{2}(y-x)+k_{3}(z-y)\\m_{3}\ddot z=-k_{3}(z-y)-k_{4}(z)\\\dot x=u\\\dot y=v\\\dot z=w\end{matrix}\right.[/math]
2.1 Método de Euler Modificado
% Euler Modificado h=0.1
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;tN=10;
% Datos
k1=4;k2=2;k3=1;k4=3;
m1=2;m2=1;m3=3;
%Condiciones iniciales
x0=0.5;y0=1;z0=0.8;
u0=0;v0=0;w0=0;
N=100; % Pasos para h=0.1
h=(tN-t0)/N; % Intervalo
% Primeros valores
xx(1)=x0;yy(1)=y0;zz(1)=z0;uu(1)=u0;vv(1)=v0;ww(1)=w0;
for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);
xa=xx(n)+h/2*k1x;
ya=yy(n)+h/2*k1y;
za=zz(n)+h/2*k1z;
ua=uu(n)+h/2*k1u;
va=vv(n)+h/2*k1v;
wa=ww(n)+h/2*k1w;
% Segunda K
k2x=ua;
k2y=va;
k2z=wa;
k2u=-(k2+k1)/m1*xa+k2/m1*ya;
k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
k2w=k3/m3*ya-((k4+k3)/m3)*za;
xx(n+1)=xx(n)+h*k2x;
yy(n+1)=yy(n)+h*k2y;
zz(n+1)=zz(n)+h*k2z;
uu(n+1)=uu(n)+h*k2u;
vv(n+1)=vv(n)+h*k2v;
ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
e=t0:0.01:tN
plot(xx+2.5,t,'-')
plot(yy+4,t,'-r')
plot(zz+8,t,'-k')
plot(2.5,e,'-g')
plot(4,e,'-g')
plot(8,e,'-g')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado')
legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio')
hold off
Aqui podemos ver los gráficos obtenidos. El de la izquierda corresponde a un paso de 0.1 y el derecho a uno de 0.025. El código es el mismo para ambos a excepción de N=100 para el primero y N=40 para el segundo.
2.2 Método de Runge-Kutta
El código de Matlab es:
% Sistema de 3 masas y 4 muelles
t0=0;tN=10;
% Datos
k1=4;k2=2;k3=1;k4=3;
m1=2;m2=1;m3=3;
%Condiciones iniciales
x0=0.5;y0=1;z0=0.8;u0=0;v0=0;w0=0;
N=100; % Pasos
h=(tN-t0)/N; % Intervalo
% Primeros valores
xx(1)=x0;yy(1)=y0;zz(1)=z0;uu(1)=u0;vv(1)=v0;ww(1)=w0;
% Metodo de Runge-Kutta
for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);
% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2z=ww(n)+1/2*h*k1w;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h)+(k3/m2)*(zz(n)+1/2*k1z*h);
k2w=k3/m3*(yy(n)+1/2*k1y*h)-((k4+k3)/m3)*(zz(n)+1/2*k1z*h);
% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3z=ww(n)+1/2*h*k2w;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h)+(k3/m2)*(zz(n)+1/2*k2z*h);
k3w=k3/m3*(yy(n)+1/2*k2y*h)-((k4+k3)/m3)*(zz(n)+1/2*k2z*h);
% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4z=ww(n)+h*k3w;
k4u=-(k2+k1)/m1*(xx(n)+h*k3x)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h*k3x)-((k2+k3)/m2)*(yy(n)+k3y*h)+(k3/m2)*(zz(n)+k3z*h);
k4w=k3/m3*(yy(n)+h*k3y)-((k4+k3)/m3)*(zz(n)+k3z*h);
xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
zz(n+1)=zz(n)+(h/6)*(k1z+2*k2z+2*k3z+k4z);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
ww(n+1)=ww(n)+(h/6)*(k1w+2*k2w+2*k3w+k4w);
end
e=t0:0.01:tN
t=t0:h:tN;
axis([0,12,0,10])
hold on
title('Runge-Kutta')
plot(xx+2.5,t,'-')
plot(yy+4,t,'-r')
plot(zz+8,t,'-k')
plot(2.5,e,'-g')
plot(4,e,'-g')
plot(8,e,'-g')
legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio')
hold offAqui podemos ver los gráficos obtenidos. El superior corresponde a un paso de 0.1 y el inferior a uno de 0.025. El código es el mismo para ambos a excepción de N=100 para el primero y N=40 para el segundo.
3 Otras condiciones iniciales
Ahora vamos a imponer unas nuevas CI para representar dos nuevas situaciones. Dichas condiciones son:
| A[math]\left\{\begin{matrix}x(0)=0.5 \\\dot x(0)=1\\y(0)=1 \\\dot y(0)=1\\z(0)=-0.5 \\\dot z=-1\\\end{matrix}\right.[/math] | B[math]\left\{\begin{matrix}x(0)=0.5 \\\dot x(0)=-1\\y(0)=1 \\\dot y(0)=0\\z(0)=-0.5 \\\dot z=0.5\\\end{matrix}\right.[/math] |
3.1 Método del trapecio
Utilizando las condiciones iniciales A y con un h=0.2 metros durante los primeros 20 segundos obtenemos:
%Sistema EDOS lineal
%Metodo de Trapecio Apartado A
%%datos del problema
t0=0;
tf=20;
y0=[-0.5 1 1 1 -0.5 1]';
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
u=-1;
v=1;
w=1;
A=[0 1 0 0 0 0; -(k2+k1)/m1 0 k2/m1 0 0 0; 0 0 0 1 0 0; k2/m2 0 -((k2+k3)/m2) 0 (k3/m2) 0; 0 0 0 0 0 1 ; 0 0 k3/m3 0 -((k4+k3)/m3) 0];
%%discretizacion
N=100;
h=(tf-t0)/N;
%%Vectores de tiempo y solucion approx.
t=t0:h:tf;
y=zeros(6,N+1);%%matriz de ceros con 3 filas (por no. de ecus)y N+1 columnas
%%Inicializacion
y(:,1)=y0;
yy=y0;
%%necesito gn y gn+1
for n=1:N
gn=[0 0 0 0 0 0]';%vector con terminos libres de y
gnp1=[0 0 0 0 0 0]';
yy=(eye(6)-h/2*A)\(((eye(6)+h/2*A)*yy)+h/2*(gn+gnp1));
y(:,n+1)=yy;
end
clf
e=t0:0.05:tf
hold on
title('Trapecio A')
plot(y(1,:)+2.5,t,'-')
plot(y(1,:)+4,t,'-r')
plot(y(1,:)+8,t,'-k')
plot(2.5,e,'-g')
plot(4,e,'-g')
plot(8,e,'-g')
legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio')
hold off
Si bien usando los mismos parámetros de h y tiempo con las condiciones B tenemos:
%Sistema EDOS lineal
%Metodo de Trapecio Apartado B
%%datos del problema
t0=0;
tf=20;
y0=[0.5 -1 1 0 -0.5 0.5]';
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
u=-1;
v=0;
w=0.5;
A=[0 1 0 0 0 0; -(k2+k1)/m1 0 k2/m1 0 0 0; 0 0 0 1 0 0; k2/m2 0 -((k2+k3)/m2) 0 (k3/m2) 0; 0 0 0 0 0 1 ; 0 0 k3/m3 0 -((k4+k3)/m3) 0];
%%discretizacion
N=100;
h=(tf-t0)/N;
%%Vectores de tiempo y solucion approx.
t=t0:h:tf;
y=zeros(6,N+1);%%matriz de ceros con 3 filas (por no. de ecus)y N+1 columnas
%%Inicializacion
y(:,1)=y0;
yy=y0;
%%necesito gn y gn+1
for n=1:N
gn=[0 0 0 0 0 0]';%vector con terminos libres de y
gnp1=[0 0 0 0 0 0]';
yy=(eye(6)-h/2*A)\(((eye(6)+h/2*A)*yy)+h/2*(gn+gnp1));
y(:,n+1)=yy;
end
clf
e=t0:0.05:tf
hold on
title('Trapecio B')
plot(y(1,:)+2.5,t,'-')
plot(y(1,:)+4,t,'-r')
plot(y(1,:)+8,t,'-k')
plot(2.5,e,'-g')
plot(4,e,'-g')
plot(8,e,'-g')
legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio')
hold off
3.2 Método de Euler Modificado
Ahora utilizaremos el método de Euler Modificado junto con las condiciones A con un h=0.05 metros y un tiempo de 10 segundos:
-->Para la primera y tercera masas tenemos:
% Euler Modificado Apartado A con Masas 1 y 3
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
%Condiciones iniciales
x0=-0.5;
y0=1;
z0=-0.5;
u0=1;
v0=1;
w0=1;
N=50; % Pasos
h=(tN-t0)/N; % Intervalo
% Primeros valores
xx(1)=x0;
yy(1)=y0;
zz(1)=z0;
uu(1)=u0;
vv(1)=v0;
ww(1)=w0;
for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);
xa=xx(n)+h/2*k1x;
ya=yy(n)+h/2*k1y;
za=zz(n)+h/2*k1z;
ua=uu(n)+h/2*k1u;
va=vv(n)+h/2*k1v;
wa=ww(n)+h/2*k1w;
% Segunda K
k2x=ua;
k2y=va;
k2z=wa;
k2u=-(k2+k1)/m1*xa+k2/m1*ya;
k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
k2w=k3/m3*ya-((k4+k3)/m3)*za;
xx(n+1)=xx(n)+h*k2x;
yy(n+1)=yy(n)+h*k2y;
zz(n+1)=zz(n)+h*k2z;
uu(n+1)=uu(n)+h*k2u;
vv(n+1)=vv(n)+h*k2v;
ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
plot(zz+8-(xx+2),t,'-k')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado Apartado A con Masas 1 y 3')
legend('Dif entre m3,m1')
hold off
-->Para todas las masas tenemos:
% Euler Modificado Apartado A con Todas las Masas
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
%Condiciones iniciales
x0=-0.5;
y0=1;
z0=-0.5;
u0=1;
v0=1;
w0=1;
N=50; % Pasos
h=(tN-t0)/N; % Intervalo
% Primeros valores
xx(1)=x0;
yy(1)=y0;
zz(1)=z0;
uu(1)=u0;
vv(1)=v0;
ww(1)=w0;
for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);
xa=xx(n)+h/2*k1x;
ya=yy(n)+h/2*k1y;
za=zz(n)+h/2*k1z;
ua=uu(n)+h/2*k1u;
va=vv(n)+h/2*k1v;
wa=ww(n)+h/2*k1w;
% Segunda K
k2x=ua;
k2y=va;
k2z=wa;
k2u=-(k2+k1)/m1*xa+k2/m1*ya;
k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
k2w=k3/m3*ya-((k4+k3)/m3)*za;
xx(n+1)=xx(n)+h*k2x;
yy(n+1)=yy(n)+h*k2y;
zz(n+1)=zz(n)+h*k2z;
uu(n+1)=uu(n)+h*k2u;
vv(n+1)=vv(n)+h*k2v;
ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
plot((yy+4)-(xx+2.5),t,'-k')
plot((zz+8)-(yy+4)+((yy+4)-(xx+2.5)),t,'-g')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado')
legend('Dif entre m1,m2','Dif entre m2,m3')
hold off
Variando a continuación solo las condiciones iniciales utilizando las condiciones de B:
-->Para la primera y tercera masas tenemos:
% Euler Modificado Apartado B con Masas 1 y 3
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
%Condiciones iniciales
x0=0.5;
y0=1;
z0=-0.5;
u0=-1;
v0=0;
w0=0.5;
N=50; % Pasos
h=(tN-t0)/N; % Intervalo
% Primeros valores
xx(1)=x0;
yy(1)=y0;
zz(1)=z0;
uu(1)=u0;
vv(1)=v0;
ww(1)=w0;
for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);
xa=xx(n)+h/2*k1x;
ya=yy(n)+h/2*k1y;
za=zz(n)+h/2*k1z;
ua=uu(n)+h/2*k1u;
va=vv(n)+h/2*k1v;
wa=ww(n)+h/2*k1w;
% Segunda K
k2x=ua;
k2y=va;
k2z=wa;
k2u=-(k2+k1)/m1*xa+k2/m1*ya;
k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
k2w=k3/m3*ya-((k4+k3)/m3)*za;
xx(n+1)=xx(n)+h*k2x;
yy(n+1)=yy(n)+h*k2y;
zz(n+1)=zz(n)+h*k2z;
uu(n+1)=uu(n)+h*k2u;
vv(n+1)=vv(n)+h*k2v;
ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
plot(zz+8-(xx+2),t,'-k')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado Apartado B con Masas 1 y 3')
legend('Dif entre m3,m1')
hold off-->Para todas las masas tenemos:
% Euler Modificado Apartado B con Todas las Masas
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
%Condiciones iniciales
x0=0.5;
y0=1;
z0=-0.5;
u0=-1;
v0=0;
w0=0.5;
N=50; % Pasos
h=(tN-t0)/N; % Intervalo
% Primeros valores
xx(1)=x0;
yy(1)=y0;
zz(1)=z0;
uu(1)=u0;
vv(1)=v0;
ww(1)=w0;
for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);
xa=xx(n)+h/2*k1x;
ya=yy(n)+h/2*k1y;
za=zz(n)+h/2*k1z;
ua=uu(n)+h/2*k1u;
va=vv(n)+h/2*k1v;
wa=ww(n)+h/2*k1w;
% Segunda K
k2x=ua;
k2y=va;
k2z=wa;
k2u=-(k2+k1)/m1*xa+k2/m1*ya;
k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
k2w=k3/m3*ya-((k4+k3)/m3)*za;
xx(n+1)=xx(n)+h*k2x;
yy(n+1)=yy(n)+h*k2y;
zz(n+1)=zz(n)+h*k2z;
uu(n+1)=uu(n)+h*k2u;
vv(n+1)=vv(n)+h*k2v;
ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
plot((yy+4)-(xx+2.5),t,'-k')
plot((zz+8)-(yy+4)+((yy+4)-(xx+2.5)),t,'-g')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado')
legend('Dif entre m1,m2','Dif entre m2,m3')
hold off
4 Estudio de un sistema en medio viscoso y con fuerza exterior
En muchos de los casos que estudian vibraciones es necesario tener en cuenta un amortiguador que va frenando las vibraciones. Vamos a estudiar el caso de una única masa (m=1kg), [math]\mu=c=1[/math] , k=4 y una fuerza exterior [math]f(t)=2·exp(-0.01t)·sin(2t)[/math]. Las condiciones iniciales son x(0)=0.3 v(0)=0.3 Aplicando las ecuaciones de Newton obtenemos la siguiente ecuación diferencial:
- [math]\ddot x+\dot x +kx -2e^{(-0.01t)·sin(2t)}=0[/math]
Para resolver numericamente este problema aplicamos el método del trapecio. Necesitamos pasar de la ecuación de segundo orden a un sistema de dos ecuaciones de primer orden, el cual queda así:
[math]\left\{\begin{matrix}\dot v+v+4x-2e^{(-0.01t)·sin(2t)}=0\\\dot x=v\end{matrix}\right.[/math]
Si operamos el sistema aplicando la fórmula del trapecio y despejamos dejando como incógnita los términos en (n+1) resolvemos con el siguiente código de paso h=0.2:
a=0
b=60
N=300
h=(b-a)/N
t=a:h:b;
t0=a
x0=0.3
v0=0.3
k=4;
c=1;
x(1)=x0;v(1)=v0;E(1)=2*(x0^2)+0.5*(v0^2)
for n=1:N
A=[x(n)-c*h/2*v(n)+h/2*k*x(n)+h/2*(2*exp(-0.01*t(n))*sin(2*t(n)))+2*exp(-0.01*t(n+1)*sin(2*t(n+1)));
c*(1+h/2)*v(n)]
B=[1-h/2*k,c*h/2;0,c]
rest=A\B;
x(n+1)=rest(1);
v(n+1)=rest(2);
end
E=2*(x.^2)+0.5*(v.^2);
semilogx(t,E,'g')
Ahora vemos la ganancia de energía del sistema a medida que transcurre el tiempo, en escala logarítmica. Primero vemos para c=1 y a continuación para c=3. Es bastante razonable que se disipe más energía con c más alta porque la viscosidad del medio es mayor y sería necesaria una mayor energía para mantener el movimiento. Sin embargo en nuestro caso se aprecia que aumenta y se puede deber a que la fuerza externa aporta más energía al sistema de la que disipa el medio viscoso.
Energía del sistema en el primer minuto para c=1
Energía del sistema en el primer minuto para c=3
