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)
    //
    // Le componenti della colonna yp sono i valori
    // in xp(1),...,xp($) del polinomio che interpola
    // i dati (x(1),y(1)),...,(x($),y($)).
    //
    k = length(x) - 1;
    N = zeros(k+1,k+1); // matrice del sist di interpolazione
    N(:,1) = ones(k+1,1);
    N(2:k+1,2) = x(2:k+1) - x(1);
    for j = 3:k+1,
        N(j:k+1,j) = N(j:k+1,j-1) .* (x(j:k+1) - x(j-1));
    end;
    // coefficienti pol interp: sol di Nx = y
    b = SA(N,y);
    // calcolo valori del pol interp con Horner
    yp = b(k+1);
    for j = k:-1:1,
        yp = b(j) + (xp - x(j)) .* yp;
    end;
endfunction
//
// Procedura PolyMQ ***************************************
//
function yp = PolyMQ(x,y,n,xp)
  //
  // Le componenti della colonna yp sono i valori
  // in xp(1),...,xp($) del polinomio di grado <= n
  // che meglio approssima i dati (x(1),y(1)),...,(x($),y($))
  // nel senso dei minimi quadrati.
  //
  k = length(x) - 1;
  N = zeros(k+1,n+1); // matrice del sist di interpolazione
  N(:,1) = ones(k+1,1);
  N(2:k+1,2) = x(2:k+1) - x(1);
  if n > 1 then
     for j = 3:n+1,
         N(j:k+1,j) = N(j:k+1,j-1) .* (x(j:k+1) - x(j-1));
     end;
  end;
  // coefficienti pol che meglio approssima i dati: sol mq di Nx = y
  b = N\y;
  // calcolo valori del pol approssimante con Horner
  yp = b(n+1);
    for j = n:-1:1,
        yp = b(j) + (xp - x(j)) .* yp;
    end;
endfunction
//
//  Dati
//
function y = f(x)
    y = 1 - x .^2;
//    y = sin(%pi * x / 10);
endfunction
//
E = 1e-2;
N = 15;
x = linspace(0,10,N)';
y = f(x) + E*rand(x);
//
// confronto interpolazione / mq
//
NN = N + 50 * (N-1);
xp = linspace(0,10,NN)'; vp = PolyMQ(x,y,4,xp); vpi = PolyInterpNewt(x,y,xp);
err_MQ = max(abs(f(xp) - vp));
err_I = max(abs(f(xp) - vpi));
printf('\nerr_MQ = %5.4e   ,   err_I = %5.4e   ,   err_I / err_MQ = %5.4e\n',err_MQ,err_I,err_I/err_MQ);
clf(); 
subplot(211); plot2d(xp,vpi,style = 2); plot2d(x,y,style = -2); xgrid(); 
legend('vpi','dati')
subplot(212); plot2d(xp,vp,style = 5); plot2d(x,y,style = -2); xgrid(); xlabel('x');
legend('vp','dati')
