Two finite difference algorithms were implemented in this simulation:
the velocity-Verlet method and the fifth-order
Gear predictor-corrector algorithm. The velocity-Verlet algorithm
involves the following steps:[19,1]
! Step 1: calculate acceleration at time t
acc = trapForce + interactionForce
! Step 2: calculate position at time t+dt, velocity at t+dt/2
pos = (pos) + (vel * dt) + (0.5 * acc * dt*dt)
vel = (vel) + (0.5 * acc * dt)
! Step 3: calculate acceleration at time t+dt
acc = trapForce + interactionForce
! Step 4: calculate velocity at time t+dt
vel = (vel) + (0.5 * acc * dt)
In the pseudocode above, acc, vel, and pos are 
arrays which store the x-, y-,
and z- components
of each ion's acceleration, velocity, and position. Note that the
velocity-Verlet algorithm needs only 3 arrays of dimension 
.
The fifth-order Gear predictor-corrector method, outlined below, is slightly more complex (but worth it, as we shall see!):
! STEP 1: predict position and up to fourth-order derivatives
! at time t+dt via Taylor expansions...
pos = pos+(vel*dt)+(d2*dt2F)+(d3*dt3F)+(d4*dt4F)+(d5*dt5F)
vel = vel+(d2*dt) +(d3*dt2F)+(d4*dt3F)+(d5*dt4F)
d2 = d2 +(d3*dt) +(d4*dt2F)+(d5*dt3F)
d3 = d3 +(d4*dt) +(d5*dt2F)
d4 = d4 +(d5*dt)
! STEP 2: evaluate forces on ions at time t+dt
acc = trapForce + interactionForce
! STEP 3: "correct" step 1's approximations by a scaled term
! which represents the discrepancy between the second derivative
! obtained from Taylor series expansion and the acceleration
! calculated explicitly from force.
del_d2=dt2F*(acc-d2)
pos=pos+a0*del_d2
vel=vel+a1*del_d2/dt
d2 = d2+a2*del_d2/dt2F
d3 = d3+a3*del_d2/dt3F
d4 = d4+a4*del_d2/dt4F
d5 = d5+a5*del_d2/dt5F
In the pseudocode for the Gear algorithm, d2, d3, d4, and
d5 are 
arrays which store the x-,
y-, and z- components of the second, third, fourth, and fifth
derivatives of ions' positions. These derivatives are used in a
Taylor series expansion:
and deld_2 , dt2F, dt3F,
dt4F, and dt5F
are factors of 
for scaling derivative terms.
Using this fifth-order Gear predictor-corrector algorithm requires
8 arrays of dimension 
, whereas the
velocity-Verlet algorithm required only 3 such arrays.
The increased memory requirement and complexity of the
Gear algorithm might, at first glance, render it less attractive
than the velocity-Verlet algorithm. However, the Gear algorithm
has an enormous advantage over the velocity-Verlet algorithm:
it requires only one calculation of the interaction force per time
step, while
the velocity-Verlet algorithm makes two calls to that function at
each update.
Recall that for this long-range interaction problem, the calculation
of interaction force is 
. Thus, the
speed of the velocity-Verlet algorithm suffers for large N.
Moreover, the Gear algorithm can achieve a higher degree of energy
conservation than the Verlet algorithm with a longer time step.
According to Allen,[1] ``based on energy conservation, Gear's
algorithm is about one order higher in accuracy than Verlet's.'' Penning
trap simulations on TrapApp illustrate the superior stability of
the Gear algorithm: one can readily construct a given configuration
and time step such that the velocity-Verlet algorithm will gain energy while
the Gear predictor-corrector algorithm remains stable.