function Knittelfelder

clc; clear; 
%% Beispiel 6 vom 3. Übungsblatt

N = 128;
N1 = N+1;
N4 = round(N/4);
V0 = [zeros(N4,1);ones(N1-2*N4,1);zeros(N4,1)];
h = 1/N;
r = 0.005; % r,k,d >= 0
k = 10; d = 0.100;

%% Matrix B und A
diag1 = -ones(N1,1);
M1 = speye(N1) + spdiags(diag1,-1,N1,N1);
M2 = speye(N1)*2 + spdiags(diag1,-1,N1,N1) + spdiags(diag1,1,N1,N1);
M2(N1,N1) = 1;
B = d/h^2*M2 + k/h*M1 + r*speye(N1);
tau = h; % tau > 0
A = speye(N1) + tau*B;
V = V0;

%% LU Zerlegung durch Gauß-Algorithmus
% Bedingung für Parameter: Sie müssen nur größer Null sein, da die Matrix
% B streng diagonal dominant ist
L = eye(N1);
A = [A, V0];
A1 = A;

for k = 1:N1-1;
    for i = k+1:N1;
        L(i,k) = A(i,k)/A(k,k);
        for j = k+1:N1+1;
            A(i,j) = A(i,j) - L(i,k)*A(k,j);
        end
    end
end

U = A(:,1:N1);
U = triu(U);
norm(A1(:,1:N1)-L*U)

%% Gauß-Seidel

x = linspace(0,1,N1);
nmax = 25;
for n = 1:nmax;
    V = vorwaerts(L,V,N1);
    V = rueckwaerts(U,V,N1);
    plot(x,V); axis([0 1 0 1]);
    drawnow; pause(0.1);
end

% Reaktion: Einzig bei der Reaktion bleibt die treppenfunktion erhalten
% Diffusion: Die Treppenfunktion wird geglättet und fügt sich den
% Randbedingungen
% Konvektion: Entspricht am Ende den Randbedingungen, wobei der Graph bis
% zu x = 0.65 bei 0 bleibt, die Konvektion muss größer 1 sein, damit der
% Graph zu Null geht.

end

function [ B ] = vorwaerts( L,V,N1 )
L = [L,V];
B = zeros(N1,1);
B(1) = L(1,N1+1)/L(1,1);
for i = 2:N1;
    s = 0;
    for j = 1:i-1;
        s = s + L(i,j)*B(j);
    end
    B(i) = (L(i,N1+1)-s)/L(i,i);
end

end

function [ B ] =rueckwaerts( U,V,N1 )
U = [U,V];
B = zeros(N1,1);
B(N1) = U(N1,N1+1)/U(N1,N1);
for i = N1-1:-1:1;
    s = 0;
    for j = i+1:N1;
        s = s + U(i,j)*B(j);
    end
    B(i) = (U(i,N1+1)-s)/U(i,i);
end

end

