%Esempio di interpolazione instabile
% tratto da Quarteroni-Saleri; Esempio 3.3 pag. 85

% Sull'intervallo [-1,1] interpoliamo la funzione
% f(x)=sin(2pi x) su 22 nodi equispaziati.
% Successivamente generiamo un insieme di valori fi(nodi)
% ottenuti perturbando in maniera casuale i valori f(nodi)
% In figura vengono confrontati i due polinomi di interpolazione
% Agli estremi dell'intervallo di interpolazione la differenza
% è molto più grande della perturbazione operata


clear all  %pulisce la memoria
close all  %chiude tutte le finestre grafuche aperte
clc        %pulisce la schermata

x=linspace(-1,1,100);  %costruisce una suddivisione equispaziata
                       %dell'intervallo [-1,1]
                       %con 100 punti
                       
y=sin(2*pi*x);         %costruisce la funzione sul vettore x

n=22;                  %n contiene il numero nodi
nodi=linspace(-1,1,n); %costruisce le ascisse dei nodi
                       %in modo equispaziato
                       
ynodi=sin(2*pi*nodi);  %costruisce le ordinate dei nodi

plot(x,y,'b',nodi,ynodi,'ro');  %disegna la funzione -blu-
                                %e i nodi -rossi-
hold on                         %mantiene attiva la finestra grafica

c=polyfit(nodi,ynodi,n-1);  %costruisce il polinomio di 
                            %interpolazione su (nodi,ynodi)
                            %Il suo grado è n-1
                            %c contiene i coefficienti 
                            %del polinomio di interpolazione
                            
yi=polyval(c,x);        %valuta il polinomio di interpolazione 
                        %sulle 100 ascisse x
                        
plot(x,yi,'g')       %disegna il polinomio di interpolazione -verde-
hold on

ynodiperturbati=ynodi+0.0005*randn(1,n);
                     %perturba le ordinate dei nodi mediante
                     %un generatore di numeri casuali randn

plot(nodi,ynodiperturbati,'k+')  %disegna i nodi pertirbati

cpert=polyfit(nodi,ynodiperturbati,n-1);  %costruisce il polinomio di 
                            %interpolazione su (nodi,ynodiperturbati)
                            %cpert contiene i coefficienti 
                            %del polinomio di interpolazione
                            
yip=polyval(cpert,x);   %valuta il polinomio di interpolazione 
                        %sulle 100 ascisse x
                        
plot(x,yip,'c')       %disegna il polinomio di interpolazione -azzurro-
grid on               %disegna la griglia nel grafico
legend('esatta','nodi','interpolata esatta','nodi perturbati','interpolata perturbata')
                      %legenda
hold off              %disattiva la finestra grafica

%Valutazione di
%1- scarto massimo tra dati esatti (ynodi) dati perturbati (ynodiperturbati)
%2- costante di Lebesgue (numero di condizionamento del problema dell'interpolazione)
%              cost_leb=2^(n+1)/(e n (log(n)+gamma))
%   con n=numero nodi, e=numero di Nepero, gamma=costante di Eulero
%3- scarto massimo tra il polinomio di interpolazione costruito a partire
%   dai dati esatti e quello costruito a partire dai dati perturbati

gamma=0.547721;       %costante di Eulero
cost_leb=2^(n+1)/(exp(1)*n*(log(n)+gamma));
scartodati=norm(ynodi-ynodiperturbati,'inf');
scartointerp=norm(yip-yi,'inf');
disp('Condizionamento:')
disp(cost_leb)
disp('Distanza tra dati reali e dati perturbati:')
disp(scartodati)
disp('Scarto tra interpolazione esatta e interpolazione perturbata:')
disp(scartointerp)