function Krainz

clc; clear;

N=10;

%% Vom Zettel kopiert 
% v ... verauschte Daten
% ustar ... gewünschte Funktion

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

%% D konstruieren

D0=-speye(N+1);
D1=sparse(diag([2;3*ones(N-1,1);2])+ diag(-ones(N,1),-1)+diag(-ones(N,1),1));
D2=sparse(diag([3;4*ones(N-1,1);3])+ diag(-ones(N,1),-1)+diag(-ones(N,1),1));
D=sparse(blkdiag(D1,kron(eye(N-1),D2),D1)+kron(diag(ones(N,1),1)+diag(ones(N,1),-1),D0));


mue=0.9;
h=1/N;

%% u Berechnen

% Backslashoperator
% u=(mue/(h^2)*D+speye((N+1)*(N+1)))\v;

% Eigener Bandgauß Algorithmus

p=N+1;
q=N+1;
u=bandgauss(((mue/(h^2)*D+speye((N+1)*(N+1)))),p,q,v);

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

%% Visualisierung vom Zettel kopiert

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

%% Vergleiche Ergebnis - Code nicht lauffähig
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);
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 = (mue*D/h^2 + speye(N1N1)) \ v;
u = reshape(u,N1,N1);
v = reshape(v,N1,N1);
subplot(1,3,3);
imagesc(u); colormap('gray'); axis image; axis off;
title('entrauscht Musterlösung Zettel')

end

function [x] = bandgauss(KoeffizientenBandMatrix,AnzahlUntereNebendiagonalen,AnzahlObereNebendiagonalen,RechteGleichungsseite)

% Bandmatrix A  p...Anzahl unter Nebendiagonalen-1 q...Anzahl obere
% Nebendiagonalen-1

%Ax=b

A=KoeffizientenBandMatrix;
p=AnzahlUntereNebendiagonalen;
q=AnzahlObereNebendiagonalen;
b=RechteGleichungsseite;

%% Gauß für Bandmatrizen beginnen

n=size(A,1); %Zeilen von A
if (size(A,2) ~= n)
  error('A soll quadratisch sein.')
end

for k=1:1:n %%n Schritte 1,...,n
    
%   for i=k:1:min(n-1,k+p+1) % Zeilen 1 bis min (n-1 , k+p+1)
    for i=(k+1):1:min(n,k+p) % Zeilen 1 bis min (n-1 , k+p+1)
%       A(i+1,k)=A(i+1,k)/A(k,k);
        A(i,k)=A(i,k)/A(k,k);
        
%       for j=k+1:1:size(A,2)%% Spalten ab Diagonale bis min {n oder Spalte+q+1)
        for j=(k+1):1:min(n,k+p)%%

%                   A(i+1,j)=A(i+1,j)-A(k,j)*A(i+1,k);
                    A(i,j)=A(i,j)-A(i,k)*A(k,j);
        end
        
%       for j=k+1:1:size(b,2)
%           b(i+1,j)=b(i+1,j)-A(k,j)*b(i+1,k);
%       end
        b(i)=b(i)-A(i,k)*b(k);
        
    end
end

%% Rückwärtssubstitution
x=zeros(size(A,1),size(b,2));

x(n) = b(n)/A(n,n);
%for k=1:1:size(b,2)
    for i=n:-1:1
%	x(i,k) = (b(i,k)-A(i,i+1:n)*x(i+1:n,k)) / A(i,i);
	x(i) = (b(i)-A(i,i+1:n)*x(i+1:n)) / A(i,i);
    end
end
