//
// Esercitazione 10
//
// Realizzazione della funzione che assegnati i dati
//
//   (x(i) , y(i))    ,    i = 1,...,K
//
// determina il polinomio interpolante p e lo calcola in
//
//   xp(1), ... ,xp(M)
//
// (1) Sia [a,b] = [0,1]. Data una funzione Fun, disegnare il grafico di Fun
//     e del polinomio p(x) che interpola i valori di Fun in x(1), ... ,x(K)
//     su uno stesso piano cartesiano. Verificare che se Fun è un polinomio
//     di grado non superiore a K-1, si ha p(x) = Fun(x) per ogni x in [0,1].
//
// (2) Disegnare il grafico degli elementi della base di Lagrange relativi a
//     x(1), ... ,x(K).
//
// (3) Per Fun(x) = exp(-x), studiare l'errore di ricostruzione con polinomio
//     interpolante per K = 5,10,20. Per K = 20, confrontare l'errore ri
//     ricostruzione che si ottiene con e senza perturbazione dei campioni.
//
//
// ****** 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 = (x - 1) .* (2*x - 1) .* (8*x - 1) .* (10*x - 8);
endfunction
// Intervallo [a,b]
a = 0; b = 1;
// Numero di dati da interpolare
K = 20;
// Numero di valori da calcolare tra due dati
n = 10;
// Dati da interpolare (colonne!)
x = linspace(a,b,K)';
y = Fun(x);
// Ascisse dove valutare il polinomio interpolante
// (includono quelle dei dati!)
xp = linspace(a,b,K + (K-1)*n)';
// Calcolo del valore del polinomio interpolante
yp = PolyInterpNewt(x, y, xp);
// disegno
clf();
// grafico del polinomio interpolante
plot2d(xp,yp, style = 2);
// grafico della funzione campionata
plot2d(xp,Fun(xp), style = 5);
// dati da interpolare
plot2d(x',y', style = -2);
xgrid(); xlabel('x');
legend('p(x)','Fun(x)','(x(i),y(i))','in_upper_right');
abort;
// ************ Errore di ricostruzione: porre Fun(x) = exp(-x)
halt();
clf();
// grafico di | Fun(x) - p(x) |
plot2d(xp, abs(yp - Fun(xp)), style = 5);
xgrid(); xlabel('x');
ylabel('| Fun(x) - p(x) |');
// errore di ricostruzione (stima)
printf('\nErr di ricostr (%d campioni):                             %.5e\n',...
       K,max(abs(Fun(xp) - yp)));
abort;
halt();
// Perturbazione massima delle ordinate dei campioni e dati perturbati
err = 1d-14;
yerr = y + err*(2*rand(x)-1);
// Calcolo del valore del polinomio interpolante
yperr = PolyInterpNewt(x,yerr,xp);
// grafico di | Fun(x) - p_err(x) |
plot2d(xp, abs(yperr - Fun(xp)), style = 5);
xgrid(); xlabel('x');
ylabel('| Fun(x) - p_err(x) |');
// errore di ricostruzione (stima)
printf('\nErr di ricostr con pert di misura %.2e (%d campioni): %.5e',...
       err,K,max(abs(yperr - Fun(xp))));
