%
% A script to do Filon quadrature of int( f(t)*cos(omega*t), t=-1..1 )
%
% Hermite interpolation of f, differentiation matrices, and the
% method of undetermined polynomials are used:
% int( p*cos(omega*t), t ) = P*cos(omega*t) + Q*sin(omega*t)
% where Q'-omega*P=0 and P'+omega*Q=p.
%
% RMC December 30, 2010
%
f = @(t) 1.0./(1+t.^2);
df= @(t) -2.0*t./(1+t.^2).^2;
%f = @(t) 1 + 2*t + 3*t.^2 + 4*t.^3 + 5*t.^4;
%df= @(t) 2 + 6*t + 12*t.^2 + 20*t.^3;
tau = [-1,0,1];
%tau = cos( pi*(0:3)/3 );
%tau = [-1, -1/2, 1/2, 1];
[gam,D] = genbarywts(tau,2);
rho = [f(tau); df(tau)];
p = rho(:);
w = 200.3; % omega=2003, Maple = -0.485623502941671423e-3, quad fails
A = D^2 + w^2*eye(size(D)); % D^2
Q = A\(w*p);
P = D*Q/w;
% use the fact that tau(1)=-1 and tau(end)=1
Filon = (P(end-1)-P(1))*cos(w) + (Q(end-1)+Q(1))*sin(w)
% Now try quad, and see how it compares
osc = @(t) f(t).*cos(w*t);
qans = quad(osc,tau(1),tau(end),1.0e-9);
reler = Filon/qans-1