%% Romberg-Quadratur
% Wir untersuchen die Romberg Extrapolation angewandt auf die summierte Trapezregel, zur Annäherung eines Integrals   I=\int_0^1 f(x)\, dx.
% Sei
%      T(h)=h\big[\frac12 f(0) +f(t_1)+\ldots+f(t_{n-1})+\frac12 f(1)\big] , mit t_j=jh, h=1/n,
% die summierte Trapezregel, und T_{i,0}:=T(2^{-i}h_0), i=0,\ldots,6, mit h_0:=0.25. Ausgehend von diesen 7 Approximationen des Integrals I wird ein  Romberg-Schema berechnet (siehe Buch):
%                 T_{i,j}:=\frac{4^j T_{i,j-1}-T_{i-1,j-1}}{4^j-1},\quad 1 \leq j  \leq i \leq 6.
% Als Beispiele betrachten wir:
%     1.   f(x)=0.25\pi^2 x\cos(0.5 \pi x)+0.5\pi e^{0.5\pi x}    (Beispiel Buch umskaliert auf [0,1])
%     2.   f(x)= x e^x      (glatte Funktion)
%     3.  f(x)= \sqrt{x}       (nicht glatt an der Stelle x=0)
%  Für diese Bespiele kann  analytisch der exakte Wert des Integrals I berechnet werden.
%  Durch Angabe von choice \in \{1,2,3\} wird eins dieser Bespiele gewählt.

choice=2
[f,Iexakt]= Beispiel(choice);
fplot(f,[0 1]);
legend('Funktion f', 'Location','best');

% Die Trapezsummen T_{i,0} und zugehörige Fehler |T_{i,0}-I|, i=0,\ldots,6, werden berechnet und gezeigt.
% Ausgehend von diesem Approximationen werden die Romberg-Approximationen T_{i,j} und entsprechende Fehler |T_{i,j}-I| berechnet.
%% Aktivität:
% Sie können den Parameter choice variieren. Beobachten Sie das Fehlerverhalten im Romberg-Schema. Weshalb liefert das Romberg-Schema für das Beispiel choice=3 keine (signifikant) genaueren Approximationen?

kmax=7;
Trapezsumme=zeros(kmax,1);
for k=1:kmax
    n=2^(k+1);
    h=1/n;
    T= 0.5*(f(0) + f(1));
    for i=1:n-1
        T= T + f(i*h);
    end
    Trapezsumme(k)=h*T;
end 
fprintf('\nErgebnisse der summierten Trapezregel:\n');
format long;
Trapezsumme
format shorte;
Fehler= abs(Trapezsumme - Iexakt)
format long;
disp('----------------');
fprintf('\nRomberg-Schema:\n');
T= extrapolation(Trapezsumme)
fprintf('\nGenauigkeit der Approximationen im Romberg-Schema:\n');
format shorte;
Fehler= Error(T,Iexakt)

%

function [f,Iexakt]= Beispiel(bsp)
switch bsp
    case 1
        f= @(x) 0.25*pi^2*x.*cos(0.5*pi.*x) + 0.5*pi*exp(0.5*pi.*x);
        Iexakt= pi/2 + exp(0.5*pi) - 2;
    case 2
        f= @(x) x.*exp(x);
        Iexakt= 1;
    case 3
        f= @(x) sqrt(x);
        Iexakt= 2/3;
end

end
 
 
function T=extrapolation(t)
    m=length(t);
    T=NaN*zeros(m,m);
    T(:,1)=t;
    for j=2:m
        for i=j:m
            T(i,j)= (4^(j-1)*T(i,j-1) - T(i-1,j-1))/(4^(j-1)-1);
        end
    end
end


function F=Error(T,Iexakt)
    m=length(T);
    F=NaN*zeros(m,m);
    for j=1:m
        for i=j:m
            F(i,j)= abs(T(i,j) - Iexakt);
        end
    end
end
 
     

%
%