//
// Esercitazione 8
//
// * 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)
//
//
// ****** 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
    // Determina la forma di Newton del polinomio che interpola
    // i dati (x(j),y(j)), j = 1:k+1 e poi ne calcola il valore
    // in yp(r), r = 1:m.
    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 in forma di Newt
    // calcolo dei valori del pol interp con metodo tipo Horner
    yp = b(k+1);
    for j = k:-1:1,
        yp = b(j) + (xp - x(j)) .* yp;
    end;
endfunction
//    *************************************************************

a = 0; b = 1; // intervallo

// *** Disegnare grafici degli elementi della base di Lagrange
//     e il grafico di
//                     sum_j=0:k  | l(j) |

//k = 4;
//n = 50; // punti aggiuntivi per il disegno
//x = linspace(a,b,k+1)';
//IdMat = eye(k+1,k+1);
//Y = [];
//for j = 1:k+1,
//    y = IdMat(:,j);
//    xp = linspace(a,b,k+1 + k*n)';
//    yp = PolyInterpNewt(x, y, xp);
//    Y(:,j) = yp;
//end;
//    legenda = 'l(1)';
//    for j = 2:k+1,
//        legenda = [legenda, 'l(' + string(j) + ')'];
//    end;
//clf();
//plot(xp,Y,'thickness',2);
//plot(x,0,'bko','thickness',3);
//xgrid();
//xlabel('$x$','fontsize',4);
//title('Base di Lagrange','fontsize',3);
//legend(legenda,'ist camp','out_upper_left');

// *** Il polinomio interpolante interpola esattamente polinomi
//     di grado non superiore a k.
//
//
//function y = Fun(x)
//    y = 2 + x - 5*x .^2 + 4*x .^3;
//endfunction
//k = 4;
//n = 50; // punti aggiuntivi per il disegno
//x = linspace(a,b,k+1)';
//y = Fun(x);
//xp = linspace(a,b,k+1 + k*n)';
//yp = PolyInterpNewt(x, y, xp);
//clf();
//plot(xp,yp, 'b','thickness',2);
//plot(xp,Fun(xp), 'r','thickness',2);
//plot(x,y, 'bkx','thickness',2);
//xgrid(); 
//xlabel('$x$','fontsize',4);
//legend('p(x)','Fun(x)','(x(i),y(i))','in_upper_right');
//title('Ricostruzione esatta di un polinomio','fontsize',3);

//
//  *** grafico sum_j=0:k  | l(j) |
//
//NC = input('numero di campioni:   ');
//while NC ~= [],
//   k = NC - 1;
//   n = 50; // punti aggiuntivi per il disegno
//   x = linspace(a,b,k+1)';
//   IdMat = eye(k+1,k+1);
//   Y = [];
//   for j = 1:k+1,
//       y = IdMat(:,j);
//       xp = linspace(a,b,k+1 + k*n)';
//       yp = PolyInterpNewt(x, y, xp);
//       Y(:,j) = abs(yp);
//   end;
//   s = sum(Y,'c')
//   clf();
//   plot(xp,s,'r','thickness',2);
//   xgrid();
//   xlabel('$x$', 'fontsize',4);
//   titolo = 'Grafico di sum_j=1:k | l(j) | con ' + string(NC)+ ' campioni.';
//   title(titolo,'fontsize',3);
//   mprintf('massimo (con %d campioni): %.4e',NC,max(s));
//   NC = input('numero di campioni:   ');
//end;

// *** Data Fun, studiare l'errore di ricostruzione con il polinomio interpolante.
//
function y = Fun(x)
    y = exp(-x) .* sin(3*x);
endfunction
mprintf('\n');
NC = input('numero di campioni:   ');
while NC ~= [],
  k = NC - 1;
  n = 50; // punti aggiuntivi per il disegno
  x = linspace(a,b,k+1)';
  y = Fun(x);
  xp = linspace(a,b,k+1 + k*n)';
  yp = PolyInterpNewt(x, y, xp);
  clf();
  plot(xp,yp, 'b','thickness',2);
  plot(xp,Fun(xp), 'r','thickness',2);
  plot(x,y, 'bkx','thickness',2);
  xgrid(); xlabel('$x$', 'fontsize',4);
  legend('p(x)','Fun(x)','(x(i),y(i))','in_upper_right');
  mprintf('err di ricostr (con %d campioni): %.4e\n\n',k+1,max(abs(yp - Fun(xp))));
  halt('premere invio');
  clf();
  plot(xp, abs(yp - Fun(xp)), 'r','thickness',2);
  title('| r(c(f)) - f |','fontsize',3);
  xgrid(); xlabel('$x$', 'fontsize',4);
  mprintf('\n');
  NC = input('numero di campioni:   ');
end;
