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

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

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

n=6; h=1/2; 
U0 = ones(n^2,1); 
I = eye(n^2,n^2);

% Dirichlet borders; applied for left and right
top = zeros(1,n); bot = top; left=ones(1,n); right=left; 
a1=[top;zeros(n-2,n);bot]; 
a2=[left' zeros(n,n-2) right'];
aa=(a1+a2)'; 
f=aa(:); 
% Neumann borders: no-flow -borders for top and bottom with exceptions
A=lapm(n);
x=[1:6 31:36]; y=[7:12 25:30]; g=0;
for i=1:2*n
  A(x(i),y(i))=2;
end
A = 1/h^2*A;
f(x)=2/h*g;
g=2; b=-5;
f([32:35])=2/h*g*[1 2 2 1];

figure(4),clf,hold on

% U' = AU+f, U vektori
% Euler explicit
if sche(1)==1
end

if sche(2) == 1
  vi=zeros(1,T-1);
  Uold = U0; 
  figure(2),clf
  'Implicit Euler', pause
  surfc([left' reshape(U0,n,n)' right']), title('Implicit')
  pause
  A1 = (I-k*A);
  for i = 2:T,
    f([3 4]) = 2*h*b*(Uold([3 4]).^4 +1);
    Unew = A1 \ (Uold +k*f);
    surfc([left' reshape(Unew,n,n)' right'])
    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([left' reshape(U0,n,n)' right']), title('Trapezoidal Rule')
  pause
  A1 = (I-k/2*A); A2 = (I+k/2*A);
  for i = 2:T,
    f([3 4]) = 2*h*b*(Uold([3 4]).^4 +1);
    Unew = A1 \ (A2*Uold +k*f);
    surfc([left' reshape(Unew,n,n)' right'])
    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

if sche(4) == 1
  m=3; vit=zeros(1,T-1);
  Uold = U0; 
  figure(1),clf
  %a few steps with Implicit, then Trapezoidal
  'Implicit + Trapezoidal'
  surfc([left' reshape(U0,n,n)' right']), title('Implicit + Trapezoidal Rule')
  pause
  A1 = (I-k*A);
  for i = 2:m,
    f([3 4]) = 2*h*b*(Uold([3 4]).^4 +1);
    Unew = A1 \ (Uold +k*f);
    surfc([left' reshape(Unew,n,n)' right'])
    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,
    f([3 4]) = 2*h*b*(Uold([3 4]).^4 +1);
    Unew = A1 \ (A2*Uold +k*f);
    surfc([left' reshape(Unew,n,n)' right'])
    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


