🚀 go-pugleaf

In article <92247.091826J0S@psuvm.psu.edu>, Jack W. Sharer <J0S@psuvm.psu.edu> writes:
|> I'm having to evaluate this integral:
|> 
|>           /\ 0.5
|>          /           N
|>         /      f(x) D (2 pi x + phi) dx
|>        /
|>      \/ -0.5
|> 
|>                        N
|> for N ~ 500 and where D (x) is the N-point Dirichlet kernel:
|> 
|> 
|>                iNx
|>        N      e    - 1
|>       D (x) = --------,  i = sqrt(-1)   (highly oscilatory for large N)
|>                 ix
|>                e  - 1
|> 
|> I'm successfully using a 2000-point Simpson integration, but I cannot
|> begin to evaluate all the needed cases (limited by computer resources).
|> 
|> Could you suggest a change of variable or another technique which
|> would dramatically improve the integration speed?
|> 
|> Thanks, Jack Sharer, Penn State, j0s@psuvm.psu.edu
Try an FFT !
Note that D(x) = 1 + exp(ix) + ...+ exp(i(N-1)x) = exp(-iMx)+...+exp(iMx)
where M=N/2.

So the integral is just the sum of the first N Fourier coeficients of f viewd
as a 1-periodic function. The only thing to do is to pick up N equaly spaced
points on [-1/2 1/2], evaluate f and then add the components of fft(f).
For a non-zero phase parameter phi, just change the variable but don't
forget to extend f. See your implementation of FFT, you might have to
multiply or divide by sqrt(N) depending on whatever normalization was done
in implementing the FFT algorithm.

I tried it in MATLAB. Up to N=2^14 i had an "instant" answer, waited
for 1 sec. with N=2^16 and a few seconds for N=2^17. (And I wasn't useing
a workstation or a fast PC, just a terminal,sharing time with 40 
other users!!!)

The FFT algorithm uses O(N log N) operations, so it is almost as 
fast as possible. It is also very stable to rounding and truncation errors.


Good luck,
Mihai Dorobantu

Dept.for Numerical Analysis and Computer Science
The Royal Institute for Technology, Stockholm - Sweden