function Blatt3

% Intervallgrenzen
a = 0.1;
b = 3;
% Funktion
f = @(x) 1./x;
% exakter Integralwert
exact_val = log(3) - log(0.1);
% Q und AnzF als globale Variablen:
% damit sind Q und AnzF in den rekursiven Funktionsaufrufen von QuadAdapt
% immer bekannt und können aktualisiert werden.
global Q
global AnzF
% Q unsere Approximation an das Integral wird mit 0 initialisiert, AnzF,
% die Anzahl an Funktionsauswertungen von f (was der Anzahl der am Ende
% verwendeten Knoten entspricht) wird mit 3 gestartet, da beim ersten
% Aufruf von QuadAdapt f 3 mal ausgewertet wird.
Q = 0;
AnzF = 3;
% Verkleinerungsschritte der Toleranz
steps = 5; 
% Initialisiere Vektoren für plot
err_Adapt = zeros(steps,1);
err_Trapez = zeros(steps,1);
AnzF_vec = zeros(steps,1);
for i = 1: steps
    tol = (1/10)^i; % Toleranz setzen
    % Rechne Approximation via Adaptive Quadratur - hier sind die ersten 3 Funktionsauswertungen
    [Approximation, Anzahl_Knoten] = QuadAdapt(f,tol,a,b,f(a),f((a+b)/2),f(b));
    % Rechne Approximation via Trapez
    n = Anzahl_Knoten - 1; % Teilintervalle der Trapezformel
    h = (b-a)/n; % Schrittweite der Trapezformel
    val = Trapez_Blatt1(f,h,a,b);
    % Zwischenspeichern
    err_Adapt(i) = abs(Approximation - exact_val);
    err_Trapez(i) = abs(val - exact_val);
    AnzF_vec(i) = Anzahl_Knoten;
    % Für den neuen Schleifendurchgang müssen natürlich Q und AnzF wieder auf die ursprünglichen Werte gesetzt werden
    Q = 0;
    AnzF = 3;
end

% plot mit labels und Legende
loglog(AnzF_vec,[err_Adapt,err_Trapez])
legend('error QuadAdapt','error Trapez')
xlabel('Anzahl Funktionsauswertungen')
ylabel('Integrationsfehler')
axis tight

end

function [Approximation, Anzahl_Knoten] = QuadAdapt(f,tol,a,b,fa,fab,fb)
% Realisiert die adaptive Quadraturformel

% Mit diesem Aufruf sind AnzF und Q innerhalb dieser Funktion bekannt
global AnzF
global Q

% Rechne T_2 und S_1 aus
T = Trapez(a,b,fa,fab,fb);
S = Simpson(a,b,fa,fab,fb);

% Prüfe nun, ob auf dem aktuellen Teilintervall halbiert werden muss oder nicht
if abs(T-S) > tol % falls Genauigkeit nicht erreicht
    tol = 1/2*tol; % halbiere Toleranz
    ab = (a+b)/2; % und das Interval
    % recycle alte Funktionswerte
    fa1 = fa; % fa linkes Intervall
    fb1 = fab; % fb linkes Intervall
    fa2 = fab; % fa rechtes Intervall
    fb2 = fb; % fa rechtes Intervall
    
    % rekursive Anwendung der Funktion auf dem linken Intervall
    fab = f((a+ab)/2); % neue Mitte des linken Intervalls
    AnzF = AnzF + 1; % hier findet 1 Funktionsauswertung (neuer Mittelpunkt) statt
    [~, ~] = QuadAdapt(f,tol,a,ab,fa1,fab,fb1);
    
    % rekursive Anwendung der Funktion auf dem rechten Intervall
    fab = f((ab+b)/2);
    AnzF = AnzF + 1; % Hier auch
    [~, ~] = QuadAdapt(f,tol,ab,b,fa2,fab,fb2);
    
    % Anmerkung: bei diesem Ansatz über Q und AnzF als globale Variablen
    % brauchen wir bei diesen rekursiven Subaufrufen von QuadAdapt keinen
    % output (mit "~" wird der output unterdrückt), denn wir interessieren
    % uns ja lediglich für die gesamte Anzahl Funktionsauswertungen und der
    % Approximation des gesamten Integrals. Diese werden aber mit AnzF und
    % Q "im Hintergrund" getrackt und können am Ende final ausgegeben
    % werden. 
else % falls Genauigkeit erreicht ist - nur dann wird auf einem Teilintervall etwas zum Gesamtintegral beigetragen
    Q = Q + S; % Addiere Wert der SimpsonFormel zur gesamten Integralapproximation
end

% Finale Ausgabe
Anzahl_Knoten = AnzF;
Approximation = Q;

end

function val = Trapez_Blatt1(f,h,startpoint,endpoint)
% Kopie der Trapezformel von der Lösung von Blatt1.

% Vektor mit Knoten
Knoten = startpoint:h:endpoint;
% f an den Knoten auswerten
eval_Knoten = f(Knoten);
% Approximation des Integralwertes bestimmen
val = h/2*(eval_Knoten(1)+eval_Knoten(end)) + h*sum(eval_Knoten(2:end-1));
end

function val = Trapez(a,b,fa,fab,fb)
% Vereinfachte Variante der zusammengesetzten Trapezregel T_2:
%
% Es ist T_2(f) = h/2 f(a) + h f((a+b)/2) + h/2 f(b) 
% und hier ist immer fa = f(a), fb = f(b) und fab = f((a+b)/2), sodass wir
% f auch nicht als Input benötigen, da wir schon alle Funktionsauswertungen
% kennen.
h = (b-a)/2;
val = h/2*fa + h*fab + h/2*fb;
end

function val = Simpson(a,b,fa,fab,fb)
% Realisierung der Simpsonformel S_1(f) = (b-a)/6 (f(a)+4f((a+b)/2)+f(b))
%
% Wie bei der Trapezformel ist fa = f(a), fb = f(b) und fab = f((a+b)/2),
% sodass f als Input nicht benötigt wird.
val = (b-a)/6*(fa+4*fab+fb);
end