% esperimenti con i polinomi ortogonali

% grafico polinomi ortogonali monici
N = 201;    % n. punti
n = 10;     % n. polinomi
a = -1;
b = 1;
t = linspace(a,b,N)';
vk = [1:n]';

% coefficienti dei polinomi di Legendre w(x)=1
alfa = zeros(n+1,1);
beta2 = [2;(vk.^2)./(4*vk.^2-1)];
beta = sqrt(beta2);

% coefficienti dei polinomi di Chebychev w(x)=1/sqrt(1-x^2)
% alfa = zeros(n+1,1);
% beta2 = [pi/2;0.25*ones(n,1)];
% beta = sqrt(beta2);

% grafico dei polinomi monici
P1 = ones(N,n+1);
P1(:,2) = (t-alfa(1));
for k = 1:n-1
    P1(:,k+2) = (t-alfa(k+1)).*P1(:,k+1)-beta2(k+1)*P1(:,k);
end
figure(1)
plot(t,P1)
title('Polinomi monici')

% grafico dei polinomi ortonormali
P2 = ones(N,n+1);
P2(:,1) = 1/beta(1);
P2(:,2) = (t-alfa(1))/beta(2);
for k = 1:n-1
    P2(:,k+2) = (t-alfa(k+1)).*P2(:,k+1)-beta(k+1)*P2(:,k);
    P2(:,k+2) = P2(:,k+2)/beta(k+2);
end
figure(2)
plot(t,P2)
title('Polinomi ortonormali')

% formula di quadratura
J = diag(alfa)+diag(beta(2:n+1),-1)+diag(beta(2:n+1),1);
[V,D] = eig(J);		% L'ERRORE ERA QUI!!!
nodi = diag(D);             % colonna
pesi = beta2(1)*V(1,:).^2;  % riga
effe = @(x) exp(x);
fquad = pesi*effe(nodi)
integralesatto = exp(1)-1/exp(1)
errore = abs(integralesatto-fquad)
