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

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 las fuerzas aplicadas transversalmente a la longitudinal de 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 [math]ρ[/math] nos indica la densidad y [math]τ[/math] la tensión. Con ello [math]\frac{τ}{ρ}=c^2[/math] nos dará la velocidad de propagación [math]c[/math] en [math]\frac{m}{s}[/math] que, en nuestro caso, será 1.

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

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)

2.2.1 Solución gráfica

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)

2.3.1 Solución gráfica

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

En este apartado vamos a realizar la aproximación del problema [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] mediante el método de Fourier con 1,3,5,10 y 20 autofunciones. Nuestra solución 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 series de Fourier.

Posición del cable mediante 5 series de Fourier.

Posición del cable mediante 10 series de Fourier.

Posición del cable mediante 20 series de Fourier.

En este apartado las figuras obtenidas mediante Fourier nos aproximan muy apropiadamente la solución. A medida que introducimos los diferentes números de autofunciones [math](1,3,5,10,20)[/math] observamos que la solución aproximada mejora considerablemente hacia la solución real.

.

4 Conservación de la Energía a lo largo del tiempo

En este apartado obtendremos la gráfica de la energía \begin{equation} E(t) = \int_{0}^{L} (u_{t}^{2}(x,t)+ u_{x}^{2}(x,t)) \ dx \end{equation} a lo largo del tiempo [math]t=40 s[/math]. Para obtenerla debemos aproximar la derivada [math]u_x[/math] mediante diferencias finitas con la aproximación [math]u_{x}\approx \frac{u(x_{n+1},t)-u(x_{n-1},t)}{2h}[/math]. Además [math]u_t[/math] será [math]V[/math] en [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]

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)
W(:,1)=W0;
for j=1:M
    W(:,j+1)=(eye(size(L))-(k/2)*L)\(W(:,j)+(k/2)*(L*W(:,j)+2*G));
end
%Condiciones Dirichlet.

U=W(1:n,:); 
UA=ua*ones(1,length(t)); UB=ub*ones(1,length(t)); U=[UA;U;UB];
 
V=W(n+1:end,:);
VA=ua*ones(1,length(t)); VB=ub*ones(1,length(t)); V=[VA;V;VB];
%====energia====
E=zeros(size(t)); 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); %Aproximación de Ux
    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);
end
plot(t,E)
%axis([0 40 0 0.6]) %Línea añadida para ver cómo se conserva la energía 
%a lo largo del tiempo 
xlabel('tiempo'); ylabel('energía');

Comportamiento de la energía a lo largo del tiempo con paso h

Comportamiento de la energía a lo largo del tiempo con paso h

.

Si dividimos el paso h por la mitad, es decir, [math]h=0.05[/math] obtenemos las siguientes gráficas:

Comportamiento de la energía a lo largo del tiempo con paso h/2

Comportamiento de la energía a lo largo del tiempo con paso h/2

.

En las gráficas vemos que se demuestra numéricamente que la energía se conserva a lo largo del tiempo. Las situadas a la izquierda tienen un zoom menor que las situadas a la derecha, así lo vemos de una forma más clara. El valor de la energía es de aproximadamente [math]0.45 J[/math]. Las dos últimas gráficas se han realizado con un paso de [math]h/2[/math] lo que permite ajustar mejor el resultado de la energía.

5 Energía en un cable sumergido en el mar

Ahora el cable estará sumergido en el mar, es decir, un medio que produce amortiguamiento. La ecuación ahora es [math]u_{tt}−u_{xx}+au_{t}=0[/math], donde a es una constante que depende del amortiguamiento que produce el medio. Dibujamos en una misma gráfica el comportamiento de la energía para [math]a=[0,1,4,10,100][/math]

[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)=\begin{cases} \frac{3x}{10} \\ \frac{3}{2}-\frac{3x}{20} \end{cases}, \; u_{t}(x,0)=0\\ \end{cases} [/math]

clear all 
clc

valores=[0,1,4,10,100];
for z=valores
    a=0; b=10;h=0.1; x=a:h:b;
    N=(b-a)/h;
    xx=x(1:N);
    xx=xx';
    tam=length(xx);
    ub=0;ua=0;
    U0=0*xx;
     for i=1:length(xx);
           if xx(i)<b/3
        U0(i)=3*xx(i)/b;
         else 
        U0(i)=1.5-1.5*xx(i)/b;
           end
     end
    V0=0*xx;

    %matriz K

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

    %Creación de F
    F=0.*xx;

    %tiempo
    t0=0;tf=40; k=h;
    t=t0:k:tf;
    M=length(t)-1;
    W=[U0;V0];

    for i=1:M   
        %Sistema
        G=[zeros(size(xx));F];
        L=[zeros(size(K)),eye(size(K));-K,-z*eye(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,:); 
    UB=ub*ones(1,length(t)); U=[U;UB];

 
    V=W(tam+1:end,:);
    VB=ub*ones(1,length(t)); V=[V;VB];
 
    %====energia====
 
    E=zeros(size(t)); Ux=zeros(size(x));
    for j=1:M+1
      for m=2:N
        Ux(m)=(U(m+1,j)-U(m-1,j))/(2*k); %Aproximación de Ux mediante diferencias finitas
      end
     Ux(1)=U(2,j)/k;
     Ux(end)=-U(end-1,j)/k;
     Ux=Ux';
        E(j)=trapz(x,V(:,j).^2)+trapz(x,Ux.^2); 
    end
  
hold on
plot(t,E)
xlabel('tiempo (s)'); ylabel('energía (J)');
clear all
end
legend ('a=0','a=1','a=4','a=10','a=100','Location','Best');
hold off

izquierda

En la gráfica adjunta vemos cómo evoluciona el valor de la energía a lo largo del tiempo. Observamos además que existe una curva para cada valor de [math]a[/math] asignado a cada una de las iteraciones. De arriba a abajo los valores de [math]a[/math] son [math][0,1,4,10,100][/math]. Vemos que, a medida que aumenta la constante que representa la fuerza de amortiguamiento, la energía se disipa a mayor velocidad en relación a los anteriores valores.

.

6 Comportamiento de la Energía. Vibraciones periódicas.

Ahora el extremo izquierdo del cable está sufriendo vibraciones periódicas sujeto a una estructura con frecuencia [math]F0[/math] Herzios. Tomamos [math]f(t)=sen(2πF0t)[/math]. [math]F0[/math] tendrá como valores [math]0.11,0.09[/math] y [math]0.1[/math]. Calculamos el comportamiento de la energía partiendo inicialmente del reposo en un tiempo [math]t∈[0,60][/math].

[math] \begin{cases} u_{tt}-u_{xx}=0, \; x∈[0,10], \; t∈[0,60], \\ u_x(0,t)=f(t), \; u(10,t)=0, \\ u(x,0)=0, \; u_{t}(x,0)=0\\ \end{cases} [/math]

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====

Ux=zeros(size(x)); E=zeros(size(t));
for l=1:M+1
    for m=2:N
        Ux(m)=(U(m+1,l)-U(m-1,l))/(2*k); %Aproximación de Ux mediante 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);
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)');

Posición del cable con F0=0.11

Comportamiento de la energía con F0=0.11

Posición del cable con F0=0.1

Comportamiento de la energía con F0=0.1

Posición del cable con F0=0.09

Comportamiento de la energía con F0=0.09

Es obvia la gran diferencia entre las gráficas. En primer lugar, con F0=0.11 y F0=0.09 la posición a la que llega el movimiento vertical del cable es mucho menor que con F0=0.1. La fuerza [math]f(t)[/math] provoca un movimiento en [math]x=0[/math] que con el tiempo se incrementa considerablemente desde el extremo izquierdo hacia el derecho aumentando con el tiempo debido a que las vibraciones dadas por [math]f(t)=sen(2πF0t)[/math] dependen de la variable [math]t[/math]. La energía en este caso no puede ser constante y así queda comprobado.

.

7 Comportamiento de la Energía. Nueva condición de contorno

El cable ahora está sujeto por el extremo izquierdo a un aparato que envía una respuesta a la vibración que recibe de manera que ahora la condición de contorno es [math]u_x (0,t)=bu(0,t)[/math]. Calculamos el comportamiento de la energía para [math]b=2,-2[/math]. Para ello utilizaremos la aproximación [math]u_{x}(0,t)−bu(0,T)\approx \frac{u_{1}(t)−u_{−1}(t)}{2h}−u_{0}(t)[/math].

Nuestro problema, por lo tanto, será [math] \begin{cases} u_{tt}-u_{xx}=0, \; x∈[0,10], \; t∈[0,40], \\ u_{x}(0,t)=bu(0,t), \; 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]

clear all 
clc
a=0; b=10;h=0.1; x=a:h:b;
N=(b-a)/h;
z=input('introducir valor elegido de b=');
xx=x(1:N);
xx=xx';
tam=length(xx);
ub=0;ua=0;
U0=0*xx;
for i=1:length(xx);
    if xx(i)<b/3
        U0(i)=3*xx(i)/b;
    else 
        U0(i)=1.5-1.5*xx(i)/b;
    end 
end
V0=0*xx;

%matriz K

K=(1/h^2)*(2*diag(ones(1,N))-diag(ones(1,N-1),-1)-diag(ones(1,N-1),1));
K(1,1)=(1/h^2)*(2*h*z+2);

%Creación de F
F=0.*xx;

%tiempo
t0=0;tf=40; k=h;
t=t0:k:tf;
M=length(t)-1;
W=[U0;V0];

for i=1:M   
    %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,:); 
 UB=ub*ones(1,length(t)); U=[U;UB];
 
V=W(tam+1:end,:);
VB=ub*ones(1,length(t)); V=[V;VB];


 
%====energia====
 
E=zeros(size(t));
Ux=zeros(size(x));
for j=1:M+1
    for m=2:N
        Ux(m)=(U(m+1,j)-U(m-1,j))/(2*k); %Aproximación de Ux
    end
    Ux(1)=U(2,j)/k;
    Ux(end)=-U(end-1,j)/k;
    Ux=Ux';
    E(j)=trapz(x,V(:,j).^2)+trapz(x,Ux.^2);
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)');

Evolución de la cuerda con b=2

Evolución de la energía de la cuerda b=2

Evolución de la cuerda con b=-2

Evolución de la energía de la cuerda b=-2

La nueva condición de contorno en [math]x=0[/math] tiene que ver con la fuerza descrita por [math]u_{x}(0,t)=bu(0,t)[/math], por lo tanto, dependiendo del valor de [math]b[/math] y de la posición de [math]u[/math] en [math]x=0[/math] a lo largo del tiempo el resultado variará. No es tan obvia en las gráficas tridimensionales de la posición pero sí en las relativas a la energía. Para [math]b=2[/math] en los primeros 3-4 segundos vemos un aumento de energía que sobrepasa los [math]0.447 Julios[/math] para posteriormente reducirse a [math]0.4445 Julios[/math] aproximadamente. Tras ello, vuelve a crecer, disminuir y aumentar hasta llegar a los 40 segundos. Para [math]b=-2[/math] la energía no crece tanto en los primeros instantes pero el cable sí mantiene una mayor cantidad de ella que para [math]b=2[/math] en el momento del descenso.