function z5_14()
%
% Kopf für jeden Übungszettel
%
ZETTEL = '5';
AUFGABE = '14';

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');
%
%
%
    
    for m = 1:1:6
        diffFormulaAndError(m, 1/2);
    end;
    
 
function diffFormulaAndError(m, point)
    % Macht alles, was in der Aufgabenstellung gefordert ist (bis auf
    % das Kommentieren) in Abh. von m und der Stelle 'point', an der
    % die Abl. approximiert werden soll.
    D_ex = Df(point, m);
    fprintf('\t\t\t ***** m = %d *****\n', m);
    fprintf('f^(%d) (%1.2f)          = %+2.16f\n', m, point, D_ex);
    for i = 0:1:8
        h(i + 1) = 10^(-8 + i); % vgl. Zettel
        % approximierte Ableitung berechnen durch Interpolation
        D_intp(i + 1) = calcDiffFormula(point, h(i + 1), m); 
        fprintf('h = 10^%+d => l_%d,%d(f) = %+2.16f\n', -8 + i, m, m, D_intp(i + 1));
    end;
    
    figure(1); subplot(3, 2, m); % Alles in eine Zeichnung plotten
    loglog(h, abs(D_intp - D_ex), '-rs', 'MarkerFaceColor', 'g');
    xlabel('Wert von h'); ylabel('Fehler im Betrag'); title(sprintf('m = %d', m));
    fprintf('\n');


function y = calcDiffFormula(point, h, m)
    % Aufgabenteil (a)
    % Berechnet die approximierte Ableitung der Funktion f an der Stelle h
    % mit Intervalllaenge 2h und m aequidistanten Stuetzstellen
    t = zeros(m + 1, 1); s = zeros(m + 1, 1);
    for i = 0:1:m
        t(i + 1) = point - h + 2*h*(i/m); % Aequidistante Aufteilung
        % des Intervalls [point - h, point + h]
        s(i + 1) = f(t(i + 1)); % Interpolationswerte
    end;
    
    a = divDiff(t, s); % Dividierte Differenz
    y = factorial(m) * a(length(a)); % = m! * a_m = l_{m,m} (f)


function y = f(x)
    y = sin(x);

    
function y = Df(x, k)
    % Gibt den Wert der k-ten Ableitung an der Stelle x zurueck
    if (k <= 0) k = 0; end;
    k = mod(k, 4); % Da sin beim Ableiten "zyklisch" ist, entspricht
    % die 5. Ableitung der 1. usw. Somit sind nur 4 Faelle zu betrachten:

    if(k == 0) % Keine Ableitung. Also f(x) = sin(x)
        y = f(x); return;
    elseif(k == 1) % 1. Ableitung ist f'(x) = cos(x)
        y = cos(x); return;
    elseif(k == 2) % f''(x) = -sin(x) = -f(x)
        y = (-1) * f(x); return;
    elseif(k == 3) % f'''(x) = -cos(x) = -f'(x)
        y = (-1) * Df(x, 1); return;
    end;


function [values] = divDiff(t, s)
    % diffDiv(t, s) berechnet die Koeffizienten
    % fuer die Tupel (t, s) und gibt diese als
    % Array zurueck
    % Berechnungsschema: Dividierte Differenzen (Newton)
    if (length(t) ~= length(s))
        % Ups, t und s sind nicht gleichlang, Fehler!
        disp('Fehlerhafte Benutzung');
        return;
    elseif (length(s) == 1)
        % Falls nur ein Wertepaar gegeben ist, ist
        % dies bereits die Loesung
        values = s;
        return;
    end;
    
    % Nullmatrix der Groesse n x n anlegen
    A = zeros(length(s), length(s));
    % erste Spalte = s setzen
    A(:, 1) = s;
    
    % nun wie im Skript...
    for j = 1:1:(length(s) - 1)
        for i = j:1:(length(s) - 1)
            A(i + 1, j + 1) = A(i + 1, j) - A(i, j);
            A(i + 1, j + 1) = A(i + 1, j + 1) / (t(i + 1) - t(i - j + 1));
        end;
    end;
    
    % Rueckgabe der Funktion bestimmen: Die Diagonale
    % der berechneten Matrix A als Vektor
    values = diag(A);
    
    
    
