%%Einführung in die numerische Mathematik
%%Beispiel 5

%Gernot Holler
%%----------------------------------------------------------

clear;clc;clf;

N=50;
my=0.001;

%%Gegebene Anfangsbedingungen
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);


%%Erzeugen der Matrix 
Z=eye((N+1)*(N-1));
T=-eye(N1^2-1);
D02=-eye(N1^2);
a=zeros(N1^2-1,1);
b=zeros(1,N1^2);
T1=[b; T a];
T2=[a T;b];
D2=T1+T2-4*D02;
D2(N1^2,N1^2)=2;
D2(1,1)=2;

m=N1;
for k=1:N  
    D2(m,m)=3;
    D2(1+m,1+m)=3;
    D2(m,1+m)=0;
    D2(1+m,m)=0;
    m=m+N1;
end
j=2;
for k=2:N
    D2(j,j)=3;
    D2(N1^2-j+1,N1^2-j+1)=3;
j=j+1;
end
D2(N1,N1)=2;
D2(N1^2-N,N1^2-N)=2;

G=eye(N1*N);
H=zeros(N1,N1^2);
O=zeros(N1,N1*N);
B=-[ O' G; H];
V=-[H; G O'];

h=1/N;

D=B+V+D2;
A=(my/h^2)*D-D02;

%%Gaußsche Elimination


%Erweitern der Matrix

A=[A v ];

tic

p=1;


while p<N1N1
    
    i=p+1;
    
    
    while i<=N1N1 && (i-p)<=N1 %% Ausnutzen der Matrixstruktur
        
        mp=A(i,p)/A(p,p);
        
          A(i,N1N1+1)=A(i,N1N1+1)-A(p,N1N1+1)*mp;
        
       j=p;
        
        while j<=(N1N1) && (j-p)<=N1 %% Ausnutzen der Matrixstruktur
            
             
          
            
            A(i,j)=A(i,j)-A(p,j)*mp;
           
            
            j=j+1;
            
        end
        
         
        i=i+1;
        
    end
    
    p=p+1;
    
end


%% Rücksubstitution
 
  u(N1N1)=A(N1N1,N1N1+1)/A(N1N1,N1N1);

  
      for j=N1N1-1:-1:1
          
         s=0;
        
          for k=N1N1:-1:j+1
              
              s=s+A(j,k)*u(k);
              
          end
          u(j)=(A(j,N1N1+1)-s)/A(j,j); 
          
      end

toc
     
 u = reshape(u,N1,N1);
 v = reshape(v,N1,N1);
 
 
 
 %%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');   
%     
% 
%Gesamtkosten der Gaußschen Elimination ohne Ausnutzung der Matrixstruktur:
%
%Divisionen:(N1N1*(N1N1-1))/2
%Multiplikationen: N1N1*(N1N1-1)/3
%Subtraktionen: N1N1*(N1N1-1)/3
%   
%     
% 
%     
%     
%     
%     
%     
%     
%     
% 


























