
function demo_integrazione4(example)

%N=[1:10 20:10:50];
N=2.^(0:20);

do_plot=0;

if nargin < 1
    example=1;
end

switch example
    case 1
        a=-5; b=5;
        f=@(x) 1./(1+x.^2);
    case 2
        a=0; b=1;
        f=@(x) exp(x.^2);
     case 3
        a=0; b=1;
        f=@(x) x.^(1/2);
   case 4
        a=0; b=1;
        f=@(x) x.^(9/2);
    case 5
        a=0; b=pi;
        f=@(x) exp(x).*cos(x);
    case 6
        a=0; b=100;
        f=@(x) exp(-x.^2);
    case 7
        a=0; b=1;
        f=@(x) (x.^3).*sqrt(x);
end
% f=@(x) x.^3+2*x.^2+0.2*x+1;
ff=chebfun(f,[a,b]);
sol=sum(ff,[a,b]);

fprintf('\n \t SOL: %1.25e',sol);


% --------------------- code starts here ---------------------


% simpson composte
fprintf(' \n \t');

for ii=1:length(N)
    [x,w]=rettangolo_pm(N(ii),a,b); 
    I(ii,1)=compute_integral(f,x,w);
    L(ii,1)=length(x);
    [x,w]=trapezi_composta(N(ii),a,b); 
    I(ii,2)=compute_integral(f,x,w);
    L(ii,2)=length(x);
    [x,w]=simpson_composta(N(ii),a,b); 
    I(ii,3)=compute_integral(f,x,w);
    L(ii,3)=length(x);
end

for ii=1:length(N)
        fprintf('\n $%2.0f$ & $%1.1e$ & $%1.1e$ & $%1.1e$ &',...
            N(ii),abs(I(ii,1)-sol),abs(I(ii,2)-sol),abs(I(ii,3)-sol));
                fprintf('$%5.0f$ & $%5.0f$ & $%5.0f$ \\\\',...
                    L(ii,1),L(ii,2),L(ii,3));
end

fprintf('\n \n');

for ii=1:length(N-1)
    rat=abs(I(ii,:)-sol)./abs(I(ii+1,:)-sol);
        fprintf('\n $%2.0f$ & $%2.2f$ & $%2.2f$ & $%2.2f$ \\\\',...
            N(ii),rat(1),rat(2),rat(3));
end

fprintf('\n \n');


if do_plot == 1
clf
semilogy(L(:,1),abs(I(:,1)-sol),'k-','LineWidth',2,...
        'MarkerEdgeColor','k',...
        'MarkerFaceColor','k',...
        'MarkerSize',10)
    hold on;
semilogy(L(:,2),abs(I(:,2)-sol),'r','LineWidth',2,...
        'MarkerEdgeColor','r',...
        'MarkerFaceColor','r',...
        'MarkerSize',10)
semilogy(L(1:11,3),abs(I(1:11,3)-sol),'c-','LineWidth',2,...
        'MarkerEdgeColor','c',...
        'MarkerFaceColor','c',...
        'MarkerSize',10)
legend('Punto di mezzo composta','Trapezi composta','Cav.-Simpson composta')
hold off;
end



function I=compute_integral(f,x,w)

fx=feval(f,x);
I=w'*fx;




function [x,w]=trapezi_composta(N,a,b)

% FORMULA DEI TRAPEZI COMPOSTA.

% INPUT:
% N: NUMERO SUBINTERVALLI.
% a, b: ESTREMI DI INTEGRAZIONE.

% OUTPUT:
% x: NODI INTEGRAZIONE.
% w: PESI INTEGRAZIONE (INCLUDE IL PASSO!).

h=(b-a)/N;             % PASSO INTEGRAZIONE.
x=a:h:b; x=x';         % NODI INTEGRAZIONE.
w=ones(N+1,1);         % PESI INTEGRAZIONE.
w(1)=0.5; w(N+1)=0.5;
w=w*h;



function [x,w]=simpson_composta(N,a,b)
% FORMULA DI SIMPSON COMPOSTA.
% INPUT:
% N: NUMERO SUBINTERVALLI.
% a, b: ESTREMI DI INTEGRAZIONE.
% OUTPUT:
% x: INTEGRAZIONE.
% w: PESI INTEGRAZIONE (INCLUDE IL PASSO!).
h=(b-a)/N;             % AMPIEZZA INTERVALLO.
x=a:(h/2):b; x=x';     % NODI INTEGRAZIONE.
w=ones(2*N+1,1);           % PESI INTEGRAZIONE.
w(3:2:2*N-1,1)=2*ones(length(3:2:2*N-1),1);
w(2:2:2*N,1)=4*ones(length(2:2:2*N),1);
w=w*h/6;




function [x,w]=rettangolo_pm(N,a,b)
% FORMULA DI SIMPSON COMPOSTA.
% INPUT:
% N: NUMERO SUBINTERVALLI.
% a, b: ESTREMI DI INTEGRAZIONE.
% OUTPUT:
% x: INTEGRAZIONE.
% w: PESI INTEGRAZIONE (INCLUDE IL PASSO!).
h=(b-a)/N;             % AMPIEZZA INTERVALLO.
hV=(h/2)*(2*(1:N)-1);
x=a+hV; x=x';     % NODI INTEGRAZIONE.
w=h*ones(length(x),1);           % PESI INTEGRAZIONE.


