Diferencia entre revisiones de «Ecuación de ondas. Grupo C2»

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