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="--")
plt.legend()
plt.show()
Notice how the solution only becomes valid at x=0.2=R, i.e. when the indicator function becomes zero. What happens here?
