function dsi = demo_simpleint()
% Vypocet urciteho integralu - naivni implementace
%   (scitani po subintervalech, nasobeni h v kazdem kroku)
%
% Metody (vse s konstantni delkou kroku h):
%   - obdelnikove pravidlo - levostranne
%   - obdelnikove pravidlo - pravostranne
%   - obdelnikove pravidlo - midpoint
%   - lichobeznikove pravidlo
%   - Simpsonovo pravidlo

%% pocet intervalu (mel by byt sudy kvuli Simpsonovi)
nint = 10;    
%nint = 100;    
%nint = 1000;    
    
%% zadani integrovane funkce a integracnich mezi
problem = 1;

if (problem==1)
    a = 1;
    b = 5;
    f = inline('exp(x)','x');
    %% analyt. integral - pro kontrolu
    fint = inline('exp(x)','x');
elseif (problem==2)
    a = pi/6;
    b = pi/2;
    f = inline('sin(x)^3','x');
    %% analyt. integral - pro kontrolu
    fint = inline('-cos(x)*(2+sin(x)^2)/3','x');
end


%% sit uzlu x_i
h  = (b-a)/nint;
x  = (a:h:b);


left  = 0.0;
righ  = 0.0;
midp  = 0.0;
trap  = 0.0;
simp  = 0.0;

%% Levostr. pravostr. midpoint a lichob: jdeme po intervalech
for i=1:nint
    left = left + h * f(x(i));
    righ = righ + h * f(x(i+1));
    midp = midp + h * f((x(i)+x(i+1))/2);
    trap = trap + h * ( f(x(i)) + f(x(i+1)) ) / 2;
end;
%% Simpson - intervaly bereme po dvou
for i=1:2:nint-1
    simp = simp + h * ( f(x(i)) + 4*f(x(i+1)) + f(x(i+2))) / 3;
end;

%% Vypocet relativnich chyb
exact = fint(b)-fint(a);
eleft = abs((left-exact)) / max(abs(exact),1e-30);
erigh = abs((righ-exact)) / max(abs(exact),1e-30);
emidp = abs((midp-exact)) / max(abs(exact),1e-30);
etrap = abs((trap-exact)) / max(abs(exact),1e-30);
esimp = abs((simp-exact)) / max(abs(exact),1e-30);


%% vypis reseni
disp(['Integral od ' num2str(a) ' do ' num2str(b) ' je presne '...
      num2str(exact) '.']);
disp(['Numericky na ' num2str(nint+1) ' bodech:']);
disp(['  Levostranne obdel.  : ' num2str(left) ',   rel.ch. ' num2str(eleft)]);
disp(['  Pravostranne obdel. : ' num2str(righ) ',   rel.ch. ' num2str(erigh)]);
disp(['  Obdelnik (midpoint) : ' num2str(midp) ',   rel.ch. ' num2str(emidp)]);
disp(['  Lichobeznikove      : ' num2str(trap) ',   rel.ch. ' num2str(etrap)]);
disp(['  Simpsonovo          : ' num2str(simp) ',   rel.ch. ' num2str(esimp)]);

