//
// Esercitazione 9
//
// * Realizzare una funzione che assegnati i dati
//
//    (x(i) , y(i))    ,    i = 1,...,K
//
//  determina i valori del polinomio interpolante in
//
//    xp(1), ... ,xp(M)
//
// * Verificare che interpola esattamente polinomi di grado non superiore
//   a K-1.
// * Usarla per disegnare grafici degli elementi della base di Lagrange.
//
// * Data Fun, studiare l'errore di ricostruzione con il polinomio interpolante.
//
//
// ****** Procedura SA ******
//
function x = SA(T,c)
//
// Per T matrice N x N triangolare inferiore invertibile e c colonna
// di N componenti, calcola la colonna x (di N numeri reali) che approssima la
// soluzione del sistema di equazioni lineari:
//
//       T z = c
//
N = size(T, 'r');
//
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
endfunction
//
// Procedura PolyInterpNewt ***************************************
//
function yp = PolyInterpNewt(x,y,xp)
    // x,y: colonne di K numeri reali
    // xp, yp: colonne di M numeri reali
    K = length(x);
    N = zeros(K,K); // matrice del sist
    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;
    c = SA(N,y); // coefficienti pol interp
    // calcolo valori del pol interp
    yp = c(K);
    for j = K-1:-1:1,
        yp = c(j) + (xp - x(j)) .* yp;
    end;
endfunction
//                          ***************************************
function y = Fun(x)
//    y = exp(-x);
//    y = exp(-x) .* sin(3*x);
//    y = 2 + x - 5*x .^2 + 4*x .^3;
endfunction
//
a = 0; b = 1;
K = 5;  // 10; 20;
n = 50; // punti aggiuntivi per il disegno
x = linspace(a,b,K)';
y = [0;0;0;1;0]; // base di Lagrange: l(4)
//y = Fun(x);
xp = linspace(a,b,K + (K-1)*n)';
yp = PolyInterpNewt(x, y, xp);
clf();
plot2d(xp,yp, style = 2);
//plot2d(xp,Fun(xp), style = 5);
plot2d(x,y, style = -2);
//plot2d(xp, abs(yp - Fun(xp)), style = 5);
xgrid(); xlabel('x');
legend('p_newt','(x(i),y(i))','in_upper_left');
//legend('p(x)','Fun(x)','(x(i),y(i))','in_upper_right');
//ylabel('| r(c(f)) - f |');
//printf('err di ricostr (con %d campioni): %.4e',K,max(abs(yp - Fun(xp))))
