clear all
close all
%------------------------------
% Solution of the 1d advection equation
% u_t+c*u_x=0
% c-velocity;
% Initial condition: u(x,0)=u_0;
% x in (0,a) with
% periodic boundary conditions: u(0,t)=u(a,t).
%-----------------------------
% Analytical solution:
% u(x,t)=u0(x-ct)
%-----------------------------
%-----------------------------
% Lax-Wendroff method
% u_ij+1=b1*u_i-1^j+b2*u_i^j+b3*u_i+1^j;
% b1=h(h+1)/2;
% b2=1-h^2;
% b3=-h(1-h)/2;
% h=c*dt/dx;
%-----------------------------
c=0.5;
a=10;
dx=0.01;
dt=0.01;
T=800;
h=(c*dt)/(dx);
n=a/dx;
tt=100;
b1=h*(h+1)/2;
b2=1-h^2;
b3=-h*(1-h)/2;
c_zahl=h
%--------------
x=0:dx:a;
u0=exp(-10*(x-2).^2);%0.5*sin(2*pi*x/a);%exp(-(x-2).^2); % Initial condition
plot(x,u0,'r','LineWidth',2)
hold on
axis([0 a 0 1.1])
% text(1.8,1.01,'t=0',...
%      'HorizontalAlignment','left')
u(1,:)=u0;
%---------------
ah=b2*ones(n+1,1);
ao=b3*ones(n,1);
au=b1*ones(n,1);
A=diag(ah)+diag(au,-1)+diag(ao,1);
A(1,n+1)=b1;
A(n+1,1)=b3;
%----------------
for t=1:T
    um=A*u(t,:)'; % 
    u(t+1,:)=um';
    if mod(t,tt)==0
       plot(x,u(t+1,:),'b','LineWidth',2);
       title(['t=' num2str(t)],'Color','b')
       axis([0 a 0 1.1])
       hold on
       pause(0.5)
       plot(x,u(t+1,:),'r','LineWidth',2);
       xlabel('x','FontName','Times New Roman','FontSize',16)
       ylabel('u(x)','FontName','Times New Roman','FontSize',16)
    end;
end;

figure(2)
       plot(x,u(1,:),'b','LineWidth',2)
       hold on
%        plot(x-c*dt*50,u(50,:),'k','LineWidth',2)
%        hold on
       plot(x-c*dt*400,u(400,:),'b--','LineWidth',2)
       hold on
%        plot(x-c*dt*150,u(150,:),'k','LineWidth',2)
%        hold on
       plot(x-c*dt*800,u(800,:),'b:','LineWidth',2)
      axis([0 4 0 1.1])
      legend('t=0','t=400','t=800')
      xlabel('x-ct','FontName','Times New Roman','FontSize',16)
       ylabel('u(x)','FontName','Times New Roman','FontSize',16)

y=0:0.01:1;
      figure(3)
      plot(y, real(i*(log(1+1^2*(cos(pi*y)-1)-i*1*sin(pi*y)))/1),'b','LineWidth',2 )
      hold on
      plot(y, real(i*(log(1+(0.5)^2*(cos(pi*y)-1)-i*0.5*sin(pi*y)))/0.5),'b--','LineWidth',2 )
      hold on
      plot(y, real(i*(log(1+(0.75)^2*(cos(pi*y)-1)-i*0.75*sin(pi*y)))/0.75),'b-.','LineWidth',2 )
      hold on
      plot(y, real(i*(log(1+(0.5)*(cos(pi*y)-1)-i*(1/sqrt(2))*sin(pi*y)))/(1/sqrt(2))),'b:','LineWidth',2 )
      xlabel('k\Delta x/\pi')
      ylabel('\omega_1\Delta t/\alpha')
      legend('\alpha=1','\alpha=0.5','\alpha=0.75', '\alpha=1/\surd2')

       figure(4)
      plot(y, -imag(i*(log(1+(0.5)^2*(cos(pi*y)-1)-i*0.5*sin(pi*y)))/0.5),'b--','LineWidth',2 )
      hold on
      plot(y, -imag(i*(log(1+(0.75)^2*(cos(pi*y)-1)-i*0.75*sin(pi*y)))/0.75),'b-.','LineWidth',2 )
      hold on
      plot(y, -imag(i*(log(1+(0.5)*(cos(pi*y)-1)-i*(1/sqrt(2))*sin(pi*y)))/(1/sqrt(2))),'b:','LineWidth',2 )
      axis([0 1 0 6])
      xlabel('k\Delta x/\pi')
      ylabel('\omega_2\Delta t/\alpha')
      legend('\alpha=0.5','\alpha=0.75', '\alpha=1/\surd2')
      