WOLFRAM|DEMONSTRATIONS PROJECT

On the Stability Limit of Leapfrog Methods

​
number of periods
10.
steps per period
8.
method
dalf
sv
mode
ψ
ϕ
errors
eigenvalues
auto range
log yMax
-1
We study two leapfrog algorithms by applying them to the simple ordinary differential equation

t
ψ(t)=-iωψ(t)
with the initial value
ψ(0)=
ψ
0
. The exact solution obviously is
ψ(t)=
ψ
0
exp(-iωt)
for all real
t
and thus represents a rotation with angular velocity
ω
. In our context the significance of this equation is that it is the most general Schrödinger equation in a one-dimensional Hilbert space. Spectral decomposition allows any Schrödinger equation in a finite-dimensional Hilbert space to reduce to a finite set of copies of this case. The numerical solutions under consideration work in a two-dimensional Hilbert space formed by pairs
(ψ(t),ϕ(t))
, where the second component is defined to have the initial value
ϕ(0)=-iωψ(0)
. In this two-dimensional Hilbert space a time-stepping evolution algorithm is given in terms of a step matrix that in this case is a complex 2×2-matrix, the matrix elements of which are specified functions of the time step
h
. In our case we have:
ψ(t+h)
ϕ(t+h)
=
1-
2
τ
2
τ-
3
τ
4
-τ
1-
2
τ
2
ψ(t)
ϕ(t)
​
Störmer–Verlet (sv),
ψ(t+h)
ϕ(t+h)
=
1-iτ-
2
τ
2
-
3
τ
8
-2τ
1+iτ-
2
τ
2
ψ(t)
ϕ(t)
densified asynchronous leapfrog (dalf),
where
τ=ωh
. Here, as always in quantum theory, the angular velocity
ω
of phase rotation is associated with an energy
ℏω
. In the Demonstration code we use the numerical values
ω=
ψ
0
=ℏ=1
. The Demonstration lets you visualize the complex quantities
ψ
and
ϕ
for a discrete trajectory. Such a trajectory is defined by a number of rotation periods and a number of steps per rotation period. You can set both numbers using the sliders; also noninteger values make sense. For less than
π
steps per period one sees singularities for
ψ
and
ϕ
. Selecting the mode "eigenvalues" shows that then one of these eigenvalues of the step matrix has a modulus larger than 1. The quantities
ν(t)=〈ψ(t)|ψ(t)〉-〈ψ(0)|ψ(0)〉
and
ϵ(t)=〈ψ(t)|ϕ(t)〉-〈ψ(0)|ϕ(0)〉
vanish for the exact solution due to unitarity and energy conservation. They do not vanish exactly for the discrete trajectories defined by the two leapfrog algorithms. They can be inspected under the method "errors". The reason for considering these error measures is that they can also be computed in cases where the exact solution is not available.
To put our main finding into perspective we notice that having
π
steps per period should be compared with "sampling at Nyquist frequency", which is two steps per period. In engineering language we have to "oversample" by a factor
π/2≈1.57
to get a stable trajectory. As we can learn from this Demonstration, a much stronger oversampling is needed to get a trajectory that is moderately accurate.