Ecuación de ondas. Grupo C2

Trabajo realizado por estudiantes
Título	Ecuación de ondas. Grupo C2
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Ana Martínez Lorente, Javier Parras Martínez, Alfredo Pazos Arjona, Antonio Pérez Mata, Javier Siguero Ginés
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Consideramos un cable de una estructura civil de longitud [math] L = 10 m [/math] sujeto por ambos extremos. Supondremos que el cable tiene una sección pequeña respecto a su longitud y que las vibraciones pueden modelizarse mediante la ecuación de ondas. Si denotamos su desplazamiento vertical por [math] u(x,t) [/math], podemos plantear el problema de su movimiento según el siguiente sistema de ecuaciones: [math] \begin{cases} u_{tt}-u_{xx}=f(x,t), \; x∈[0,10], \; t∈[0,T],\\ u(0,t)=g_{0}(t), \; u(L,t)=g_{1}(t),\\ u(x,0)=h_{0}(x), \; u_{t}(x,0)=h_{1}(x)\\ \end{cases} [/math]

2 Vibración sin amortiguamiento. Condiciones Dirichlet. Resolución por el método de líneas

Vamos a tratar el problema de vibración de un cable de longitud [math] L=10 [/math], en la cual los dos extremos de la misma se encuentran fijos a lo largo del tiempo y con una desplazamiento nulo. Al inicio, en [math] t=0 [/math], sujetamos el cable desde el punto [math] x=L/3 [/math], y lo desplazamos 1 m en la dirección perpendicular a la recta que une sus extremos.

El problema viene modelizado por la siguente ecuación: [math] \begin{cases} u_{tt}-u_{xx}=0, x∈[0,L], t∈[0,40], \\ u(0,t)=0, u(L,t)=0, \\ u(x,0)=\begin{cases} \frac{3x}{L} \\ \frac{3}{2}-\frac{3x}{2L} \end{cases}, u_{t}(x,0)=0\\ \end{cases} [/math] A continuación se resuelve el problema por el método de diferencias finitas, aplicando para la resolución de la ecuación matricial que aparece los métodos del trapecio, de Euler explícito y de Heun.

2.1 Método del trapecio

Gráfica de la posición de cada punto del cable.

clear all
%Datos en x
a=0; b=10; %Longitud del cable L=10.
h=0.1; %Paso.
x=a:h:b; %Discretización espacial del cable.
N=round((b-a)/h);
%Definimos las matrices de la ecuación
xx=x(2:N);
xx=xx';
ua=0;ub=0; %Condiciones de contorno.
U0=zeros(size(xx)); %Preasignación de U0.
%Recorremos mediante un bucle U0, y añadimos los valores que correspondan.
for j=1:length(xx);
    if xx(j)<b/3
        U0(j)=3*xx(j)/b;
    else
        U0(j)=1.5-1.5*xx(j)/b;
    end
end
V0=zeros(size(xx)); %Preasignación de V0.
%Matriz K
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
%Término F y valor inicial
F=0*xx;
F(1)=F(1)+ua/h^2;
F(end)=F(end)+ub/h^2;
%Resolución del sistema de ecuaciones de EDO de orden 1.
t0=0;tM=40;
k=h; %Paso en t.
t=t0:k:tM; %Discretización del vector de tiempos.
M=length(t)-1; %Número de subintervalos.
%Añadimos en la primera columna las condiciones iniciales.
U(:,1)=U0;
V(:,1)=V0;
for i=1:M
    %Resolución del sistema de ecuaciones por el método del trapecio
    U(:,i+1)=(eye(size(K))+0.25*(k^2)*K)\(U(:,i)+0.5*k*(2*V(:,i)+0.5*k*(-K*U(:,i)+2*F)));
    V(:,i+1)=(eye(size(K))+0.25*(k^2)*K)\(V(:,i)+0.5*k*(-K*U(:,i)+2*F)-0.5*k*K*(U(:,i)+0.5*k*V(:,i)));
end
%Incluimos condiciones Dirichlet.
UA=ua*ones(1,length(t));
UB=ub*ones(1,length(t));
U=[UA;U;UB];
%Dibujamos el gráfico.
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,U)
xlabel('Longitud (m)'); ylabel('Tiempo (s)'); zlabel('Posición (m)');

2.2 Método de Euler explícito

Gráfica de la posición de cada punto del cable.

clear all
%Datos en x
a=0; b=10; %Longitud del cable L=10.
h=0.1; %Paso.
x=a:h:b; %Discretización espacial del cable.
N=round((b-a)/h);
%Definimos las matrices de la ecuación
xx=x(2:N);
xx=xx';
ua=0;ub=0; %Condiciones de contorno.
U0=zeros(size(xx)); %Preasignación de U0.
%Recorremos mediante un bucle U0, y añadimos los valores que correspondan.
for j=1:length(xx);
    if xx(j)<b/3
        U0(j)=3*xx(j)/b;
    else
        U0(j)=1.5-1.5*xx(j)/b;
    end
end
V0=zeros(size(xx)); %Preasignación de V0.
%Matriz K
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
%Término F y valor inicial
F=0*xx;
F(1)=F(1)+ua/h^2;
F(end)=F(end)+ub/h^2;
%Resolución del sistema de ecuaciones de EDO de orden 1.
t0=0;tM=40;
k=h; %Paso en t.
t=t0:k:tM; %Discretización del vector de tiempos.
M=length(t)-1; %Número de subintervalos.
%Añadimos en la primera columna las condiciones iniciales.
U(:,1)=U0;
V(:,1)=V0;
for i=1:M
    %Resolución del sistema de ecuaciones por el método de Euler (explícito).
    U(:,i+1)=U(:,i)+k*V(:,i);
    V(:,i+1)=V(:,i)+k*(-K*U(:,i)+F);
end
%Incluimos condiciones Dirichlet.
UA=ua*ones(1,length(t));
UB=ub*ones(1,length(t));
U=[UA;U;UB];
%Dibujamos el gráfico.
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,U)
xlabel('Longitud (m)'); ylabel('Tiempo (s)'); zlabel('Posición (m)');

2.3 Método de Heun

Gráfica de la posición de cada punto del cable.

clear all
%Datos en x
a=0; b=10; %Longitud del cable L=10.
h=0.1; %Paso.
x=a:h:b; %Discretización espacial del cable.
N=round((b-a)/h);
%Definimos las matrices de la ecuación
xx=x(2:N);
xx=xx';
ua=0;ub=0; %Condiciones de contorno.
U0=zeros(size(xx)); %Preasignación de U0.
%Recorremos mediante un bucle U0, y añadimos los valores que correspondan.
for j=1:length(xx);
    if xx(j)<b/3
        U0(j)=3*xx(j)/b;
    else
        U0(j)=1.5-1.5*xx(j)/b;
    end
end
V0=zeros(size(xx)); %Preasignación de V0.
%Matriz K
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
%Término F y valor inicial
F=0*xx;
F(1)=F(1)+ua/h^2;
F(end)=F(end)+ub/h^2;
%Resolución del sistema de ecuaciones de EDO de orden 1.
t0=0;tM=40;
k=h; %Paso en t.
t=t0:k:tM; %Discretización del vector de tiempos.
M=length(t)-1; %Número de subintervalos.
%Añadimos en la primera columna las condiciones iniciales.
U(:,1)=U0;
V(:,1)=V0;
for i=1:M
    %Resolución del sistema de ecuaciones por el método de Heun
    K1U=V(:,i);
    K2U=V(:,i)+K1U*k;
    K1V=-K*U(:,i)+F;
    K2V=-K*(U(:,i)+K1V*k)+F;
    U(:,i+1)=U(:,i)+(k/2)*(K1U+K2U);
    V(:,i+1)=V(:,i)+(k/2)*(K1V+K2V);
end
%Incluimos condiciones Dirichlet.
UA=ua*ones(1,length(t));
UB=ub*ones(1,length(t));
U=[UA;U;UB];
%Dibujamos el gráfico.
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,U)
xlabel('Longitud (m)'); ylabel('Tiempo (s)'); zlabel('Posición (m)');

3 Vibración con amortiguamiento. Cálculo de la energía

Gráfica de la energía para distintos amortiguamientos

clear all, clf
%Coeficiente de amortiguamiento
am=[0,1,4,10,100];
%Hacemos un bucle donde calcular la energía para cada coeficiente.
for n=am
    %Datos en x
    a=0; b=10; %Longitud del cable L=10.
    h=0.1; %Paso.
    x=a:h:b; %Discretización espacial del cable.
    N=round((b-a)/h);
    %Definimos las matrices de la ecuación
    xx=x(2:N);
    xx=xx';
    ua=0;ub=0; %Condiciones de contorno.
    U0=zeros(size(xx)); %Preasignación de U0.
    %Recorremos mediante un bucle U0, y añadimos los valores que correspondan.
    for j=1:length(xx);
        if xx(j)<b/3
            U0(j)=3*xx(j)/b;
        else
            U0(j)=1.5-1.5*xx(j)/b;
        end
    end
    V0=zeros(size(xx)); %Preasignación de V0.
    %Matriz K
    K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
    %Término F y valor inicial
    F=0*xx;
    F(1)=F(1)+ua/h^2;
    F(end)=F(end)+ub/h^2;
    %Resolución del sistema de ecuaciones de EDO de orden 1
    t0=0;tM=40;
    k=h; %Paso en t.
    t=t0:k:tM; %Discretización del vector de tiempos.
    M=length(t)-1; %Número de subintervalos.
    %Añadimos en la primera columna las condiciones iniciales.
    U(:,1)=U0;
    V(:,1)=V0;
    for i=1:M
        %Sistema de ecuaciones por el método del trapecio
        V(:,i+1)=((1+0.5*k*n)*eye(size(K))+0.25*(k^2)*K)\(V(:,i)+0.5*k*(-n*V(:,i)-K*U(:,i)+2*F)-0.5*k*K*(U(:,i)+0.5*k*V(:,i)));
        U(:,i+1)=U(:,i)+0.5*h*(V(:,i)+V(:,i+1));
    end
    %Incluimos condiciones Dirichlet.
    UA=ua*ones(1,length(t));
    UB=ub*ones(1,length(t));
    U=[UA;U;UB];
    %Como las condiciones Dirichlet son nulas, las velocidades de estos puntos
    %también lo serán
    V=[UA;V;UB];
    %Energía
    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=Ux';
        E(l)=trapz(x,V(:,l).^2)+trapz(x,Ux.^2); %Cálculo de la energía.
    end
    %Dibujamos la gráfica de la energía.
    hold on
    plot(t,E)
    xlabel('Tiempo (s)'); ylabel('Energía (J)');
    %Borramos todos los datos para realizar el bucle de nuevo.
    clear all
end
legend('a=0','a=1','a=4','a=10','a=100','Location','Best');
hold off

4 Cambio de condiciones en los extremos. Cálculo de la energía

Consideramos que nuestro cable está sujeto a una estructura que sufre vibraciones periódicas con frecuencia F0 Herzios. Vamos a tomar la función [math]f(t)=\sin(2\pi*F0*t)[/math] que será la que defina la posición del extremo izquierdo, que está sujeto a la estructura, en función del tiempo. Por tanto, nuestro problema queda: [math] \begin{cases} u_{tt}-u_{xx}=0, x∈[0,10], t∈[0,60],\\ u(0,t)=\sin(2\pi*F0*t), u(10,t)=0,\\ u(x,0)=0, u_{t}(x,0)=0.\\ \end{cases} [/math]

clear all
%Datos en x
a=0; b=10; %Longitud del cable L=10.
h=0.1; %Paso.
x=a:h:b; %Discretización espacial del cable.
N=round((b-a)/h);
%Definimos las matrices de la ecuación
xx=x(2:N);
xx=xx';
ub=0; %condición de contorno en el extremo derecho.
%Preasignación de la posición y la velocidad incial.
U0=zeros(size(xx));
V0=U0;
%Matriz K
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
%Término F y valor inicial
F=0*xx;
%Resolución del sistema de ecuaciones de EDO de orden 1
t0=0;tM=60; %Discretización del vector de tiempos.
k=h; %Paso en t.
t=t0:k:tM;
M=length(t)-1; %Número de subintervalos.
%Añadimos en la primera columna las condiciones iniciales.
U(:,1)=U0;
V(:,1)=V0;
%Pedimos por teclado al usuario los distintos valores de la frecuencia que
%le transmite la estructura al cable. Estas serán F0=1/L+0.01 Hz, 
%F0=1/L-0.01 Hz y F0=1/L Hz, siendo L=b=10.
F0=input('Introduzca la frecuencia (Hz) transmitida al cable: ');
for i=1:M
    F(1)=F(1)+sin(2*pi*F0*t(i))/h^2;
    F(end)=F(end)+ub/h^2;
    %Resolución del sistema de ecuaciones por el método del trapecio.
    U(:,i+1)=(eye(size(K))+0.25*(k^2)*K)\(U(:,i)+0.5*k*(2*V(:,i)+0.5*k*(-K*U(:,i)+2*F)));
    V(:,i+1)=(eye(size(K))+0.25*(k^2)*K)\(V(:,i)+0.5*k*(-K*U(:,i)+2*F)-0.5*k*K*(U(:,i)+0.5*k*V(:,i)));
end
%Incluimos las condiciones Dirichlet en nuestra solución. 
UA=sin(2*pi*F0*t);
UB=ub*ones(1,length(t));
U=[UA;U;UB];
VB=zeros(1,length(t));
VA=VB;
V=[VA;V;VB];
%Dibujamos el gráfico.
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,U)
xlabel('Longitud (m)'); ylabel('Tiempo (s)'); zlabel('Posición (m)');
%Cadena de texto con la frecuencia introducida, para el título de la gráfica.
Frec=sprintf('Frecuencia = %.2f Hz',F0); 
title(Frec);
%Energía
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=Ux';
    E(l)=trapz(x,V(:,l).^2)+trapz(x,Ux.^2); %Cálculo de la energía.
end
%Dibujamos la gráfica de la energía.
figure(2)
plot(t,E)
xlabel('Tiempo (s)'); ylabel('Energía (J)');

Gráfica de la posición del cable en cada punto y en cada instante.

Gráfica de la posición del cable en cada punto y en cada instante.

Gráfica de la posición del cable en cada punto y en cada instante.