% preqcn.m
% fd, iacs, mai, 08 juillet 2004, 2 janvier  2005, 25 fevrier 2025 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xlabel ('preqcn. Precision du schema de Crank Nicolson pour un retour a l equilibre') ;
% title ('Cnam Paris. Travaux pratiques Iacs.') ;
m=9
pas=zeros(m,1);
err=zeros(m,1);
%
  n=8;
  for j=1:m;
dt=1/n;
t=zeros(n+1,1);
u=zeros(n+1,1);
u;
u(1:1)=1;
t(1:1)=0;
  for i=1:n;
t(i+1)=t(i)+dt;
u(i+1)= (1-dt/2)*u(i) / (1+dt/2) ;
end ;
% plot(t,u(:,1),'*')
pas(j)=log(n)/log(2);
err(j)=log ( abs(u(n+1,1) - exp(-1)) ) / log(2) ;
n=2*n; end;

x=zeros(m,1); x =pas ;
y=zeros(m,1); y= err(:,1);
xm=0;
ym=0;
x2m=0;
y2m=0;
xym=0;
  for j=1:m;
xm=xm+x(j);
ym=ym+y(j);
x2m=x2m+x(j)*x(j);
y2m=y2m+y(j)*y(j);
xym=xym+x(j)*y(j);
end ;
xm=xm / m ;
ym=ym / m ;
x2m=x2m / m ;
y2m=y2m / m ;
xym=xym / m ;
alpha = ( xym -xm*ym ) / (x2m - xm*xm )


plot(pas,err(:,1),'*;logarithme des ereurs;') ; 
axis ([ 2.7 ,  11.3 ,  -27.6 , -8.9 ]) 
title('retour equilibre        schema de Crank-Nicolson : ordre 2')
print -dpdf 09-preqcn.pdf
print -djpg 09-preqcn.jpg

