%*****************************************************************
% 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


%*****************************************************************
% 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


%*****************************************************************
% File:  QUADCHEB.M
%
% Scopo: Quadratura di Chebyshev in [-1,1]
%
% Uso:	inte = quadcheb(f,n)
%
% Input:   f	 funzione integranda  (stringa, es 'x.^2.*cos(x)')
%	   n	 numero nodi per quadratura.
%
% Output:  inte  integrale calcolato.
%
%*****************************************************************

function inte=quadcheb(f,n)
j=[0:n-1];
pi2n=pi/(2*n);
x=cos((2*j+1)*pi2n);
yf=eval(f);
inte=pi*sum(yf)/n;


%*****************************************************************
% File:  ZPLEGE.M
%
% Scopo:  Calcolo degli zeri e pesi di Legendre in [-1,1].
%
% Uso:	 [z,p] = zplege(n)
%
% Input:   n   grado del polinomio
%
% Output:  z   nodi della formula di quadratura
%	   p   pesi della formula di quadratura
%
% Richiami:
% Commenti:   p(k+1) = x*p(k) - k^2/(4*k-1)*p(k-1)
% Autori: mf
%*****************************************************************

function [z,p] = zplege(n)
if n<=1
  z=0;p=2;
  return
end
jac=zeros(n);
k=[1:n-1];
v=(k)./(sqrt(4*(k.^2)-1));
jac=jac+diag(v,1)+diag(v,-1);
[p,z]=eig(jac);
norm2=sqrt(diag(p'*p));    % normalizzazione per pesi
p=(2*p(1,:)'.^2)./norm2;   % 2=mu0; ' per aver p colonna
z=diag(z);		   % z colonna


%*****************************************************************
% File:  QUADLEGE.M
%
% Scopo: Quadratura di Legendre in [-1,1]
%
% Uso:	inte = quadlege(f,n)
%
% Input:    n	numero punti (ordine formula)
%	    f	funzione da integrare  (stringa, es 'x.^2.*cos(x)')
%
% Output:   t	integrale calcolato.
%
% Richiami: ZPLEGE
%*****************************************************************

function  inte = quadlege (f,n)
%
%  Calcola zeri e pesi
%
[x,p]=zplege(n);
%
y=eval(f);
inte=p'*y;
%
return


%*****************************************************************
% File:  QUADLEGETOLL.M
%
% Scopo: Quadratura di Legendre in [-1,1]
%
% Uso:	inte = quadlege(f,n)
%
% Input:    f     funzione da integrare  (stringa, es 'x.^2.*cos(x)')
%           toll  tolleranza richiesta
%           N_max numero massimo di nodi
%
% Output:   t     integrale calcolato.
%
% Richiami: ZPLEGE
%*****************************************************************

function  [inte2,m] = quadlegetoll(f,toll,N_max)
%
%  Calcola zeri e pesi
%
err = inf;
n   = 5;
while(n<=N_max & err > toll)
[x,p]=zplege(n);
y=eval(f);
inte1=p'*y;
%
m = n+2;
[x,p]=zplege(m);
y=eval(f);
inte2=p'*y;
%
err = abs(inte2-inte1);
n = n+4;
%
end

if n>=N_max
disp('Be careull! Maximal number of nodes reached')
end

return