function Blatt5
close all

% die verschiedenen Funktionen der Aufgabe
f1 = @(x) sin(7/4*pi*x + 1/4*pi);
f2 = @(x) sin(x);
% Start- und Endwert der Integration
a = -1; b = 1;
% Die exakten Werte der Integrale
exact_val_f1 = -4/(7*pi);
exact_val_f2 = 0;

%% Aufgabenteil (a)
a_err = zeros(3,1);
for i = 1:3
    a_err(i) = abs(Golub_Welsh(f1,5*i) - exact_val_f1);
end
sprintf('Für %i Knoten ergibt sich als Fehler %0.16f.\n Für %i Knoten ergibt sich als Fehler %0.16f.\n Für %i Knoten ergibt sich als Fehler %0.16f.\n', 5, a_err(1), 10, a_err(2), 15, a_err(3))

%% Aufgabenteil (b)
b_Knoten = 100:100:10000;
b_err = zeros(length(b_Knoten),1);
for i = 1:length(b_Knoten)
    b_err(i) = abs(Monte_Carlo(f2,b_Knoten(i),a,b) - exact_val_f2);
end
% plot des Monte-Carlo Fehlers inkl. der Monte-Carlo-Konvergenzgeschwindigkeit
figure;
loglog(b_Knoten,b_err, 'r', b_Knoten,1./sqrt(b_Knoten),'g')
legend('error Monte Carlo','1/sqrt(n)')
xlabel('Anzahl Stützstellen')
ylabel('absoluter Integrationsfehler')

%% Aufgabenteil (c)
c_Knoten = 4:2:1000;
c_vals = zeros(length(c_Knoten),3); % haben 3 verschiedene Verfahren
c_err = zeros(length(c_Knoten),3);
for i = 1:length(c_Knoten)
    m = c_Knoten(i);
    c_vals(i,1) = Golub_Welsh(f1,m);
    c_vals(i,2) = Monte_Carlo(f1,m,a,b);
    c_vals(i,3) = Trapez_Blatt1(f1,(b-a)/(m-1),a,b);
    c_err(i,:) = abs(c_vals(i,:) - exact_val_f1);
end
% plot
figure;
loglog(c_Knoten,[c_err(:,1),c_err(:,2),c_err(:,3)])
axis tight
legend('error_Golub_welsh','error_Monte_Carlo','error_Trapez')
xlabel('Anzahl Stützstellen')
ylabel('absoluter Integrationsfehler')

end

%% Die verschiedenen Verfahren

function val = Golub_Welsh(f,m)
% Realisiert die Approximation von \Int_-1^1 f(x)dx mit Hilfe des Golub Welsh Algorithmus
% m ist die Anzahl der Stützstellen, f die zu integrierende Funktion

% erstelle die betas
betas = (1:m-1) ./ sqrt(4*(1:m-1).^2 -1);
% baue A
A = diag(betas,-1)+diag(betas,1);
% berechne Eigenwerte/Eigenvectoren
[EV,EW] = eig(A);
% die Quadratur
val = 0;
for i = 1:m
    val = val + 2*EV(1,i)^2/norm(EV(:,i))^2*f(EW(i,i));
end
end


function val = Monte_Carlo(f,m,startpoint,endpoint)
% Realisiert die Approximation von \Int_a^b f(x)dx mit Hilfe des vereinfachten Monte-Carlo Verfahrens
% m ist die Anzahl der Stützstellen, f die zu integrierende Funktion, a = startpoint, b = endpoint

% erstelle gleichverteilte Stützstellen
Knoten = startpoint + (endpoint-startpoint)*rand(m,1); % Zufallszahlen aus [a,b] - Standardeinstellung von rand ist uniform - also gleichverteilt
% die Quadratur
val = (endpoint-startpoint)/m*sum(f(Knoten));
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