function dr = demo_romberg()
   
    
%% integrovana funkce
f = inline('sin(x)*exp(cos(x))','x');
%% integracni meze
a = 0.0;
b = pi;

%% primitivni funkce (jen pro vykresleni)
fint = inline('-exp(cos(x))','x');


%%  slozene lichobeznikove pravidlo
function [integral] = lichobeznik(funkce,odkud,kam,krok)
    xArr     = odkud:krok:kam;
    integral = 0;
    for x=xArr
        % secteme vsechny funkcni hodnoty, ktere se opakuji dvakrat
    	integral = integral + funkce(x);
    end;
    % a odecteme poloviny kraju
    integral = integral - 0.5*funkce(xArr(1)) - 0.5*funkce(xArr(size(xArr,2)));
    integral = integral*krok;
end


%%  postupne spocteme integraly na 4 rozlisenich
%%  s kroky h0=0.1, h0/2, h0/4, h0/8
nint  = 4;
nref  = 4;
integ = zeros(nref,1);
h0    = (b-a)/nint;
disp(['Zakladni rozliseni ' num2str(nint) ...
      ' interv., krok h0= ' num2str(h0) ]);
for i=1:nref,
  h = h0/(2^(i-1));
  integ(i) = lichobeznik(f,a,b,h);
  disp(['Krok h= ' num2str(h) ' : int= ' num2str(integ(i))]);
end;


%%  hodnota integralu analyticky (jen pro vykresleni)
exact = fint(b) - fint(a);
disp(['Presna hodnota integralu je ' num2str(exact) ]);


%%  vykresleni jednotlivych integralu a extrapolaci
figure;

xmin = 1-nref/20;
xmax = nref+nref/20;
plot([xmin xmax], [exact exact], 'k--');
hold on;
plot(1:nref,integ,'*', 'Linewidth', 2);
xlim([xmin xmax]);
legend({'Presne (analyticke) reseni';...
        'Lichob. integral s krokem h_0/2^{n-1}'},'Location','SE')
xlabel('Zjemneni  n  (n=1 odpovida h_0)')
ylabel('(Vypoctena) hodnota integralu')
title('... stisknete ENTER');
pause

r1=4/3*integ(2)-1/3*integ(1);
disp(['Extrap. ze 2 bodu: int= ' num2str(r1) ]);
plot([2 2], [r1 r1], 'rd' , 'Linewidth', 2);
legend({'Presne (analyt.) reseni'; ...
        'Lichob. integral s krokem h_0/2^{n-1}'; ...
        'Romberg extrap. z n bodu'},'Location','SE')
pause;

r2=64/45*integ(3)-20/45*integ(2)+1/45*integ(1);
disp(['Extrap. ze 3 bodu: int= ' num2str(r2) ]);
plot([2 3], [r1,r2], 'rd' , 'Linewidth', 2);
r21 = 4/3*integ(3)-1/3*integ(2);
err_odh = abs(r2-r21);
disp(['  Odhad chyby extrap. ze 3 bodu: err= ' num2str(err_odh) ]);
pause;

r3 = 4096/2835*integ(4)-1344/2835*integ(3)+84/2835*integ(2)-integ(1)/2835;
disp(['Extrap. ze 4 bodu: int= ' num2str(r3,10) ]);
plot([2 3 4],[r1 r2 r3],'rd', 'Linewidth', 2)
hold on;
r32=64/45*integ(4)-20/45*integ(3)+1/45*integ(2);
err_od4 = abs(r3-r32);
disp(['  Odhad chyby extrap. ze 4 bodu: err= ' num2str(err_od4) ]);
err_skut = r3 - exact;
disp(['  Skutecna chyba extrap. ze 4 bodu: errskut= ' num2str(err_skut) ]);
title(' ');

hold off;

end
