Convolution error when using DFT for non-periodic functions

Question

I need help understanding why FFT-based computation of convolution on a finite domain, between two non-periodic functions often gets it wrong. Specific questions:

(1) why does FFT-based computation have error when interval [a, b] is short in the specific example below but no error when interval [a, b] is wider ;

(2) how can I anticipate this error when applying FFT-based approach to convolve functions?;

(3) Can this error be mitigated/reduced, e.g. different zero-padding options? Increasing discretisation size does not help at all.

Let's convolve two Gaussians on $x\in [a, b]$ for two choices of $[a, b]$: (1) $[a,b]=[0,5]$ and (2) $[a,b]=[-1, 5]$. We will compute convolution via DFT and compare against an exact result and a real-space evaluation.

Specific functions:

$$ f_1(x)=\exp{(-x^2)} $$

$$ f_2(x)=\exp{(-x^2)} $$

I want to compute the convolution on $[a, b]$ $$ f(x)=\int\limits_a^b f_1(x) f_2(x-y) dx $$

Implementation:

import numpy as np

import matplotlib.pyplot as plt

import scipy.fft as fft
from scipy.special import erf

N = 64
a = 0
b = 5
x = np.linspace(a, b, N)
dx = np.diff(x)[0]

f1 = np.exp(-x**2)
f2 = np.exp(-x**2)

# exact convolution
fe = .5*np.exp(-x**2/2)*\
        np.sqrt(np.pi/2)*\
        (-erf((2*a-x)/np.sqrt(2))+erf((2*b-x)/np.sqrt(2)))

# computation on real domain
c = np.zeros(x.size)
for i in range(x.size):
    c[i]=np.trapz(np.exp(-x**2)*np.exp(-(x-x[i])**2), x)


f3 = fft.ifft(fft.fft(f1)*fft.fft(f2)).real*dx

f31 = fft.ifft(fft.fft(np.concatenate((f1, np.zeros(N-1))))*\
               fft.fft(np.concatenate((f2, np.zeros(N-1))))
               ).real*dx

x1 = np.arange(x[0], x[0] + f31.size*dx, dx)

plt.figure()

plt.subplot(2,1,1)
plt.plot(x, f1, label = 'f1(x)');
plt.plot(x, f2, '.', label = 'f2(x)');
plt.plot(x, fe, label = 'Exact f1*f2 (x)');
plt.plot(x, c, '.', label = 'f1*f2 (x) computed directly');
plt.legend()
plt.title(f"Exact convolution on [a,b]=[{a}, {b}]");

plt.subplot(2,1,2)
plt.plot(x, fe,'.', label = 'f1*f2: exact result');
plt.plot(x1+a, f31, label = 'f1*f2: computed with FFT')
plt.axvline(a)
plt.axvline(b)
plt.legend()
plt.title('FFT computation of convolution');
plt.tight_layout();

Results for $[a, b]=[0, 5]$: NOTICE THAT FFT APPROACH HAS HUGE ERROR

Results for $[a, b]=[-2, 5]$: NOTICE THAT FFT APPROACH HAS ALMOST NO ERROR

Answer 1

This seems to be because the significant non-zero component of the Gaussian before the region of interest $[0,5]$ is neglected from the fft calculation but not in the direct calculations.

If I do the calculation over $[-10,10]$ but just look at the range of interest, it works fine.

Instead of using

$$f(x) = e^{-x^2}$$

in the DFT calculation, this windowed version

$$f(x) = e^{-x^2}(u[x-a] - u[x-b])$$

is being used (where $u[n]$ is the unit step).

---

Updated code below

import numpy as np

import matplotlib.pyplot as plt

import scipy.fft as fft
from scipy.special import erf

N = 64
a = 0
b = 5
aa = -10
bb = 10
x = np.linspace(aa, bb, 2*N)
dx = np.diff(x)[0]

f1 = np.exp(-x**2)
f2 = np.exp(-x**2)

# exact convolution
fe = .5*np.exp(-x**2/2)*\
        np.sqrt(np.pi/2)*\
        (-erf((2*aa-x)/np.sqrt(2))+erf((2*bb-x)/np.sqrt(2)))

# computation on real domain
c = np.zeros(x.size)
for i in range(x.size):
    c[i]=np.trapz(np.exp(-x**2)*np.exp(-(x-x[i])**2), x)


f3 = fft.ifft(fft.fft(f1)*fft.fft(f2)).real*dx

#f31 = fft.ifft(fft.fft(np.concatenate((f1, np.zeros(N-1))))*\
#               fft.fft(np.concatenate((f2, np.zeros(N-1))))
#               ).real*dx

f31 = fft.ifft(np.multiply(fft.fft(np.concatenate((f1, np.zeros(2*N)))),
               fft.fft(np.concatenate((f2, np.zeros(2*N)))))
               ).real*dx


x1 = np.arange(x[0], x[0] + f31.size*dx, dx)

plt.figure()

plt.subplot(3,1,1)
plt.plot(x, f1, label = 'f1(x)');
plt.plot(x, f2, '.', label = 'f2(x)');
plt.plot(x, fe, label = 'Exact f1*f2 (x)');
plt.plot(x, c, '.', label = 'f1*f2 (x) computed directly');
plt.legend()
plt.title(f"Exact convolution on [a,b]=[{a}, {b}]");

plt.subplot(3,1,2)
plt.plot(x, fe,'.', label = 'f1*f2: exact result');
plt.plot(x1+aa, f31, label = 'f1*f2: computed with FFT')
plt.axvline(a)
plt.axvline(b)
plt.legend()
plt.title('FFT computation of convolution');
plt.tight_layout();


plt.subplot(3,1,3)
plt.plot(x, fe,'.', label = 'f1*f2: exact result');
plt.plot(x1+aa, f31, label = 'f1*f2: computed with FFT')
plt.axvline(a)
plt.axvline(b)
plt.legend()
plt.title('FFT computation of convolution');
plt.xlim([a,b])
plt.tight_layout();