function [int_fce,I_O,I_L,I_S]=integral_vypocet_newton(f, a, b, N, zobraz);
%
% Vypocte numericky urcity integral symbolicky zadane funkce f 
% na intervalu <a,b> s poctem deleni N pomoci Newton-Cotesovych vzorcu
% I_O ... obdelnikove pravidlo
% I_L ... lichobeznikove pravidlo
% I_S ... Sipmsonovo pravidlo


if mod(N,2)
  N=N+1;
end;
krok=(b-a)/N;
deleni=krok/100;

if zobraz==1
 clc;
 disp('------------------------------------------------------------------');
 disp(sprintf('  Numericky vypocet urciteho integralu funkce f(x)=%s',f));
 disp(sprintf('  na intervalu <%f,%f> s poctem deleni N=%d',a,b,N));
 disp('------------------------------------------------------------------');
end;

int_fce=eval(int(f,a,b));
I_O=0;
I_L=0;
I_S=0;

% obdelnikove pravidlo
for i=1:N
  x=a+(i-0.5)*krok;
  ff(i)=eval(f); 
  I_O=I_O+ff(i)*krok;  
end;

% lichobeznikove pravidlo
for i=1:N
  x=a+(i-1)*krok;
  levy=eval(f);
  x=a+i*krok;
  pravy=eval(f);
  ff(i)=(pravy+levy)/2;
  I_L=I_L+ff(i)*krok;  
end;

% Simpsonovo pravidlo
for i=1:N/2
  x=a+2*(i-1)*krok;
  levy=eval(f);
  x=a+(2*i-1)*krok;
  stred=eval(f);
  x=a+2*i*krok;
  pravy=eval(f);
  ff(i)=(pravy+4*stred+levy)/3;
  I_S=I_S+ff(i)*krok;  
end;


if zobraz==1
 figure(1);
 clf;
 hold on;
 zoom on;
 xx=a:deleni:b;
 i=1;
 for x=a:deleni:b
   fcef(i)=eval(f);
   i=i+1;
 end;

 max_y=max(fcef);
 min_y=min(fcef);
 sy=max_y-min_y;
 sx=b-a;

 axis([a-sx/10 b+sx/10 min_y-.1*sy max_y+.1*sy]);
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 texty(1,:)=sprintf('zadana funkce                   ');
 legend(texty,0);
 pause;
 clf;
 hold on;
 zoom on;
 axis([a-sx/10 b+sx/10 min_y-.1*sy max_y+.1*sy]);

 disp(' ');
 disp(sprintf(' Presna hodnota integralu je  %f ',int_fce));
 disp(' ');
 barva=[0.75 0.75 0.85];
 fill([xx,b,a],[fcef,0,0],barva)
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,0);
 title('Presna hodnota integralu');
 pause;

 disp(sprintf([' Priblizna hodnota integralu',...
                '                                  chyba ']));
 disp(sprintf(...
 '     - pomoci obdelnikoveho pravidla       %f          %f',I_O,I_O-int_fce));
 clf;
 hold on;
 zoom on;
 axis([a-sx/10 b+sx/10 min_y-.1*sy max_y+.1*sy]);
 barva=[0.75 0.85 0.75];
 for i=1:N
   levy=a+(i-1)*krok;
   pravy=a+i*krok;
   x=(pravy+levy)/2;
   pom=eval(f); 
   fill([levy,levy,pravy,pravy],[0 pom pom 0], barva);
 end;
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,0);
 title('Obdelnikove pravidlo');
 pause;              

 disp(sprintf(...
 '     - pomoci lichobeznikoveho pravidla    %f          %f',I_L,I_L-int_fce));
 clf;
 hold on;
 zoom on;
 axis([a-sx/10 b+sx/10 min_y-.1*sy max_y+.1*sy]);
 barva=[0.85 0.75 0.75];
 for i=1:N
   levy=a+(i-1)*krok;
   pravy=a+i*krok;
   x=levy;
   pom1=eval(f); 
   x=pravy;
   pom2=eval(f); 
   fill([levy,levy,pravy,pravy],[0 pom1 pom2 0], barva);
 end;
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,0);
 title('Lichobeznikove pravidlo');
 pause;              

 disp(sprintf(...
 '     - pomoci Simpsonova pravidla          %f          %f',I_S,I_S-int_fce));
 clf;
 hold on;
 zoom on;
 axis([a-sx/10 b+sx/10 min_y-.1*sy max_y+.1*sy]);
 barva=[0.85 0.65 0.75];
 for i=1:N/2
   x=a+2*(i-1)*krok;
   xlevy=x;
   levy=eval(f);
   x=a+(2*i-1)*krok;
   xstred=x;
   stred=eval(f);
   x=a+2*i*krok;
   xpravy=x;
   pravy=eval(f);
   [xparabola,parabola,koef]=polynom([xlevy,xstred,xpravy],[levy,stred,pravy],deleni,0);
   fill([xparabola,xpravy,xlevy],[parabola,0,0], barva);
 end;
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,0);
 title('Simpsonovo pravidlo');
 disp(' ');
 pause;
 close(1);
end;









