%%%%%%%%%%%%%%%%%MAIN5

close all
clear all

xp = -1 ;	% poczatek i 
xk =   1   ;		% koniec przedzialu w ktorym calkujemy
n =  3   ;		% ilosc podprzedzialow calkowania


x = linspace(xp,xk, 100) ;    %generuje wierszowy wektor 100 liczb rozłożonych równomiernie w przedziale od xp do xk (podział przedziału na 99 równych odcinków),

y = f1(x) ;
plot(x,y, 'r')


% % % Kwadratury Newtona-Cotesa
[I2] = f_rectI('f1', xp, xk, n) 
[I3] = f_trapI('f1', xp, xk, n) 
[I4] = f_SimpI('f1', xp, xk, n) 


tol = 10 ;	% zadana dokladnosc

I1a = quad('f1', xp, xk,tol,1); 
I1b = quadl('f1', xp, xk,tol,1);
I1c = trapz(x,y) 	


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%f_rectl

function [I] = f_rectI(f, xp, xk, n)

%%%%%%%%%%%%%%%%%%%%% ARGUMENTY %%%%%%%%%%%%%
% f 	    - zadana funkcja
% xp,xk 	- przedzial calkowania
% n         - ilosc podprzedzialow, na ktore dzielimy odcinek <a,b>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%<<<<<<<<<<<< ZWRACANE WARTOSCI <<<<<<<<<<<<<
% I  - przyblizenie calki
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

% Copyright MB
% Revision: 0.00		Date: April 20, 2011 		Time: 12:38:46 PM 
% Matlab ver: 7.5.0.342 (R2007b)


% % % szerokosc pojedynczego przedzialu
dx =  (xk-xp)/n  ;   % wzór 8

% % % wezly kwadratury; wzór 7
%linspace generuje wierszowy wektor n liczb rozłożonych równomiernie w zadanym przedziale
x_i = linspace(xp+0.5*dx, xk-0.5*dx, n) ;

% % % wartosci funkcji w wezlach kwadratury
f_i = eval([f '(x_i)']) ;   % wzór 9
 

I = 0 ;		% wartosc calki
for ii=1:n		% dla kazdego podprzedzialu: dodaj do tymczasowej wartosci calki, 
				% pole prostokata o wysokosci rownej wartosci funkcji w srodku ii-tego podprzedzialu
	I = I + f_i  ;   % wzór 10, 11
end





% % % wizualizacja metody prostokatow
xx= linspace(xp,xk,n+1) ;
figure, hold on
for ii=1:n
	plot( [xx(ii) xx(ii+1)], [f_i(ii) f_i(ii)] )
	plot( [xx(ii) xx(ii)], [0  f_i(ii)] )
	plot( [xx(ii+1) xx(ii+1)], [0  f_i(ii)] )
end

xx= linspace(xp,xk,10*n) ;
plot(xx, eval([f '(xx)']), 'r' )
title('Metoda prostokatow')
hold off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%f_Simpl

function [I] = f_SimpI(f, xp, xk, n)

%%%%%%%%%%%%%%%%%%%%% ARGUMENTY %%%%%%%%%%%%%
% f 	    - zadana funkcja
% xp, xk 	- przedzial calkowania
% n         - ilosc podprzedzialow, na ktore dzielimy odcinek <xp,xk>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%<<<<<<<<<<<< ZWRACANE WARTOSCI <<<<<<<<<<<<<
% I  - przyblizenie calki
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

% Copyright MB
% Revision: 0.00		Date: April 20, 2011 		Time: 12:38:46 PM 
% Matlab ver: 7.5.0.342 (R2007b)


% % % szerokosc pojedynczego przedzialu
dx = (xk-xp)/n ;
% % % wezly kwadratury
x_i = linspace(xp,xk, n+1+n) ;
% % % wartosci funkcji w wezlach kwadratury
f_i = eval([f ' (x_i)']) ;

% ilosc wezlow
%nn = n+n+1 ;

%    1     2     3     4      - numeracja przedzialow
% |--.--|--.--|--.--|--.--|
% 1  2  3  4  5  6  7  8  9   - numeracja wezlow


I = 0 ;
for ii=1:n	% dla kazdego podprzedzialu: dodaj do tymczasowej wartosci calki, 
				% pole pod parabola opisana na 3 wezlach ii-tego przedzialu
                          % wzór 28
	I = I + ( f_i(2*ii-2 +1) + 4*f_i(2*ii-1 +1) + f_i( 2*ii ) ) * dx/6 ;
end
  




% wizualizacja metody Simpsona
figure, hold on
for ii=1:n
	ind = (2*ii-1) ;
	xx= linspace(x_i(ind  ),x_i(ind+2), 10) ;
	
	A = [x_i(ind  )^2 x_i(ind  ) 1
		 x_i(ind+1)^2 x_i(ind+1) 1
		 x_i(ind+2)^2 x_i(ind+2) 1 ] ;
	d = [f_i(ind) f_i(ind+1) f_i(ind+2) ]' ;
	c = A\d ;
	
	
	plot( xx, c(1)*xx.^2 + c(2)*xx + c(3)  )
	plot( [x_i(ind)   x_i(ind)], [0  f_i(ind)] )
	plot( [x_i(ind+2) x_i(ind+2)], [0  f_i(ind+2)] )
end

zz = linspace(xp,xk,10*n) ;
plot(zz, eval([f '(zz)']), 'r' )
title('Metoda Simpsona')
hold off