function z8_21_22()
%
% Kopf für jeden Übungszettel
%
ZETTEL = '8';
AUFGABE = '21,22';

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

    % Aufg. 21
    A = [1 0 2; 1 10^(-20) 3; 2 1 1];
    b = [15 17 21];
    x = zeros(3, 4);
    x(:,1) = LoeseLGS(A, b, 'none')';
    x(:,2) = LoeseLGS(A, b, 'abscol')';
    x(:,3) = LoeseLGS(A, b)';
    x(:,4) = (A \ b')';
    
    fprintf('Aufgabe 21: A und b wie auf Zettel. Berechne Ax = b ...\n');
    fprintf('... ohne Pivotisierung => x = (%1.3f, %1.3f, %1.3f)\n', x(1, 1), x(2, 1), x(3, 1));
    fprintf('... mit abs. Spalt.-P. => x = (%1.3f, %1.3f, %1.3f)\n', x(1, 2), x(2, 2), x(3, 2));
    fprintf('... mit rel. Spalt.-P. => x = (%1.3f, %1.3f, %1.3f)\n', x(1, 3), x(2, 3), x(3, 3));
    fprintf('... via A \\ b          => x = (%1.3f, %1.3f, %1.3f)\n', x(1, 4), x(2, 4), x(3, 4));
    fprintf('\nNun Aufgabe 22 (nur Plots)\n');    
    
    e = zeros(41, 1);
    for m = 2:2:80
        if(m == 26 || m == 50) fig = m;
        else fig = 0; end;
        e(m) = Aufg22(m, @funcA, fig);
    end;
    
    figure(42);
    m = 2:2:80;
    semilogy(m, e(m), 'k-');
    title('Maximaler Fehler');
    xlabel('m'); ylabel('e_m = max |f(t) - p(t)|');
    
end

function emax = Aufg22(m, f, fig)
    t = zeros(m + 1, 1);
    for j = 0:1:m
        t(j + 1) = (2*pi*j) / m;
    end;
    s = feval(f, t);
    
    K = TrigInterpKoeff(t, s');
    x = 0 : 0.001 : 2*pi;
    x = sort([x t']); % t dazu packen
    y_f = feval(f, x);
    y_p = TrigInterpolationVector(K, x);
    e = abs(y_p' - y_f);
    emax = max(e);
    
    if(fig > 0)
        figure(42*fig);
        subplot(2, 1, 1);
        plot(x, y_f, 'k-', t, s, 'ks'); 
        legend('f(x)', 'Stuetzstln.'); xlabel('x'); ylabel('f(x)');
        title(sprintf(' m = %d ', m));
    
        subplot(2, 1, 2);
        plot(x, y_p, 'k-', t, s, 'ks'); 
        legend('p(x)', 'Stuetzstln.'); xlabel('x'); ylabel('p(x)');
    
        figure(42*fig + 1);
        subplot(2, 1, 1);
        plot(x, y_f, 'r-', x, y_p, 'g-', t, s, 'ks');
        legend('f(x)', 'p(x)', 'Stuetzstln.'); xlabel('x'); ylabel('y');
        title(sprintf(' m = %d ', m));
    
        subplot(2, 1, 2);
        semilogy(x, e, 'r-');
        xlabel('x'); ylabel('Fehler im Betrag');
    end;
end

function x = TrigInterpolation(K, t)
    % Berechnet das trigonometrische Interpolationspolynom
    % zu den Koeffizieten K an der Stelle t
    if(mod(length(K), 2) == 0)
        disp('Falsche Benutzung');
        return
    end;

    c = (2*pi) / (2*pi + 1);
    for i = 1:1:((length(K) - 1) / 2)
        K(2*i) = K(2*i) * sin( (c*i) * t);
        K(2*i + 1) = K(2*i + 1) * cos( (c*i) * t);
    end;
    x = sum(K);
end

function [x] = TrigInterpolationVector(K, t)
    % Wie oben, nur diesmal fuer Vektoren t
    % (multiple Auswertung)
    x = ones(length(t), 1); x = x .* NaN;
    for i = 1:1:length(t)
        x(i) = TrigInterpolation(K, t(i));
    end;
end

function [x] = TrigInterpKoeff(t, s)
    % Bestimmt die Koeffizienten fuer die trigonometrische Interpolation
    m = length(t);
    if(m ~= length(s) || mod(m, 2) == 0)
        disp('Falsche Benutzung');
        return
    end;
    
    c = (2*pi) / (2*pi + 1);
    
    A = zeros(m, m);
    A(:, 1) = 1;
    for j = 1:1:((m - 1)/2)
        A(:, 2*j + 1) = cos((c .* j) .* t);
        A(:, 2*j) = sin((c .* j) .* t);
    end;
    
    x = LoeseLGS(A, s, 'abscol');
end

function [b] = LoeseLGS(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
    
    % 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;

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 y = funcA(x)
    y = (8/5) .* ((x - pi) .^ 2) .* sin(16 .* (x - pi));
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
