

disp("testy kwadratur zlozonych");

function y=trap(f,a,b,N=200)
  %f uchwyt do funkcji (dla wektora x zwraca wektor y(k)=[f(x(1),..,f(x(N))]
  %a,b - konce odcinka calkowania
  %N - ilosc punktow minus jeden (zawieraja konce)
  x=linspace(a,b,N+1);
  fx=f(x);
  y=(0.5*(fx(1)+fx(N+1))+sum(fx(2:N)))*(b-a)/N;
end

function [y,tf]=testtrap(f,a,b,cf=quad(f,a,b),N0=10,M=8)
%testujemy kwad zl trapezow
%INPUT
%f uchwyt od funckji calkowanej
%a,b - konce odcinka calkowania
%cf - dokladna wartosc calki (default policzona quad()
N=10;
M=10;tf=h=zeros(M,1);
for k=1:M,
  h=(b-a)/N;
  tf(k)=trap(f,a,b,N);
  %y(k)/(h*h)
  N*=2;
end  
y=tf-cf;
end

a=0;b=1;
cf=-cos(1)+cos(0);
f=@sin; 
trap(f,a,b,N=200);
y=testtrap(f,a,b,cf)
y(1:M-1)./y(2:M)

a=0;b=1;
f=@(x) sqrt(x)
cf=(sqrt(b)*b-sqrt(a)*a)*2/3;
[y,tf]=testtrap(f,a,b,cf)
y(1:M-1)./y(2:M)


