clear all; close all;
%-------------------------- demo_DRQ.m ---------------------------
% This script returns the DRQ error, the exact integral and its
% approssimation with Direct Rational Quadrature (DRQ method).
% 
% Function called: weights, gauss_jacobi_integration_xw, eval_bi,
%                  gauss_jacobi_integration_xw
% 
%-----------------------------------------------------------------
% Authors: 
%      Stefano De Marchi and Cinzia Bandiziol
%      Department of Pure and Applied Math.
%      University of Padova
%           
% Version: May 2015		
%------------------------------------------------------------------			

a=-5;
b=5;
n=21;
d=3;

%f=inline('abs(x)','x');
f=inline('1./(1+x.^(2))','x');
%f=inline('sin(x)','x');

% Use equidistant nodes
x=linspace(a,b,n)';
y=feval(f,x);

% Compute the F-H's weights
w_fh=weights(d,x);
[x_gauss w_gauss]=gauss_jacobi_integration_xw(500,0,0,a,b);
C=eval_bi(x,w_fh,x_gauss);
Int_vect=C*w_gauss;
Int_DRQ=Int_vect'*y;

% Compute the integral of f with Gauss-Legendre quadrature formula
[Int_exact]=gauss_jacobi_integration_I(1000,0,0,a,b,f);


Err=abs(Int_exact-Int_DRQ);
fprintf('Exact Integral:%f\nDRQ_Integral:%f\nDRQ_Error:%f\n',Int_exact,Int_DRQ,Err);