next up previous contents
Next: About this document ... Up: Appendix Previous: Appendix   Contents

Program


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);


Julien Delbove
2000-11-23