I need to solve system of two coupled partial differential equations numerically.

\begin{align} \frac{\partial x_1}{\partial t} &= c_1\nabla ^2 x_1 + f_1(x_1,x_2)\\ \frac{\partial x_2}{\partial t} &= c_2\nabla ^2 x_2 + K\frac{\partial x_1}{\partial t} \end{align}

The domain of system is a square region.

Boundary condition:

\begin{align} x &= \text{constant} \implies \frac{\partial x_1}{\partial x} = \frac{\partial x_2}{\partial x} = 0\\ y &= \text{constant} \implies \frac{\partial x_1}{\partial y} = \frac{\partial x_2}{\partial y} = 0 \end{align}

I tried to solve this system with Fourier transform. Solution becomes unstable after few iterations. I have solved this system earlier with finite difference scheme and it worked well so I know that constants of system are perfectly fine.

  • My question is can Fourier transform be used to solve these equations?
  • I read somewhere that it because of Neumann boundary condition one cannot apply Fourier transform. Is this correct?
  • If yes what is alternative?(I have read that cosine transform should be used but want to confirm).
  • $\begingroup$ How do you define $f_1(x_1,x_2)$ ? Is $f_1$ a periodic function ? $\endgroup$
    – ucsky
    Sep 10, 2012 at 16:35
  • $\begingroup$ $f_1(x_1,x_2) = P(x_1)arctan(x_2)$, where P(x) is a polynomial. $\endgroup$
    – chatur
    Sep 10, 2012 at 19:36
  • 1
    $\begingroup$ Maybe the numerical unstability come from the function $f_1$. If $P$ is not equal to zero at the limite of the domain you can have a jump at the boundary condition when reconstructing the periodic signal. $\endgroup$
    – ucsky
    Sep 10, 2012 at 20:07
  • $\begingroup$ @aberration, Thanks for comment but do not understand what it means. Probably I should study FFT more thoroughly. But If I am solving system in all the changes are limited to center part only(roughly) so that values at the boundary are not affected can I use FFT here? $\endgroup$
    – chatur
    Sep 11, 2012 at 6:01
  • 1
    $\begingroup$ Well maybe what I said is wrong, I think it's may work with some smooth function $f_1$ but anyway it's will give probably a wrong solution because near boundary process will not be solved properly. As you said you probably need to use a cosine modes or Chebyshev modes. If you really want to use Fourier transform you can use some kind of penalization method. $\endgroup$
    – ucsky
    Sep 11, 2012 at 16:08

1 Answer 1


The FFT can be used for periodic boundary conditions. Because von Neumann boundary conditions are effectively "mirror" boundary conditions, you have to do a "mirrored continuation", before you can apply an FFT. One drawback of this approach is that you will blow up the data volume by a factor 4 (which is not important if you are just interested in experimenting a bit). The use of the cosine transform implicitly does the "mirrored continuation" and avoids the factor 4 overhead.

Note that depending on where the grid points near the boundary are located, there are two different ways to do a "discrete mirrored continuation". Hence you will find that libraries like FFTW offer different variants of the cosine transform (corresponding to these different "discrete mirrored continuations").


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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