%   Author: Reija Ristilä <reija@hut.fi>
%   Last modified: 1999-03-01 16:03:56 EET

% define sche, T and k first
% sche=[1 1 1 1], T=10, k=1/4

% calculating diffusion for a square with zero borders
% comparision of different numerical schemes

n=6; h=1/2; 
U0 = zeros(n^2,1); U0([15 16 21 22]) = [ 1 1; 1 1];
A = 1/h^2*lapm(n);
I = eye(n^2,n^2);

% Dirichlet borders
top = zeros(1,n); left = top; right = top; bot = top; 
a1=[top;zeros(n-2,n);bot]; 
a2=[left' zeros(n,n-2) right'];
aa=(a1+a2)'; 
f=aa(:); 
NW=0.5*(left(n)+top(1));NE=0.5*(top(n)+right(n));  
SE=0.5*(bot(n)+right(1));SW=0.5*(bot(1)+left(1));

figure(4),clf,hold on

% U' = AU+f, U vektori
% Euler explicit
if sche(1)==1
  ve=zeros(1,T-1);
  Uold = U0;
  figure(1),clf
  'Explicit Euler', pause
  surfc([NW top NE; left' reshape(U0,n,n)' right'; SW bot SE;])
  title('Explicit')
  pause
  for i = 2:T,
    Unew = (I+k*A)*Uold +k*f;
    surfc([NW top NE; left' reshape(Unew,n,n)' right'; SW bot SE;])
    Uold=Unew;
    fprintf('%6.2f ',max(Unew))
    ve(i-1)=max(Unew);
    pause
  end
  fprintf('\n')
  Uee = reshape(Unew,n,n)'; title('Explicit')
  figure(4),semilogy(ve,'g--p')
end

if sche(2) == 1
  vi=zeros(1,T-1);
  Uold = U0; 
  figure(2),clf
  'Implicit Euler', pause
  surfc([NW top NE; left' reshape(U0,n,n)' right'; SW bot SE;])
  title('Implicit')
  pause
  A1 = (I-k*A);
  for i = 2:T,
    Unew = A1 \ (Uold +k*f);
    surfc([NW top NE; left' reshape(Unew,n,n)' right'; SW bot SE;])
    Uold=Unew;
    fprintf('%6.3f ',max(Unew))
    vi(i-1)=max(Unew);
    pause
  end
  fprintf('\n')
  Uie = reshape(Unew,n,n)'; title('Implicit')
  figure(4),semilogy(vi,'-o')%,legend('Impl',1)
end

if sche(3) == 1
  vt=zeros(1,T-1);
  Uold = U0; 
  figure(3),clf
  'Trapezoidal rule', pause
  surfc([NW top NE; left' reshape(U0,n,n)' right'; SW bot SE;])
  title('Trapezoidal Rule')
  pause
  A1 = (I-k/2*A); A2 = (I+k/2*A);
  for i = 2:T,
    Unew = A1 \ (A2*Uold +k*f);
    surfc([NW top NE; left' reshape(Unew,n,n)' right'; SW bot SE;])
    Uold=Unew;
    fprintf('%6.3f ',max(Unew))
    vt(i-1)=max(Unew);
    pause
  end
  fprintf('\n')
  Utr = reshape(Unew,n,n)'; title('Trapezoidal Rule')
  figure(4),semilogy(vt,'r-^')%,legend('Trap',2)
end

% For Euler Explicit to be stable, k < 2/abs(min(eig(A)))
%k=2/abs(min(eig(A)))
%sche=[1 1 0 0];

% Trapezoidal rule is more stable, if k is smaller
%k=1/4/2
%%sche=[0 1 1 0];

if sche(4) == 1
  m=4; vit=zeros(1,T-1);
  Uold = U0; 
  figure(1),clf
  %a few steps with Implicit, then Trapezoidal
  'Implicit + Trapezoidal'
  surfc([NW top NE; left' reshape(U0,n,n)' right'; SW bot SE;])
  title('Implicit + Trapezoidal Rule')
  pause
  A1 = (I-k*A);
  for i = 2:m,
    Unew = A1 \ (Uold +k*f);
    surfc([NW top NE; left' reshape(Unew,n,n)' right'; SW bot SE;])
    Uold=Unew;
    fprintf('%6.3f ',max(Unew))
    vit(i-1)=max(Unew);
    pause
  end
  'Switching to Trapezoidal Rule'
  pause
  A1 = (I-k/2*A); A2 = (I+k/2*A);
  for i = (m+1):T,
    Unew = A1 \ (A2*Uold +k*f);
    surfc([NW top NE; left' reshape(Unew,n,n)' right'; SW bot SE;])
    Uold=Unew;
    fprintf('%6.3f ',max(Unew))
    vit(i-1)=max(Unew);
    pause
  end
  fprintf('\n')
  Uit = reshape(Unew,n,n)';
  figure(4),semilogy(vit,'k-x'),%legend('Impl + Trap',3)
end
pause
figure(4),legend

% semilogy(vi,'b-o'),hold on,semilogy(vt,'r-^'),semilogy(vit,'k-x')

