%% Numerik Blatt 2 Beispiel 5
clear; clc; clf;

N=20;
h=1/N;
mu=0.012;
N1 = N+1;
N4 = round(N/4);
ustar = [zeros(N4,1);ones(N1-2*N4,1);zeros(N4,1)];
ustar = kron(ustar,ustar');
v = ustar + 0.1*randn(N1,N1);
N1N1=N1*N1;
v=reshape(v,N1N1,1);


%% Vektoren um die Matrix D mit spdiags zu erzeugen
d=3*ones(N-1,1);
d1 = [2;d;2];
N2= N1N1-2*(N+1);
d2 = 4*ones(N2,1);

for k=1:N2/(N+1)
    d2((k*N)+k) = 3;
end

for k=1:N2/(N+1)-1
    d2((k*N)+k+1) = 3;
end

d2(1) = 3;
D1  = [d1;d2;d1];                   % Vektor fuer die Hauptdiagonale


d1_neben1 = -1*ones(N1N1,1);
d2_neben = d1_neben1;                   % Vektor fuer Nebendiagonalen
d1_neben2 = d1_neben1;

for k=1:N1N1/(N+1)
    d1_neben1((N+1)*k) = 0;             % Vektor fuer die erste obere Nebendiagonale
end

for k=1:(N1N1-1)
    d1_neben2(k)=d1_neben1(k+1);        % Vektor fuer die erste untere Nebendiagonale
end

a=[D1 d1_neben2 d1_neben1 d2_neben d2_neben];
b=[0 1 -1 N+1 -(N+1)];
D = spdiags(a,b,N1N1,N1N1);

%% Gaussalgorithmus

A = mu*D/h^2 + speye(N1N1);
b = v;

A=full([A,b]);

m=size(A,1);                %Zeilen
n=size(A,2);                %Spalten
tic
for k=1:(m-1)
    for i=(k+1):m
        A(i,k) = A(i,k) / A(k,k);
        A(i,(k+1):n) = A(i,(k+1):n)-A(i,k)*A(k,(k+1):n);
%         for j=(k+1):n
%             A(i,j) = A(i,j) - A(i,k)*A(k,j);
%         end
    end
end

u=zeros(m,1);
u(m) = A(m,n) / A(m,m);
for i=(m-1):-1:1
    s=0;
    for j=(i+1):m
        s = s + A(i,j)*u(j);
    end
    u(i) = (A(i,m+1)-s) / A(i,i);
end
toc

%% Benoetigte Rechenzeit mit Gauss fuer
%     N=10        0.542427 seconds
%     N=20        19.015835 seconds
%     N=30        200.390454 seconds
%     N=50        3549.749848 seconds (fast eine Stunde!)
%       ->fuer große N extrem langsam (Groessenordnung der Rechenschritte O(n^3)!)
     
u = reshape(u,N1,N1);
v = reshape(v,N1,N1);

%% Graphische Darstellung
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');
