%====================================================================
% Script per determinare il polinomio al piu’ di quarto grado
% approssimante ai minimi quadrati. Pesi unitari.
% Necessita delle function triinf e trisup
%====================================================================
clear all
% Definizione dati e controllo
x = [-3.490 -2.948 -2.574 -2.157 -1.377 -1.234 -0.861 -0.116 ...
0.235 0.558 1.036 1.318 2.139 2.566 2.736 3.312]';
y = [27.200 4.720 -0.978 4.100 16.013 19.656 22.498 21.650 ...
16.770 12.671 4.042 -2.158 -16.901 -11.437 -13.449 31.1840]';
n = 4;
m1 = length(x);
if length(y) ~= m1
error('Le dimensioni di x e di y non sono compatibili')
end
% Creazione della matrice A
A = [ones(m1,1) x x.^2 x.^3 x.^4];
% oppure
% A = fliplr(vander(x));
% A = A(:,1:n+1);
% Creazione del vettore b
b = y;

% Calcolo coefficienti con risoluzione sistema delle equazioni normali
% I coefficienti vengono calcolati rispetto alla
% base canonica, ordinati per potenze crescenti della x.
M=A'*A;
z=A'*b;
[L,flag] = chol(M,'lower');
yy=L\z;
a = L'\yy;
% Calcola il baricentro
disp('Baricentro')
xbari=sum(x)/m1;
ybari=sum(y)/m1;
disp([xbari,ybari])
% Grafico
minx=min(x);
maxx=max(x);
xs = linspace(minx, maxx, 200);
xs = xs(:);
ys = polyval(flipud(a), xs);
plot([minx maxx], [0 0], 'k', x, y, 'b*', xs, ys, 'r', xbari,ybari,'go')
title(['Approssimazione ai minimi quadrati con polinomio di grado ',num2str(n)])