%%Einführung in die numerische Mathematik 


%%Blatt 4
%%Gernot Holler


%%Beispiel 6

clear; clc;clf;
M=20;
N=30;
N1=N+1;
N1N1=N1*N1;

h=1/N;
tau=1/M;

%%k,r,d müssen größer als null sein. Das Verfahren der konjugierten
%%Gradienten funktioniert garantiert bei SPD-Matritzen

k=1;
r=1;
d=1;

%% Die Schadstoffe bewegen sich aufgrund der Anfangsdaten nur in eine
%% Richtung und die Bewegung ist am Rand begrenzt


%% Erstellen der Matrix A
N4=round(N/4);
V0=[zeros(N4,1);ones(N1-2*N4,1);zeros(N4,1)];
V0=kron(V0,V0');
V0=reshape(V0,N1N1,1);

D=diag([4*ones(N1N1-1,1);3])-diag(ones(N1N1-1,1),1)-diag(ones(N1N1-1,1),-1)-diag(ones(N1N1-N1,1),N1)-diag(ones(N1N1-N1,1),-N1);
F=eye(N1N1)-diag(ones(N1N1-1,1),-1);



for kg=N1:N1:N1N1-N1
    
    D(kg:kg+1,kg)=[3,0];
    D(kg,kg+1)=0;
    F(kg+1,kg)=0;
    
end

B=(d/(h^2))*D+(k/h)*F+r*eye(N1N1);

A=eye(N1N1)+tau*B;

%%Implementierungen zur Präkonditionierung

%K=diag(diag(A))^(1/2); %%Langsamer
%kinv=eye(N1N1)/K;
u=sqrt(diag(A));
K=diag(u); 
kinv=diag(1./u);
A=kinv*A*kinv;
       

    for j=0:M
       
xm=zeros(N1N1,1);

    x=K*V0;

v=kinv*V0-A*x;
r1=v;
i=1;
  
   
   while i<=N*N  && norm(x-xm)>1e-13*norm(x) %%mit geeignetem my   Bei exakter Arithmetik für die Methode der konjug iierten Gradienten für 1<=k<=N1N1 zum Ziel
    xm=x;
       
    t=(r1'*r1)/(v'*A*v);
    x=x+t*v;
    r2=r1-t*A*v;
    s=(r2'*r2)/(r1'*r1);
    v=r2+s*v;
    r1=r2;

    i=i+1;
   
     end

     V0=kinv*x;

    %% x=reshape(x,N1,N1);
  
    surf(reshape(V0,N1,N1));
    
    axis([1 N1 1 N1 0 1]);
    view ([1 1 1]); drawnow;
    pause(0.1);
    
    
    end
