//---------------------------------------
// Vibration d'une corde. 
// Schéma aux différences finies
// comparaison avec approximation par modes propres
//----------------------------------------
/////////////////////////////////////////////
//paramètres
L=1;
a=0.7;b=0.8;
N=60;
M=610;
omega=3/2*(%pi/L);
Tfin=10;
dx=L/(N+1);
dt=Tfin/(M+1);
umin=-0.05;umax=0.05 //pour la fenêtre de tracé
////////////////////////////////////////////////
//initialisation vecteurs
u=zeros(N,1);//On place ici les données calculées pour le temps n+1
ua=zeros(N,1);//positions au temps n
ub=zeros(N,1);//positions au temps n-1
Fext=zeros(N,1);
////////////////////////////////////
// matrice laplacien
/////
A=toeplitz([2*(1-(dt^2)/(dx^2)),(dt^2)/(dx^2),zeros(1,N-2)]);//
////////////////////////////////////////
//calcul et affichage
figure(1)
xtitle('vibrations de la corde')
//
for i=1:(M+1);
    for j=1:N
      if ((j*dx)<b)&(a<(j*dx)) then  
      Fext(j)=cos(omega*i*dt)
      end
    end  
    u=A*ua-ub+(dt)^2*Fext;
    drawlater
    clf()
    plot2d([0:dx:1],[0;u;0],rect=[0,umin,L,umax])
    drawnow
    ub=ua;
    ua=u;
end
//////////////////////////////////////////////////////////
// decompostion suivant les modes propres
///////////////////////////////////////////////////////////
gk=zeros(20,1)
//
for i=1:20
  gk(i)=sqrt(2)/(i*%pi)*(cos(i*%pi*a)- cos(i*%pi*b))
end
/////////////////
valu3=zeros(N,1);valu20=zeros(N,1);
for j=1:N
 u3=0;u20=0;
 for i=1:20
    u20=u20+gk(i)/((i*%pi)^2-(3*%pi/2)^2)*(cos(3*%pi/2*Tfin)-cos(i*%pi*Tfin))*sqrt(2)*sin(i*%pi*j*dx)
  end
  for i=1:3
    u3=u3+gk(i)/((i*%pi)^2-(3*%pi/2)^2)*(cos(3*%pi/2*Tfin)-cos(i*%pi*Tfin))*sqrt(2)*sin(i*%pi*j*dx)
  end
  valu3(j)=u3;valu20(j)=u20;
end
figure(2)
clf()
plot2d([0:dx:1],[0;u;0],1,rect=[0,umin,L,umax])  
plot2d([0:dx:1],[0;valu3;0],2,rect=[0,umin,L,umax])  
plot2d([0:dx:1],[0;valu20;0],3,rect=[0,umin,L,umax])
legends(['différences finies','approximation 3 modes','approximation 20 modes'],[1,2,3])
