8

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).
Gilles
  • 253
  • 1
  • 2
  • 10
chatur
  • 181
  • 3
  • How do you define $f_1(x_1,x_2)$ ? Is $f_1$ a periodic function ? – ucsky Sep 10 '12 at 16:35
  • $f_1(x_1,x_2) = P(x_1)arctan(x_2)$, where P(x) is a polynomial. – chatur Sep 10 '12 at 19:36
  • 1
    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. – ucsky Sep 10 '12 at 20:07
  • @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? – chatur Sep 11 '12 at 06:01
  • 1
    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. – ucsky Sep 11 '12 at 16:08
  • Just to clarify, for the boundary conditions do you mean that along the lines of constant $x$, the $x$ derivative of $x_1$ and $x_2$ is zero? – David Z Sep 13 '12 at 06:47
  • @DavidZaslavsky, Yes. So that flux of $x_1, x_2$ is zero in direction normal to surface at boundary. – chatur Sep 13 '12 at 08:04

1 Answers1

4

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

Thomas Klimpel
  • 2,141
  • 15
  • 35