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

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Ecuación de ondas en el cable de una estructura civil (Grupo 9-B)
Asignatura Ecuaciones Diferenciales
Curso Curso 2013-14
Autores

María García Fernández

José Manuel Corbí Garrido

Estefanía Jiménez Ocampo

Diego García Vaquero

José Francisco Aguilera Corral

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 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 trata de ser inextensible, es decir, que intenta ser 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 \(u_{tt}-u_{xx}=0\). 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 \(u(x,t)\), dependiente de dos variables; la \(x\) que indica la distancia al extremo de la cuerda tomado como referencia, y la \(t\) que define el tiempo.

Las condiciones de frontera en este caso nos indican que el cable no sufrirá desplazamiento vertical en sus extremos \(u(0,t)=0; u(10,t)=0\). Por otra parte las condiciones iniciales \(u(x,0)=0; u_{t}(x,0)=0\) 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\[ \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} \]

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).

2 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 \(u(x,t)\). 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\[ \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} \]

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.

2.1 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 nuestro problema de ondas 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:

 1 %% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x?[0,10]; t>0 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.1;
15 N=L/h;
16 x=0:h:L;
17 xint=h:h:L-h;
18 
19 dt=h;
20 t=0:dt:T;
21 
22 k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
23 k=(1/h^2)*k;;
24 
25 F=zeros(N-1,1);
26 
27 u0=(2-abs(2-xint*2/5))';
28 v0=0*xint';
29 
30 M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
31 G=[zeros(N-1,1); F];
32 
33 %% Aproximación de u mediante el método del trapecio
34 
35 w0=[u0 ; v0];
36 ww=w0;
37 sol(1,:)=[ 0 , w0(1:N-1)' , 0];
38 
39 for j=1:length(t)-1
40     ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);
41     
42     sol(j+1,:)=[0 , ww(1:N-1)' , 0];
43 end
44 
45 %% Representación gráfica
46 
47 figure
48 [xx,tt]=meshgrid(x,t);
49 surf(xx,tt,sol)
50 xlabel('Cable')
51 ylabel('Tiempo')
52 zlabel('Desplazamiento vertical')
53 title('u aproximado mediante diferencias finitas con el método del trapecio')


Obteniendo la gráfica:

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 para que pueda ser o no factible de ejecutar. El siguiente código matlab nos permite ejecutarlo:

 1 %% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler explícito
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x?[0,10]; t>0 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.1;
15 N=L/h;
16 x=0:h:L;
17 xint=h:h:L-h;
18 
19 dt=h;
20 t=0:dt:T;
21 
22 k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
23 k=(1/h^2)*k;
24 
25 F=zeros(N-1,1);
26 
27 u0=(2-abs(2-xint*2/5))';
28 v0=0*xint';
29 
30 %% Aproximación de u mediante el método de Euler explícito
31 
32 sol(1,:)=[ 0 , u0' , 0];
33 u=u0;
34 v=v0;
35 for j=1:length(t)-1
36     u=u+dt*v+F;
37     v=v-dt*k*u+F;
38     sol(j+1,:)=[0 , u' , 0];
39 end
40 
41 %% Representación gráfica
42 
43 [xx,tt]=meshgrid(x,t);
44 surf(xx,tt,sol)
45 xlabel('Cable')
46 ylabel('Tiempo')
47 zlabel('Desplazamiento vertical')
48 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:

Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler explícito

El último método que empleamos para resolver el problema planteado mediante diferencias finitas es el método de Euler modificado, el cual tiene más precisión que el método de Euler explícito. El código que se muestra a continuación ejecuta dicho método y nos proporciona una gráfica del cable:

 1 %% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método de Euler modificado
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x?[0,10]; t>0 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.1;
15 N=L/h;
16 x=0:h:L;
17 xint=h:h:L-h;
18 
19 dt=h;
20 t=0:dt:T;
21 
22 k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
23 k=(1/h^2)*k;
24 
25 F=zeros(N-1,1);
26 
27 u0=(2-abs(2-xint*2/5))';
28 v0=0*xint';
29 
30 %% Aproximación de u mediante el método de Euler modificado
31 
32 sol(1,:)=[ 0 , u0' , 0];
33 u=u0;
34 v=v0;
35 for n=1:length(t)-1
36   l=k*u;
37     % calculamos k1=f(t_n,y_n)
38     k1=[ v l]';
39     % Calculamos k2=f(t_n+h/2,y_n+1/2k1*h)
40     %xp=x(n)+dt; %argumento de tiempo
41     %solp=sol+dt*k1; %argumento y 
42     k2=[v l]';
43     %Calculamos k3=f(t_n+h/2,y_n+1/2k2*h)
44     %mismo argumento tiempo que antes
45     
46     
47     u=u+dt/2*(k1(1,:)'+k2(1,:)')+F;
48     l=+k*u;
49     k1=[v l]';
50     k2=[v l]';
51     v=v-dt/2*(k1(2,:)'+k2(2,:)')+F;
52    % v=v-dt/2*(k*u+k*u)+F;
53    sol(n+1,:)=[0 , u' , 0];
54      
55 end  
56 
57 %% Representación gráfica
58 
59 [xx,tt]=meshgrid(x,t);
60 surf(xx,tt,sol)
61 xlabel('Cable')
62 ylabel('Tiempo')
63 zlabel('Desplazamiento vertical')
64 title('u aproximado mediante diferencias finitas con el método de Euler modificado')

La gráfica que se obtiene es:

Ecuación de ondas para el cable de una estructura civil. Aproximación mediante diferencias finitas con el método de Euler modificado

2.2 Aproximación mediante el método de Fourier

2.2.1 Con un término de la serie

 1 %% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con un término
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x?[0,10]; t>0 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.1;
15 N=L/h;
16 dt=0.1; 
17 M=T/dt;
18 x=0:h:10;
19 posicion=(2-abs(2-x*2/5));
20 velocidad=0*x;
21 t=0:dt:T;
22 
23 %% Aproximación de u mediante Fourier con un término (Q=1)
24 
25 Q=1;
26 for k=1:Q
27 p=sin(k*pi/10*x) ;
28 fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
29 
30 end
31 [XX,TT]=meshgrid(x,t);
32 U=0;
33 
34 for k=1:Q
35 U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
36 end
37 
38 %% Representación gráfica
39 
40 surf(XX,TT,U)
41 xlabel('Cable')
42 ylabel('Tiempo')
43 zlabel('Desplazamiento vertical')
44 title('u aproximado mediante Fourier con un término')

Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con un término

2.2.2 Con tres términos de la serie

 1 %% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con tres términos
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x?[0,10]; t>0 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.1;
15 N=L/h;
16 dt=0.1; 
17 M=T/dt;
18 x=0:h:10;
19 posicion=(2-abs(2-x*2/5));
20 velocidad=0*x;
21 t=0:dt:T;
22 
23 %% Aproximación de u mediante Fourier con tres términos (Q=3)
24 
25 Q=3;
26 
27 for k=1:Q
28     p=sin(k*pi/10*x) ;
29     fourier1(k)=trapz(x,posicion.*p)/trapz(x,p.*p);
30  end
31 
32 [XX,TT]=meshgrid(x,t);
33 U=0;
34 
35 for k=1:Q
36     U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
37 end
38 
39 %% Representación gráfica
40 
41 surf(XX,TT,U)
42 xlabel('Cable')
43 ylabel('Tiempo')
44 zlabel('Desplazamiento vertical')
45 title('u aproximado mediante Fourier con tres términos')

Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con tres términos

2.2.3 Con cinco términos de la serie

 1 %% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con cinco términos
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x?[0,10]; t>0 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.1;
15 N=L/h;
16 dt=0.1; 
17 M=T/dt;
18 x=0:h:10;
19 posicion=(2-abs(2-x*2/5));
20 velocidad=0*x;
21 t=0:dt:T;
22 
23 %% Aproximación de u mediante Fourier con cinco términos (Q=5)
24 
25 Q=5;
26 
27 for k=1:Q
28     p=sin(k*pi/10*x) ;
29     fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
30 end
31 
32 [XX,TT]=meshgrid(x,t);
33 U=0;
34 
35 for k=1:Q
36     U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
37 end
38 
39 %% Representación gráfica
40 
41 surf(XX,TT,U)
42 xlabel('Cable')
43 ylabel('Tiempo')
44 zlabel('Desplazamiento vertical')
45 title('u aproximado mediante Fourier con cinco términos')

Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con cinco términos

2.2.4 Con diez términos de la serie

 1 %% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con diez términos
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x?[0,10]; t>0 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.1;
15 N=L/h;
16 dt=0.1; 
17 M=T/dt;
18 x=0:h:10;
19 posicion=(2-abs(2-x*2/5));
20 velocidad=0*x;
21 t=0:dt:T;
22 
23 %% Aproximación de u mediante Fourier con diez términos (Q=10)
24 
25 Q=10;
26 
27 for k=1:Q
28     p=sin(k*pi/10*x) ;
29     fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
30 end
31 
32 [XX,TT]=meshgrid(x,t);
33 U=0;
34 
35 for k=1:Q
36     U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
37 end
38 
39 %% Representación gráfica
40 
41 surf(XX,TT,U)
42 xlabel('Cable')
43 ylabel('Tiempo')
44 zlabel('Desplazamiento vertical')
45 title('u aproximado mediante Fourier con diez términos')

Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con diez términos

2.2.5 Con veinte términos de la serie

 1 %% Aproximación del desplazamiento vertical del cable (u) mediante Fourier con veinte términos
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x?[0,10]; t>0 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.1;
15 N=L/h;
16 dt=0.1; 
17 M=T/dt;
18 x=0:h:10;
19 posicion=(2-abs(2-x*2/5));
20 velocidad=0*x;
21 t=0:dt:T;
22 
23 %% Aproximación de u mediante Fourier con diez términos (Q=20)
24 
25 Q=20;
26 
27 for k=1:Q
28     p=sin(k*pi/10*x) ;
29     fourier1(k)=trapz(x,posicion.*p)/trapz(x, p.*p);
30     
31 end
32 
33 [XX,TT]=meshgrid(x,t);
34 U=0;
35 
36 for k=1:Q
37     U=U+(fourier1(k)*cos(k*pi/10*TT)).*sin(k*pi/10*XX);
38 end
39 
40 %% Representación gráfica
41 
42 surf(XX,TT,U)
43 xlabel('Cable')
44 ylabel('Tiempo')
45 zlabel('Desplazamiento vertical')
46 title('u aproximado mediante Fourier con veinte términos')

Ecuación de ondas para el cable de una estructura civil. Aproximación mediante Fourier con veinte términos

2.3 Comparativa de diferencias finitas y Fourier. Interpretación de las gráficas anteriores

Todas las gráficas obtenidas mediante diferencias finitas (con sus distintos métodos) y Fourier representan la aproximación numérica de la solución del problema planteado, es decir, el desplazamiento vertical de todos los puntos de nuestro cable a lo largo del tiempo.

Considerando los órdenes de precisión de cada uno de los métodos empleados en los programas de diferencias finitas, para aproximar la solución de nuestro problema (la función \(u(x,t)\) que representa el desplazamiento vertical): Podemos concluir que el método del trapecio (orden 2) es el más preciso, seguido del método de Euler modificado y por último del método de Euler explícito (orden 1).

En cuanto a la resolución del problema mediante Fourier, 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. El aumento del número de términos considerados en Fourier hace el método más estable.

3 Energía del cable

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

 1 %% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x?[0,10]; t>0 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.05;
15 N=L/h;
16 x=0:h:L;
17 xint=h:h:L-h;
18 
19 dt=h;
20 t=0:dt:T;
21 
22 k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
23 k=(1/h^2)*k;
24 
25 F=zeros(N-1,1);
26 
27 u0=(2-abs(2-xint*2/5))';
28 v0=0*xint';
29 
30 M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
31 G=[zeros(N-1,1); F];
32 
33 %% Aproximación de u mediante el método del trapecio
34 
35 w0=[u0 ; v0];
36 ww=w0;
37 sol(1,:)=[ 0 , w0(1:N-1)' , 0];
38 
39 for j=1:length(t)-1
40     ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);
41     sol(j+1,:)=[0 , ww(1:N-1)' , 0];
42 end
43 
44 [xx,tt]=meshgrid(x,t);
45 surf(xx,tt,sol)
46 
47 %% Aproximación de u_t y u_x basándonos en el concepto de derivada
48 
49 v=size(sol);
50 for i=2:v(1)-1
51   for j=2:v(2)-1
52     Ux(i,j)=(sol(i+1,j)-sol(i-1,j))/(2*h);
53     Ut(i,j)=(sol(i,j+1)-sol(i,j-1))/(2*h);
54   end
55 end
56 
57 %% Cálculo de la energía
58 
59 e=Ux.^2+Ut.^2;
60 xaux=[0,xint];
61 
62 for k=1:v(1)-2
63     E(k)=trapz(xaux,e(k,:));
64 end
65 t(401)=[ ];
66 
67 %% Representación gráfica de la energía a lo largo del tiempo
68 
69 plot(t(2:end-1),E(2:end))
70 xlabel('Tiempo')
71 ylabel('Energía')
72 title('Energía a lo largo del tiempo')


Energía del cable a lo largo del tiempo)


Se observa que el valor de la energía total del cable, para un paso del tiempo de 0.05 segundos, oscila aproximadamente entre 1.560 y 1.590. Este rango se reduce si elegimos un paso más pequeño para el tiempo y el espacio. Por tanto, aunque la aproximación no es perfecta, sí podemos confirmar que la energía mecánica del sistema se conserva.

4 Ejemplos de aplicación de cables en estructuras civiles

4.1 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\[ \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} \]

Vamos a estudiar lo que ocurre con la energía en estas condiciones, dependiendo del parámetro \(a\). Damos a \(a\) los valores de 0, 1 ,4, 10 y 100. Para ello utilizamos el método numérico de las diferencias finitas empleando Euler con un paso muy pequeño. En el caso de que a valga cero estamos ante el mismo caso que si no estuviese sumergido.

 1 %% Cable sumergido en un medio viscoso (p.ej el mar)con a=0
 2 % Cambiando el valor de a definido más adelante obtenemos todas las gráficas
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}+a*u_{t}=0; x[0,10]; t[0,40]; 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=0; u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.1;
15 N=L/h;
16 x=0:h:L;
17 xint=h:h:L-h;
18 
19 dt=h/2;
20 t=0:dt:T;
21 
22 k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
23 
24 k=(1/h^2)*k;
25 
26 F=zeros(N-1,1);
27 
28 u0=(2-abs(2-xint*2/5))';
29 v0=0*xint';
30 
31 
32 %% Aproximación de u mediante el método de Euler y obtención de las derivadas
33 
34 sol(1,:)=[ 0 , u0' , 0];
35 u=u0;
36 v=v0;
37 a=0;
38 for j=1:length(t)-1
39     u=u+dt*v;
40     v=v-dt*k*u-dt*a*v;
41     sol(j+1,:)=[0 , u' , 0];
42 end
43 
44 %% Gráficas
45 
46 [xx,tt]=meshgrid(x,t);
47 
48 figure
49 
50 surf(xx,tt,sol);
51 xlabel('Cable')
52 ylabel('Tiempo')
53 zlabel('Desplazamiento vertical')
54 title('u aproximado mediante diferencias finitas con el método de Euler explícito')
55 
56 figure
57 
58 ffff=sol(1,:);
59 plot(x,ffff);
60 xlabel('Cable')
61 ylabel('Desplazamiento inicial del cable')
62 title('Desplazamiento inicial de los puntos del cable')
63 
64 v=size(sol);
65 for i=1:v(1)-1
66   for j=1:v(2)-1
67     Ux(i,j)=(sol(i+1,j)-sol(i,j))/h;
68     Ut(i,j)=(sol(i,j+1)-sol(i,j))/h;
69   end
70 end
71 
72 e=Ux.^2+Ut.^2;
73 xaux=[0,xint];
74 laux=length(xaux);
75 otro=length(e(2,:));
76 for k=1:v(1)-1
77   E(k)=trapz(xaux,e(k,:));
78  end
79 t(401)=[ ];
80 lE=length(E);
81 lt=length(t);
82 
83 figure
84 
85 plot(t,E)
86 xlabel('Tiempo')
87 ylabel('Energía')
88 title ('Energía con a=0')


Cable sumergido en un medio viscoso (p.ej el mar)con a=0 Cable sumergido en un medio viscoso (p.ej el mar)con a=1 Cable sumergido en un medio viscoso (p.ej el mar)con a=4 Cable sumergido en un medio viscoso (p.ej el mar)con a=10 Cable sumergido en un medio viscoso (p.ej el mar)con a=100

Observamos que, al sumergir el cable en un medio viscoso, el desplazamiento vertical de las ondas es cada vez más pequeño. Esto se debe a que la energía disminuye con el tiempo. Además, al aumentar el valor de \(a\), la energía se disipa más rápidamente. Como el medio es más viscoso, la oposición al desplazamiento vertical del cable es mayor. Para el caso de a=100, el método no es estable y da problemas. Sería necesario reducir aún más el paso temporal y espacial o utilizar otro método de mayor precisión.

4.2 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 \(F_{0}\) Herzios. Tomamos por ejemplo \(f(t) = sin(2 \pi F_{0})\). 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 \(F_{0}\) para observar si se produce alguna diferencia.

Tomamos \(F_{0}\) 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\[ \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} \]

Lo resolvemos numéricamente por diferencias finitas utilizando Euler.

  1 %% Cable sujeto a una estructura sometida a vibraciones periódicas
  2 
  3 %% Enunciado del problema
  4 
  5 % u_{tt}-u_{xx}=0; x[0,10]; t[0,60] 
  6 % u(0,t)=0; u(10,t)=sin(2*pi*f0)
  7 % u(x,0)=0; u(5,0)=2\\  u_{t}(x,0)=0
  8 
  9 %% Planteamiento y discretización (f0=0.11)
 10  
 11 clear all
 12 L=10;
 13 T=60;
 14 h=0.1;
 15 N=L/h;
 16 x=0:h:L;
 17 xint=h:h:L;
 18  
 19 dt=h;
 20 t=0:dt:T;
 21  
 22 k=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
 23 k(N,N-1)=-2;
 24  
 25 k=(1/h^2)*k;
 26  
 27 F=zeros(N-1,1);
 28  
 29 u0=0*xint';
 30 v0=0*xint';
 31  
 32 %% Aproximación de u mediante el método de Euler explícito y obtención de las derivadas (f0=0.11)
 33  
 34 f0=0.11;
 35 sol(1,:)=[ 0 , u0'];
 36 u=u0;
 37 v=v0;
 38 for j=1:length(t)-1
 39     u=u+dt*v+[F;sin(2*pi*f0*j)];
 40     v=v-dt*k*u+[F;sin(2*pi*f0*j)];
 41     sol(j+1,:)=[0 , u' ];
 42 end
 43  
 44 [xx,tt]=meshgrid(x,t);
 45 figure(1)
 46 surf(xx,tt,sol)
 47  
 48 v=size(sol);
 49 for i=1:v(1)-1
 50   for j=1:v(2)-1
 51     Ux(i,j)=(sol(i+1,j)-sol(i,j))/h;
 52     Ut(i,j)=(sol(i,j+1)-sol(i,j))/h;
 53   end
 54 end
 55 
 56 %% Energía (f0=0.11)
 57 
 58 e=Ux.^2+Ut.^2;
 59 xaux=[xint];
 60 laux=length(xaux);
 61 otro=length(e(2,:));
 62 
 63 for k=1:v(1)-1
 64     E(k)=trapz(xaux,e(k,:));
 65 end
 66 
 67 t(401)=[ ];
 68 
 69 %% Representación Energía-Tiempo (f0=0.11)
 70 
 71 figure(2)
 72 plot(t,E)
 73 xlabel('Tiempo')
 74 ylabel('Energía')
 75 title('Gráfica energía-tiempo f0=0.11')
 76 
 77 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 78  
 79 %% Planteamiento y discretización (f0=0.09)
 80  
 81 clear all
 82 L=10;
 83 T=60;
 84 h=0.1;
 85 N=L/h;
 86 x=0:h:L;
 87 xint=h:h:L;
 88  
 89 dt=h;
 90 t=0:dt:T;
 91  
 92 k=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
 93 k(N,N-1)=-2;
 94  
 95 k=(1/h^2)*k;
 96  
 97 F=zeros(N-1,1);
 98  
 99 u0=0*xint';
100 v0=0*xint';
101  
102 %% Aproximación de u mediante el método de Euler explícito y obtención de las derivadas (f0=0.09)
103  
104 f0=0.09;
105 sol(1,:)=[ 0 , u0'];
106 u=u0;
107 v=v0;
108 
109 for j=1:length(t)-1
110     u=u+dt*v+[F;sin(2*pi*f0*j)];
111     v=v-dt*k*u+[F;sin(2*pi*f0*j)];
112     sol(j+1,:)=[0 , u' ];
113 end
114  
115 [xx,tt]=meshgrid(x,t);
116 figure(3)
117 surf(xx,tt,sol)
118 
119 
120 v=size(sol);
121 for i=1:v(1)-1
122   for j=1:v(2)-1
123     Ux(i,j)=(sol(i+1,j)-sol(i,j))/h;
124     Ut(i,j)=(sol(i,j+1)-sol(i,j))/h;
125   end
126 end
127 
128 %% Energía (f0=0.09)
129 
130 e=Ux.^2+Ut.^2;
131 xaux=[xint];
132 laux=length(xaux);
133 otro=length(e(2,:));
134 for k=1:v(1)-1
135     E(k)=trapz(xaux,e(k,:));
136 end
137 
138 t(401)=[ ];
139 
140 %% Representación Energía-Tiempo (f0=0.09)
141 
142 figure(4)
143 plot(t,E)
144 xlabel('Tiempo')
145 ylabel('Energía')
146 title('Gráfica energía-tiempo f0=0.09')
147 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
148 
149 %% Planteamiento y discretización (f0=0.1)
150 
151 clear all
152 L=10;
153 T=60;
154 h=0.1;
155 N=L/h;
156 x=0:h:L;
157 xint=h:h:L;
158  
159 dt=h;
160 t=0:dt:T;
161  
162 k=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
163 k(N,N-1)=-2;
164  
165 k=(1/h^2)*k;
166  
167 F=zeros(N-1,1);
168 
169 u0=0*xint';
170 v0=0*xint';
171  
172 %% Aproximación de u mediante el método de Euler explícito y obtención de las derivadas (f0=0.1)
173  
174 f0=0.1;
175 sol(1,:)=[ 0 , u0'];
176 u=u0;
177 v=v0;
178 for j=1:length(t)-1
179     u=u+dt*v+[F;sin(2*pi*f0*j)];
180     v=v-dt*k*u+[F;sin(2*pi*f0*j)];
181     sol(j+1,:)=[0 , u' ];
182 end
183  
184 [xx,tt]=meshgrid(x,t);
185 figure(5)
186 surf(xx,tt,sol)
187 v=size(sol);
188 
189 for i=1:v(1)-1
190   for j=1:v(2)-1
191     Ux(i,j)=(sol(i+1,j)-sol(i,j))/h;
192     Ut(i,j)=(sol(i,j+1)-sol(i,j))/h;
193   end
194 end
195 
196 %% Energía (f0=0.1)
197 
198 e=Ux.^2+Ut.^2;
199 xaux=[xint];
200 laux=length(xaux);
201 otro=length(e(2,:));
202 
203 for k=1:v(1)-1
204     E(k)=trapz(xaux,e(k,:));
205 end
206 
207 t(401)=[ ];
208 
209 %% Representación Energía-Tiempo (f0=0.1)
210 
211 figure(6)
212 plot(t,E)
213 xlabel('Tiempo')
214 ylabel('Energía')
215 title('Gráfica energía-tiempo f0=0.1')

Cable sujeto a una estructura sometida a vibraciones periódicas f0=0.11 Cable sujeto a una estructura sometida a vibraciones periódicas f0=0.09 Cable sujeto a una estructura sometida a vibraciones periódicas f0=0.1

Se observa que cuanto mayor es el valor de \(F_{0}\), la energía alcanza un máximo más elevado. Esto es debido a que las vibraciones en el extremo derecho de la estructura dependen de una función senoidal que está en función de \(F_{0}\). A medida que aumenta \(F_{0}\) crece el seno y, con ello, la derivada parcial en x, esto se traduce en un incremento de la energía. El valor de \(F_{0}\) únicamente afecta a los máximos de la energía pero no a su carácter temporal. Se trata de una función periódica que depende del tiempo. Para cualquier \(F_{0}\), el periodo es de 40 segundos.

4.3 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 \(u_{x}(L,t)=bu(L,t)\). Tomamos como condición inicial \(u(x,0)=2-|2-\frac{2}{5}x|\) y observando en qué cambia la energía en función de distintos valores de b. estudiamos dos casos concretos: b=-2 y b=2.

\( \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} \)

Resolvemos este problema por diferencias finitas utilizando el método de Euler, con un paso pequeño para que la aproximación sea más exacta.


 1 %% Cable sujeto a un aparato que envía una respuesta a las vibraciones (b=2)
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x[0,10]; t[0,40]; 
 6 % u(0,t)=0; u(10,t)=bu(10,0) 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 b=2;
13 L=10;
14 T=40;
15 h=0.01;
16 N=L/h;
17 x=0:h:L;
18 xint=h:h:L;
19 
20 dt=h/2;
21 t=0:dt:T;
22 
23 k=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
24 k(N,N-1)=-2;
25 k(N,N)=2+b*(2*h);
26 
27 k=(1/h^2)*k;
28 
29 F=zeros(N,1);
30 
31 u0=(2-abs(2-xint*2/5))';
32 v0=0*xint';
33 
34 %% Aproximación de u mediante el método de Euler y obtención de las derivadas
35 
36 f0=0.1;
37 sol(1,:)=[ 0 , u0'];
38 u=u0;
39 v=v0;
40 for j=1:length(t)-1
41     u=u+dt*v;
42     v=v-dt*k*u;
43     sol(j+1,:)=[0 , u' ];
44 end
45 
46 [xx,tt]=meshgrid(x,t);
47 figure(1)
48 surf(xx,tt,sol)
49   
50 v=size(sol);
51 for i=2:v(1)-1
52   for j=2:v(2)-1
53     Ux(i,j)=(sol(i+1,j)-sol(i-1,j))/(2*h);
54     Ut(i,j)=(sol(i,j+1)-sol(i,j-1))/(2*h);
55   end
56 end
57 e=Ux.^2+Ut.^2;
58 xaux=[0,xint];
59 
60 %% Energía
61 
62 for k=1:v(1)-2
63     E(k)=trapz(xaux(2:end),e(k,:));
64 end
65 t(401)=[ ];
66 plot(t(2:end-1),E(2:end))
67 xlabel('Tiempo')
68 ylabel('Energía')
69 title('Gráfica energía-tiempo para b=2')

Cable sujeto a un aparato que envía una respuesta a las vibraciones (b=2)

No podemos apreciar para valores de b positivos (como el estudiado, b=2) cambios importantes en la energía. Esta no aumenta ni disminuye de manera significativa.

 1 %% Cable sujeto a un aparato que envía una respuesta a las vibraciones (b=-2)
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x[0,10]; t[0,40]; 
 6 % u(0,t)=0; u(10,t)=bu(10,0) 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 b=-2;
13 L=10;
14 T=40;
15 h=0.01;
16 N=L/h;
17 x=0:h:L;
18 xint=h:h:L;
19 
20 dt=h/2;
21 t=0:dt:T;
22 
23 k=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
24 k(N,N-1)=-2;
25 k(N,N)=2+b*(2*h);
26 
27 k=(1/h^2)*k;
28 
29 F=zeros(N,1);
30 
31 u0=(2-abs(2-xint*2/5))';
32 v0=0*xint';
33 
34 %% Aproximación de u mediante el método de Euler y obtención de las derivadas
35 
36 f0=0.1;
37 sol(1,:)=[ 0 , u0'];
38 u=u0;
39 v=v0;
40 for j=1:length(t)-1
41     u=u+dt*v;
42     v=v-dt*k*u;
43     sol(j+1,:)=[0 , u' ];
44 end
45 
46 [xx,tt]=meshgrid(x,t);
47 figure(1)
48 surf(xx,tt,sol)
49   
50 v=size(sol);
51 for i=2:v(1)-1
52   for j=2:v(2)-1
53     Ux(i,j)=(sol(i+1,j)-sol(i-1,j))/(2*h);
54     Ut(i,j)=(sol(i,j+1)-sol(i,j-1))/(2*h);
55   end
56 end
57 e=Ux.^2+Ut.^2;
58 xaux=[0,xint];
59 
60 %% Energía
61 
62 for k=1:v(1)-2
63     E(k)=trapz(xaux(2:end),e(k,:));
64 end
65 t(401)=[ ];
66 plot(t(2:end-1),E(2:end))
67 xlabel('Tiempo')
68 ylabel('Energía')
69 title('Gráfica energía-tiempo para b=-2')

Cable sujeto a un aparato que envía una respuesta a las vibraciones (b=-2)

Se observa que cuanto mayor es la tensión del cable, mayor es el desplazamiento vertical. Se está otorgando energía al sistema. De está forma, para valores de b negativos (por ejemplo b=-2), la energía crece de manera exponencial y tiende a infinito. Se produce un efecto de resonancia, incrementándose la energía.

5 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.

 1 %% Aproximación del desplazamiento vertical del cable (u) mediante diferencias finitas con el método del trapecio
 2 
 3 %% Enunciado del problema
 4 
 5 % u_{tt}-u_{xx}=0; x?[0,10]; t>0 
 6 % u(0,t)=0; u(10,t)=0 
 7 % u(x,0)=2-abs(2-\frac{2}{5}x); u(5,0)=2\\  u_{t}(x,0)=0\\
 8 
 9 %% Planteamiento y discretización
10 
11 clear all
12 L=10;
13 T=40;
14 h=0.1;
15 N=L/h;
16 x=0:h:L;
17 xint=h:h:L-h;
18 
19 dt=h;
20 t=0:dt:T;
21 
22 k=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
23 k=(1/h^2)*k;;
24 
25 F=zeros(N-1,1);
26 
27 u0=(2-abs(2-xint*2/5))';
28 v0=0*xint';
29 
30 M=[zeros(N-1) eye(N-1); -k zeros(N-1) ];
31 G=[zeros(N-1,1); F];
32 
33 %% Aproximación de u mediante el método del trapecio
34 
35 w0=[u0 ; v0];
36 ww=w0;
37 sol(1,:)=[ 0 , w0(1:N-1)' , 0];
38 
39 for j=1:length(t)-1
40     ww=(eye(2*N-2)-dt/2*M)\((eye(2*N-2)+dt/2*M)*ww+dt*G);
41     
42     sol(j+1,:)=[0 , ww(1:N-1)' , 0];
43 end
44 
45 %% Representación gráfica
46 
47 figure
48 [xx,tt]=meshgrid(x,t);
49 surf(xx,tt,sol)
50 xlabel('Cable')
51 ylabel('Tiempo')
52 zlabel('Desplazamiento vertical')
53 title('u aproximado mediante diferencias finitas con el método del trapecio')
54 
55 %% Animación.
56 
57 % Muestra el desplazamiento vertical de cada punto del cable en cada instante
58 
59 figure
60 [n,m]=size(sol);
61 aux=zeros(size(sol));
62 for j=1:n
63     % Se va analizando el desplazamiento vertical en cada instante (j) de
64     % todos los puntos (i) del cable
65     for i=1:m
66         aux(j,i)=sol(j,i);
67     end
68     surf(xx,tt,aux); % Es la representación de todos los puntos (i) del cable en un instante dado
69     xlabel('Cable')
70     ylabel('Tiempo')
71     zlabel('Desplazamiento vertical')
72     zlim([-2 2])
73     MFRAMES=getframe;   % Guarda una foto del cable en cada instante
74 end
75 movie=MFRAMES;
76 
77 figure
78 for j=1:n
79     % Se va analizando el desplazamiento vertical en cada instante (j) de
80     % todos los puntos (i) del cable
81     for i=1:m
82         aux(j,i)=sol(j,i);
83     end
84     plot(x,aux(j,:)); % Es la representación de todos los puntos (i) del cable en un instante dado
85     xlabel('Cable')
86     ylabel('Desplazamiento vertical')
87     ylim([-2 2])
88     MFRAMES=getframe;   % Guarda una foto del cable en cada instante
89 end
90 movie=MFRAMES;