Numerical Instability

Hopefully we have provided you with an understanding of the foundations behind numerical methods on the GPU. We even provided a straightforward example application that generated some pretty cool visualizations. It all seems too good to be true. Unfortunately, when something seems too good to be true, it probably is. This simple finite difference implementation is numerically unstable. Even with the most carefully chosen initial conditions, as the simulation runs over time, errors will become evident and grow without bound. If we aren't careful with the initial conditions, these errors will quickly dominate.

For example, if we double the xResolution from 800 to 1600 grid points, even with a flat potential, we get some, shall we say, spectacular results.

We can make significant improvements by revisiting our first approximation for the time derivative of the wave function.

$$\frac{\partial }{\partial t}\Psi \left(x,t\right)\approx \frac{\Psi \left(x,t+\Delta t\right)-\Psi \left(x,t\right)}{\Delta t}$$

This approximates the time derivative using the difference between the forward step and the current timestep, a technique known as a forward difference approximation. We can improve on this significantly with a central difference approximation. We will see, though that this comes at a cost in memory.

The central difference approximates the time derivative using the difference between the previous and the next timestep.

$$\frac{\partial }{\partial t}\Psi \left(x,t\right)\approx \frac{\Psi \left(x,t+\Delta t\right)-\Psi \left(x,t-\Delta t\right)}{2\Delta t}$$

Generally, we expect this to yield more reasonable results. We can think of these difference equations as approximating the derivative at the midpoint between the two points. So the forward difference really approximates the derivative at $t+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\Delta t$ when we really want the derivative at $t$ . This is what the central difference provides.

This time, we put the central difference approximation into the the Schrödinger equation and isolate the $t+\Delta t$ term on the left.

$$i\Psi \left(x,t+\Delta t\right)\approx i\Psi \left(x,t-\Delta t\right)-\frac{{\partial }^2}{\partial x^2}\Psi \left(x,t\right)\Delta t+V\left(x\right)\Psi \left(x,t\right)2\Delta t$$

This looks almost exactly like our first path through the schrodinger equation, except for the first term on the right side, where instead of $\Psi \left(x,t\right)$ we have $\Psi \left(x,t-\Delta t\right)$

Remember that the terms on the right all became texture reads in the code. This equation shows that we will now need to read, and therefore keep, wavefunction values for $t-\Delta t$ . The wave function for $t+\Delta t$ becomes the wavefunction from $t-\Delta t$ plus adjustments from $t$ . This gives raise to the name of the method, leapfrog, which we will cover next.

• Back
• Next