function Wegger

%% Numerik Blatt 03 Bsp 6

clear; clc; clf;
%% Initialisiere beliebige Daten
N = 100;
N1 = N+1;
N4 = round(N/4);
V0 = [zeros(N4,1);ones(N1-2*N4,1);zeros(N4,1)];
r = 1;          % r kleiner -> groeßere Schwingung
k = 5;          % Konvergenzgeschwindigkeit (schneller fuer großes k)
d = 1;        
h = 1/N;
tau = 0.01;

%% Erzeuge Matrix B
d1 = -ones(N1,1);

B1 = speye(N1)*2 + spdiags(d1,-1,N1,N1) + spdiags(d1,1,N1,N1);
B1(N1,N1) = 1;
B2 = speye(N1) + spdiags(d1,-1,N1,N1);

B = (d/h^2)*B1 + k/h*B2 + r*speye(N1);
A = speye(N1) + tau*B;
V = V0;

%% LU Zerlegung von A

L = eye(N1);
A = [A, V0];

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

%% Graphische Darstellung
x = linspace(0,1,N1);
nmax = 50;
for n = 1:nmax;
    V = vorwaerts(L,V);
    V = rueckwaerts(U,V);
    plot(x,V);
    axis equal;
    drawnow; pause(0.1);
end

end

%% Funktion fuer Vorwaertssubstitution
function [v] = vorwaerts(L,V)

L = [L,V];
b = size(V,1);
v = zeros(b,1);
v(1) = L(1,b+1)/L(1,1);

for i = 2:b
    s = 0;
    for j = 1:i-1
        s = s + L(i,j)*v(j);
    end
    v(i) = (L(i,b+1)-s)/L(i,i);
end

end

%% Funktion fuer Ruekwaertssubstitution

function [v] = rueckwaerts(U,V)
U = [U,V];
a = size(V,1);
v = zeros(a,1);
v(a) = U(a,a+1)/U(a,a);

for i = a-1:-1:1
    s = 0;
    for j = i+1:a
        s = s + U(i,j)*v(j);
    end
    v(i) = (U(i,a+1)-s)/U(i,i);
end

end
