%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lagrange.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mostra il grafico dei polinomi di Lagrange

n = 6;
xnodes = linspace(0,1,n);
Y = eye(n); % la riga i descrive la delta di Dirac D_ij sui nodi xnodes;

for i=1:n
  pol_coeff(i,:) = polyfit(xnodes,Y(i,:),n-1); % coefficienti del polinomio
end;

% N.B. ciascuna riga riporta i coefficienti del polinomio interpolante
% nell ordine: dal coefficiente del termine di grado massimo
% al coefficiente del termine di grado 0

x = linspace(0,1,100); % nodi di valutazione dei grafici

for i=1:n
  figure(i)
  pol_values = polyval(pol_coeff(i,:),x);
  plot(x,pol_values,'b-',xnodes,Y(i,:),'ro')
  grid on
%  axis([0 1 0 1]);
  title(sprintf('l_%d (x)',i));
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% func.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val = func(x)

val = abs(x);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dfuncdx.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val = func(x)

n = length(x);

for i=1:n
  if x(i) < 0
     val(i) = -1;
  else
     val(i) = 1;
  end;
end;

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% d2funcdx2.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val = func(x)

val = 0;

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spline_der_conv.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mostra la convergenza delle derivate dell'interpolata spline
% nei punti in cui la funzione interpolata e` sufficientemente regolare.

xnodes = [-2 -1 0 1 2];

y = feval('func',xnodes);

x = linspace(-2, 2, 100);

spline_values = spline(xnodes,y,x);

poly_coeff = polyfit(xnodes,y,4);
poly_values = polyval(poly_coeff,x);

figure(1)
plot(x,spline_values,'b-',x,poly_values,'r-',xnodes,y,'g-o')

disp('Strike any key to continue')
pause

% Ora confrontiamo le derivate sull'intervallo [1 2]

x = linspace(1, 2, 100);

spline_piecewise_polynomial = spline(xnodes,y);
[xbreak, coeff, piece, ord] = unmkpp(spline_piecewise_polynomial);

% xbreak sono i punti xnodes
% coeff sono i coefficienti dei tratti polinomiali
% spl_i(x) = coeff(i,1)*(x-x_i)^3 + coeff(i,2)*(x-x_i)^2 +
%            coeff(i,3)*(x-x_i)   + coeff(i,4)

dspldx_coeff = [3*coeff(:,1) 2*coeff(:,2) coeff(:,3)];
d2spldx2_coeff = [6*coeff(:,1) 2*coeff(:,2)];

dspldx_pp = mkpp(xnodes,dspldx_coeff);
d2spldx2_pp = mkpp(xnodes,d2spldx2_coeff);

dspldx_values = ppval(dspldx_pp,x);
d2spldx2_values = ppval(d2spldx2_pp,x);


poly_coeff = polyfit(xnodes,y,4);
dpolydx = [4*poly_coeff(1) 3*poly_coeff(2) 2*poly_coeff(3) poly_coeff(4)]; 
d2polydx2 = [12*poly_coeff(1) 6*poly_coeff(2) 2*poly_coeff(3)]; 

dpolydx_values = polyval(dpolydx,x);
d2polydx2_values = polyval(d2polydx2,x);

spline_values = spline(xnodes,y,x);

poly_coeff = polyfit(xnodes,y,4);
poly_values = polyval(poly_coeff,x);

yv = feval('func',x);
dyvdx = feval('dfuncdx',x);
d2yvdx2 = feval('d2funcdx2',x);

figure(1)
plot(x,spline_values,'b-',x,poly_values,'r-',x,yv,'g-.')

figure(2)
plot(x,dspldx_values,'b-',x,dpolydx_values,'r-',x,dyvdx,'g-.')

figure(3)
plot(x,d2spldx2_values,'b-',x,d2polydx2_values,'r-',x,d2yvdx2,'g-.')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Strike any key to continue')
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = 15;

xnodes = linspace(-2, 2, N);

y = feval('func',xnodes);

x = linspace(-2, 2, 100);

spline_values = spline(xnodes,y,x);

poly_coeff = polyfit(xnodes,y,N-1);
poly_values = polyval(poly_coeff,x);

figure(4)
plot(x,spline_values,'b-',x,poly_values,'r-',xnodes,y,'g-o')

disp('Strike any key to continue')
pause

% Ora confrontiamo le derivate sull'intervallo [1 2]

x = linspace(1, 2, 100);

spline_piecewise_polynomial = spline(xnodes,y);
[xbreak, coeff, piece, ord] = unmkpp(spline_piecewise_polynomial);

% xbreak sono i punti xnodes
% coeff sono i coefficienti dei tratti polinomiali
% spl_i(x) = coeff(i,1)*(x-x_i)^3 + coeff(i,2)*(x-x_i)^2 +
%            coeff(i,3)*(x-x_i)   + coeff(i,4)

dspldx_coeff = [3*coeff(:,1) 2*coeff(:,2) coeff(:,3)];
d2spldx2_coeff = [6*coeff(:,1) 2*coeff(:,2)];

dspldx_pp = mkpp(xnodes,dspldx_coeff);
d2spldx2_pp = mkpp(xnodes,d2spldx2_coeff);

dspldx_values = ppval(dspldx_pp,x);
d2spldx2_values = ppval(d2spldx2_pp,x);


poly_coeff = polyfit(xnodes,y,N-1);
for i=1:N-1
    dpolydx(i)   = (N-i)*poly_coeff(1);
end;
for i=1:N-2
    d2polydx2(i) = (N-i)*(N-i-1)*poly_coeff(i); 
end;

dpolydx_values = polyval(dpolydx,x);
d2polydx2_values = polyval(d2polydx2,x);

spline_values = spline(xnodes,y,x);

poly_coeff = polyfit(xnodes,y,N-1);
poly_values = polyval(poly_coeff,x);

yv = feval('func',x);
dyvdx = feval('dfuncdx',x);
d2yvdx2 = feval('d2funcdx2',x);

figure(5)
plot(x,spline_values,'b-',x,poly_values,'r-',x,yv,'g-.')

figure(6)
plot(x,dspldx_values,'b-',x,dpolydx_values,'r-',x,dyvdx,'g-.')

figure(7)
plot(x,d2spldx2_values,'b-',x,d2polydx2_values,'r-',x,d2yvdx2,'g-.')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% difdiv.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dd = difdiv(x,y)

% dd: vettore delle differenze divise
% x:  vettore dei nodi di interpolazione
% y:  vettore dei valori della funzione da interpolare

n = length(x);

if length(y) ~= n
  disp('Errore, il vettore dei nodi e il vettore dei valori della funzione')
  disp('non hanno la stessa dimensione!')
  return
end;

for i = 1:n-1
  for j = n:-1:i+1
y(j) = (y(j)-y(j-1))/(x(j)-x(j-i));
  end;
end;

dd = y;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% interp.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function p = interp(x,dd,t)

% p:  valore del polinomio in t (puo` essere un vettore)
% dd: dettore delle differenze divise (Output di difdiv)
% x:  vettore dei nodi di interpolazione
% t:  valore dei nodi in cui valutare il polinomio (puo` essere un vettore)

     n = length(x);

if length(dd) ~= n
     disp('Errore, il vettore dei nodi e il vettore delle differenze divise')
     disp('non hanno la stessa dimensione!')
     return
end;

p = dd(n)*ones(size(t));
for i=n-1:-1:1
  p = p.*(t-x(i))+dd(i);
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% interpola.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [] = interpola(N)

% Questa funzione mostra come usare le funzioni difdiv e interp
% applicate alla funzione di Runge
% Grado dei polinomi

x = linspace(-5,5,N+1);

y = 1./(1+x.^2);

dd = difdiv(x,y);

xx = linspace(-5,5);
yy = 1./(1+xx.^2);
pp = interp(x,dd,xx);

plot(xx,yy,'r-',xx,pp,'g-')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% trapezi.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%*****************************************************************
% File: TRAPEZI.M
%
% Uso:	t = trapezi(a,b,n,f)
%
% Scopo: Calcolo integrali con formula composta dei trapezi
%
% Input:  a,b	intervallo di integrazionee
%	   n	numero punti
%	   f	funzione integranda (stringa, es 'x.^2.*cos(x)')
%
% Output:  t	approssimazione dell' integrale
%*****************************************************************

function  t=trapezi(a,b,n,f)
h=(b-a)/(n-1);
x=[a:h:b];
% per certi valori di n, h non fa arrivare a b
if size(x)~=n;x(n)=b;end
%
%
y=eval(f);
t=2*sum(y(2:n-1));
t=.5*h*(t+y(1)+y(n));
%
return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% simpson.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%*****************************************************************
% File: SIMPSON.M
%
% Uso:	t = simpson(a,b,n,f)
%
% Scopo: Calcolo integrali con con formula composta di Simpson
%
% Input:  a,b	intervallo di integrazionee
%	   n	numero punti (dispari)
%	   f	funzione integranda (stringa, es 'x.^2.*cos(x)')
%
% Output:  t	approssimazione dell' integrale
%*****************************************************************

function  t=simpson(a,b,n,f)
if n == fix(n/2)*2
   t=' ** ERRORE ** numero pari di punti';
   return
end
%
h=(b-a)/(n-1);
x=[a:h:b];
% per certi valori di n, h non fa arrivare a b
if max(size(x))~=n;x(n)=b;end
%
y=eval(f);
t=4*sum(y(2:2:n));
t=t+2*sum(y(3:2:n-1));
t=h*(t+y(1)+y(n))/3;
%
return