# Solving a simple Schroedinger equation with Fast Fourier Transforms

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:

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

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?

• Your closed-form solution isn't true in general — it's only valid for $k$ that are integer multiples of $2\pi/L$, so it's not the right solution if your initial condition is non-periodic (i.e., has a discontinuity at the boundary). Perhaps try deriving the correct solution accounting for the discontinuity and comparing the numerical solution to that? Even then, you'd still get numerical artifacts from the discontinuous IC. – Kirill Jan 14 '16 at 21:08
• @Kirill is right. The issue is that your initial conditions have a discontinuity at the (periodic) boundary. The discontinuity is represented by large-magnitude high-frequency spectral components which eventually spread and corrupt your entire solution. You can't really say it doesn't match the "correct" solution because the solution you were expecting doesn't have discontinuous initial conditions. – Doug Lipinski Jan 14 '16 at 21:17
• This question was already answered on MO. If you want it here, it should be migrated, not reposted. – David Ketcheson Jan 15 '16 at 20:09

## 1 Answer

The fundamental solution of the heat equation $$u_t = iu_{xx}$$ with periodic boundary conditions on $[0,1]$ can be expressed in closed form in terms of a Jacobi theta function, $$\vartheta_3(x, -4\pi t) = 1 + 2\sum_{n\geq1} e^{-4\mathrm{i}\pi^2n^2t}\cos 2\pi n x.$$

So if the initial conditions are discontinuous in the domain (as in the case of a non-periodic IC $e^{ikx}$ with $k$ not an integer multiple of $2\pi$), the solution you should be comparing your numerical results to is not the single wave $e^{i k x-ik^2t}$, which only satisfies the boundary conditions for those specific $k$, but the convolution with the fundamental solution, $$\int_0^1 u_0(y)\vartheta_3(x-y, -4\pi t)\,dy.$$