%%Lance Simms                                               2/13/06
%%AA214B-MacCormack
%
%This Script will attempt to solve the model hyperbolic equation 
%
%       du    du
%       -- +c -- = 0
%       dt    dx
%
%using several implicit time marching methods.
%
%The initial condition is specified and the wave speed is assumed to be
%positive

%The number of mesh points and so on
I=41;
delx=2/(I-1);
i=1:I;
x=(i-1).*delx;
cfl=.9;

%Initial condition has first half of wave at 1 and second half at 1/2
uo=zeros(1,I);
uo(i(find(i < I/4+1)))=1;
uo(i(find(i > I/4+1)))=1/2;

nsteps=10;
%METHOD 1-----------------------Non-linear backwards difference
u=uo;
delt=delx;

for i=1:nsteps
    f=u.^2./2;                  %u du/dx=.5* d(u^2)/dx
    u(2:I)=u(2:I)-(delt./delx).*(f(2:I)-f(1:I-1));
    plot(x,u);
end

%METHOD 4-----------------------Central Implicit Method
u=uo;
alpha=1;
a=zeros(1,I);   a(1:I-1)=1;      a(I)=1+cfl;
b(1:I-2)=-.5*alpha*cfl;          b(I-1)=-cfl;
c(2:I-1)=.5*alpha*cfl;          c(1)=0;
[L,U,alph,gamm]=LU_decomposition(a,b,c,I);
   
for i=1:nsteps
   f(1)=uo(1);                                 f(I)=u(I);
   f(2:I-1)=u(2:I-1)-.5*(1-alpha)*cfl.*(u(3:I)-u(1:I-2));
   u=LU_solve(alph,b,gamm,f,I);
end

subplot(2,1,1);
plot(x,u);
set(gca,'Ylim',[.25,1.25]);
text(.1,.9,'Central Implicit','units','normalized');
text(.8,.9,sprintf('CFL= %6.2f', cfl),'units','normalized');

%METHOD 5------------------------Crank-Nikolson
u=uo;
alpha=.5;
a=zeros(1,I);   a(1:I-1)=1;      a(I)=1+cfl;
b(1:I-2)=-.5*alpha*cfl;          b(I-1)=-cfl;
c(2:I-1)=.5*alpha*cfl;           c(1)=0;
[L,U,alph,gamm]=LU_decomposition(a,b,c,I);
    
for i=1:nsteps
   f(1)=uo(1);                              f(I)=u(I);
   f(2:I-1)=u(2:I-1)-.5*(1-alpha)*cfl.*(u(3:I)-u(1:I-2));
   u=LU_solve(alph,b,gamm,f,I);
end

subplot(2,1,2);
plot(x,u);
set(gca,'Ylim',[.25,1.25]);
text(.1,.9,'Crank-Nicolson Implicit','units','normalized');
text(.8,.9,sprintf('CFL= %6.2f', cfl),'units','normalized')