// Esercitazione 11
//
// Confronto numerico dell'errore di ricostruzione con interpolazione polinomiale
// e con funzioni continue lineari a tratti.
//
// ****** Procedura SA ******
//
function [x,info] = SA(T,c)
//
// Funzione SA: per T matrice N x N triangolare inferiore invertibile e c colonna
// di N componenti, calcola la colonna soluzione del sistema di equazioni lineari
//
//       T x = c
//
// x: se T invertibile, la colonna x (di N numeri reali) approssima la soluzione
//    del sistema, altrimenti x = [].
// info = 0 se T non invertibile,
//      = 1 se T invertibile.
//
N = size(T, 'r');
//
// *** Verifica se T invertibile. 
//
if find(diag(T)==0) ~= [] then
   info = 0; x = [];
   error("SA: elementi nulli sulla diagonale.");
   else
       info = 1;
       //
       // *** Calcolo soluzione
       //
       x = zeros(N,1);
       //
       x(1) = c(1) / T(1,1);
       for k = 2:N,
           x(k) = (c(k) - T(k,1:k-1)*x(1:k-1)) / T(k,k);
       end;
end
endfunction
//
// ****** Procedura PolyInterpNewt ****** 
//
function yp = PolyInterpNewt(x,y,xp)
    // x,y: colonne di K numeri reali.
    // xp: colonna di M numeri reali.
    //
    // yp: colonna di M numeri reali. La componente yp(i) è il valore
    //     del polinomio interpolante in xp(i).
    //
    // Controllo argomenti
    if size(x,'c') > 1 then error('il primo argomento non è una colonna'); end;
    if ~isequal(size(x),size(y)) then
       error('i primi due argomenti non hanno la stessa dimensione'); end;
    if size(xp,'c') > 1 then error('il terzo argomento non è una colonna'); end;
    // Costruzione della matrice del sistema N per determinare la colonna b dei
    // coefficienti della forma di Newton.
    K = size(x,'r');
    N(:,1) = ones(K,1);
    for j = 2:K,
        N(j:K,j) = N(j:K,j-1) .* (x(j:K) - x(j-1));
    end;
    b = SA(N,y);
    // Calcolo dei valori del polinomio interpolante in xp(i).
    yp = b(K) * ones(xp);
    for j = K-1:-1:1,
        yp = b(j) + (xp - x(j)) .* yp;
    end;
endfunction
// ***************************************************************
function y = Fun(x)
    //y = exp(-x);
    y = sin(%pi * x);
endfunction
// numero campioni
K = 5; // K = 5,15,35,45,50,...
// per il grafico: n punti tra t(i) e t(i+1)
n = 10;
// dati
x = linspace(0,1,K)';
y = Fun(x);
// punti dove calcolare Fun e le ricostruzioni
xp = linspace(0,1,K + n*(K-1))';
// ricostruzione polinomiale
yp = PolyInterpNewt(x,y,xp);
// ricostruzione lineare a tratti
ypwl = interp1(x,y,xp);
// grafici dell'errore
clf(); plot2d(xp,[abs(Fun(xp)-ypwl),abs(Fun(xp)-yp)],style = [2,5]);
xlabel('x'); legend('lineare a tratti','polinomiale');
xgrid();
//abort;
// errore di ricostruzione max
printf('\n**** K = %d',K);
printf('\ninterp polin e(f) <= %.5e',%pi^(K) / factorial(K));
printf('\ninterp lin a tratti: e(f) <= %.5e', (%pi^2) / (8 * (K-1)^2));
//abort;
// precisione di macchina
u = 2^(-53);
// condizionamento ricostruzione polinomiale (caso camp uniforme):
//
//                                          2^K
//  max(|L(0,K)| + ... + |L(K,K)|) ~ ----------------
//                                   e (K-1) log(K-1)
//
printf('\nmax err interp polin = %.5e',u * (2^(K) / (%e * (K-1)*log(K-1))) );
