% Einfuehrung in die Numerik 1  SS12
% Bsp_5 Blatt_2

% Feichtel Ernst
% 0630363


% Kosten zur Loesung des linearen Gleichungssystems 
% mit Gaußscher Elimination


clear
clc
clf


N=20;
N1 = N+1;
N4 = round(N/4);

% verrauschte Daten v

ustar = [zeros(N4,1);ones(N1-2*N4,1);zeros(N4,1)];
ustar = kron(ustar,transpose(ustar));
v = ustar + 0.1*randn(N1,N1);

N1N1 = N1*N1;
v = reshape(v,N1N1,1);


% Matrix D (mit Eintraegen D0, D1, D0) bilden
 
Dxdiag = -ones(N1,N1); Dxdiag(N1,:) = 0;
Dxsup = ones(N1,N1); Dxsup(1,:) = 0;
Dx = spdiags(Dxdiag(:), 0,N1N1,N1N1)+ spdiags( Dxsup(:), 1,N1N1,N1N1);
Dydiag = -ones(N1,N1); Dydiag(:,N1) = 0;
Dysup = ones(N1,N1); Dysup(:,1) = 0;
Dy = spdiags(Dydiag(:), 0,N1N1,N1N1)+ spdiags( Dysup(:),N1,N1N1,N1N1);
D = transpose(Dx)*Dx + transpose(Dy)*Dy;

% Gleichungssystem erstellen

h=1/N;
muy=0.002;
A=(muy/h^2)*D+eye(N1*N1);
Av=[A,v];
n=N1*N1;

% Gauss, Vorwaerts Substitution
 
for k=1:(n-1)
   for i=(k+1):n
       muv(i,k)=Av(i,k)/Av(k,k);
       %Av(i,k:(n+1))=Av(i,k:(n+1))-mv(i,k)*Av(k,k:(n+1));
       for j=k:(n+1)
           Av(i,j)=Av(i,j)-muv(i,k)*Av(k,j);
       end
   end
end

% Gauss, Rueckwaerts Substitution

u=zeros(n,1);
u(n)=Av(n,n+1)/Av(n,n);
for i=(n-1):(-1):1
   s=0;
   for j=1:n
       if (i+j) < (n+1)
           s=s+Av(i,i+j)*u(i+j);
       end
   end
   u(i)=(Av(i,n+1)-s)/Av(i,i);
end

% u,v in Bilder umwandeln

u = reshape(u,N1,N1);
v = reshape(v,N1,N1);

% ustar, u, v graphisch darstellen

h1 = figure(1); clf; set(h1,'Position',[10 10 900 300]);
subplot(1,3,1);
imagesc(ustar); colormap('gray'); axis image; axis off;
title('exakt');
subplot(1,3,2);
imagesc(v); colormap('gray'); axis image; axis off;
title('verrauscht');
subplot(1,3,3);
imagesc(u); colormap('gray'); axis image; axis off;
title('entrauscht');

% Waehlt man N > 35 steigt die Laufzeit des Programmes an,
% im Vergleich zu (8 < N < 35)
