function [int_fce,I_G]=integral_vypocet_gauss(f, a, b, N, pocet, zobraz);
%
% Vypocte numericky urcity integral symbolicky zadane funkce f 
% na intervalu <a,b> s poctem deleni N pomoci Gaussovych vzorcu
% I_G ... vektor s vysledky 
%         (Gaussovy zakladni vzorce s uzly 1...pocet, (max 3)

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_G=zeros(min([pocet,3]),1);

% Gaussovy kvadraturni vzorce s 1 uzlem na zakladanim intervalu
for i=1:N
  xlevy=a+(i-1)*krok;
  xpravy=a+i*krok;
  w1=krok;
  x=(xlevy+xpravy)/2;
  f1=eval(f);
  I_G(1)=I_G(1)+w1*f1;
end;

if pocet>1
 % Gaussovy kvadraturni vzorce s 2 uzly na zakladanim intervalu
 for i=1:N
  xlevy=a+(i-1)*krok;
  xpravy=a+i*krok;
  w1=krok/2;
  w2=krok/2;
  x=(xlevy+xpravy)/2-(1/3)^(0.5)*(xpravy-xlevy)/2;
  f1=eval(f);
  x=(xlevy+xpravy)/2+(1/3)^(0.5)*(xpravy-xlevy)/2;
  f2=eval(f);
  I_G(2)=I_G(2)+w1*f1+w2*f2;
 end;
end;

if pocet>2
 % Gaussovy kvadraturni vzorce s 3 uzly na zakladanim intervalu
 for i=1:N
  xlevy=a+(i-1)*krok;
  xpravy=a+i*krok;
  w1=5/18*krok;
  w2=8/18*krok;
  w3=5/18*krok;
  x=(xlevy+xpravy)/2-(3/5)^(0.5)*(xpravy-xlevy)/2;
  f1=eval(f);
  x=(xlevy+xpravy)/2;
  f2=eval(f);
  x=(xlevy+xpravy)/2+(3/5)^(0.5)*(xpravy-xlevy)/2;
  f3=eval(f);
  I_G(3)=I_G(3)+w1*f1+w2*f2+w3*f3;
 end;
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 Gausova vzorce s 1 uzlem       %f       %f',I_G(1),I_G(1)-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
   xlevy=a+(i-1)*krok;
   xpravy=a+i*krok;
   w1=krok;
   x=(xlevy+xpravy)/2;
   f1=eval(f);
   fill([xlevy,xlevy,xpravy,xpravy],[0 f1 f1 0], barva);
 end;
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,0);
 title('Gaussuv kvadraturni vzorec s  1  uzlem');
 for i=1:N
   x=a+(i-0.5)*krok;
   plot(x,0,'bx');
 end;
 pause;              

if pocet >1
 disp(sprintf(...
 '     - pomoci Gausova vzorce s 2 uzly        %f       %f',I_G(2),I_G(2)-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
   xlevy=a+(i-1)*krok;
   xpravy=a+i*krok;
   w1=krok/2;
   w2=krok/2;
   x=(xlevy+xpravy)/2-(1/3)^(0.5)*(xpravy-xlevy)/2;
   f1=eval(f);
   x=(xlevy+xpravy)/2+(1/3)^(0.5)*(xpravy-xlevy)/2;
   f2=eval(f);
   fill([xlevy,xlevy,(xlevy+xpravy)/2,(xlevy+xpravy)/2],...
        [0 f1 f1 0], barva);
   fill([(xlevy+xpravy)/2,(xlevy+xpravy)/2,xpravy,xpravy],...
        [0 f2 f2 0], barva);
 end;
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,0);
 title('Gaussuv kvadraturni vzorec s  2  uzly');
 for i=1:N
   xlevy=a+(i-1)*krok;
   xpravy=a+i*krok;
   w1=krok/2;
   w2=krok/2;
   x=(xlevy+xpravy)/2-(1/3)^(0.5)*(xpravy-xlevy)/2;
   f1=eval(f);
   x=(xlevy+xpravy)/2+(1/3)^(0.5)*(xpravy-xlevy)/2;
   f2=eval(f);
   [fmin,ind]=min(abs([f1,f2]));
   if ind==1
     fmin=sign(f1)*fmin;
   else
     fmin=sign(f2)*fmin;
   end;
   if f1*f2<0 
     fmin=0;
   end;
   obj=line([(xlevy+xpravy)/2,(xlevy+xpravy)/2],[0 fmin]);
   set(obj,'Color',barva);
   set(obj,'LineStyle','--');
   x=(xlevy+xpravy)/2-(1/3)^(0.5)*(xpravy-xlevy)/2;
   plot(x,0,'bx');
   x=(xlevy+xpravy)/2+(1/3)^(0.5)*(xpravy-xlevy)/2;
   plot(x,0,'bx');
 end;
 pause;              
end;

if pocet>2
 disp(sprintf(...
 '     - pomoci Gausova vzorce s 3 uzly        %f       %f',I_G(3),I_G(3)-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
   xlevy=a+(i-1)*krok;
   xpravy=a+i*krok;
   w1=5/18*krok;
   w2=8/18*krok;
   w3=5/18*krok;
   x=(xlevy+xpravy)/2-(3/5)^(0.5)*(xpravy-xlevy)/2;
   f1=eval(f);
   x=(xlevy+xpravy)/2;
   f2=eval(f);
   x=(xlevy+xpravy)/2+(3/5)^(0.5)*(xpravy-xlevy)/2;
   f3=eval(f);
   fill([xlevy,xlevy,xlevy+w1,xlevy+w1],...
        [0 f1 f1 0], barva);
   fill([xlevy+w1,xlevy+w1,xpravy-w3,xpravy-w3],...
        [0 f2 f2 0], barva);
   fill([xpravy-w3,xpravy-w3,xpravy,xpravy],...
        [0 f3 f3 0], barva);
 end;
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,0);
 title('Gaussuv kvadraturni vzorec s  3  uzly');
 for i=1:N
   xlevy=a+(i-1)*krok;
   xpravy=a+i*krok;
   w1=5/18*krok;
   w2=8/18*krok;
   w3=5/18*krok;
   x=(xlevy+xpravy)/2-(3/5)^(0.5)*(xpravy-xlevy)/2;
   f1=eval(f);
   x=(xlevy+xpravy)/2;
   f2=eval(f);
   x=(xlevy+xpravy)/2+(3/5)^(0.5)*(xpravy-xlevy)/2;
   f3=eval(f);
   [fmin,ind]=min(abs([f1,f2]));
   if ind==1
     fmin=sign(f1)*fmin;
   else
     fmin=sign(f2)*fmin;
   end;
   if f1*f2<0 
     fmin=0;
   end;
   obj=line([xlevy+w1,xlevy+w1],[0 fmin]);
   set(obj,'Color',barva);
   set(obj,'LineStyle','--');
   [fmin,ind]=min(abs([f2,f3]));
   if ind==1
     fmin=sign(f2)*fmin;
   else
     fmin=sign(f3)*fmin;
   end;
   if f2*f3<0 
     fmin=0;
   end;
   obj=line([xpravy-w3,xpravy-w3],[0 fmin]);
   set(obj,'Color',barva);
   set(obj,'LineStyle','--');
   x=(xlevy+xpravy)/2-(3/5)^(0.5)*(xpravy-xlevy)/2;
   plot(x,0,'bx');
   x=(xlevy+xpravy)/2;
   plot(x,0,'bx');
   x=(xlevy+xpravy)/2+(3/5)^(0.5)*(xpravy-xlevy)/2;
   plot(x,0,'bx');
 end;
 disp(' ');
 pause;
end;
close(1);
end;









