1
$\begingroup$

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()

enter image description here

$\endgroup$
6
  • $\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, 2019 at 1:11
  • $\begingroup$ Hint: $\sin(x) = \frac{e^{ix}-e^{-ix}}{2}$ $\endgroup$ Nov 27, 2019 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$
    – reloh100
    Nov 27, 2019 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, 2019 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$
    – reloh100
    Nov 28, 2019 at 16:12

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.