%% Übungsblatt 4 - Beispiel 6
clc; clear; clf;
%% Anfangsdaten
N = 50;
N1 = N+1;
h = 1/N;
N1N1 = N1*N1;
M = 20;
tau = 1/M;
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);
r = 2;
k = 0.0; % konvergiert bei großem k
% damit A SPD ist, müsste ja k = 0 sein, aber mit Methode 1 konvergiert es
% auch so.
d = 1; % sollte > 1 für CG sein

%% Matrix B
D0 = -ones(N1N1,1); 
D = spdiags(D0,-N1,N1N1,N1N1) + spdiags(D0,N1,N1N1,N1N1) + spdiags(D0,1,N1N1,N1N1) ...
    + spdiags(D0,-1,N1N1,N1N1) + speye(N1N1)*4; D(N1N1,N1N1) = 3;
F = speye(N1N1) + spdiags(D0,-1,N1N1,N1N1);
for i = 1:N1-1;
        D(i*N1+1,i*N1) = 0;
        F(i*N1+1,i*N1) = 0;
        D(i*N1,i*N1+1) = 0;
        D(i*N1,i*N1) = 3;
end

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

%% Grafische Darstellung
V = V0;

for n = 0:M;
% Methode 1 von Kelley (bei dieser Methode sind die Transportparameter egal,
% da diag(diag(A)) eine symmetrische, positiv definite Matrix ist
    imax = N1N1;

    eps = 1.0e-5; 
 
    r = V-A*V; 
    rho = norm(r,2)^2; i = 2;
    d = diag(A);
    while sqrt(rho)>eps*norm(V,2) && i<imax; % Abbruchskriterium
        z = r./d;
        if i == 2;
            t1 = z'*r;
            b = 0; p = z;
        else
            t0 = t1;
            t1 = z'*r;
            b = t1/t0; p = z + b*p;
        end
        w = A*p; 
        a = t1/(p'*w); 
        V = V + a*p; 
        r = r-a*w; 
        rho = r'*r;
        i = i+1;
    end

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

end

V = V0;

for n = 0:M;
% Methode 2: Backslash konvergiert sehr schnell und sogar für meinen
% Problemfall d = 0
    V = A\V;

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

end

V = V0;
[L,U] = lu(A);

for n = 0:M
% Methode 3: LU Zerlegung von A
     V = L\V;
     V = U\V;
    
% Diese Methode konvergiert, jedoch muss man weitaus öfters iterieren, z.B.
% 2*M mal, wenn vor allem k nahe an 0 ist.

    surf(reshape(V,N1,N1)); axis([1 N1 1 N1 0 1]);
    title('Methode der LU Zerlegung');
    view([1 1 1]); drawnow; pause(0.1);

end

V = V0;
L = tril(A,-1);
U = triu(A,1);
D = diag(diag(A));
imax = 100;
tol  = 1.0e-5;

for n = 0:M
% Methode 3: Symmetrisch Gauß-Seidel Methode
% (für streng diagonal dominante Matrizen)
    for i=1:imax
        Vi = V;
        V  = (D + L)\(V - U*V);
        V  = (D + U)\(V - L*V);
        if (norm(V-Vi) < tol*norm(V))
            break;
        end
    end
    
% Diese Methode konvergiert, jedoch muss man weitaus öfters iterieren, z.B.
% 2*M mal, wenn vor allem k nahe an 0 ist.

    surf(reshape(V,N1,N1)); axis([1 N1 1 N1 0 1]);
    title('Gauss-Seidel Methode');
    view([1 1 1]); drawnow; pause(0.1);

end
