Pseudo-spectral method (Split-step FFT)

There exist some physically interesting problems that have parts easily solvable in position space and parts easily solvable in momentum space. For these kinds of problems we can use a technique known as the pseudo-spectral method. Essentially, all this technique does is break the problem down into parts involving the bits easily solvable in each space and then merely makes the relevant transformations to move between spaces. Breaking the problem into parts usually involves an approximation to the actual evolution desired, but since numerical computation of the evolution PDEs is an approximation, it doesn't hurt to do another if it allows us to solve a problem, and we can be reasonably accurate as it turns out anyway. In this discussion I will again follow the text of DeVries.

Just to be different, we're going to analyse the one-dimensional Schrödinger equation instead of the diffusion equation^3.12. The reason for this is because the diffusion equation would reduce back to the standard FFT method described in the last section, and therefore wouldn't be very informative or interesting. The 1-D Schrödinger equation is

$$i \hbar \frac{\partial }{\partial t} \psi(x,t) = - \frac{\hbar^2}{2 m} \frac{\partial ^2}{\partial x^2} \psi(x,t) + V(x) \psi(x,t),$$ (3.63)

where $m$ is the mass of the quantum mechanical particle, $\psi(x,t)$ is its wavefunction and $V(x)$ is some potential dependent only upon the spatial position. We can write this in terms of a kinetic energy term $\mbf{T}$ and a potential energy term $\mbf{V}$ as

$$i \hbar \frac{\partial }{\partial t} \psi(x,t) = (\mbf{T}+ \mbf{V}) \psi(x,t).$$ (3.64)

Note that $\mbf{T}$ is basically the $\mbf{D}$ matrix (in general, the derivative operator) that we have been using in most of the examples to date. The potential $\mbf{V}$ merely multiplies $\psi$. One should mention at this point that the time dependence of Schrödinger's equation in Equation (3.63) can be solved formally, but this is really boring, and we want to do this on a computer anyway. Nevertheless, the formal solution can lead us to a practical method of solving Schrödinger's equation computationally. The solution of Schrödinger's equation in the form of Equation (3.64) is

$$\psi(x,t) = e^{-i(\mbf{T}+ \mbf{V}) \delta_t/\hbar} \psi(x,t_0),$$ (3.65)

where $\delta_t = t - t_0$ and $t_0$ is the initial point in time.

Note that this solution involves the exponential of operators. How do we take exponentials of operators? In much the same way we would take the exponential of a scalar;

$$e^\mbf{A}= 1 + \mbf{A}+ \frac{\mbf{A}\cdot \mbf{A}}{2} + \frac{\mbf{A}\cdot \mbf{A}\cdot \mbf{A}}{3!} + \cdots,$$ (3.66)

When we realise that $\mbf{A}$ could have a representation as a matrix, one can imagine that this could get nasty to calculate fairly quickly.

It would be really cool if we could write the exponential in Equation (3.65) as

$$e^{-i (\mbf{T}+ \mbf{V}) \delta_t/\hbar} = e^{-i \mbf{T}\delta_t/\hbar} e^{-i \mbf{V}\delta_t/\hbar},$$ (3.67)

but we can't because $\mbf{T}$ and $\mbf{A}$ are operators and in general don't commute. However, even though the above relation isn't true in a strictly mathematical sense, it doesn't mean it isn't useful as an approximation for what we want to do. Performing this decomposition of the exponential into two factors is known as the split-operator approach and is why this technique is also known as the split-step FFT method. Although it may seem a bit bizarre, this method is actually widely applicable, but we won't go into that here.

Equation (3.67) is only accurate to $O(\delta_t)$, nevertheless if we take the symmetric decomposition

$$e^{-i (\mbf{T}+ \mbf{V}) \delta_t/\hbar} \approx e^{-i \mbf{V}\delta_t/2\hbar} e^{-i \mbf{T}\delta_t/\hbar} e^{-i \mbf{V}\delta_t/2\hbar},$$ (3.68)

then the split-operator approximation is accurate to $O(\delta_t^2)$ (which is much better). What Equation (3.68) is saying is that we will work out the solution for half of the $\mbf{V}$ step, then all of the $\mbf{T}$ step then the other half of the $\mbf{V}$ step. The process we will follow is then:

1. Solve half of the $\mbf{V}$ step in $x$-space.
2. Transform into $k$-space.
3. Solve all of the $\mbf{T}$ step in $k$-space.
4. Transform back to $x$-space.
5. Solve the other half of $\mbf{V}$ in $x$-space.
6. Giving the final solution.

Right, let's do it then.

The wavefunction $\psi$ at time $t$ is approximated by

$$\psi(x,t) \approx e^{-i \mbf{V}\delta_t/2\hbar} e^{-i \mbf{T}\delta_t/\hbar} e^{-i \mbf{V}\delta_t/2\hbar} \psi(x,t_0).$$ (3.69)

Since $\mbf{V}$ is evaluated in $x$-space, we define the function

$$\phi(x) = e^{-i \mbf{V}\delta_t/2\hbar} \psi(x,t),$$ (3.70)

and so we want to evaluate the effect of the kinetic exponential acting on $\phi$. In other words

$$e^{-i \mbf{T}\delta_t/\hbar} \phi(x).$$ (3.71)

This is nasty in $x$-space, but easy in $k$-space, so we Fourier transform:

$$\Phi(k) = \mathcal{F} \bigl[ \phi(x) \bigr] = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty \phi(x) e^{-ikx} dx.$$ (3.72)

After a lot of hard work, one can eventually show that

$$e^{-i \mbf{T}\delta_t/\hbar} \phi(x) = \mathcal{F}^{-1} \Bigl[ e^{-i T(k) \delta_t/\hbar} \Phi \Bigr],$$ (3.73)

where we have defined the quantity

$$T(k) = \frac{\hbar^2 k^2}{2 m}.$$ (3.74)

Therefore, the full propagation is described by:

$$\psi(x,t) \approx e^{-i V(x) \delta_t/2\hbar} \mathcal{F}^{-1} \B... ...hbar} \mathcal{F} \Bigl[ e^{-i V(x) \delta_t/2\hbar} \psi(x,t_0) \Bigr] \Bigr].$$ (3.75)

And that's it, the solution evolves after jumping back and forth between $x$- and $k$-space, but nevertheless solves a reasonably difficult problem in an efficient and accurate manner.

Footnotes

... equation^3.12

Yay!