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	Nuestros nombres
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 L = 10 m 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 u(x,t), podemos plantear el problema de su movimiento según el siguiente sistema de ecuaciones:

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

clear all, clf
%Coeficiente de amortiguamiento
am=[0,1,4,10,100];
for n=am
    %Datos en x
    a=0; b=10;
    h=0.1;
    x=a:h:b;
    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));
    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));
    %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;
    M=length(t)-1; %número de subintervalos
    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));
    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);
        end
        Ux=Ux';
        E(l)=trapz(x,V(:,l).^2)+trapz(x,Ux.^2);
    end
    %Dibujamos la gráfica de la energía
    hold on
    plot(t,E)
    %Borramos todos los datos para realizar el bucle de nuevo
    clear all
end
hold off