I have the following functions: $$\phi(x) = \exp{\left(-i\frac{kx^2}{2f}\right)}$$ and $$f(x) = \exp\left(-\frac{x^2}{2\sigma^2}\right)$$ for some given values of $k$, $f$, and $\sigma$.
I am interested in computing the product $f \cdot \phi$ in a numerical solver. I have 2 approaches to do this. The first one is the obvious one:
- take the two functions and multiply them for some values on a discrete equally spaced numerical grid if there is no alias.
- if there is some alias for "large" values of $x$ the go to Fourier space and compute using $F^{-1}(F[\phi](\xi)\star F[f](\xi))$ where $F$ is the Fourier transform $F[f](\xi) = \int_{-\infty}^{\infty} f \cdot \exp(-2i\pi x \xi)$, $\star$ is the convolution.
The first case works, naturally. The second one has issues.
Since the 2 functions allow for analytic computation of the Fourier Transform, I assumed that computing the analytic functions $F[\phi]$ and $F[f]$ and then evaluating it on a discrete domain would be the way to go. In order to have the correct answer for the Fourier Transform, I have used Wolfram Alpha, which returned $$F[\phi](\xi) = \sqrt{\frac{2\pi f}{ik}}\exp\left(\frac{2i\pi^2 f \xi^2}{k}\right)$$
All good so far, except when I what to compare the outputs for both approaches. For the direct product approach, the $L_2$ norm remains the same as for $f$ since $|\phi(x)|=1, \forall x\in \mathbb{R}$. For the second approach, the resulting product is at least multiplied by a constant. I say at least because I am not certain at this point what the difference is.
In order to check myself I have considered computing the FFT of $\phi$ and use that in the calculation. The amplitude of the analytic and numeric Fourier Transforms are different.
As it can be seen in the figure (plot of the amplitude of the analytic and numeric Fourier transforms), the blue (analytic) and orange (numeric) lines are not the same. The phases are similar as it can be seen in the second figure where the same color legend has been used.
The strange part, for me, is that the numerical approach gives a better answer than the analytic implementation. Why does this happen?
Apparently the problem is related to the discrete and finite spatial domain. Given the shape of $\phi$, when the FFT is computed numerically the resulting function goes to the theoretical one, but since the domain is a closed interval it is only an approximation.
I don't think that there is a solution to this problem as long as a numerical approach is used.