function final2(nu,anim,nbite) close all; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Constant values % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nbpoint=100; L=1; lambda=1; deltat=0.00001; deltax=L/nbpoint; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Calcul of the control parameter % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% D=lambda/(L*L*nu); s=num2str(D); strcat('la valeur de D est : ',s) n=1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Array initialisation % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:nbpoint+4 sol(i,1)=0; u(i,1)=0; nsol(i,1)=0; for j=1:nbpoint+4 mat(i,j)=0; A(i,j)=0; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Type of perturbation % % f->file sol.txt % % something else ->perturbation % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R=input('Quel type de perturbation ?','s') if R=='f' fid=fopen('sol2.txt','rt'); for i=1:nbpoint+4 resul=fscanf(fid,'%f %f',2); abscisse(i)=resul(1); sol(i,1)=resul(2); end fclose(fid); else for i=1:nbpoint+4 abscisse(i)=(i-3)*deltax; sol(i,1)=0.1*sin(n*2*pi*abscisse(i)); end %sol(50,1)=1; %sol(49,1)=1; %sol(51,1)=1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Save of the first picture of the animation % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% M=moviein(anim+1); plot(abscisse,sol,'r'); %axis([0 1 -0.1 1]); grid; az=num2str(0); nom=strcat('im',az,'.jpg'); %print(1,'-djpeg90',nom); M(:,1)=getframe; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Begin of the time iteration % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=1:nbite % Fill of the matrix which have to be inversed mat(1,1)=1; mat(1,nbpoint+1)=-1; mat(2,2)=1; mat(2,nbpoint+2)=-1; for i=3:nbpoint+2 mat(i,i)=1-2*nu*deltat/(deltax*deltax)+6*lambda*deltat/(deltax*deltax*... deltax*deltax); mat(i,i-1)=nu*deltat/(deltax*deltax)-4*lambda*deltat/(deltax*deltax*... deltax*deltax); mat(i,i-2)=lambda*deltat/(deltax*deltax*deltax*deltax); mat(i,i+1)=nu*deltat/(deltax*deltax)-4*lambda*deltat/(deltax*deltax*... deltax*deltax); mat(i,i+2)=mat(i,i-2); end mat(nbpoint+3,nbpoint+3)=1; mat(nbpoint+3,3)=-1; mat(nbpoint+4,nbpoint+4)=1; mat(nbpoint+4,4)=-1; % Fill of the matrix represeted the explicit part of the modelisation for i=3:nbpoint+2 u(i,1)=sol(i,1)-1/2*deltat*((sol(i+1,1)*sol(i+1,1)-sol(i-1,1)*sol(i-1,1))... /(4*deltax)+sol(i,1)*... (sol(i+1,1)-sol(i-1,1))/(2*deltax)); end % Inversion of matrix determination of the solution inm=inv(mat); nsol=inm*u; sol=nsol; % Selection of the picture if mod(j,fix(nbite/anim))==0 plot(abscisse,sol,'r'); %axis([0 1 -0.1 1 ]); grid; az=num2str(j); nom=strcat('im',az,'.jpg'); % print(1,'-djpeg90',nom); M(:,j/fix(nbite/anim))=getframe; end sol(1,1)=0; sol(2,1)=0; sol(nbpoint+4,1)=0; sol(nbpoint+3,1)=0; te(j)=deltat*j; no(j)=norme(sol,nbpoint,deltax); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Post processoring % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Evolution of the norm of the solution figure; plot(te,no,'r.') % Save of the solution to restart the simulation R1=input('enregistrement ?','s'); if R1=='o' fid=fopen('sol2.txt','wt'); for i=1:nbpoint+4 fprintf(fid,' %f %f \n',[abscisse(i),nsol(i)]); end fclose(fid); end function y=norme(vec,int,deltax) y=0; for i=3:int+2 y=y+deltax*vec(i,1)*vec(i,1); end y=sqrt(y);