MateWiki - Contribuciones del usuario [es]

Modelo de vibraciones con muelles masas (Grupo A8)

2016-05-02T18:56:53Z

Danieldelpotro: /* Supuesto 2 */

{{ TrabajoED | Modelo de vibraciones con muelles masas. Grupo 8-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Daniel del Potro Gabin 1423
Humberto del Castillo 1281
Luis de la Fuente Puig 939
Ignacio del Río 1693
}}
<big>
== '''Enunciado y modelización del sistema''' ==
</big>
Muchos modelos de vibraciones se modelizan con sistemas de muelles masas. Consideramos aquí un ejemplo simple de cuatro muelles con constantes kiN/m; i= 1;2;3;4, (fuerzas restauradoras), y tres masas mi; i= 1;2;3, unidos a dos paredes que se deslizan libremente sobre un plano horizontal.

Se modeliza el sistema de ecuaciones diferenciales para el desplazamiento de 3 masas desde la posición de equilibrio. Los grados de libertad que se utilizan son x1, x2 y x3, siendo estos la distancia de cada masa a su posición de equilibrio (Se toma el desplazamiento positivo hacia la derecha). 
Se desprecia la fuerza de rozamiento, por lo que se supone que deslizan sobre una superficie horizontal y lisa, y por tanto las únicas fuerzas a tener en cuenta son las ejercidas por los muelles. Según la ley de Hooke, la fuerza producida por un muelle es proporcional a su elongación. 
Por tanto las fuerzas aplicadas a la primera masa serán: :
<math> F_1=-k_1x_1 </math> :
<math> F_2=k_2(x_2-x_1) </math>
Sobre la segunda masa: :
<math> F_3=-k_2(x_2-x_1) </math> :
<math> F_4=-k_3x_2 </math>
Y las fuerzas que actúan en la tercera masa: :
<math> F_5=-k_3(x_3-x_2) </math> :
<math> F_6=-k_4x_3 </math>

Donde k1, k2, k3 y k4 son las constantes elásticas del movimiento de cada muelle.

Aplicando la segunda Ley de Newton, <math> F=mx'' </math>, se obtiene el siguiente sistema de ecuaciones diferenciales:
<math> m_1x_1''=-k_1x_1+k_2(x_2-x_1)</math>:
<math> m_2x_2''=-k_2(x_2-x_1)-k_3(x_2-x_3)</math>:
<math> m_3x_3''=-k_3(x_3-x_2)-k_4x_3</math>

=='''RESOLUCIÓN DEL SISTEMA DE ECUACIONES'''==

En el instante t=0, las tres masas están desplazadas 0.5,1 y 0.8 metros hacia la derecha respecto de la posición de equilibrio y se sueltan repentinamente,sin velocidad.Se resuelve el sistema de ecuaciones diferenciales utilizando dos métodos numéricos: 
* Método de Euler modificado. 
* Método de Runge-Kutta. 

Primero planteamos el problema en Matlab:

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=0.8;x11=0;x22=0;x33=0;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z=[z1;z2];
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=input('Introduce instante final:');
h=input('Introduce tamaño de paso:');%Probar para h=0.1 y h=0.025
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%f=M*z son las solucines del sistema
%Selecciono la segunda mitad del vector solución, ya que contiene las
%posiciones y la parte superior las derivadas
}}

Y lo resolvemos con ambos métodos:

{{matlab| codigo=
%Método Euler modificado
zEM=z;
for i=1:N
K1=M*zEM(:,i);
zEM(:,i+1)=zEM(:,i)+h*(M*(zEM(:,i)+K1*h/2));
end
zEM=zEM([4 5 6],:);%Selecciono las soluciones de posición
zEM1=zEM(1,:)+2.5;%Posición masa 1
zEM2=zEM(2,:)+4;%Posición masa 2
zEM3=zEM(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(zEM1,t,zEM2,t,zEM3,t);
title('Método Euler Modificado');
hold off

%Método Runge-Kutta
zRK=z;
for j=1:N
K1=M*zRK(:,j);
K2=M*(zRK(:,j)+(1/2)*K1*h);
K3=M*(zRK(:,j)+(1/2)*K2*h);
K4=M*(zRK(:,j)+K3*h);
zRK(:,j+1)=zRK(:,j)+h/6*(K1+2*K2+2*K3+K4);
end
zRK=zRK([4 5 6],:);%Selecciono las soluciones de posición
zRK1=zRK(1,:)+2.5;%Posición masa 1
zRK2=zRK(2,:)+4;%Posición masa 2
zRK3=zRK(3,:)+8;%Posición masa 3
ylabel('Tiempo');xlabel('Posición');
hold on
subplot(1,2,2)
plot(zRK1,t,zRK2,t,zRK3,t);
title('Método Runge-Kutta orden 4');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo para Euler-Modificado y Runge-Kutta con h=0.1:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones h=0.1]]

Y lo comparamos con la gráfica obtenida con Euler-Modificado y Runge-Kutta con h=0.025:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones h=0.025]]

=='''EFECTO DE LA VELOCIDAD INICIAL EN LA TRAYECTORIA'''==

Del mismo modo que en el caso anterior, las tres masas parten de una posición diferente a la del equilibrio, pero en este caso se impone una velocidad inicial. Se analiza la trayectoria de las masas en dos supuestos distintos. 

===Supuesto 1===

La primera masa se desplaza a su izquierda 0.5 metros, la segunda 1 metro a su derecha, y la tercera 0.5 metros hacia la izquierda, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha.. 

Se obtienen los desplazamientos utilizando el método del trapecio.
{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=-0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=20;%Tiempo final
h=0.2;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método del Trapecio
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
z(:,i+1)=inv(I-h/2*M)*(z(:,i)+h/2*(M*z(:,i)));%Despejando
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
plot(z1,t,z2,t,z3,t);
title('Método del trapecio');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

===Supuesto 2===

La primera masa se desplaza a su derecha 0.5 metros con una velocidad de 1 m/s en dirección hacia la izquierda, la segunda masa se suelta sin velocidad inicial con un desplazamiento de 1 metro a la derecha, y la tercera masa se desplaza a su izquierda 0.5 metros y se suelta con una velocidad de 0.5 metros hacia la derecha.

Para obtener los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio, y para hallar las posiciones relativas de las masas mediante el uso del método de Euler Modificado, el código de MATLAB para ambos casos sería el mismo que en el supuesto 1, cambiando las condiciones iniciales

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=-0.5;x11=-1;x22=0;x33=0.5;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=10;%Tiempo final
h=0.05;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método de Euler Modificado
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
K1=M*z(:,i);
z(:,i+1)=z(:,i)+h*(M*(z(:,i)+K1*h/2));
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(z1,t,z3,t);
title('Masas 1 y 3');
subplot(1,2,2);
plot(z1,t,z2,t,z3,t);
title('Masas 1,2 y 3');
hold off
}}

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''SISTEMA DE UNA MASA EN MEDIO VISCOSO CON FUERZA EXTERIOR APLICADA'''==

{{matlab| codigo=
%Datos iniciales
m=1;%Masa
k=4;%Constante del muelle
x1=0.3;x11=0.3;%Datos en t0
x0=[x1;x11]';%Condiciones iniciales
t0=0;%Instante inicial
tN=60;%Instante final
nu=1;%Constante de amortiguamiento
h=0.2;%Tamaño de paso
N=(tN-t0)/h;
%Matriz de coeficientes
M=[-nu/m,-k/m;1,0];
%vector tiempo y matriz posicion
t=linspace(t0,tN,N+1);
x=zeros(2,N+1);
x(:,1)=x0;
%Para conocer el valor de la fuerza en cada instante t
f=[];
for i=1:N+1
f=[f;2*exp(-0.01*t(i))*sin(2*t(i))];
end
%Método del Trapecio
I=eye(size(M));
for i=1:N
x(:,i+1)=inv(I-h/2*M)*(x(:,i)+h/2*(M*x(:,i)+f(i)+f(i+1)));
end
subplot(2,1,1);
plot(t,x(2,:));%Tiempo respecto a velocidad
ylabel('Tiempo');xlabel('Velocidad');
E=k/2*(x(1,:).^2)+m/2*(x(2,:).^2);%Ecuación de la energía
subplot(2,1,2);
semilogx(t,E);
title('Cantidad de energía');
}}

La gráfica obtenida es la siguiente:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Modelo de vibraciones con muelles masas (Grupo A8)

2016-05-02T18:54:05Z

Danieldelpotro: /* Supuesto 1 */

{{ TrabajoED | Modelo de vibraciones con muelles masas. Grupo 8-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Daniel del Potro Gabin 1423
Humberto del Castillo 1281
Luis de la Fuente Puig 939
Ignacio del Río 1693
}}
<big>
== '''Enunciado y modelización del sistema''' ==
</big>
Muchos modelos de vibraciones se modelizan con sistemas de muelles masas. Consideramos aquí un ejemplo simple de cuatro muelles con constantes kiN/m; i= 1;2;3;4, (fuerzas restauradoras), y tres masas mi; i= 1;2;3, unidos a dos paredes que se deslizan libremente sobre un plano horizontal.

Se modeliza el sistema de ecuaciones diferenciales para el desplazamiento de 3 masas desde la posición de equilibrio. Los grados de libertad que se utilizan son x1, x2 y x3, siendo estos la distancia de cada masa a su posición de equilibrio (Se toma el desplazamiento positivo hacia la derecha). 
Se desprecia la fuerza de rozamiento, por lo que se supone que deslizan sobre una superficie horizontal y lisa, y por tanto las únicas fuerzas a tener en cuenta son las ejercidas por los muelles. Según la ley de Hooke, la fuerza producida por un muelle es proporcional a su elongación. 
Por tanto las fuerzas aplicadas a la primera masa serán: :
<math> F_1=-k_1x_1 </math> :
<math> F_2=k_2(x_2-x_1) </math>
Sobre la segunda masa: :
<math> F_3=-k_2(x_2-x_1) </math> :
<math> F_4=-k_3x_2 </math>
Y las fuerzas que actúan en la tercera masa: :
<math> F_5=-k_3(x_3-x_2) </math> :
<math> F_6=-k_4x_3 </math>

Donde k1, k2, k3 y k4 son las constantes elásticas del movimiento de cada muelle.

Aplicando la segunda Ley de Newton, <math> F=mx'' </math>, se obtiene el siguiente sistema de ecuaciones diferenciales:
<math> m_1x_1''=-k_1x_1+k_2(x_2-x_1)</math>:
<math> m_2x_2''=-k_2(x_2-x_1)-k_3(x_2-x_3)</math>:
<math> m_3x_3''=-k_3(x_3-x_2)-k_4x_3</math>

=='''RESOLUCIÓN DEL SISTEMA DE ECUACIONES'''==

En el instante t=0, las tres masas están desplazadas 0.5,1 y 0.8 metros hacia la derecha respecto de la posición de equilibrio y se sueltan repentinamente,sin velocidad.Se resuelve el sistema de ecuaciones diferenciales utilizando dos métodos numéricos: 
* Método de Euler modificado. 
* Método de Runge-Kutta. 

Primero planteamos el problema en Matlab:

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=0.8;x11=0;x22=0;x33=0;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z=[z1;z2];
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=input('Introduce instante final:');
h=input('Introduce tamaño de paso:');%Probar para h=0.1 y h=0.025
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%f=M*z son las solucines del sistema
%Selecciono la segunda mitad del vector solución, ya que contiene las
%posiciones y la parte superior las derivadas
}}

Y lo resolvemos con ambos métodos:

{{matlab| codigo=
%Método Euler modificado
zEM=z;
for i=1:N
K1=M*zEM(:,i);
zEM(:,i+1)=zEM(:,i)+h*(M*(zEM(:,i)+K1*h/2));
end
zEM=zEM([4 5 6],:);%Selecciono las soluciones de posición
zEM1=zEM(1,:)+2.5;%Posición masa 1
zEM2=zEM(2,:)+4;%Posición masa 2
zEM3=zEM(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(zEM1,t,zEM2,t,zEM3,t);
title('Método Euler Modificado');
hold off

%Método Runge-Kutta
zRK=z;
for j=1:N
K1=M*zRK(:,j);
K2=M*(zRK(:,j)+(1/2)*K1*h);
K3=M*(zRK(:,j)+(1/2)*K2*h);
K4=M*(zRK(:,j)+K3*h);
zRK(:,j+1)=zRK(:,j)+h/6*(K1+2*K2+2*K3+K4);
end
zRK=zRK([4 5 6],:);%Selecciono las soluciones de posición
zRK1=zRK(1,:)+2.5;%Posición masa 1
zRK2=zRK(2,:)+4;%Posición masa 2
zRK3=zRK(3,:)+8;%Posición masa 3
ylabel('Tiempo');xlabel('Posición');
hold on
subplot(1,2,2)
plot(zRK1,t,zRK2,t,zRK3,t);
title('Método Runge-Kutta orden 4');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo para Euler-Modificado y Runge-Kutta con h=0.1:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones h=0.1]]

Y lo comparamos con la gráfica obtenida con Euler-Modificado y Runge-Kutta con h=0.025:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones h=0.025]]

=='''EFECTO DE LA VELOCIDAD INICIAL EN LA TRAYECTORIA'''==

Del mismo modo que en el caso anterior, las tres masas parten de una posición diferente a la del equilibrio, pero en este caso se impone una velocidad inicial. Se analiza la trayectoria de las masas en dos supuestos distintos. 

===Supuesto 1===

La primera masa se desplaza a su izquierda 0.5 metros, la segunda 1 metro a su derecha, y la tercera 0.5 metros hacia la izquierda, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha.. 

Se obtienen los desplazamientos utilizando el método del trapecio.
{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=-0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=20;%Tiempo final
h=0.2;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método del Trapecio
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
z(:,i+1)=inv(I-h/2*M)*(z(:,i)+h/2*(M*z(:,i)));%Despejando
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
plot(z1,t,z2,t,z3,t);
title('Método del trapecio');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

===Supuesto 2===

La primera masa se desplaza a su derecha 0.5 metros con una velocidad de 1 m/s en dirección hacia la izquierda, la segunda masa se suelta sin velocidad inicial con un desplazamiento de 1 metro a la derecha, y la tercera masa se desplaza a su izquierda 0.5 metros y se suelta con una velocidad de 0.5 metros hacia la derecha.

Para obtener los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio, y para hallar las posiciones relativas de las masas mediante el uso del método de Euler Modificado, el código de MATLAB para ambos casos sería el mismo que en el supuesto 1, cambiando las condiciones iniciales

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=10;%Tiempo final
h=0.05;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método de Euler Modificado
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
K1=M*z(:,i);
z(:,i+1)=z(:,i)+h*(M*(z(:,i)+K1*h/2));
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(z1,t,z3,t);
title('Masas 1 y 3');
subplot(1,2,2);
plot(z1,t,z2,t,z3,t);
title('Masas 1,2 y 3');
hold off
}}

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''SISTEMA DE UNA MASA EN MEDIO VISCOSO CON FUERZA EXTERIOR APLICADA'''==

{{matlab| codigo=
%Datos iniciales
m=1;%Masa
k=4;%Constante del muelle
x1=0.3;x11=0.3;%Datos en t0
x0=[x1;x11]';%Condiciones iniciales
t0=0;%Instante inicial
tN=60;%Instante final
nu=1;%Constante de amortiguamiento
h=0.2;%Tamaño de paso
N=(tN-t0)/h;
%Matriz de coeficientes
M=[-nu/m,-k/m;1,0];
%vector tiempo y matriz posicion
t=linspace(t0,tN,N+1);
x=zeros(2,N+1);
x(:,1)=x0;
%Para conocer el valor de la fuerza en cada instante t
f=[];
for i=1:N+1
f=[f;2*exp(-0.01*t(i))*sin(2*t(i))];
end
%Método del Trapecio
I=eye(size(M));
for i=1:N
x(:,i+1)=inv(I-h/2*M)*(x(:,i)+h/2*(M*x(:,i)+f(i)+f(i+1)));
end
subplot(2,1,1);
plot(t,x(2,:));%Tiempo respecto a velocidad
ylabel('Tiempo');xlabel('Velocidad');
E=k/2*(x(1,:).^2)+m/2*(x(2,:).^2);%Ecuación de la energía
subplot(2,1,2);
semilogx(t,E);
title('Cantidad de energía');
}}

La gráfica obtenida es la siguiente:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Modelo de vibraciones con muelles masas (Grupo A8)

2016-05-02T18:53:03Z

Danieldelpotro: /* RESOLUCIÓN DEL SISTEMA DE ECUACIONES */

{{ TrabajoED | Modelo de vibraciones con muelles masas. Grupo 8-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Daniel del Potro Gabin 1423
Humberto del Castillo 1281
Luis de la Fuente Puig 939
Ignacio del Río 1693
}}
<big>
== '''Enunciado y modelización del sistema''' ==
</big>
Muchos modelos de vibraciones se modelizan con sistemas de muelles masas. Consideramos aquí un ejemplo simple de cuatro muelles con constantes kiN/m; i= 1;2;3;4, (fuerzas restauradoras), y tres masas mi; i= 1;2;3, unidos a dos paredes que se deslizan libremente sobre un plano horizontal.

Se modeliza el sistema de ecuaciones diferenciales para el desplazamiento de 3 masas desde la posición de equilibrio. Los grados de libertad que se utilizan son x1, x2 y x3, siendo estos la distancia de cada masa a su posición de equilibrio (Se toma el desplazamiento positivo hacia la derecha). 
Se desprecia la fuerza de rozamiento, por lo que se supone que deslizan sobre una superficie horizontal y lisa, y por tanto las únicas fuerzas a tener en cuenta son las ejercidas por los muelles. Según la ley de Hooke, la fuerza producida por un muelle es proporcional a su elongación. 
Por tanto las fuerzas aplicadas a la primera masa serán: :
<math> F_1=-k_1x_1 </math> :
<math> F_2=k_2(x_2-x_1) </math>
Sobre la segunda masa: :
<math> F_3=-k_2(x_2-x_1) </math> :
<math> F_4=-k_3x_2 </math>
Y las fuerzas que actúan en la tercera masa: :
<math> F_5=-k_3(x_3-x_2) </math> :
<math> F_6=-k_4x_3 </math>

Donde k1, k2, k3 y k4 son las constantes elásticas del movimiento de cada muelle.

Aplicando la segunda Ley de Newton, <math> F=mx'' </math>, se obtiene el siguiente sistema de ecuaciones diferenciales:
<math> m_1x_1''=-k_1x_1+k_2(x_2-x_1)</math>:
<math> m_2x_2''=-k_2(x_2-x_1)-k_3(x_2-x_3)</math>:
<math> m_3x_3''=-k_3(x_3-x_2)-k_4x_3</math>

=='''RESOLUCIÓN DEL SISTEMA DE ECUACIONES'''==

En el instante t=0, las tres masas están desplazadas 0.5,1 y 0.8 metros hacia la derecha respecto de la posición de equilibrio y se sueltan repentinamente,sin velocidad.Se resuelve el sistema de ecuaciones diferenciales utilizando dos métodos numéricos: 
* Método de Euler modificado. 
* Método de Runge-Kutta. 

Primero planteamos el problema en Matlab:

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=0.8;x11=0;x22=0;x33=0;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z=[z1;z2];
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=input('Introduce instante final:');
h=input('Introduce tamaño de paso:');%Probar para h=0.1 y h=0.025
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%f=M*z son las solucines del sistema
%Selecciono la segunda mitad del vector solución, ya que contiene las
%posiciones y la parte superior las derivadas
}}

Y lo resolvemos con ambos métodos:

{{matlab| codigo=
%Método Euler modificado
zEM=z;
for i=1:N
K1=M*zEM(:,i);
zEM(:,i+1)=zEM(:,i)+h*(M*(zEM(:,i)+K1*h/2));
end
zEM=zEM([4 5 6],:);%Selecciono las soluciones de posición
zEM1=zEM(1,:)+2.5;%Posición masa 1
zEM2=zEM(2,:)+4;%Posición masa 2
zEM3=zEM(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(zEM1,t,zEM2,t,zEM3,t);
title('Método Euler Modificado');
hold off

%Método Runge-Kutta
zRK=z;
for j=1:N
K1=M*zRK(:,j);
K2=M*(zRK(:,j)+(1/2)*K1*h);
K3=M*(zRK(:,j)+(1/2)*K2*h);
K4=M*(zRK(:,j)+K3*h);
zRK(:,j+1)=zRK(:,j)+h/6*(K1+2*K2+2*K3+K4);
end
zRK=zRK([4 5 6],:);%Selecciono las soluciones de posición
zRK1=zRK(1,:)+2.5;%Posición masa 1
zRK2=zRK(2,:)+4;%Posición masa 2
zRK3=zRK(3,:)+8;%Posición masa 3
ylabel('Tiempo');xlabel('Posición');
hold on
subplot(1,2,2)
plot(zRK1,t,zRK2,t,zRK3,t);
title('Método Runge-Kutta orden 4');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo para Euler-Modificado y Runge-Kutta con h=0.1:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones h=0.1]]

Y lo comparamos con la gráfica obtenida con Euler-Modificado y Runge-Kutta con h=0.025:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones h=0.025]]

=='''EFECTO DE LA VELOCIDAD INICIAL EN LA TRAYECTORIA'''==

Del mismo modo que en el caso anterior, las tres masas parten de una posición diferente a la del equilibrio, pero en este caso se impone una velocidad inicial. Se analiza la trayectoria de las masas en dos supuestos distintos. 

===Supuesto 1===

La primera masa se desplaza a su izquierda 0.5 metros, la segunda 1 metro a su derecha, y la tercera 0.5 metros hacia la derecha, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha.. 

Se obtienen los desplazamientos utilizando el método del trapecio.
{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=-0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=20;%Tiempo final
h=0.2;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método del Trapecio
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
z(:,i+1)=inv(I-h/2*M)*(z(:,i)+h/2*(M*z(:,i)));%Despejando
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
plot(z1,t,z2,t,z3,t);
title('Método del trapecio');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

===Supuesto 2===

La primera masa se desplaza a su derecha 0.5 metros con una velocidad de 1 m/s en dirección hacia la izquierda, la segunda masa se suelta sin velocidad inicial con un desplazamiento de 1 metro a la derecha, y la tercera masa se desplaza a su izquierda 0.5 metros y se suelta con una velocidad de 0.5 metros hacia la derecha.

Para obtener los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio, y para hallar las posiciones relativas de las masas mediante el uso del método de Euler Modificado, el código de MATLAB para ambos casos sería el mismo que en el supuesto 1, cambiando las condiciones iniciales

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=10;%Tiempo final
h=0.05;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método de Euler Modificado
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
K1=M*z(:,i);
z(:,i+1)=z(:,i)+h*(M*(z(:,i)+K1*h/2));
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(z1,t,z3,t);
title('Masas 1 y 3');
subplot(1,2,2);
plot(z1,t,z2,t,z3,t);
title('Masas 1,2 y 3');
hold off
}}

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''SISTEMA DE UNA MASA EN MEDIO VISCOSO CON FUERZA EXTERIOR APLICADA'''==

{{matlab| codigo=
%Datos iniciales
m=1;%Masa
k=4;%Constante del muelle
x1=0.3;x11=0.3;%Datos en t0
x0=[x1;x11]';%Condiciones iniciales
t0=0;%Instante inicial
tN=60;%Instante final
nu=1;%Constante de amortiguamiento
h=0.2;%Tamaño de paso
N=(tN-t0)/h;
%Matriz de coeficientes
M=[-nu/m,-k/m;1,0];
%vector tiempo y matriz posicion
t=linspace(t0,tN,N+1);
x=zeros(2,N+1);
x(:,1)=x0;
%Para conocer el valor de la fuerza en cada instante t
f=[];
for i=1:N+1
f=[f;2*exp(-0.01*t(i))*sin(2*t(i))];
end
%Método del Trapecio
I=eye(size(M));
for i=1:N
x(:,i+1)=inv(I-h/2*M)*(x(:,i)+h/2*(M*x(:,i)+f(i)+f(i+1)));
end
subplot(2,1,1);
plot(t,x(2,:));%Tiempo respecto a velocidad
ylabel('Tiempo');xlabel('Velocidad');
E=k/2*(x(1,:).^2)+m/2*(x(2,:).^2);%Ecuación de la energía
subplot(2,1,2);
semilogx(t,E);
title('Cantidad de energía');
}}

La gráfica obtenida es la siguiente:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Modelo de vibraciones con muelles masas (Grupo A8)

2016-05-02T18:45:52Z

Danieldelpotro: /* RESOLUCIÓN DEL SISTEMA DE ECUACIONES */

{{ TrabajoED | Modelo de vibraciones con muelles masas. Grupo 8-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Daniel del Potro Gabin 1423
Humberto del Castillo 1281
Luis de la Fuente Puig 939
Ignacio del Río 1693
}}
<big>
== '''Enunciado y modelización del sistema''' ==
</big>
Muchos modelos de vibraciones se modelizan con sistemas de muelles masas. Consideramos aquí un ejemplo simple de cuatro muelles con constantes kiN/m; i= 1;2;3;4, (fuerzas restauradoras), y tres masas mi; i= 1;2;3, unidos a dos paredes que se deslizan libremente sobre un plano horizontal.

Se modeliza el sistema de ecuaciones diferenciales para el desplazamiento de 3 masas desde la posición de equilibrio. Los grados de libertad que se utilizan son x1, x2 y x3, siendo estos la distancia de cada masa a su posición de equilibrio (Se toma el desplazamiento positivo hacia la derecha). 
Se desprecia la fuerza de rozamiento, por lo que se supone que deslizan sobre una superficie horizontal y lisa, y por tanto las únicas fuerzas a tener en cuenta son las ejercidas por los muelles. Según la ley de Hooke, la fuerza producida por un muelle es proporcional a su elongación. 
Por tanto las fuerzas aplicadas a la primera masa serán: :
<math> F_1=-k_1x_1 </math> :
<math> F_2=k_2(x_2-x_1) </math>
Sobre la segunda masa: :
<math> F_3=-k_2(x_2-x_1) </math> :
<math> F_4=-k_3x_2 </math>
Y las fuerzas que actúan en la tercera masa: :
<math> F_5=-k_3(x_3-x_2) </math> :
<math> F_6=-k_4x_3 </math>

Donde k1, k2, k3 y k4 son las constantes elásticas del movimiento de cada muelle.

Aplicando la segunda Ley de Newton, <math> F=mx'' </math>, se obtiene el siguiente sistema de ecuaciones diferenciales:
<math> m_1x_1''=-k_1x_1+k_2(x_2-x_1)</math>:
<math> m_2x_2''=-k_2(x_2-x_1)-k_3(x_2-x_3)</math>:
<math> m_3x_3''=-k_3(x_3-x_2)-k_4x_3</math>

=='''RESOLUCIÓN DEL SISTEMA DE ECUACIONES'''==

En el instante t=0, las tres masas están desplazadas 0.5,1 y 0.8 metros hacia la derecha respecto de la posición de equilibrio y se sueltan repentinamente,sin velocidad.Se resuelve el sistema de ecuaciones diferenciales utilizando dos métodos numéricos: 
* Método de Euler modificado. 
* Método de Runge-Kutta. 

Primero planteamos el problema en Matlab:

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=0.8;x11=0;x22=0;x33=0;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z=[z1;z2];
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=input('Introduce instante final:');
h=input('Introduce tamaño de paso:');%Probar para h=0.1 y h=0.025
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%f=M*z son las solucines del sistema
%Selecciono la segunda mitad del vector solución, ya que contiene las
%posiciones y la parte superior las derivadas
}}

Y lo resolvemos con ambos métodos:

{{matlab| codigo=
%Método Euler modificado
zEM=z;
for i=1:N
K1=M*zEM(:,i);
zEM(:,i+1)=zEM(:,i)+h*(M*zEM(:,i)+K1*h/2);
end
zEM=zEM([4 5 6],:);%Selecciono las soluciones de posición
zEM1=zEM(1,:)+2.5;%Posición masa 1
zEM2=zEM(2,:)+4;%Posición masa 2
zEM3=zEM(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(zEM1,t,zEM2,t,zEM3,t);
title('Método Euler Modificado');
hold off

%Método Runge-Kutta
zRK=z;
for j=1:N
K1=M*zRK(:,j);
K2=M*(zRK(:,j)+(1/2)*K1*h);
K3=M*(zRK(:,j)+(1/2)*K2*h);
K4=M*(zRK(:,j)+K3*h);
zRK(:,j+1)=zRK(:,j)+h/6*(K1+2*K2+2*K3+K4);
end
zRK=zRK([4 5 6],:);%Selecciono las soluciones de posición
zRK1=zRK(1,:)+2.5;%Posición masa 1
zRK2=zRK(2,:)+4;%Posición masa 2
zRK3=zRK(3,:)+8;%Posición masa 3
ylabel('Tiempo');xlabel('Posición');
hold on
subplot(1,2,2)
plot(zRK1,t,zRK2,t,zRK3,t);
title('Método Runge-Kutta orden 4');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo para Euler-Modificado y Runge-Kutta con h=0.1:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones h=0.1]]

Y lo comparamos con la gráfica obtenida con Euler-Modificado y Runge-Kutta con h=0.025:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones h=0.025]]

=='''EFECTO DE LA VELOCIDAD INICIAL EN LA TRAYECTORIA'''==

Del mismo modo que en el caso anterior, las tres masas parten de una posición diferente a la del equilibrio, pero en este caso se impone una velocidad inicial. Se analiza la trayectoria de las masas en dos supuestos distintos. 

===Supuesto 1===

La primera masa se desplaza a su izquierda 0.5 metros, la segunda 1 metro a su derecha, y la tercera 0.5 metros hacia la derecha, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha.. 

Se obtienen los desplazamientos utilizando el método del trapecio.
{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=-0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=20;%Tiempo final
h=0.2;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método del Trapecio
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
z(:,i+1)=inv(I-h/2*M)*(z(:,i)+h/2*(M*z(:,i)));%Despejando
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
plot(z1,t,z2,t,z3,t);
title('Método del trapecio');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

===Supuesto 2===

La primera masa se desplaza a su derecha 0.5 metros con una velocidad de 1 m/s en dirección hacia la izquierda, la segunda masa se suelta sin velocidad inicial con un desplazamiento de 1 metro a la derecha, y la tercera masa se desplaza a su izquierda 0.5 metros y se suelta con una velocidad de 0.5 metros hacia la derecha.

Para obtener los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio, y para hallar las posiciones relativas de las masas mediante el uso del método de Euler Modificado, el código de MATLAB para ambos casos sería el mismo que en el supuesto 1, cambiando las condiciones iniciales

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=10;%Tiempo final
h=0.05;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método de Euler Modificado
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
K1=M*z(:,i);
z(:,i+1)=z(:,i)+h*(M*(z(:,i)+K1*h/2));
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(z1,t,z3,t);
title('Masas 1 y 3');
subplot(1,2,2);
plot(z1,t,z2,t,z3,t);
title('Masas 1,2 y 3');
hold off
}}

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''SISTEMA DE UNA MASA EN MEDIO VISCOSO CON FUERZA EXTERIOR APLICADA'''==

{{matlab| codigo=
%Datos iniciales
m=1;%Masa
k=4;%Constante del muelle
x1=0.3;x11=0.3;%Datos en t0
x0=[x1;x11]';%Condiciones iniciales
t0=0;%Instante inicial
tN=60;%Instante final
nu=1;%Constante de amortiguamiento
h=0.2;%Tamaño de paso
N=(tN-t0)/h;
%Matriz de coeficientes
M=[-nu/m,-k/m;1,0];
%vector tiempo y matriz posicion
t=linspace(t0,tN,N+1);
x=zeros(2,N+1);
x(:,1)=x0;
%Para conocer el valor de la fuerza en cada instante t
f=[];
for i=1:N+1
f=[f;2*exp(-0.01*t(i))*sin(2*t(i))];
end
%Método del Trapecio
I=eye(size(M));
for i=1:N
x(:,i+1)=inv(I-h/2*M)*(x(:,i)+h/2*(M*x(:,i)+f(i)+f(i+1)));
end
subplot(2,1,1);
plot(t,x(2,:));%Tiempo respecto a velocidad
ylabel('Tiempo');xlabel('Velocidad');
E=k/2*(x(1,:).^2)+m/2*(x(2,:).^2);%Ecuación de la energía
subplot(2,1,2);
semilogx(t,E);
title('Cantidad de energía');
}}

La gráfica obtenida es la siguiente:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Modelo de vibraciones con muelles masas (Grupo A8)

2016-05-02T18:43:28Z

Danieldelpotro: /* RESOLUCIÓN DEL SISTEMA DE ECUACIONES */

{{ TrabajoED | Modelo de vibraciones con muelles masas. Grupo 8-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Daniel del Potro Gabin 1423
Humberto del Castillo 1281
Luis de la Fuente Puig 939
Ignacio del Río 1693
}}
<big>
== '''Enunciado y modelización del sistema''' ==
</big>
Muchos modelos de vibraciones se modelizan con sistemas de muelles masas. Consideramos aquí un ejemplo simple de cuatro muelles con constantes kiN/m; i= 1;2;3;4, (fuerzas restauradoras), y tres masas mi; i= 1;2;3, unidos a dos paredes que se deslizan libremente sobre un plano horizontal.

Se modeliza el sistema de ecuaciones diferenciales para el desplazamiento de 3 masas desde la posición de equilibrio. Los grados de libertad que se utilizan son x1, x2 y x3, siendo estos la distancia de cada masa a su posición de equilibrio (Se toma el desplazamiento positivo hacia la derecha). 
Se desprecia la fuerza de rozamiento, por lo que se supone que deslizan sobre una superficie horizontal y lisa, y por tanto las únicas fuerzas a tener en cuenta son las ejercidas por los muelles. Según la ley de Hooke, la fuerza producida por un muelle es proporcional a su elongación. 
Por tanto las fuerzas aplicadas a la primera masa serán: :
<math> F_1=-k_1x_1 </math> :
<math> F_2=k_2(x_2-x_1) </math>
Sobre la segunda masa: :
<math> F_3=-k_2(x_2-x_1) </math> :
<math> F_4=-k_3x_2 </math>
Y las fuerzas que actúan en la tercera masa: :
<math> F_5=-k_3(x_3-x_2) </math> :
<math> F_6=-k_4x_3 </math>

Donde k1, k2, k3 y k4 son las constantes elásticas del movimiento de cada muelle.

Aplicando la segunda Ley de Newton, <math> F=mx'' </math>, se obtiene el siguiente sistema de ecuaciones diferenciales:
<math> m_1x_1''=-k_1x_1+k_2(x_2-x_1)</math>:
<math> m_2x_2''=-k_2(x_2-x_1)-k_3(x_2-x_3)</math>:
<math> m_3x_3''=-k_3(x_3-x_2)-k_4x_3</math>

=='''RESOLUCIÓN DEL SISTEMA DE ECUACIONES'''==

En el instante t=0, las tres masas están desplazadas 0.5,1 y 0.8 metros hacia la derecha respecto de la posición de equilibrio y se sueltan repentinamente,sin velocidad.Se resuelve el sistema de ecuaciones diferenciales utilizando dos métodos numéricos: 
* Método de Euler modificado. 
* Método de Runge-Kutta. 

Primero planteamos el problema en Matlab:

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=0.8;x11=0;x22=0;x33=0;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z=[z1;z2];
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=input('Introduce instante final:');
h=input('Introduce tamaño de paso:');%Probar para h=0.1 y h=0.025
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%f=M*z son las solucines del sistema
%Selecciono la segunda mitad del vector solución, ya que contiene las
%posiciones y la parte superior las derivadas
}}

Y lo resolvemos con ambos métodos:

{{matlab| codigo=
%Método Euler modificado
zEM=z;
for i=1:N
K1=M*zEM(:,i);
zEM(:,i+1)=zEM(:,i)+h*(M*zEM(:,i)+K1*h/2);
end
zEM=zEM([4 5 6],:);%Selecciono las soluciones de posición
zEM1=zEM(1,:)+2.5;%Posición masa 1
zEM2=zEM(2,:)+4;%Posición masa 2
zEM3=zEM(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(zEM1,t,zEM2,t,zEM3,t);
title('Método Euler Modificado');
hold off

%Método Runge-Kutta
zRK=z;
for j=1:N
K1=M*zRK(:,j);
K2=M*(zRK(:,j)+(1/2)*K1*h);
K3=M*(zRK(:,j)+(1/2)*K2*h);
K4=M*(zRK(:,j)+K3*h);
zRK(:,j+1)=zRK(:,j)+h/6*(K1+2*K2+2*K3+K4);
end
zRK=zRK([4 5 6],:);%Selecciono las soluciones de posición
zRK1=zRK(1,:)+2.5;%Posición masa 1
zRK2=zRK(2,:)+4;%Posición masa 2
zRK3=zRK(3,:)+8;%Posición masa 3
ylabel('Tiempo');xlabel('Posición');
hold on
subplot(1,2,2)
plot(zRK1,t,zRK2,t,zRK3,t);
title('Método Runge-Kutta orden 4');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo para Euler-Modificado y Runge-Kutta con h=0.1:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Runge-Kutta]]

Y lo comparamos con la gráfica obtenida con Euler-Modificado y Runge-Kutta con h=0.025:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''EFECTO DE LA VELOCIDAD INICIAL EN LA TRAYECTORIA'''==

Del mismo modo que en el caso anterior, las tres masas parten de una posición diferente a la del equilibrio, pero en este caso se impone una velocidad inicial. Se analiza la trayectoria de las masas en dos supuestos distintos. 

===Supuesto 1===

La primera masa se desplaza a su izquierda 0.5 metros, la segunda 1 metro a su derecha, y la tercera 0.5 metros hacia la derecha, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha.. 

Se obtienen los desplazamientos utilizando el método del trapecio.
{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=-0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=20;%Tiempo final
h=0.2;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método del Trapecio
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
z(:,i+1)=inv(I-h/2*M)*(z(:,i)+h/2*(M*z(:,i)));%Despejando
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
plot(z1,t,z2,t,z3,t);
title('Método del trapecio');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

===Supuesto 2===

La primera masa se desplaza a su derecha 0.5 metros con una velocidad de 1 m/s en dirección hacia la izquierda, la segunda masa se suelta sin velocidad inicial con un desplazamiento de 1 metro a la derecha, y la tercera masa se desplaza a su izquierda 0.5 metros y se suelta con una velocidad de 0.5 metros hacia la derecha.

Para obtener los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio, y para hallar las posiciones relativas de las masas mediante el uso del método de Euler Modificado, el código de MATLAB para ambos casos sería el mismo que en el supuesto 1, cambiando las condiciones iniciales

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=10;%Tiempo final
h=0.05;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método de Euler Modificado
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
K1=M*z(:,i);
z(:,i+1)=z(:,i)+h*(M*(z(:,i)+K1*h/2));
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(z1,t,z3,t);
title('Masas 1 y 3');
subplot(1,2,2);
plot(z1,t,z2,t,z3,t);
title('Masas 1,2 y 3');
hold off
}}

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''SISTEMA DE UNA MASA EN MEDIO VISCOSO CON FUERZA EXTERIOR APLICADA'''==

{{matlab| codigo=
%Datos iniciales
m=1;%Masa
k=4;%Constante del muelle
x1=0.3;x11=0.3;%Datos en t0
x0=[x1;x11]';%Condiciones iniciales
t0=0;%Instante inicial
tN=60;%Instante final
nu=1;%Constante de amortiguamiento
h=0.2;%Tamaño de paso
N=(tN-t0)/h;
%Matriz de coeficientes
M=[-nu/m,-k/m;1,0];
%vector tiempo y matriz posicion
t=linspace(t0,tN,N+1);
x=zeros(2,N+1);
x(:,1)=x0;
%Para conocer el valor de la fuerza en cada instante t
f=[];
for i=1:N+1
f=[f;2*exp(-0.01*t(i))*sin(2*t(i))];
end
%Método del Trapecio
I=eye(size(M));
for i=1:N
x(:,i+1)=inv(I-h/2*M)*(x(:,i)+h/2*(M*x(:,i)+f(i)+f(i+1)));
end
subplot(2,1,1);
plot(t,x(2,:));%Tiempo respecto a velocidad
ylabel('Tiempo');xlabel('Velocidad');
E=k/2*(x(1,:).^2)+m/2*(x(2,:).^2);%Ecuación de la energía
subplot(2,1,2);
semilogx(t,E);
title('Cantidad de energía');
}}

La gráfica obtenida es la siguiente:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Modelo de vibraciones con muelles masas (Grupo A8)

2016-05-02T18:41:40Z

Danieldelpotro: /* Supuesto 2 */

{{ TrabajoED | Modelo de vibraciones con muelles masas. Grupo 8-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Daniel del Potro Gabin 1423
Humberto del Castillo 1281
Luis de la Fuente Puig 939
Ignacio del Río 1693
}}
<big>
== '''Enunciado y modelización del sistema''' ==
</big>
Muchos modelos de vibraciones se modelizan con sistemas de muelles masas. Consideramos aquí un ejemplo simple de cuatro muelles con constantes kiN/m; i= 1;2;3;4, (fuerzas restauradoras), y tres masas mi; i= 1;2;3, unidos a dos paredes que se deslizan libremente sobre un plano horizontal.

Se modeliza el sistema de ecuaciones diferenciales para el desplazamiento de 3 masas desde la posición de equilibrio. Los grados de libertad que se utilizan son x1, x2 y x3, siendo estos la distancia de cada masa a su posición de equilibrio (Se toma el desplazamiento positivo hacia la derecha). 
Se desprecia la fuerza de rozamiento, por lo que se supone que deslizan sobre una superficie horizontal y lisa, y por tanto las únicas fuerzas a tener en cuenta son las ejercidas por los muelles. Según la ley de Hooke, la fuerza producida por un muelle es proporcional a su elongación. 
Por tanto las fuerzas aplicadas a la primera masa serán: :
<math> F_1=-k_1x_1 </math> :
<math> F_2=k_2(x_2-x_1) </math>
Sobre la segunda masa: :
<math> F_3=-k_2(x_2-x_1) </math> :
<math> F_4=-k_3x_2 </math>
Y las fuerzas que actúan en la tercera masa: :
<math> F_5=-k_3(x_3-x_2) </math> :
<math> F_6=-k_4x_3 </math>

Donde k1, k2, k3 y k4 son las constantes elásticas del movimiento de cada muelle.

Aplicando la segunda Ley de Newton, <math> F=mx'' </math>, se obtiene el siguiente sistema de ecuaciones diferenciales:
<math> m_1x_1''=-k_1x_1+k_2(x_2-x_1)</math>:
<math> m_2x_2''=-k_2(x_2-x_1)-k_3(x_2-x_3)</math>:
<math> m_3x_3''=-k_3(x_3-x_2)-k_4x_3</math>

=='''RESOLUCIÓN DEL SISTEMA DE ECUACIONES'''==

En el instante t=0, las tres masas están desplazadas 0.5,1 y 0.8 metros hacia la derecha respecto de la posición de equilibrio y se sueltan repentinamente,sin velocidad.Se resuelve el sistema de ecuaciones diferenciales utilizando dos métodos numéricos: 
* Método de Euler modificado. 
* Método de Runge-Kutta. 

Primero planteamos el problema en Matlab:

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=0.8;x11=0;x22=0;x33=0;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z=[z1;z2];
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=input('Introduce instante final:');
h=input('Introduce tamaño de paso:');%Probar para h=0.1 y h=0.025
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%f=M*z son las solucines del sistema
%Selecciono la segunda mitad del vector solución, ya que contiene las
%posiciones y la parte superior las derivadas
}}

Y lo resolvemos con ambos métodos:

{{matlab| codigo=
%Método Euler modificado
zEM=z;
for i=1:N
K1=M*zEM(:,i);
zEM(:,i+1)=zEM(:,i)+h*(M*zEM(:,i)+K1*h/2);
end
zEM=zEM([4 5 6],:);%Selecciono las soluciones de posición
zEM1=zEM(1,:)+2.5;%Posición masa 1
zEM2=zEM(2,:)+4;%Posición masa 2
zEM3=zEM(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(zEM1,t,zEM2,t,zEM3,t);
title('Método Euler Modificado');
hold off

%Método Runge-Kutta
zRK=z;
for j=1:N
K1=M*zRK(:,j);
K2=M*(zRK(:,j)+(1/2)*K1*h);
K3=M*(zRK(:,j)+(1/2)*K2*h);
K4=M*(zRK(:,j)+K3*h);
zRK(:,j+1)=zRK(:,j)+h/6*(K1+2*K2+2*K3+K4);
end
zRK=zRK([4 5 6],:);%Selecciono las soluciones de posición
zRK1=zRK(1,:)+2.5;%Posición masa 1
zRK2=zRK(2,:)+4;%Posición masa 2
zRK3=zRK(3,:)+8;%Posición masa 3
ylabel('Tiempo');xlabel('Posición');
hold on
subplot(1,2,2)
plot(zRK1,t,zRK2,t,zRK3,t);
title('Método Runge-Kutta orden 4');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo para Euler-Modificado:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Runge-Kutta]]

Y lo comparamos con la gráfica obtenida con Runge-Kutta:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''EFECTO DE LA VELOCIDAD INICIAL EN LA TRAYECTORIA'''==

Del mismo modo que en el caso anterior, las tres masas parten de una posición diferente a la del equilibrio, pero en este caso se impone una velocidad inicial. Se analiza la trayectoria de las masas en dos supuestos distintos. 

===Supuesto 1===

La primera masa se desplaza a su izquierda 0.5 metros, la segunda 1 metro a su derecha, y la tercera 0.5 metros hacia la derecha, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha.. 

Se obtienen los desplazamientos utilizando el método del trapecio.
{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=-0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=20;%Tiempo final
h=0.2;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método del Trapecio
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
z(:,i+1)=inv(I-h/2*M)*(z(:,i)+h/2*(M*z(:,i)));%Despejando
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
plot(z1,t,z2,t,z3,t);
title('Método del trapecio');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

===Supuesto 2===

La primera masa se desplaza a su derecha 0.5 metros con una velocidad de 1 m/s en dirección hacia la izquierda, la segunda masa se suelta sin velocidad inicial con un desplazamiento de 1 metro a la derecha, y la tercera masa se desplaza a su izquierda 0.5 metros y se suelta con una velocidad de 0.5 metros hacia la derecha.

Para obtener los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio, y para hallar las posiciones relativas de las masas mediante el uso del método de Euler Modificado, el código de MATLAB para ambos casos sería el mismo que en el supuesto 1, cambiando las condiciones iniciales

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=10;%Tiempo final
h=0.05;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método de Euler Modificado
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
K1=M*z(:,i);
z(:,i+1)=z(:,i)+h*(M*(z(:,i)+K1*h/2));
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(z1,t,z3,t);
title('Masas 1 y 3');
subplot(1,2,2);
plot(z1,t,z2,t,z3,t);
title('Masas 1,2 y 3');
hold off
}}

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''SISTEMA DE UNA MASA EN MEDIO VISCOSO CON FUERZA EXTERIOR APLICADA'''==

{{matlab| codigo=
%Datos iniciales
m=1;%Masa
k=4;%Constante del muelle
x1=0.3;x11=0.3;%Datos en t0
x0=[x1;x11]';%Condiciones iniciales
t0=0;%Instante inicial
tN=60;%Instante final
nu=1;%Constante de amortiguamiento
h=0.2;%Tamaño de paso
N=(tN-t0)/h;
%Matriz de coeficientes
M=[-nu/m,-k/m;1,0];
%vector tiempo y matriz posicion
t=linspace(t0,tN,N+1);
x=zeros(2,N+1);
x(:,1)=x0;
%Para conocer el valor de la fuerza en cada instante t
f=[];
for i=1:N+1
f=[f;2*exp(-0.01*t(i))*sin(2*t(i))];
end
%Método del Trapecio
I=eye(size(M));
for i=1:N
x(:,i+1)=inv(I-h/2*M)*(x(:,i)+h/2*(M*x(:,i)+f(i)+f(i+1)));
end
subplot(2,1,1);
plot(t,x(2,:));%Tiempo respecto a velocidad
ylabel('Tiempo');xlabel('Velocidad');
E=k/2*(x(1,:).^2)+m/2*(x(2,:).^2);%Ecuación de la energía
subplot(2,1,2);
semilogx(t,E);
title('Cantidad de energía');
}}

La gráfica obtenida es la siguiente:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Modelo de vibraciones con muelles masas (Grupo A8)

2016-05-02T18:35:08Z

Danieldelpotro: /* Supuesto 2 */

{{ TrabajoED | Modelo de vibraciones con muelles masas. Grupo 8-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Daniel del Potro Gabin 1423
Humberto del Castillo 1281
Luis de la Fuente Puig 939
Ignacio del Río 1693
}}
<big>
== '''Enunciado y modelización del sistema''' ==
</big>
Muchos modelos de vibraciones se modelizan con sistemas de muelles masas. Consideramos aquí un ejemplo simple de cuatro muelles con constantes kiN/m; i= 1;2;3;4, (fuerzas restauradoras), y tres masas mi; i= 1;2;3, unidos a dos paredes que se deslizan libremente sobre un plano horizontal.

Se modeliza el sistema de ecuaciones diferenciales para el desplazamiento de 3 masas desde la posición de equilibrio. Los grados de libertad que se utilizan son x1, x2 y x3, siendo estos la distancia de cada masa a su posición de equilibrio (Se toma el desplazamiento positivo hacia la derecha). 
Se desprecia la fuerza de rozamiento, por lo que se supone que deslizan sobre una superficie horizontal y lisa, y por tanto las únicas fuerzas a tener en cuenta son las ejercidas por los muelles. Según la ley de Hooke, la fuerza producida por un muelle es proporcional a su elongación. 
Por tanto las fuerzas aplicadas a la primera masa serán: :
<math> F_1=-k_1x_1 </math> :
<math> F_2=k_2(x_2-x_1) </math>
Sobre la segunda masa: :
<math> F_3=-k_2(x_2-x_1) </math> :
<math> F_4=-k_3x_2 </math>
Y las fuerzas que actúan en la tercera masa: :
<math> F_5=-k_3(x_3-x_2) </math> :
<math> F_6=-k_4x_3 </math>

Donde k1, k2, k3 y k4 son las constantes elásticas del movimiento de cada muelle.

Aplicando la segunda Ley de Newton, <math> F=mx'' </math>, se obtiene el siguiente sistema de ecuaciones diferenciales:
<math> m_1x_1''=-k_1x_1+k_2(x_2-x_1)</math>:
<math> m_2x_2''=-k_2(x_2-x_1)-k_3(x_2-x_3)</math>:
<math> m_3x_3''=-k_3(x_3-x_2)-k_4x_3</math>

=='''RESOLUCIÓN DEL SISTEMA DE ECUACIONES'''==

En el instante t=0, las tres masas están desplazadas 0.5,1 y 0.8 metros hacia la derecha respecto de la posición de equilibrio y se sueltan repentinamente,sin velocidad.Se resuelve el sistema de ecuaciones diferenciales utilizando dos métodos numéricos: 
* Método de Euler modificado. 
* Método de Runge-Kutta. 

Primero planteamos el problema en Matlab:

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=0.8;x11=0;x22=0;x33=0;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z=[z1;z2];
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=input('Introduce instante final:');
h=input('Introduce tamaño de paso:');%Probar para h=0.1 y h=0.025
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%f=M*z son las solucines del sistema
%Selecciono la segunda mitad del vector solución, ya que contiene las
%posiciones y la parte superior las derivadas
}}

Y lo resolvemos con ambos métodos:

{{matlab| codigo=
%Método Euler modificado
zEM=z;
for i=1:N
K1=M*zEM(:,i);
zEM(:,i+1)=zEM(:,i)+h*(M*zEM(:,i)+K1*h/2);
end
zEM=zEM([4 5 6],:);%Selecciono las soluciones de posición
zEM1=zEM(1,:)+2.5;%Posición masa 1
zEM2=zEM(2,:)+4;%Posición masa 2
zEM3=zEM(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(zEM1,t,zEM2,t,zEM3,t);
title('Método Euler Modificado');
hold off

%Método Runge-Kutta
zRK=z;
for j=1:N
K1=M*zRK(:,j);
K2=M*(zRK(:,j)+(1/2)*K1*h);
K3=M*(zRK(:,j)+(1/2)*K2*h);
K4=M*(zRK(:,j)+K3*h);
zRK(:,j+1)=zRK(:,j)+h/6*(K1+2*K2+2*K3+K4);
end
zRK=zRK([4 5 6],:);%Selecciono las soluciones de posición
zRK1=zRK(1,:)+2.5;%Posición masa 1
zRK2=zRK(2,:)+4;%Posición masa 2
zRK3=zRK(3,:)+8;%Posición masa 3
ylabel('Tiempo');xlabel('Posición');
hold on
subplot(1,2,2)
plot(zRK1,t,zRK2,t,zRK3,t);
title('Método Runge-Kutta orden 4');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo para Euler-Modificado:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Runge-Kutta]]

Y lo comparamos con la gráfica obtenida con Runge-Kutta:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''EFECTO DE LA VELOCIDAD INICIAL EN LA TRAYECTORIA'''==

Del mismo modo que en el caso anterior, las tres masas parten de una posición diferente a la del equilibrio, pero en este caso se impone una velocidad inicial. Se analiza la trayectoria de las masas en dos supuestos distintos. 

===Supuesto 1===

La primera masa se desplaza a su izquierda 0.5 metros, la segunda 1 metro a su derecha, y la tercera 0.5 metros hacia la derecha, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha.. 

Se obtienen los desplazamientos utilizando el método del trapecio.
{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=-0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=20;%Tiempo final
h=0.2;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método del Trapecio
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
z(:,i+1)=inv(I-h/2*M)*(z(:,i)+h/2*(M*z(:,i)));%Despejando
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
plot(z1,t,z2,t,z3,t);
title('Método del trapecio');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

===Supuesto 2===

La primera masa se desplaza a su derecha 0.5 metros con una velocidad de 1 m/s en dirección hacia la izquierda, la segunda masa se suelta sin velocidad inicial con un desplazamiento de 1 metro a la derecha, y la tercera masa se desplaza a su izquierda 0.5 metros y se suelta con una velocidad de 0.5 metros hacia la derecha.

Para obtener los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio, y para hallar las posiciones relativas de las masas mediante el uso del método de Euler Modificado, el código de MATLAB para ambos casos sería el mismo que en el supuesto 1, cambiando las condiciones iniciales

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=10;%Tiempo final
h=0.05;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método de Euler Modificado
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
K1=M*z(:,i);
z(:,i+1)=z(:,i)+h*(M*(z(:,i)+K1*h/2));
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(z1,t,z3,t);
title('Masas 1 y 3');
subplot(1,2,2);
plot(z1,t,z2,t,z3,t);
title('Masas 1,2 y 3');
hold off
}}

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Modelo de vibraciones con muelles masas (Grupo A8)

2016-05-02T18:34:06Z

Danieldelpotro: /* Supuesto 2 */

{{ TrabajoED | Modelo de vibraciones con muelles masas. Grupo 8-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Daniel del Potro Gabin 1423
Humberto del Castillo 1281
Luis de la Fuente Puig 939
Ignacio del Río 1693
}}
<big>
== '''Enunciado y modelización del sistema''' ==
</big>
Muchos modelos de vibraciones se modelizan con sistemas de muelles masas. Consideramos aquí un ejemplo simple de cuatro muelles con constantes kiN/m; i= 1;2;3;4, (fuerzas restauradoras), y tres masas mi; i= 1;2;3, unidos a dos paredes que se deslizan libremente sobre un plano horizontal.

Se modeliza el sistema de ecuaciones diferenciales para el desplazamiento de 3 masas desde la posición de equilibrio. Los grados de libertad que se utilizan son x1, x2 y x3, siendo estos la distancia de cada masa a su posición de equilibrio (Se toma el desplazamiento positivo hacia la derecha). 
Se desprecia la fuerza de rozamiento, por lo que se supone que deslizan sobre una superficie horizontal y lisa, y por tanto las únicas fuerzas a tener en cuenta son las ejercidas por los muelles. Según la ley de Hooke, la fuerza producida por un muelle es proporcional a su elongación. 
Por tanto las fuerzas aplicadas a la primera masa serán: :
<math> F_1=-k_1x_1 </math> :
<math> F_2=k_2(x_2-x_1) </math>
Sobre la segunda masa: :
<math> F_3=-k_2(x_2-x_1) </math> :
<math> F_4=-k_3x_2 </math>
Y las fuerzas que actúan en la tercera masa: :
<math> F_5=-k_3(x_3-x_2) </math> :
<math> F_6=-k_4x_3 </math>

Donde k1, k2, k3 y k4 son las constantes elásticas del movimiento de cada muelle.

Aplicando la segunda Ley de Newton, <math> F=mx'' </math>, se obtiene el siguiente sistema de ecuaciones diferenciales:
<math> m_1x_1''=-k_1x_1+k_2(x_2-x_1)</math>:
<math> m_2x_2''=-k_2(x_2-x_1)-k_3(x_2-x_3)</math>:
<math> m_3x_3''=-k_3(x_3-x_2)-k_4x_3</math>

=='''RESOLUCIÓN DEL SISTEMA DE ECUACIONES'''==

En el instante t=0, las tres masas están desplazadas 0.5,1 y 0.8 metros hacia la derecha respecto de la posición de equilibrio y se sueltan repentinamente,sin velocidad.Se resuelve el sistema de ecuaciones diferenciales utilizando dos métodos numéricos: 
* Método de Euler modificado. 
* Método de Runge-Kutta. 

Primero planteamos el problema en Matlab:

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=0.8;x11=0;x22=0;x33=0;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z=[z1;z2];
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=input('Introduce instante final:');
h=input('Introduce tamaño de paso:');%Probar para h=0.1 y h=0.025
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%f=M*z son las solucines del sistema
%Selecciono la segunda mitad del vector solución, ya que contiene las
%posiciones y la parte superior las derivadas
}}

Y lo resolvemos con ambos métodos:

{{matlab| codigo=
%Método Euler modificado
zEM=z;
for i=1:N
K1=M*zEM(:,i);
zEM(:,i+1)=zEM(:,i)+h*(M*zEM(:,i)+K1*h/2);
end
zEM=zEM([4 5 6],:);%Selecciono las soluciones de posición
zEM1=zEM(1,:)+2.5;%Posición masa 1
zEM2=zEM(2,:)+4;%Posición masa 2
zEM3=zEM(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(zEM1,t,zEM2,t,zEM3,t);
title('Método Euler Modificado');
hold off

%Método Runge-Kutta
zRK=z;
for j=1:N
K1=M*zRK(:,j);
K2=M*(zRK(:,j)+(1/2)*K1*h);
K3=M*(zRK(:,j)+(1/2)*K2*h);
K4=M*(zRK(:,j)+K3*h);
zRK(:,j+1)=zRK(:,j)+h/6*(K1+2*K2+2*K3+K4);
end
zRK=zRK([4 5 6],:);%Selecciono las soluciones de posición
zRK1=zRK(1,:)+2.5;%Posición masa 1
zRK2=zRK(2,:)+4;%Posición masa 2
zRK3=zRK(3,:)+8;%Posición masa 3
ylabel('Tiempo');xlabel('Posición');
hold on
subplot(1,2,2)
plot(zRK1,t,zRK2,t,zRK3,t);
title('Método Runge-Kutta orden 4');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo para Euler-Modificado:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Runge-Kutta]]

Y lo comparamos con la gráfica obtenida con Runge-Kutta:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''EFECTO DE LA VELOCIDAD INICIAL EN LA TRAYECTORIA'''==

Del mismo modo que en el caso anterior, las tres masas parten de una posición diferente a la del equilibrio, pero en este caso se impone una velocidad inicial. Se analiza la trayectoria de las masas en dos supuestos distintos. 

===Supuesto 1===

La primera masa se desplaza a su izquierda 0.5 metros, la segunda 1 metro a su derecha, y la tercera 0.5 metros hacia la derecha, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha.. 

Se obtienen los desplazamientos utilizando el método del trapecio.
{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=-0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=20;%Tiempo final
h=0.2;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método del Trapecio
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
z(:,i+1)=inv(I-h/2*M)*(z(:,i)+h/2*(M*z(:,i)));%Despejando
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
plot(z1,t,z2,t,z3,t);
title('Método del trapecio');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

===Supuesto 2===

La primera masa se desplaza a su derecha 0.5 metros con una velocidad de 1 m/s en dirección hacia la izquierda, la segunda masa se suelta sin velocidad inicial con un desplazamiento de 1 metro a la derecha, y la tercera masa se desplaza a su izquierda 0.5 metros y se suelta con una velocidad de 0.5 metros hacia la derecha.

Para obtener los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio, y para hallar las posiciones relativas de las masas mediante el uso del método de Euler Modificado, el código de MATLAB para ambos casos sería el mismo que en el supuesto 1, cambiando las condiciones iniciales

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=10;%Tiempo final
h=0.05;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método de Euler Modificado
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
K1=M*z(:,i);
z(:,i+1)=z(:,i)+h*(M*(z(:,i)+K1*h/2));
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(z1,t,z3,t);
title('Masas 1 y 3');
subplot(1,2,2);
plot(z1,t,z2,t,z3,t);
title('Masas 1,2 y 3');
hold off
}}

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2015-16]]

Modelo de vibraciones con muelles masas (Grupo A8)

2016-05-02T18:31:17Z

Danieldelpotro: Página creada con «{{ TrabajoED | Modelo de vibraciones con muelles masas. Grupo 8-A | Ecuaciones Diferenciales|Curso 2015-16...»

{{ TrabajoED | Modelo de vibraciones con muelles masas. Grupo 8-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Daniel del Potro Gabin 1423
Humberto del Castillo 1281
Luis de la Fuente Puig 939
Ignacio del Río 1693
}}
<big>
== '''Enunciado y modelización del sistema''' ==
</big>
Muchos modelos de vibraciones se modelizan con sistemas de muelles masas. Consideramos aquí un ejemplo simple de cuatro muelles con constantes kiN/m; i= 1;2;3;4, (fuerzas restauradoras), y tres masas mi; i= 1;2;3, unidos a dos paredes que se deslizan libremente sobre un plano horizontal.

Se modeliza el sistema de ecuaciones diferenciales para el desplazamiento de 3 masas desde la posición de equilibrio. Los grados de libertad que se utilizan son x1, x2 y x3, siendo estos la distancia de cada masa a su posición de equilibrio (Se toma el desplazamiento positivo hacia la derecha). 
Se desprecia la fuerza de rozamiento, por lo que se supone que deslizan sobre una superficie horizontal y lisa, y por tanto las únicas fuerzas a tener en cuenta son las ejercidas por los muelles. Según la ley de Hooke, la fuerza producida por un muelle es proporcional a su elongación. 
Por tanto las fuerzas aplicadas a la primera masa serán: :
<math> F_1=-k_1x_1 </math> :
<math> F_2=k_2(x_2-x_1) </math>
Sobre la segunda masa: :
<math> F_3=-k_2(x_2-x_1) </math> :
<math> F_4=-k_3x_2 </math>
Y las fuerzas que actúan en la tercera masa: :
<math> F_5=-k_3(x_3-x_2) </math> :
<math> F_6=-k_4x_3 </math>

Donde k1, k2, k3 y k4 son las constantes elásticas del movimiento de cada muelle.

Aplicando la segunda Ley de Newton, <math> F=mx'' </math>, se obtiene el siguiente sistema de ecuaciones diferenciales:
<math> m_1x_1''=-k_1x_1+k_2(x_2-x_1)</math>:
<math> m_2x_2''=-k_2(x_2-x_1)-k_3(x_2-x_3)</math>:
<math> m_3x_3''=-k_3(x_3-x_2)-k_4x_3</math>

=='''RESOLUCIÓN DEL SISTEMA DE ECUACIONES'''==

En el instante t=0, las tres masas están desplazadas 0.5,1 y 0.8 metros hacia la derecha respecto de la posición de equilibrio y se sueltan repentinamente,sin velocidad.Se resuelve el sistema de ecuaciones diferenciales utilizando dos métodos numéricos: 
* Método de Euler modificado. 
* Método de Runge-Kutta. 

Primero planteamos el problema en Matlab:

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=0.8;x11=0;x22=0;x33=0;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z=[z1;z2];
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=input('Introduce instante final:');
h=input('Introduce tamaño de paso:');%Probar para h=0.1 y h=0.025
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%f=M*z son las solucines del sistema
%Selecciono la segunda mitad del vector solución, ya que contiene las
%posiciones y la parte superior las derivadas
}}

Y lo resolvemos con ambos métodos:

{{matlab| codigo=
%Método Euler modificado
zEM=z;
for i=1:N
K1=M*zEM(:,i);
zEM(:,i+1)=zEM(:,i)+h*(M*zEM(:,i)+K1*h/2);
end
zEM=zEM([4 5 6],:);%Selecciono las soluciones de posición
zEM1=zEM(1,:)+2.5;%Posición masa 1
zEM2=zEM(2,:)+4;%Posición masa 2
zEM3=zEM(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(zEM1,t,zEM2,t,zEM3,t);
title('Método Euler Modificado');
hold off

%Método Runge-Kutta
zRK=z;
for j=1:N
K1=M*zRK(:,j);
K2=M*(zRK(:,j)+(1/2)*K1*h);
K3=M*(zRK(:,j)+(1/2)*K2*h);
K4=M*(zRK(:,j)+K3*h);
zRK(:,j+1)=zRK(:,j)+h/6*(K1+2*K2+2*K3+K4);
end
zRK=zRK([4 5 6],:);%Selecciono las soluciones de posición
zRK1=zRK(1,:)+2.5;%Posición masa 1
zRK2=zRK(2,:)+4;%Posición masa 2
zRK3=zRK(3,:)+8;%Posición masa 3
ylabel('Tiempo');xlabel('Posición');
hold on
subplot(1,2,2)
plot(zRK1,t,zRK2,t,zRK3,t);
title('Método Runge-Kutta orden 4');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo para Euler-Modificado:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Runge-Kutta]]

Y lo comparamos con la gráfica obtenida con Runge-Kutta:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

=='''EFECTO DE LA VELOCIDAD INICIAL EN LA TRAYECTORIA'''==

Del mismo modo que en el caso anterior, las tres masas parten de una posición diferente a la del equilibrio, pero en este caso se impone una velocidad inicial. Se analiza la trayectoria de las masas en dos supuestos distintos. 

===Supuesto 1===

La primera masa se desplaza a su izquierda 0.5 metros, la segunda 1 metro a su derecha, y la tercera 0.5 metros hacia la derecha, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha.. 

Se obtienen los desplazamientos utilizando el método del trapecio.
{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=-0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=20;%Tiempo final
h=0.2;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método del Trapecio
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
z(:,i+1)=inv(I-h/2*M)*(z(:,i)+h/2*(M*z(:,i)));%Despejando
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
plot(z1,t,z2,t,z3,t);
title('Método del trapecio');
hold off
}}

Se obtiene la siguiente gráfica posición-tiempo:

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

===Supuesto 2===

La primera masa se desplaza a su derecha 0.5 metros con una velocidad de 1 m/s en dirección hacia la izquierda, la segunda masa se suelta sin velocidad inicial con un desplazamiento de 1 metro a la derecha, y la tercera masa se desplaza a su izquierda 0.5 metros y se suelta con una velocidad de 0.5 metros hacia la derecha.

Para obtener los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio, y para hallar las posiciones relativas de las masas mediante el uso del método de Euler Modificado, el código de MATLAB para ambos casos sería el mismo que en el supuesto 1, cambiando las condiciones iniciales

{{matlab| codigo=
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=10;%Tiempo final
h=0.05;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método de Euler Modificado
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
K1=M*z(:,i);
z(:,i+1)=z(:,i)+h*(M*(z(:,i)+K1*h/2));
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(z1,t,z3,t);
title('Masas 1 y 3');
subplot(1,2,2);
plot(z1,t,z2,t,z3,t);
title('Masas 1,2 y 3');
hold off
}}

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por el método del trapecio]]

[[Archivo:INSERTAR-IMAGEN.jpg|marco|centro|Aproximación numérica del sistema de ecuaciones por Euler Modificado]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Comportamiento de un fluido incompresible alrededor de un obstáculo circular

2015-12-03T22:16:20Z

Danieldelpotro:

{{ TrabajoED | Comportamiento de un fluido incompresible. Grupo 4-C | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC15/16|2015-16]] | Adrián Coronado López 
Humberto del Castillo Montes de Oca 
Carlos Nieto Egido 
Jorge García Flores 
Daniel del Potro Gabín 
Gonzalo Fernández Arroyo}} 

En este articulo nos centraremos en analizar el comportamiento de un fluido incompresible al interponerse en su camino un obstáculo.

== Visualización del recinto==
Inicialmente consideramos una región con forma de corona circular de radio interior dos y radio exterior seis del plano Z=0 ocupada por un fluido.
El obstáculo esta constituido por la circunferencia interior.
Para visualizar el campo vectorial creamos un mallado plano que ocupe la región de estudio del fluido, que sería, en los ejes cartesianos [-5,5]x[-5,5].
[[Archivo:Región ocupada oor el fluido.jpg|1000px|thumb|right|Región ocupada por el fluido]]
{{matlab|codigo=
clear
%Declarar paso y variables
h=0.1;
ro=2:h:6;
th=0:h:2*pi+h;
%generamos la retícula
[U,V]=meshgrid(ro,th);

%hacemos un cambio de coordenadas Polares-Cartesianas
X=U.*cos(V);
Y=U.*sin(V);

%dibujo del mallado
mesh(X,Y,0*X)
axis ([-5,5,-5,5])
axis equal
view(2) %para verlo en 2D
}}

==Campo de velocidad del fluido==
Ahora veremos como ese fluido se mueve en nuestro mallado, para ello veremos la variación de la velocidad de las partículas del fluido.
Esta velocidad viene dada por el gradiente de la función potencial: \[\varphi=2cos\theta(\rho+\frac{4}{\rho})\] , la cual nos indica la dirección más probable que toma el fluido.
 \[\vec{u}=\nabla\varphi=\frac{\partial\varphi}{\partial\rho}\vec{g^\rho}+\frac{\partial\varphi}{\partial\theta}\vec{g^\theta}+\frac{\partial\varphi}{\partial z}\vec{g^z}\]
\[\frac{\partial\varphi}{\partial z}\vec{g^z}=0\]
\[\vec{u}=2cos\theta(1-\frac{4}{\rho^2})\vec{g^\rho}-2sin\theta(\rho+\frac{4}{\rho})\vec{g^\theta}\]
\[\vec{g^\rho}=\vec{g_\rho}; \vec{g^\theta}=\frac{1}{\rho^2}\vec{g_\theta}\]
\[\vec{u}=2cos\theta(1-\frac{4}{\rho^2})\vec{g_\rho}-2sin\theta(\frac{1}{\rho}+\frac{4}{\rho^3})\vec{g_\theta}\]

Si '''\(\vec{n}\)''' es el vector normal a los puntos del obstáculo, <math>\vec{u}</math>•\(\vec{n}\)=0, lo que significa que \(\vec{n}\) es ortogonal a <math>\vec{u}</math>

===Representación gráfica===
El programa para representar la función potencial del fluído:

[[Archivo:Grafico2 C19.jpg|400px|miniaturadeimagen|derecha|Función potencial <math>\varphi</math> del fluido]]{{matlab|codigo=
%Declaramos el paso y las variables
h=0.1;
rho=2:h:6;
th=0:h:2*pi+h;

%Generamos la retícula
X=U.*cos(V);
Y=U.*sin(V);

%Función potencial
f=2*(U+(4./U)).*cos(V);

%Representación de la función
surf(X,Y,f)
axis([-5,5,-5,5])
axis equal
view(2)
}}

El programa para representar el campo de velocidades <math>\vec u</math> del fluido:

[[Archivo:Grafico3 C19.jpg|400px|miniaturadeimagen|derecha|Campo de velocidades <math>\vec u</math> del fluido sobre las curvas de nivel de <math>\varphi</math>]]{{matlab|codigo=
%Declarar paso y variables
h=0.3; %tomamos un paso diferente para que los vectores se vean más claramente
rho=2:h:6;
th=0:h:2*pi+h;

%Generamos la retícula
[U,V]=meshgrid(rho,th);

X=U.*cos(V);
Y=U.*sin(V);

%Funciones
fx=((2*(1./U-4./U.^3).*cos(V)).*X)+((2*(1./U+4./U.^3).*sin(V)).*Y); %componente en i
fy=((2*(1./U-4./U.^3).*cos(V)).*Y)-((2*(1./U+4./U.^3).*sin(V)).*X); %componente en j

%Representación de las curvas de nivel
contour(X,Y,f,50,'LineWidth',1.5) %utilizamos un grosor mayor para poder luego hacer
axis([-5,5,-5,5]) %una observación más precisa
axis equal
view(2)
hold on

%Representación del campo de velocidades
quiver(X,Y,fx,fy)
axis([-5,5,-5,5])
axis equal
view(2)
}}

Como podemos observar en el gráfico las lineas de nivel de la función potencial son ortogonales a los vectores, velocidad del fluido,<math>\vec{u}</math>. Al existir un campo vectorial que es igual al gradiente de un campo escalar, ese campo vectorial es un campo conservativo. Como veremos más adelante la circulación es nula.
[[Archivo:Zoom.jpg|600px|miniaturadeimagen|centro|Detalle del campo de velocidades del fluido perpendicular a la función potencial <math>\varphi</math>]]

==Suponemos <math>\rho</math> lo suficientemente grande como para despreciar <math>\frac{1}{\rho}</math>==

Calculamos aproximadamente el campo de velocidades:

<math>\displaystyle\lim_{\rho \to {\infty}}{\displaystyle(2cos\theta(1-\frac{4}{\rho^2})cos\theta+2sin\theta(\frac{1}{\rho}+\frac{4}{\rho^3}\rho sin\theta))\vec{i}+(2cos\theta sin\theta(1-\frac{4}{\rho^2})-2sin\theta\rho cos\theta(\frac{1}{\rho}+\frac{4}{\rho^3}))\vec{j}}=(2(cos^2\theta)+2(sin^2\theta))\vec{i}</math>

==Rotacional y Divergencia==

===Rotacional===
En el cálculo vectorial, el rotacional es un operador vectorial que muestra la tendencia de un campo vectorial a inducir rotación alrededor de un punto.
<math>\nabla \times\bar { u } =\frac { 1 }{ \sqrt { g } } \begin{vmatrix} \bar { { g }_{ \rho } } & \bar { { g }_{ \theta } } & \bar { { g }_{ z } } \\ \frac { \partial }{ \partial \rho } & \frac { \partial }{ \partial \theta } & \frac { \partial }{ \partial z } \\ { u }_{ \rho } & { u }_{ \theta } & { u }_{ z } \end{vmatrix}=\frac { 1 }{ \rho } (\bar { { g }_{ \rho } } (\frac { \partial { u }_{ z } }{ \partial \theta } -\frac { \partial { u }_{ \theta } }{ \partial z } )+\bar { { g }_{ \theta } } (\frac { \partial { u }_{ \rho } }{ \partial z } -\frac { \partial { u }_{ z } }{ \partial \rho } )+\bar { { g }_{ z } } (\frac { \partial { u }_{ \theta } }{ \partial \rho } -\frac { \partial { u }_{ \rho } }{ \partial \theta } ))=\frac { 1 }{ \rho } (\frac { \partial { u }_{ \theta } }{ \partial \rho } -\frac { \partial { u }_{ \rho } }{ \partial \theta } )\bar { { g }_{ z } }</math> 

<math>{ u }_{ z }=0</math> 

<math>{ u }_{ \theta }=\frac { \partial \varphi }{ \partial \theta } =-2(\rho +\frac { 4 }{ \rho } )sen(\theta )\qquad \qquad \frac { \partial { u }_{ \theta } }{ \partial \rho } =-2sen(\theta )(1-\frac { 4 }{ { \rho }^{ 2 } } )</math> 

<math>{ u }_{ \rho }=\frac { \partial \varphi }{ \partial \rho } =2cos(\theta )(1-\frac { 4 }{ { \rho }^{ 2 } } ) \qquad \qquad \frac { \partial { u }_{ \rho } }{ \partial \theta } =-2sen(\theta )(1-\frac { 4 }{ { \rho }^{ 2 } } )</math> 

<math>\boxed { \nabla \times (\bar { u } )=0 } </math>

===Divergencia===

Un fluido incompresible permanece constante con el tiempo y tiene la capacidad de oponerse a la compresión del mismo bajo cualquier condición. Esto quiere decir que ni la masa ni el volumen del fluido puede cambiar. 
El agua es un fluido incompresible, es decir, la cantidad de volumen y de masa permanecerán prácticamente iguales, aún bajo presión. 
Para comprobar la incompresibilidad debemos comprobar que la divergencia sea nula.
Calculando para nuestra situación estudiada::

<math>div(\bar { u } )=\frac { 1 }{ \sqrt { g } } (\frac { \partial }{ \partial { x }^{ i } } (\sqrt { g } { u }^{ i }))=\frac { 1 }{ \rho } \left( \frac { \partial (\rho (2cos(\theta )(1-\frac { 4 }{ { \rho }^{ 2 } } ))) }{ \partial \rho } +\frac { \partial (\rho (-2\left( \frac { { \rho }^{ 2 }+4 }{ { \rho }^{ 3 } } \right) sen\theta ) }{ \partial \theta } \right) </math> :

<math>\frac { 1 }{ \rho } (2(cos(\theta )(1+\frac { 4 }{ { \rho }^{ 2 } } )-2cos(\theta )(1+\frac { 4 }{ { \rho }^{ 2 } } ))=0</math> :
<math>\boxed { div(\bar { u } )=0 } </math> 

===Conclusiones===
Sabiendo que <math>\vec u </math> es el campo de velocidades de nuestro fluido, hemos calculado que su rotacional es nulo, por lo que las partículas no giran. En general, este rotacional indica la cantidad neta de giro del fluido en dirección contraria a las manecillas del reloj alrededor de nuestra superficie. Indica por tanto su circulación, que como ya hemos indicado es nula.
La divergencia del campo vectorial <math>\vec u </math> en un punto P es el flujo que sale del punto por unidad de volumen. En este caso la divergencia representa la razón de expansión por unidad de volumen bajo el flujo del fluido. Como la divergencia también nos da nula, significa que el fluido no varía, se mantiene constante.
Todo esto nos lleva a afirmar que nuestro campo ni gira ni se expande, concluyendo que se trata de un fluido incompresible.

==Puntos de la frontera==

Vamos a analizar los puntos de la frontera en los valores 2 y 6. Para ello primero sustituiremos los valores de <math>( \rho =2,\quad \rho =6 )</math> en el campo <math>\overline { u }</math> y después derivamos el campo <math>\overline { u }</math> y sacamos los puntos en los que se anula la derivada de <math>\overline { u }</math> para obtener los puntos en los que encontraremos máximos y mínimos. Para saber si los puntos son máximos o mínimos lo evaluaremos en la función, no en la derivada.

Campo <math>\bar { u }</math>::
<math>\overline { u } =2cos(\theta )(1-\frac { 4 }{ \rho ^{ 2 } } ){ \overline { g } }_{ \rho }-2sen(\theta )\left( \frac { { \rho }^{ 2 }+4 }{ { \rho }^{ 3 } } \right) { \overline { g } }_{ \theta }</math> 

Para <math>\rho =2</math> 
<math>\bar { u } =-2sen(\theta )</math> 
<math>{ \bar { u } }'=-2cos(\theta )=0 </math> 

Igualando a 0 obtenemos los valores 
<math>\theta =\frac { \pi }{ 2 } </math>Evaluando en u sabemos que es un Mínimo 
Análogamente se haría lo mismo en los otros valores de <math>\theta</math> para saber si son máximos o mínimos 
<math>\theta =3\frac { \pi }{ 2 } </math> Máximo 

Para el valor <math>\underline { \rho =6 }:</math> 
<math>\bar { u } ={ \left( \frac { 16 }{ 9 } cos\theta \right) }^{ 2 }-{ \left( \frac { 10 }{ 2z } sen\theta \right) }^{ 2 }</math> 

<math>{ \bar { u } }'=2{ \left( \frac { 16 }{ 9 } cos\theta \right) }{ \left( \frac { -16 }{ 9 } sen\theta \right) }-{ 2\left( \frac { 10 }{ 2z } sen\theta \right) }{ \left( \frac { 10 }{ 2z } cos\theta \right) }=0</math> 

Al igualar a 0 obtenemos los siguientes valores: 
<math>\theta =0 </math> maximo 
<math>\theta =\frac { \pi }{ 2 } </math>Mínimo 
<math>\theta = 3\frac { \pi }{ 2 } </math>Mínimo 
<math>\theta = { \pi } </math>Máximo 

A continuación vamos a obtener numericamente los máximos y los mínimos con MATLAB:
{{matlab|codigo=
%Creamos el mallado.
u=linspace(2,6,100);
v=linspace(0,2*pi,100);
[U,V]=meshgrid(u,v);
%Definimos la parametrizacion y la función velocidad en ro=2.
X=U.*cos(V);
Y=U.*sin(V);
fu=(2.*sin(V));
%Calculamos la velocidad maxima.
maximo=max(max(fu))
}}

Vamos a ver los puntos del fluido que tienen módulo de velocidad menor: 
{{matlab|codigo=
%Creamos el mallado.
u=linspace(2,6,100);
v=linspace(0,2*pi,100);
[U,V]=meshgrid(u,v);
%Definimos la parametrizacion y la funcion velocidad en ro=2.
X=U.*cos(V);
Y=U.*sin(V);
fu=(2.*sin(V));
%Calculamos el minimo valor absoluto de esta.
minimo=min(min(abs(fu)))
}}

[[Archivo:LineasDeCorrienteYCurvasDeNivel.jpg|miniaturadeimagen|derecha]]
[[Archivo:CorrientePotencial|derecha]]
[[Archivo:Abd|sinmarco|izquierda|kdfdn]]
[[Archivo:CorrientePotencial.jpg|miniaturadeimagen]]
==Líneas de Corriente==
Las '''líneas de corriente''' del campo \(\vec{u}\) son tangentes a la velocidad, entendida como el gradiente de la función potencial '''φ'''.
Para poder dibujar las líneas de corriente procedemos al cálculo de un vector perpendicular a \(\vec{u}\), el vector <math>\vec{v}=\vec{k}\times\vec{u}</math>

<math>\vec{v}=\vec{k}\times\vec{u}= \sqrt { g }\begin{vmatrix} \bar { { g }^{ \rho } } & \bar { { g }^{ \theta } } & \bar { { g }^{ z } } \\ 0 & 0 & 1 \\ { u }^{ \rho } & { u }^{ \theta } & { u }^{ z } \end{vmatrix}= ρ (\vec{u^ρ} \vec{g^θ} - \vec{u^θ} \vec{g^ρ})</math>

→ <math>\vec{v}={2(\frac{1}{ρ}-\frac{4}{ρ^3})cosθ \vec{g_θ} + 2(1 + \frac{4}{ρ^2}) sen θ \vec{g_ρ}}</math>

<math>\nabla\times\vec{v}=\frac{1}{\sqrt { g }}\begin{vmatrix} \bar { { g }_{ \rho } } & \bar { { g }_{ \theta } } & \bar { { g }_{ z } } \\ \frac { \partial }{ \partial \rho } & \frac { \partial }{ \partial \theta } & \frac { \partial }{ \partial z } \\ { v }_{ \rho } & { v }_{ \theta } & { v }_{ z } \end{vmatrix}=\frac{2}{ρ}\begin{vmatrix} \bar { { g }_{ \rho } } & \bar { { g }_{ \theta } } & \bar { { g }_{ z } } \\ \frac { \partial }{ \partial \rho } & \frac { \partial }{ \partial \theta } & \frac { \partial }{ \partial z } \\ (1+\frac{4}{ρ^2}) senθ & ρ(1-\frac{4}{ρ^2}) cosθ & 0 \end{vmatrix}=\frac{1}{ρ}((2+\frac{8}{ρ^2})cosθ- (2+\frac{8}{ρ^2})cosθ)\vec{g_z}= \frac{2}{ρ}((1+\frac{4}{ρ^2}) - (1+\frac{4}{ρ^2}))cosθ\vec{g_z} = 0</math> 

→ <math> { \bar { rot } (\bar { v } )=\nabla\times\vec{v}=0 } </math>

Además \(\vec{v}\) tendrá un '''potencial escalar ψ''' que se conoce como la función de corriente de \(\vec{u}\). Dicho potencial lo calcularemos gracias al campo \(\vec{v}\) pues es el gradiente de dicho potencial.

<math>\vec{v}=\vec{grad}ψ= \frac{∂ψ}{∂ρ}\vec{g^ρ} + \frac{∂ψ}{∂θ}\vec{g^θ} + \frac{∂ψ}{∂z}\vec{g^z} = 2(1 + \frac{4}{ρ^2})senθ \vec{g^ρ} + 2ρ(1 - \frac{4}{ρ^2})cosθ \vec{g^θ}</math>

<math>\frac{∂ψ}{∂ρ}= 2(1 + \frac{4}{ρ^2})senθ → ψ=\int 2(1 + \frac{4}{ρ^2})senθ dρ= 2(ρ-\frac{4}{ρ})senθ + h(θ) \\ \frac{∂ψ}{∂θ}=2ρ(1-\frac{4}{ρ^2})cosθ=2(ρ-\frac{4}{ρ})cosθ+h'(θ) → h'(θ) = 0 → h(θ) = cte </math> 

→ <math> {ψ=2(ρ-\frac{4}{ρ})senθ}</math>.

Dibujamos las líneas de corriente del campo de velocidades \(\vec{u}\) sobre las que añadiremos las curvas de nivel del potencial escalar ψ.

[[Archivo:CorrienteYVelocidad.png|800x500px|miniaturadeimagen|derecha|Campo de velocidad y Líneas De Corriente ]]

{{matlab|codigo=
clear
%Declara paso y variables
h = 0.3;
ro = 2:h:2*pi+h;
tg= 0:h:2*pi+h;
%Generar retícula

[U,V] = meshgrid(ro,tg);

%Cambio de polares a cartesianas
X = U.*cos(V);
Y = U.*sin(V);

%Función de corriente
f = 2*(U-(4./U)).*sin(V);

%Funciones
fx=((2*(1./U-4./U.^3).*cos(V)).*X)+((2*(1./U+4./U.^3).*sin(V)).*Y);
fy=((2*(1./U-4./U.^3).*cos(V)).*Y)-((2*(1./U+4./U.^3).*sin(V)).*X);

%Dibujar lineas de corriente
contour(X,Y,f,50,'Linewidth',1.3)
axis([-5,5,-5,5])
axis equal
view(2)
hold on

%Dibujar campo de velocidades
quiver(X,Y,fx,fy)
axis([-5,5,-5,5])
axis equal
view(2)
}}

Finalmente, representamos las líneas de corriente y las líneas equipotenciales. Se puede comprobar que son ortogonales.

[[Archivo:CorrienteYEquipotenciales.png|500x500px|miniaturadeimagen|derecha|Líneas de corriente y equipotenciales ]]

{{matlab|codigo=
clear
%Declara paso y variables
h = 0.3;
ro = 2:h:2*pi+h;
tg= 0:h:2*pi+h;

%Generar retícula
[U,V] = meshgrid(ro,tg);

%Cambio de polares a cartesianas
X = U.*cos(V);
Y = U.*sin(V);

%Función de corriente
f = 2*(U-(4./U)).*sin(V);

%Funciones
fx=((2*(1./U-4./U.^3).*cos(V)).*X)+((2*(1./U+4./U.^3).*sin(V)).*Y);
fy=((2*(1./U-4./U.^3).*cos(V)).*Y)-((2*(1./U+4./U.^3).*sin(V)).*X);

%Dibujar lineas de corriente
contour(X,Y,f,50,'Linewidth',1.3)
axis([-5,5,-5,5])
axis equal
view(2)
hold on

%Declara paso y variables
h=0.3;
ro=2:h:6;
tg=0:h:2*pi+h;

%Generar retícula
[U,V]=meshgrid(ro,tg);

%Cambio de polares a cartesianas
X = U.*cos(V);
Y = U.*sin(V);

%Función potencial
f=2*(U+(4./U)).*cos(V);

%Funciones
fx=((2*(1./U-4./U.^3).*cos(V)).*X)+((2*(1./U+4./U.^3).*sin(V)).*Y);
fy=((2*(1./U-4./U.^3).*cos(V)).*Y)-((2*(1./U+4./U.^3).*sin(V)).*X);

%Dibujar curvas de nivel
contour(X,Y,f,50,'Linewidth',1.6)
axis([-5,5,-5,5])
axis equal
view(2)
}}

==Distribución de Presiones en el fluido==
===La ecuación de Bernouilli===
Esta ecuación describe la distribución de presiones en el fluido. \[\frac{1}{2}d\left | \vec{u} \right |^2+p=15\]

Siendo:
d: densidad del fluido incomprensible
<math>\vec{u}</math>: campo vectorial de velocidades del fluido
p: presión estática a la que esta sometido el fluido debido a sus moléculas

Para poder aplicar esta ecuación debemos suponer que:
:1.el fluido(gas o liquido) se mueve en un régimen permanente, es decir, <math> \frac{\partial}{\partial t}=0 </math>
:2.se desprecia la viscosidad (<math>\mu=0</math>)

La ecuación nos indica que a medida que se aumenta la velocidad del flujo, la presión estática en el fluido disminuye. Aplicamos esto a nuestros datos para obtener la ecuación de la presión: \[|\vec{u}|^2+p=15\]
\[|\vec{u}|^2=4(cos\theta)^2(1-\frac{4}{\rho^2})^2+4(sin\theta)^2(\frac{1}{\rho}+\frac{4}{\rho^3})^2\]

\[p=15-4(cos\theta)^2(1-\frac{4}{\rho^2})^2+4(sin\theta)^2(\frac{1}{\rho}+\frac{4}{\rho^3})^2\]

===Representación gráfica===
[[Archivo:Presiones2D.jpg|800px|miniaturadeimagen|derecha|Representación en 2D del comportamiento de las presiones]]

{{matlab|codigo=
clear
%Declarar paso y variables
h=0.1;
ro=2:h:6;
th=0:h:2*pi+h;
%generamos la retícula
[U,V]=meshgrid(ro,th);

%hacemos un cambio de coordenadas Polares-Cartesianas
X=U.*cos(V);
Y=U.*sin(V);
%la presion viene dada por la función
p=15-(4*cos(V).^2.*(1-4./U.^2).^2+4*sin(V).^2.*(1./U+4./U.^3).^2)

%dibujamos la superficie de presiones en el fluido
figure(1)
surf(X,Y,p)
axis[-5,5,-5,5]
axis equal
view(2)
figure(2)
surf(X,Y,p)
axis[-5,5,-5,5]
axis equal
view(3)

Pmax=max(max(p));%presión máxima
Pmin=min(min(p));%presión mínima

%la velocidad será mínima cuando la presión sea máxima y viceversa
}}
[[Archivo:Presiones(3D).jpg|800px|miniaturadeimagen|centro|Representación en 3D del comportamiento de las presiones]]

===Interpretación de una linea de corriente cercana al objeto===
Como podemos observar en las imágenes, a medida nos acercamos al objeto, aumenta la presión y disminuye la velocidad del líquido.
Incrementando dicha velocidad paulatinamente a medida que el objeto es rodeado, lo cual indica que la presión es inversamente proporcional a la velocidad.

===Presiones medias===
Para ello aproximaremos la integral de la presión en todo el fluido y dividirlo por el área total en el anillo 2 < <math>\rho</math> < 6.
:Para el cálculo de estas presiones medias podemos utilizar la aproximación trapezoidal:

{{matlab|codigo=
clear
%lo primero que haremos sera definir el mallado y para ello necesito
%definir las variables

h=0.01;
ro=2:h:6-h;
th=0:h:2*pi-h;

%el mallado

[U,V]=meshgrid(ro,th);

%haremos un cambio de coordenadas
X=U.*cos(V);
Y=U.*sin(V);
p=15-(4*(cos(V).^2).*((1-4./U.^2).^2)+4*(sin(V).^2).*((1./U+4./U.^3).^2));

%realizamos una integral aproximada
P=U.*p;

%resolución de la integral
volumen=h^2.*P;
resultado=sum(sum(volumen));
%area del anillo
area=pi*(6^2-2^2);
Pmedia=resultado/area;
}}
el resultado de la <math>\textbf{Presión media=13.574}</math> que al fin y al cabo no es mas que obrservar la representación gráfica de la distribución de presiones.

==La Paradoja de D'Alambert==

Se trata de una contradicción a la que llegó D'Alembert luego de estudiar matemáticamente el fenómeno de la resistencia producida sobre un cuerpo cuando una corriente de fluido (líquido o gas) circula sobre él.
D'Alembert aplicó la teoría de flujo potencial para modelar el fenómeno, y concluyó que la fuerza resultante sobre el cuerpo sobre el cual fluye el aire es cero, lo cual se contradice con la observación.
Para saber que esa fuerza es nula aplicamos el Teorema de Kutta-Joukowski que sabiendo que la circulación de <math>\overrightarrow {u} </math> es nula y que es proporcional a la fuerza que se ejerce sobre el cuerpo. Por lo tanto, vamos a demostrar que la circulación de <math> \overrightarrow {u} </math> es nula.

Para ello utilizamos el Teorema de Stokes el cual nos dice que: \[\int \vec{u}\vec{t}ds=\int_{S}(\nabla\times \vec{u})\vec{\partial S}\] 

\[\nabla\times \vec{u}=\begin{vmatrix} \vec{g_\rho} & \vec{g_\theta} & \vec{g_z}\\ \frac{\partial}{\partial\rho} & \frac{\partial}{\partial\theta} & \frac{\partial}{\partial z}\\ 2cos\theta(1-\frac{4}{\rho^2}) & -2sin\theta(\rho+\frac{4}{\rho}) & 0 \end{vmatrix}=\frac{1}{\rho}(-2sin\theta(1-\frac{4}{\rho^2})\vec{g_z}+2sin\theta(1-\frac{4}{\rho^2})\vec{g_z})=0\]

Pero como, al fin y al cabo, el objeto si que se mueve si un fluido lo rodea, tendríamos que considerar la fuerza de la viscosidad, es decir, en un fluido perfecto e incompresible como el que tenemos la fuerza de rozamiento debido a la viscosidad es nula. Por ello si tomamos el fluido como un fluido mas sofisticado (<math>\mu>0</math>), esa fuerza sera diferente de 0.

== Ecuación de '''Navier-Stokes''' ==

Las '''ecuaciones de Navier-Stokes''' reciben su nombre de Louis Navier y George Gabriel Stokes. Se trata de un conjunto de ecuaciones en derivadas parciales no lineales que describen el movimiento de un fluido.

\(\vec{u} \cdot \nabla \vec{u} + \nabla p=μ \Delta\vec{u} \)

=== Comprobación de la ecuación ===

A continuación, probamos que nuestro campo de presiones y velocidades satisfacen la ecuación y la condición de incompresibilidad (el agua ocupa siempre el mismo volumen).

[[Archivo:apartado1formu.PNG|500px|centro]]

[[Categoría:Teoría de Campos]]
[[Categoría:TC15/16]]