1
$\begingroup$

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

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?

$\endgroup$

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.