matlab

matlab script.m

% POLNEWTON Calcola i coefficienti del polinomio interpolatore
% utilizzando la forma di Newton con le differenze
% divise
%
% Uso:
% c = polnewton (x,y)
%
% Dati di ingresso:
% x vettore dei nodi
% y vettore dei valori della funzione da interpolare nei nodi
%
% Dati di uscita:
% c vettore colonna dei coefficienti ordinati per indici
% crescenti (c_0, c_1, ... )   
function [c] = polnewton(x,y)

n1=length(x);
if n1~= length(y)
   error('tab errata') 
end
    
    ctab=zeros(n1,n1);

    ctab(:,1)=y;
    %inizializza la prima colonna
    
  for j=2:n1
      for i=1:n1-j+1
          
          ctab(i,j)=(ctab(i,j-1)-ctab(i+1,j-1))/(x(i)-x(i+j-1));
      
      end
  end
c=ctab(1,:);
c=c(:);
end

%-------------------------------------------------------------

% HORNER Calcola il valore del polinomio interpolatore in x^*
% utilizzando la forma di Newton e l'algoritmo di Horner
%
% Uso:
% fxstar = horner (x,c,xstar)
%
% Dati di ingresso:
% x vettore dei nodi
% c vettore dei coefficienti della forma di Newton
% ordinati per indici crescenti (c_0, c_1, ... )
% xstar valore in cui si vuole valutare il polinomio
%
% Dati di uscita:
% fxstar valore di P(x^*)
function fxstar = horner (x,c,xstar)

n1=length(x);

if(n1~=length(c))
    error('error')
end

fxstar=c(n1);

 for i=n1-1:-1:1
   fxstar=fxstar*(xstar-x(i))+c(i);
 end


end
%---------------------------------------------------------------------


%Si scriva poi uno script che calcoli il valore del polinomio di interpolazione in 201 valori dell'intervallo
%I, per poterlo rappresentare graficamente (linea tratteggiata rossa).
%Sullo stesso grafico si rappresentino anche i punti della tabulazione (con un cerchietto rosso) e la funzione
%data (linea intera nera).


claer all

f=inline('1./(1+x.^2)');
a=-5; b=5;
nnodi=11;

xn=linspace(a,b,nnodi);
fxn=feval(f,xn);

c=polnewton(xn fxn);

xint= linspace(a,b,201);
xint=xint(:);

y=feval(f,xint);

for i=1:201
    fxstar(i)=horner (xn,c,xint(i));
end

fxstar=fxstar(:);
figure(1)
plot(xn,fxn,'ro')
hold on
plot(xint,fxstar,'r--')
hold off

legend('punti dati','Funzione','Polinomio','Location','Best')
hold off

%----------------------------------------------------------

%BISEZFUN Metodo di Bisezione 
% 
% Uso: 
% [xv, fxv, n] = bisezione(f, a, b, toll, nmax) 
% % Dati di ingresso: 
% f: funzione (inline function) 
% a: estremo sinistro 
% b: estremo destro 
% toll: tolleranza richiesta per l’ampiezza 
% dell’intervallo 
% nmax: massimo indice dell’iterata permesso 
% % Dati di uscita: 
% xv: vettore contenente le iterate 
% fxv: vettore contenente i corrispondenti residui 
% n: indice dell’iterata finale calcolata

function [xv, fxv, n] = Bisezione (f, a, b, toll, nmax)
    n = -1;
    fa = f(a);
    amp = toll + 1;
    xv = [];
    fxv = [];
    while (amp >= toll && n < nmax)
        n = n + 1;
        amp = abs(b-a);
        x = a + amp/2;
        fx = f(x);
        xv = [xv : x];
        fxv = [fxv : fx];
        if (fa * fx < 0)
            b = x;
        else if (fa * fx > 0)
                a = x;
                fa = fx;
            else
                amp = 0;
            end
        end
    end
end

%-----------------------------------------
%NEWTONFUN Metodo di Newton
% Uso:
% [xv, fxv, n, flag] = newtonfun (f, f1, x0, toll, nmax)
%
% Dati di ingresso:
% f: funzione
% f1: derivata prima
% x0: valore iniziale
% toll: tolleranza richiesta per il valore assoluto
% della differenza di due iterate successive
% nmax: massimo numero di iterazioni permesse
%
% Dati di uscita:
% xv: vettore contenente le iterate
% fxv: vettore contenente i corrispondenti residui
% n: numero di iterazioni effettuate
% flag: Se flag = 1 la derivata prima si e' annullata

function [xv, fxv, n, flag] = Newton (f, f1, x0, toll, nmax)
    n = 0;
    x = x0;
    diff = toll;
    flag = 0;
    xv = [];
    fxv = [];
    while (abs(diff) >= toll && n < nmax)
        n = n + 1;
        fx = f(x);
        f1x = f1(x);
        if (f1x == 0)
            diff = 0;
            flag = 1;
        else
            diff = -fx/f1x;
            x = x + diff;
            xv = [xv , x];
            fxv = [fxv , fx];
        end
    end
end

%---------------------------------------------------------------------------------------------

% Script minimiquadrati

clear all;

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]';%Vettore Xn Trasposto colonna!
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.184]';%Vettore f(Xn) Trasposto!

m1=length(x);
A = [ones(m1,1) x x.^2 x.^3 x.^4 ];%polinomio al piu di quarto grado
 % A' trasposta di A
 
 M=A'*A; 
 z=A'*y;
 a=M\z;
 
 xs= linspace(min(x),max(x),200);
 xs=xs';
 xz=polyval(flipud(a),xs);
 
 xbari=sum(x)/m1;
 ybari=sum(y)/m1;
 
 
 plot([min(x),max(x)],[0,0],'k',x,y,'b*',xs,xz,'r-',xbari,ybari,'go');
 
 
 %------------------------------------------------------------------------------------------