Discrete convolution

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

• In principle you are right. The details of the method depend on the details of the problem (especially boundary treatment). You may add more details to get a good answer.
– Dirk
Feb 8 '13 at 9:29

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).$$