// Equation de la chaleur - schema de Crank-Nicholson (theta schema avec theta=1/2)
//clear; close,
//da=gda(); da.auto_clear='on';
//
V=1;	// vitesse de convection
T=0.01;	// temps final
theta=0.5;	
cfl=1.;	
N=[50,100,200,400,800,1600]';	// Nb de point de discretisation

erreur=zeros(N);
for i=1:size(N,'r')
	dx=1/N(i);
// on peut se placer a cfl fixée (ordre 2 en espace) 
//	dt=cfl*dx^2/2. 
// ou faire dependre lineairement dt de dx (ordre 2 en temps)
	dt=dx;
	x=linspace(0,1,N(i)+1);
//	u0=(1.-4.*((x-0.5)^2))';
	u0=cos(4*%pi*x)';
	u=u0;
	plot(x,u);
	A=diag(sparse(2.*ones(1,N(i)+1)))-diag(sparse(ones(1,N(i))),1)-diag(sparse(ones(1,N(i))),-1);
	A(1,2)=-2.;
	A(N(i)+1,N(i))=-2.;
	A=A/(dx)^2;
	B=speye(N(i)+1,N(i)+1)+theta*dt*A;
	B2=speye(N(i)+1,N(i)+1)-(1.-theta)*dt*A;
	for t=dt:dt:T,
		v=B2*u;
		u=B\v;// theta-schema
		//	plot2d(x,[u u0]);
	end;
	uexact=exp(-16*%pi^2*T)*cos(4*%pi*x)';
	erreur(i)=sqrt(((u-uexact)'*(u-uexact)*dx));
	plot2d(x,[u u0 uexact]);
scf
end;
dN=ones(size(N,'r'),1)./N;
dN2=ones(size(N,'r'),1)./(N^2);
plot2d(N,[erreur dN  dN2],logflag='ll');

