Filon quadrature (oscillatory integrands) in maxima.
Subject: Filon quadrature (oscillatory integrands) in maxima.
From: Richard Fateman
Date: Tue, 24 Apr 2007 21:51:31 -0700
In case someone is interested in numerical integration of functions of the
form f(x)*sin(x), Filon quadrature does a good job. (Especially if it is
like f(x)*sin(1000*x) if you get the drift.)
Why not just use quadpack? Because the function below can use bigfloats.
I attach only the shortest version, from a paper in progress...
filons0(f,m,p):= /*integrate f(x)*sin(m pi x) from 0 to 1, using p panels
*/
block([h:bfloat(1/(2*p),
k:m*%pi, theta,s2p,s2pm1,alf,bet,gam],
theta:bfloat(k*h),
s2p:
block([sum:0],
for i:0 thru p do sum:sum+f(2*h*i)*sin(2*k*h*i),
sum-1/2*(f(1)*sin(k))),
s2pm1:
block([sum:0],
for i:1 thru p do sum:sum+f(h*(2*i-1))*sin (k*h*(2*i-1)), sum),
alf: 1/theta + sin(2*theta)/2/theta^2-2*sin(theta)^2/theta^3,
bet: 2*((1+cos(theta)^2)/theta^2 -sin(2*theta)/theta^3),
gam: 4*(sin(theta)/theta^3-cos(theta)/theta^2),
h*(alf*(f(0)-f(1))*cos(k) +bet*s2p+ gam*s2pm1))$
Reference: CACM algorithm 353, by Chase and Fosdick. Find it with google..