% calcolare il polinomio di
% interpolazione nei nodi
% xn di una funzione data

% definisco la funzione
f = inline('1./(1+x.^2)');

% disegno la funzione 
% nell'intervallo [a, b]
a = -5;
b =  5;
x = linspace(a,b,500);
y = f(x);
plot(x,y,'k','lineWidth',2)
xlabel('x','fontSize',18)
ylabel('y','fontSize',18)

% specifico i nodi dove
% interpolare la funzione
numNodi = 4;
xnEq = linspace(a,b,numNodi);
ynEq = f(xnEq);

% aggiungo i nodi al grafico
% precedente
hold on
plot(xnEq,ynEq,'or','lineWidth',3)

% calcolo il polinomio di
% interpolazione secondo Lagrange
yPeq = lagpoly(xnEq,ynEq,x);

% aggiungo al grafico precedente
plot(x,yPeq,'b','lineWidth',2)
hold off

% aggiungo legenda al grafico
legend('f(x)','nodi','P(x) equispaziati')

%=========================================


% disegno la funzione 
% nell'intervallo [a, b]
figure(2)
plot(x,y,'k','lineWidth',2)
xlabel('x','fontSize',18)
ylabel('y','fontSize',18)

% specifico i nodi di Chebyshev
% dove interpolare la funzione
%            2k+1     pi
% tk = cos(--------* ---),  k=0,....n
%            n+1      2
%
% xn = (a+b)/2 +tk*(b-a)/2, k=0,...,n
%
k    = 0:numNodi;
phi  = ( (2*k+1)*pi )/(2*numNodi +2);
tk   = cos(phi);
xnCh = (a+b)/2+(b-a)/2*cos(phi);
ynCh = f(xnCh);

% aggiungo i nodi al grafico
% precedente
hold on
plot(xnCh,ynCh,'or','lineWidth',3)

% calcolo il polinomio di
% interpolazione secondo Lagrange
yPCh = lagpoly(xnCh,ynCh,x);

% aggiungo al grafico precedente
plot(x,yPCh,'b','lineWidth',2)
hold off

% aggiungo legenda al grafico
legend('f(x)','nodi','P(x) Chebyshev')


==============================================================

function y = lagpoly(xn,yn,x)
%LAGPOLY polinomio interpolante secondo Lagrange
%   Y = LAGPOLY(XN,YN,X) calcola le ordinate
%   del polinomio di Lagrange per le ascisse
%   X. I nodi sono contenuti in XN e i valori
%   nodali in YN.
%   L'uscita e' il vettore Y che ha le stesse
%   dimensioni di X.

y = zeros( size(x) );
for k=1:length(xn)
    % aggiungo il contributo del k-esimo
    % polinomio di Lagrange
    y = y + yn(k)*lk(xn,x,k);
end

=============================================================

function y = lk(xn,x,k)
%LK k-esimo polinomio di Lagrange
%   Y = LK(X) restituisce le ordinate
%   del k-esimo polinomio di Lagrange
%   valutate in corrispondenza del 
%   vettore delle ascisse X.
%            (x-xn(1))*(x-xn(2))*...
%   LK(x) = ---------------------------------
%            (xn(k)-xn(1))*(xn(k)-xn(2))*...
%
%

y = ones( size(x) );
for i=1:k-1
    % aggiungo contributo nodo i-esimo
    y = y.*( x-xn(i) )/( xn(k) - xn(i) );
end

for i=k+1:length(xn)
    % aggiungo contributo nodo i-esimo
    y = y.*( x-xn(i) )/( xn(k) - xn(i) );    
end

=============================================================

ESERCIZIO PROPOSTO

function nt = somma(s,nmax)
% SOMMA da' numero termini del tipo 1/(n(n+1))
% necessari per superare soglia s
%

somma = 0;
n = 1;
while(somma<s && n<=nmax)
    somma = somma + 1/(n*(n+1));
    n = n+1;
end
nt = n-1;



















