function z6_16()
%
% Kopf für jeden Übungszettel
%
ZETTEL = '6';
AUFGABE = '16';

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');
%
%
%
    fprintf('Bitte etwas warten...\n');

    times = zeros(7, 1); times(1) = 0; 
    values = zeros(7, 1);

    % Berechnung aller Werte
    values(1) = abs(F1(-4) - F1(2)); % "echter Wert"

    % TEIL (A)
        TIME = clock; % Zeitmessung
    values(2) = integrateNumerical(@f1, -4, 2, 30, 10^(-2), 0);
        times(2) = etime(clock, TIME);
        TIME = clock;
    values(3) = integrateNumerical(@f1, -4, 2, 30, 10^(-7), 0);
        times(3) = etime(clock, TIME);
 
    % TEIL (B)
    % Die Werte 200 und 55000 fuer die Allokierung der Arrays habe ich
    % vorher grob ermittelt. Sie sind in jedem Fall groesser als die
    % eigentliche Groesse von p = 168 bzw. p = 50238, welche am Ende der
    % Berechnung vorliegt
        TIME = clock;
    values(4) = integrateNumerical(@f1, -4, 2, 30, 10^(-2), 200);
        times(4) = etime(clock, TIME);
        TIME = clock;
    values(5) = integrateNumerical(@f1, -4, 2, 30, 10^(-7), 55000);
        times(5) = etime(clock, TIME);
    
    % TEIL (C)
        TIME = clock;
    values(6) = recursivelyIntegrate(@f1, -4, 2, 30, 10^(-2), 6);
        times(6) = etime(clock, TIME);
        TIME = clock;
    values(7) = recursivelyIntegrate(@f1, -4, 2, 30, 10^(-7), 6);
        times(7) = etime(clock, TIME);

    % Ausgabe
    fprintf('Geschaetzter Wert des Integrals:\n    30\n\n');
    fprintf('Durch Stammfkt numerisch errechnet:\n    %2.16f\n\n', values(1));
    fprintf('Aufg. (a) mit accuracy = 10^(-2)\n    %2.16f\n', values(2));
    fprintf(' =>  %2.16f   Abweichung\n    Gebrauchte Zeit: %02.3fs\n\n', ...
        abs(values(1) - values(2)), times(2));
    fprintf('Aufg. (a) mit accuracy = 10^(-7)\n    %2.16f\n', values(3));
    fprintf(' =>  %2.16f   Abweichung\n    Gebrauchte Zeit: %02.3fs\n\n', ...
        abs(values(1) - values(3)), times(3));
    fprintf('Aufg. (b) mit accuracy = 10^(-2) und Speicherallokierung\n    %2.16f\n', values(4));
    fprintf(' =>  %2.16f   Abweichung\n    Gebrauchte Zeit: %02.3fs\n\n', ...
        abs(values(1) - values(4)), times(4));
    fprintf('Aufg. (b) mit accuracy = 10^(-7) und Speicherallokierung\n    %2.16f\n', values(5));
    fprintf(' =>  %2.16f   Abweichung\n    Gebrauchte Zeit: %02.3fs\n\n', ...
        abs(values(1) - values(5)), times(5));
    fprintf('Aufg. (c) mit accuracy = 10^(-2) rekursiv\n    %2.16f\n', values(6));
    fprintf(' =>  %2.16f   Abweichung\n    Gebrauchte Zeit: %02.3fs\n\n', ...
        abs(values(1) - values(6)), times(6));
    fprintf('Aufg. (c) mit accuracy = 10^(-7) rekursiv\n    %2.16f\n', values(7));
    fprintf(' =>  %2.16f   Abweichung\n    Gebrauchte Zeit: %02.3fs\n\n', ...
        abs(values(1) - values(7)), times(7));

end


function I = integrateNumerical(func, a, b, W, accuracy, allocate)
    % Diese Funktion integriert numerisch die Funktion func (mittels feval)
    % auf [a, b]. Der Algorithmus entstammt Seite 103.
    % Der Schaetzwert W muss gegeben sein, sowie die Genauigkeit accuracy =
    % epsilon.
    % allocate bezeichnet die Groesse des zu allokierenden Arrays. Ist dies
    % nicht bekannt, sollte man es einfach auf 0 setzen.
    % Moeglicher Aufruf: integrateNumerical(@myFunction, -1, 1, 0.5, 10^-8, 0);
    
    I = 0;
    j = 0; k = 1; p = 1; % Initialisierung
    if (allocate < 2) allocate = 2; end;
    t = zeros(allocate, 1); f = zeros(allocate, 1); i = zeros(allocate, 1);
    t(1) = a; t(2) = b; i(1) = 1;
    f(1) = feval(func, a);
    f(2) = feval(func, b);
    
    while 1
        
        % Berechnung
        h = t(k + 1) - t(j + 1);
        T = (h / 2) * (f(k + 1) + f(j + 1));
        tau = (t(k + 1) + t(j + 1)) / 2;
        y = feval(func, tau);
        S = (T + (2 * h * y)) / 3;
        
        r = abs((T - S) / W);
        
        if (r <= accuracy * (h / (b - a)))
            
            I = I + S;
            
            if (k == 1)
                break; % Exitpoint der Schleife!
            else
                j = k;
                k = i(k + 1);
            end;
            
        else
            
            p = p + 1;
            i(p + 1) = k;
            i(j + 1) = p;
            k = p;
            f(p + 1) = y;
            t(p + 1) = tau;
            
        end;
        
    end;
    

end


function I = recursivelyIntegrate(func, a, b, W, accuracy, diff)
    % Berechnet rekursiv das Integral fuer f auf dem Interval [a,b] wie in
    % Aufgabenteil (c) und damit speichersparend.
    % W ist der geschaetzte Wert des Integrals und accuracy die Genauigkeit
    % und diff ist die Differenz b - a (muss beim Aufruf gesetzt werden)
    % Bsp: recusivelyIntegrate(@myFunc, -5, 10, 1, 10^(-2), 10 + 5);
    I = 0;
    simpsonsum = calcSimpsonSum(func, a, b);
    e = abs( (calcTrapezeoidSum(func, a, b) - simpsonsum) / W );
    
    if(e <= accuracy * ((b - a) / diff))
        I = simpsonsum;
    else
        I = recursivelyIntegrate(func, a, (a + b) / 2, W, accuracy, diff) ...
            + recursivelyIntegrate(func, (a + b) / 2, b, W, accuracy, diff);
    end;
end


function y = calcSimpsonSum(f, a, b)
    % vgl. Seite 94
    y = ((b - a) / 6) * (feval(f, a) + feval(f, b) + 4 * feval(f, (a + b) / 2 ));
end

function y = calcTrapezeoidSum(f, a, b)
    % vgl. Seite 94
    y = ((b - a) / 2) * (feval(f, a) + feval(f, b));
end


function y = f1(x)
    % f(x) wie in der Aufgabe
    y = 5 - ((pi^3 * x) / (x^2 + 1)^2) * cos(25 / (x^2 + 1));
end

function y = F1(x)
    % Stammfunktion zu f1(x)
    % Die Stammfkt. habe ich mit Mathematica ermittelt, aber natuerlich mit
    % Ableiten ihre Richtigkeit bewiesen.
    y = 5*x + (1/50) * pi^3 * sin(25 / (1 + x^2));
end

