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.

Details

Snapshot 1: for only a few steps per rotation the trajectory is polygonal; the corresponding graphics for method sv give a more oblate, less circular trajectory
Snapshot 2: error curves for the same trajectory as snapshot 1
Snapshot 3: many steps per period create smooth error curves and trajectories
Snapshot 4: eigenvalues of the step matrix for a stable trajectory
Snapshot 5: eigenvalues of the step matrix for an exploding trajectory
This Demonstration builds trajectories by iterated application of the step matrix on states, starting from some initial state. In the present simple case, however, it is also possible to compute the generic
th
n
power of the step matrix as a symbolic expression. For the Störmer–Verlet method this was done in[1, 2].
For the densified asynchronous leapfrog method, corresponding results follow from the fact that the step matrices
U
sv
and
U
dalf
of the two methods, as given explicitly earlier, satisfy
U
dalf
=
-1
M
U
sv
M
with the invertible matrix
M=
2
-
0
1

. In particular, the eigenvalues of the step matrices are the same. Let us collect the main formulas that hold for
hω<2
:
n
U
sv
=
cos(nα)
βsin(nα)
-1
β
sin(nα)
cos(nα)
where
α=2arctan
γ
β
and
β=1-
2
γ
with
γ=
ωh
2
.
For the error functions
ν
and
ϵ
(for
ψ
0
=1
) this implies
ν
sv
=-
2
γ
2
sin(nα)
,
ϵ
sv
=-
1
2
2
γ
sin(2nα)
.
The similarity of the step matrices by matrix
M
implies the relations of the error functions for the two integration methods:
ν
dalf
=-
2
γ
β
ν
sv
,
ϵ
dalf
=
1
2
ϵ
sv
-2i
-2
γ
ν
sv
.
As carried out in[1, 2], all these equations remain valid when the real number
ω
is replaced by the Hamiltonian operator and the definition of the
ν
and
ϵ
functions is understood as implying the scalar products with respect to the initial state of a trajectory. If you want to see the reasons behind the particular step matrices, inspect the code that defines them as a product of simple building blocks in accordance with the basic leapfrog idea.

References

[1] U. Mutze, "The Direct Midpoint Method as a Quantum Mechanical Integrator." http://www.ma.utexas.edu/mp_arc/c/06/06-356.pdf.
[2] U. Mutze, "The Direct Midpoint Method as a Quantum Mechanical Integrator II." http://www.ma.utexas.edu/mp_arc/c/07/07-176.pdf.

Permanent Citation

Ulrich Mutze
​
​"On the Stability Limit of Leapfrog Methods"​
​http://demonstrations.wolfram.com/OnTheStabilityLimitOfLeapfrogMethods/​
​Wolfram Demonstrations Project​
​Published: September 7, 2011