
%kwadratura prostokatow zlozona podzial rownomierny
1==1;
function [c,h]=prost(f,a,b,N=1000)
%obliczamy przyblizenie calki zf na [a,b] zlozona kawdratura
%prostokatow \sum_{k=0}^(N-1) f(a+(k+0.5)*h)h h=(b-a)/N
%INPUT
%f uchwyt do funkcji (zakladamy ze f(x) dla x wektora jest wektorem f(x(k))
%a,b konce odcinka a<b
%N - ilosc pktow (domyslnie 1000)
%OUTPUT
%c - przyblizenie calki zlozona kwad. prost.
x=linspace(a,b,N+1);
h=(b-a)/N;
xh=x(1:N)+0.5*h;
fxh=f(xh);
c=sum(fxh)*h;
end

prost(@(x) x.*x,0,1)-1/3
disp("testowanie rzedu kwadratury zlozonej prostokatow dla x^(3/2) f kalcy C^1 ale nie C^2");

function [er,rat]=testprost(f,a,b,c=quad(f,a,b),M=7,N0=10)
 %testujemy rzad zlozonej kwad prostokatow
 %c=quad(f,a,b); %if we do not know the exact value of the integral
%it may approximated by a octave black box integral solver
 N=N0;
 er=zeros(M,4);

 for k=1:M,
  [p,h]=prost(f,a,b,N);
  ih=1/h;
  er(k,:)=[N,c-p,(c-p)*ih,(c-p)*ih*ih];
  N*=2;
 end  
  rat=er(1:end-1,2)./er(2:end,2);   

end

[er,rat]=testprost(@(x) x.*x,0,1,1/3,M=7,N0=10);
[rat,log2(rat)]
%sin na [0,1]
prost(@sin,0,1)+cos(1)-1


[er,rat]=testprost(@sin,0,b=1,1-cos(b),M=7,N0=10)

disp("testowanie rzedu kwadratury zlozonej prostokatow dla sin(x) na [0,1]");
disp("pierwsza kol: rat(N)=er(N)/er(2N); 2ga: log2(rat(N)) -rzad zboeznosci"); 
[rat,log2(rat)]

%sqrt(x) na [0,1] \int_0^1 sqrt(x) = 2/3
[er,rat]=testprost(@(x) sqrt(x),0,1,2/3,M=7,N0=10);
disp("testowanie rzedu kwad zlozonej prost dla x^(1/2) na [0,1]  funk klasy C");
[rat,log2(rat)]

%teraz x^(3/2) calak na [0,1] = 2/5
[er,rat]=testprost(@(x) x.*sqrt(x),0,1,2/5,M=7,N0=10);
disp("testowanie rzedu kwad. zlozonej prostokatow dla x^(3/2) na [0,1] ");
[rat,log2(rat)]


%x^(1/10) na [0,1] \int_0^1 sqrt(x) = 1/1.1
[er,rat]=testprost(@(x) x.^0.1,0,1,1/(1.1),M=7,N0=10);
disp("testowanie rzedu kwad zlozonej prost dla x^(1/10) na [0,1]  funk klasy C");
[rat,log2(rat)]