# Calculate integrals using numpy.fft

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")
plt.legend()
plt.show()


• 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? – Anton Menshov Nov 27 '19 at 1:11
• Hint: $\sin(x) = \frac{e^{ix}-e^{-ix}}{2}$ – Alone Programmer Nov 27 '19 at 2:40
• 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}.$ – SchroedingersLion Nov 27 '19 at 17:10
• 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? – nicoguaro Nov 28 '19 at 3:34
• 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. – SchroedingersLion Nov 28 '19 at 16:12