%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Master Pro II SNM
%2006 / 2007
%Simulation de l'èqwuation de la corde vibrante
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
        %
        % 1) données physiques
        %
            c=1;%vitesse
            mu=1;%%amortissement
        %données initiales (corde d'une guitare)
            u0=inline('2*(0.5-abs(x-0.5))')
        %
        %2) Données numériques
        %
            dt=0.001;%pas de temps
            Npdt=4000;%nb pas de temps
            Nper=10;%périodicité des calculs
            N=255;%nb de points intérieurs
            h=1/(N+1);%Pas de discrétisation en espace
            x=h:h:1-h;
            % les matrices
            %dérivée seconde
            A=gallery('tridiag',N)/h^2;
            
            % Partie implicite  du schéma
            
            MI=(1+0.5*mu*dt)*speye(N,N)+0.5*c^2*dt^2*A;
            ME=(0.5*mu*dt-1)*speye(N,N)-0.5*c^2*dt^2*A;


            
        %
        % 3) Initialisation des itérations en temps
        %

            U0=u0(x)';
            U1=U0;%la vitesse initiale est nulle
            
            figure(1)
            %Energie discrète
            EN=[];
            for k=2:Npdt
                %
                % un pas en temps
                %
                F=2*U1+ME*U0;
                U=MI\F;
                U0=U1;
                U1=U;
                %
                % représentation graphique périodique
                %
                if mod(k,Nper)==0
                    plot(x,U)
                    axis([0 1 -1 1])
                    drawnow
                    %pause(0.2)
                    e=0.5*((norm(U1-U0)/dt)^2+c^2*0.5*((A*U1)'*U1+(A*U0)'*U0));
                    EN=[EN e];
                end
                
            end
            figure(2)
            semilogy(EN)
