//
// Esercitazione 7
//
// * Realizzare una funzione che assegnati i dati
//
//    (x(i) , y(i))    ,    i = 1,...,k+1
//
//  determina i valori del polinomio interpolante in
//
//    xp(1), ... ,xp(m)
//
// * Verificare che interpola esattamente polinomi di grado non superiore
//   a k.
// * 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+1 numeri reali
    // xp, yp: colonne di m numeri reali
    k = length(x) - 1;
    N = zeros(k+1,k+1); // matrice del sist
    N(:,1) = ones(k+1,1);
    for j = 2:k+1,
        N(j:k+1,j) = N(j:k+1,j-1) .* (x(j:k+1) - x(j-1));
    end;
    b = SA(N,y); // coefficienti pol interp
    // calcolo valori del pol interp
    yp = b(k+1);
    for j = k:-1:1,
        yp = b(j) + (xp - x(j)) .* yp;
    end;
endfunction
//                          ***************************************
function y = Fun(x)
    y = exp(-x) .* sin(3*x);
//    y = 2 + x - 5*x .^2 + 4*x .^3;
endfunction
//
a = 0; b = 1; // intervallo
k = 15;
n = 50; // punti aggiuntivi per il disegno
x = linspace(a,b,k+1)';
//y = [0;0;0;1;0]; // base di Lagrange: l(4)
y = Fun(x);
xp = linspace(a,b,k+1 + k*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+1,max(abs(yp - Fun(xp))))
