clc; clear; clf;

%% Beispiel 1_4

% Daten
N = 9; % N > 8
N1 = N+1;
N4 = round(N/4);
h = 1/N;

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);

mu = 0.012; % mu > 0

%% Matrix D
D0 = -ones(N1N1,1);
for k = 1:N1;
    D0(k*N1) = 0;
end
D_main = 4*ones(N1N1,1);
for k = [1 N1 N1N1-N1+1 N1N1];
    D_main(k) = 2;
end
for j = [2:N N1N1-N1+2:N1N1-1];
    D_main(j) = 3;
end
for i = 1:N1-2;
    D_main(N1*i + 1) = 3;
    D_main(N1*(i+1)) = 3;
end

D_leo = spdiags(D0,1,N1N1,N1N1) + spdiags(D0,-1,N1N1,N1N1) + spdiags(D0,N1,N1N1,N1N1) ...
    + spdiags(D0,-N1,N1N1,N1N1) + spdiags(D_main,0,N1N1,N1N1);

%% Gauß Elimination
G = (mu*D_leo/h^2 + speye(N1N1));
G = [G,v];

for k = 1:N1N1-1;
    for i = k+1:N1N1;
        M = G(i,k)/G(k,k);
        for j = k+1:N1N1+1;
            G(i,j) = G(i,j) - G(k,j)*M;
        end
    end
end

%% Rücksubstitution
u(N1N1) = G(N1N1,N1N1+1)/G(N1N1,N1N1);
for i = N1N1-1:-1:1;
    s = 0;
    for j = i+1:N1N1;
        s = s + G(i,j)*u(j);
    end
    u(i) = (G(i,N1N1+1)-s)/G(i,i);
end

%% Vergleichsdaten

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 = Dx'*Dx + Dy'*Dy;

u_2 = (mu*D/h^2 + speye(N1N1))\v;

u_2 = reshape(u_2,N1,N1);
u = reshape(u,N1,N1);
v = reshape(v,N1,N1);
%% Darstellung
h1 = figure(1); clf; set(h1,'Position',[10 10 900 300]);
subplot(1,4,1);
imagesc(ustar); colormap('gray'); axis image; axis off;
title('exact');
subplot(1,4,2);
imagesc(v); colormap('gray'); axis image; axis off;
title('verrauscht');
subplot(1,4,3);
imagesc(u); colormap('gray'); axis image; axis off;
title('entrauscht mit Gauß');
subplot(1,4,4);
imagesc(u_2); colormap('gray'); axis image; axis off;
title('entrauscht mit "/"');

% Bei großen N rechnet Matlab schon ziemlich lange. Der Gauß-Algorithmus
% müsste in seiner Arithmetik verbessert werden.
