clc;clear;clf;

% Dieser Code ist nur für kleine N sinnvoll. 
% Für N > 25 braucht er zu lange. (wegen der Erzeugung der Matrix D)
% BSP: N=60, mu=0.1 -> sec: 1962.1482

N=input('Größe der Matrix (N>8) =' );
N1 = N+1;
N4 = round(N/4);
N1N1 = N1*N1;

ustar = [zeros(N4,1);ones(N1-2*N4,1);zeros(N4,1)];
ustar = kron(ustar,ustar');

v = ustar + 0.1*randn(N1,N1);

v = reshape(v,N1N1,1);

%% Matrizen erzeugen


mu=input('Gib mu ein (mu>0) = ');

t=tic;
D0=-eye(N+1);

D1=zeros(N+1);
D2=zeros(N+1);



i=1;
while i<=N
    D1(i+1,i)=-1;
    D1(i,i+1)=-1;
    D1(i,i)=3;
    
    D2(i+1,i)=-1;
    D2(i,i+1)=-1;
    D2(i,i)=4;
    i=i+1;
end

    D1(1,1)=2;
    D1(N+1,N+1)=2;
    
    D2(1,1)=3;
    D2(N+1,N+1)=3;
    
    D3=zeros(N+1);
   
    % Matrizen zusammenstückeln (code ist 
    % sehr schwerfällig/ langsam)
    
    D=zeros(1,N1N1);
    
ii=1; %Block-Zeile
while ii<=N+1
    
    if  (ii==1)  %Block 1,1
        DZ=D1;

        else if (ii==2) %Block 2,1
        DZ=D0;
        
    else 
        DZ=D3;
            end
    end
       
 
    
    j=2;    %Spalte
    while j<=N+1
        
            if j==ii+1  %Nebendiagonale
                 DZ=[DZ D0];
                 
            else if (j==ii-1 && ii>2) %Nebendiagonale
                 DZ=[DZ D0];
                 
            else if (j==ii)&&(j~=N+1)
                 DZ=[DZ D2];
            else if (j==ii)&& (ii==N+1)
                    DZ=[DZ D1];
                 
             else              %alles sonst
                  DZ=[DZ D3];
                end
                end
                end
            end
            
    j=j+1;        
    end
    
    D=[D;DZ];
    
    
    ii=ii+1;
end
                    
                  
    D(1,:)=[];
    I=eye(N1N1);
    h=1/N;
    A=(mu/h^2*D+I);

   
 %% Gaußsche ELimination
 vneu=v;  
 
    jj=1;  %zeilen von A
    
    while jj<N1N1
    
        m=1/A(jj,jj);
        ij=jj+1;
        
        while ij<=N1N1
        vneu(ij)=vneu(ij)-(m*A(ij,jj)*vneu(jj));
        A(ij,:)=A(ij,:)-(m*A(ij,jj)*A(jj,:));
        ij=ij+1;
        end
        
    
    jj=jj+1;
    end
    
    % Gaußsche Elimination kostet: 
    % N1N1(N1N1 - 1)/2 + 2N1N1(N1N1^2 - 1)/3

    % Rückwärssubstitution:
    
    u=zeros(N1N1,1);
    
  

    for t1=N1N1-1:-1:1  %Zeile
       sum=0;
       u(N1N1)=vneu(N1N1)/A(N1N1,N1N1);
       
        for t2=N1N1:-1:t1+1
            sum=sum+A(t1,t2)*u(t2);
        end
 
           
        u(t1)=(vneu(t1)-sum)/A(t1,t1);
        
    end
    
       % Rücksubstitution kostet: (N1N1^2+N1N1-1)
        
    t=toc(t);
    
    
    %Bildliche Darstellung
vneu=v;    
u = reshape(u,N1,N1);
v = reshape(v,N1,N1);


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


%% Vergleichscode
tt=tic;

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 = (mu*D/h^2 + speye(N1N1)) \ vneu;

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

    tt=toc(tt);


disp(['Wir brauchen für die erste Berechnung (in sec) : ', num2str(t)]);

disp(['Wir brauchen für die zweite Berechnung (in sec): ', num2str(tt)]);

 
   
