Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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

share|improve this question
    
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
add comment

1 Answer 1

up vote 4 down vote accepted

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

share|improve this answer
    
many thanks, Christian! –  lovis Feb 10 '13 at 16:01
    
If one is computing a linear convolution of a Fourier series, then it can be advantageous to put the origin in the middle of the array, in which case the array of added zeroes need only be half the length of the input data. The library fftw++ (fftwpp.sf.net) provides some convolution routines which take care of the padding for you, which might be worth checking out. –  Malcolm Nov 23 '13 at 17:37
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.