0

As my question states, I want to calculate the Fourier transform F(q) of a radial function f(r) (defined on [0,infinity[ and which decays like an exponential exp(-Ar +b) at large r) as accurately as possible in Fortran. The function values come from a data file (which I can easily interpolate through cubic interpolation for example and extrapolate since the behaviour at large r is known).

I'm using the "physics" definition of the Fourier transform in 3D, which gives (because f is radial) : enter image description here

I first tried to calculate this integral for some chosen values of q by using Gauss-Legendre quadrature, by generating some 60 or 100 abscissas and weights via the NAG routine D01BCF (D01BCF link). In the case of Gauss Legendre quadrature, the problem is to choose the interval [0,B] on which to integrate. While the function f loses 4 to 5 orders of magnitude from r=10 to r=20 (example), the choice of B as a strong influence on the result of the calculation... When I compared the result I get to a "nearly exact" calculation (made with matlab but with a veeeery long computation time), I saw that in fact this was only valid for small values of q (of the order of 5, when I have to deal with values as large as 150). A Gauss-Laguerre quadrature does not give any better result, probably because of the oscillatory part of the integrand.

I then tried to compute this Fourier transform for some given values of q with the routine D01ASF (D01ASF link). It is a "one-dimensional quadrature, adaptive, semi-infinite interval, weight function cos(ωx) or sin(ωx) ", which is exactly what I need. The results are quite convincing for q up to 80 or 100 if I input absolute error tolerances of 10E-5. Problems are : I would need to go at larger q, and the Fourier transform F(q) oscillates with a magnitude of ~ 10E-6 at such q's. Lowering the tolerance to 10E-5 already takes some time and even makes the whole thing to output some error message from the subroutine so I don't know if 10E-6 would be feasible.

I'm thus currently wondering if trying to calculate this Fourier transform with FFT wouldn't be a good idea ? The problems I face are that I don't know how to calculate radial wave functions with FFT (and also that I don't even know how to use FFT properly either since the definition of the transform is not even the same (exponent sign and argument) and that I never used it before).

Would you have ideas ? :)

EDIT 2 : I tried by FFT (using the routine C06FAF from NAG library). It works quite well up to some large values of q. The problem I face is that there is always some constant normalising factor to account for. I don't get why. This normalising factor evolves with the number N of points used in the mesh. It has the for of a power law : Normalising Factor F = N^(-0.5) x exp(9.9) approximately (see figure figure where the black line is the "exact" Fourier Transform and the green, magenta, blue, red and yellow lines are the FFT calculated for different values of N)

EDIT3 : I found the factor to be A*N^(-0.5) where A is the length of the integration mesh

mwoua
  • 403
  • 5
  • 13
  • one bit of advice is to test your quadrature schemes using an analytic `f(r)` for which you have an exact solution. ( `f(r)=Exp(-r)` for example ) – agentp Aug 24 '16 at 15:58
  • 3
    My idea is to post this question at http://scicomp.stackexchange.com and remove it from here. – High Performance Mark Aug 24 '16 at 16:21
  • @agentp I tried by FFT and I get some weird result when I compare the analytical FT(exp(-r)) to the computed one. There seems to be a normalising factor to take care of (look also at the EDIT) – mwoua Aug 30 '16 at 16:45
  • @agentp Also, at some point, the Fourier transform calculated by FFT (for exp(-r)) drops (at the end of the frequencies mesh). I tested it for some fermi dirac distribution (1/exp((r-R)/a)) and it also diverges at some point (and quite soon...) – mwoua Aug 30 '16 at 16:51
  • A scale factor of `n**(-0.5)` in the definition of the transform is mentioned in the documentation. Are you saying that you see a factor beyond that? – francescalus Aug 31 '16 at 10:38
  • I found the factor. It is actually A*N**(-0.5) where [0,A] is the interval of the calculation. – mwoua Aug 31 '16 at 12:49
  • The FFT algorithm calculates the discrete Fourier transform (DFT) which is a poor approximation to the continuous Fourier transform (FT). It calculates the integral with a rectangular rule or Rieman sum so it is not vey accurate. Quadpack from [link](http://nines.cs.kuleuven.be/software/QUADPACK/) has a quadrature subroutine for oscillating integrands (I think it's called QDAG). Unfortunately, speed and accuracy conflict with each other. – JAC Jan 08 '19 at 17:02

0 Answers0