home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #19 / NN_1992_19.iso / spool / sci / math / numanal / 2632 < prev    next >
Encoding:
Text File  |  1992-09-03  |  2.3 KB  |  68 lines
  1. Newsgroups: sci.math.num-analysis
  2. Path: sparky!uunet!mcsun!sunic!kth.se!nada.kth.se!mihai
  3. From: mihai@nada.kth.se (Mihai Dorobantu)
  4. Subject: Re: Oscillatory Integration?
  5. Message-ID: <1992Sep3.163359.6239@kth.se>
  6. Sender: usenet@kth.se (Usenet)
  7. Nntp-Posting-Host: cyklop.nada.kth.se
  8. Organization: Royal Institute of Technology, Stockholm, Sweden
  9. References:  <92247.091826J0S@psuvm.psu.edu>
  10. Date: Thu, 3 Sep 1992 16:33:59 GMT
  11. Lines: 55
  12.  
  13. In article <92247.091826J0S@psuvm.psu.edu>, Jack W. Sharer <J0S@psuvm.psu.edu> writes:
  14. |> I'm having to evaluate this integral:
  15. |> 
  16. |>           /\ 0.5
  17. |>          /           N
  18. |>         /      f(x) D (2 pi x + phi) dx
  19. |>        /
  20. |>      \/ -0.5
  21. |> 
  22. |>                        N
  23. |> for N ~ 500 and where D (x) is the N-point Dirichlet kernel:
  24. |> 
  25. |> 
  26. |>                iNx
  27. |>        N      e    - 1
  28. |>       D (x) = --------,  i = sqrt(-1)   (highly oscilatory for large N)
  29. |>                 ix
  30. |>                e  - 1
  31. |> 
  32. |> I'm successfully using a 2000-point Simpson integration, but I cannot
  33. |> begin to evaluate all the needed cases (limited by computer resources).
  34. |> 
  35. |> Could you suggest a change of variable or another technique which
  36. |> would dramatically improve the integration speed?
  37. |> 
  38. |> Thanks, Jack Sharer, Penn State, j0s@psuvm.psu.edu
  39. Try an FFT !
  40. Note that D(x) = 1 + exp(ix) + ...+ exp(i(N-1)x) = exp(-iMx)+...+exp(iMx)
  41. where M=N/2.
  42.  
  43. So the integral is just the sum of the first N Fourier coeficients of f viewd
  44. as a 1-periodic function. The only thing to do is to pick up N equaly spaced
  45. points on [-1/2 1/2], evaluate f and then add the components of fft(f).
  46. For a non-zero phase parameter phi, just change the variable but don't
  47. forget to extend f. See your implementation of FFT, you might have to
  48. multiply or divide by sqrt(N) depending on whatever normalization was done
  49. in implementing the FFT algorithm.
  50.  
  51. I tried it in MATLAB. Up to N=2^14 i had an "instant" answer, waited
  52. for 1 sec. with N=2^16 and a few seconds for N=2^17. (And I wasn't useing
  53. a workstation or a fast PC, just a terminal,sharing time with 40 
  54. other users!!!)
  55.  
  56. The FFT algorithm uses O(N log N) operations, so it is almost as 
  57. fast as possible. It is also very stable to rounding and truncation errors.
  58.  
  59.  
  60. Good luck,
  61. Mihai Dorobantu
  62.  
  63. Dept.for Numerical Analysis and Computer Science
  64. The Royal Institute for Technology, Stockholm - Sweden
  65.  
  66.  
  67.  
  68.