// Esercitazione 10
//
clear;
clf();
//
function J = integrale(a,b,f,k)
    pesi = ones(1,k+1);
    pesi(1) = 1/2; pesi($) = 1/2;
    nodi = linspace(a,b,k+1)';
    J = (b-a) * pesi * f(nodi) / k;
endfunction
//
function [k, M] = grafico(f,a,b,E,d2f)
    // Stima iniziale di M
    k = 10;
    nodi = linspace(a,b,k)';
    M = max(abs(d2f(nodi)));
    // Valore iniziale di k
    k = ceil((b - a) * sqrt(M / (8 * E)));
    nodi = linspace(a,b,k)';
    // Nuova stima di M
    nuovo_M = max(abs(d2f(nodi)));
    while nuovo_M > M, // finché M cresce...
          M = [M, nuovo_M];
          k = ceil((b - a) * sqrt(M($) / (8 * E)));
          nodi = linspace(a,b,k)';
          nuovo_M = max(abs(d2f(nodi)));
    end;
    plot2d(nodi,f(nodi),style = 2);
    xgrid(); xlabel('x');
    xtitle('Grafico di f(x) con errore < ' + msprintf('%3.2e',E) + '  (' + msprintf('%d',k) + ' punti)');
endfunction
//
// Approssimazione dell'integrale da 0 a pi
//
function y = f(x)
    y = sin(x);
endfunction
//
function y = d2f(x)
    y = - sin(x);
endfunction
//
a = 0;
b = %pi;
for k = 50:50:500,
    nodi = linspace(a,b,k+1)';
    J = integrale(a,b,f,k);
    M = max(abs(d2f(nodi)));
    L = M * (b-a)^3 / (8 * k^2);
    printf('\nJ = %10.9e   ,   errore = %3.2e < %3.2e\n',J,abs(2 - J),L);
end;
//abort;
//
// Grafico con E dato.
//
E = 1d-3;
subplot(211);
[k,M] = grafico(f,a,b,E);
//
nodi = linspace(a,b,k)';
x = linspace(a,b,10*k)';
errore = max(abs( f(x) - interp1(nodi,f(nodi),x) ));
printf('\nErrore nel grafico: %3.2e', errore);
subplot(212);
plot2d(x,abs( f(x) - interp1(nodi,f(nodi),x) ));
xgrid(); xlabel('x'); ylabel('errore');
