%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ATTENZIONE!!!!! NON MODIFICARE NULLA DI QUESTO SCRIPT
%====================================================================
% Script per il Metodo dei trapezi composto
% Necessita della Function trapezi.m
%====================================================================

clear all
file = fopen('tabella_trapezi.txt','w');

disp('=======================================================')
fun = 'sin(x.^2)';
fprintf('Funzione integranda = %s \n', fun);
fprintf(file,'Funzione integranda = %s \n', fun);
f = inline(fun);
a = 0;
b = pi/4;
fprintf('Intervallo di integrazione = [ %d , %18.16f ] \n', a,b);
fprintf(file,'Intervallo di integrazione = [ %d , %18.16f ] \n', a,b);
toll = 1.e-8;
fprintf('Tolleranza toll = %2.0e \n', toll);
fprintf(file,'Tolleranza toll = %2.0e \n', toll);
f = inline(fun);
m = 3;
fprintf('Suddivisione iniziale m = %d \n',m);
fprintf(file,'Suddivisione iniziale m = %d \n',m);
h = (b-a)/3;
fprintf('Passo iniziale h = %10.6f \n', h);
fprintf(file,'Passo iniziale h = %10.6f \n', h);
nmax = 15;
fprintf('Numero massimo di iterazioni nmax = %d \n', nmax);
fprintf(file,'Numero massimo di iterazioni nmax = %d \n', nmax);
disp('=======================================================')
% Parte esecutiva

[intvec(1),h] = trapezi (f,a,b,m);
fprintf('%s \t %s \t\t\t %s \t\t %s \n\n','   m','h','  Integrale',' Test');
fprintf(file,'\n %s \t %s \t\t\t %s \t\t %s \n\n','   m','h','  Integrale',' Test');
fprintf('%4.0f %16.5e %19.16f \n',[m,h,intvec(1)]);
fprintf(file,'%4.0f %16.5e %19.16f \n',[m,h,intvec(1)]);
inold=intvec(1);
err = toll+1;
n = 1;
while (err > toll) && (n<nmax)
    n = n+1;
    m = m*2;
    [intvec(n),h] = trapezi (f,a,b,m);
    err = abs(inold-intvec(n));
    fprintf('%4.0f %16.5e %19.16f %11.4e \n',[m,h,intvec(n),err]);
    fprintf(file,'%4.0f %16.5e %19.16f %11.4e \n',[m,h,intvec(n),err]);
    inold = intvec(n);
end   
if n==nmax
    disp(' ')
    disp('*** ATTENZIONE RAGGIUNTO IL NUMERO MASSIMO DI APPROSSIMAZIONI ***')
end
q = quadl(f,a,b,1e-15);

fprintf('\nValore integrale di riferimento q = %20.16f \n',q);
fprintf(file,'\nValore integrale di riferimento q = %20.16f \n',q);
fclose(file);


% figura
semilogy(1:n,abs(intvec-q),'r-')
title('Metodo dei trapezi')
xlabel('numero di iterazioni (raddoppio dei nodi)')
ylabel('norma dell''errore')



