function Hrassnigg

clc;
clear;
clf;

%% Anfangsdaten
N=30;
M=N;
N1 = N+1;
N1N1=N1^2;
x=linspace(0,1,N+1);
y=x;
t=linspace(0,1,M+1);
h=1/N;
tau=1/M;

d = 1.2;
k = 1;
r = 3;

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

aa=linspace(4,4,(N+1)^2);
for i=1:N+1
    aa(i*(N+1))=3;
end

D=diag(aa)-diag(ones((N+1)^2-1,1),1)-diag(ones((N+1)^2-1,1),-1)...
     -diag(ones((N+1)^2-N-1,1),N+1) -diag(ones((N+1)^2-N-1,1),-N-1);
 
for i=linspace(N1,N1N1-N1,N1-1)
     D(i+1,i)=0;
     D(i,i+1)=0;
end
 
aaa=ones((N+1)^2,1);

F=diag(ones((N+1)^2,1))-diag(ones((N+1)^2-1,1),-1);

for i=linspace(N1,N1N1-N1,N1-1)
    F(i+1,i)=0;
end

I2=eye(N1N1);

B=d/h^2*D+k/h*F+r*I2;
A=(I2+tau*B);

%% Lösung des Problems  1. ohne Präkonditionierung!!!

V=V0;

for n=0:M

    if (k == 0)
        V=LoeseSystem(A,V);
    else
        V=LoeseSystem(A'*A,A'*V);
    end
 
    surf(reshape(V,N1,N1)); axis([1 N1 1 N1 0 1]);
    view([1 1 1]); drawnow;
    title('ohne Praekonditionierung') 
    pause(0.1);
 
end
 
% Anhand der graphisch Dargestellten Ergebnisse würde ich vermuten,
% dass wir einen Rand bei 0 festhalten und beim gegenüberliegenden 
% Rand die Schadstoffe entweichen. 

%% Lösung des Problems 2. mit Präkonditionierung
 
clear V;
 
C=diag(diag(A));
Cinv=zeros(N1N1);
 
for s=1:N1N1
    Cinv(s,s)=1/C(s,s);
end

A_praekond=Cinv*A*Cinv;
V0_praekond=Cinv*V0;
V_praekond=V0_praekond;
 
for ss=1:M

    V_praekond=LoeseSystem(A_praekond,V_praekond);
 
    V=Cinv*V_praekond;
    V_praekond=Cinv*V;

    surf(reshape(V,N1,N1));
    axis([1 N1 1 N1 0 1]);
    view([1 1 1]); drawnow;
    title('mit Praekonditionierung') 
    pause(0.1);

end 

% Für welche Transportparameter 
% ist die Konvergenz der Methode garantiert?

% Solange A SPD ist funktioniert die Methode
% des konjugierten Gradienten gut.(Matrix nur bei k=0 symmetrisch)

%%  Wenn A *nicht* SPD ist, welche Methoden haben wir?
% Jacobi und  Gauss-Seidel könnten funktionieren, 
% Matrix ist aber nicht diagonal dominant. 
% Gauß-Elimination/backslash funktioniert.

%% Lösung des Problems 3. mit backslash
clear V;

V_n=V0;
for i=0:M
    
    V_nplus1=A\V_n;
    V_n=V_nplus1;
    V=V_n;
    
    surf(reshape(V,N1,N1)); axis([1 N1 1 N1 0 1]);
    view([1 1 1]); drawnow;
    title('mit backslash') 
    pause(0.1);

end

end

function[V]=LoeseSystem(A,M)

    ny=0.01;
    x_kminus1=M;
    xk=M;
    v_kplus1=M-(A*M);
    rk=v_kplus1;

    while norm(M-A*x_kminus1,2)>(ny*norm(M,2))

        r_kminus1=rk;
        vk=v_kplus1;
        x_kminus1=xk;
        
        tk=r_kminus1'*r_kminus1/(vk'*A*vk);
        xk=x_kminus1+tk*vk;
        rk=r_kminus1-tk*A*vk;
        sk=rk'*rk/(r_kminus1'*r_kminus1);
        v_kplus1=rk+sk*vk;
        
%norm(M-A*x_kminus1,2)-(ny*norm(M,2))

    end

    V=xk;

end
