in my quest to understand how I can use FFT to compute integrals (see my other question click, still no answer there), I came across the fact that a convolution of two functions can be calculated by multiplying the Fourier transforms of the functions and taking the inverse transform of that product.

Now suppose I want to calculate $$\int_{0}^R x \text{sin}(y-x)dx. $$ I can write this as the convolution of $\text{sin}(x)$ with $x\cdot \text{Ind}_{x<R}$ with $\text{Ind}$ as an indicator function that gives the proper integral boundaries.

My $x$ vector is given on N data points, equally spaced between 0 and some L, and to these N points I also know $\text{sin}(x)$.

The strange thing is, if $L=2\pi$, the method seems to work. But with another L, say $L=1$ it does not work. What is the reason for this, and what can I do?

I have read about the method of zero-padding where I add N-1 zeros to both function vectors. Would that help? And if so, I want to have a vector of the original size N after the ifft (which will give me a vector of 2N-1), so that I have the convolution value to each of my N $x$ values. Is there a way to do this?

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

pi2 = 2 * np.pi 

N = 2001   ## should be odd!!
R = 0.2
L = 1

x = np.linspace(0, L-L/N, N)
f = np.sin(x)
g = x.copy()
g[g>R]=0      ## indicator function

# f=np.pad(f, (0, N-1), 'constant')   ## zero padding??
# g=np.pad(g, (0, N-1), 'constant')

F = fft(f)
G = fft(g)
C = (L/N)*ifft(F*G) 

analytic_sol = -np.sin(R-x)+R*np.cos(R-x)-np.sin(x)   ## analytic result

plt.plot(x, C, label="test")
plt.plot(x, analytic_sol, label="analytic",linestyle="--")

enter image description here

Notice how the solution only becomes valid at x=0.2=R, i.e. when the indicator function becomes zero. What happens here?



Your Answer

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