Good evening,

I would like to understand why I do not get the correct result: I assume that I know my function on discrete data points and expand it as a discrete Fourier transform: $\text{sin}(x)=\sum_{k=-m}^m A_ke^{i2\pi kx}$

This I use to calculate an integral: $$\int_{-R}^{R}\text{sin}(x+y)\text{d}x = \sum_{k=-m}^m A_ke^{i2\pi ky}\int_{-R}^R e^{i2\pi kx} \text{d}x = \sum_{k=-m}^m\frac{A_k}{i2\pi k}\left( e^{i2\pi kR}-e^{-i2\pi kR} \right) e^{i2\pi ky}$$

In other words, I calculate $A_k$ from FFT(sin(x)), multiply each $A_k$ with the term $\frac{1}{i2\pi k}\left( e^{i2\pi kR}-e^{-i2\pi kR} \right) $ and go back into real space with inverse FFT of these modified $A_k$.

Yet I do not get the correct result (the analyitc solution is $-\text{cos}(R+y)+\text{cos}(-R+y)$ ). What am I doing wrong?

import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, fftfreq, ifft

pi2 = 2*np.pi

L = pi2
N = 1001   ## should be odd!!
R = 0.2

x = np.linspace(0, L-L/N, N)
f = np.sin(x)
FOU_f = np.zeros(N, dtype=complex)
phi = np.zeros(N, dtype=complex)
FOU_phi = np.zeros(N, dtype=complex)

FOU_f = fft(f, N)    ## function to fourier space

for k in range(1, N):
    FOU_phi[k] = ( FOU_f[k]/(1j*pi2*k) ) * ( np.exp(1j*pi2*k*R) - np.exp(-1j*pi2*k*R) )

phi = ifft(FOU_phi, N)  ## back to real space

analytic_sol = -np.cos(x+R) + np.cos(-R+x)

plt.plot(x, phi, label="test")
plt.plot(x, analytic_sol, label="analytic")

enter image description here

  • $\begingroup$ Currently, your ifft returns a complex result. A sign that something in your derivation is wrong. Are you sure in the correctness of step 1 in your formula of going from $\sin(x+y)$ to the sum? $\endgroup$
    – Anton Menshov
    Nov 27 '19 at 1:11
  • $\begingroup$ Hint: $\sin(x) = \frac{e^{ix}-e^{-ix}}{2}$ $\endgroup$ Nov 27 '19 at 2:40
  • $\begingroup$ I don't know where the theoretical error would be. You can always write down the fourier transform the way I did: If I have 2m+1 data points, I can write $f(x_i)=\sum_{k=-m}^m A_k e^{i2\pi kx_i}.$ $\endgroup$ Nov 27 '19 at 17:10
  • $\begingroup$ It is recommended that the numbers of points for FFT is a power of 2. Also, did you check the normalization of the FFT that you are using? $\endgroup$
    – nicoguaro
    Nov 28 '19 at 3:34
  • $\begingroup$ There is no mention of the power of 2 thing in the documentation: docs.scipy.org/doc/numpy/reference/routines.fft.html And the normalization is just 1/N in the inverse transform. $\endgroup$ Nov 28 '19 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.