Ecuación de Ondas aplicada a un Cable (Grupo 2A)

De MateWiki
Revisión del 15:12 14 may 2015 de JM Alonso de Caso (Discusión | contribuciones) (→‎Aproximación mediante el Método de Fourier)

Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Ecuación de Ondas aplicada a un Cable. Grupo 2-A
Asignatura Ecuaciones Diferenciales
Curso Curso 2014-15
Autores Juan Raúl Ruíz Méndez (531); Jaime Enrech Martínez (532); Manuel Mudarra Hernández (551); Jose Manuel Alonso de Caso Gilsanz (618); Guillermo Díaz Rivera (649); Iago Rodriguez Romero (824)
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Ecuación de ondas y Cable a considerar

Se nos presenta un cable utilizado para estructura civiles de longitud L=10m, con una sección que se nos permite despreciar con respecto a su longitud. Además, las vibraciones que este sufre se pueden modelizar mediante la ecuación de ondas. Expresaremos como el desplazamiento vertical la función [math]u(x,t)[/math] que sufren los distintos puntos del cable en función de su posición a lo largo del mismo (discretizados en [math]x ∈ [0,L] [/math]) y en distintos instantes de tiempo t.

Inicialmente, el cable se encuentra en un medio de viscosidad nula. Por ello, en la ecuación ya mencionada no aparece ningún término que amortigüe el desplazamiento de los puntos del cable. El cable se encuentra sujeto por ambos extremos, lo que nos facilita las primeras condiciones de frontera (condiciones Dirichlet). Concretamente, establecemos que tanto el desplazamiento del extremo izquierdo (x=0) como el extremo derecho (x=10) son nulos.

Se considera que el cable es perfectamente homogéneo en toda su longitud. Por tanto, no existen imperfecciones que alteren su comportamiento. Consecuentemente, la función a la que igualamos nuestra ecuación de ondas f(x) , la cual indica la fuerza aplicada transversal a la cuerda, es nula. Si las vibraciones que sufre el cable son lo suficientemente pequeñas, la función [math]u(x,t)[/math] satisface la ecuación de ondas [math]u_{xx} - u_{xx}=f(x)[/math]. El sistema completo que muestra la situación del cable es el siguiente:

[math] \begin{cases} ρu_{tt} - τu_{xx}=0), \; x∈[0,10], \; t∈[0,T],\\ u(0,t)=0, \; u(10,t)=0,\\ u(x,0)=h_0(x), \; u_{t}(x,0)=h_1(x)\\ \end{cases} [/math]

Donde ρ nos indica la densidad y τ la tensión. Con ello [math]\frac{τ}{ρ}=c^2[/math] nos dará la velocidad de propagación en [math]\frac{m}{s}[/math] que, en nuestro caso, será 1.

2 Aproximación mediante el Método de líneas. Trapecio. Euler. Heun

El cable estará sujeto en [math]x=10/3[/math] y desplazado 1 m en dirección perpendicular. Empezará a vibrar una vez que lo soltemos. La posición del cable la aproximaremos mediante el método de líneas con un paso de longitud de [math]∆x = 0.1[/math] y un paso de tiempo [math]∆t = ∆x[/math]. Para ello utilizaremos los métodos del trapecio, euler explícito y heun para el tiempo [math]t ∈ [0, 40][/math]. Nuestro problema vendrá dado por: [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)=\begin{cases} \frac{3x}{10} \\ \frac{3}{2}-\frac{3x}{20} \end{cases}, \; u_{t}(x,0)=0\\ \end{cases} [/math]

Para el método de líneas realizamos la aproximación de [math]u_{xx}\approx \frac{u(x_{n-1},t)-2u(x_n,t)+u(x_{n+1},t)}{h^2}[/math]. Mediante ello llegaremos al problema dado por [math] \begin{cases} U''=-KU+F\\ U(0)=u^{0}\\ U'(0)=v^{0}\\ \end{cases} [/math] el cual reduciremos a un sistema de orden 1 con el cambio de variable [math] \begin{cases} U'=V\\ V'=-KU+F\\ \end{cases} [/math]. Tendremos


[math] \begin{cases} \begin{bmatrix} U \\ V \end{bmatrix}' = \begin{bmatrix} 0 & 1 \\ -K & 0 \end{bmatrix} \begin{bmatrix} U \\ V \end{bmatrix} + \begin{bmatrix} 0 \\ F \end{bmatrix}\\ \begin{bmatrix} U(0) \\ V(0) \end{bmatrix} = \begin{bmatrix} U^{0} \\ V^{0} \end{bmatrix} \end{cases} [/math]que sustituiremos por [math] \begin{cases} W'=LW+G\\ W(0)=W^{0}\\ \end{cases} [/math].

2.1 Método del Trapecio

Este método viene definido por [math]W^{m+1}=W^{m}+\frac{k}{2}({L^{m}W^{m}+G^{m}+L^{m+1}W^{m+1}+G^{m+1}})[/math] que despejando [math]W^{m+1}[/math] tendremos [math]W^{m+1}=\frac{W^{m}+\frac{k}{2}(L^{m}W^{m}+G^{m}+G^{m+1})}{I-\frac{k}{2}(L^{m+1})}[/math] 

clear all
a=0; b=10;
h=0.1;
x=a:h:b;
N=round(b-a)/h;
%Condiciones de Contorno
ua=0; ub=0;
xx=x(2:N); %Extraemos los elementos interiores del vector x
xx=xx';%Trasponemos el vector para poder utilizarlo en el sistema
U0=zeros(size(xx)); %Generamos un vector de ceros al que luego introduciremos los valores de la función según el punto del cable en el que nos encontremos
n=length(xx);
for i=1:n
    if xx(i)<=10/3
        U0(i)=(3/10)*xx(i); %Condición Inicial Posición: Función a ejecutar si x<=10/3
    else
        U0(i)=1.5-(3/20)*xx(i);  %Condición Inicial Posición: Función a ejecutar si x>10/3
    end
end
V0=0*xx; %Condición Inicial Velocidad
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1)); %Matriz de coeficientes del sistema
F=0*xx;
%Añadimos el primer y último elemento de las condiciones iniciales
F(1)=F(1)+ua/h^2; %Extremo izquierdo (x=0)
F(end)=F(end)+ub/h^2; %Extremo derecho (x=10)
%Procedemos a la discretización del tiempo
t0=0; tM=40;
k=h;
t=t0:k:tM;
M=(tM-t0)/k;
%El sistema a resolver es W'=LW+G, W(0)=W(super)0
W0=[U0;V0]; %Valor Inicial Sistema de orden 1
G=[zeros(size(xx));F]; %Vector término independiente Sistema de orden 1
L=[zeros(size(K)),eye(size(K));-K,zeros(size(K))];
%RESOLUCIÓN DEL SISTEMA (Por TRAPECIO)
sol(:,1)=W0;
for j=1:M
    sol(:,j+1)=(eye(size(L))-(k/2)*L)\(sol(:,j)+(k/2)*(L*sol(:,j)+2*G));
end
sol=sol(1:N-1,:);
%Tomamos las N-1 primeras filas e incluimos las condiciones de contorno
UA=ua*ones(1,length(t));
UB=ub*ones(1,length(t));
sol=[UA;sol;UB];
%Hay que tomar la parte que afecta a los primeros L/3 del cable de sol1 y
%la parte que afecta al resto del cable de sol2
[Mt,Mx]=meshgrid(t,x);
set(gcf,'renderer','painters');
mesh(Mx,Mt,sol)

2.1.1 Solución gráfica

Posición del Cable mediante Método del Trapecio

En la figura con ejes [math]Longitud-Tiempo-Posición[/math] se puede apreciar el desplazamiento del cable a lo largo de [math]L=10 m[/math] durante 40 segundos. Vemos que inicialmente el cable se sitúa 1 metro por encima de su posición de equilibrio en [math]x=10/3 m[/math]. Tras ello, su distancia vertical oscilará entre [math][1,-1][/math] a lo largo del tiempo. Esto es debido a la no presencia de fuerzas externas. Este método es muy adecuado para estos problemas.














.

2.2 Método de Euler explícito

La expresión a tratar de este método es [math]W^{m+1}=W^{m}+k(L^{m}W^{m}+G^{m})[/math] 

clear all
a=0; b=10;
h=0.1;
x=a:h:b;
N=round(b-a)/h;
%Condiciones de Contorno
ua=0; ub=0;
xx=x(2:N); %Extraemos los elementos interiores del vector x
xx=xx';%Trasponemos el vector para poder utilizarlo en el sistema
U0=zeros(size(xx)); %Generamos un vector de ceros al que luego introduciremos los valores de la función según el punto del cable en el que nos encontremos
n=length(xx);
for i=1:n
    if xx(i)<=10/3
        U0(i)=(3/10)*xx(i); %Condición Inicial Posición: Función a ejecutar si x<=10/3
    else
        U0(i)=1.5-(3/20)*xx(i);  %Condición Inicial Posición: Función a ejecutar si x>10/3
    end
end
V0=0*xx; %Condición Inicial Velocidad
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1)); %Matriz de coeficientes del sistema
F=0*xx;
%Añadimos el primer y último elemento de las condiciones iniciales
F(1)=F(1)+ua/h^2; %Extremo izquierdo (x=0)
F(end)=F(end)+ub/h^2; %Extremo derecho (x=10)
%Procedemos a la discretización del tiempo
t0=0; tM=40;
k=h;
t=t0:k:tM;
M=(tM-t0)/k;
%El sistema a resolver es W'=LW+G, W(0)=W(super)0
W0=[U0;V0]; %Valor Inicial Sistema de orden 1
G=[zeros(size(xx));F]; %Vector término independiente Sistema de orden 1
L=[zeros(size(K)),eye(size(K));-K,zeros(size(K))];
%RESOLUCIÓN DEL SISTEMA (Por Euler Explícito)
sol(:,1)=W0;
for j=1:M
    sol(:,j+1)=sol(:,j)+k*(L*sol(:,j)+G);
end
sol=sol(1:N-1,:);
%Tomamos las N-1 primeras filas e incluimos las condiciones de contorno
UA=ua*ones(1,length(t));
UB=ub*ones(1,length(t));
sol=[UA;sol;UB];
%Hay que tomar la parte que afecta a los primeros L/3 del cable de sol1 y
%la parte que afecta al resto del cable de sol2
[Mt,Mx]=meshgrid(t,x);
set(gcf,'renderer','painters');
mesh(Mx,Mt,sol)


Posición del cable por el método de Euler explícito

En la figura con ejes [math]Longitud−Tiempo−Posición[/math] se puede apreciar el desplazamiento del cable a lo largo de [math]L=10m[/math] durante 40 segundos. Vemos que este método no es especialmente adecuado para aproximar el problema obteniendo una posición de [math]10^{136}[/math] lo cual no es correcto. Con un paso más pequeño podría obtenerse una solución más cercana a la real.













.

2.3 Método de Heun

[math]W^{m+1}=W^{m}+\frac{k}{2}(K1+K2)[/math] es la expresión a tratar donde [math] \begin{cases} K1=L^{m}W^{m}+G^{m}\\ K2=W^{m}+kK1\\ \end{cases} [/math] 

clear all
a=0; b=10;
h=0.1;
x=a:h:b;
N=round(b-a)/h;
%Condiciones de Contorno
ua=0; ub=0;
xx=x(2:N); %Extraemos los elementos interiores del vector x
xx=xx';%Trasponemos el vector para poder utilizarlo en el sistema
U0=zeros(size(xx)); %Generamos un vector de ceros al que luego introduciremos los valores de la función según el punto del cable en el que nos encontremos
n=length(xx);
for i=1:n
    if xx(i)<=10/3
        U0(i)=(3/10)*xx(i); %Condición Inicial Posición: Función a ejecutar si x<=10/3
    else
        U0(i)=1.5-(3/20)*xx(i);  %Condición Inicial Posición: Función a ejecutar si x>10/3
    end
end
V0=0*xx; %Condición Inicial Velocidad
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1)); %Matriz de coeficientes del sistema
F=0*xx;
%Añadimos el primer y último elemento de las condiciones iniciales
F(1)=F(1)+ua/h^2; %Extremo izquierdo (x=0)
F(end)=F(end)+ub/h^2; %Extremo derecho (x=10)
%Procedemos a la discretización del tiempo
t0=0; tM=40;
k=h;
t=t0:k:tM;
M=(tM-t0)/k;
%El sistema a resolver es W'=LW+G, W(0)=W(super)0
W0=[U0;V0]; %Valor Inicial Sistema de orden 1
G=[zeros(size(xx));F]; %Vector término independiente Sistema de orden 1
L=[zeros(size(K)),eye(size(K));-K,zeros(size(K))];
%RESOLUCIÓN DEL SISTEMA (Por HEUN)
sol(:,1)=W0;
for j=1:M
    K1=L*sol(:,j)+G;
    K2=sol(:,j)+K1*k;
    sol(:,j+1)=sol(:,j)+k/2*(K1+K2);
end
sol=sol(1:N-1,:);
%Tomamos las N-1 primeras filas e incluimos las condiciones de contorno
UA=ua*ones(1,length(t));
UB=ub*ones(1,length(t));
sol=[UA;sol;UB];
%Hay que tomar la parte que afecta a los primeros L/3 del cable de sol1 y
%la parte que afecta al resto del cable de sol2
[Mt,Mx]=meshgrid(t,x);
set(gcf,'renderer','painters');
mesh(Mx,Mt,sol)
Posición del cable mediante el método de Heun

En la figura con ejes [math]Longitud−Tiempo−Posición[/math] se puede apreciar el desplazamiento del cable a lo largo de [math]L=10m[/math] durante 40 segundos. Al igual que el método de Euler, con este método no se obtiene una aproximación correcta del problema ya que la figura marca una posición de [math]10^{69}[/math] lo cual es imposible que suceda. Con un paso más pequeño podría obtenerse una solución más cercana a la real.













.

3 Aproximación mediante el Método de Fourier

[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)=\begin{cases} \frac{3x}{10} \\ \frac{3}{2}-\frac{3x}{20} \end{cases}, \; u_{t}(x,0)=0\\ \end{cases} [/math] Nuestra solución mediante este método tendrá la forma [math]u(x,t)=\sum_{i=1}^k X_k(x)T_k(t)[/math]. Sustituyendo en la ecuación obtenemos [math]T(t)''X(x)-T(t)X(x)''=-λ[/math] y [math]\frac{T(t)''}{T(t)}=\frac{X(x)''}{X(x)}=-λ[/math] . Nuestro [math] PA \begin{cases} X(x)''+λX(x)=0\\ X(0)=0\\ X(10)=0\\ \end{cases} [/math] tendrá como autovalores [math]λ_k=\frac{k^{2}π^2}{10^{2}}[/math] para [math]k=1,2,3,...[/math] y [math]X_k=sen(\frac{kπx}{10})[/math] como autofunciones. Queremos resolver ahora \begin{cases} T(t)+λT(t)=0\\ T(0)=\sum_{i=1}^k c_k sen(\frac{kπx}{10})\\ T'(0)=0\\ \end{cases}

[math]c_k= \frac{\int_a^b f(x)sen(\frac{kπx}{10})\,dx}{\int_a^b sen^2(\frac{kπx}{10})\,dx} [/math] donde [math]f(x)=\begin{cases} \frac{3x}{10} \\ \frac{3}{2}-\frac{3x}{20} \end{cases}[/math]. Nuestra solución será [math]u_k(x,t)=\sum_{i=1}^k c_k cos(\frac{kπt}{10})sen(\frac{kπx}{10})[/math]. 

clear all
%datos del problema
a=0;b=10; %espacio, extremos de la cuerda
h=0.1;%En x--------Paso espacial=h
Q=input('Introduzca número de autofunciones: ');
%discretización en x
x=a:h:b;
t=0:0.1:40;
[Mx,Mt]=meshgrid(x,t);
u0=zeros(size(x)); %primera función valor inicial
  for i=1:length(x)
      if x(i)<10/3 
       u0(i)=(3*x(i))/10;
      else 
       u0(i)=1.5-1.5*x(i)/10;
      end
  end
u0t=0; %segunda función valor inicial
U=0;
 for j=1:Q
    p=sin((j*pi/10)*x);
    cj=trapz(x,u0.*p)/trapz(x,p.*p);
    U=U+(cj.*cos((j*pi*Mt)/10)).*sin((j*pi*Mx)/10);
 end
mesh(Mx,Mt,U)


Posición del cable mediante 1 serie de Fourier.
Posición del cable mediante 3 serie de Fourier.
Posición del cable mediante 5 serie de Fourier.
Posición del cable mediante 10 serie de Fourier.
Posición del cable mediante 20 serie de Fourier.

4 Energía

clear all 
clc
a=0; b=10; h=0.1; x=a:h:b; N=(b-a)/h;
xx=x(2:N); xx=xx'; tam=length(xx);
ub=0; %contorno x=L
U0=0*xx; V0=0*xx; %cond. iniciales
%Matriz K
K=(1/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
K(1,2)=0; 

%Creación de F 
F=0.*xx;
F0=input('introduzca la frecuencia F0 en Hz aplicada al cable ');

%tiempo
t0=0; tf=60; k=0.1; t=t0:k:tf; M=length(t)-1;

W=[U0; V0];
 
for i=1:M
    F(1)=-2.*sin(2*pi*F0*t(i))/h^2;
    
    %Sistema 
    G=[zeros(size(xx));F] ; 
    L=[zeros(size(K)),eye(size(K));-K,zeros(size(K))];

    %EULER 
    %W(:,i+1)=W(:,i)+k*(L*W(:,i)+G);  
    %TRAPECIO 
    W(:,i+1)=(eye(size(L))-(k/2)*L)\(W(:,i)+(k/2)*(L*W(:,i)+2*G));
    
end

U=W(1:tam,:); 
UA=sin(2*pi*F0*t); UB=ub*ones(1,length(t)); U=[UA;U;UB];

V=W(tam+1:end,:);
VA=2*pi*F0*cos(2*pi*F0*t); VB=ub*ones(1,length(t)); V=[VA;V;VB];

%====energia====

E=zeros(size(t)); %Preasignación.
Ux=zeros(size(x));
for l=1:M+1
    for m=2:N
        Ux(m)=(U(m+1,l)-U(m-1,l))/(2*k); %Cálculo de la derivada Ux mediante la aproximación por diferencias finitas.
    end
    Ux(1)=U(2,l)/k;
    Ux(end)=-U(end-1,l)/k;
    Ux=Ux';
    E(l)=trapz(x,V(:,l).^2)+trapz(x,Ux.^2); %Cálculo de la energía.
end

figure(1)
[Mt,Mx]=meshgrid(t,x);
set(gcf,'renderer','painters');
mesh(Mx,Mt,U);
xlabel('Longitud (m)'); ylabel('tiempo (s)'); zlabel('altura de los puntos (m)');

figure(2)
plot(t,E)
xlabel('tiempo (s)'); ylabel('energía (J/kg)');
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_de_Ondas_aplicada_a_un_Cable_(Grupo_2A)&oldid=29994»