Discrete convolution

Question

please can I ask a bit stupid question? Let say I need to solve an equation in a form $\frac{\partial X}{\partial t}=\sum_k M_k * X_{n-k}$ How can I do the discrete convolution numerically? I will say I will expand X into series (FFT) and then, please? I have a layer of free surfaces on the top and bottom and from the sides the material can inflow/outflow.. Many thanks

Answer 1

In the continuous case on the whole real line, the convolution of two functions is equal to the inverse Fourier transform of the pointwise product of their Fourier transforms.

In the finite (discrete) case, it depends on the how you define the behavior of the convolution at the end points (as Dirk pointed out), since $n-k$ can be outside the range of indices for your vectors. There are two possibilities to prevent this (but since you didn't specify limits for your sum, I don't know which one you need):

1. You can "wrap around" the indices, i.e., for $n<k$ take $X_{N+n-k}$, where $N$ is the length of $X$ and $M$. This is the cyclic convolution, defined as $$ (M*X)_n = \sum_{k=1}^N M_k X_{(n-k)\ \textrm{mod}\, N},\qquad 1\leq n\leq N.$$ In this case, you can compute the convolution as in the continuous case: $$ M*X = FFT_N^{-1}\left(FFT_N(M)FFT_N(X)\right), $$ where $FFT_N$ is the (fast) discrete Fourier transform of length $N$.

2. You can "truncate" the indices, i.e., take $\max(1,n-k)$ instead of $n-k$. Without going into details, you now need to sum over a larger range; the linear convolution is defined as $$ (M*X)_n = \sum_{k=\max(1,n-N)}^{\min(n,N)} M_k X_{n-k}\qquad 1\leq n\leq 2N-1.$$ Note that the linear convolution of two vectors of length $N$ has therefore length $2N-1$. To compute this using the FFT, you need to pad the vectors: Let $(X;0_m)$ denote the vector $X$ padded by $m$ zeros. Then, $$ (M*X;0_1) = FFT_{2N}^{-1}\left(FFT_{2N}((M;0_N))FFT_{2N}((X;0_N))\right).$$