MateWiki - Contribuciones del usuario [es]

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T19:18:12Z

Jose Manuel Corbí Garrido: /* Con veinte términos de la serie */

{{ Beta }}
{{ TrabajoED | Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]]|

María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

=Introducción y modelización del problema=

En este artículo analizaremos el comportamiento de un cable tenso de 10 metros de longitud en una estructura civil. El cable se encontrará sujeto por los extremos y será sometido a pequeñas vibraciones. Para estudiar el problema haremos ciertas hipótesis que simplificarán el problema, y que en muchos casos dan una buena aproximación del comportamiento real que tendría un cable como el que aquí estudiamos.

En cuanto a las '''características mecánicas del cable''' haremos las siguientes '''hipótesis''':

* El cable tendrá una sección despreciable frente a su longitud.
* Será un cable homogéneo, es decir, la masa por unidad de volumen será constante.
* Dicho cable también será perfectamente flexible, es decir, solo ofrecerá resistencia a esfuerzos tangenciales (en su dirección longitudinal), siendo por ello nula su resistencia frente a esfuerzos de flexión y corte.
* El cable al someterlo a tracción es inextensible, es decir, que es lo suficientemente rígido (en dirección longitudinal) como para poder despreciar su extensibilidad, por el contrario al someterlo a compresión no ofrece resistencia y se arruga.

Por otra parte el cable será sometido a '''pequeñas vibraciones''', las cuales podemos modelizar mediante la '''ecuación de ondas''' <math>u_{tt}-u_{xx}=0</math>. Estas vibraciones, al ser pequeñas, son despreciables frente a su longitud total; preocupándonos únicamente por los desplazamientos verticales (fruto de la tracción) que sufre el cable. Estos desplazamientos verticales serán modelizados mediante la función <math>u(x,t)</math>, dependiente de dos variables; la <math>x</math> que indica la distancia al extremo de la cuerda tomado como referencia, y la <math>t</math> que define el tiempo.

Las condiciones de frontera en este caso nos indican que el cable no sufrirá desplazamiento vertical en sus extremos <math>u(0,t)=0; u(10,t)=0</math>. Por otra parte las condiciones iniciales <math>u(x,0)=0; u_{t}(x,0)=0</math> se traducen en que inicialmente el desplazamiento vertical y la velocidad del cable en todos sus puntos es cero.

El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

En los sucesivos apartados estudiaremos el comportamiento del cable en diferentes situaciones, que se traducirán en diferentes condiciones para el problema (cambios en las condiciones y en la ecuación de ondas con la que modelizamos las pequeñas vibraciones).

=Desplazamiento vertical del cable=
Sujetando el cable desde el centro y desplazándolo dos metros en la dirección perpendicular, estudiaremos el comportamiento del cable una vez lo soltemos. Estará sometido a unas pequeñas vibraciones, que originarán desplazamientos verticales a lo largo del cable que aproximaremos mediante la función <math>u(x,t)</math>. Teniendo en cuenta el desplazamiento vertical de dos metros y que el cable es soltado con una velocidad nula, el problema que vamos a estudiar es el siguiente:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

A continuación estudiaremos el desplazamiento vertical al que se ve cada punto del cable expuesto y para cada instante de tiempo. Para ello aproximaremos la ecuación diferencial que modeliza el comportamiento de dicho cable mediante los métodos de diferencias finitas y mediante series de Fourier, obteniendo gráficas que reflejan la posición de cada punto del cable en cada instante.
==Aproximación mediante diferencias finitas==
Dentro del método de diferencias finitas se nos plantea un sistema de ecuaciones diferenciales, que resolveremos aplicando diferentes métodos de los cuales luego compararemos su precisión a la hora de aproximar el comportamiento del cable y su estabilidad.

El primer método que utilizaremos para resolver el sistema de ecuaciones diferenciales será el '''método del trapecio''', el cual es '''incondicionalmente estable'''. Para implementar este método y observar el comportamiento del cable en una gráfica usaremos el siguiente código en matlab:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')
}}
Obteniendo la gráfica:

[[Archivo:ImagenCodigo1GRP9.jpeg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método del trapecio]]

Usaremos también el '''método de Euler''', el cual es un método explícito y por tanto dependemos de la ecuación en que pueda ser o no factible de ejecutar. El siguiente código matlab nos permite ejecutarlo:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler explícito

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler explícito

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for j=1:length(t)-1
u=u+dt*v+F;
v=v-dt*k*u+F;
sol(j+1,:)=[0 , u' , 0];
end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler explícito')
}}
Obteniendo el siguiente gráfico del comportamiento del cable respecto al tiempo:

[[Archivo:ImagenCodigo2GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler explícito]]

Por último método para resolver el sistema planteado mediante diferencias finitas, resolvemos numéricamente mediante el '''método de Euler modificado''', el cual tiene más precisión que el método de Euler. El código que se muestra a continuación ejecuta dicho método y nos proporciona una gráfica del cable:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler modificado

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler modificado

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for n=1:length(t)-1
l=k*u;
% calculamos k1=f(t_n,y_n)
k1=[ v l]';
% Calculamos k2=f(t_n+h/2,y_n+1/2k1*h)
%xp=x(n)+dt; %argumento de tiempo
%solp=sol+dt*k1; %argumento y
k2=[v l]';
%Calculamos k3=f(t_n+h/2,y_n+1/2k2*h)
%mismo argumento tiempo que antes

u=u+dt/2*(k1(1,:)'+k2(1,:)')+F;
l=+k*u;
k1=[v l]';
k2=[v l]';
v=v-dt/2*(k1(2,:)'+k2(2,:)')+F;
% v=v-dt/2*(k*u+k*u)+F;
sol(n+1,:)=[0 , u' , 0];

end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler modificado')
}}
La gráfica que se obtiene es:

[[Archivo:ImagenCodigo3GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler modificado]]
===Comparativa de los métodos e interpretación===
Dados los órdenes de precisión de cada uno de los métodos para aproximar la función solución del sistema de ecuaciones planteado por el método de las diferencias finitas, el más preciso es el método del trapecio (orden 2), seguido por Euler modificado y por último Euler explicito (orden 1).

==Aproximación mediante el método de Fourier==
===Con un término de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con un término

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con un término (Q=1)

Q=1;
for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);

end
[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con un término')

}}
[[Archivo:ImagenCodigo4GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con un término]]

===Con tres términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con tres términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con tres términos (Q=3)

Q=3;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x,p.*p);
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con tres términos')
}}
[[Archivo:ImagenCodigo5GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con tres términos]]

===Con cinco términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con cinco términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con cinco términos (Q=5)

Q=5;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con cinco términos')
}}
[[Archivo:ImagenCodigo6GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con cinco términos]]

===Con diez términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con diez términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con diez términos (Q=10)

Q=10;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con diez términos')
}}
[[Archivo:ImagenCodigo7GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con diez términos]]

===Con veinte términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con veinte términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con diez términos (Q=20)

Q=20;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);

end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con veinte términos')
}}
[[Archivo:ImagenCodigo8.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con veinte términos]]

===Comparativa de la diferencia de términos e interpretación===
Podemos observar que a medida que aumentamos el número de términos de la serie de Fourier, la forma de la solución pasa de ser más redondeada a más recta, pareciéndose cada vez más a la solución obtenida por el método de las diferencias finitas.

==Comparativa de diferencias finitas y Fourier==
Se puede obsevar que al aumentar el numéro de términos con Fourier, la solución se va aproximando cada vez más a la obtenida mediante el método de las diferencias finitas.

=Energía del cable=

La energía del cable viene definida por la siguiente función:: <math> E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx </math>
Mediante el programa que se muestra a continuación procederemos a dibujar la gráfica de la función <math>E(t)</math>, que cuantifica la energía del cable. Previamente, y recurriendo al concepto de derivada, aproximaremos los valores que toman <math>u_{t}</math> y <math>u_{x}</math>, para así después poder calcular la integral que define <math>E(t)</math>.
{{matlab|codigo=Aquí hay que subir CodigoEnergia.m
}}

[[Archivo:ImagenCodigoEnergiaGRP9.jpg|border|Variación de la energía del cable a lo largo del tiempo]]

Como se puede apreciar en la gráfica, esta es de apariencia trigonométrica, periódica en el tiempo (el periodo es de 10 segundos).

'''¿Qué quiere decir esto?'''

Esta gráfica no es otra cosa si no la '''representación de la variación de la energía potencial en valor absoluto a lo largo del tiempo''', tomando como referencia de los potenciales un desplazamiento vertical nulo <math>u=0</math>. Para una interpretación más clara solo '''nos fijaremos en el punto situado a la mitad del cable''', en el que para <math>t=0</math> el desplazamiento vertical es de dos metros.

Inicialmente soltamos el cable desde una altura de dos metros y con velocidad nula, es decir, que en este instante la energía potencial será máxima mientras que la energía cinética será inexistente. A medida que desciende el cable, este incrementa su energía cinética, contrarrestando así la energía potencial que decrece, hasta que en <math>t=5</math> la energía potencial es nula (al igual que el desplazamiento vertical); esto quiere decir que la energía cinética es máxima en este instante. La energía potencial continúa aumentando su valor hasta que en <math>t=10</math> se vuelve a hacer casi máxima (el desplazamiento vertical del cable esta cerca de los dos metros hacia abajo) mientras que la cinética nula. Comienza nuevamente a aumentar la cinética y decrecer la potencial, de forma que para <math>t=15</math> tanto el desplazamiento vertical como la energía potencial son nulos (la cinética vuelve a ser casi máxima); vuelve a aumentar nuevamente y en <math>t=20</math> el cable vuelve a estar próximo a los dos metros de altura con una energía potencial cercana a la máxima.

La '''energía mecánica''' (suma de la cinética y la potencial) '''tiende a conservarse, pero esto no ocurre''', puesto que hay pérdidas (puede ser por el rozamiento con el aire, etc), a simple vista no se aprecia en la gráfica, pero si se ejecuta el programa y se comprueba la matriz “E” (almacena la energía en cada instante) se observa que esta va decreciendo a lo largo del tiempo.

El siguiente código auxiliar muestra el desplazamiento vertical del punto central en cada instante, a simple vista no se aprecia, pero si se ejecuta y comprueba la matriz “sol” se verá que el desplazamiento vertical es máximo al principio (2 metros), no volviendo a recuperar el cable este valor (aunque se quede próximo a él en ciertos instantes, para ello hay que observar la columna 51 de la matriz), '''el movimiento se va “amortiguando” reduciéndose lentamente el desplazamiento vertical y la energía'''. Si estudiáramos un intervalo de tiempo mayor esto se apreciaría mucho mejor, pruebe a cambiar en el siguiente código T=40 por T=500, es evidente, cada vez se aleja más de los dos metros iniciales que marcan los límites de la gráfica (observe la segunda imagen).
{{matlab|codigo=Aquí hay que subir CodigoAuxEnergia.m
}}
[[Archivo:ImagenCodigoAuxEnergia.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=40]]
[[Archivo:ImagenCodigoAuxEnergiaT500.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=500]]

=Ejemplos de aplicación de cables en estructuras civiles=
==Cable sumergido en el mar (medio viscoso)==
En este caso el cable está sumergido en un medio viscoso que produce amortiguamiento.
El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}+au_{t}=0; x∈[0,10]; t∈[0,40] \\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

Vamos a estudiar lo que ocurre con la energía en estas condiciones, dependiendo del parámetro a.Vamos a dar a a los valores de 0, 1 ,4, 10 y 100.
En el caso de que a valga cero estamos ante el mismo caso que si no estuviese sumergido, el cable permanecería en reposo.

Se ve que, al aumentar la a, la energía se disipa más rápidamente.

==Cable sujeto a una estructura sometida a vibraciones periódicas==
Ahora nos preguntamos qué ocurriría si el extremo derecho del cable está sujeto a una estructura que sufre vibraciones periódicas. La frecuencia de estas vibraciones es de <math>F_{0}</math> Herzios. Tomamos por ejemplo <math>f(t) = sin(2 \pi F_{0})</math>. Vamos a calcular el comportamiento de la energía cuando el cable parte inicialmente del reposo y tomamos un tiempo t. Repetiremos el proceso con distintos valores de <math>F_{0}</math> para observar si se produce alguna diferencia.

Tomamos <math>F_{0}</math> igual a 1/L+0,01 Hz, 1/L-0,01 Hz y 1/L Hz.

El problema nos quedaría planteado de la siguiente manera:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,60]\\
u(0,t)=0; u_{x}(10,t)=sin(2 \pi F_{0})\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

==Cable sujeto a un aparato que envía una respuesta a las vibraciones==
En este apartado vamos a calcular el comportamiento de la energía para el valor b suponiendo ahora que el extremo derecho del cable está sujeto a un aparato que envía una respuesta a la vibración que recibe de manera que la condición de contorno ahora es, tomando como condición inicial la del apartado 2 y observando en qué cambia

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u_{x}(10,t)=b u(10,0) \\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

=Visualización del movimiento del cable=
Si el siguiente código basado en el método de diferencias finitas (mediante el método del trapecio) analizado anteriormente (en el apartado 2) se ejecuta en MATLAB (Octave no admite animaciones), se puede apreciar en la primera imagen que proporciona el programa la solución al problema enunciado en el apartado 2. Lo siguiente en salir por pantalla es una animación de la solución obtenida a lo largo del tiempo, de forma que cuando se agota el intervalo de tiempo en el que estudiamos el problema, queda la misma imagen que la primera que nos proporciona este programa; con este programa vemos como la onda va "viajando". La tercera imagen en salir por pantalla es la representación del desplazamiento vertical que sufre el cable a lo largo del tiempo, es decir, la vibración de todos los puntos de nuestro cable.
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')

%% Animación.

% Muestra el desplazamiento vertical de cada punto del cable en cada instante

figure
[n,m]=size(sol);
aux=zeros(size(sol));
for j=1:n
% Se va analizando el desplazamiento vertical en cada instante (j) de
% todos los puntos (i) del cable
for i=1:m
aux(j,i)=sol(j,i);
end
surf(xx,tt,aux); % Es la representación de todos los puntos (i) del cable en un instante dado
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
zlim([-2 2])
MFRAMES=getframe; % Guarda una foto del cable en cada instante
end
movie=MFRAMES;

figure
for j=1:n
% Se va analizando el desplazamiento vertical en cada instante (j) de
% todos los puntos (i) del cable
for i=1:m
aux(j,i)=sol(j,i);
end
plot(x,aux(j,:)); % Es la representación de todos los puntos (i) del cable en un instante dado
xlabel('Cable')
ylabel('Desplazamiento vertical')
ylim([-2 2])
MFRAMES=getframe; % Guarda una foto del cable en cada instante
end
movie=MFRAMES;
}}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T19:16:32Z

Jose Manuel Corbí Garrido: /* Con diez términos de la serie */

{{ Beta }}
{{ TrabajoED | Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]]|

María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

=Introducción y modelización del problema=

En este artículo analizaremos el comportamiento de un cable tenso de 10 metros de longitud en una estructura civil. El cable se encontrará sujeto por los extremos y será sometido a pequeñas vibraciones. Para estudiar el problema haremos ciertas hipótesis que simplificarán el problema, y que en muchos casos dan una buena aproximación del comportamiento real que tendría un cable como el que aquí estudiamos.

En cuanto a las '''características mecánicas del cable''' haremos las siguientes '''hipótesis''':

* El cable tendrá una sección despreciable frente a su longitud.
* Será un cable homogéneo, es decir, la masa por unidad de volumen será constante.
* Dicho cable también será perfectamente flexible, es decir, solo ofrecerá resistencia a esfuerzos tangenciales (en su dirección longitudinal), siendo por ello nula su resistencia frente a esfuerzos de flexión y corte.
* El cable al someterlo a tracción es inextensible, es decir, que es lo suficientemente rígido (en dirección longitudinal) como para poder despreciar su extensibilidad, por el contrario al someterlo a compresión no ofrece resistencia y se arruga.

Por otra parte el cable será sometido a '''pequeñas vibraciones''', las cuales podemos modelizar mediante la '''ecuación de ondas''' <math>u_{tt}-u_{xx}=0</math>. Estas vibraciones, al ser pequeñas, son despreciables frente a su longitud total; preocupándonos únicamente por los desplazamientos verticales (fruto de la tracción) que sufre el cable. Estos desplazamientos verticales serán modelizados mediante la función <math>u(x,t)</math>, dependiente de dos variables; la <math>x</math> que indica la distancia al extremo de la cuerda tomado como referencia, y la <math>t</math> que define el tiempo.

Las condiciones de frontera en este caso nos indican que el cable no sufrirá desplazamiento vertical en sus extremos <math>u(0,t)=0; u(10,t)=0</math>. Por otra parte las condiciones iniciales <math>u(x,0)=0; u_{t}(x,0)=0</math> se traducen en que inicialmente el desplazamiento vertical y la velocidad del cable en todos sus puntos es cero.

El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

En los sucesivos apartados estudiaremos el comportamiento del cable en diferentes situaciones, que se traducirán en diferentes condiciones para el problema (cambios en las condiciones y en la ecuación de ondas con la que modelizamos las pequeñas vibraciones).

=Desplazamiento vertical del cable=
Sujetando el cable desde el centro y desplazándolo dos metros en la dirección perpendicular, estudiaremos el comportamiento del cable una vez lo soltemos. Estará sometido a unas pequeñas vibraciones, que originarán desplazamientos verticales a lo largo del cable que aproximaremos mediante la función <math>u(x,t)</math>. Teniendo en cuenta el desplazamiento vertical de dos metros y que el cable es soltado con una velocidad nula, el problema que vamos a estudiar es el siguiente:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

A continuación estudiaremos el desplazamiento vertical al que se ve cada punto del cable expuesto y para cada instante de tiempo. Para ello aproximaremos la ecuación diferencial que modeliza el comportamiento de dicho cable mediante los métodos de diferencias finitas y mediante series de Fourier, obteniendo gráficas que reflejan la posición de cada punto del cable en cada instante.
==Aproximación mediante diferencias finitas==
Dentro del método de diferencias finitas se nos plantea un sistema de ecuaciones diferenciales, que resolveremos aplicando diferentes métodos de los cuales luego compararemos su precisión a la hora de aproximar el comportamiento del cable y su estabilidad.

El primer método que utilizaremos para resolver el sistema de ecuaciones diferenciales será el '''método del trapecio''', el cual es '''incondicionalmente estable'''. Para implementar este método y observar el comportamiento del cable en una gráfica usaremos el siguiente código en matlab:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')
}}
Obteniendo la gráfica:

[[Archivo:ImagenCodigo1GRP9.jpeg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método del trapecio]]

Usaremos también el '''método de Euler''', el cual es un método explícito y por tanto dependemos de la ecuación en que pueda ser o no factible de ejecutar. El siguiente código matlab nos permite ejecutarlo:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler explícito

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler explícito

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for j=1:length(t)-1
u=u+dt*v+F;
v=v-dt*k*u+F;
sol(j+1,:)=[0 , u' , 0];
end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler explícito')
}}
Obteniendo el siguiente gráfico del comportamiento del cable respecto al tiempo:

[[Archivo:ImagenCodigo2GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler explícito]]

Por último método para resolver el sistema planteado mediante diferencias finitas, resolvemos numéricamente mediante el '''método de Euler modificado''', el cual tiene más precisión que el método de Euler. El código que se muestra a continuación ejecuta dicho método y nos proporciona una gráfica del cable:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler modificado

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler modificado

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for n=1:length(t)-1
l=k*u;
% calculamos k1=f(t_n,y_n)
k1=[ v l]';
% Calculamos k2=f(t_n+h/2,y_n+1/2k1*h)
%xp=x(n)+dt; %argumento de tiempo
%solp=sol+dt*k1; %argumento y
k2=[v l]';
%Calculamos k3=f(t_n+h/2,y_n+1/2k2*h)
%mismo argumento tiempo que antes

u=u+dt/2*(k1(1,:)'+k2(1,:)')+F;
l=+k*u;
k1=[v l]';
k2=[v l]';
v=v-dt/2*(k1(2,:)'+k2(2,:)')+F;
% v=v-dt/2*(k*u+k*u)+F;
sol(n+1,:)=[0 , u' , 0];

end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler modificado')
}}
La gráfica que se obtiene es:

[[Archivo:ImagenCodigo3GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler modificado]]
===Comparativa de los métodos e interpretación===
Dados los órdenes de precisión de cada uno de los métodos para aproximar la función solución del sistema de ecuaciones planteado por el método de las diferencias finitas, el más preciso es el método del trapecio (orden 2), seguido por Euler modificado y por último Euler explicito (orden 1).

==Aproximación mediante el método de Fourier==
===Con un término de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con un término

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con un término (Q=1)

Q=1;
for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);

end
[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con un término')

}}
[[Archivo:ImagenCodigo4GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con un término]]

===Con tres términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con tres términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con tres términos (Q=3)

Q=3;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x,p.*p);
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con tres términos')
}}
[[Archivo:ImagenCodigo5GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con tres términos]]

===Con cinco términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con cinco términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con cinco términos (Q=5)

Q=5;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con cinco términos')
}}
[[Archivo:ImagenCodigo6GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con cinco términos]]

===Con diez términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con diez términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con diez términos (Q=10)

Q=10;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con diez términos')
}}
[[Archivo:ImagenCodigo7GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con diez términos]]

===Con veinte términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo8.m
}}
[[Archivo:ImagenCodigo8.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con veinte términos]]
===Comparativa de la diferencia de términos e interpretación===
Podemos observar que a medida que aumentamos el número de términos de la serie de Fourier, la forma de la solución pasa de ser más redondeada a más recta, pareciéndose cada vez más a la solución obtenida por el método de las diferencias finitas.

==Comparativa de diferencias finitas y Fourier==
Se puede obsevar que al aumentar el numéro de términos con Fourier, la solución se va aproximando cada vez más a la obtenida mediante el método de las diferencias finitas.

=Energía del cable=

La energía del cable viene definida por la siguiente función:: <math> E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx </math>
Mediante el programa que se muestra a continuación procederemos a dibujar la gráfica de la función <math>E(t)</math>, que cuantifica la energía del cable. Previamente, y recurriendo al concepto de derivada, aproximaremos los valores que toman <math>u_{t}</math> y <math>u_{x}</math>, para así después poder calcular la integral que define <math>E(t)</math>.
{{matlab|codigo=Aquí hay que subir CodigoEnergia.m
}}

[[Archivo:ImagenCodigoEnergiaGRP9.jpg|border|Variación de la energía del cable a lo largo del tiempo]]

Como se puede apreciar en la gráfica, esta es de apariencia trigonométrica, periódica en el tiempo (el periodo es de 10 segundos).

'''¿Qué quiere decir esto?'''

Esta gráfica no es otra cosa si no la '''representación de la variación de la energía potencial en valor absoluto a lo largo del tiempo''', tomando como referencia de los potenciales un desplazamiento vertical nulo <math>u=0</math>. Para una interpretación más clara solo '''nos fijaremos en el punto situado a la mitad del cable''', en el que para <math>t=0</math> el desplazamiento vertical es de dos metros.

Inicialmente soltamos el cable desde una altura de dos metros y con velocidad nula, es decir, que en este instante la energía potencial será máxima mientras que la energía cinética será inexistente. A medida que desciende el cable, este incrementa su energía cinética, contrarrestando así la energía potencial que decrece, hasta que en <math>t=5</math> la energía potencial es nula (al igual que el desplazamiento vertical); esto quiere decir que la energía cinética es máxima en este instante. La energía potencial continúa aumentando su valor hasta que en <math>t=10</math> se vuelve a hacer casi máxima (el desplazamiento vertical del cable esta cerca de los dos metros hacia abajo) mientras que la cinética nula. Comienza nuevamente a aumentar la cinética y decrecer la potencial, de forma que para <math>t=15</math> tanto el desplazamiento vertical como la energía potencial son nulos (la cinética vuelve a ser casi máxima); vuelve a aumentar nuevamente y en <math>t=20</math> el cable vuelve a estar próximo a los dos metros de altura con una energía potencial cercana a la máxima.

La '''energía mecánica''' (suma de la cinética y la potencial) '''tiende a conservarse, pero esto no ocurre''', puesto que hay pérdidas (puede ser por el rozamiento con el aire, etc), a simple vista no se aprecia en la gráfica, pero si se ejecuta el programa y se comprueba la matriz “E” (almacena la energía en cada instante) se observa que esta va decreciendo a lo largo del tiempo.

El siguiente código auxiliar muestra el desplazamiento vertical del punto central en cada instante, a simple vista no se aprecia, pero si se ejecuta y comprueba la matriz “sol” se verá que el desplazamiento vertical es máximo al principio (2 metros), no volviendo a recuperar el cable este valor (aunque se quede próximo a él en ciertos instantes, para ello hay que observar la columna 51 de la matriz), '''el movimiento se va “amortiguando” reduciéndose lentamente el desplazamiento vertical y la energía'''. Si estudiáramos un intervalo de tiempo mayor esto se apreciaría mucho mejor, pruebe a cambiar en el siguiente código T=40 por T=500, es evidente, cada vez se aleja más de los dos metros iniciales que marcan los límites de la gráfica (observe la segunda imagen).
{{matlab|codigo=Aquí hay que subir CodigoAuxEnergia.m
}}
[[Archivo:ImagenCodigoAuxEnergia.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=40]]
[[Archivo:ImagenCodigoAuxEnergiaT500.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=500]]

=Ejemplos de aplicación de cables en estructuras civiles=
==Cable sumergido en el mar (medio viscoso)==
En este caso el cable está sumergido en un medio viscoso que produce amortiguamiento.
El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}+au_{t}=0; x∈[0,10]; t∈[0,40] \\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

Vamos a estudiar lo que ocurre con la energía en estas condiciones, dependiendo del parámetro a.Vamos a dar a a los valores de 0, 1 ,4, 10 y 100.
En el caso de que a valga cero estamos ante el mismo caso que si no estuviese sumergido, el cable permanecería en reposo.

Se ve que, al aumentar la a, la energía se disipa más rápidamente.

==Cable sujeto a una estructura sometida a vibraciones periódicas==
Ahora nos preguntamos qué ocurriría si el extremo derecho del cable está sujeto a una estructura que sufre vibraciones periódicas. La frecuencia de estas vibraciones es de <math>F_{0}</math> Herzios. Tomamos por ejemplo <math>f(t) = sin(2 \pi F_{0})</math>. Vamos a calcular el comportamiento de la energía cuando el cable parte inicialmente del reposo y tomamos un tiempo t. Repetiremos el proceso con distintos valores de <math>F_{0}</math> para observar si se produce alguna diferencia.

Tomamos <math>F_{0}</math> igual a 1/L+0,01 Hz, 1/L-0,01 Hz y 1/L Hz.

El problema nos quedaría planteado de la siguiente manera:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,60]\\
u(0,t)=0; u_{x}(10,t)=sin(2 \pi F_{0})\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

==Cable sujeto a un aparato que envía una respuesta a las vibraciones==
En este apartado vamos a calcular el comportamiento de la energía para el valor b suponiendo ahora que el extremo derecho del cable está sujeto a un aparato que envía una respuesta a la vibración que recibe de manera que la condición de contorno ahora es, tomando como condición inicial la del apartado 2 y observando en qué cambia

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u_{x}(10,t)=b u(10,0) \\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

=Visualización del movimiento del cable=
Si el siguiente código basado en el método de diferencias finitas (mediante el método del trapecio) analizado anteriormente (en el apartado 2) se ejecuta en MATLAB (Octave no admite animaciones), se puede apreciar en la primera imagen que proporciona el programa la solución al problema enunciado en el apartado 2. Lo siguiente en salir por pantalla es una animación de la solución obtenida a lo largo del tiempo, de forma que cuando se agota el intervalo de tiempo en el que estudiamos el problema, queda la misma imagen que la primera que nos proporciona este programa; con este programa vemos como la onda va "viajando". La tercera imagen en salir por pantalla es la representación del desplazamiento vertical que sufre el cable a lo largo del tiempo, es decir, la vibración de todos los puntos de nuestro cable.
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')

%% Animación.

% Muestra el desplazamiento vertical de cada punto del cable en cada instante

figure
[n,m]=size(sol);
aux=zeros(size(sol));
for j=1:n
% Se va analizando el desplazamiento vertical en cada instante (j) de
% todos los puntos (i) del cable
for i=1:m
aux(j,i)=sol(j,i);
end
surf(xx,tt,aux); % Es la representación de todos los puntos (i) del cable en un instante dado
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
zlim([-2 2])
MFRAMES=getframe; % Guarda una foto del cable en cada instante
end
movie=MFRAMES;

figure
for j=1:n
% Se va analizando el desplazamiento vertical en cada instante (j) de
% todos los puntos (i) del cable
for i=1:m
aux(j,i)=sol(j,i);
end
plot(x,aux(j,:)); % Es la representación de todos los puntos (i) del cable en un instante dado
xlabel('Cable')
ylabel('Desplazamiento vertical')
ylim([-2 2])
MFRAMES=getframe; % Guarda una foto del cable en cada instante
end
movie=MFRAMES;
}}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T19:16:15Z

Jose Manuel Corbí Garrido: /* Con cinco términos de la serie */

{{ Beta }}
{{ TrabajoED | Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]]|

María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

=Introducción y modelización del problema=

En este artículo analizaremos el comportamiento de un cable tenso de 10 metros de longitud en una estructura civil. El cable se encontrará sujeto por los extremos y será sometido a pequeñas vibraciones. Para estudiar el problema haremos ciertas hipótesis que simplificarán el problema, y que en muchos casos dan una buena aproximación del comportamiento real que tendría un cable como el que aquí estudiamos.

En cuanto a las '''características mecánicas del cable''' haremos las siguientes '''hipótesis''':

* El cable tendrá una sección despreciable frente a su longitud.
* Será un cable homogéneo, es decir, la masa por unidad de volumen será constante.
* Dicho cable también será perfectamente flexible, es decir, solo ofrecerá resistencia a esfuerzos tangenciales (en su dirección longitudinal), siendo por ello nula su resistencia frente a esfuerzos de flexión y corte.
* El cable al someterlo a tracción es inextensible, es decir, que es lo suficientemente rígido (en dirección longitudinal) como para poder despreciar su extensibilidad, por el contrario al someterlo a compresión no ofrece resistencia y se arruga.

Por otra parte el cable será sometido a '''pequeñas vibraciones''', las cuales podemos modelizar mediante la '''ecuación de ondas''' <math>u_{tt}-u_{xx}=0</math>. Estas vibraciones, al ser pequeñas, son despreciables frente a su longitud total; preocupándonos únicamente por los desplazamientos verticales (fruto de la tracción) que sufre el cable. Estos desplazamientos verticales serán modelizados mediante la función <math>u(x,t)</math>, dependiente de dos variables; la <math>x</math> que indica la distancia al extremo de la cuerda tomado como referencia, y la <math>t</math> que define el tiempo.

Las condiciones de frontera en este caso nos indican que el cable no sufrirá desplazamiento vertical en sus extremos <math>u(0,t)=0; u(10,t)=0</math>. Por otra parte las condiciones iniciales <math>u(x,0)=0; u_{t}(x,0)=0</math> se traducen en que inicialmente el desplazamiento vertical y la velocidad del cable en todos sus puntos es cero.

El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

En los sucesivos apartados estudiaremos el comportamiento del cable en diferentes situaciones, que se traducirán en diferentes condiciones para el problema (cambios en las condiciones y en la ecuación de ondas con la que modelizamos las pequeñas vibraciones).

=Desplazamiento vertical del cable=
Sujetando el cable desde el centro y desplazándolo dos metros en la dirección perpendicular, estudiaremos el comportamiento del cable una vez lo soltemos. Estará sometido a unas pequeñas vibraciones, que originarán desplazamientos verticales a lo largo del cable que aproximaremos mediante la función <math>u(x,t)</math>. Teniendo en cuenta el desplazamiento vertical de dos metros y que el cable es soltado con una velocidad nula, el problema que vamos a estudiar es el siguiente:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

A continuación estudiaremos el desplazamiento vertical al que se ve cada punto del cable expuesto y para cada instante de tiempo. Para ello aproximaremos la ecuación diferencial que modeliza el comportamiento de dicho cable mediante los métodos de diferencias finitas y mediante series de Fourier, obteniendo gráficas que reflejan la posición de cada punto del cable en cada instante.
==Aproximación mediante diferencias finitas==
Dentro del método de diferencias finitas se nos plantea un sistema de ecuaciones diferenciales, que resolveremos aplicando diferentes métodos de los cuales luego compararemos su precisión a la hora de aproximar el comportamiento del cable y su estabilidad.

El primer método que utilizaremos para resolver el sistema de ecuaciones diferenciales será el '''método del trapecio''', el cual es '''incondicionalmente estable'''. Para implementar este método y observar el comportamiento del cable en una gráfica usaremos el siguiente código en matlab:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')
}}
Obteniendo la gráfica:

[[Archivo:ImagenCodigo1GRP9.jpeg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método del trapecio]]

Usaremos también el '''método de Euler''', el cual es un método explícito y por tanto dependemos de la ecuación en que pueda ser o no factible de ejecutar. El siguiente código matlab nos permite ejecutarlo:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler explícito

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler explícito

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for j=1:length(t)-1
u=u+dt*v+F;
v=v-dt*k*u+F;
sol(j+1,:)=[0 , u' , 0];
end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler explícito')
}}
Obteniendo el siguiente gráfico del comportamiento del cable respecto al tiempo:

[[Archivo:ImagenCodigo2GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler explícito]]

Por último método para resolver el sistema planteado mediante diferencias finitas, resolvemos numéricamente mediante el '''método de Euler modificado''', el cual tiene más precisión que el método de Euler. El código que se muestra a continuación ejecuta dicho método y nos proporciona una gráfica del cable:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler modificado

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler modificado

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for n=1:length(t)-1
l=k*u;
% calculamos k1=f(t_n,y_n)
k1=[ v l]';
% Calculamos k2=f(t_n+h/2,y_n+1/2k1*h)
%xp=x(n)+dt; %argumento de tiempo
%solp=sol+dt*k1; %argumento y
k2=[v l]';
%Calculamos k3=f(t_n+h/2,y_n+1/2k2*h)
%mismo argumento tiempo que antes

u=u+dt/2*(k1(1,:)'+k2(1,:)')+F;
l=+k*u;
k1=[v l]';
k2=[v l]';
v=v-dt/2*(k1(2,:)'+k2(2,:)')+F;
% v=v-dt/2*(k*u+k*u)+F;
sol(n+1,:)=[0 , u' , 0];

end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler modificado')
}}
La gráfica que se obtiene es:

[[Archivo:ImagenCodigo3GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler modificado]]
===Comparativa de los métodos e interpretación===
Dados los órdenes de precisión de cada uno de los métodos para aproximar la función solución del sistema de ecuaciones planteado por el método de las diferencias finitas, el más preciso es el método del trapecio (orden 2), seguido por Euler modificado y por último Euler explicito (orden 1).

==Aproximación mediante el método de Fourier==
===Con un término de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con un término

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con un término (Q=1)

Q=1;
for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);

end
[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con un término')

}}
[[Archivo:ImagenCodigo4GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con un término]]

===Con tres términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con tres términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con tres términos (Q=3)

Q=3;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x,p.*p);
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con tres términos')
}}
[[Archivo:ImagenCodigo5GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con tres términos]]

===Con cinco términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con cinco términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con cinco términos (Q=5)

Q=5;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con cinco términos')
}}
[[Archivo:ImagenCodigo6GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con cinco términos]]

===Con diez términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con diez términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con diez términos (Q=10)

Q=10;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
%fourier2(k)=(1/(k�pi ) ) �trapz (x , ve l o c idad . � p) /
%trapz (x , p . � p) ;
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con diez términos')
}}
[[Archivo:ImagenCodigo7GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con diez términos]]
===Con veinte términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo8.m
}}
[[Archivo:ImagenCodigo8.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con veinte términos]]
===Comparativa de la diferencia de términos e interpretación===
Podemos observar que a medida que aumentamos el número de términos de la serie de Fourier, la forma de la solución pasa de ser más redondeada a más recta, pareciéndose cada vez más a la solución obtenida por el método de las diferencias finitas.

==Comparativa de diferencias finitas y Fourier==
Se puede obsevar que al aumentar el numéro de términos con Fourier, la solución se va aproximando cada vez más a la obtenida mediante el método de las diferencias finitas.

=Energía del cable=

La energía del cable viene definida por la siguiente función:: <math> E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx </math>
Mediante el programa que se muestra a continuación procederemos a dibujar la gráfica de la función <math>E(t)</math>, que cuantifica la energía del cable. Previamente, y recurriendo al concepto de derivada, aproximaremos los valores que toman <math>u_{t}</math> y <math>u_{x}</math>, para así después poder calcular la integral que define <math>E(t)</math>.
{{matlab|codigo=Aquí hay que subir CodigoEnergia.m
}}

[[Archivo:ImagenCodigoEnergiaGRP9.jpg|border|Variación de la energía del cable a lo largo del tiempo]]

Como se puede apreciar en la gráfica, esta es de apariencia trigonométrica, periódica en el tiempo (el periodo es de 10 segundos).

'''¿Qué quiere decir esto?'''

Esta gráfica no es otra cosa si no la '''representación de la variación de la energía potencial en valor absoluto a lo largo del tiempo''', tomando como referencia de los potenciales un desplazamiento vertical nulo <math>u=0</math>. Para una interpretación más clara solo '''nos fijaremos en el punto situado a la mitad del cable''', en el que para <math>t=0</math> el desplazamiento vertical es de dos metros.

Inicialmente soltamos el cable desde una altura de dos metros y con velocidad nula, es decir, que en este instante la energía potencial será máxima mientras que la energía cinética será inexistente. A medida que desciende el cable, este incrementa su energía cinética, contrarrestando así la energía potencial que decrece, hasta que en <math>t=5</math> la energía potencial es nula (al igual que el desplazamiento vertical); esto quiere decir que la energía cinética es máxima en este instante. La energía potencial continúa aumentando su valor hasta que en <math>t=10</math> se vuelve a hacer casi máxima (el desplazamiento vertical del cable esta cerca de los dos metros hacia abajo) mientras que la cinética nula. Comienza nuevamente a aumentar la cinética y decrecer la potencial, de forma que para <math>t=15</math> tanto el desplazamiento vertical como la energía potencial son nulos (la cinética vuelve a ser casi máxima); vuelve a aumentar nuevamente y en <math>t=20</math> el cable vuelve a estar próximo a los dos metros de altura con una energía potencial cercana a la máxima.

La '''energía mecánica''' (suma de la cinética y la potencial) '''tiende a conservarse, pero esto no ocurre''', puesto que hay pérdidas (puede ser por el rozamiento con el aire, etc), a simple vista no se aprecia en la gráfica, pero si se ejecuta el programa y se comprueba la matriz “E” (almacena la energía en cada instante) se observa que esta va decreciendo a lo largo del tiempo.

El siguiente código auxiliar muestra el desplazamiento vertical del punto central en cada instante, a simple vista no se aprecia, pero si se ejecuta y comprueba la matriz “sol” se verá que el desplazamiento vertical es máximo al principio (2 metros), no volviendo a recuperar el cable este valor (aunque se quede próximo a él en ciertos instantes, para ello hay que observar la columna 51 de la matriz), '''el movimiento se va “amortiguando” reduciéndose lentamente el desplazamiento vertical y la energía'''. Si estudiáramos un intervalo de tiempo mayor esto se apreciaría mucho mejor, pruebe a cambiar en el siguiente código T=40 por T=500, es evidente, cada vez se aleja más de los dos metros iniciales que marcan los límites de la gráfica (observe la segunda imagen).
{{matlab|codigo=Aquí hay que subir CodigoAuxEnergia.m
}}
[[Archivo:ImagenCodigoAuxEnergia.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=40]]
[[Archivo:ImagenCodigoAuxEnergiaT500.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=500]]

=Ejemplos de aplicación de cables en estructuras civiles=
==Cable sumergido en el mar (medio viscoso)==
En este caso el cable está sumergido en un medio viscoso que produce amortiguamiento.
El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}+au_{t}=0; x∈[0,10]; t∈[0,40] \\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

Vamos a estudiar lo que ocurre con la energía en estas condiciones, dependiendo del parámetro a.Vamos a dar a a los valores de 0, 1 ,4, 10 y 100.
En el caso de que a valga cero estamos ante el mismo caso que si no estuviese sumergido, el cable permanecería en reposo.

Se ve que, al aumentar la a, la energía se disipa más rápidamente.

==Cable sujeto a una estructura sometida a vibraciones periódicas==
Ahora nos preguntamos qué ocurriría si el extremo derecho del cable está sujeto a una estructura que sufre vibraciones periódicas. La frecuencia de estas vibraciones es de <math>F_{0}</math> Herzios. Tomamos por ejemplo <math>f(t) = sin(2 \pi F_{0})</math>. Vamos a calcular el comportamiento de la energía cuando el cable parte inicialmente del reposo y tomamos un tiempo t. Repetiremos el proceso con distintos valores de <math>F_{0}</math> para observar si se produce alguna diferencia.

Tomamos <math>F_{0}</math> igual a 1/L+0,01 Hz, 1/L-0,01 Hz y 1/L Hz.

El problema nos quedaría planteado de la siguiente manera:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,60]\\
u(0,t)=0; u_{x}(10,t)=sin(2 \pi F_{0})\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

==Cable sujeto a un aparato que envía una respuesta a las vibraciones==
En este apartado vamos a calcular el comportamiento de la energía para el valor b suponiendo ahora que el extremo derecho del cable está sujeto a un aparato que envía una respuesta a la vibración que recibe de manera que la condición de contorno ahora es, tomando como condición inicial la del apartado 2 y observando en qué cambia

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u_{x}(10,t)=b u(10,0) \\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

=Visualización del movimiento del cable=
Si el siguiente código basado en el método de diferencias finitas (mediante el método del trapecio) analizado anteriormente (en el apartado 2) se ejecuta en MATLAB (Octave no admite animaciones), se puede apreciar en la primera imagen que proporciona el programa la solución al problema enunciado en el apartado 2. Lo siguiente en salir por pantalla es una animación de la solución obtenida a lo largo del tiempo, de forma que cuando se agota el intervalo de tiempo en el que estudiamos el problema, queda la misma imagen que la primera que nos proporciona este programa; con este programa vemos como la onda va "viajando". La tercera imagen en salir por pantalla es la representación del desplazamiento vertical que sufre el cable a lo largo del tiempo, es decir, la vibración de todos los puntos de nuestro cable.
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')

%% Animación.

% Muestra el desplazamiento vertical de cada punto del cable en cada instante

figure
[n,m]=size(sol);
aux=zeros(size(sol));
for j=1:n
% Se va analizando el desplazamiento vertical en cada instante (j) de
% todos los puntos (i) del cable
for i=1:m
aux(j,i)=sol(j,i);
end
surf(xx,tt,aux); % Es la representación de todos los puntos (i) del cable en un instante dado
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
zlim([-2 2])
MFRAMES=getframe; % Guarda una foto del cable en cada instante
end
movie=MFRAMES;

figure
for j=1:n
% Se va analizando el desplazamiento vertical en cada instante (j) de
% todos los puntos (i) del cable
for i=1:m
aux(j,i)=sol(j,i);
end
plot(x,aux(j,:)); % Es la representación de todos los puntos (i) del cable en un instante dado
xlabel('Cable')
ylabel('Desplazamiento vertical')
ylim([-2 2])
MFRAMES=getframe; % Guarda una foto del cable en cada instante
end
movie=MFRAMES;
}}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T19:15:42Z

Jose Manuel Corbí Garrido: /* Con un término de la serie */

{{ Beta }}
{{ TrabajoED | Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]]|

María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

=Introducción y modelización del problema=

En este artículo analizaremos el comportamiento de un cable tenso de 10 metros de longitud en una estructura civil. El cable se encontrará sujeto por los extremos y será sometido a pequeñas vibraciones. Para estudiar el problema haremos ciertas hipótesis que simplificarán el problema, y que en muchos casos dan una buena aproximación del comportamiento real que tendría un cable como el que aquí estudiamos.

En cuanto a las '''características mecánicas del cable''' haremos las siguientes '''hipótesis''':

* El cable tendrá una sección despreciable frente a su longitud.
* Será un cable homogéneo, es decir, la masa por unidad de volumen será constante.
* Dicho cable también será perfectamente flexible, es decir, solo ofrecerá resistencia a esfuerzos tangenciales (en su dirección longitudinal), siendo por ello nula su resistencia frente a esfuerzos de flexión y corte.
* El cable al someterlo a tracción es inextensible, es decir, que es lo suficientemente rígido (en dirección longitudinal) como para poder despreciar su extensibilidad, por el contrario al someterlo a compresión no ofrece resistencia y se arruga.

Por otra parte el cable será sometido a '''pequeñas vibraciones''', las cuales podemos modelizar mediante la '''ecuación de ondas''' <math>u_{tt}-u_{xx}=0</math>. Estas vibraciones, al ser pequeñas, son despreciables frente a su longitud total; preocupándonos únicamente por los desplazamientos verticales (fruto de la tracción) que sufre el cable. Estos desplazamientos verticales serán modelizados mediante la función <math>u(x,t)</math>, dependiente de dos variables; la <math>x</math> que indica la distancia al extremo de la cuerda tomado como referencia, y la <math>t</math> que define el tiempo.

Las condiciones de frontera en este caso nos indican que el cable no sufrirá desplazamiento vertical en sus extremos <math>u(0,t)=0; u(10,t)=0</math>. Por otra parte las condiciones iniciales <math>u(x,0)=0; u_{t}(x,0)=0</math> se traducen en que inicialmente el desplazamiento vertical y la velocidad del cable en todos sus puntos es cero.

El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

En los sucesivos apartados estudiaremos el comportamiento del cable en diferentes situaciones, que se traducirán en diferentes condiciones para el problema (cambios en las condiciones y en la ecuación de ondas con la que modelizamos las pequeñas vibraciones).

=Desplazamiento vertical del cable=
Sujetando el cable desde el centro y desplazándolo dos metros en la dirección perpendicular, estudiaremos el comportamiento del cable una vez lo soltemos. Estará sometido a unas pequeñas vibraciones, que originarán desplazamientos verticales a lo largo del cable que aproximaremos mediante la función <math>u(x,t)</math>. Teniendo en cuenta el desplazamiento vertical de dos metros y que el cable es soltado con una velocidad nula, el problema que vamos a estudiar es el siguiente:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

A continuación estudiaremos el desplazamiento vertical al que se ve cada punto del cable expuesto y para cada instante de tiempo. Para ello aproximaremos la ecuación diferencial que modeliza el comportamiento de dicho cable mediante los métodos de diferencias finitas y mediante series de Fourier, obteniendo gráficas que reflejan la posición de cada punto del cable en cada instante.
==Aproximación mediante diferencias finitas==
Dentro del método de diferencias finitas se nos plantea un sistema de ecuaciones diferenciales, que resolveremos aplicando diferentes métodos de los cuales luego compararemos su precisión a la hora de aproximar el comportamiento del cable y su estabilidad.

El primer método que utilizaremos para resolver el sistema de ecuaciones diferenciales será el '''método del trapecio''', el cual es '''incondicionalmente estable'''. Para implementar este método y observar el comportamiento del cable en una gráfica usaremos el siguiente código en matlab:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')
}}
Obteniendo la gráfica:

[[Archivo:ImagenCodigo1GRP9.jpeg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método del trapecio]]

Usaremos también el '''método de Euler''', el cual es un método explícito y por tanto dependemos de la ecuación en que pueda ser o no factible de ejecutar. El siguiente código matlab nos permite ejecutarlo:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler explícito

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler explícito

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for j=1:length(t)-1
u=u+dt*v+F;
v=v-dt*k*u+F;
sol(j+1,:)=[0 , u' , 0];
end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler explícito')
}}
Obteniendo el siguiente gráfico del comportamiento del cable respecto al tiempo:

[[Archivo:ImagenCodigo2GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler explícito]]

Por último método para resolver el sistema planteado mediante diferencias finitas, resolvemos numéricamente mediante el '''método de Euler modificado''', el cual tiene más precisión que el método de Euler. El código que se muestra a continuación ejecuta dicho método y nos proporciona una gráfica del cable:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler modificado

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler modificado

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for n=1:length(t)-1
l=k*u;
% calculamos k1=f(t_n,y_n)
k1=[ v l]';
% Calculamos k2=f(t_n+h/2,y_n+1/2k1*h)
%xp=x(n)+dt; %argumento de tiempo
%solp=sol+dt*k1; %argumento y
k2=[v l]';
%Calculamos k3=f(t_n+h/2,y_n+1/2k2*h)
%mismo argumento tiempo que antes

u=u+dt/2*(k1(1,:)'+k2(1,:)')+F;
l=+k*u;
k1=[v l]';
k2=[v l]';
v=v-dt/2*(k1(2,:)'+k2(2,:)')+F;
% v=v-dt/2*(k*u+k*u)+F;
sol(n+1,:)=[0 , u' , 0];

end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler modificado')
}}
La gráfica que se obtiene es:

[[Archivo:ImagenCodigo3GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler modificado]]
===Comparativa de los métodos e interpretación===
Dados los órdenes de precisión de cada uno de los métodos para aproximar la función solución del sistema de ecuaciones planteado por el método de las diferencias finitas, el más preciso es el método del trapecio (orden 2), seguido por Euler modificado y por último Euler explicito (orden 1).

==Aproximación mediante el método de Fourier==
===Con un término de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con un término

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con un término (Q=1)

Q=1;
for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);

end
[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con un término')

}}
[[Archivo:ImagenCodigo4GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con un término]]

===Con tres términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con tres términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con tres términos (Q=3)

Q=3;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x,p.*p);
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con tres términos')
}}
[[Archivo:ImagenCodigo5GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con tres términos]]

===Con cinco términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con cinco términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con cinco términos (Q=5)

Q=5;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
%fourier2(k)=(1/(k�pi ) ) �trapz (x , ve l o c idad . � p) /
%trapz (x , p . � p) ;
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con cinco términos')
}}
[[Archivo:ImagenCodigo6GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con cinco términos]]
===Con diez términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con diez términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con diez términos (Q=10)

Q=10;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
%fourier2(k)=(1/(k�pi ) ) �trapz (x , ve l o c idad . � p) /
%trapz (x , p . � p) ;
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con diez términos')
}}
[[Archivo:ImagenCodigo7GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con diez términos]]
===Con veinte términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo8.m
}}
[[Archivo:ImagenCodigo8.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con veinte términos]]
===Comparativa de la diferencia de términos e interpretación===
Podemos observar que a medida que aumentamos el número de términos de la serie de Fourier, la forma de la solución pasa de ser más redondeada a más recta, pareciéndose cada vez más a la solución obtenida por el método de las diferencias finitas.

==Comparativa de diferencias finitas y Fourier==
Se puede obsevar que al aumentar el numéro de términos con Fourier, la solución se va aproximando cada vez más a la obtenida mediante el método de las diferencias finitas.

=Energía del cable=

La energía del cable viene definida por la siguiente función:: <math> E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx </math>
Mediante el programa que se muestra a continuación procederemos a dibujar la gráfica de la función <math>E(t)</math>, que cuantifica la energía del cable. Previamente, y recurriendo al concepto de derivada, aproximaremos los valores que toman <math>u_{t}</math> y <math>u_{x}</math>, para así después poder calcular la integral que define <math>E(t)</math>.
{{matlab|codigo=Aquí hay que subir CodigoEnergia.m
}}

[[Archivo:ImagenCodigoEnergiaGRP9.jpg|border|Variación de la energía del cable a lo largo del tiempo]]

Como se puede apreciar en la gráfica, esta es de apariencia trigonométrica, periódica en el tiempo (el periodo es de 10 segundos).

'''¿Qué quiere decir esto?'''

Esta gráfica no es otra cosa si no la '''representación de la variación de la energía potencial en valor absoluto a lo largo del tiempo''', tomando como referencia de los potenciales un desplazamiento vertical nulo <math>u=0</math>. Para una interpretación más clara solo '''nos fijaremos en el punto situado a la mitad del cable''', en el que para <math>t=0</math> el desplazamiento vertical es de dos metros.

Inicialmente soltamos el cable desde una altura de dos metros y con velocidad nula, es decir, que en este instante la energía potencial será máxima mientras que la energía cinética será inexistente. A medida que desciende el cable, este incrementa su energía cinética, contrarrestando así la energía potencial que decrece, hasta que en <math>t=5</math> la energía potencial es nula (al igual que el desplazamiento vertical); esto quiere decir que la energía cinética es máxima en este instante. La energía potencial continúa aumentando su valor hasta que en <math>t=10</math> se vuelve a hacer casi máxima (el desplazamiento vertical del cable esta cerca de los dos metros hacia abajo) mientras que la cinética nula. Comienza nuevamente a aumentar la cinética y decrecer la potencial, de forma que para <math>t=15</math> tanto el desplazamiento vertical como la energía potencial son nulos (la cinética vuelve a ser casi máxima); vuelve a aumentar nuevamente y en <math>t=20</math> el cable vuelve a estar próximo a los dos metros de altura con una energía potencial cercana a la máxima.

La '''energía mecánica''' (suma de la cinética y la potencial) '''tiende a conservarse, pero esto no ocurre''', puesto que hay pérdidas (puede ser por el rozamiento con el aire, etc), a simple vista no se aprecia en la gráfica, pero si se ejecuta el programa y se comprueba la matriz “E” (almacena la energía en cada instante) se observa que esta va decreciendo a lo largo del tiempo.

El siguiente código auxiliar muestra el desplazamiento vertical del punto central en cada instante, a simple vista no se aprecia, pero si se ejecuta y comprueba la matriz “sol” se verá que el desplazamiento vertical es máximo al principio (2 metros), no volviendo a recuperar el cable este valor (aunque se quede próximo a él en ciertos instantes, para ello hay que observar la columna 51 de la matriz), '''el movimiento se va “amortiguando” reduciéndose lentamente el desplazamiento vertical y la energía'''. Si estudiáramos un intervalo de tiempo mayor esto se apreciaría mucho mejor, pruebe a cambiar en el siguiente código T=40 por T=500, es evidente, cada vez se aleja más de los dos metros iniciales que marcan los límites de la gráfica (observe la segunda imagen).
{{matlab|codigo=Aquí hay que subir CodigoAuxEnergia.m
}}
[[Archivo:ImagenCodigoAuxEnergia.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=40]]
[[Archivo:ImagenCodigoAuxEnergiaT500.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=500]]

=Ejemplos de aplicación de cables en estructuras civiles=
==Cable sumergido en el mar (medio viscoso)==
En este caso el cable está sumergido en un medio viscoso que produce amortiguamiento.
El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}+au_{t}=0; x∈[0,10]; t∈[0,40] \\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

Vamos a estudiar lo que ocurre con la energía en estas condiciones, dependiendo del parámetro a.Vamos a dar a a los valores de 0, 1 ,4, 10 y 100.
En el caso de que a valga cero estamos ante el mismo caso que si no estuviese sumergido, el cable permanecería en reposo.

Se ve que, al aumentar la a, la energía se disipa más rápidamente.

==Cable sujeto a una estructura sometida a vibraciones periódicas==
Ahora nos preguntamos qué ocurriría si el extremo derecho del cable está sujeto a una estructura que sufre vibraciones periódicas. La frecuencia de estas vibraciones es de <math>F_{0}</math> Herzios. Tomamos por ejemplo <math>f(t) = sin(2 \pi F_{0})</math>. Vamos a calcular el comportamiento de la energía cuando el cable parte inicialmente del reposo y tomamos un tiempo t. Repetiremos el proceso con distintos valores de <math>F_{0}</math> para observar si se produce alguna diferencia.

Tomamos <math>F_{0}</math> igual a 1/L+0,01 Hz, 1/L-0,01 Hz y 1/L Hz.

El problema nos quedaría planteado de la siguiente manera:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,60]\\
u(0,t)=0; u_{x}(10,t)=sin(2 \pi F_{0})\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

==Cable sujeto a un aparato que envía una respuesta a las vibraciones==
En este apartado vamos a calcular el comportamiento de la energía para el valor b suponiendo ahora que el extremo derecho del cable está sujeto a un aparato que envía una respuesta a la vibración que recibe de manera que la condición de contorno ahora es, tomando como condición inicial la del apartado 2 y observando en qué cambia

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u_{x}(10,t)=b u(10,0) \\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

=Visualización del movimiento del cable=
Si el siguiente código basado en el método de diferencias finitas (mediante el método del trapecio) analizado anteriormente (en el apartado 2) se ejecuta en MATLAB (Octave no admite animaciones), se puede apreciar en la primera imagen que proporciona el programa la solución al problema enunciado en el apartado 2. Lo siguiente en salir por pantalla es una animación de la solución obtenida a lo largo del tiempo, de forma que cuando se agota el intervalo de tiempo en el que estudiamos el problema, queda la misma imagen que la primera que nos proporciona este programa; con este programa vemos como la onda va "viajando". La tercera imagen en salir por pantalla es la representación del desplazamiento vertical que sufre el cable a lo largo del tiempo, es decir, la vibración de todos los puntos de nuestro cable.
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')

%% Animación.

% Muestra el desplazamiento vertical de cada punto del cable en cada instante

figure
[n,m]=size(sol);
aux=zeros(size(sol));
for j=1:n
% Se va analizando el desplazamiento vertical en cada instante (j) de
% todos los puntos (i) del cable
for i=1:m
aux(j,i)=sol(j,i);
end
surf(xx,tt,aux); % Es la representación de todos los puntos (i) del cable en un instante dado
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
zlim([-2 2])
MFRAMES=getframe; % Guarda una foto del cable en cada instante
end
movie=MFRAMES;

figure
for j=1:n
% Se va analizando el desplazamiento vertical en cada instante (j) de
% todos los puntos (i) del cable
for i=1:m
aux(j,i)=sol(j,i);
end
plot(x,aux(j,:)); % Es la representación de todos los puntos (i) del cable en un instante dado
xlabel('Cable')
ylabel('Desplazamiento vertical')
ylim([-2 2])
MFRAMES=getframe; % Guarda una foto del cable en cada instante
end
movie=MFRAMES;
}}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T19:15:11Z

Jose Manuel Corbí Garrido: /* Con tres términos de la serie */

{{ Beta }}
{{ TrabajoED | Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]]|

María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

=Introducción y modelización del problema=

En este artículo analizaremos el comportamiento de un cable tenso de 10 metros de longitud en una estructura civil. El cable se encontrará sujeto por los extremos y será sometido a pequeñas vibraciones. Para estudiar el problema haremos ciertas hipótesis que simplificarán el problema, y que en muchos casos dan una buena aproximación del comportamiento real que tendría un cable como el que aquí estudiamos.

En cuanto a las '''características mecánicas del cable''' haremos las siguientes '''hipótesis''':

* El cable tendrá una sección despreciable frente a su longitud.
* Será un cable homogéneo, es decir, la masa por unidad de volumen será constante.
* Dicho cable también será perfectamente flexible, es decir, solo ofrecerá resistencia a esfuerzos tangenciales (en su dirección longitudinal), siendo por ello nula su resistencia frente a esfuerzos de flexión y corte.
* El cable al someterlo a tracción es inextensible, es decir, que es lo suficientemente rígido (en dirección longitudinal) como para poder despreciar su extensibilidad, por el contrario al someterlo a compresión no ofrece resistencia y se arruga.

Por otra parte el cable será sometido a '''pequeñas vibraciones''', las cuales podemos modelizar mediante la '''ecuación de ondas''' <math>u_{tt}-u_{xx}=0</math>. Estas vibraciones, al ser pequeñas, son despreciables frente a su longitud total; preocupándonos únicamente por los desplazamientos verticales (fruto de la tracción) que sufre el cable. Estos desplazamientos verticales serán modelizados mediante la función <math>u(x,t)</math>, dependiente de dos variables; la <math>x</math> que indica la distancia al extremo de la cuerda tomado como referencia, y la <math>t</math> que define el tiempo.

Las condiciones de frontera en este caso nos indican que el cable no sufrirá desplazamiento vertical en sus extremos <math>u(0,t)=0; u(10,t)=0</math>. Por otra parte las condiciones iniciales <math>u(x,0)=0; u_{t}(x,0)=0</math> se traducen en que inicialmente el desplazamiento vertical y la velocidad del cable en todos sus puntos es cero.

El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

En los sucesivos apartados estudiaremos el comportamiento del cable en diferentes situaciones, que se traducirán en diferentes condiciones para el problema (cambios en las condiciones y en la ecuación de ondas con la que modelizamos las pequeñas vibraciones).

=Desplazamiento vertical del cable=
Sujetando el cable desde el centro y desplazándolo dos metros en la dirección perpendicular, estudiaremos el comportamiento del cable una vez lo soltemos. Estará sometido a unas pequeñas vibraciones, que originarán desplazamientos verticales a lo largo del cable que aproximaremos mediante la función <math>u(x,t)</math>. Teniendo en cuenta el desplazamiento vertical de dos metros y que el cable es soltado con una velocidad nula, el problema que vamos a estudiar es el siguiente:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

A continuación estudiaremos el desplazamiento vertical al que se ve cada punto del cable expuesto y para cada instante de tiempo. Para ello aproximaremos la ecuación diferencial que modeliza el comportamiento de dicho cable mediante los métodos de diferencias finitas y mediante series de Fourier, obteniendo gráficas que reflejan la posición de cada punto del cable en cada instante.
==Aproximación mediante diferencias finitas==
Dentro del método de diferencias finitas se nos plantea un sistema de ecuaciones diferenciales, que resolveremos aplicando diferentes métodos de los cuales luego compararemos su precisión a la hora de aproximar el comportamiento del cable y su estabilidad.

El primer método que utilizaremos para resolver el sistema de ecuaciones diferenciales será el '''método del trapecio''', el cual es '''incondicionalmente estable'''. Para implementar este método y observar el comportamiento del cable en una gráfica usaremos el siguiente código en matlab:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')
}}
Obteniendo la gráfica:

[[Archivo:ImagenCodigo1GRP9.jpeg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método del trapecio]]

Usaremos también el '''método de Euler''', el cual es un método explícito y por tanto dependemos de la ecuación en que pueda ser o no factible de ejecutar. El siguiente código matlab nos permite ejecutarlo:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler explícito

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler explícito

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for j=1:length(t)-1
u=u+dt*v+F;
v=v-dt*k*u+F;
sol(j+1,:)=[0 , u' , 0];
end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler explícito')
}}
Obteniendo el siguiente gráfico del comportamiento del cable respecto al tiempo:

[[Archivo:ImagenCodigo2GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler explícito]]

Por último método para resolver el sistema planteado mediante diferencias finitas, resolvemos numéricamente mediante el '''método de Euler modificado''', el cual tiene más precisión que el método de Euler. El código que se muestra a continuación ejecuta dicho método y nos proporciona una gráfica del cable:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler modificado

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler modificado

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for n=1:length(t)-1
l=k*u;
% calculamos k1=f(t_n,y_n)
k1=[ v l]';
% Calculamos k2=f(t_n+h/2,y_n+1/2k1*h)
%xp=x(n)+dt; %argumento de tiempo
%solp=sol+dt*k1; %argumento y
k2=[v l]';
%Calculamos k3=f(t_n+h/2,y_n+1/2k2*h)
%mismo argumento tiempo que antes

u=u+dt/2*(k1(1,:)'+k2(1,:)')+F;
l=+k*u;
k1=[v l]';
k2=[v l]';
v=v-dt/2*(k1(2,:)'+k2(2,:)')+F;
% v=v-dt/2*(k*u+k*u)+F;
sol(n+1,:)=[0 , u' , 0];

end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler modificado')
}}
La gráfica que se obtiene es:

[[Archivo:ImagenCodigo3GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler modificado]]
===Comparativa de los métodos e interpretación===
Dados los órdenes de precisión de cada uno de los métodos para aproximar la función solución del sistema de ecuaciones planteado por el método de las diferencias finitas, el más preciso es el método del trapecio (orden 2), seguido por Euler modificado y por último Euler explicito (orden 1).

==Aproximación mediante el método de Fourier==
===Con un término de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con un término

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con un término (Q=1)

Q=1;
for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
%fourier2(k)=(1/(k*pi))*trapz(x,velocidad.*p) /
%trapz (x,p.*p) ;
end
[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con un término')

}}
[[Archivo:ImagenCodigo4GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con un término]]

===Con tres términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con tres términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con tres términos (Q=3)

Q=3;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x,p.*p);
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con tres términos')
}}
[[Archivo:ImagenCodigo5GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con tres términos]]

===Con cinco términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con cinco términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con cinco términos (Q=5)

Q=5;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
%fourier2(k)=(1/(k�pi ) ) �trapz (x , ve l o c idad . � p) /
%trapz (x , p . � p) ;
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con cinco términos')
}}
[[Archivo:ImagenCodigo6GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con cinco términos]]
===Con diez términos de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con diez términos

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con diez términos (Q=10)

Q=10;

for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
%fourier2(k)=(1/(k�pi ) ) �trapz (x , ve l o c idad . � p) /
%trapz (x , p . � p) ;
end

[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con diez términos')
}}
[[Archivo:ImagenCodigo7GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con diez términos]]
===Con veinte términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo8.m
}}
[[Archivo:ImagenCodigo8.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con veinte términos]]
===Comparativa de la diferencia de términos e interpretación===
Podemos observar que a medida que aumentamos el número de términos de la serie de Fourier, la forma de la solución pasa de ser más redondeada a más recta, pareciéndose cada vez más a la solución obtenida por el método de las diferencias finitas.

==Comparativa de diferencias finitas y Fourier==
Se puede obsevar que al aumentar el numéro de términos con Fourier, la solución se va aproximando cada vez más a la obtenida mediante el método de las diferencias finitas.

=Energía del cable=

La energía del cable viene definida por la siguiente función:: <math> E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx </math>
Mediante el programa que se muestra a continuación procederemos a dibujar la gráfica de la función <math>E(t)</math>, que cuantifica la energía del cable. Previamente, y recurriendo al concepto de derivada, aproximaremos los valores que toman <math>u_{t}</math> y <math>u_{x}</math>, para así después poder calcular la integral que define <math>E(t)</math>.
{{matlab|codigo=Aquí hay que subir CodigoEnergia.m
}}

[[Archivo:ImagenCodigoEnergiaGRP9.jpg|border|Variación de la energía del cable a lo largo del tiempo]]

Como se puede apreciar en la gráfica, esta es de apariencia trigonométrica, periódica en el tiempo (el periodo es de 10 segundos).

'''¿Qué quiere decir esto?'''

Esta gráfica no es otra cosa si no la '''representación de la variación de la energía potencial en valor absoluto a lo largo del tiempo''', tomando como referencia de los potenciales un desplazamiento vertical nulo <math>u=0</math>. Para una interpretación más clara solo '''nos fijaremos en el punto situado a la mitad del cable''', en el que para <math>t=0</math> el desplazamiento vertical es de dos metros.

Inicialmente soltamos el cable desde una altura de dos metros y con velocidad nula, es decir, que en este instante la energía potencial será máxima mientras que la energía cinética será inexistente. A medida que desciende el cable, este incrementa su energía cinética, contrarrestando así la energía potencial que decrece, hasta que en <math>t=5</math> la energía potencial es nula (al igual que el desplazamiento vertical); esto quiere decir que la energía cinética es máxima en este instante. La energía potencial continúa aumentando su valor hasta que en <math>t=10</math> se vuelve a hacer casi máxima (el desplazamiento vertical del cable esta cerca de los dos metros hacia abajo) mientras que la cinética nula. Comienza nuevamente a aumentar la cinética y decrecer la potencial, de forma que para <math>t=15</math> tanto el desplazamiento vertical como la energía potencial son nulos (la cinética vuelve a ser casi máxima); vuelve a aumentar nuevamente y en <math>t=20</math> el cable vuelve a estar próximo a los dos metros de altura con una energía potencial cercana a la máxima.

La '''energía mecánica''' (suma de la cinética y la potencial) '''tiende a conservarse, pero esto no ocurre''', puesto que hay pérdidas (puede ser por el rozamiento con el aire, etc), a simple vista no se aprecia en la gráfica, pero si se ejecuta el programa y se comprueba la matriz “E” (almacena la energía en cada instante) se observa que esta va decreciendo a lo largo del tiempo.

El siguiente código auxiliar muestra el desplazamiento vertical del punto central en cada instante, a simple vista no se aprecia, pero si se ejecuta y comprueba la matriz “sol” se verá que el desplazamiento vertical es máximo al principio (2 metros), no volviendo a recuperar el cable este valor (aunque se quede próximo a él en ciertos instantes, para ello hay que observar la columna 51 de la matriz), '''el movimiento se va “amortiguando” reduciéndose lentamente el desplazamiento vertical y la energía'''. Si estudiáramos un intervalo de tiempo mayor esto se apreciaría mucho mejor, pruebe a cambiar en el siguiente código T=40 por T=500, es evidente, cada vez se aleja más de los dos metros iniciales que marcan los límites de la gráfica (observe la segunda imagen).
{{matlab|codigo=Aquí hay que subir CodigoAuxEnergia.m
}}
[[Archivo:ImagenCodigoAuxEnergia.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=40]]
[[Archivo:ImagenCodigoAuxEnergiaT500.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=500]]

=Ejemplos de aplicación de cables en estructuras civiles=
==Cable sumergido en el mar (medio viscoso)==
En este caso el cable está sumergido en un medio viscoso que produce amortiguamiento.
El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}+au_{t}=0; x∈[0,10]; t∈[0,40] \\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

Vamos a estudiar lo que ocurre con la energía en estas condiciones, dependiendo del parámetro a.Vamos a dar a a los valores de 0, 1 ,4, 10 y 100.
En el caso de que a valga cero estamos ante el mismo caso que si no estuviese sumergido, el cable permanecería en reposo.

Se ve que, al aumentar la a, la energía se disipa más rápidamente.

==Cable sujeto a una estructura sometida a vibraciones periódicas==
Ahora nos preguntamos qué ocurriría si el extremo derecho del cable está sujeto a una estructura que sufre vibraciones periódicas. La frecuencia de estas vibraciones es de <math>F_{0}</math> Herzios. Tomamos por ejemplo <math>f(t) = sin(2 \pi F_{0})</math>. Vamos a calcular el comportamiento de la energía cuando el cable parte inicialmente del reposo y tomamos un tiempo t. Repetiremos el proceso con distintos valores de <math>F_{0}</math> para observar si se produce alguna diferencia.

Tomamos <math>F_{0}</math> igual a 1/L+0,01 Hz, 1/L-0,01 Hz y 1/L Hz.

El problema nos quedaría planteado de la siguiente manera:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,60]\\
u(0,t)=0; u_{x}(10,t)=sin(2 \pi F_{0})\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

==Cable sujeto a un aparato que envía una respuesta a las vibraciones==
En este apartado vamos a calcular el comportamiento de la energía para el valor b suponiendo ahora que el extremo derecho del cable está sujeto a un aparato que envía una respuesta a la vibración que recibe de manera que la condición de contorno ahora es, tomando como condición inicial la del apartado 2 y observando en qué cambia

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u_{x}(10,t)=b u(10,0) \\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

=Visualización del movimiento del cable=
Si el siguiente código basado en el método de diferencias finitas (mediante el método del trapecio) analizado anteriormente (en el apartado 2) se ejecuta en MATLAB (Octave no admite animaciones), se puede apreciar en la primera imagen que proporciona el programa la solución al problema enunciado en el apartado 2. Lo siguiente en salir por pantalla es una animación de la solución obtenida a lo largo del tiempo, de forma que cuando se agota el intervalo de tiempo en el que estudiamos el problema, queda la misma imagen que la primera que nos proporciona este programa; con este programa vemos como la onda va "viajando". La tercera imagen en salir por pantalla es la representación del desplazamiento vertical que sufre el cable a lo largo del tiempo, es decir, la vibración de todos los puntos de nuestro cable.
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')

%% Animación.

% Muestra el desplazamiento vertical de cada punto del cable en cada instante

figure
[n,m]=size(sol);
aux=zeros(size(sol));
for j=1:n
% Se va analizando el desplazamiento vertical en cada instante (j) de
% todos los puntos (i) del cable
for i=1:m
aux(j,i)=sol(j,i);
end
surf(xx,tt,aux); % Es la representación de todos los puntos (i) del cable en un instante dado
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
zlim([-2 2])
MFRAMES=getframe; % Guarda una foto del cable en cada instante
end
movie=MFRAMES;

figure
for j=1:n
% Se va analizando el desplazamiento vertical en cada instante (j) de
% todos los puntos (i) del cable
for i=1:m
aux(j,i)=sol(j,i);
end
plot(x,aux(j,:)); % Es la representación de todos los puntos (i) del cable en un instante dado
xlabel('Cable')
ylabel('Desplazamiento vertical')
ylim([-2 2])
MFRAMES=getframe; % Guarda una foto del cable en cada instante
end
movie=MFRAMES;
}}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T19:11:54Z

Jose Manuel Corbí Garrido: /* Con un término de la serie */

{{ Beta }}
{{ TrabajoED | Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]]|

María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

=Introducción y modelización del problema=

En este artículo analizaremos el comportamiento de un cable tenso de 10 metros de longitud en una estructura civil. El cable se encontrará sujeto por los extremos y será sometido a pequeñas vibraciones. Para estudiar el problema haremos ciertas hipótesis que simplificarán el problema, y que en muchos casos dan una buena aproximación del comportamiento real que tendría un cable como el que aquí estudiamos.

En cuanto a las '''características mecánicas del cable''' haremos las siguientes '''hipótesis''':

* El cable tendrá una sección despreciable frente a su longitud.
* Será un cable homogéneo, es decir, la masa por unidad de volumen será constante.
* Dicho cable también será perfectamente flexible, es decir, solo ofrecerá resistencia a esfuerzos tangenciales (en su dirección longitudinal), siendo por ello nula su resistencia frente a esfuerzos de flexión y corte.
* El cable al someterlo a tracción es inextensible, es decir, que es lo suficientemente rígido (en dirección longitudinal) como para poder despreciar su extensibilidad, por el contrario al someterlo a compresión no ofrece resistencia y se arruga.

Por otra parte el cable será sometido a '''pequeñas vibraciones''', las cuales podemos modelizar mediante la '''ecuación de ondas''' <math>u_{tt}-u_{xx}=0</math>. Estas vibraciones, al ser pequeñas, son despreciables frente a su longitud total; preocupándonos únicamente por los desplazamientos verticales (fruto de la tracción) que sufre el cable. Estos desplazamientos verticales serán modelizados mediante la función <math>u(x,t)</math>, dependiente de dos variables; la <math>x</math> que indica la distancia al extremo de la cuerda tomado como referencia, y la <math>t</math> que define el tiempo.

Las condiciones de frontera en este caso nos indican que el cable no sufrirá desplazamiento vertical en sus extremos <math>u(0,t)=0; u(10,t)=0</math>. Por otra parte las condiciones iniciales <math>u(x,0)=0; u_{t}(x,0)=0</math> se traducen en que inicialmente el desplazamiento vertical y la velocidad del cable en todos sus puntos es cero.

El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

En los sucesivos apartados estudiaremos el comportamiento del cable en diferentes situaciones, que se traducirán en diferentes condiciones para el problema (cambios en las condiciones y en la ecuación de ondas con la que modelizamos las pequeñas vibraciones).

=Desplazamiento vertical del cable=
Sujetando el cable desde el centro y desplazándolo dos metros en la dirección perpendicular, estudiaremos el comportamiento del cable una vez lo soltemos. Estará sometido a unas pequeñas vibraciones, que originarán desplazamientos verticales a lo largo del cable que aproximaremos mediante la función <math>u(x,t)</math>. Teniendo en cuenta el desplazamiento vertical de dos metros y que el cable es soltado con una velocidad nula, el problema que vamos a estudiar es el siguiente:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

A continuación estudiaremos el desplazamiento vertical al que se ve cada punto del cable expuesto y para cada instante de tiempo. Para ello aproximaremos la ecuación diferencial que modeliza el comportamiento de dicho cable mediante los métodos de diferencias finitas y mediante series de Fourier, obteniendo gráficas que reflejan la posición de cada punto del cable en cada instante.
==Aproximación mediante diferencias finitas==
Dentro del método de diferencias finitas se nos plantea un sistema de ecuaciones diferenciales, que resolveremos aplicando diferentes métodos de los cuales luego compararemos su precisión a la hora de aproximar el comportamiento del cable y su estabilidad.

El primer método que utilizaremos para resolver el sistema de ecuaciones diferenciales será el '''método del trapecio''', el cual es '''incondicionalmente estable'''. Para implementar este método y observar el comportamiento del cable en una gráfica usaremos el siguiente código en matlab:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')
}}
Obteniendo la gráfica:

[[Archivo:ImagenCodigo1GRP9.jpeg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método del trapecio]]

Usaremos también el '''método de Euler''', el cual es un método explícito y por tanto dependemos de la ecuación en que pueda ser o no factible de ejecutar. El siguiente código matlab nos permite ejecutarlo:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler explícito

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler explícito

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for j=1:length(t)-1
u=u+dt*v+F;
v=v-dt*k*u+F;
sol(j+1,:)=[0 , u' , 0];
end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler explícito')
}}
Obteniendo el siguiente gráfico del comportamiento del cable respecto al tiempo:

[[Archivo:ImagenCodigo2GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler explícito]]

Por último método para resolver el sistema planteado mediante diferencias finitas, resolvemos numéricamente mediante el '''método de Euler''' modificado, el cual tiene más precisión que el método de Euler. El código que se muestra a continuación ejecuta dicho método y nos proporciona una gráfica del cable:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler modificado

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler modificado

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for n=1:length(t)-1
l=k*u;
% calculamos k1=f(t_n,y_n)
k1=[ v l]';
% Calculamos k2=f(t_n+h/2,y_n+1/2k1*h)
%xp=x(n)+dt; %argumento de tiempo
%solp=sol+dt*k1; %argumento y
k2=[v l]';
%Calculamos k3=f(t_n+h/2,y_n+1/2k2*h)
%mismo argumento tiempo que antes

u=u+dt/2*(k1(1,:)'+k2(1,:)')+F;
l=+k*u;
k1=[v l]';
k2=[v l]';
v=v-dt/2*(k1(2,:)'+k2(2,:)')+F;
% v=v-dt/2*(k*u+k*u)+F;
sol(n+1,:)=[0 , u' , 0];

end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler modificado')
}}
La gráfica que se obtiene es:

[[Archivo:ImagenCodigo3GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler modificado]]
===Comparativa de los métodos e interpretación===
Dados los órdenes de precisión de cada uno de los métodos para aproximar la función solución del sistema de ecuaciones planteado por el método de las diferencias finitas, el más preciso es el método del trapecio (orden 2), seguido por Euler modificado y por último Euler explicito (orden 1).

==Aproximación mediante el método de Fourier==
===Con un término de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con un término

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con un término (Q=1)

Q=1;
for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
%fourier2(k)=(1/(k*pi))*trapz(x,velocidad.*p) /
%trapz (x,p.*p) ;
end
[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con un término')

}}
[[Archivo:ImagenCodigo4GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con un término]]

===Con tres términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo5.m
}}
[[Archivo:ImagenCodigo5GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con tres términos]]
===Con cinco términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo6.m
}}
[[Archivo:ImagenCodigo6GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con cinco términos]]
===Con diez términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo7.m
}}
[[Archivo:ImagenCodigo7GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con diez términos]]
===Con veinte términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo8.m
}}
[[Archivo:ImagenCodigo8.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con veinte términos]]
===Comparativa de la diferencia de términos e interpretación===
Podemos observar que a medida que aumentamos el número de términos de la serie de Fourier, la forma de la solución pasa de ser más redondeada a más recta, pareciéndose cada vez más a la solución obtenida por el método de las diferencias finitas.

==Comparativa de diferencias finitas y Fourier==
Se puede obsevar que al aumentar el numéro de términos con Fourier, la solución se va aproximando cada vez más a la obtenida mediante el método de las diferencias finitas.

=Energía del cable=

La energía del cable viene definida por la siguiente función:: <math> E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx </math>
Mediante el programa que se muestra a continuación procederemos a dibujar la gráfica de la función <math>E(t)</math>, que cuantifica la energía del cable. Previamente, y recurriendo al concepto de derivada, aproximaremos los valores que toman <math>u_{t}</math> y <math>u_{x}</math>, para así después poder calcular la integral que define <math>E(t)</math>.
{{matlab|codigo=Aquí hay que subir CodigoEnergia.m
}}

[[Archivo:ImagenCodigoEnergiaGRP9.jpg|border|Variación de la energía del cable a lo largo del tiempo]]

Como se puede apreciar en la gráfica, esta es de apariencia trigonométrica, periódica en el tiempo (el periodo es de 10 segundos).

'''¿Qué quiere decir esto?'''

Esta gráfica no es otra cosa si no la '''representación de la variación de la energía potencial en valor absoluto a lo largo del tiempo''', tomando como referencia de los potenciales un desplazamiento vertical nulo <math>u=0</math>. Para una interpretación más clara solo '''nos fijaremos en el punto situado a la mitad del cable''', en el que para <math>t=0</math> el desplazamiento vertical es de dos metros.

Inicialmente soltamos el cable desde una altura de dos metros y con velocidad nula, es decir, que en este instante la energía potencial será máxima mientras que la energía cinética será inexistente. A medida que desciende el cable, este incrementa su energía cinética, contrarrestando así la energía potencial que decrece, hasta que en <math>t=5</math> la energía potencial es nula (al igual que el desplazamiento vertical); esto quiere decir que la energía cinética es máxima en este instante. La energía potencial continúa aumentando su valor hasta que en <math>t=10</math> se vuelve a hacer casi máxima (el desplazamiento vertical del cable esta cerca de los dos metros hacia abajo) mientras que la cinética nula. Comienza nuevamente a aumentar la cinética y decrecer la potencial, de forma que para <math>t=15</math> tanto el desplazamiento vertical como la energía potencial son nulos (la cinética vuelve a ser casi máxima); vuelve a aumentar nuevamente y en <math>t=20</math> el cable vuelve a estar próximo a los dos metros de altura con una energía potencial cercana a la máxima.

La '''energía mecánica''' (suma de la cinética y la potencial) '''tiende a conservarse, pero esto no ocurre''', puesto que hay pérdidas (puede ser por el rozamiento con el aire, etc), a simple vista no se aprecia en la gráfica, pero si se ejecuta el programa y se comprueba la matriz “E” (almacena la energía en cada instante) se observa que esta va decreciendo a lo largo del tiempo.

El siguiente código auxiliar muestra el desplazamiento vertical del punto central en cada instante, a simple vista no se aprecia, pero si se ejecuta y comprueba la matriz “sol” se verá que el desplazamiento vertical es máximo al principio (2 metros), no volviendo a recuperar el cable este valor (aunque se quede próximo a él en ciertos instantes, para ello hay que observar la columna 51 de la matriz), '''el movimiento se va “amortiguando” reduciéndose lentamente el desplazamiento vertical y la energía'''. Si estudiáramos un intervalo de tiempo mayor esto se apreciaría mucho mejor, pruebe a cambiar en el siguiente código T=40 por T=500, es evidente, cada vez se aleja más de los dos metros iniciales que marcan los límites de la gráfica (observe la segunda imagen).
{{matlab|codigo=Aquí hay que subir CodigoAuxEnergia.m
}}
[[Archivo:ImagenCodigoAuxEnergia.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=40]]
[[Archivo:ImagenCodigoAuxEnergiaT500.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=500]]

=Ejemplos de aplicación de cables en estructuras civiles=
==Cable sumergido en el mar (medio viscoso)==
En este caso el cable está sumergido en un medio viscoso que produce amortiguamiento.
El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}+au_{t}=0; x∈[0,10]; t∈[0,40] \\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

Vamos a estudiar lo que ocurre con la energía en estas condiciones, dependiendo del parámetro a.Vamos a dar a a los valores de 0, 1 ,4, 10 y 100.
En el caso de que a valga cero estamos ante el mismo caso que si no estuviese sumergido, el cable permanecería en reposo.

Se ve que, al aumentar la a, la energía se disipa más rápidamente.

==Cable sujeto a una estructura sometida a vibraciones periódicas==
Ahora nos preguntamos qué ocurriría si el extremo derecho del cable está sujeto a una estructura que sufre vibraciones periódicas. La frecuencia de estas vibraciones es de <math>F_{0}</math> Herzios. Tomamos por ejemplo <math>f(t) = sin(2 \pi F_{0})</math>. Vamos a calcular el comportamiento de la energía cuando el cable parte inicialmente del reposo y tomamos un tiempo t. Repetiremos el proceso con distintos valores de <math>F_{0}</math> para observar si se produce alguna diferencia.

Tomamos <math>F_{0}</math> igual a 1/L+0,01 Hz, 1/L-0,01 Hz y 1/L Hz.

El problema nos quedaría planteado de la siguiente manera:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,60]\\
u(0,t)=0; u_{x}(10,t)=sin(2 \pi F_{0})\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

==Cable sujeto a un aparato que envía una respuesta a las vibraciones==
En este apartado vamos a calcular el comportamiento de la energía para el valor b suponiendo ahora que el extremo derecho del cable está sujeto a un aparato que envía una respuesta a la vibración que recibe de manera que la condición de contorno ahora es, tomando como condición inicial la del apartado 2 y observando en qué cambia

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u_{x}(10,t)=b u(10,0) \\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

=Visualización del movimiento del cable=
Si el siguiente código basado en el método de diferencias finitas (mediante el método del trapecio) analizado anteriormente (en el apartado 2) se ejecuta en MATLAB (Octave no admite animaciones), se puede apreciar en la primera imagen que proporciona el programa la solución al problema enunciado en el apartado 2. Lo siguiente en salir por pantalla es una animación de la solución obtenida a lo largo del tiempo, de forma que cuando se agota el intervalo de tiempo en el que estudiamos el problema, queda la misma imagen que la primera que nos proporciona este programa; con este programa vemos como la onda va "viajando". La tercera imagen en salir por pantalla es la representación del desplazamiento vertical que sufre el cable a lo largo del tiempo, es decir, la vibración de todos los puntos de nuestro cable.
{{matlab|codigo=Aquí hay que subir CodigoAnimacion.m
}}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T19:11:30Z

Jose Manuel Corbí Garrido: /* Con un término de la serie */

{{ Beta }}
{{ TrabajoED | Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]]|

María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

=Introducción y modelización del problema=

En este artículo analizaremos el comportamiento de un cable tenso de 10 metros de longitud en una estructura civil. El cable se encontrará sujeto por los extremos y será sometido a pequeñas vibraciones. Para estudiar el problema haremos ciertas hipótesis que simplificarán el problema, y que en muchos casos dan una buena aproximación del comportamiento real que tendría un cable como el que aquí estudiamos.

En cuanto a las '''características mecánicas del cable''' haremos las siguientes '''hipótesis''':

* El cable tendrá una sección despreciable frente a su longitud.
* Será un cable homogéneo, es decir, la masa por unidad de volumen será constante.
* Dicho cable también será perfectamente flexible, es decir, solo ofrecerá resistencia a esfuerzos tangenciales (en su dirección longitudinal), siendo por ello nula su resistencia frente a esfuerzos de flexión y corte.
* El cable al someterlo a tracción es inextensible, es decir, que es lo suficientemente rígido (en dirección longitudinal) como para poder despreciar su extensibilidad, por el contrario al someterlo a compresión no ofrece resistencia y se arruga.

Por otra parte el cable será sometido a '''pequeñas vibraciones''', las cuales podemos modelizar mediante la '''ecuación de ondas''' <math>u_{tt}-u_{xx}=0</math>. Estas vibraciones, al ser pequeñas, son despreciables frente a su longitud total; preocupándonos únicamente por los desplazamientos verticales (fruto de la tracción) que sufre el cable. Estos desplazamientos verticales serán modelizados mediante la función <math>u(x,t)</math>, dependiente de dos variables; la <math>x</math> que indica la distancia al extremo de la cuerda tomado como referencia, y la <math>t</math> que define el tiempo.

Las condiciones de frontera en este caso nos indican que el cable no sufrirá desplazamiento vertical en sus extremos <math>u(0,t)=0; u(10,t)=0</math>. Por otra parte las condiciones iniciales <math>u(x,0)=0; u_{t}(x,0)=0</math> se traducen en que inicialmente el desplazamiento vertical y la velocidad del cable en todos sus puntos es cero.

El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

En los sucesivos apartados estudiaremos el comportamiento del cable en diferentes situaciones, que se traducirán en diferentes condiciones para el problema (cambios en las condiciones y en la ecuación de ondas con la que modelizamos las pequeñas vibraciones).

=Desplazamiento vertical del cable=
Sujetando el cable desde el centro y desplazándolo dos metros en la dirección perpendicular, estudiaremos el comportamiento del cable una vez lo soltemos. Estará sometido a unas pequeñas vibraciones, que originarán desplazamientos verticales a lo largo del cable que aproximaremos mediante la función <math>u(x,t)</math>. Teniendo en cuenta el desplazamiento vertical de dos metros y que el cable es soltado con una velocidad nula, el problema que vamos a estudiar es el siguiente:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

A continuación estudiaremos el desplazamiento vertical al que se ve cada punto del cable expuesto y para cada instante de tiempo. Para ello aproximaremos la ecuación diferencial que modeliza el comportamiento de dicho cable mediante los métodos de diferencias finitas y mediante series de Fourier, obteniendo gráficas que reflejan la posición de cada punto del cable en cada instante.
==Aproximación mediante diferencias finitas==
Dentro del método de diferencias finitas se nos plantea un sistema de ecuaciones diferenciales, que resolveremos aplicando diferentes métodos de los cuales luego compararemos su precisión a la hora de aproximar el comportamiento del cable y su estabilidad.

El primer método que utilizaremos para resolver el sistema de ecuaciones diferenciales será el '''método del trapecio''', el cual es '''incondicionalmente estable'''. Para implementar este método y observar el comportamiento del cable en una gráfica usaremos el siguiente código en matlab:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')
}}
Obteniendo la gráfica:

[[Archivo:ImagenCodigo1GRP9.jpeg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método del trapecio]]

Usaremos también el '''método de Euler''', el cual es un método explícito y por tanto dependemos de la ecuación en que pueda ser o no factible de ejecutar. El siguiente código matlab nos permite ejecutarlo:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler explícito

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler explícito

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for j=1:length(t)-1
u=u+dt*v+F;
v=v-dt*k*u+F;
sol(j+1,:)=[0 , u' , 0];
end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler explícito')
}}
Obteniendo el siguiente gráfico del comportamiento del cable respecto al tiempo:

[[Archivo:ImagenCodigo2GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler explícito]]

Por último método para resolver el sistema planteado mediante diferencias finitas, resolvemos numéricamente mediante el '''método de Euler''' modificado, el cual tiene más precisión que el método de Euler. El código que se muestra a continuación ejecuta dicho método y nos proporciona una gráfica del cable:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler modificado

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler modificado

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for n=1:length(t)-1
l=k*u;
% calculamos k1=f(t_n,y_n)
k1=[ v l]';
% Calculamos k2=f(t_n+h/2,y_n+1/2k1*h)
%xp=x(n)+dt; %argumento de tiempo
%solp=sol+dt*k1; %argumento y
k2=[v l]';
%Calculamos k3=f(t_n+h/2,y_n+1/2k2*h)
%mismo argumento tiempo que antes

u=u+dt/2*(k1(1,:)'+k2(1,:)')+F;
l=+k*u;
k1=[v l]';
k2=[v l]';
v=v-dt/2*(k1(2,:)'+k2(2,:)')+F;
% v=v-dt/2*(k*u+k*u)+F;
sol(n+1,:)=[0 , u' , 0];

end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler modificado')
}}
La gráfica que se obtiene es:

[[Archivo:ImagenCodigo3GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler modificado]]
===Comparativa de los métodos e interpretación===
Dados los órdenes de precisión de cada uno de los métodos para aproximar la función solución del sistema de ecuaciones planteado por el método de las diferencias finitas, el más preciso es el método del trapecio (orden 2), seguido por Euler modificado y por último Euler explicito (orden 1).

==Aproximación mediante el método de Fourier==
===Con un término de la serie===
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con un término

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con un término (Q=1)

Q=1;
for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
%fourier2(k)=(1/(k*pi))*trapz(x,velocidad.*p) /
%trapz (x,p.*p) ;
end
[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con un término')

}}
}}
[[Archivo:ImagenCodigo4GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con un término]]
{{matlab|codigo=

===Con tres términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo5.m
}}
[[Archivo:ImagenCodigo5GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con tres términos]]
===Con cinco términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo6.m
}}
[[Archivo:ImagenCodigo6GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con cinco términos]]
===Con diez términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo7.m
}}
[[Archivo:ImagenCodigo7GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con diez términos]]
===Con veinte términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo8.m
}}
[[Archivo:ImagenCodigo8.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con veinte términos]]
===Comparativa de la diferencia de términos e interpretación===
Podemos observar que a medida que aumentamos el número de términos de la serie de Fourier, la forma de la solución pasa de ser más redondeada a más recta, pareciéndose cada vez más a la solución obtenida por el método de las diferencias finitas.

==Comparativa de diferencias finitas y Fourier==
Se puede obsevar que al aumentar el numéro de términos con Fourier, la solución se va aproximando cada vez más a la obtenida mediante el método de las diferencias finitas.

=Energía del cable=

La energía del cable viene definida por la siguiente función:: <math> E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx </math>
Mediante el programa que se muestra a continuación procederemos a dibujar la gráfica de la función <math>E(t)</math>, que cuantifica la energía del cable. Previamente, y recurriendo al concepto de derivada, aproximaremos los valores que toman <math>u_{t}</math> y <math>u_{x}</math>, para así después poder calcular la integral que define <math>E(t)</math>.
{{matlab|codigo=Aquí hay que subir CodigoEnergia.m
}}

[[Archivo:ImagenCodigoEnergiaGRP9.jpg|border|Variación de la energía del cable a lo largo del tiempo]]

Como se puede apreciar en la gráfica, esta es de apariencia trigonométrica, periódica en el tiempo (el periodo es de 10 segundos).

'''¿Qué quiere decir esto?'''

Esta gráfica no es otra cosa si no la '''representación de la variación de la energía potencial en valor absoluto a lo largo del tiempo''', tomando como referencia de los potenciales un desplazamiento vertical nulo <math>u=0</math>. Para una interpretación más clara solo '''nos fijaremos en el punto situado a la mitad del cable''', en el que para <math>t=0</math> el desplazamiento vertical es de dos metros.

Inicialmente soltamos el cable desde una altura de dos metros y con velocidad nula, es decir, que en este instante la energía potencial será máxima mientras que la energía cinética será inexistente. A medida que desciende el cable, este incrementa su energía cinética, contrarrestando así la energía potencial que decrece, hasta que en <math>t=5</math> la energía potencial es nula (al igual que el desplazamiento vertical); esto quiere decir que la energía cinética es máxima en este instante. La energía potencial continúa aumentando su valor hasta que en <math>t=10</math> se vuelve a hacer casi máxima (el desplazamiento vertical del cable esta cerca de los dos metros hacia abajo) mientras que la cinética nula. Comienza nuevamente a aumentar la cinética y decrecer la potencial, de forma que para <math>t=15</math> tanto el desplazamiento vertical como la energía potencial son nulos (la cinética vuelve a ser casi máxima); vuelve a aumentar nuevamente y en <math>t=20</math> el cable vuelve a estar próximo a los dos metros de altura con una energía potencial cercana a la máxima.

La '''energía mecánica''' (suma de la cinética y la potencial) '''tiende a conservarse, pero esto no ocurre''', puesto que hay pérdidas (puede ser por el rozamiento con el aire, etc), a simple vista no se aprecia en la gráfica, pero si se ejecuta el programa y se comprueba la matriz “E” (almacena la energía en cada instante) se observa que esta va decreciendo a lo largo del tiempo.

El siguiente código auxiliar muestra el desplazamiento vertical del punto central en cada instante, a simple vista no se aprecia, pero si se ejecuta y comprueba la matriz “sol” se verá que el desplazamiento vertical es máximo al principio (2 metros), no volviendo a recuperar el cable este valor (aunque se quede próximo a él en ciertos instantes, para ello hay que observar la columna 51 de la matriz), '''el movimiento se va “amortiguando” reduciéndose lentamente el desplazamiento vertical y la energía'''. Si estudiáramos un intervalo de tiempo mayor esto se apreciaría mucho mejor, pruebe a cambiar en el siguiente código T=40 por T=500, es evidente, cada vez se aleja más de los dos metros iniciales que marcan los límites de la gráfica (observe la segunda imagen).
{{matlab|codigo=Aquí hay que subir CodigoAuxEnergia.m
}}
[[Archivo:ImagenCodigoAuxEnergia.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=40]]
[[Archivo:ImagenCodigoAuxEnergiaT500.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=500]]

=Ejemplos de aplicación de cables en estructuras civiles=
==Cable sumergido en el mar (medio viscoso)==
En este caso el cable está sumergido en un medio viscoso que produce amortiguamiento.
El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}+au_{t}=0; x∈[0,10]; t∈[0,40] \\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

Vamos a estudiar lo que ocurre con la energía en estas condiciones, dependiendo del parámetro a.Vamos a dar a a los valores de 0, 1 ,4, 10 y 100.
En el caso de que a valga cero estamos ante el mismo caso que si no estuviese sumergido, el cable permanecería en reposo.

Se ve que, al aumentar la a, la energía se disipa más rápidamente.

==Cable sujeto a una estructura sometida a vibraciones periódicas==
Ahora nos preguntamos qué ocurriría si el extremo derecho del cable está sujeto a una estructura que sufre vibraciones periódicas. La frecuencia de estas vibraciones es de <math>F_{0}</math> Herzios. Tomamos por ejemplo <math>f(t) = sin(2 \pi F_{0})</math>. Vamos a calcular el comportamiento de la energía cuando el cable parte inicialmente del reposo y tomamos un tiempo t. Repetiremos el proceso con distintos valores de <math>F_{0}</math> para observar si se produce alguna diferencia.

Tomamos <math>F_{0}</math> igual a 1/L+0,01 Hz, 1/L-0,01 Hz y 1/L Hz.

El problema nos quedaría planteado de la siguiente manera:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,60]\\
u(0,t)=0; u_{x}(10,t)=sin(2 \pi F_{0})\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

==Cable sujeto a un aparato que envía una respuesta a las vibraciones==
En este apartado vamos a calcular el comportamiento de la energía para el valor b suponiendo ahora que el extremo derecho del cable está sujeto a un aparato que envía una respuesta a la vibración que recibe de manera que la condición de contorno ahora es, tomando como condición inicial la del apartado 2 y observando en qué cambia

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u_{x}(10,t)=b u(10,0) \\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

=Visualización del movimiento del cable=
Si el siguiente código basado en el método de diferencias finitas (mediante el método del trapecio) analizado anteriormente (en el apartado 2) se ejecuta en MATLAB (Octave no admite animaciones), se puede apreciar en la primera imagen que proporciona el programa la solución al problema enunciado en el apartado 2. Lo siguiente en salir por pantalla es una animación de la solución obtenida a lo largo del tiempo, de forma que cuando se agota el intervalo de tiempo en el que estudiamos el problema, queda la misma imagen que la primera que nos proporciona este programa; con este programa vemos como la onda va "viajando". La tercera imagen en salir por pantalla es la representación del desplazamiento vertical que sufre el cable a lo largo del tiempo, es decir, la vibración de todos los puntos de nuestro cable.
{{matlab|codigo=Aquí hay que subir CodigoAnimacion.m
}}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T19:10:44Z

Jose Manuel Corbí Garrido: /* Con un término de la serie */

{{ Beta }}
{{ TrabajoED | Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]]|

María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

=Introducción y modelización del problema=

En este artículo analizaremos el comportamiento de un cable tenso de 10 metros de longitud en una estructura civil. El cable se encontrará sujeto por los extremos y será sometido a pequeñas vibraciones. Para estudiar el problema haremos ciertas hipótesis que simplificarán el problema, y que en muchos casos dan una buena aproximación del comportamiento real que tendría un cable como el que aquí estudiamos.

En cuanto a las '''características mecánicas del cable''' haremos las siguientes '''hipótesis''':

* El cable tendrá una sección despreciable frente a su longitud.
* Será un cable homogéneo, es decir, la masa por unidad de volumen será constante.
* Dicho cable también será perfectamente flexible, es decir, solo ofrecerá resistencia a esfuerzos tangenciales (en su dirección longitudinal), siendo por ello nula su resistencia frente a esfuerzos de flexión y corte.
* El cable al someterlo a tracción es inextensible, es decir, que es lo suficientemente rígido (en dirección longitudinal) como para poder despreciar su extensibilidad, por el contrario al someterlo a compresión no ofrece resistencia y se arruga.

Por otra parte el cable será sometido a '''pequeñas vibraciones''', las cuales podemos modelizar mediante la '''ecuación de ondas''' <math>u_{tt}-u_{xx}=0</math>. Estas vibraciones, al ser pequeñas, son despreciables frente a su longitud total; preocupándonos únicamente por los desplazamientos verticales (fruto de la tracción) que sufre el cable. Estos desplazamientos verticales serán modelizados mediante la función <math>u(x,t)</math>, dependiente de dos variables; la <math>x</math> que indica la distancia al extremo de la cuerda tomado como referencia, y la <math>t</math> que define el tiempo.

Las condiciones de frontera en este caso nos indican que el cable no sufrirá desplazamiento vertical en sus extremos <math>u(0,t)=0; u(10,t)=0</math>. Por otra parte las condiciones iniciales <math>u(x,0)=0; u_{t}(x,0)=0</math> se traducen en que inicialmente el desplazamiento vertical y la velocidad del cable en todos sus puntos es cero.

El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

En los sucesivos apartados estudiaremos el comportamiento del cable en diferentes situaciones, que se traducirán en diferentes condiciones para el problema (cambios en las condiciones y en la ecuación de ondas con la que modelizamos las pequeñas vibraciones).

=Desplazamiento vertical del cable=
Sujetando el cable desde el centro y desplazándolo dos metros en la dirección perpendicular, estudiaremos el comportamiento del cable una vez lo soltemos. Estará sometido a unas pequeñas vibraciones, que originarán desplazamientos verticales a lo largo del cable que aproximaremos mediante la función <math>u(x,t)</math>. Teniendo en cuenta el desplazamiento vertical de dos metros y que el cable es soltado con una velocidad nula, el problema que vamos a estudiar es el siguiente:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u(10,t)=0\\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

A continuación estudiaremos el desplazamiento vertical al que se ve cada punto del cable expuesto y para cada instante de tiempo. Para ello aproximaremos la ecuación diferencial que modeliza el comportamiento de dicho cable mediante los métodos de diferencias finitas y mediante series de Fourier, obteniendo gráficas que reflejan la posición de cada punto del cable en cada instante.
==Aproximación mediante diferencias finitas==
Dentro del método de diferencias finitas se nos plantea un sistema de ecuaciones diferenciales, que resolveremos aplicando diferentes métodos de los cuales luego compararemos su precisión a la hora de aproximar el comportamiento del cable y su estabilidad.

El primer método que utilizaremos para resolver el sistema de ecuaciones diferenciales será el '''método del trapecio''', el cual es '''incondicionalmente estable'''. Para implementar este método y observar el comportamiento del cable en una gráfica usaremos el siguiente código en matlab:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
G=[zeros(N-1,1); F];

%% Aproximación de u mediante el método del trapecio

w0=[u0 ; v0];
ww=w0;
sol(1,:)=[ 0 , w0(1:N-1)' , 0];

for j=1:length(t)-1
ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);

sol(j+1,:)=[0 , ww(1:N-1)' , 0];
end

%% Representación gráfica

figure
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método del trapecio')
}}
Obteniendo la gráfica:

[[Archivo:ImagenCodigo1GRP9.jpeg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método del trapecio]]

Usaremos también el '''método de Euler''', el cual es un método explícito y por tanto dependemos de la ecuación en que pueda ser o no factible de ejecutar. El siguiente código matlab nos permite ejecutarlo:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler explícito

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler explícito

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for j=1:length(t)-1
u=u+dt*v+F;
v=v-dt*k*u+F;
sol(j+1,:)=[0 , u' , 0];
end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler explícito')
}}
Obteniendo el siguiente gráfico del comportamiento del cable respecto al tiempo:

[[Archivo:ImagenCodigo2GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler explícito]]

Por último método para resolver el sistema planteado mediante diferencias finitas, resolvemos numéricamente mediante el '''método de Euler''' modificado, el cual tiene más precisión que el método de Euler. El código que se muestra a continuación ejecuta dicho método y nos proporciona una gráfica del cable:
{{matlab|codigo=%% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler modificado

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L-h;

dt=h;
t=0:dt:T;

k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
k=(1/h^2)*k;

F=zeros(N-1,1);

u0=(2-abs(2-xint*2/5))';
v0=0*xint';

%% Aproximación de u mediante el método de Euler modificado

sol(1,:)=[ 0 , u0' , 0];
u=u0;
v=v0;
for n=1:length(t)-1
l=k*u;
% calculamos k1=f(t_n,y_n)
k1=[ v l]';
% Calculamos k2=f(t_n+h/2,y_n+1/2k1*h)
%xp=x(n)+dt; %argumento de tiempo
%solp=sol+dt*k1; %argumento y
k2=[v l]';
%Calculamos k3=f(t_n+h/2,y_n+1/2k2*h)
%mismo argumento tiempo que antes

u=u+dt/2*(k1(1,:)'+k2(1,:)')+F;
l=+k*u;
k1=[v l]';
k2=[v l]';
v=v-dt/2*(k1(2,:)'+k2(2,:)')+F;
% v=v-dt/2*(k*u+k*u)+F;
sol(n+1,:)=[0 , u' , 0];

end

%% Representación gráfica

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante diferencias finitas con el método de Euler modificado')
}}
La gráfica que se obtiene es:

[[Archivo:ImagenCodigo3GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler modificado]]
===Comparativa de los métodos e interpretación===
Dados los órdenes de precisión de cada uno de los métodos para aproximar la función solución del sistema de ecuaciones planteado por el método de las diferencias finitas, el más preciso es el método del trapecio (orden 2), seguido por Euler modificado y por último Euler explicito (orden 1).

==Aproximación mediante el método de Fourier==
===Con un término de la serie===
{{matlab|codigo=Aquí hay que subir Codigo4.m
}}
[[Archivo:ImagenCodigo4GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con un término]]
{{matlab|codigo=
%% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con un término

%% Enunciado del problema

% u_{tt}-u_{xx}=0; x?[0,10]; t>0
% u(0,t)=0; u(10,t)=0
% u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\ u_{t}(x,0)=0\\

%% Planteamiento y discretización

clear all
L=10;
T=40;
h=0.1;
N=L/h;
dt=0.1;
M=T/dt;
x=0:h:10;
posicion=(2-abs(2-x*2/5));
velocidad=0*x;
t=0:dt:T;

%% Aproximación de u mediante Fourier con un término (Q=1)

Q=1;
for k=1:Q
p=sin(k*pi/10*x) ;
fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
%fourier2(k)=(1/(k*pi))*trapz(x,velocidad.*p) /
%trapz (x,p.*p) ;
end
[XX,TT]=meshgrid(x,t);
U=0;

for k=1:Q
U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
end

%% Representación gráfica

surf(XX,TT,U)
xlabel('Cable')
ylabel('Tiempo')
zlabel('Desplazamiento vertical')
title('u aproximado mediante Fourier con un término')

}}

===Con tres términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo5.m
}}
[[Archivo:ImagenCodigo5GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con tres términos]]
===Con cinco términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo6.m
}}
[[Archivo:ImagenCodigo6GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con cinco términos]]
===Con diez términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo7.m
}}
[[Archivo:ImagenCodigo7GRP9.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con diez términos]]
===Con veinte términos de la serie===
{{matlab|codigo=Aquí hay que subir Codigo8.m
}}
[[Archivo:ImagenCodigo8.jpg|border|Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con veinte términos]]
===Comparativa de la diferencia de términos e interpretación===
Podemos observar que a medida que aumentamos el número de términos de la serie de Fourier, la forma de la solución pasa de ser más redondeada a más recta, pareciéndose cada vez más a la solución obtenida por el método de las diferencias finitas.

==Comparativa de diferencias finitas y Fourier==
Se puede obsevar que al aumentar el numéro de términos con Fourier, la solución se va aproximando cada vez más a la obtenida mediante el método de las diferencias finitas.

=Energía del cable=

La energía del cable viene definida por la siguiente función:: <math> E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \cdot dx </math>
Mediante el programa que se muestra a continuación procederemos a dibujar la gráfica de la función <math>E(t)</math>, que cuantifica la energía del cable. Previamente, y recurriendo al concepto de derivada, aproximaremos los valores que toman <math>u_{t}</math> y <math>u_{x}</math>, para así después poder calcular la integral que define <math>E(t)</math>.
{{matlab|codigo=Aquí hay que subir CodigoEnergia.m
}}

[[Archivo:ImagenCodigoEnergiaGRP9.jpg|border|Variación de la energía del cable a lo largo del tiempo]]

Como se puede apreciar en la gráfica, esta es de apariencia trigonométrica, periódica en el tiempo (el periodo es de 10 segundos).

'''¿Qué quiere decir esto?'''

Esta gráfica no es otra cosa si no la '''representación de la variación de la energía potencial en valor absoluto a lo largo del tiempo''', tomando como referencia de los potenciales un desplazamiento vertical nulo <math>u=0</math>. Para una interpretación más clara solo '''nos fijaremos en el punto situado a la mitad del cable''', en el que para <math>t=0</math> el desplazamiento vertical es de dos metros.

Inicialmente soltamos el cable desde una altura de dos metros y con velocidad nula, es decir, que en este instante la energía potencial será máxima mientras que la energía cinética será inexistente. A medida que desciende el cable, este incrementa su energía cinética, contrarrestando así la energía potencial que decrece, hasta que en <math>t=5</math> la energía potencial es nula (al igual que el desplazamiento vertical); esto quiere decir que la energía cinética es máxima en este instante. La energía potencial continúa aumentando su valor hasta que en <math>t=10</math> se vuelve a hacer casi máxima (el desplazamiento vertical del cable esta cerca de los dos metros hacia abajo) mientras que la cinética nula. Comienza nuevamente a aumentar la cinética y decrecer la potencial, de forma que para <math>t=15</math> tanto el desplazamiento vertical como la energía potencial son nulos (la cinética vuelve a ser casi máxima); vuelve a aumentar nuevamente y en <math>t=20</math> el cable vuelve a estar próximo a los dos metros de altura con una energía potencial cercana a la máxima.

La '''energía mecánica''' (suma de la cinética y la potencial) '''tiende a conservarse, pero esto no ocurre''', puesto que hay pérdidas (puede ser por el rozamiento con el aire, etc), a simple vista no se aprecia en la gráfica, pero si se ejecuta el programa y se comprueba la matriz “E” (almacena la energía en cada instante) se observa que esta va decreciendo a lo largo del tiempo.

El siguiente código auxiliar muestra el desplazamiento vertical del punto central en cada instante, a simple vista no se aprecia, pero si se ejecuta y comprueba la matriz “sol” se verá que el desplazamiento vertical es máximo al principio (2 metros), no volviendo a recuperar el cable este valor (aunque se quede próximo a él en ciertos instantes, para ello hay que observar la columna 51 de la matriz), '''el movimiento se va “amortiguando” reduciéndose lentamente el desplazamiento vertical y la energía'''. Si estudiáramos un intervalo de tiempo mayor esto se apreciaría mucho mejor, pruebe a cambiar en el siguiente código T=40 por T=500, es evidente, cada vez se aleja más de los dos metros iniciales que marcan los límites de la gráfica (observe la segunda imagen).
{{matlab|codigo=Aquí hay que subir CodigoAuxEnergia.m
}}
[[Archivo:ImagenCodigoAuxEnergia.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=40]]
[[Archivo:ImagenCodigoAuxEnergiaT500.jpg|border|Desplazamiento vertical del cable a lo largo del tiempo en el punto medio T=500]]

=Ejemplos de aplicación de cables en estructuras civiles=
==Cable sumergido en el mar (medio viscoso)==
En este caso el cable está sumergido en un medio viscoso que produce amortiguamiento.
El problema que se nos plantea es el siguiente:
<math>
\begin{cases}
u_{tt}-u_{xx}+au_{t}=0; x∈[0,10]; t∈[0,40] \\
u(0,t)=0; u(10,t)=0\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

Vamos a estudiar lo que ocurre con la energía en estas condiciones, dependiendo del parámetro a.Vamos a dar a a los valores de 0, 1 ,4, 10 y 100.
En el caso de que a valga cero estamos ante el mismo caso que si no estuviese sumergido, el cable permanecería en reposo.

Se ve que, al aumentar la a, la energía se disipa más rápidamente.

==Cable sujeto a una estructura sometida a vibraciones periódicas==
Ahora nos preguntamos qué ocurriría si el extremo derecho del cable está sujeto a una estructura que sufre vibraciones periódicas. La frecuencia de estas vibraciones es de <math>F_{0}</math> Herzios. Tomamos por ejemplo <math>f(t) = sin(2 \pi F_{0})</math>. Vamos a calcular el comportamiento de la energía cuando el cable parte inicialmente del reposo y tomamos un tiempo t. Repetiremos el proceso con distintos valores de <math>F_{0}</math> para observar si se produce alguna diferencia.

Tomamos <math>F_{0}</math> igual a 1/L+0,01 Hz, 1/L-0,01 Hz y 1/L Hz.

El problema nos quedaría planteado de la siguiente manera:

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,60]\\
u(0,t)=0; u_{x}(10,t)=sin(2 \pi F_{0})\\
u(x,0)=0; u_{t}(x,0)=0\\
\end{cases}
</math>

==Cable sujeto a un aparato que envía una respuesta a las vibraciones==
En este apartado vamos a calcular el comportamiento de la energía para el valor b suponiendo ahora que el extremo derecho del cable está sujeto a un aparato que envía una respuesta a la vibración que recibe de manera que la condición de contorno ahora es, tomando como condición inicial la del apartado 2 y observando en qué cambia

<math>
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\
u(0,t)=0; u_{x}(10,t)=b u(10,0) \\
u(x,0)=2-|2-\frac{2}{5}x|; u(5,0)=2\\
u_{t}(x,0)=0\\
\end{cases}
</math>

=Visualización del movimiento del cable=
Si el siguiente código basado en el método de diferencias finitas (mediante el método del trapecio) analizado anteriormente (en el apartado 2) se ejecuta en MATLAB (Octave no admite animaciones), se puede apreciar en la primera imagen que proporciona el programa la solución al problema enunciado en el apartado 2. Lo siguiente en salir por pantalla es una animación de la solución obtenida a lo largo del tiempo, de forma que cuando se agota el intervalo de tiempo en el que estudiamos el problema, queda la misma imagen que la primera que nos proporciona este programa; con este programa vemos como la onda va "viajando". La tercera imagen en salir por pantalla es la representación del desplazamiento vertical que sufre el cable a lo largo del tiempo, es decir, la vibración de todos los puntos de nuestro cable.
{{matlab|codigo=Aquí hay que subir CodigoAnimacion.m
}}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T18:54:53Z

Jose Manuel Corbí Garrido: /* Comparativa de la diferencia de términos e interpretación */

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T18:50:28Z

Jose Manuel Corbí Garrido: /* Cable sujeto a una estructura sometida a vibraciones periódicas */

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T18:50:10Z

Jose Manuel Corbí Garrido: /* Cable sumergido en el mar (medio viscoso) */

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T18:49:40Z

Jose Manuel Corbí Garrido: /* Cable sumergido en el mar (medio viscoso) */

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T18:48:09Z

Jose Manuel Corbí Garrido: /* Cable sujeto a un aparato que envía una respuesta a las vibraciones */

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T18:47:55Z

Jose Manuel Corbí Garrido: /* Cable sujeto a un aparato que envía una respuesta a las vibraciones */

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T18:30:08Z

Jose Manuel Corbí Garrido: /* Cable sujeto a una estructura sometida a vibraciones periódicas */

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-19T18:29:49Z

Jose Manuel Corbí Garrido: /* Cable sujeto a una estructura sometida a vibraciones periódicas */

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-18T19:26:26Z

Jose Manuel Corbí Garrido: /* Cable sumergido en el mar (medio viscoso) */

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-18T19:01:49Z

Jose Manuel Corbí Garrido: /* Cable sumergido en el mar (medio viscoso) */

Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)

2014-05-18T12:50:41Z

Jose Manuel Corbí Garrido: /* Comparativa de diferencias finitas y Fourier */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-05T20:27:38Z

Jose Manuel Corbí Garrido: /* Cálculo y representación de u_1 */

{{ TrabajoED | Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

{{ Beta }}
=Introducción y enunciado del problema=
En este artículo, se desarrollará el estudio de la capa límite de un fluido laminar sobre una placa plana semi-infinita mediante métodos numéricos programados en lenguaje M. De esta forma el lector podrá reproducir (si así lo desea) los algoritmos y códigos aquí escritos en MATLAB u otras alternativas de software libre como Octave.

[[Archivo:Grp9BIMAG1.png|miniaturadeimagen|Fluido laminar por encima de una placa plana]]

El tratamiento de este problema mediante métodos numéricos requiere una definición precisa de las condiciones iniciales y de contorno. En nuestro caso definiremos la semi-placa en el plano cartesiano como la semirecta <math>y=0</math>, con <math>x∈(0,∞)</math>; y estudiaremos lo que pasa por encima de ella <math>(x>0, y>0)</math>. El fluido laminar que se encuentra con la placa, antes de llegar a ella tiene una velocidad constante de <math>\overrightarrow{u} = u_0\cdot\overrightarrow{i}, u_0 = 2 </math>, iniciándose una fase de transición en su encuentro con la placa (la velocidad es nula en los puntos de la placa), y lejos del punto de encuentro, una vez finalizada la transición, volviendo a fluir a la velocidad constante inicial.
El campo de velocidades del fluido se define como::
<math>\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})</math>
donde <math>\psi</math> es la corriente del fluido::
<math>\psi(x,y) = \sqrt[]{ \nu \cdot\ u_0 \cdot x} f(\eta)</math> y <math>ν</math> es su viscosidad (que supondremos <math>ν = 1</math>), <math> f(\eta) </math> con <math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math> satisface la '''ecuación de Blasius'''::

<math>f''' + \frac{1}{2}f f'' =0</math>

En los siguientes apartados resolveremos numéricamente esta ecuación diferencial ordinaria no lineal, para unas condiciones iniciales::
<math>f(0)=f(0)'=0; f(0)''=k</math> <math> f'\rightarrow 1</math> cuando <math>\eta\rightarrow \infty</math>
tomando <math>k</math> valores <math> k \in \mbox{(0,1;1)}</math> y <math>dk=0,01</math>
A continuación estudiaremos la resolución de la ecuación mediante diferentes métodos numéricos, y estudiaremos el <math>k</math> "óptimo" para las condiciones iniciales enunciadas.

=Resolución numérica de la ecuación de Blasius=
Para poder resolver esta ecuación independientemente del posterior procedimiento numérico que vayamos a aplicar, primero hemos tenido que convertirla en un sistema, aplicando el método analítico que detallamos a continuación::

<math>y_{1}=f(x)</math>:

<math>y_{2}=f'(x)</math>:

<math>y_{3}=f''(x)</math>

Así <math>y_{3}'=f'''(x)</math>, si despejamos de nuestra ecuación diferencial obtenemos que::

<math>y_{2}=y_{1}'</math>:

<math>y_{3}=y_{2}'</math>:

<math>y_{3}'=-\frac{1}{2}y_{1}y_{3}</math>

Como resultado de esta transformación obtenemos un sistema para el cual podemos obtener su solución numéricamente.

==Mediante el método de Euler modificado==
A continuación se muestra el código en lenguaje matlab que proporciona las soluciones numéricas a la ecuación de Blasius para los diferentes valores de k (recordemos que tomaba valores de 0,1 hasta 1 con un paso de 0,01) en un intervalo <math> \eta \in \mbox{(0,20)}</math>. Hemos tomado una <math>h=0,5</math>.

Implementando el método de Euler en este sistema:
{{matlab|codigo=
%% Método de Euler modificado

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
%% Datos del problema y discretización
x0=0; % las x son las eta de las que depende f
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta, k y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada x
%% Iteraciones y(x_n)-->y(x_n+1)
for n=1:N
% calculamos k1=f(x_n,y_n)
k1=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
% Calculamos k2=f(x_n+h/2,y_n+1/2k1*h)
xp=x(n)+h; %argumento de x
fp=ff+h*k1; %argumento y
k2=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';

ff=ff+h/2*(k1+k2);
f(:,n+1)=ff;

end
laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solución
f2=laf2paraescogerk;
figure(1)
clf
plot(x,f)
legend('f(x)','f´(x)','f´´(x)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f, f´,f´´ en intervalo de eta, siendo k el que mejor aproxima f´(20) a 1')
corte=1;
figure (2)
hold on
plot(k,laf2paraescogerk)
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')
hold off
%Para calcular la k de forma exacta utilizamos este bucle que nos dice con qué k la diferencia en valor absoluto entre la f´(n) y uno es menor
for n=1:K
l=n;
p=abs(1-f2(n));
if min(abs(1-f2))==p

break
end
end
kfinal=k(l)
}}

[[Archivo:Grp9BIMAG2eulermodif1.jpg|miniaturadeimagen|f,f´,f´´ con k=0,33 para intervalo de eta de 0 a 20]]
[[Archivo:Grp9BIMAG2eulermodif2.jpg|border|f´(20) para los distintos valores de k]]

Como se puede apreciar en el gráfico, según este método la k que mejor se ajusta a las condiciones iniciales es <math>k=0.330</math>

==Mediante el método de Runge-Kutta de cuarto orden==
La siguiente secuencia de código da solución a la ecuación de Blasius mediante el método de Runge-Kutta, siendo aplicado al sistema de ecuaciones diferenciales obtenido en el primer apartado mediante procedimientos analíticos:

{{matlab|codigo=
% Resolver sistema con un método RK4
%% Método de Runge-Kutta orden 4

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
% Datos del problema y discretización
% La variable eta recibe el nombre x
x0=0;
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada tiempo

for n=1:N
% calculamos k1=f(x_n,f_n)
k1=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
% Calculamos k2=f(x_n+h/2,f_n+1/2k1*h)
xp=x(n)+h/2; %argumento de eta
fp=ff+1/2*h*k1; %argumento f
k2=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
%Calculamos k3=f(x_n+h/2,f_n+1/2k2*h)
%mismo argumento eta que antes
fp=ff+1/2*h*k2;%argumento y
k3=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
%calculamos k4=f(x_n+h/2,f_n+1/2k3*h)
xp=x(n)+h;%argumento eta
fp=ff+h*k3; %argumento f
k4=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
ff=ff+h/6*(k1+2*k2+2*k3+k4);
f(:,n+1)=ff;
end
laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solucion
f2=laf2paraescogerk;
figure(1)
clf
plot(x,f)
legend('f1(x)','f2(x)','f3(x)')
corte=1;
figure (2)
plot(k,laf2paraescogerk)
hold on
plot(k,corte,'red')
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')

for n=1:K
l=n;
p=abs(1-f2(n));
if min(abs(1-f2))==p

break
end
end

kfinal=k(l)
}}
[[Archivo:Grafipartado2d.jpeg|border]]

Según este método la k que mejor se ajusta a las condiciones iniciales es k=0.330, como podemos apreciar en el gráfico.

==Mediante el método de Euler==
Al igual que en los otros apartados resolvemos el sistema asociado a la ecuación de Blasius numéricamente, en este caso mediante el método de Euler:
{{matlab|codigo=
%% Método de Euler

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
% Datos del problema y discretización
% x es la variable eta del problema
x0=0;
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada valor de eta
%% Iteraciones f(x_n)-->f(x_n+1)
for n=1:N

y=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
ff=ff+h*y;
f(:,n+1)=ff;

end
laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solucion
f2=laf2paraescogerk;
clf
corte=1;
figure (1)
plot(k,laf2paraescogerk)
hold on
plot(k,corte,'red')
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')

for n=1:K
l=n;
p=abs(1-f2(n));
if min(abs(1-f2))==p

break
end
end

kfinal=k(l)
figure(2)
%ya con la k=0,33

x0=0;
xf=20;
k=0.32
f0=[0 0 k]';
%% Datos discretizacion
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
f=zeros(3,N+1);
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada eta
%% Iteraciones f(x_n)-->f(x_n+1)
for n=1:N

y=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
ff=ff+h*y;
f(:,n+1)=ff;

end
plot(x,f)
legend('f(x)','f´(x)','f´´(x)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f, f´,f´´ en intervalo de eta, con Euler f´(20) a 1, ya con la k escogida')
}}

[[Archivo:Grafipartado2c.jpeg|border]]

Como se puede apreciar en el gráfico, según este método la k que mejor se ajusta a las condiciones iniciales es k=0.320

==Comparativa entre los tres métodos, elección de k==
La diferencia entre los tres métodos numéricos empleados para resolver la ecuacion de Blasius (Euler modificado, Runge-Kutta y Euler), es muy pequeña. La forma de las gráficas es completamente similar y el valor que se obtiene para que f'(20) sea lo más proxima a uno es en el caso de Euler modificado y Runge-Kutta igual a 0.33 y en Euler es 0.32. Por la tanto la diferencia entre los tres métodos es pequeña.

[[Archivo:Grp9BIMAG2eulermodif2.jpg|border|f´(20) para los distintos valores de k]]

[[Archivo:Grafipartado2d.jpeg|border|Con RK4]]

[[Archivo:Grafipartado2c.jpeg|border|Con Euler]]

'''De hecho a simple vista no se aprecia la diferencia.'''

=Representación e interpretación de <math>f’(\eta)</math>=
Vamos a estudiar a partir de qué valor de <math>\eta</math> la gráfica difiere 0,01 o menos de 1.

Para esto se puede escoger cualquiera de los métodos numéricos anteriormente utilizados. Y después, para que nos muestre la gráfica de <math>f’(\eta)</math>, se escriben las siguientes líneas de código:

{{matlab|codigo=
plot(eta,f(2,:),'g')
legend('f2(eta)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f´ para k=0,33')
}}:

[[Archivo:grafipartado3.jpg|grafipartado3.jpg]]

Observamos como tiende rápidamente al valor 1, por lo que confirma que la k escogida en los apartados anteriores es correcta.

Ahora para encontrar a partir de qué valor de <math>\eta</math> la gráfica difiere 0,01 o menos de 1; introducimos la siguiente secuencia de código:

{{matlab|codigo=
figure(2)
gg=abs(f(2,:)-1);
for n=1:N
l=n;
if gg(n)<0.01
break
end
end
nu0=l-1
loquequermos=eta0*h

plot(eta0,gg)

}}
[[Archivo:grafipartado3b.jpg|grafipartado3b.jpg]]

Visualmente vemos que el valor buscado está más o menos entre 5 y 6. Con la condición del código de que la diferencia en valor absoluto sea menor que 0.01 llegamos a que el valor que buscamos concretamente es 5,2.

=Campo de velocidades del fluido=
Definimos el campo de velocidades del fluido <math>\vec{u}</math> como::

<math>\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})</math>

==Cálculo y representación de <math>u_1</math>==
A continuación vamos a calcular analíticamente la componente de <math> u_1 </math> de la velocidad del fluido en función de <math>f( \eta)</math>

Para ello, desarrollamos la componente hasta la obtención de su expresión final::

<math>u_1=\frac{\partial \psi}{\partial y}=\frac{\partial }{\partial y}(\sqrt[]{\nu u_0 x} f(\eta))= \sqrt[]{\nu u_0 x} \frac{\partial f(\eta)}{\partial y}=\sqrt[]{\nu u_0 x} \frac{\partial f(\eta)}{\partial \eta} \frac{\partial \eta}{\partial y}</math>:

Definimos la función <math>\eta(x,y)</math> como::

<math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math>:

<math>u_1=\sqrt[]{\nu u_0 x} \sqrt[]{\frac{u_0}{\nu x}} f’(\eta)=u_0 f’(\eta)</math>

Resolviendo esta ecuación mediante el método de Rungge-Kutta para un <math>y \in (0,3)</math> y unos valores de <math>x \in \{0.05, 0.2, 0.4, 0.6, 0.8\}</math> mediante el siguiente código, que nos permite calcular las gráficas para cada una de las <math>x_{i}</math> particulares:

{{matlab|codigo=

% Componente u1 de la velocidad
% Introducimos las variables
clear all
u0=2;
nu=1;
xk=[0.05 0.2 0.4 0.6 0.8];
y0=0;
yf=3;
h=0.05;
y=y0:h:yf;
N=(yf-y0)/h;
% Iteraciones
for i=1:5
eta(i,:)=y*sqrt(u0/(nu*xk(i)));
% H es el paso para resolver el sistema
H=0.05*sqrt(u0/(nu*xk(i)));
k=0.33;
f0=[0 0 k]';

ff=f0;
for n=1:N
A=[0 1 0;0 0 1;(-ff(3)/2) 0 0];
j=ff+H*A*ff;
ff=ff+(H/2)*(A*ff+A*j);
f(:,n+1)=[ff(1);ff(2);ff(3)];
end
G(i,:)=f(2,:);
end
u1=u0*G;

%% Dibujamos la solución
hold on
plot(y,u1)

legend('para x=0.05','para x=0.2','para x=0.4','para x=0.6','para x=0.8')
xlabel('y')
ylabel('u1')
title ('Componente u1 de la velocidad')
hold off

}}
[[Archivo:Gráfica del apartado 4.JPG|Gráfica del apartado 4.JPG]]

==Conclusión==
Como vemos en la gráfica, podemos deducir que el fluido cuando se mueve sobre el eje de abscisas, necesita una altura mayor para poder tener una velocidad límite más estabilizada. Observamos que la primera línea es muy próxima a una recta pero que las siguiente según avanzamos en las x se van acercando a una figura más parabólica. Estas curvas representan la transición de la velocidad del fluido desde la capa (donde es nula) hasta que alcanza u0 aumenta a medida que avanzamos por la placa.

=Capa límite, transición entre la velocidad nula al comienzo de la placa y lejos de la misma=
Hemos comprobado que la altura (<math>y</math>) a la que el fluido alcanza la velocidad constante (<math>u_0=2</math>) depende de la distancia recorrida sobre la placa (<math>x</math>). Ahora vamos a dibujar la función <math>g(x)</math>, donde <math>g(x)</math> es el valor <math>y</math> para el cual <math>\eta=\eta_0</math>, para observar como es ese cambio de <math>y</math> en función de la <math>x</math>.
La relación entre estas variables viene dada por la expresión: <math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math>
Despejamos la <math>y</math> y suponemos <math>u_0=2</math>, <math>\nu=1</math> y <math>\eta=\eta_0=5.2</math>
<center><math>y=\eta_0 \sqrt[]{ \frac{x}{2}}=g(x)</math></center>
Representamos esta función en Matlab.
{{matlab|codigo=
% Dibujamos la función g(x), y en función eta0(cte.) y de x.
% Introducimos variables
eta0=5.2;
u0=2;
v=1;
x=0:0.05:10;
% Calculamos y
y=eta0*sqrt(v*x/u0);
% Dibujamos la gráfica
plot(x,y)
legend('g(x)')
xlabel('x')
ylabel('y')
title('Capa límite de transición de la velocidad')
}}
[[Archivo:Gráficaapartado5.JPG|Gráficaapartado5.JPG]]

Interpretamos la función <math>g(x)</math> como el borde de la capa límite que hace la transición entre la velocidad
nula en la placa y la velocidad lejos de la misma del fluido. Se confirma lo que ya podíamos ver en el apartado anterior. A medida que el fluido se aleja del origen de la placa (<math>x=0</math>), es decir, a medida que aumenta el valor de la variable <math>x</math>, el fluido alcanza la velocidad inicial constante a una mayor altura (<math>y</math>).
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-05T20:08:37Z

Jose Manuel Corbí Garrido: /* Conclusión */

{{ TrabajoED | Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

{{ Beta }}
=Introducción y enunciado del problema=
En este artículo, se desarrollará el estudio de la capa límite de un fluido laminar sobre una placa plana semi-infinita mediante métodos numéricos programados en lenguaje M. De esta forma el lector podrá reproducir (si así lo desea) los algoritmos y códigos aquí escritos en MATLAB u otras alternativas de software libre como Octave.

[[Archivo:Grp9BIMAG1.png|miniaturadeimagen|Fluido laminar por encima de una placa plana]]

El tratamiento de este problema mediante métodos numéricos requiere una definición precisa de las condiciones iniciales y de contorno. En nuestro caso definiremos la semi-placa en el plano cartesiano como la semirecta <math>y=0</math>, con <math>x∈(0,∞)</math>; y estudiaremos lo que pasa por encima de ella <math>(x>0, y>0)</math>. El fluido laminar que se encuentra con la placa, antes de llegar a ella tiene una velocidad constante de <math>\overrightarrow{u} = u_0\cdot\overrightarrow{i}, u_0 = 2 </math>, iniciándose una fase de transición en su encuentro con la placa (la velocidad es nula en los puntos de la placa), y lejos del punto de encuentro, una vez finalizada la transición, volviendo a fluir a la velocidad constante inicial.
El campo de velocidades del fluido se define como::
<math>\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})</math>
donde <math>\psi</math> es la corriente del fluido::
<math>\psi(x,y) = \sqrt[]{ \nu \cdot\ u_0 \cdot x} f(\eta)</math> y <math>ν</math> es su viscosidad (que supondremos <math>ν = 1</math>), <math> f(\eta) </math> con <math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math> satisface la '''ecuación de Blasius'''::

<math>f''' + \frac{1}{2}f f'' =0</math>

En los siguientes apartados resolveremos numéricamente esta ecuación diferencial ordinaria no lineal, para unas condiciones iniciales::
<math>f(0)=f(0)'=0; f(0)''=k</math> <math> f'\rightarrow 1</math> cuando <math>\eta\rightarrow \infty</math>
tomando <math>k</math> valores <math> k \in \mbox{(0,1;1)}</math> y <math>dk=0,01</math>
A continuación estudiaremos la resolución de la ecuación mediante diferentes métodos numéricos, y estudiaremos el <math>k</math> "óptimo" para las condiciones iniciales enunciadas.

=Resolución numérica de la ecuación de Blasius=
Para poder resolver esta ecuación independientemente del posterior procedimiento numérico que vayamos a aplicar, primero hemos tenido que convertirla en un sistema, aplicando el método analítico que detallamos a continuación::

<math>y_{1}=f(x)</math>:

<math>y_{2}=f'(x)</math>:

<math>y_{3}=f''(x)</math>

Así <math>y_{3}'=f'''(x)</math>, si despejamos de nuestra ecuación diferencial obtenemos que::

<math>y_{2}=y_{1}'</math>:

<math>y_{3}=y_{2}'</math>:

<math>y_{3}'=-\frac{1}{2}y_{1}y_{3}</math>

Como resultado de esta transformación obtenemos un sistema para el cual podemos obtener su solución numéricamente.

==Mediante el método de Euler modificado==
A continuación se muestra el código en lenguaje matlab que proporciona las soluciones numéricas a la ecuación de Blasius para los diferentes valores de k (recordemos que tomaba valores de 0,1 hasta 1 con un paso de 0,01) en un intervalo <math> \eta \in \mbox{(0,20)}</math>. Hemos tomado una <math>h=0,5</math>.

Implementando el método de Euler en este sistema:
{{matlab|codigo=
%% Método de Euler modificado

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
%% Datos del problema y discretización
x0=0; % las x son las eta de las que depende f
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta, k y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada x
%% Iteraciones y(x_n)-->y(x_n+1)
for n=1:N
% calculamos k1=f(x_n,y_n)
k1=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
% Calculamos k2=f(x_n+h/2,y_n+1/2k1*h)
xp=x(n)+h; %argumento de x
fp=ff+h*k1; %argumento y
k2=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';

ff=ff+h/2*(k1+k2);
f(:,n+1)=ff;

end
laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solución
f2=laf2paraescogerk;
figure(1)
clf
plot(x,f)
legend('f(x)','f´(x)','f´´(x)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f, f´,f´´ en intervalo de eta, siendo k el que mejor aproxima f´(20) a 1')
corte=1;
figure (2)
hold on
plot(k,laf2paraescogerk)
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')
hold off
%Para calcular la k de forma exacta utilizamos este bucle que nos dice con qué k la diferencia en valor absoluto entre la f´(n) y uno es menor
for n=1:K
l=n;
p=abs(1-f2(n));
if min(abs(1-f2))==p

break
end
end
kfinal=k(l)
}}

[[Archivo:Grp9BIMAG2eulermodif1.jpg|miniaturadeimagen|f,f´,f´´ con k=0,33 para intervalo de eta de 0 a 20]]
[[Archivo:Grp9BIMAG2eulermodif2.jpg|border|f´(20) para los distintos valores de k]]

Como se puede apreciar en el gráfico, según este método la k que mejor se ajusta a las condiciones iniciales es <math>k=0.330</math>

==Mediante el método de Runge-Kutta de cuarto orden==
La siguiente secuencia de código da solución a la ecuación de Blasius mediante el método de Runge-Kutta, siendo aplicado al sistema de ecuaciones diferenciales obtenido en el primer apartado mediante procedimientos analíticos:

{{matlab|codigo=
% Resolver sistema con un método RK4
%% Método de Runge-Kutta orden 4

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
% Datos del problema y discretización
% La variable eta recibe el nombre x
x0=0;
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada tiempo

for n=1:N
% calculamos k1=f(x_n,f_n)
k1=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
% Calculamos k2=f(x_n+h/2,f_n+1/2k1*h)
xp=x(n)+h/2; %argumento de eta
fp=ff+1/2*h*k1; %argumento f
k2=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
%Calculamos k3=f(x_n+h/2,f_n+1/2k2*h)
%mismo argumento eta que antes
fp=ff+1/2*h*k2;%argumento y
k3=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
%calculamos k4=f(x_n+h/2,f_n+1/2k3*h)
xp=x(n)+h;%argumento eta
fp=ff+h*k3; %argumento f
k4=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
ff=ff+h/6*(k1+2*k2+2*k3+k4);
f(:,n+1)=ff;
end
laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solucion
f2=laf2paraescogerk;
figure(1)
clf
plot(x,f)
legend('f1(x)','f2(x)','f3(x)')
corte=1;
figure (2)
plot(k,laf2paraescogerk)
hold on
plot(k,corte,'red')
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')

for n=1:K
l=n;
p=abs(1-f2(n));
if min(abs(1-f2))==p

break
end
end

kfinal=k(l)
}}
[[Archivo:Grafipartado2d.jpeg|border]]

Según este método la k que mejor se ajusta a las condiciones iniciales es k=0.330, como podemos apreciar en el gráfico.

==Mediante el método de Euler==
Al igual que en los otros apartados resolvemos el sistema asociado a la ecuación de Blasius numéricamente, en este caso mediante el método de Euler:
{{matlab|codigo=
%% Método de Euler

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
% Datos del problema y discretización
% x es la variable eta del problema
x0=0;
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada valor de eta
%% Iteraciones f(x_n)-->f(x_n+1)
for n=1:N

y=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
ff=ff+h*y;
f(:,n+1)=ff;

end
laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solucion
f2=laf2paraescogerk;
clf
corte=1;
figure (1)
plot(k,laf2paraescogerk)
hold on
plot(k,corte,'red')
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')

for n=1:K
l=n;
p=abs(1-f2(n));
if min(abs(1-f2))==p

break
end
end

kfinal=k(l)
figure(2)
%ya con la k=0,33

x0=0;
xf=20;
k=0.32
f0=[0 0 k]';
%% Datos discretizacion
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
f=zeros(3,N+1);
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada eta
%% Iteraciones f(x_n)-->f(x_n+1)
for n=1:N

y=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
ff=ff+h*y;
f(:,n+1)=ff;

end
plot(x,f)
legend('f(x)','f´(x)','f´´(x)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f, f´,f´´ en intervalo de eta, con Euler f´(20) a 1, ya con la k escogida')
}}

[[Archivo:Grafipartado2c.jpeg|border]]

Como se puede apreciar en el gráfico, según este método la k que mejor se ajusta a las condiciones iniciales es k=0.320

==Comparativa entre los tres métodos, elección de k==
La diferencia entre los tres métodos numéricos empleados para resolver la ecuacion de Blasius (Euler modificado, Runge-Kutta y Euler), es muy pequeña. La forma de las gráficas es completamente similar y el valor que se obtiene para que f'(20) sea lo más proxima a uno es en el caso de Euler modificado y Runge-Kutta igual a 0.33 y en Euler es 0.32. Por la tanto la diferencia entre los tres métodos es pequeña.

[[Archivo:Grp9BIMAG2eulermodif2.jpg|border|f´(20) para los distintos valores de k]]

[[Archivo:Grafipartado2d.jpeg|border|Con RK4]]

[[Archivo:Grafipartado2c.jpeg|border|Con Euler]]

'''De hecho a simple vista no se aprecia la diferencia.'''

=Representación e interpretación de <math>f’(\eta)</math>=
Vamos a estudiar a partir de qué valor de <math>\eta</math> la gráfica difiere 0,01 o menos de 1.

Para esto se puede escoger cualquiera de los métodos numéricos anteriormente utilizados. Y después, para que nos muestre la gráfica de <math>f’(\eta)</math>, se escriben las siguientes líneas de código:

{{matlab|codigo=
plot(eta,f(2,:),'g')
legend('f2(eta)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f´ para k=0,33')
}}:

[[Archivo:grafipartado3.jpg|grafipartado3.jpg]]

Observamos como tiende rápidamente al valor 1, por lo que confirma que la k escogida en los apartados anteriores es correcta.

Ahora para encontrar a partir de qué valor de <math>\eta</math> la gráfica difiere 0,01 o menos de 1; introducimos la siguiente secuencia de código:

{{matlab|codigo=
figure(2)
gg=abs(f(2,:)-1);
for n=1:N
l=n;
if gg(n)<0.01
break
end
end
nu0=l-1
loquequermos=eta0*h

plot(eta0,gg)

}}
[[Archivo:grafipartado3b.jpg|grafipartado3b.jpg]]

Visualmente vemos que el valor buscado está más o menos entre 5 y 6. Con la condición del código de que la diferencia en valor absoluto sea menor que 0.01 llegamos a que el valor que buscamos concretamente es 5,2.

=Campo de velocidades del fluido=
Definimos el campo de velocidades del fluido <math>\vec{u}</math> como::

<math>\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})</math>

==Cálculo y representación de <math>u_1</math>==
A continuación vamos a calcular analíticamente la componente de <math> u_1 </math> de la velocidad del fluido en función de <math>f( \eta)</math>

Para ello, desarrollamos la componente hasta la obtención de su expresión final::

<math>u_1=\frac{\partial \psi}{\partial y}=\frac{\partial }{\partial y}(\sqrt[]{\nu u_0 x} f(\eta))= \sqrt[]{\nu u_0 x} \frac{\partial f(\eta)}{\partial y}=\sqrt[]{\nu u_0 x} \frac{\partial f(\eta)}{\partial \eta} \frac{\partial \eta}{\partial y}</math>:

Definimos la función <math>\eta(x,y)</math> como::

<math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math>:

<math>u_1=\sqrt[]{\nu u_0 x} \sqrt[]{\frac{u_0}{\nu x}} f’(\eta)=u_0 f’(\eta)</math>

Resolviendo esta ecuación mediante el método de Rungge-Kutta para un <math>y \in (0,3)</math> y unos valores de <math>x \in \{0.05, 0.2, 0.4, 0.6, 0.8\}</math> mediante el siguiente código, que nos permite calcular las gráficas para cada una de las <math>x_{i}</math> particulares:

{{matlab|codigo=

% Componente u1 de la velocidad
% Introducimos las variables
clear all
u0=2;
nu=1;
xk=[0.05 0.2 0.4 0.6 0.8];
y0=0;
yf=3;
h=0.05;
y=y0:h:yf;
N=(yf-y0)/h;
% Iteraciones
for i=1:5
eta(i,:)=y*sqrt(u0/(nu*xk(i)));
% H es el paso para resolver el sistema
H=0.05*sqrt(u0/(nu*xk(i)));
k=0.33;
f0=[0 0 k]';
f=[f0(1);f0(2);f0(3)];
ff=f0;
for n=1:N
A=[0 1 0;0 0 1;(-ff(3)/2) 0 0];
j=ff+H*A*ff;
ff=ff+(H/2)*(A*ff+A*j);
f(:,n+1)=[ff(1);ff(2);ff(3)];
end
G(i,:)=f(2,:);
end
u1=u0*G;

%% Dibujamos la solución
hold on
plot(y,u1(1,:),'m')
plot(y,u1(2,:),'b')
plot(y,u1(3,:),'r')
plot(y,u1(4,:),'g')
plot(y,u1(5,:),'c')
legend('para x=0.05','para x=0.2','para x=0.4','para x=0.6','para x=0.8')
xlabel('y')
ylabel('u1')
title ('Componente u1 de la velocidad')
hold off

}}
[[Archivo:Gráfica del apartado 4.JPG|Gráfica del apartado 4.JPG]]

==Conclusión==
Como vemos en la gráfica, podemos deducir que el fluido cuando se mueve sobre el eje de abscisas, necesita una altura mayor para poder tener una velocidad límite más estabilizada. Observamos que la primera línea es muy próxima a una recta pero que las siguiente según avanzamos en las x se van acercando a una figura más parabólica. Estas curvas representan la transición de la velocidad del fluido desde la capa (donde es nula) hasta que alcanza u0 aumenta a medida que avanzamos por la placa.

=Capa límite, transición entre la velocidad nula al comienzo de la placa y lejos de la misma=
Hemos comprobado que la altura (<math>y</math>) a la que el fluido alcanza la velocidad constante (<math>u_0=2</math>) depende de la distancia recorrida sobre la placa (<math>x</math>). Ahora vamos a dibujar la función <math>g(x)</math>, donde <math>g(x)</math> es el valor <math>y</math> para el cual <math>\eta=\eta_0</math>, para observar como es ese cambio de <math>y</math> en función de la <math>x</math>.
La relación entre estas variables viene dada por la expresión: <math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math>
Despejamos la <math>y</math> y suponemos <math>u_0=2</math>, <math>\nu=1</math> y <math>\eta=\eta_0=5.2</math>
<center><math>y=\eta_0 \sqrt[]{ \frac{x}{2}}=g(x)</math></center>
Representamos esta función en Matlab.
{{matlab|codigo=
% Dibujamos la función g(x), y en función eta0(cte.) y de x.
% Introducimos variables
eta0=5.2;
u0=2;
v=1;
x=0:0.05:10;
% Calculamos y
y=eta0*sqrt(v*x/u0);
% Dibujamos la gráfica
plot(x,y)
legend('g(x)')
xlabel('x')
ylabel('y')
title('Capa límite de transición de la velocidad')
}}
[[Archivo:Gráficaapartado5.JPG|Gráficaapartado5.JPG]]

Interpretamos la función <math>g(x)</math> como el borde de la capa límite que hace la transición entre la velocidad
nula en la placa y la velocidad lejos de la misma del fluido. Se confirma lo que ya podíamos ver en el apartado anterior. A medida que el fluido se aleja del origen de la placa (<math>x=0</math>), es decir, a medida que aumenta el valor de la variable <math>x</math>, el fluido alcanza la velocidad inicial constante a una mayor altura (<math>y</math>).
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-05T20:04:09Z

Jose Manuel Corbí Garrido: /* Mediante el método de Runge-Kutta de cuarto orden */

{{ TrabajoED | Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral}}

{{ Beta }}
=Introducción y enunciado del problema=
En este artículo, se desarrollará el estudio de la capa límite de un fluido laminar sobre una placa plana semi-infinita mediante métodos numéricos programados en lenguaje M. De esta forma el lector podrá reproducir (si así lo desea) los algoritmos y códigos aquí escritos en MATLAB u otras alternativas de software libre como Octave.

[[Archivo:Grp9BIMAG1.png|miniaturadeimagen|Fluido laminar por encima de una placa plana]]

El tratamiento de este problema mediante métodos numéricos requiere una definición precisa de las condiciones iniciales y de contorno. En nuestro caso definiremos la semi-placa en el plano cartesiano como la semirecta <math>y=0</math>, con <math>x∈(0,∞)</math>; y estudiaremos lo que pasa por encima de ella <math>(x>0, y>0)</math>. El fluido laminar que se encuentra con la placa, antes de llegar a ella tiene una velocidad constante de <math>\overrightarrow{u} = u_0\cdot\overrightarrow{i}, u_0 = 2 </math>, iniciándose una fase de transición en su encuentro con la placa (la velocidad es nula en los puntos de la placa), y lejos del punto de encuentro, una vez finalizada la transición, volviendo a fluir a la velocidad constante inicial.
El campo de velocidades del fluido se define como::
<math>\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})</math>
donde <math>\psi</math> es la corriente del fluido::
<math>\psi(x,y) = \sqrt[]{ \nu \cdot\ u_0 \cdot x} f(\eta)</math> y <math>ν</math> es su viscosidad (que supondremos <math>ν = 1</math>), <math> f(\eta) </math> con <math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math> satisface la '''ecuación de Blasius'''::

<math>f''' + \frac{1}{2}f f'' =0</math>

En los siguientes apartados resolveremos numéricamente esta ecuación diferencial ordinaria no lineal, para unas condiciones iniciales::
<math>f(0)=f(0)'=0; f(0)''=k</math> <math> f'\rightarrow 1</math> cuando <math>\eta\rightarrow \infty</math>
tomando <math>k</math> valores <math> k \in \mbox{(0,1;1)}</math> y <math>dk=0,01</math>
A continuación estudiaremos la resolución de la ecuación mediante diferentes métodos numéricos, y estudiaremos el <math>k</math> "óptimo" para las condiciones iniciales enunciadas.

=Resolución numérica de la ecuación de Blasius=
Para poder resolver esta ecuación independientemente del posterior procedimiento numérico que vayamos a aplicar, primero hemos tenido que convertirla en un sistema, aplicando el método analítico que detallamos a continuación::

<math>y_{1}=f(x)</math>:

<math>y_{2}=f'(x)</math>:

<math>y_{3}=f''(x)</math>

Así <math>y_{3}'=f'''(x)</math>, si despejamos de nuestra ecuación diferencial obtenemos que::

<math>y_{2}=y_{1}'</math>:

<math>y_{3}=y_{2}'</math>:

<math>y_{3}'=-\frac{1}{2}y_{1}y_{3}</math>

Como resultado de esta transformación obtenemos un sistema para el cual podemos obtener su solución numéricamente.

==Mediante el método de Euler modificado==
A continuación se muestra el código en lenguaje matlab que proporciona las soluciones numéricas a la ecuación de Blasius para los diferentes valores de k (recordemos que tomaba valores de 0,1 hasta 1 con un paso de 0,01) en un intervalo <math> \eta \in \mbox{(0,20)}</math>. Hemos tomado una <math>h=0,5</math>.

Implementando el método de Euler en este sistema:
{{matlab|codigo=
%% Método de Euler modificado

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
%% Datos del problema y discretización
x0=0; % las x son las eta de las que depende f
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta, k y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada x
%% Iteraciones y(x_n)-->y(x_n+1)
for n=1:N
% calculamos k1=f(x_n,y_n)
k1=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
% Calculamos k2=f(x_n+h/2,y_n+1/2k1*h)
xp=x(n)+h; %argumento de x
fp=ff+h*k1; %argumento y
k2=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';

ff=ff+h/2*(k1+k2);
f(:,n+1)=ff;

end
laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solución
f2=laf2paraescogerk;
figure(1)
clf
plot(x,f)
legend('f(x)','f´(x)','f´´(x)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f, f´,f´´ en intervalo de eta, siendo k el que mejor aproxima f´(20) a 1')
corte=1;
figure (2)
hold on
plot(k,laf2paraescogerk)
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')
hold off
%Para calcular la k de forma exacta utilizamos este bucle que nos dice con qué k la diferencia en valor absoluto entre la f´(n) y uno es menor
for n=1:K
l=n;
p=abs(1-f2(n));
if min(abs(1-f2))==p

break
end
end
kfinal=k(l)
}}

[[Archivo:Grp9BIMAG2eulermodif1.jpg|miniaturadeimagen|f,f´,f´´ con k=0,33 para intervalo de eta de 0 a 20]]
[[Archivo:Grp9BIMAG2eulermodif2.jpg|border|f´(20) para los distintos valores de k]]

Como se puede apreciar en el gráfico, según este método la k que mejor se ajusta a las condiciones iniciales es <math>k=0.330</math>

==Mediante el método de Runge-Kutta de cuarto orden==
La siguiente secuencia de código da solución a la ecuación de Blasius mediante el método de Runge-Kutta, siendo aplicado al sistema de ecuaciones diferenciales obtenido en el primer apartado mediante procedimientos analíticos:

{{matlab|codigo=
% Resolver sistema con un método RK4
%% Método de Runge-Kutta orden 4

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
% Datos del problema y discretización
% La variable eta recibe el nombre x
x0=0;
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada tiempo

for n=1:N
% calculamos k1=f(x_n,f_n)
k1=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
% Calculamos k2=f(x_n+h/2,f_n+1/2k1*h)
xp=x(n)+h/2; %argumento de eta
fp=ff+1/2*h*k1; %argumento f
k2=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
%Calculamos k3=f(x_n+h/2,f_n+1/2k2*h)
%mismo argumento eta que antes
fp=ff+1/2*h*k2;%argumento y
k3=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
%calculamos k4=f(x_n+h/2,f_n+1/2k3*h)
xp=x(n)+h;%argumento eta
fp=ff+h*k3; %argumento f
k4=[fp(2) fp(3) -1/2*fp(1)*fp(3)]';
ff=ff+h/6*(k1+2*k2+2*k3+k4);
f(:,n+1)=ff;
end
laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solucion
f2=laf2paraescogerk;
figure(1)
clf
plot(x,f)
legend('f1(x)','f2(x)','f3(x)')
corte=1;
figure (2)
plot(k,laf2paraescogerk)
hold on
plot(k,corte,'red')
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')

for n=1:K
l=n;
p=abs(1-f2(n));
if min(abs(1-f2))==p

break
end
end

kfinal=k(l)
}}
[[Archivo:Grafipartado2d.jpeg|border]]

Según este método la k que mejor se ajusta a las condiciones iniciales es k=0.330, como podemos apreciar en el gráfico.

==Mediante el método de Euler==
Al igual que en los otros apartados resolvemos el sistema asociado a la ecuación de Blasius numéricamente, en este caso mediante el método de Euler:
{{matlab|codigo=
%% Método de Euler

% y1=f;
% y2=f';
% y3=f'';
% y3'=f'''=-1/2*f*f''=-1/2*y1*y3
% Datos del problema y discretización
% x es la variable eta del problema
x0=0;
xf=20;
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
k0=0.1;
kf=1;
k=0.1:0.01:1;
K=length(k);
for j=1:K
f0=[0 0 k(j)]';
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada valor de eta
%% Iteraciones f(x_n)-->f(x_n+1)
for n=1:N

y=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
ff=ff+h*y;
f(:,n+1)=ff;

end
laf2paraescogerk(j)=f(2,N);
end
%% dibujo la solucion
f2=laf2paraescogerk;
clf
corte=1;
figure (1)
plot(k,laf2paraescogerk)
hold on
plot(k,corte,'red')
legend('f´(20)')
xlabel('k')
ylabel('valores de f´(20)')
title('Representación de f(20)´ para los posibles k')
plot(k,corte,'red')

for n=1:K
l=n;
p=abs(1-f2(n));
if min(abs(1-f2))==p

break
end
end

kfinal=k(l)
figure(2)
%ya con la k=0,33

x0=0;
xf=20;
k=0.32
f0=[0 0 k]';
%% Datos discretizacion
h=0.05;
N=(xf-x0)/h;
%% Vectores de eta y aproximaciones
x=x0:h:xf;
f=zeros(3,N+1);
%% Inicialización
f(:,1)=f0;
ff=f0; % Valor de la aprox. en cada eta
%% Iteraciones f(x_n)-->f(x_n+1)
for n=1:N

y=[ff(2) ff(3) -1/2*ff(1)*ff(3)]';
ff=ff+h*y;
f(:,n+1)=ff;

end
plot(x,f)
legend('f(x)','f´(x)','f´´(x)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f, f´,f´´ en intervalo de eta, con Euler f´(20) a 1, ya con la k escogida')
}}

[[Archivo:Grafipartado2c.jpeg|border]]

Como se puede apreciar en el gráfico, según este método la k que mejor se ajusta a las condiciones iniciales es k=0.320

==Comparativa entre los tres métodos, elección de k==
La diferencia entre los tres métodos numéricos empleados para resolver la ecuacion de Blasius (Euler modificado, Runge-Kutta y Euler), es muy pequeña. La forma de las gráficas es completamente similar y el valor que se obtiene para que f'(20) sea lo más proxima a uno es en el caso de Euler modificado y Runge-Kutta igual a 0.33 y en Euler es 0.32. Por la tanto la diferencia entre los tres métodos es pequeña.

[[Archivo:Grp9BIMAG2eulermodif2.jpg|border|f´(20) para los distintos valores de k]]

[[Archivo:Grafipartado2d.jpeg|border|Con RK4]]

[[Archivo:Grafipartado2c.jpeg|border|Con Euler]]

'''De hecho a simple vista no se aprecia la diferencia.'''

=Representación e interpretación de <math>f’(\eta)</math>=
Vamos a estudiar a partir de qué valor de <math>\eta</math> la gráfica difiere 0,01 o menos de 1.

Para esto se puede escoger cualquiera de los métodos numéricos anteriormente utilizados. Y después, para que nos muestre la gráfica de <math>f’(\eta)</math>, se escriben las siguientes líneas de código:

{{matlab|codigo=
plot(eta,f(2,:),'g')
legend('f2(eta)')
xlabel('eta')
ylabel('valores de f, f´, f´´')
title('Representación de f´ para k=0,33')
}}:

[[Archivo:grafipartado3.jpg|grafipartado3.jpg]]

Observamos como tiende rápidamente al valor 1, por lo que confirma que la k escogida en los apartados anteriores es correcta.

Ahora para encontrar a partir de qué valor de <math>\eta</math> la gráfica difiere 0,01 o menos de 1; introducimos la siguiente secuencia de código:

{{matlab|codigo=
figure(2)
gg=abs(f(2,:)-1);
for n=1:N
l=n;
if gg(n)<0.01
break
end
end
nu0=l-1
loquequermos=eta0*h

plot(eta0,gg)

}}
[[Archivo:grafipartado3b.jpg|grafipartado3b.jpg]]

Visualmente vemos que el valor buscado está más o menos entre 5 y 6. Con la condición del código de que la diferencia en valor absoluto sea menor que 0.01 llegamos a que el valor que buscamos concretamente es 5,2.

=Campo de velocidades del fluido=
Definimos el campo de velocidades del fluido <math>\vec{u}</math> como::

<math>\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})</math>

==Cálculo y representación de <math>u_1</math>==
A continuación vamos a calcular analíticamente la componente de <math> u_1 </math> de la velocidad del fluido en función de <math>f( \eta)</math>

Para ello, desarrollamos la componente hasta la obtención de su expresión final::

<math>u_1=\frac{\partial \psi}{\partial y}=\frac{\partial }{\partial y}(\sqrt[]{\nu u_0 x} f(\eta))= \sqrt[]{\nu u_0 x} \frac{\partial f(\eta)}{\partial y}=\sqrt[]{\nu u_0 x} \frac{\partial f(\eta)}{\partial \eta} \frac{\partial \eta}{\partial y}</math>:

Definimos la función <math>\eta(x,y)</math> como::

<math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math>:

<math>u_1=\sqrt[]{\nu u_0 x} \sqrt[]{\frac{u_0}{\nu x}} f’(\eta)=u_0 f’(\eta)</math>

Resolviendo esta ecuación mediante el método de Rungge-Kutta para un <math>y \in (0,3)</math> y unos valores de <math>x \in \{0.05, 0.2, 0.4, 0.6, 0.8\}</math> mediante el siguiente código, que nos permite calcular las gráficas para cada una de las <math>x_{i}</math> particulares:

{{matlab|codigo=

% Componente u1 de la velocidad
% Introducimos las variables
clear all
u0=2;
nu=1;
xk=[0.05 0.2 0.4 0.6 0.8];
y0=0;
yf=3;
h=0.05;
y=y0:h:yf;
N=(yf-y0)/h;
% Iteraciones
for i=1:5
eta(i,:)=y*sqrt(u0/(nu*xk(i)));
% H es el paso para resolver el sistema
H=0.05*sqrt(u0/(nu*xk(i)));
k=0.33;
f0=[0 0 k]';
f=[f0(1);f0(2);f0(3)];
ff=f0;
for n=1:N
A=[0 1 0;0 0 1;(-ff(3)/2) 0 0];
j=ff+H*A*ff;
ff=ff+(H/2)*(A*ff+A*j);
f(:,n+1)=[ff(1);ff(2);ff(3)];
end
G(i,:)=f(2,:);
end
u1=u0*G;

%% Dibujamos la solución
hold on
plot(y,u1(1,:),'m')
plot(y,u1(2,:),'b')
plot(y,u1(3,:),'r')
plot(y,u1(4,:),'g')
plot(y,u1(5,:),'c')
legend('para x=0.05','para x=0.2','para x=0.4','para x=0.6','para x=0.8')
xlabel('y')
ylabel('u1')
title ('Componente u1 de la velocidad')
hold off

}}
[[Archivo:Gráfica del apartado 4.JPG|Gráfica del apartado 4.JPG]]

==Conclusión==
Como vemos en la gráfica, podemos deducir que el fluido cuando se mueve sobre el eje de abscisas, necesita una altura mayor para poder tener una velocidad limite mas estabilizada. Observamos que la primera línea es muy próxima a una recta pero que las siguiente según avanzamos en las x se van acercando a una figura más parabólica. Estas curvas representan la transición de la velocidad del fluido desde la capa (donde es nula) hasta que alcanza u0 aumenta a medida que avanzamos por la placa.

=Capa límite, transición entre la velocidad nula al comienzo de la placa y lejos de la misma=
Hemos comprobado que la altura (<math>y</math>) a la que el fluido alcanza la velocidad constante (<math>u_0=2</math>) depende de la distancia recorrida sobre la placa (<math>x</math>). Ahora vamos a dibujar la función <math>g(x)</math>, donde <math>g(x)</math> es el valor <math>y</math> para el cual <math>\eta=\eta_0</math>, para observar como es ese cambio de <math>y</math> en función de la <math>x</math>.
La relación entre estas variables viene dada por la expresión: <math>\eta = y \sqrt[]{ \frac{u_0}{\nu x}}</math>
Despejamos la <math>y</math> y suponemos <math>u_0=2</math>, <math>\nu=1</math> y <math>\eta=\eta_0=5.2</math>
<center><math>y=\eta_0 \sqrt[]{ \frac{x}{2}}=g(x)</math></center>
Representamos esta función en Matlab.
{{matlab|codigo=
% Dibujamos la función g(x), y en función eta0(cte.) y de x.
% Introducimos variables
eta0=5.2;
u0=2;
v=1;
x=0:0.05:10;
% Calculamos y
y=eta0*sqrt(v*x/u0);
% Dibujamos la gráfica
plot(x,y)
legend('g(x)')
xlabel('x')
ylabel('y')
title('Capa límite de transición de la velocidad')
}}
[[Archivo:Gráficaapartado5.JPG|Gráficaapartado5.JPG]]

Interpretamos la función <math>g(x)</math> como el borde de la capa límite que hace la transición entre la velocidad
nula en la placa y la velocidad lejos de la misma del fluido. Se confirma lo que ya podíamos ver en el apartado anterior. A medida que el fluido se aleja del origen de la placa (<math>x=0</math>), es decir, a medida que aumenta el valor de la variable <math>x</math>, el fluido alcanza la velocidad inicial constante a una mayor altura (<math>y</math>).
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T18:31:23Z

Jose Manuel Corbí Garrido: /* Representación e interpretación de f’(\eta) */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T18:30:42Z

Jose Manuel Corbí Garrido: /* Representación e interpretación de f’(\eta) */

Archivo:Grafipartado3b.jpg

2014-03-04T18:27:16Z

Jose Manuel Corbí Garrido:

Archivo:Grafipartado3.jpg

2014-03-04T18:26:47Z

Jose Manuel Corbí Garrido:

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T18:26:27Z

Jose Manuel Corbí Garrido: /* Representación e interpretación de f’(\eta) */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T18:07:30Z

Jose Manuel Corbí Garrido: /* Comparativa entre los tres métodos, elección de k */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T18:06:45Z

Jose Manuel Corbí Garrido: /* Comparativa entre los tres métodos, elección de k */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T18:06:20Z

Jose Manuel Corbí Garrido: /* Comparativa entre los tres métodos, elección de k */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T18:00:57Z

Jose Manuel Corbí Garrido: /* Mediante el método de Euler */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T18:00:28Z

Jose Manuel Corbí Garrido: /* Mediante el método de Runge-Kutta de cuarto orden */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:59:52Z

Jose Manuel Corbí Garrido: /* Mediante el método de Runge-Kutta de cuarto orden */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:58:23Z

Jose Manuel Corbí Garrido: /* Mediante el método de Runge-Kutta de cuarto orden */

Archivo:Grafipartado2d.jpeg

2014-03-04T17:56:48Z

Jose Manuel Corbí Garrido:

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:56:35Z

Jose Manuel Corbí Garrido: /* Mediante el método de Runge-Kutta de cuarto orden */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:54:15Z

Jose Manuel Corbí Garrido: /* Mediante el método de Runge-Kutta de cuarto orden */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:52:52Z

Jose Manuel Corbí Garrido: /* Mediante el método de Euler */

Archivo:Grafipartado2c.jpeg

2014-03-04T17:52:26Z

Jose Manuel Corbí Garrido:

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:51:42Z

Jose Manuel Corbí Garrido: /* Mediante el método de Euler */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:50:07Z

Jose Manuel Corbí Garrido: /* Comparativa entre los tres métodos, elección de k */

Archivo:Grafipartado2b.jpeg

2014-03-04T17:44:52Z

Jose Manuel Corbí Garrido:

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:44:33Z

Jose Manuel Corbí Garrido: /* Mediante el método de Euler */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:42:15Z

Jose Manuel Corbí Garrido: /* Mediante el método de Euler */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:41:08Z

Jose Manuel Corbí Garrido: /* Mediante el método de Euler */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:39:58Z

Jose Manuel Corbí Garrido: /* Mediante el método de Euler */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:30:36Z

Jose Manuel Corbí Garrido: /* Mediante el método de Runge-Kutta de cuarto orden */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:30:22Z

Jose Manuel Corbí Garrido: /* Mediante el método de Runge-Kutta de cuarto orden */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:29:58Z

Jose Manuel Corbí Garrido: /* Mediante el método de Runge-Kutta de cuarto orden */

Capa límite de un fluido laminar sobre una placa plana (Grupo 9-B)

2014-03-04T17:28:07Z

Jose Manuel Corbí Garrido: /* Mediante el método de Runge-Kutta de cuarto orden */