Split-step FFT method

Another method of achieving the goal of stable propagation with high accuracy, is to use the so-called split-step FFT method, in which the differential operator is calculated in Fourier space. This technique involves the use of periodic boundary conditions, but has the advantage that the inverse differential operator is completely trivial in Fourier space!

The Fourier transform on a computer is a DISCRETE Fourier transform, as obviously one cannot Fourier transform over a continuous domain. These have many applications in computational physics, experimental data processing and image processing.

It is an especially useful method in higher dimensions, as inverting a Fourier transform is much faster than inverting a large multi-dimensional matrix. The FFT operator is usually defined for these discrete Fourier transforms, as

$$[\mathcal{F} A]_j = \sum^{L-1}_{l=0} e^{-2 \pi ilj/L} A_l$$ (B.6)

In the Matlab implementation (A(t_n)fft), where indices start at $l = 1$, the definition is actually:

$$\widetilde{A}_j = [\mathcal{F} A]_j = \sum^L_{l=1} e^{-2 \pi i(l-1)(j-1)/L} A_l$$ (B.7)

Note that the inverse FFT operation may require normalization by the number of lattice sites $L$, since the inverse FFT operator is defined for these discrete Fourier transforms, as

$$A_j = [\mathcal{F}^{-1} \widetilde{A}]_j = \frac{1}{L} \sum^{L-1}_{l=0} e^{2 \pi ilj/L} \widetilde{A}_l$$ (B.8)

The additional normalization is not required in the Matlab implementation of A(t_n+1/2)ifft since the normalization is carried out by the inverse Fourier transform subroutine, but is sometimes required in other implementations. The main point is that $[\mathcal{F}^{-1} \mathcal{F}] = 1$. Note, however, that the normalization of the FFT for a discrete Fourier transformation is not the same as the usual convention for a continuous Fourier transform, given earlier. This does not matter, as long as we only plot the results in $x$-space, not in $k$-space!

A more subtle point is the indexing convention with FFTs, which typically identifies $k=0$ (zero momentum) with the first lattice site in reciprocal space, say at $j = 1$ in Matlab. Positive $k$ values then correspond to $k = \Delta k (j-1)$, where the definition of $\Delta k$ is given by examination of the relationship between the discrete approximation to the Fourier transform, and the usual integral form; this leads to the results that: $\Delta k = 2 \pi /X$. The first negative value of $k = -\Delta k$ is then obtained by using the $l = L$ site, and reducing the lattice number from its highest value, to give increasingly negative $k$ values! The value of $k$ at $l = L/2$ is ambiguous, and can be given either a positive or negative value.

Thus, the total method for propagating a distance $\Delta t/2$ for the diffusion problem, can be written in just one step, as:

$$\mbf{A}' = \mathcal{P} [\mbf{A}(t_n)]$$ (B.9)

$$= F^{-1} \left[ \mbf{e}^{-\Delta t \kappa k^2/2} \mbf{F}\ \right]$$ (B.10)

To use this method more generally with nonlinear terms, requires one Fourier step like this, then a nonlinear step in the $x$-domain, then a further Fourier step. The total method can then be written in two steps, as:

$$\mbf{A}(t_{n+1/2}) = \mathcal{P} [\mbf{A}(t_n)] + \frac{\Delta t}{2} \mbf{U}[\mbf{A}(t_{n+1/2})]$$

$$\mbf{A}(t_{n+1}) = \mathcal{P} \left[ \mbf{A}(t_{n+1/2}) + \frac{\Delta t}{2} \mbf{U}\ \right],$$ (B.11)

Here, in general, the midpoint value of $\mbf{A}(t_{n+1})$ is to be solved for implicitly, since $\mbf{A}(t_{n+1})$ is present on both sides of the equation. This can be achieved simply in two or three iterations. In the Schrödinger equation--either linear or nonlinear--iteration isn't needed, since we can write:

$$\mbf{U}[\mbf{A}] = i \mbf{A}\Phi[\vert\mbf{A}\vert].$$ (B.12)

In this case, the evolution under the local term $\mbf{U}$--without the dispersion (Laplacian) term--is exact, because $\mbf{U}$ only affects the phase of $\mbf{A}$, and not its amplitude. Thus, one can treat the problem as three exactly soluble stages--i.e., Fourier propagation followed by local time-evolution, then Fourier propagation. A local discretisation error of size $\Delta t^3$ per step occurs as previously, because the three stages don't commute with each other:

$$\mbf{A}' = \mathcal{P} [\mbf{A}(t_{n})]$$

$$\mbf{A}(t_{n+1}) = \mathcal{P} \left[ e^{i \Delta t \Phi[\vert\mbf{A}'\vert]} \cdot \mbf{A}' \right].$$ (B.13)

There are alternative schemes to this one, based on Runge-Kutta methods, with better convergence properties--although these require more intermediate storage, resulting in caching problems with large lattices. Note that the most efficient way to use the split-step FFT is to only transform to the x-space variables at the mid-points, since this is the only time that x-space values are needed. A potential problem with this method is `aliasing' error, caused by the fact that periodic boundaries result in feed-through from one boundary to the other. This is usually not physically appropriate, though it can be eliminated by including a large absorption term (`apodisation') at the boundaries.