xdel(winsid());


disp('Exercice 2')
alpha=0.9;

upwind=struct('v',[0,1-alpha,alpha],'nom','Upwind');
centre=struct('v',[-alpha/2,1,alpha/2],'nom','centré')
laxfriedrichs=struct('v',[(1-alpha)/2,0,(1+alpha)/2],'nom','Lax-Friedrichs')
laxwendroff=struct('v',[(-alpha+alpha^2)/2,1-alpha^2,(alpha+alpha^2)/2],'nom','Lax-Wendroff')


function[y]=heaviside(x)
  if length(x)==0 then
    y='H(x)';
    else
  y=zeros(x);
  I=find(x>1/2);
  y(I)=1;
  I=find(x<=1/2);
  y(I)=-1;
  end
endfunction


function[y]=heaviside2(x)
  if length(x)==0 then
   y='H2(x)';
    else
  y=zeros(x);
  I=find(x>1/2);
  y(I)=.908734643304;
  I=find(x<=1/2);
  y(I)=0.012938732984;
  end
endfunction

function[y]=sinusoide(x)
  if length(x)==0 then
        y='\sin(8\pi x)';

  else
         y=sin(4*2*%pi*x);
  end
endfunction  


function[y]=sinusoide2(x)
  if length(x)==0 then
        y='\sin(4\pi x)^2';

  else
         y=sin(2*2*%pi*x)^2;
  end
endfunction  


function[y]=fonctionreguliere(x)
  if length(x)==0 then
    y='|\sin(2\pi x)|^{6.2}';
  else
    y=(sin(2*%pi*x).^2).^(3.1);
  end
endfunction



function[]=schema_periodique(u0,schema)

N=100;
J=50;
x=(0:J-1)/J;
erreur=zeros(2,N);
u=u0(x);
f=figure();
f.pixmap='on';

for n=1:N

u=schema.v(2)*u+schema.v(3)*u([J,1:J-1])+schema.v(1)*u([2:J,1]);

clf;
subplot(2,1,1)
xmct=pmodulo(x-alpha*n/J,1);
I=pmodulo((1:J)+ceil(alpha*n)-1,J)+1;

plot2d(x,[u;u0(xmct)]');

legend(['Solution approchée','Solution exacte']);
title(schema.nom+' , n='+string(n)+', J='+string(J))
g=gca();
g.data_bounds=[0,1,-2,2];
xlabel('x-ct')
ylabel('u')
show_pixmap();
erreur(:,n)=[norm(u-u0(xmct),2)/(J^(1/2));norm(u-u0(xmct),%inf)];
end
f.pixmap='off';

subplot(2,1,2)
plot((1:n)/J,erreur')
xtitle('$\mbox{Erreur pour '+schema.nom+', } u_0(x)='+u0([])+'$')
legend('$\mbox{Erreur } l^2$','$\mbox{Erreur } l^\infty$')
xlabel('t')
endfunction

schema_periodique(sinusoide,upwind);
schema_periodique(fonctionreguliere,upwind);
schema_periodique(heaviside,upwind);

schema_periodique(sinusoide,centre);

schema_periodique(sinusoide,laxfriedrichs);
schema_periodique(fonctionreguliere,laxfriedrichs);
schema_periodique(heaviside,laxfriedrichs);

schema_periodique(sinusoide,laxwendroff);
schema_periodique(fonctionreguliere,laxwendroff);
schema_periodique(heaviside,laxwendroff);




disp('Exercice 3')


function[]=schema_bords(u0,schema)
//On impose u_0^n=0 et u_(J+1)^n=u_J^n.
N=60;
J=50;
x=(0:J-1)/J;
erreur=zeros(2,N);
u=u0(x);

f=figure();
f.pixmap='on';

for n=1:N

u=schema.v(2)*u+schema.v(3)*[0,u([1:J-1])]+schema.v(1)*u([2:J,J]);

clf;
subplot(2,1,1)
xmct=max(x-alpha*n/J,0);

plot2d(x,[u;u0(xmct)]');

legend(['Solution approchée','Solution exacte']);
xtitle(schema.nom+' , n='+string(n)+', J='+string(J)+', Dirichlet en 0, Neumann en 1')
g=gca();
g.data_bounds=[0,1,-2,2];
xlabel('x')
ylabel('u')
show_pixmap();
erreur(:,n)=[norm(u-u0(xmct),2)/(J^(1/2));norm(u-u0(xmct),%inf)];
end
f.pixmap='off';
subplot(2,1,2)
plot2d((1:n)/J,erreur')
xtitle('$\mbox{Erreur pour '+schema.nom+', } u_0(x)='+u0([])+'$')
legend('$\mbox{Erreur } l^2$','$\mbox{Erreur } l^\infty$')
xlabel('t')
endfunction
//
schema_bords(sinusoide2,upwind)
schema_bords(sinusoide2,laxfriedrichs)
schema_bords(sinusoide2,laxwendroff)



disp('Exercice 3')


function[]=schema_implicite(u0)
N=60;
J=100;
x=(0:J-1)/J;
erreur=zeros(2,N);
u=u0(x);

//D est une matrice creuse telle que u*D=u([J,1:J-1])
e=sparse(ones(J-1,1));
D=diag(e,1);
D(J,1)=1;

// Le schema, avec conditions aux bords periodiques peut s ecrire
// sous la forme Au^(n+1)=u^n.
A=speye(J,J)-alpha*D/2+alpha*D'/2;

f=figure();
f.pixmap='on';
for n=1:N


u=u/A;


clf;
subplot(2,1,1)
xmct=pmodulo(x-alpha*n/J,1);
I=pmodulo((1:J)+ceil(alpha*n)-1,J)+1;

plot2d([xmct(I);x]',[u(I);u0(x)]');

legend(['Solution approchée','Solution exacte']);
title('Schéma implicite centré , n='+string(n)+', J='+string(J))
g=gca();
g.data_bounds=[0,1,-2,2];
xlabel('x-ct')
ylabel('u')
show_pixmap()
erreur(:,n)=[norm(u-u0(xmct),2)/(J^(1/2));norm(u-u0(xmct),%inf)];
end
f.pixmap='off';
subplot(2,1,2)
plot2d((1:n)/J,erreur')
xtitle('$\mbox{Erreur pour le schema implicite, } u_0(x)='+u0([])+'$')
legend('$\mbox{Erreur } l^2$','$\mbox{Erreur } l^\infty$')
xlabel('t')
endfunction

schema_implicite(fonctionreguliere);


disp('Exercice 4')

function[]=schema_second_membre(u0)
  b=1;
N=100;
J=100;
dx=1/J;
dt=dx;
c=alpha*dx/dt;
bdt=b*dt;

x=(0:J-1)/J;
erreur=zeros(2,N);
u=u0(x);
f=figure();
f.pixmap='on';
for n=1:N
u=(1+alpha+bdt)*u([J,1:J-1])/2+(1-alpha+bdt)*u([2:J,1])/2;

clf;
subplot(2,1,1)
xmct=pmodulo(x-alpha*n/J,1);
I=pmodulo((1:J)+ceil(alpha*n)-1,J)+1;

plot2d([xmct(I);x]',[u(I);u0(x)*exp(bdt*n)]');

legend(['Solution approchée','Solution exacte']);
title('Transport avec second membre, Lax-Friedrichs, n='+string(n)+', J='+string(J))
g=gca();
g.data_bounds=[0,1,-2*exp(bdt*n),2*exp(bdt*n)];
xlabel('x-ct')
ylabel('u')
show_pixmap();
erreur(:,n)=[norm(u-u0(xmct)*exp(bdt*n),2)/(J^(1/2));norm(u-u0(xmct),%inf)];
end
f.pixmap='off'
subplot(2,1,2)
plot((1:n)/J,erreur')
xtitle('$\mbox{Erreur pour Transport avec second membre, Lax-Friedrichs, } u_0(x)='+u0([])+'$')
legend('$\mbox{Erreur } l^2$','$\mbox{Erreur } l^\infty$')
xlabel('t')
endfunction

schema_second_membre(fonctionreguliere);


disp('Exercice 5')


function[y]=F(x)
  y=x.*(1-x);
endfunction

function[y]=DF(x)
  y=1-2*x;
endfunction


function[]=schema_trafic_routier(u0)
  b=1;
N=200;
J=200;
dx=1/J;

dt=alpha*dx;
bdt=b*dt;

x=(0:J-1)/J;
erreur=zeros(2,N);
u=u0(x);
f=figure();
f.pixmap='on';
for n=1:N
  // Calcul du flux de voitures.
  phi=zeros(u);
  up=u([2:J,1]);
  
  phi=(F(u)+F(up))/2;
  u=(up+u([J,1:J-1]))/2+(dt/dx)*(phi([J,1:J-1])-phi);

clf;
xmct=pmodulo(x+alpha*n/J*DF(u0(x)),1);

plot2d([x;xmct]',[u;u0(x)]',style=[2,0]);
legend(['Solution approchée','Solution exacte (tant que les caractéristiques ne se sont pas croisées)']);
title('Trafic routier, n='+string(n)+', J='+string(J))
g=gca();
g.data_bounds=[0,1,-1,2];
xlabel('x')
ylabel('u')
show_pixmap();
end
f.pixmap='off';
endfunction

alpha=0.6
schema_trafic_routier(fonctionreguliere);
schema_trafic_routier(heaviside2);