While trying to solve a stochastic Gross-Piaevskii equation I have found a problem that can be tracked down to something buggy occuring in the simplest Schrodinger equation possible:

$\partial_t \psi = i \partial^2_x \psi$

This has the very simple solution of a plane wave:

$\psi(t) = A e^{i k x - i\omega t}$, with $\omega = k^2$

and my problem is that in some cases I cannot recover this very simple solution numerically. This is due to the fact that my code uses a very standar Fast Fourier Transform method to go to k-space, where the evolution is trivial, and then does the inverse FFT to go back to real space. This method naturally enforces periodic boundary conditions in your simulation grid since you are expanding your solution as a sum of periodic functions. What occurs is that, when the initial condition is a plane wave with a periodicity that does not match the domain size, the algorithm is not able to provide the correct solution.

I show here an example. The problem is clear in $|\psi(t)|$, which should remain constant:

Evolution of the wavefunction amplitude

The effect is less dramatic in the phase but still present. One can see some ripples that come from the bottom edges.

Evolution of the phase

I'm not surprised by the fact that, if I enforce periodic boundary conditions, computing the evolution of a wave which is not periodic in the simulation domain yields problems. However, the method is very standard and use extensively, with packages like XMDS (xmds.org) employing it by default. Therefore, I am surprised by the fact that I did not found any mention of this method failing to solve such an extremely simple example. My question is, is there something that I'm missing? Should I just get over it and assume that this is not a good method if I expect a solution of this kind? Does someone know a reference where this is documented?

  • 1
    $\begingroup$ Maybe you could use some kind of regularization on the space you are already working. Or perform the FFT on another (resampled) space. $\endgroup$ – mathreadler Jan 14 '16 at 21:07
  • $\begingroup$ This is a very interesting question and I am sorry it did not find attention. AFAIK, the numerical simulations of "oscillating" equations is hard even in simple cases, like you showed in this question. I saw in a course of Edriss Titi that he used to make fun of the audience, asking them to numerically find out whether a certain simple initial value problem blew up or not. Everybody would answer that the solution blows up, because they simulate it with periodic boundary conditions, but actually the problem posed on free space does not blow up. $\endgroup$ – Giuseppe Negro Sep 21 '18 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.