function z3_8a()
%
% Kopf für jeden Übungszettel
%
ZETTEL = '3';
AUFGABE = '8a';

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');
%
%
%
    m = 19;
    a = -2; b = 2;
    t = zeros(m + 1, 1); % t und s als Vektor der Länge 20 setzen
    s = zeros(m + 1, 1);
    tau = zeros(m, 1); % Zwischenstellen "Tau"
    for i = 0:1:m
        % Werte von t und s nach Definition setzen
        t(i + 1) = a + ( ( (b - a) / m ) * i );
        s(i + 1) = f( t(i + 1) );
        if (i ~= m) % tau laueft nur von 0 bis 18 !!!
                    % tau(19) also nicht definiert
            tau(i + 1) = a + ( (b - a) / (2*m) ) + ( (b - a) / m) * i;
        end;
    end;
    
    % Koeffizienten fuer das Polynom bestimmen
    K = divDiff(t, s);


    % Polynom und f an den Zwischenstellen t auswerten
    % und gleichzeitig den "echten" Wert f(t) ausgeben
    for i = 1:1:(m+1)
        fprintf('p(t) = %1.16f fuer t = %1.16f (i = %d)\n', NewtonPolynom(t(i), K, t), t(i), i - 1);
        fprintf('f(t) = %1.16f\n', s(i));
    end;
    
    
    fprintf('\n');
    % Fehler bestimmen: | f(tau) - p(tau) |
    fehler = zeros(m, 1);
    for i = 1:1:length(fehler)
        fehler(i) = abs( f(tau(i)) - NewtonPolynom(tau(i), K, t) );
        fprintf('Fehler %1.16f bei tau = %1.16f\n', fehler(i), tau(i));
    end;
    
    
    
    % Ein paar Werte fuer das Zeichnen... Uninteressant
    u = t(1) : 0.001 : t(length(t));
    tmp = zeros(length(u), 1);
    for i = 1:1:length(tmp)
        tmp(i) = NewtonPolynom(u(i), K, t); % Auswertung des Polynoms speichern
    end;
    
    figure(1);
    subplot(2, 1, 1);
    % Polynom und Funktion f plotten mit den Stuetzstellen
    plot(t, s, 'ms', u, f(u), 'g-', u, tmp, 'b-');
    xlabel('x'); ylabel('y');
    title('f(t) und p(t)');
    legend('Stuetzstelle t_i', 'f(t)', 'p(t)');
    
    subplot(2, 1, 2);
    % Fehler plotten
    plot(tau, fehler, 'r-s');
    xlabel('x'); ylabel('y');
    title('Fehler zwischen f und p');
    
    
    % AB HIER:
    % Berechnung des m's, für welches der Fehler klein wird...
    fprintf('\n');
    m_tmp = 1; e_m_tmp = 42;
    for m = 1:1:99
        e(m) = maxFehler(m, a, b);
        if (e(m) < e_m_tmp) e_m_tmp = e(m); m_tmp = m; end;
        fprintf('m = %02d => e_m = %1.16f\n', m, e(m));
    end;
    fprintf('=> kleinster Fehler bei m = %d\n\n', m_tmp);
    
    % und das ganze plotten:
    u = [1:1:99];
    figure(2);
    semilogy(u, e);
    title('Fehler in Abh. von m, log. skaliert');
    xlabel('m'); ylabel('e_m');
    
    clear all;

function y = f(x)
    % f(t) aus Aufgabenteil a)
    y = exp( -(x.^2) );

function y = NewtonPolynom(x, K, t)
    % Bestimmt das Interpolationspolynom in der Newtonschen
    % Darstellung abhaengig von den Interpol.-Wertepaaren (t, s)
    % und wertet es an der Stelle x aus. Die Auswertung wird
    % durch das Horner-Schema vorgenommen.
    % Die Koeffizienten muessen durch K gegeben sein!
    
    % Mit dem Horner-Schema auswerten
    y = K(length(K)); % a(length(a)) ist der letzte Koeffizient
    for i = 1:1:length(t)
        j = length(t) - i + 1; % Index fuer die Rueckwaertssummierung
        y = y * (x - t(j)) + K(j);
    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);
    
 
 function e = maxFehler(m, a, b)
    % Weil's so schoen ist mal ohne Kommentare...
     
    t = zeros(m + 1, 1); s = zeros(m + 1, 1);
    for i = 0:1:m
        t(i + 1) = a + ( ( (b - a) / m ) * i );
        s(i + 1) = f( t(i + 1) );
    end;

    K = divDiff(t, s);
    e = 0;
    
    for i = 1:1:2000
        c = a + ( (b - a) / 2001 ) * i;
        f = abs( f(c) - NewtonPolynom(c, K, t) ); % begannener Fehler
        if (f >= e) e = f; end;
    end;
    
    
    
    
    
