function z9_23()
%
% Kopf für jeden Übungszettel
%
ZETTEL = '9';
AUFGABE = '23';

fprintf('\n\n###\tNumerik I SS 2010\t\t###\n###\tZettel %s, Aufgabe %s', ZETTEL, AUFGABE);
fprintf('\t\t###\n###\t\t\t\t\t###\n###\tvon\tFelix Nyenhuis\t\t###\n');
fprintf('###\t\tLouisa Marie Puls\t###\n###\t\tReidar Janssen\t\t###\n\n');
%
%
%

    Aufg23a(0, 10, 1/2, 0, 1/10);
    Aufg23a(0, 10, 1/2, 0, 1);
    Aufg23a(0, 10, 1/2, 1/2, 1/10);
    Aufg23a(0, 10, 1/2, 1/2, 1);
    Aufg23a(0, 10, 1/200, 0, 1/10);
    Aufg23a(0, 10, 1/200, 0, 1);
    Aufg23a(0, 10, 1/200, 1/2, 1/10);
    Aufg23a(0, 10, 1/200, 1/2, 1);

end


function Aufg23a(a, b, h, alpha, beta)

    function z = K(x1, x2)
        z = max(0, 1 - abs( beta * (x1 - x2) ));
    end

    function z = r(x1)
        z = 1 + cos((5.*x1) ./ pi);
    end

    N = abs(b - a) / h;
    x = a + h .* [0:1:N];
    
    A = zeros(N + 1, N + 1);
    for i = 0:1:N
       for j = 0:1:N
           if(i == j) A(i + 1, j + 1) = alpha; end;
           c = 1;
           if(j == 0 || j == N) c = 1/2; end;
           
           A(i + 1, j + 1) = A(i + 1, j + 1) + h * c * K(x(i + 1), x(j + 1));
       end
    end
    
    y = LoeseLGS_mitPivot(A, r(x), 'abscol');
    
    f = mod(round(1/h + 2 + alpha + 10 * beta), 42);
    
    fprintf('### Bild %d ### \n', f);
    fprintf('h = %1.3f, alpha = %1.1f, beta = %1.1f\n', h, alpha, beta);
    fprintf('=>\tcond(A) = %f\n', cond(A));
    fprintf('\tgroesstes Pivotel. = %f\n', y(N + 2));
    fprintf('\tkleinstes Pivotel. = %f\n', y(N + 3));
    fprintf('\n');
    
    figure(f);
    subplot(2, 1, 1);
    plot(x, y(1:1:N+1), 'k-');    
    title(sprintf('Bild %d\nh = %1.3f, alpha = %1.1f, beta = %1.1f\n', ...
        f, h, alpha, beta));
    xlabel('x'); ylabel('y(x)');
    
    subplot(2, 1, 2);
    spy(A);

end

function [b] = LoeseLGS_mitPivot(A, b, pivot)
    % Loest das lin. Gleichungssystem Ax = b
    % mit quadratischer Matrix A
    % mithilfe der LR-Zerlegung und der Pivotiersierung
    % pivot, wobei pivot folgendes sein kann:
    %   'none' (keine Pivotisierung), 
    %   'abscol' (absolute Spaltenpivot.), 
    %   'relcol' (relative Spaltenpivot., Standard)
    % ist b eine Matrix, so werden die Spalten als einzelne
    % Vektoren zur Ausrechnung verwendet
    
    [n, m] = size(A);
    [n_b, m_b] = size(b);
    if(m ~= n || n == 0)
        disp('Bitte nur quadratische Matrizen!');
        return;
    end
    if(n_b == 1) b = b'; end; % damit wir A*b machen koennen
    if(m_b ~= m)
        disp('Vektoren b muessen von der Laenge zu A passen!');
        return;
    end;
    
    % Pivotstategie pruefen
    if(nargin < 3) pivot = 'relcol'; end;
    
    [LR, P] = LR_P_Zerlegung(A, pivot); % LR, P-Zerlegung berechnen
    %LR, P
    
    % biggest or smalles pivot element
    bigP = -1; smallP = -1;
    for j = 1:1:length(A)
        if(abs(A(P(j), j)) > bigP || bigP < 0) bigP = A(P(j), j); end;
        if(abs(A(P(j), j)) < smallP || smallP < 0) smallP = A(P(j), j); end;
    end;
    
    % Berechnung mit b
    for j = 1 : n-1
        b = swapRows(b, P(j), j); % bei den b's die Zeilen tauschen
        % Quasi also Pb berechnen fuer jedes einzelne b
    end;
    
    for k = 1 : 1 : n_b % Alle b's durchgehen und die Loesung berechnen
        
        d = b(:, k);
        % Lu = Pb
        for j = 2 : n
            d(j) = d(j) - LR( j, 1 : j-1 ) * d( 1 : j-1 );
        end;
        
        % Rx = u
        for i = n : -1 : 1
            d(i) = d(i) - LR( i, i+1 : n ) * d( i+1 : n );
            d(i) = d(i) / LR(i, i);
        end;
        
        % Und das Ergebnis zurueckschreiben
        b(:, k) = d;
        
    end;
    
    b(size(b) + 1) = bigP;
    b(size(b) + 1) = smallP;

end

function [A, P] = LR_P_Zerlegung(A, pivot)
    n = length(A);
    P = zeros(length(A), 1);
    
    for j = 1 : 1 : n-1
        
        if(strcmp(pivot, 'abscol')) 
            % Absolute Spaltenpivotisierung
            P(j) = getAbsColPivot(A, j);
        elseif(strcmp(pivot, 'none'))
            P(j) = j; % Keine Pivotisierung
        else
            % Relative Spaltenpivotisierung (Standard)
            P(j) = getRelColPivot(A, j);
        end
        
        % Zeilen vertauschen
        A = swapRows(A, P(j), j);
        
        % Nun die eigentliche Berechnung
        for i = j+1 : n
            A(i, j) = A(i, j) / A(j, j); % q_i,k setzen
            A(i, j+1 : n) = A(i, j+1 : n) - ...
                A(i, j) * A(j, j+1 : n); % Zeilenumformung
        end
        
    end
    
    P(n) = n; % noch den letzten Eintrag setzen
end

function A = swapRows(A, k, l)
    % Vertauscht bei der Matrix A die Zeilen k und l,
    % falls l != k
    if(k ~= l)
        tmp = A(k, :); A(k, :) = A(l, :); A(l, :) = tmp;
    end
end

function p = getRelColPivot(A, j)
    % Geht sicher etwas besser und schneller, aber zur Demonstration 
    % einmal exakt wie die Definition
    p = 0; tmp_p = -2; tmp = -1;
    for i = j : length(A)
        tmp = abs(A(i, j)) / sum(abs(A(i, j : length(A))));
        if(tmp > tmp_p) 
            p = i; tmp_p = tmp;
        end;
    end;
end

function p = getAbsColPivot(A, j)
    % Geht sicher etwas besser und schneller, aber zur Demonstration 
    % einmal exakt wie die Definition
    p = 0; tmp_p = -2; tmp = -1;
    for i = j : length(A)
        tmp = abs(A(i, j));
        if(tmp > tmp_p)
            p = i; tmp_p = tmp;
        end;
    end;
end
