%=====================================================================
% Solving PDE's in spectral space
%
% Computer Lab for Math 221, Fall 2002
% Sorin Mitran 09/24/2002
%=====================================================================

% Clear previous work
clc; clf; clear;

% I. Advection equation
%
%  q_t + u q_x = 0
%  q(x,t=0)=q0(x) 0<=x<=2 pi
%  u=1
n=64; dx=2*pi/n;
x=(0:n-1)*dx+dx/2;
t=0.; tfinal=0.5*pi;
nmodesmax=2; u=1;

% Initial condition is built up as a sum of Fourier modes
for nmodes=1:nmodesmax
 q0=zeros(size(x)); a0=zeros([nmodes 1]); b0=zeros([nmodes 1]);
 for k=1:nmodes
   a0(k)=rand; b0(k)=rand;  
   q0 = q0 + a0(k) * cos(k*x) + b0(k) * sin(k*x);
 end
end

figure(1); plot(x,q0);
in_ax=axis;
disp('Initial condition. Strike a key');
title('Initial condition');
pause

cfl=0.5;

% Upwind scheme
q1=q0; qinit=q0;
nt=1; dt=cfl*abs(dx/u); sigma=u*dt/dx; I=sqrt(-1);
while t<tfinal  
  for i=2:n
    q1(i) = q0(i) - sigma*(q0(i)-q0(i-1));
  end
  q1(1) = q0(1) - sigma*(q0(1)-q0(n));   
  clc
  t=t+dt
  figure(1); plot(x,qinit,'.g',x,q1,'b');
  axis(in_ax); title('Upwind scheme');  
  pause(0.01)
  q0=q1;
end
qupwind=q1;
figure(1); plot(x,qinit,'.g',x,q1,'b');
axis(in_ax); title('Upwind scheme');  
clc;
disp('Upwind solution. Strike a key to continue ...');  
pause 

% Lax-Friedrichs scheme
q1=q0; t=0.
nt=1; dt=cfl*abs(dx/u); sigma=0.5*u*dt/dx; I=sqrt(-1);
while t<tfinal  
  for i=2:n-1
    q1(i) = 0.5*(q0(i+1)+q0(i-1)) - sigma*(q0(i+1)-q0(i-1));
  end
  q1(1) = 0.5*(q0(1)+q0(n)) - sigma*(q0(2)-q0(n));    
  q1(n) = 0.5*(q0(n-1)+q0(1)) - sigma*(q0(1)-q0(n-1));  
  clc
  t=t+dt
  figure(1); plot(x,qinit,'.g',x,q1,'b');
  axis(in_ax); title('Lax Friedrichs scheme');  
  pause(0.01)
  q0=q1;
end
qLaxF=q1;
figure(1); plot(x,qinit,'.g',x,q1,'b');
axis(in_ax); title('Lax Friedrichs scheme');  
clc;
disp('Lax Friedrichs solution. Strike a key to continue ...');  
pause

% Lax-Wendroff scheme
q1=q0; t=0.
nt=1; dt=cfl*abs(dx/u); sigma=0.5*u*dt/dx; sigma2=sigma*u*dt/dx; I=sqrt(-1);
while t<tfinal  
  for i=2:n-1
    q1(i) = q0(i) - sigma*(q0(i+1)-q0(i-1)) + sigma2*(q0(i+1)-2*q0(i)+q0(i-1));
  end
  q1(1) = q0(1) - sigma*(q0(2)-q0(n)) + sigma2*(q0(2)-2*q0(1)+q0(n));
  q1(n) = q0(n) - sigma*(q0(1)-q0(n-1)) + sigma2*(q0(1)-2*q0(n)+q0(n-1));
  clc
  t=t+dt
  figure(1); plot(x,qinit,'.g',x,q1,'b');
  axis(in_ax); title('Lax Wendroff scheme');  
  pause(0.01)
  q0=q1;
end
qLaxF=q1;
figure(1); plot(x,qinit,'.g',x,q1,'b');
axis(in_ax); title('Lax Wendroff scheme');  
clc;
disp('Lax Wendroff solution. Strike a key to continue ...');  
pause