Boundary Conditions

So for we have chosen simulations that never reach the edge or the boundary. We have either limited the time or, as with the particle in a box, setup a potential that forces the wave function to zero at the boundary. For many examples, this is sufficient. However, in general, we have to deal with the wave function at the boundary of our simulation. This example shows a wave packet that interacts with the boundary almost immediately.

The interactions at the boundary clearly do not reflect reality.

Clearly the behavior at the boundary does not reflect physical reality. We can find the cause if we look at the update equations, and consider their evaluation at the ends of the grid.

To update the real part of the wave function
$Re (Ψ (x, t + Δ t)) \approx Re (Ψ (x, t - Δ t)) - \frac{Im (Ψ (x + Δ x, t)) - 2 Im (Ψ (x, t)) + Im (Ψ (x - Δ x, t))}{{(Δ x)}^{2}} Δ t + 2 V (x) Im (Ψ (x, t)) Δ t$ and to update the imaginary part we have
$Im (Ψ (x, t + Δ t)) \approx Im (Ψ (x, t - Δ t)) + \frac{Re (Ψ (x + Δ x, t)) - 2 Re (Ψ (x, t)) + Re (Ψ (x - Δ x, t))}{{(Δ x)}^{2}} Δ t - 2 V (x) Re (Ψ (x, t)) Δ t$

Remember that the terms like $Ψ (x, t)$ on the right side correspond to texture reads.

For simplicity focus on the right side of the grid. On the right edge, the problematic term is $Ψ (x + Δ x, t)$ , which reads beyond the edge of the grid.

Texture reads for computing

Ψ (x, t + Δ t)

at the right edge a 1D grid.

Remember that when we setup the textures, we set


  gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE);

So that when we read beyond the edge of the grid, we get the value at the edge. Clearly we need to do better for a realistic physical simulation.

This issue of what to do at the boundary of a simulation is one that recurs again and again in scientific computing. Indeed, the solution we will put forward borrows directly from the methods developed for electromagnetic fields.

We impose a wave equation at the boundary ¹, and numerically only treat the part that propagates from the interior region outward across the boundary. This is referred to as a one way wave equation.

Of course the wave equation is not an exact match for the Schrödinger equation. However, the errors introduced by using the wave equation along the boundary are small compare to the errors introduced without this treatment of the boundary. We will see an obvious improvement.

The one way wave equation is a close kin to the standard wave equation. We start with a standard wave equation, then factor it into operators. $\frac{\partial^{2}}{\partial x^{2}} Ψ (x, t) - \frac{1}{{v_{p}}^{2}} \frac{\partial^{2}}{\partial t^{2}} Ψ (x, t) = (\frac{\partial^{2}}{\partial x^{2}} - \frac{1}{{v_{p}}^{2}} \frac{\partial^{2}}{\partial t^{2}}) Ψ (x, t) = (\frac{\partial}{\partial x^{}} - \frac{1}{v_{p}} \frac{\partial}{\partial t}) (\frac{\partial}{\partial x} + \frac{1}{v_{p}} \frac{\partial}{\partial t}) Ψ (x, t) = 0$

Now, we can look at each of these differential operators separately. $(\frac{\partial}{\partial x} - \frac{1}{v_{p}} \frac{\partial}{\partial t}) Ψ (x, t) = 0$ which has solutions of the form $f (x + v_{p} t)$ , representing a leftward propagating wave. Impose this condition on the left edge of the grid, allowing only outward bound waves.

$(\frac{\partial}{\partial x} + \frac{1}{v_{p}} \frac{\partial}{\partial t}) Ψ (x, t) = 0$ has solutions of the form $f (x - v_{p} t)$ , representing a rightward propagating wave. Impose this condition at the right edge of the grid restricting the solution to rightward, outward, propagating waves.

Work through the application along the right edge in detail then carry what we learn over to the left edge. $\frac{\partial}{\partial x} Ψ (x, t) = - \frac{1}{v_{p}} \frac{\partial}{\partial t} Ψ (x, t)$

Write out expressions for the space and time derivatives at the edge of the grid as finite differences. $\frac{\partial}{\partial t} Ψ (x, t) = \frac{Ψ (x, t + Δ t) - Ψ (x, t)}{Δ t}$ $\frac{\partial}{\partial x} Ψ (x, t) = \frac{Ψ (x t) - Ψ (x - Δ x, t)}{Δ x}$

We might be tempted substitute these expressions directly into the one way wave equation, but we can't quite do that yet. Something a little bit subtle has happened here. These expressions are defined at different locations, and we need to do a little more work before we have an optimal solution.

The naive approach yields $\frac{\partial}{\partial x} Ψ (x, t)$ and $\frac{\partial}{\partial t} Ψ (x, t)$ at different locations.

Luckily there is a simple process to shift the derivatives onto a common location. Average the x derivative at t with a similar derivative at t+Δt. This shifts the derivative from $(x - \frac{Δ x}{2}, t)$ to $(x - \frac{Δ x}{2}, t + \frac{Δ t}{2})$ .

Similarly, evaluate the $t$ derivative at $x - Δ x$ . Averaging this with the derivative at $x$ shifts the derivative from $(x, t + \frac{Δ t}{2})$ to $(x - \frac{Δ x}{2}, t + \frac{Δ t}{2})$ .

This gives us new expressions for the derivatives. $\frac{\partial}{\partial t} Ψ (x, t) = \frac{1}{2} (\frac{Ψ (x, t + Δ t) - Ψ (x, t)}{Δ t} + \frac{Ψ (x - Δ x, t + Δ t) - Ψ (x - Δ x, t)}{Δ t})$ $\frac{\partial}{\partial x} Ψ (x, t) = \frac{1}{2} (\frac{Ψ (x t) - Ψ (x - Δ x, t)}{Δ x} + \frac{Ψ (x, t + Δ t) - Ψ (x - Δ x, t + Δ t)}{Δ x})$ .

Put these expressions into the one way wave equation for the right edge, giving a slightly long, but luckily not very complicated, expression. $\frac{Ψ (L t) - Ψ (L - Δ x, t)}{Δ x} + \frac{Ψ (L, t + Δ t) - Ψ (L - Δ x, t + Δ t)}{Δ x} = \frac{1}{v_{p}} (\frac{Ψ (L, t + Δ t) - Ψ (L, t)}{Δ t} + \frac{Ψ (L - Δ x, t + Δ t) - Ψ (L - Δ x, t)}{Δ t})$

Rearrange terms to get an expression for $Ψ$ at the right edge of the grid at the most advanced time in the simulation. $Ψ (L, t + Δ t) = Ψ (L - Δ x, t) + \frac{v_{p} Δ t - Δ x}{v_{p} Δ t + Δ x} (Ψ (L - Δ x, t + Δ t) - Ψ (L, t))$

A similar process gives us an expression for $Ψ$ at the left ( $x = 0$ ) edge. $Ψ (0, t + Δ t) = Ψ (Δ x, t) + \frac{v_{p} Δ t - Δ x}{v_{p} Δ t + Δ x} (Ψ (Δ x, t + Δ t) - Ψ (0, t))$

These update equations are significantly different from those we have seen before. In the next section we will see that they can be expressed elegantly with OpenGL geometry and shaders.

G. Mur, "Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations," IEEE Transactions on Electromagnetic Compatibility, vol. 23, no. 4, pp. 337–382, November 1981.