1

I want to numerically compute following integral (as fast as possible):

\begin{equation} f(x) = \left| \int e^{ \mathrm{i} p x / 2} f(p) \ \mathrm{d}p\right|^2 , \end{equation}

where $f(p)$ is given as a table on a discrete grid in p.

Example function:

f[p_] := Exp[-p^2] Exp[I p] Cos[p]^2 p^2

Generate the table on a discrete grid in p:

fp = Transpose[{Range[-5, 5, 0.01], 
   f[p] /. p -> Range[-5, 5, 0.01]}]

How I do it:

fx[fp_] := 
  ParallelTable[{x, ((fp[[All, 1]][[-1]] - fp[[All, 1]][[1]])/
        Length[fp[[All, 1]]])^2 Abs[
       Sum[fp[[All, 2]][[n]] Exp[I fp[[All, 1]][[n]] x/2], {n, 1, 
         Length[fp[[All, 1]]]}]]^2}, {x, -20, 20, 0.1}];

fx[fp] works fine, but takes 16 seconds on my machine.

Is there a way to do this much faster?

Edit: I have an answer code to my own question. I think it might be helpful for others.

Luke
  • 838
  • 4
  • 11
  • 2
    Fourier performs a fast Fourier transform, perhaps that's what you are looking for. – Ulrich Neumann Jun 22 '20 at 11:38
  • I can give Fourier only a list of values and it doesn't know about the grid, or am I missing something? – Luke Jun 22 '20 at 11:46
  • 2
    It is good so that it does not know about the grid. That is why it is so fast. – yarchik Jun 22 '20 at 11:56
  • 1
    As long as the grid is uniform and the wave length is a multiple of the interval's length, Fourier should work. Of course, you have to adjust some constants. – Henrik Schumacher Jun 22 '20 at 12:05
  • Could someone please write down an example code how I can map the result from 'Fourier' back to a proper grid? The results I get with 'Fourier' do not make much sense for me. I get large values at the boundaries and extremely small values in between. – Luke Jun 22 '20 at 12:33
  • 1
    This has already been asked many times here https://mathematica.stackexchange.com/questions/85139/what-do-the-x-and-y-axis-stand-for-in-the-fourier-transform-domain/85167#85167 – yarchik Jun 22 '20 at 14:20
  • Thank you for the link. I tried this approach but it seems to me that for Fourier there is a fundamental trade off between Range of my grid and the resolution I can get in Fourierspace. . – Luke Jun 22 '20 at 16:29
  • Thank you for the link. I tried this approach but it seems to me that for Fourier there is a fundamental trade off between the range of my p-space grid (max-value - min-value) and the resolution I can get in x-space. This is not the case for the code I gave above. Am I correct? Unfortunately, I know my numerical $f(p)$ only for a quite narrow grid and it has quite narrow features, so my grid spacing in p-space is quite small. The resolution in x-space is therefore really bad. Is there a way around that? – Luke Jun 22 '20 at 16:34
  • 1
    You cannot get more information than it is already contained in FFT. – yarchik Jun 22 '20 at 17:28
  • 1
    Then I guess it comes down to carefully choosing the discretization range, right? I have no experience with FFT, sorry :) – Luke Jun 22 '20 at 17:36

1 Answers1

0

We want to reproduce the result of following code with a FFT

Abs[FourierTransform[f[p], p, x, FourierParameters -> {a, b}]]^2

This code does it correctly for all positive values of the x-space:

 Clear[fft2]
fft2[fp_(*f(p) in Matrix form*), dp_(*grid resolution*), a_, b_] := 
 Module[{pgrid, samples, n, factor, transform, xgrid},
  pgrid = Range[fp[[All, 1]][[1]], fp[[All, 1]][[-1]], dp];
  samples = Interpolation[fp][p] /. p -> pgrid;
  n = Length[pgrid];
  factor = 
   Sqrt[Abs[b] n/(2 Pi)^(1 - a)] Exp[
     2 Pi I pgrid[[1]]/(n*dp) Range[0, n - 1]]*dp;
  transform = 
   Chop[Fourier[samples, FourierParameters -> {0, b}]*factor];
  xgrid = 
   RotateRight[Range[-Ceiling[n/2] + 1, Floor[n/2]]/(n*dp), 
     Floor[n/2] + 1]*2*Pi;
  Interpolation[Transpose[{xgrid, Abs[transform]^2}]]
  ]

It takes as argument $f(p)$ stored in the following form $fp=\{\{p_1,fp_1\},\{p_2,fp_2\}...\}$ and a proper grid spacing $dp$ depending on what resolution of the FFT you require, $a$ and $b$ as the Fourier Parameters and returns an interpolation function.

Not sure how I can make this work for negative x-values as well but for me they are not necessary.

Luke
  • 838
  • 4
  • 11
  • Negative x values should give the same result as their positive counterparts, since the integrands (summands, using Fourier) are complex conjugates of one another and an absolute value is taken at the end. – Daniel Lichtblau Jun 26 '20 at 15:56