clear all;
%   Programma che calcola la soluzione di
%   w_t+alfa w_x=0 in [0,1]x[0,T]
%   w(x,0)=uinit(x)
%   con vari metodi
%
%   la soluzione viene calcolata in w usando il precedente valore w0
%
%   la soluzione esatta e' u(x,t)=pulse(x-alfa t)
%
%   icase=1 --> pulse=sin, icase=2 --> pulse funzione a campana
%
icase=input('Che dato iniziale? (1=sin, 2=campana) ');
%
%   scelta del metodo
%
%  imetodo=1  --> differenze centrate
%  imetodo=2  --> differenze upwind
%
imetodo=input('Con che metodo? (1=diff. centrate, 2=diff. upwind, 3=LW) ');
%
T=1;
N1=input('Quanti intervalli spaziali? ');
N=N1-1;
Nstep=input('Quanti intervalli temporali? ');
%N=199;
%Nstep=200;
h=1/N1;
dt=T/Nstep;
dt2=dt^2;
lambda=dt/h;
alfa=1;
la=lambda*alfa;
la2=la^2;
figure(1);
clf;
set(1,'position',[407 12 612 684]);
set(1,'defaultlinelinewidth',3);
%
%
%  costruzione dati iniziali
%
%  w(x,0)
%
xxx=linspace(0,1,200);
xx=linspace(0,1,N+2);
x=linspace(h,1,N+1);
w0=pulse(x,icase);
plot([0 x],[pulse(0,icase) w0])
axis([0 1 -1.1 1.1])
pause(0.4)
%
%  ciclo temporale
%
if imetodo==1
%
t=0;
for m=1:Nstep
  w(1)=w0(2)-la*(w0(1)-pulse(-alfa*t,icase));
  for i=2:N
    w(i)=w0(i)-la*(w0(i+1)-w0(i-1))/2;
  end
  w(N+1)=w0(N+1)-la*(w0(N+1)-w0(N));
  t=t+dt;
  y=pulse(xx-alfa*t,icase);
  plot(xx, [pulse(-alfa*t,icase) w],'*--',xx,y,'-')
  axis([0 1 -1.1 1.1])
  pause(0.2)
  w0=w;
end
%
elseif imetodo==2
%
t=0;
for m=1:Nstep
  w(1)=w0(1)-la*(w0(1)-pulse(-alfa*t,icase));
  for i=2:N+1
    w(i)=w0(i)-la*(w0(i)-w0(i-1));
  end
  t=t+dt;
  y=pulse(xx-alfa*t,icase);
  yy=pulse(xxx-alfa*t,icase);
  plot(xxx,yy,'r')
  hold on;
  plot(xx, [pulse(-alfa*t,icase) w],'*--',xx,y,'-')
  hold off;
  axis([0 1 -1.1 1.1])
  pause(0.2)
  w0=w;
end
%
elseif imetodo==3
%
t=0;
for m=1:Nstep
w(1)=w0(1)-la*(w0(2)-pulse(-alfa*t,icase))/2+ ...
     la2*(w0(2)-2*w0(1)+pulse(-alfa*t,icase))/2;
for i=2:N
  w(i)=w0(i)-la*(w0(i+1)-w0(i-1))/2+la2*(w0(i+1)-2*w0(i)+w0(i-1))/2;
end
dw=(3*w0(N+1)-4*w0(N)+w0(N-1))/2;
d2w=(-5*w0(N+1)+12*w0(N)-9*w0(N-1)+2*w0(N-2))/6;
w(N+1)=w0(N+1)-la*dw+la2*d2w/2;
t=t+dt;
y=pulse(xx-alfa*t,icase);
yy=pulse(xxx-alfa*t,icase);
plot(xxx,yy,'r')
hold on;
plot(xx, [pulse(-alfa*t,icase) w],'*--',xx,y,'-')
hold off;
axis([0 1 -1.1 1.1])
pause(0.4)
%
w0=w;
end
%
end
